THEORETICAL

POPULATION

BIOLOGY

11,

23-36

(1977)

On a General Stochastic Epidemic Model NIELS

Department

G.

BECKER*

of Mathematics, La Trobe University, Bundoora, Victoria 3083, Australia Received May 5, 1976

A non-Markovian epidemic model is proposed for which a stochastic epidemic threshold theorem, like that by Whittle (1955) for a simpler model, is shown to hold. The threshold theorem can be a useful guide in determining public health measures aimed at preventing major outbreaks of a communicable disease. Bounds are obtained for the mean size of minor epidemics. A parameter of the model, whose value is crucial to applications of these results, is the product of the infection rate and the mean duration of the infectious period. Estimators are suggested for this parameter by identifying martingales associated with counting processes of the epidemic model and by constructing other martingales which can be written as stochastic integrals.

1.

INTRODUCTION

According to mathematical theory, the epidemic threshold theorem of Kermack and McKendrick (1927), and in particular its stochastic form as given by Whittle (1955), is fundamental to the prevention of major epidemics. The theorem indicates that major epidemics can be prevented in homogeneously mixing communities, by keeping the product of (i) the size of the susceptible population, (ii) the infection rate, and (iii) the mean duration of the infectious period, sufficiently small. The form of the stochastic threshold theorem can be deduced by using a suitable simple birth-death process to describe the size of the infective population and interpreting extinction of the process as a minor epidemic. It is tempting to retain this interpretation while describing the size of the infective population by a general branching process and thereby obtain a threshold theorem under more realistic assumptions. The objection to this is that a finite number of cases does not necessarily imply a minor epidemic. A minor epidemic is more properly defined as one in which at most a specified fraction of the susceptible community is affected. Using this definition, it is here verified that the stochastic epidemic threshold theorem does hold under more * This work was conducted at the Center for Environmental Quality Management, Cornell University, and was supported in part with funds provided by the Zurn Founda-

tion.

23 Copyright All rights

0 1917 by Academic Press, Inc. of reproduction in any form reserved.

ISSN

00404809

24

NIELS

G. BECKER

realistic assumptions. By “sandwiching” the epidemic model between two general branching processes, bounds are obtained for the mean size of minor epidemics. The estimation of certain parameter values needed for the application of the theorem is then considered.

2. AN EPIDEWC

MODEL

FOR INFINITELY

LARGE COMMUNITIES

At a given time, chosen as the time origin, a individuals who have just been infected with a certain communicable disease join a community of individuals who are susceptible to the disease. For the initial formulation of the model the susceptible community is assumed to be infinitely large. It is also assumed that no further infected immigrants arrive during the course of the epidemic. Each infective has the potential to spread the disease until he is removed as a result of death, recovery, or isolation. After his infection, an individual is latent for a random duration S. He is then infectious until his removal after a further random duration T. The probability distribution of (S, T) is assumed to be the same for each of the infectives. During his infectious period, the individual infects susceptibles, independently of other infectives, according to a Poisson process K(v) with intensity /3. The infection intensity /3 is also assumed to be the same for each of the infectives. Hence, with each infective Y there is associated a random triple (S, , T, , K,.(.)) and the triples are assumed to be stochastically independent. This epidemic model is therefore a general branching process. For such a process it is known (see Theorem (6.5.1) of Jagers (1975)) that the extinction probability @) is the smallest nonnegative root of the equation s = qsm+-j

= zqw-l)]

= MT@(S -

I)],

where MT is the moment generating function of the infectious period. Except in the degenerate case &‘#(s - l)] = s, it follows that q(p) < 1 if and only if VP > 1, where Y = n/r,‘(O) is the mean duration of the infectious period. The distribution of T is seen to play a major role in this result and it is appropriate to indicate how it can be specified in more detail. Firstly, Pr(T = 0) = p might be positive due to the chance of an infective being removed during his latent period. Conditional on T > 0, an infective can be removed in two possible ways, The infection within an infective, who is not discovered by public health measures, is assumed to evolve naturally and the infective then becomes a removal by death or recovery after a random time U. However, an infective might be removed prior to his “natural” removal by being discovered and isolated. Given the duration U of an infective’s “natural” infectious period, the intensity of removal by public health measures is taken to be y(s ) U), if the individual has

