Contact Matrices for Multipopulation Epidemic Models: How to Build a Consistent Matrix Close to Data ANDREA PUGLIESE
Dipartimento di Matematica, Universitd di Trento, 38050 Povo (TN), Italy [Received 23 January 1991 and in revised form 26 June 1991]
Keywords: multigroup epidemic models; mixing matrices; contact networks.
1. Introduction
In models of the dynamics of sexually transmitted diseases in several heterogeneous populations, it is necessary to specify the matrix of the contacts between populations. Until recently, the proportionate mixing model was the most common description of the mixing between populations, mostly because of its simplicity. In fact, in this case the computation of many quantities of interest, such as the reproduction number Ro, can be performed using a single population model with a contact rate that keeps track both of the mean contact rate and of its variance [2] (for a more general discussion, see also [9]). Nold [16] and Hethcote & Yorke [12] introduced another form of mixing model, the so-called 'preferred mixing'; in this model, it is assumed that a fixed fraction of contacts takes place within the population to which an individual belongs, while the remaining fraction follows the proportionate mixing model. It became apparent [14] that the dynamics of the model with preferred mixing is rather different from the dynamics of models with proportionate mixing. It is also clear that the preferred mixing model is inadequate to describe several patterns of contacts, and that other types of preference models may be necessary. 249 © Oxford Univenily Pren 1991
Downloaded from http://imammb.oxfordjournals.org/ at Cambridge University on August 23, 2015
In models of the dynamics of sexually transmitted diseases between and within N interacting populations, it is necessary to specify the matrix of contacts between populations. Mixing matrices have to satisfy a consistency condition that generally will not be satisfied by empirically obtained matrices. The problem of inferring a mixing matrix from data is phrased here as the problem of finding a matrix in a prescribed set that minimizes an opportune distance from the matrix of data. Two different distances that attempt to measure relative errors are suggested. When no constraints are posed on the activity rates of the populations, the author shows that the minimum distance from data is attained at Knox's (1986) matrix. When activity rates are to be preserved, the minimum cannot be found explicitly for N > 2. An algorithm proposed by Area, Perucci, Spadea, and Rossi (1990) is then investigated, and shown always to converge to a consistent matrix. Through several examples, it is shown that this limiting matrix does not minimize distance from data, but is generally close to the minimum. Finally, the author simulated the collection of data with sampling errors and possible bias and evaluated the performance of this algorithm in approximating the 'true' contact matrix starting from the simulated 'data' matrix.
250
ANDREA PUGLIESE
The problem, however, is that mixing matrices cannot be arbitrary. In fact, once the number of individuals, Nh and the mean level of sexual activity, a(, are assigned for each population i, the mixing matrix ptJ must satisfy (1)
*UTJ>
(2a)
'~l~T^T
0}, and let Aj = {C e M(n xn):Ctj>0
if (i,j) e / ,
CtJ = 0 if (i,j)
Suppose now that C and D are in Aj and let
da(C,D) = ( X llogCCy/zyrY.
(5)
Here dx is a distance in Aj for all a ^ 1. Formally, we also consider a = oo intending the maximum over (ij) e J.
Downloaded from http://imammb.oxfordjournals.org/ at Cambridge University on August 23, 2015
J. (4b) They found empirically that this procedure always converged to a limiting matrix p that satisfied (1); then, at each time step of their simulations, they updated p according to the current values of Nh using algorithm (4). Note that the problem of constructing consistent mixing matrices from inconsistent ones presents itself at each step of a dynamic simulation of a multipopulation epidemic model. Even if the initial mixing matrix p is consistent, when population sizes TV, change with time, it will no longer be consistent. To perform a simulation, therefore, one needs a rule to make consistent, at each time step, a mixing matrix p, which will itself possibly vary according to predictable behavioural changes. An empirical procedure for achieving this aim, also with modifications of activity rates, has been used by Anderson and co-workers [1, 11]. Another option, chosen by Blythe et al. [6], is to assume that the matrix
0 for i = 1 ,..., N.
This assumption is further discussed in Section 4. We also change the distance used in Aj, in agreement with the procedure of Area et al. [3] and with the greater role that the rows of C play, since we force the sum of each row to a prescribed constant. More precisely, let Jt = {j: (i,j) e J) and note that, because of (A), Jt is nonempty for i = 1 ,..., N. Then, if C and D are in Ay, let
d'x(C, D) = \Y, (maxllogCCyD^lYT. Note that dx = d'x for a = oo only. Finally, let =
{CeAy:CiJ=CJi}.
(6)
Downloaded from http://imammb.oxfordjournals.org/ at Cambridge University on August 23, 2015
PROBLEM 1 Find the symmetric matrix C e Aj such that dJjC, C°) = min da(C, C°) over all symmetric matrices C e Aj.
CONTACT MATRICES FOR EPIDEMIC MODELS
253
We can then state the problem of finding the consistent contact matrix closest to the assigned one, as follows. PROBLEM
matrices
2 Find the matrix CeAnB CeAnB.
such that d'a(C, C°) = min d'x(C, C°) over all
Finding a solution to Problem 2 appears to be rather difficult. When n = 2, the solution (unique for a > 1) can be found, but its explicit expression depends on the value of several inequalities. For instance, setting p = p°l2 < 1 and q = p21 < 1, if T\V > T2q, p ^ i and T2(\ -q)^
Txp, then C 1 2 = {TYpT2q)112
for all values of a > 1.
PROBLEM 2a Given a matrix C e B, find the matrix De A such that dx(C, D) = min da(C, D) over all matrices D e A.
2b Given a matrix D e A such that Da > 0 for i = 1 ,..., N, find the matrix CeB such that d'x(D, C) = mind'a(D, C) over all matrices C e B.
PROBLEM
Problem 2a is just Problem 1; therefore we again have Knox's solution. As for Problem 2b, we have the following result. PROPOSITION
1 / / 1 ^ a < oo, there exists a unique solution of Problem 2b:
^ (7) Proof. Note first that all elements in Aj, and therefore in A and B, have strictly positive diagonal elements, because of assumption (A). Clearly the value of C,y influences only the ith term in the summation in (6). Therefore, for any fixed i, the values of CiJt with j e Jt, will be such as to minimize max|log(Cy/Dy)|
(8)
under the constraint that Ci} are nonnegative and sum up to Tt. To simplify the notation, let Jt = {i,jl ,...,],}, where 0 ^ r < n — 1, and rename xk = ClJk, yk = DlJk (k = 1 ,..., r), y0 = Dih and T = T,. The minimization of (8) can be restated as the problem of finding positive x, ,..., xr, such that X*-i x*
2, there does not seem to be any explicit formula. It is certainly possible to use numerical minimization techniques; however, they only guarantee convergence to a local minimum, not necessarily a global one. The procedure of Area et al. [3] provides an hybrid solution to this problem. In fact, each step can be read as follows: given a matrix C e B find the matrix DHe A that minimizes dx(C, £>„); then, find the matrix C + 1 eB that minimizes d'x(D", C + 1). One may hope that this procedure will converge to a matrix C = D that will be in A n B (this will indeed be proved in the following section) and that this limit C will conserve some minimizing property (which does not seem to be generally true). In order to show that indeed algorithm (4) has the above property, we study the solution to the following subproblems.
254
ANDREA PUGLIESE
minimize F(xi,..., xr) = max] max |log(xj/.y,)|, log 0-1
(9)
r
The matrix C of the thesis then corresponds to
Note that we have log
= F(Xl ,., xr)
for all j = 1,..., r.
Now take x ^ x, choose _/ with Xj ¥= Xj, and suppose, for instance, X*=o ^t ^ T. If x^ < Xj, we have
,..., xr)
= F(x t ,..., x r ).
If Xj > Xj and x t > xk for all /c = 1 ,..., r, we have T — £ l = j x t < T — ^ ^
t
xk, and
so log
= log log
>log
= F(x!,..., x r ).
In any case, we have therefore F(xx ,..., x r ) > F(x t ,..., x r ), that is, the thesis. When a = 00, C,j will minimize (8) only for j such that log
=
max
/
T.
When i does not have this property, Ctt can be chosen somewhat arbitrarily. Note that if in Problem 2a we substitute d'x for da, or in Problem 2b dx for d'a, the solution is not necessarily unique, and certainly has no simple expression. The choice of this hybrid step is therefore mainly dictated by simplicity. 3. Convergence of the approximating algorithm We now consider the convergence of algorithm (4). It can be rewritten as follows: given a probability matrix p° satisfying assumption (A), we define the symmetric matrix C° by "0 _ y —
(10)
Then, given a symmetric matrix C, let (ii)
Downloaded from http://imammb.oxfordjournals.org/ at Cambridge University on August 23, 2015
ii-oy
CONTACT MATRICES FOR EPIDEMIC MODELS
255
and, using (7) and (10), we obtain
a;1 = (F(c-)) y .= c&o/d? dp1*2.
(12)
Let S = {C symmetric nonnegative matrices} and n = {CeS
such that dk > 0 for all / = 1,..., N } ,
THEOREM
2
Matrix C converges (as n -* oo) to a matrix C such that
£ Cy = 7J V i = 1 ,..., JV. For its proof, we need several lemmas. Most of the results can be extended to the case where the geometric mean (xy)112 is replaced by a generalized mean m(x, y). We start with some simple but useful properties of F. LEMMA 3 Let CeS such that Cu > Ofor i = 1,..., N; let C = F(C). Without loss of generality, assume dl = maX|_j N dt and dN = min ( = 1 N dt. Then 12
'.
12
The equality {djd^ only if dt = dN.
(13) u2
= 3, holds only ifdt = dx; and the equality (d,/dN)
max