STOCHASTIC

25

EPIDEMIC

been infectious for a time s. The distribution of U is determined by the development of the infection within an individual, while the intensity function y(. 1 U) is determined by the type of symptoms infectious individuals display and the nature of public health measures. The distribution of the duration V from the time an individual becomes infectious until he is discovered is given by Pr(V < a 1 U) = 1 - exp I-

joV y(s 1 U) ds( ,

and the distribution for U. An infective might never be discovered and so V need Conditional on T > 0, the duration of the not have a proper distribution. infectious period is given by T = U A V, the smaller of U and V.

TABLE

I

The Number of Smallpox Importations into Europe over the Period January 1961 to September 1970, Which Resulted in k Cases ”1

k Observed

number

2

3

4

5

>5

12

3

2

1

1

9

Estimated expected number (geometric)

14.5

3.6

1.8

1.1

0.8

6.1

Estimated expected number (Poisson)

11.1

4.1

2.2

1.5

1.0

8.1

One might expect this model to have some practical relevance to large communities. That the branching process can be used to describe minor epidemics and the early stages of major epidemics in large populations has been suggested previously, for example by Bharucha-Reid (1956) Neyman and Scott (1964) and Bartoszynski (1967). Use of this idea has been made by Becker (1972, 1975) in proposing methods for assessing vaccination policies. It seems appropriate to present some empirical evidence to support the adequacy of a branching process as an epidemic model for small outbreaks in large communities. Table I summarizes the sizes of smallpox outbreaks as a result of 28 importations into Europe over a period of 10 years, as given by the World Health Organization (1970). We assume that no infectives are removed during their latent period and consider two different distributions for the duration T of the infectious period. First, it is assumed that T has a negative exponential distribution, in which case

MT[P(S- 111= l/(1 + 4 - $4

26

NIELS G. BECKER

and the offspring distribution is seen to be a geometric distribution. assumed that T = Y with probability one, in which case M&3(s

-

Second, it is

l)] = eVE+l)

and the offspring distribution is seen to be a Poisson distribution. Both of these offspring distributions belong to the family of one parameter power series distributions. The probability distribution for the total number of infections follows, in each case, from Eq. (2) o f section 4 and the results of Patil(l963) on the sum of independent variates distributed according to a power series distribution. Each of the two resulting distributions involve one parameter, whose value needs to be estimated. For this purpose we note that under the assumption of a homogeneous general branching process, it follows that the 28 separate outbreaks are equivalent to one outbreak started by 28 infectives. Now the results of Becker (1974) show that the maximum likelihood estimator of the unknown parameter is, for both cases, obtained by equating v/3 to 1 - a/R@), where R(p) denotes the total number of cases. Here a = 28 and the observation on R(b) is 391. The resulting (estimated) expected frequencies, corresponding to each of the geometric and Poisson offspring distributions, are also given in Table I. The agreement is seen to be good. Specifically, chi-square goodness-of-fit tests based on the classes {l}, (2, 3,4, 5}, and (6, 7,...} lead to test values 1.8 (geometric) and 0.5 (Poisson) which are not significant when compared with the critical value

xl”(.05) = 3.8.

3. THE EPIDEMIC

MODEL

FOR A FINITE

COMMUNITY

Consider now a finite closed population. At time t, the community consists of X, + Yt infectives, and 2, removals. Here X, denotes the number of infectives who are in their latent period, while Yt denotes the number who are in their infectious period, at time t. The initial conditions are W,, = n, X0 = a, Y,, = 2, = 0 and the assumption of a closed population gives

W, susceptibles,

WS + Xt + Yt + 2, = n + a for all t 2 0. The only change made to the formulation of Section 2 is that the infection intensity is now made to depend on the number of susceptibles remaining in the community. This adjustment is made to allow for the depletion of the susceptible population. When W, = w the infection intensity is now written as fiW and we assume /I to be a nondecreasing function with &, = 0. By allowing the infection intensity to depend on W, , the infection processes K,(e) become stochastically dependent and the epidemic model is no longer a general branching process. The infection process K,(e) is now what is called a compound regular point process in the literature; see Rubin (1972).

27

STOCHASTIC EPIDEMIC

4. THE STOCHASTIC EPIDEMIC

THRESHOLD

THEOREM

For a specified fraction E, 0 < E < 1, we say that the community suffered a minor epidemic if the total number of secondary infections, that is all cases except the a initial cases,is less than no. Let 7rCdenote the probability of a minor epidemic. Since fi is a nondecreasing function it follows for all minor epidemics that

for all t > 0. It is intuitively clear, therefore, that the probability of a minor epidemic has the bounds given by

where R(b) denotes the total number of progeny in the general branching process of Section 2 with infection intensity fi. Hence

We indicate that the term Pr{ne + a < I?@,) < co} becomes negligible for large communities. Dwass (1969) has shown that WR(l%) = 4 = (u/k) WQl + Q2 + *-a + Qk = k - 4,

(2)

where Q1 , Q2,..., Qkare independent and identically distributed discrete random variables, which, in the present application, have the distribution corresponding to probability generating function AIr[j3,(s - I)]. This result of Dwass was in fact derived earlier by Lloyd (1963) and Mott (1963) in connection with an equivalent problem in dam theory. Now using the central limit theorem for sums of independent and identically distributed random variables we have to a firstorder approximation Pr k-u--~ 1, thenMBnlla< T G MBn-nc)la.

6)

If &-,,

(ii)

If v/$-~~ -=I 1 < v/3, , then [q(&#

(iii)

If I&

< rc < I.

< 1, then 7~~= 1.

It seems reasonable to assume that is,, ‘v j3n-,B for large communities and we see that the probability of a minor epidemic is roughly [q@J]“. Under the assumption that infectives do not have a latent period and that the infectious period is ewonentially distributed, we find the above theorem to be equivalent to the stochastic epidemic threshold theorem of Whittle (1955).

5. THE MEAN SIZE OF MINOR EPIDEMICS A similar argument leads to bounds on the mean size of minor Let 7 denote the duration of the epidemic. That is

T = i${t: x, + I’, = O},

epidemics.

(3)

so that 2, denotes the total number of cases. If G denotes the probability generating function for the total number of cases in the general branching model of Section 2 with one initial infective, then G(s) = s%[B{G(s) see Harris (1963, p. 32). This leads to E[sR’“’ / I?(@) < co] = [G(s),/q(@]” and

-

111;

29

STOCHASTIC EPIDEMIC

which is taken as +co if the denominator is zero. For large communities we then find the expected size of minor epidemics to lie between the bounds given by

l/(1 - $L-m) < EL-GI -T < m + 4 < l/(1 - $n), when v/3, < I, and between the bounds given by

when I&,-,, > 1. Here we have ignored terms exp(--nc(m, - 1)2/2an2}]. Assuming that Pn N&+,, see that the mean size of a minor epidemic is roughly

6. PREVENTING

MAJOR

of the order O[~l/~a;’ for large communities we

EPIDEMICS

The quantity $3, , which we shall call the initial infection rate of the community, is seen to play a major role in determining the outcome of epidemics. According to the threshold theorem, it is possible to prevent major epidemics by introducing public health measures which are adequate to make VP, < 1. The initial infection rate v& is also involved in the upper bound on the mean of the total number of cases, when $7, < 1. The limited amount of data available will generally force us to make more specific assumptions in our model and we assume from here on that

t% = k%?(W)? w

= 0, 1) 2 )..., n,

where g is a known nondecreasing function with g(0) = 0. When g is the identity function, that is g(w) = w, we say that the community is homogeneously mixi%. It is possible to make the initial infection rate less than unity in several ways: (i) measures leading decrease in v,

to earlier

(ii) measures which decrease individuals will decrease /3, while (iii) immunizing g(n - no).

a fraction

discovery

of infectives

the rate of “close”

19 of the susceptibles

will

will

contacts

effect a between

decrease g(n) to

30

NIELS

G. BECKER

In order to determine how effective these measures must be in order to prevent major epidemics, it is necessary to know the value of v/3. It is then, for example, possible to determine what fraction 6 of the susceptibles must be immunized in order to prevent major epidemics by computing the smallest value of 0 such that g(n - nt9) < I/$. A ccordingly, we now consider the estimation of VP, or equivalently the estimation of the initial infection rate @g(n).

7. ESTIMATING

THE INITIAL

INFECTION

RATE

Estimating the parameters of communicable disease models presents difficulties, because only a limited part of the epidemic process is observable. The quantitative data on an epidemic often consists only of the total size and the duration of the epidemic. Occasionally, the times when cases are detected are also available. However, such data can be adequate for estimating the initial infection rate, even when these data are available only for a single outbreak of the disease. This estimating problem has been considered previously for a particular form of the above epidemic model; see Bailey (1975) for details. The approach is, in general, to use the method of maximum likelihood for the Markovian model obtained by assuming that there is no latent period and that the duration T of the infectious period has a negative exponential distribution. In order to avoid these unrealistic simplifying assumptions and the considerable amount of computation involved in these methods, Becker (1976) uses a modified embedded Galton-Watson process as an approximate description of the spread of a disease and considers estimation for this approximating model. The present approach is a more direct alternative, based on martingales associated with the epidemic model described in Sections 2 and 3. This approach is suggested by the study of Aalen (1975). In order to allow some comparisons along the way, we first consider the maximum likelihood estimation under the assumption that the process ( W, , X, , Yt , Z,} is completely observable. Let F denote the joint distribution function of (S, T), the durations of the latent and infectious periods. Split the time interval [0, T], r as given by (3) into the 32, - a intervals between the “jumps” of { W, , X, , Yt , Z,> and consider each of these intervals separately. Now recall that for a Poisson process with time-dependent intensity p, the probability that there are no jumps in the interval [0, t] is exp(-Ji pzl du), while the wait until the first jump has probability density given by pt exp(-ji pu du). Using these facts for each of the 3.2, - a subintervals of [0,7] leads to the likelihood function

where Yj* infection.

denotes the number

of infectious

infectives

just prior

to the jth

STOCHASTIC

31

EPIDEMIC

It follows that the maximum likelihood estimator of /3 is given by

(4) By specifying F, the maximum likelihood estimates of p = E[S], Y = E[T], and other parameters of interest can be found. In particular, suppose that S and T are independent and let T have an exponential distribution. Then the maximum likelihood estimator of Yis given by

which is just the average duration of the Z, infectious periods. The Fisher information associated with ,6 and v are E[Z, - a]/,P and v2E[Z,], respectively. This suggests that for large epidemics, that is when si Yt dt and si g( IV,) Yt dt are large, we use and

(a [

Y, dt)-’

to estimate the variance of the asymptotic distribution of fl and P, respectively. Actually, in the present application, the usual asymptotic theory of maximum likelihood estimation would seem to need modifications, similar to those arising for the birth-death process and counting processes in general; see Keiding (1975a) and Aalen (1975). More specifically, a random norming is required to get asymptotic normality, but this does not seem to effect any changes to the usual procedures in applications with data. Using the large-sample formula var($) N G2var@) + 2i$ cov(6, fi) + p var($) with cov($, 8) = 0 gives Q$)2/(ZT - 4 + B”/P2Z) as an estimator for var(@). These estimators are only of practical importance when the process is completely observable. This is the case, for example, with diseases for which the durations of the latent and infectious periods may be taken as known constants. Generally, however, the process {IV, , Y,} must be considered unobservable and we adopt the approach of Aalen (1975) based on martingales. We retain the assumption that T is exponentially distributed and so our model differs from that used by Bailey (1975) only in that we allow the duration of the latent period to have an arbitrary distribution. We also let g be an arbitrary, but known, nondecreasing function with g(0) = 0, rather than taking the identity function for g. 653/11/1-3

32

NIELS

G. BECKER

Under the assumptions of the model {n - W,} is a compound regular point process with intensity ~g(W,) Y, , while (2,) is a compound regular point process with intensity YJw. We now take W, and 2, to be right-continuous, without loss of generality. However, without indicating this explicitly by a change of notation, we consider that every integrand occurring from here on is modified, if necessary, to make it left-continuous and have right-hand limits. It is known (see, for example, Segall and Kailath (1975)), that M, = n - W, -B

jot g(Wu)Yu du

and Nt = vZ, -

t Y, du I0

are zero mean martingales. By equating each of these two martingales to zero, setting t = 7, and solving the resulting estimating equations, we obtain the maximum likelihood estimators for /3 and Ygiven by (4) and (5). In order to arrive at estimators that are more useful when (W, , Yt} is not observed, we construct martingales from n/r, and Nt which do not involve Yt . The theory of stochastic integrals with respect to martingales, as introduced by Kunita and Watanabe (1967), implies that the stochastic integral

Bjot gW’J dlv,

(6)

is a martingale. Since the community size is finite and g(W,) < g(n) < co, it is clear that E J6 j g( W,)[ 1dN% ( < co, where the integral is a Stieltjes integral. It follows from Proposition 3 of DolCans-Dade and Meyer (1970) that the martingale (6) can be written 4 jt g(Wu) dzu - B jt gWu)Yu dus 0

0

(7)

where the integrals are Stieltjes integrals. By subtracting the martingale (7) from Mt , it is clear that n-

W,-v/3

tg(Wu)dZu s0

(8)

is a zero mean martingale. By equating this martingale to zero and setting t = 7, the estimator

STOCHASTIC

EPIDEMIC

33

of V/I emerges, where Wit denotes the number of susceptibles at the time just prior to thejth removalj = 1, 2,..., Z, . Note that the expression for ~7 still depends on {W,}, which is generally not observed. However, the behavior of iwd, given Z, , is considerably restricted and it is possibIe to obtain bounds for the Wj+ based on Z, . Each removal after the ath removal is known to decrease the size of the susceptible population by one. It is therefore clear that, with probability one, n+a-Z,=W, 0. For a homogeneously munity $ is seen to be similar to an estimator suggested in Startsev (1970), for the epidemic model with no latent period. from the abstract how Startsev arrived at his estimator. Note that of the form

(1 - 4qg-‘,

This estimator mixing coman abstract by It is not clear $ is essentially

(15)

where g-1 denotes the average of the reciprocals of the values of g( W,) at the times when the infections occur. The difference between the estimators $ and V$ becomes more transparent when (15) is compared with (11) and (12).

8. NUMERICAL

ILLUSTRATION

Bailey and Thomas (1971) quote data of D.M. Thompson and W.H. Foege on an outbreak of smallpox in a closed community in southeastern Nigeria. The data are summarized by 1z = 119, a = 1, 2, = 30. The bounds on ~7 specified by (9) lead to 1.11 < &z < 1.28. These bounds are sufficiently close to each other to be useful. The corresponding bounds as obtained from (lo), which were obtained by a completely different approach, give I .lO and 1.26. The estimator $ti gives the estimate 1.10, under the assumption of homogeneous mixing. The approach of Bailey and Thomas (1971) gives the estimate 1.15 for v@n; see Bailey (1975, p. 125).

REFERENCES AALEN, 0. 0. 1975. “Statistical Inference for a Family of Counting Processes,” Ph.D. Thesis, University of California, Berkeley. BAILEY, N. T. J. 1975. “The Mathematical Theory of Infectious Diseases and Its Application,” Griffin, London. BAILEY, N. T. J. AND THOMAS, A. S. 1971. The estimation of parameters from population data on the general stochastic epidemic, Theor. Popul. Bid. 2, 253-270. BARTOSZYNSKI, R. 1967. Branching processes and the theory of epidemics, Proc. 5th Berkeley Symp. 4, 259-269. BECKER, N. G. 1972. Vaccination programmes for rare infectious diseases, Biometrika 59, 443-453.

36

NIELS

G. BECKER

BECKER, N. G. 1974. On parametric estimation for mortal branching processes, BiometriRa 61, 393-399. BECKER, N. G. 1975. The use of mathematical models in determining vaccination policies, Bull. Int. Statist. Inst. 46. BECKER, N. G. 1976. Estimation for an epidemic model, Biometrics 32. BHARUCIIA-REID, A. T. 1956. On the stochastic theory of epidemics, Proc. 3~d Berkeley symp. 4, 1 I l-l 19. DOL%ANS-DAD& C. AND MEYER, P. A. 1970. Integrales stochastiques par rapport aux martingales locales, Lecture Notes in Mathematics, Vol. 124, pp. 77-107, SpringerVerlag, Berlin. DOOB, J. L. 1953. “Stochastic Processes,” Wiley, New York. DWASS, M. 1969. The total progeny in a branching process and a related random walk, J. Appl. Probability 6, 682-686. HARRIS, T. E. 1963. “The Theory of Branching Processes,” Springer-Verlag, Berlin. JAGERS,P. 1975. “Branching Processes with Biological Applications,” Wiley, New York. KEIDING, N. 1975a. Maximum likelihood estimation in the birth-and-death process, Ann. Statist. 3, 363-372. KEIDING, N. 1975b. Estimation theory for branching processes, Bull. Int. Statist. Inst. 46. KERMACK, W. 0. AND MCKENDRICK, A. G. 1927. A contribution to the mathematical theory of epidemics: 1, PYOC.Roy. Sot. Ser. A 115, 700-721. KCNITA, H. AND WATANABE, S. 1967. On square integrable martingales, Nagoya Math. J. 30, 209-245. LLOYD, E. H. 1963. The epochs of emptiness of a semi-infinite discrete reservoir, J. Roy. Statist. Sot. Ser. B 25, 131-136. MOTT, J. L. 1963. The distribution of the time-to-emptiness of a discrete dam under steady demand, J. Roy. Statist. Sot. Ser. B 25, 137-139. NEYMAN, J. AND SCOTT, E. L. 1964. lit “Stochastic Models in Medicine and Biology” (J. Gurland, Ed.), pp. 45-55, University of Wisconsin Press, Madison. PATIL, G. P. 1963. Minimum variance umbiased estimation and certain problems of additive number theory, Ann. Math. Statist. 34, 1050-1056. RUBIN, I. 1972. Regular point processes and their detection, IEEE Trans. Information Theory IT-18, 547-557. SEGALL, A. AND KAILATH, T. 1975. The modeling of randomly modulated jump processes, IEEE Trans. Information Theory IT-21 (2), 135-143. STARTSEV, A. N. 1970. On estimation of the regulating parameter in the stochastic model for epidemics, (In Russia) Kratk. Nauch. Soobstch. No. 6. Izv. Akad. Nauk Uzb. SSR 1970. WHITTLE, P. 1955. The outcome of a stochastic epidemic-a note on Bailey’s paper, Biometrika 42, I 16-l 22. WORLD HEALTH ORGANIZATION 1970. Smallpox surveillance, Weekly Epidemiological Record 45, 452-468.