- Email: [email protected]
Reconstruction of measles dynamics in a vaccinated population
Vaccine 21 (2003) 2643–2650
Reconstruction of measles dynamics in a vaccinated population J. Wallinga∗ , P. Teunis, M. Kretzschmar National Institute for Public Health and the Environment, P.O. Box 1, 3720 BA, Bilthoven, The Netherlands Received 25 June 2002; received in revised form 6 December 2002; accepted 23 December 2002
Abstract The evaluation of measles vaccination programmes has been problematic because the change in actual numbers of infections and susceptibles over time cannot be directly observed. In this paper, we present a method for estimating the time series of number of susceptibles and infections, as well as the critical vaccination coverage in a vaccinated population. The proposed method is applied to data on measles outbreaks in The Netherlands. We show that the results are self-consistent and in line with available independent estimates. A potential application of the proposed method lies in detecting the loss of herd immunity and assessing the risk of major outbreaks in vaccinated populations. © 2003 Elsevier Science Ltd. All rights reserved. Keywords: Host–pathogen dynamics; Time series analysis; Basic reproduction ratio
1. Introduction Though measuring the effectiveness of measles vaccination of individuals is a relatively straightforward affair, measuring the effectiveness of measles vaccination programmes for entire populations is fraught with difficulties. An obvious indicator of a programme’s effectiveness is its ability to reduce the number of measles infections. However, the number of measles infections is usually incompletely observed in the sense that not all measles cases are presented at a general practitioner and not all general practitioners notify the presented cases. Other indicators of a programme’s effectiveness are the vaccination coverage and the resulting proportion of susceptibles among the population. However, the temporal changes in proportion of susceptibles are incompletely observed, in the sense that there is no direct clue on the susceptibility in the population before or after collecting serological data. Therefore, the incomplete observation of dynamics of measles infections and susceptibles presents a serious problem to measuring the effectiveness of measles vaccination programmes. The problem of incomplete observation of dynamics of measles infections can be overcome in part by simply assuming that the notification ratio remains constant over time. Fine and Clarkson [1] applied this assumption in a thorough study of measles case-notification data for England and Wales. They showed that before mass vaccination there was a ∗
Corresponding author. Tel.: +31-30-274-2553; fax: +31-30-274-4409. E-mail address: [email protected] (J. Wallinga).
predominant biannual pattern of major and minor epidemics. After introduction of mass vaccination there were less regular fluctuations. Building on Fine and Clarkson’s work, Finkenstädt and Grenfell [2] showed that it is possible to estimate time-varying notification ratios from case notifications, but they could not validate their estimates against independent observations of notification ratios. Recently, Bjørnstad et al. [3] showed that it is possible to estimate the mean proportion of susceptibles from a set of case-notification data for various populations, but they could not validate their estimates against independent measurements of the proportion of susceptibles. The problem of incomplete observation of temporal changes in proportion of susceptibles can be overcome in part by simply assuming that the proportion of susceptibles does not change over time, that is, by assuming that the observed proportion of susceptibles reflects an equilibrium situation. This assumption allows for sophisticated statistical analysis of serological data, and estimation of variables that are of clear epidemiological interest, such as the basic reproduction ratio [4]. This basic reproduction ratio determines the value of the critical vaccination coverage—the coverage required for a vaccination programme to achieve herd immunity and eliminate the infectious agent [5]. The value of the critical vaccination coverage for measles is estimated to be near 96% [5]. This estimate is based on data collected in England and Wales before introduction of mass vaccination and relies on the equilibrium assumption. The equilibrium assumption is clearly at odds with the biannual pattern of epidemics that is inferred from case notification
0264-410X/03/$ – see front matter © 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0264-410X(03)00051-3
2644
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
data [1–3]. There is currently no consistent way to overcome the problems of incomplete observation of both dynamics of measles infections and dynamics of susceptibles. To underline the implications of these problems of incomplete observation in the context of measuring the effectiveness of measles vaccination programmes, we consider a large measles outbreak in The Netherlands in 1999–2000 [6]. This outbreak resulted in the largest number of notified cases since introduction of mass vaccination against measles, despite the high vaccination coverage of around 95% during the two decades preceding the outbreak. The question arises how effective the measles vaccination programme is. To what extent can the large number of notifications be attributed to a higher rate of reporting measles infections? And to what extent can it be attributed to a high number of susceptibles accumulated before the outbreak? It is impossible to provide a consistent answer to these questions by looking separately at the available case-notification data and serological survey data. This paper presents a method for estimating the time series of susceptibles and measles infections—that is, a method to reconstruct the measles dynamics—by integrating information from case-notification data, serological survey data and vaccination coverage data. This reconstruction also yields various other estimates of immediate epidemiological interest, such as notification ratios and critical vaccination coverage. Moreover, it allows the identification of the periods when a vaccination programme fails to provide herd immunity and the population becomes susceptible to a major outbreak of measles. The practical usefulness of the proposed method for evaluating a vaccination programme’s effectiveness is illustrated by an application to data on measles outbreaks in The Netherlands. 2. Material and methods 2.1. Data sources Measles case notifications for The Netherlands are available from 1976 until 2001 in four-weekly intervals according to date of notification (source: Medical Inspectorate of Health, Fig. 1a). Serological data on prevalence of measles-specific antibodies are available for the years 1987 and 1996. The 1987 data comprise an age-stratified sample of 2041 individuals from the town of Nijmegen, of which a proportion of 0.053 (S.E. 0.015) tested negative. The 1996 data concern an age-stratified sample of 8303 individuals from the Dutch population, of which a proportion of 0.043 (S.E. 0.002) tested negative [6]. Annual birth rates are divided into 4-week intervals, assuming a constant birth rate within the year (source: Statistics Netherlands, Fig. 1b). Population sizes at the start of the year are available from Statistics Netherlands (Fig. 1c). Mass vaccination against measles was started in 1976. From 1976 until 1987, a plain live attenuated measles vaccine was administered to
Fig. 1. Data used for reconstructing the time series of measles in The Netherlands: (a) four weekly measles notifications from 1976 to 2001; (b) yearly births from 1975 to 2000; (c) population size from 1976 to 2001; (d) vaccine uptake from 1976 to 2001. The required serological data is indicated in Fig. 2b.
children at the age of 14 months. Since 1987 a combined measles, mumps, and rubella vaccine has been administered to children at the ages of 14 months and 9 years. Vaccine uptake rates for each age cohort are divided into 4-week intervals, assuming a constant uptake rate within
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
each age cohort (source: Medical Inspectorate of Health, Fig. 1d). 2.2. A balance equation for susceptibles
2645
B, case reports C, and vaccine uptake V , once we know the initial number of susceptibles S0 and the notification ratios ηt . 2.3. Reconstruction
The fact that susceptibles are recruited by birth and depleted by infections or vaccinations suggests a balance equation for susceptibles [1–3,7]: St = St−1 + Bt−m − ηt Ct − Vt + et ,
(1)
where the index t = 0, 1, 2, . . . refers to time. Because the measles case notifications are available in 4-week intervals, the time step in the model is 4 weeks. St is the number of susceptibles at time t. Bt−m is the number of children born in the earlier interval t − m − 1, t − m] that become susceptible over the time interval t − 1, t], m is the duration of maternally derived immunity, which is estimated to be about five 4-week intervals for measles [8]. The term ηt Ct denotes the number of individuals who are infected over the interval t − 1, t], where ηt is the notification ratio over interval t − 1, t] and Ct is the number of reported cases over that interval. The number of cases is fully reported if ηt equals unity, and under-reported if ηt is larger than unity. Vt is the number of individuals that are immunised by vaccination over the interval t − 1, t]. The book-keeping error et is assumed to follow a Gaussian distribution with a mean 0 and standard deviation σe . Note that the book-keeping error et only serves to close the balance of susceptibles; it accounts only for the occasional imbalance between immigration and emigration or death of susceptibles and eventual errors in reported births and vaccine uptake. The number of individuals who are successfully vaccinated over the 4-week interval t − 1, t] is related to the number of births at an earlier time Vt = wut Bt−v + (1 − w)ut+v −v Bt−v .
(2)
Here, v and v are the ages at which the first and second dose of vaccine are administered. The ages are measured in 4-week intervals. The first term in Eq. (2) is the number of susceptibles that are immunised by a first dose of vaccine, w is the probability that vaccination results in effective protection, which is estimated to be 0.97 [8], ut is the uptake of the first dose of vaccine administered to individuals of age v at time t. The second term in Eq. (2) denotes the number of susceptibles that are immunised by a second dose of vaccine. This number is approximated by the number of individuals of age v at time t that did receive the first dose (at time t + v − v , when they were at age v ) but that did not develop effective protection. Iteration of Eq. (1) for an initial number of susceptibles S0 yields a cumulative balance equation for susceptibles: S t = S0 +
t i=1
Bi−m −
t i=1
ηi C i −
t i=1
Vi +
t
ei .
(3)
i=1
This equation shows that the numbers of susceptibles can be reconstructed from the available data on reported births
The first unknown in the cumulative balance equation for susceptibles (Eq. (3)) is the initial number of susceptibles S0 . Its value can be estimated, together with the standard deviation of the book-keeping errors σe , from serological survey data according to the method of maximum likelihood. A condition for this maximum likelihood estimation is that we know the values of the notification ratios ηt . Appendix A shows how the likelihood function is derived from the balance equation for susceptibles (Eq. (1)). With the maximum likelihood estimate for σe we can check whether the book-keeping errors are small enough to make use of the balance equation. With the maximum likelihood estimate for S0 and Eq. (3), we can reconstruct the time series of the absolute numbers of susceptibles St . Dividing the number of susceptibles St by population size Nt gives the proportion of susceptibles. Averaging over the whole time series gives the average proportion of susceptibles, which we denote by E(S/N) (see Appendix A). The remaining unknowns in the balance equation are the notification ratios ηt . If we know the average proportion of susceptibles, and if we assume this proportion is stationary over each epidemic cycle, we know what the expected number of infections over the epidemic cycle must be (see Appendix B). The notification ratio at time t is then estimated as the ratio of the expected number of infections and the number of notified cases over the epidemic cycle. Appendix B shows how this estimator of notification ratio is derived from the balance equation for susceptibles (Eq. (1)). We take as an initial estimate for the notification ratios ηt = 1 and apply the following steps: 1. Estimate the initial number of susceptibles from the serological data, given the estimated notification ratios (using Eq. (A.4)). 2. Calculate the corresponding average proportion of susceptibles over the whole time series (using Eq. (A.5)). 3. Estimate the notification ratios, given the calculated average proportion of susceptibles (using Eq. (B.4)). We iterate these steps until the estimates have converged. The resulting estimates are substituted into Eq. (3) to reconstruct the time series of susceptibles St . The estimated notification ratios are multiplied with the observed number of cases per 4-week period to reconstruct the time series of infections ηt Ct . The convergence and accuracy of the estimation algorithm are tested by a parametric bootstrap method. The bootstrap method uses for each time interval t−1, t] a binomial distribution of notified cases with parameters ηt Ct and 1/ηt . Successive sampling from these binomial distributions yields a bootstrapped time series of case notifications. Bootstrapped
2646
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
serological data are obtained by sampling from Gaussian distributions with mean and standard error as observed in serological surveys. To each of the bootstrapped data sets the same reconstruction procedure is applied as for the original data set. The resulting 95% bootstrap confidence intervals on the estimates indicate the uncertainty due to the incomplete reporting of cases and sampling error in serological data. 2.4. Analysis of reconstructed time series Spectral analysis is used to identify periodic cycles in the reconstructed time series of measles infections ηt Ct . The power spectral density is calculated using the Lomb algorithm [9]. Here, we report results from analysing the complete time series of measles infections. If only a very small proportion of all measles infections is notified, there will be many short periods with successive zero case notifications, which might alter the power spectral densities. To check for such an artefact, we also carried out spectral analysis for a time series where 4-week periods with zero case notifications were excluded. The results for the non-zero time series were found to be almost identical to the results for the complete time series. 2.5. Estimation of the threshold level for susceptibles and critical vaccination coverage The basic reproduction ratio R0 is defined as the expected number of secondary cases produced by an infectious individual in a completely susceptible population [5,10–12]. For a homogeneously mixing population that consists of susceptibles and immunes, the expected number of cases produced per infectious individual at time t is R0 (St /Nt ), where St is the number of susceptibles and Nt is the population size. Over a long period of time, the geometric mean of the numbers of secondary cases produced per infectious individual must be one if the number of infections remains larger than zero but bounded. In mathematical terms this reads R0 (St /Nt ) = 1, where · represents the geometric mean. Rewriting this expression shows that the basic reproduction ratio is estimated by the geometric mean of the inverse proportion of susceptibles [12]: Nt R0 = . (4) St The threshold level for susceptibles sc is defined as the proportion of susceptibles where each infectious individual produces on average one secondary case. A major outbreak is only possible if the proportion of susceptibles exceeds the threshold level. In a homogeneously mixing population, this level is equal to the proportion of susceptibles in an endemic equilibrium [5]: sc =
1 . R0
(5)
The critical vaccination coverage is defined for a hypothetical vaccine that is 100% effective and administered soon after birth. It is defined as the coverage that maintains the proportion of susceptibles at the threshold level [5]: pc = 1 −
1 . R0
(6)
All three equations in this section rely on the assumption that the population is mixing homogeneously. Of course, this assumption can only be true by approximation.
3. Results Fig. 2a shows the estimated proportion of cases actually notified. The estimated proportion of cases that were notified over the years 1976–1980 is on average 0.024, which is in line with an independent estimate of 0.028 obtained by comparing nation-wide notifications with measles incidence as reported in sentinel studies (van den Hof et al., 1998, RIVM report 213676 008). The increased proportion of cases notified around 1999 coincides with the set-up of a case-register and more active searching for measles cases. The estimated proportion of cases that were notified in 1999–2000 is on average 0.067, which is in line with an independent estimate of 0.09 based on comparing questionnaires that revealed 164 cases among a study population with 15 case reports among the same population during a measles outbreak (van Isterdael et al., in preparation). Therefore, the estimated time series of notification ratios agrees with all available independent information about notification ratios. The reconstructed time series of infections, as shown in Fig. 2b, reveals a slightly different picture than the time series of case notifications (compare Fig. 1a) in that the largest outbreak occurred just after the start of vaccination in 1976–1977. Other major outbreaks occurred in 1987–1988 and 1999–2000, and in between these major outbreaks were episodes with smaller outbreaks in 1980–1984 and 1992–1996. Time series analysis reveals a significant periodicity at a frequency of 0.17 per year, corresponding to an interepidemic period of 5.9 years (the null-hypothesis that this periodicity is produced by uncorrelated white noise is rejected with P 0.05). As predicted by epidemic theory [13,14], the estimated interepidemic period is much larger for a vaccinated population than for unvaccinated populations. The book-keeping errors in the balance equation for susceptibles are estimated to have an extremely small standard deviation of σe = 3 × 10−13 . Therefore, we can almost exactly reconstruct the time series of susceptibles from births, vaccinations, and the estimated number of infections. The estimated number of susceptibles at onset of the mass vaccination programme is estimated to be 704 × 103 (95% bootstrap confidence interval (634 − 803) × 103 ); this number corresponds to a proportion 0.05 of the entire population being susceptible at onset of the mass vaccination
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
2647
tervals is of the order of 0.01 and much of this uncertainty is due to sampling error in serological survey data. The basic reproduction ratio R0 is estimated as 23.01 (95% bootstrap confidence interval 20.81–25.33). The corresponding critical vaccination coverage pc is estimated as 0.957 (95% bootstrap confidence interval 0.952–0.960). As far as we are aware of, this is the first estimate of measles critical vaccination coverage based on data obtained from a population after introduction of mass vaccination. The corresponding threshold level of susceptibles sc is estimated as 0.043 (95% bootstrap confidence interval 0.037–0.049). Fig. 2d shows the difference between the proportion of susceptibles and the threshold level. The excess proportion of susceptibles can be estimated accurately; the 95% bootstrap confidence limits almost overlap the estimate. Fig. 2d reveals the episodes during which the proportion of susceptibles exceeds the threshold level, that is, when the vaccination programme fails to provide herd immunity and the vaccinated population is susceptible to major measles outbreaks. Comparison with the reconstructed times series of infections (Fig. 2b) confirms that major outbreaks indeed occur only when the proportion of susceptibles has reached or exceeded the threshold level.
4. Discussion
Fig. 2. Reconstructed dynamics of measles in The Netherlands after introduction of mass vaccination: (a) proportion of cases actually notified; (b) number of infections in 4-weekly intervals; (c) proportion of susceptibles among the population, and the observed proportions of susceptibles with 95% confidence intervals (markers); (d) proportion of susceptibles in excess of the threshold level for major outbreaks. Dashed lines indicate the 95% bootstrap confidence limits.
programme. Fig. 2c shows that the reconstructed proportion of susceptibles is stationary around a mean value of 0.043 and that temporal fluctuations remain within a limited range of 0.04–0.05. The width of the 95% bootstrap confidence in-
This study reveals the values of hitherto inaccessible variables that characterize measles in a vaccinated population, such as proportion of susceptibles, number of measles infections, notification ratio, and critical vaccination coverage. The values are inferred from combined information on vaccination coverage, serological survey data and case notifications. The framework for combining information from the various data sources is provided by the balance equation for susceptibles (Eq. (1)). The use of such a balance equation to analyse trends in susceptibles is not entirely new; the use of an explicit balance equation dates back to 1933 [16] and the more implicit use of balance considerations dates back even further to 1906 [17]. The contribution of the present study to the rapidly expanding literature on this topic [1–3,7] is to use the balance equation to estimate the temporal changes in the absolute number of susceptibles from serological data (Appendix A) and to use the balance equation to estimate the time-varying notification ratios (Appendix B). We would like to point out that the use of the balance equation can be extended even further by keeping track of susceptibles per age-cohort [18] or per region [3]. This makes it possible to infer age-specific or region-specific notification ratios and identify the age groups or regions that are least protected. A reason for not including this refinement here is that the Dutch case-notification data is structured by age and by geographical location only from 1988 onward: the time series is too short to obtain reliable estimates of the notification ratios.
2648
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
The estimation of the time-varying notification ratios relies on the assumption that the proportion of susceptibles remains stationary over epidemic cycles. That is, we assume similarity of expected proportions of susceptibles that are in the same phase of the epidemic cycle. This is a considerable improvement over the assumption of an equilibrium. The proposed method for estimating notification ratios is at first sight more or less equivalent to the local regression approach proposed by Finkenstädt and Grenfell [2]. This approach implicitly assumes stationarity of absolute number of susceptibles over an arbitrarily chosen period. However, applying the local regression approach to a set of case notification data for subpopulations of varying sizes shows that it is the proportion of susceptibles, rather than the absolute number, that remains stationary [3]. In this sense, the assumptions underlying the local regression approach are inconsistent with its results, whereas our approach does not have such a problem. As in the local regression approach, the assumption of stationarity over epidemic cycles does not hold in all situations. For example, it does not hold when the vaccination programme has eliminated measles and an epidemic cycle is absent. The excellent agreement between the estimated notification ratios and the available independent information suggests that at least for measles in The Netherlands the assumption of stationarity holds. The estimation of the basic reproduction ratio, the critical vaccination coverage, and the threshold for susceptibles, relies on the assumption of a homogeneously mixing population. This assumption is motivated by parsimony, as homogeneous mixing is the simplest assumption that produces meaningful and useful estimates. However, heterogeneous mixing may be more realistic as it accounts for the socio-geographical clustering of unvaccinated individuals in The Netherlands [8] and for the observed distribution of immunity across ages [15]. We are currently analyzing the estimated time series of measles infections under various mixing assumptions, and the preliminary results suggest that some features of the time series are better explained by heterogeneous mixing than homogeneous mixing, although this has little effect on estimated values of the basic reproduction ratio and the critical vaccination coverage. Notwithstanding these reservations concerning the homogeneous mixing assumption, the practical usefulness of the resulting estimates should not be underestimated. The pattern of major outbreaks is consistent with the pattern of excess of susceptibles over the threshold value, that is, with the pattern of losing herd immunity according to the homogeneous mixing assumption. For practical purposes, the assumption of homogeneous mixing allows us to estimate when a vaccinated population is protected by herd immunity and when it is not. Our estimate of the basic reproduction ratio R0 for measles in a vaccinated population is 23.01. This value is higher than earlier estimates based on the average age at infection as observed before introduction of mass vaccination (range: 14–18) [5], but lower than the estimate of 29.9 based on time series of case notifications before introduction of
mass vaccination [3]. The slight discrepancies between the estimated values could result from differences between the estimation methods, differences between the observed populations, and differences between the times of observation. The fact that our estimate for a vaccinated population lies within the range of previously published values for unvaccinated populations suggests that the basic reproduction ratio for measles has not been drastically altered by vaccination. Our estimates of the proportion of measles cases that are notified are remarkably low and show fluctuations over time. Apparently, the number of cases notified by practicing physicians is a rather poor indicator of the number of actual measles infections among their patients. Therefore case-notification data alone gives a distorted picture of the impact of the measles vaccination programme on the number of measles infections. Time series reconstruction, as proposed here, is necessary to correct for the low and fluctuating notification ratios. But the statistical algorithms for time series reconstruction cannot retrieve all the information that is lost due to low reporting levels, and a higher level of reporting would certainly improve the estimates of the basic reproduction ratio and the critical vaccination coverage. An approach that has proven to increase the reporting level is to search actively for measles cases by public health services, another approach might be to increase the motivation for reporting cases by communicating the insights obtained from notification data back to the reporting physicians. The reconstruction of measles dynamics in the vaccinating Dutch population makes it possible to identify some of the causes of the large measles outbreak in 1999–2000. Firstly, the actual vaccination coverage remained near or below the critical level needed for elimination, and the proportion of susceptibles increased beyond the threshold level, such that from 1997 onwards a measles outbreak could occur. Secondly, the rate of import of measles infections is such that previous outbreaks started on average 2 years after the proportion of susceptibles increased beyond the threshold level, and therefore an outbreak could be expected around 1999, approximately 2 years after the threshold was exceeded in 1997. Thirdly, there has been a drastic increase in the proportion of infections that were notified in 1999. Together, these three causes account for the outbreak with the largest number of notified cases since introduction of vaccination. It is of course not implied that these three causes exhaust all possible relevant causes. The reconstruction of measles dynamics also reveals the epidemiology of measles in a vaccinated population in much greater detail than ever before. A striking feature is that regular epidemics are maintained, even when measles transmission is interrupted after each major outbreak. Already in 1956, Bartlett [13] proposed a mathematical explanation for regularity in such recurrent epidemic outbreaks: the epidemic period consists of a fixed duration required for replenishing the proportion of susceptibles up to the threshold level, and a variable duration for waiting until measles is reintroduced and another outbreak starts. Bartlett assumed
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
that all susceptibles are infected during a major outbreak. The reconstructed time series of proportion of susceptibles, as presented here, shows this is not true: the proportion of susceptibles varies during an epidemic cycle only within a very small range. We are currently adapting Bartlett’s mathematical argument to allow for smaller, more realistic, fluctuations in proportion of susceptibles. In conclusion, this work demonstrates that we can overcome the problems of incomplete observation of dynamics of measles infections and susceptibles by combining case-notification data and serological data. Much more information relevant to the evaluation of vaccination programmes can be extracted from the combined data sources than from each source in separation. We hope that the proposed methods of analysis will be applied to other data sets for other infectious diseases, other vaccination programmes and other populations. Reconstructing the dynamics of infectious diseases in vaccinated populations greatly increases our under− (S0 , σe |xk , sk , ηt ) =
kobs 1 k=1
2
log(2π(sk2 + tk σe2 )) +
2649
where σk is the standard deviation of proportion of susceptibles at time tk . The expected number of susceptibles at time tk follows from Eq. (3), and dividing this number by the population size Nk gives the proportion of susceptibles: tk tk tk Bi−m − i=1 ηi Ci − i=1 S0 + i=1 Vi Sk = . (A.2) Nk Nk The variance of proportion of susceptibles at time tk can be split into two factors, one factor is due to the book-keeping error and follows from Eq. (3), the other factor is the observation error sk : σk2 = sk2 + tk σe2 .
(A.3)
Substituting Eqs. (A.2) and (A.3) into Eq. (A.1) and summing over the number of serological surveys, gives the log-likelihood function for initial number of susceptibles S0 and the standard deviation of the book-keeping error σe
2 tk tk tk kobs x − S + B − /N η C − V 0 i−m i i i k k i=1 i=1 i=1 2(sk2 + tk σe2 )
k=1
standing of how vaccination influences the epidemiology of the infectious disease, it indicates when the vaccination programme fails to provide herd immunity, and thus enhances our capability to measure the actual impact of vaccination.
Acknowledgements This manuscript has profited from discussion with Bryan Grenfell at the University of Cambridge, and from comments on an earlier version by Susan van den Hof and Nico Nagelkerke at the National Institute for Public Health and the Environment.
.
(A.4) Maximizing the log-likelihood for S0 and σe gives the maximum likelihood estimates. With the maximum likelihood estimate for S0 and Eq. (3), we can reconstruct the time series of the absolute numbers of susceptibles St . The average proportion of susceptibles is then calculated as T S 1 St E , (A.5) = N T Nt t=0
where T is the length of the time series of case notifications.
Appendix B. Estimation of notification ratios Appendix A. Maximum likelihood estimation of initial number of susceptibles The initial number of susceptibles S0 and the standard deviation of the book-keeping error σe are estimated from serological data. We use the method of maximum likelihood and assume that the notification ratios ηt are given. We have data on serological surveys conducted at times tk (k = 1, . . . , kobs ; kobs denotes the number of serological surveys available). The population sizes at time of each survey are indicated as Nk . For each serological survey we have the mean proportion of susceptibles xk and standard error of the mean sk . For a Gaussian distribution of mean proportion of susceptibles, the log-likelihood function for number of susceptibles at the time of the serological survey Sk takes the following form: −(Sk , σk |xk ) =
1 (xk − Sk /Nk )2 , log(2πσk2 ) + 2 2σk2
(A.1)
The notification ratio at time t, denoted as ηt , is estimated as the ratio of expected number of infections over the number of reported cases. Both numbers are taken over the interval t−d, t+d, where d is a positive integer. We derive this ratio from the balance equation for susceptibles (Eq. (1)) by first computing St−d − St+d and then solving for ηt , assuming that the cumulative error over this interval is negligible. This yields: ηt ≈
(St−d − St+d ) +
t+d
i=t−d+1 Bi−m t+d i=t−d+1 Ci
−
t+d
i=t−d+1 Vi
.
(B.1) The length of the interval t − d, t + d is 2d. We choose values of 2d close to the length of an epidemic interval or a multiple thereof. We find such values by calculating a weight ρd that is approximately proportional to the autocorrelation function of case reports at lag 2d while
2650
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
remaining positive: ρd = q
T
Ct−d Ct+d ,
(B.2)
t=1
where T is the length of the time series of case notifications and q is a normalization constant, chosen such that all ρd sum to one. For the specific time series of case notifications studied here, the weights are highest around a length of about 6 years, which is the mean epidemic interval. We exclude lengths shorter than 4 years to ensure that the denominator in Eq. (B.1) is large enough, and lengths longer than 12 years because it is unlikely that the notification ratio remains constant over such a long interval. Thus, we consider only lengths between 4 and 12 years (that is, 26 ≤ d ≤ 78). Furthermore, we have to express the unknown term for change in number of susceptibles in Eq. (B.1), (St−d −St+d ), in terms of known variables. We rely on the approximate stationarity of proportion of susceptibles over epidemic cycles. That is, the change in number of susceptibles is approximately the change in population size times the average proportion of susceptibles in the population, E(S/N), for values of interval length 2d close to an epidemic interval or a multiple thereof: S St−d − St+d ≈ E (Nt−d − Nt+d ). (B.3) N The expected value of proportion of susceptibles E(S/N) is calculated according to Eq. (A.5). The estimator for notification ratio is obtained by substitution of Eq. (B.3) into Eq. (B.1) and taking a weighted average over the appropriate values of d: ηt ≈
d=78 d=26
ρd
t+d E(S/N)(Nt−d −Nt+d )+ t+d i=t−d+1 Bi−m − i=t−d+1 Vi . t+d i=t−d+1 Ci (B.4)
References [1] Fine PEM, Clarkson JA. Measles in England and Wales. I. An analysis of factors underlying seasonal patterns. Int J Epidemiol 1982;11:5–14.
[2] Finkenstädt BF, Grenfell BT. Time series modelling of childhood diseases: a dynamical systems approach. Appl Statist 2000;49:187– 205. [3] Bjørnstad ON, Finkenstädt BF, Grenfell BT. Dynamics of measles epidemics: estimating scaling of transmission rates using a time series SIR model. Ecol Monogr 2002;72:169–84. [4] Farrington CP, Kanaan MN, Gay NJ. Estimation of the basic reproduction number for infectious diseases from age-stratified serological survey data. Appl Statist 2001;50:251–9. [5] Anderson RM, May RM. Directly transmitted infectious diseases: control by vaccination. Science 1982;215:1053–60. [6] vanden Hof S, Berbers GAM, de Melker H, Conyn-van Spaendonck MAE. Sero-epidemiology of measles antibodies in The Netherlands: a cross-sectional study in a national sample and in communities with low vaccine coverage. Vaccine 2000;18:931–40. [7] Grenfell BT, Bjørnstad ON, Finkenstädt BF. Dynamics of measles epidemics: scaling noise, determinism, and predictability with the TSIR model. Ecol Monogr 2002;72:185–202. [8] van den Hof S, Wallinga J, Widdowson M-A, Conyn-van Spaendonck MAE. Protecting the vaccinating population in the face of a measles epidemic: assessing the impact of adjusted vaccination schedules. Epidemiol Infect 2002;128:47–57. [9] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C. Cambridge: Cambridge University Press; 1992. [10] Diekmann O, Heesterbeek JAP, Metz JAJ. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J Math Biol 1990;28:365–82. [11] Dietz K. The estimation of the basic reproduction number for infectious diseases. Stat Methods Med Res 1993;2:23–41. [12] Keeling MJ, Grenfell BT. Individual-based perspectives on R0 . J Theor Biol 2000;203:51–61. [13] Bartlett MS. Deterministic and stochastic models for recurrent epidemics. In: Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability. Berkeley: University of California Press; 1956. p. 81–109. [14] Griffiths DA. The effect of measles vaccination on the incidence of measles in the community. J R Stat Soc A 1973;136:441–9. [15] Wallinga J, Lévy-Bruhl D. Gay NJ, Wachmann CH. Estimation of measles reproduction ratios and prospects for elimination of measles by vaccination in some Western European countries. Epidemiol Infect 2001;127:281–95. [16] Hedrich AW. Monthly estimates of the child population ‘susceptible’ to measles, 1900–1931, Baltimore, Md. Am J Hyg 1933;17:613– 36. [17] Hamer WH. Epidemic disease in England. Lancet 1906;1:733– 9. [18] Fine PEM, Clarkson JA. Measles in England and Wales. II. The impact of the measles vaccination programme on the distribution of immunity in the population. Int J Epidemiol 1982;11:15– 25.
Reconstruction of measles dynamics in a vaccinated population J. Wallinga∗ , P. Teunis, M. Kretzschmar National Institute for Public Health and the Environment, P.O. Box 1, 3720 BA, Bilthoven, The Netherlands Received 25 June 2002; received in revised form 6 December 2002; accepted 23 December 2002
Abstract The evaluation of measles vaccination programmes has been problematic because the change in actual numbers of infections and susceptibles over time cannot be directly observed. In this paper, we present a method for estimating the time series of number of susceptibles and infections, as well as the critical vaccination coverage in a vaccinated population. The proposed method is applied to data on measles outbreaks in The Netherlands. We show that the results are self-consistent and in line with available independent estimates. A potential application of the proposed method lies in detecting the loss of herd immunity and assessing the risk of major outbreaks in vaccinated populations. © 2003 Elsevier Science Ltd. All rights reserved. Keywords: Host–pathogen dynamics; Time series analysis; Basic reproduction ratio
1. Introduction Though measuring the effectiveness of measles vaccination of individuals is a relatively straightforward affair, measuring the effectiveness of measles vaccination programmes for entire populations is fraught with difficulties. An obvious indicator of a programme’s effectiveness is its ability to reduce the number of measles infections. However, the number of measles infections is usually incompletely observed in the sense that not all measles cases are presented at a general practitioner and not all general practitioners notify the presented cases. Other indicators of a programme’s effectiveness are the vaccination coverage and the resulting proportion of susceptibles among the population. However, the temporal changes in proportion of susceptibles are incompletely observed, in the sense that there is no direct clue on the susceptibility in the population before or after collecting serological data. Therefore, the incomplete observation of dynamics of measles infections and susceptibles presents a serious problem to measuring the effectiveness of measles vaccination programmes. The problem of incomplete observation of dynamics of measles infections can be overcome in part by simply assuming that the notification ratio remains constant over time. Fine and Clarkson [1] applied this assumption in a thorough study of measles case-notification data for England and Wales. They showed that before mass vaccination there was a ∗
Corresponding author. Tel.: +31-30-274-2553; fax: +31-30-274-4409. E-mail address: [email protected] (J. Wallinga).
predominant biannual pattern of major and minor epidemics. After introduction of mass vaccination there were less regular fluctuations. Building on Fine and Clarkson’s work, Finkenstädt and Grenfell [2] showed that it is possible to estimate time-varying notification ratios from case notifications, but they could not validate their estimates against independent observations of notification ratios. Recently, Bjørnstad et al. [3] showed that it is possible to estimate the mean proportion of susceptibles from a set of case-notification data for various populations, but they could not validate their estimates against independent measurements of the proportion of susceptibles. The problem of incomplete observation of temporal changes in proportion of susceptibles can be overcome in part by simply assuming that the proportion of susceptibles does not change over time, that is, by assuming that the observed proportion of susceptibles reflects an equilibrium situation. This assumption allows for sophisticated statistical analysis of serological data, and estimation of variables that are of clear epidemiological interest, such as the basic reproduction ratio [4]. This basic reproduction ratio determines the value of the critical vaccination coverage—the coverage required for a vaccination programme to achieve herd immunity and eliminate the infectious agent [5]. The value of the critical vaccination coverage for measles is estimated to be near 96% [5]. This estimate is based on data collected in England and Wales before introduction of mass vaccination and relies on the equilibrium assumption. The equilibrium assumption is clearly at odds with the biannual pattern of epidemics that is inferred from case notification
0264-410X/03/$ – see front matter © 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0264-410X(03)00051-3
2644
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
data [1–3]. There is currently no consistent way to overcome the problems of incomplete observation of both dynamics of measles infections and dynamics of susceptibles. To underline the implications of these problems of incomplete observation in the context of measuring the effectiveness of measles vaccination programmes, we consider a large measles outbreak in The Netherlands in 1999–2000 [6]. This outbreak resulted in the largest number of notified cases since introduction of mass vaccination against measles, despite the high vaccination coverage of around 95% during the two decades preceding the outbreak. The question arises how effective the measles vaccination programme is. To what extent can the large number of notifications be attributed to a higher rate of reporting measles infections? And to what extent can it be attributed to a high number of susceptibles accumulated before the outbreak? It is impossible to provide a consistent answer to these questions by looking separately at the available case-notification data and serological survey data. This paper presents a method for estimating the time series of susceptibles and measles infections—that is, a method to reconstruct the measles dynamics—by integrating information from case-notification data, serological survey data and vaccination coverage data. This reconstruction also yields various other estimates of immediate epidemiological interest, such as notification ratios and critical vaccination coverage. Moreover, it allows the identification of the periods when a vaccination programme fails to provide herd immunity and the population becomes susceptible to a major outbreak of measles. The practical usefulness of the proposed method for evaluating a vaccination programme’s effectiveness is illustrated by an application to data on measles outbreaks in The Netherlands. 2. Material and methods 2.1. Data sources Measles case notifications for The Netherlands are available from 1976 until 2001 in four-weekly intervals according to date of notification (source: Medical Inspectorate of Health, Fig. 1a). Serological data on prevalence of measles-specific antibodies are available for the years 1987 and 1996. The 1987 data comprise an age-stratified sample of 2041 individuals from the town of Nijmegen, of which a proportion of 0.053 (S.E. 0.015) tested negative. The 1996 data concern an age-stratified sample of 8303 individuals from the Dutch population, of which a proportion of 0.043 (S.E. 0.002) tested negative [6]. Annual birth rates are divided into 4-week intervals, assuming a constant birth rate within the year (source: Statistics Netherlands, Fig. 1b). Population sizes at the start of the year are available from Statistics Netherlands (Fig. 1c). Mass vaccination against measles was started in 1976. From 1976 until 1987, a plain live attenuated measles vaccine was administered to
Fig. 1. Data used for reconstructing the time series of measles in The Netherlands: (a) four weekly measles notifications from 1976 to 2001; (b) yearly births from 1975 to 2000; (c) population size from 1976 to 2001; (d) vaccine uptake from 1976 to 2001. The required serological data is indicated in Fig. 2b.
children at the age of 14 months. Since 1987 a combined measles, mumps, and rubella vaccine has been administered to children at the ages of 14 months and 9 years. Vaccine uptake rates for each age cohort are divided into 4-week intervals, assuming a constant uptake rate within
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
each age cohort (source: Medical Inspectorate of Health, Fig. 1d). 2.2. A balance equation for susceptibles
2645
B, case reports C, and vaccine uptake V , once we know the initial number of susceptibles S0 and the notification ratios ηt . 2.3. Reconstruction
The fact that susceptibles are recruited by birth and depleted by infections or vaccinations suggests a balance equation for susceptibles [1–3,7]: St = St−1 + Bt−m − ηt Ct − Vt + et ,
(1)
where the index t = 0, 1, 2, . . . refers to time. Because the measles case notifications are available in 4-week intervals, the time step in the model is 4 weeks. St is the number of susceptibles at time t. Bt−m is the number of children born in the earlier interval t − m − 1, t − m] that become susceptible over the time interval t − 1, t], m is the duration of maternally derived immunity, which is estimated to be about five 4-week intervals for measles [8]. The term ηt Ct denotes the number of individuals who are infected over the interval t − 1, t], where ηt is the notification ratio over interval t − 1, t] and Ct is the number of reported cases over that interval. The number of cases is fully reported if ηt equals unity, and under-reported if ηt is larger than unity. Vt is the number of individuals that are immunised by vaccination over the interval t − 1, t]. The book-keeping error et is assumed to follow a Gaussian distribution with a mean 0 and standard deviation σe . Note that the book-keeping error et only serves to close the balance of susceptibles; it accounts only for the occasional imbalance between immigration and emigration or death of susceptibles and eventual errors in reported births and vaccine uptake. The number of individuals who are successfully vaccinated over the 4-week interval t − 1, t] is related to the number of births at an earlier time Vt = wut Bt−v + (1 − w)ut+v −v Bt−v .
(2)
Here, v and v are the ages at which the first and second dose of vaccine are administered. The ages are measured in 4-week intervals. The first term in Eq. (2) is the number of susceptibles that are immunised by a first dose of vaccine, w is the probability that vaccination results in effective protection, which is estimated to be 0.97 [8], ut is the uptake of the first dose of vaccine administered to individuals of age v at time t. The second term in Eq. (2) denotes the number of susceptibles that are immunised by a second dose of vaccine. This number is approximated by the number of individuals of age v at time t that did receive the first dose (at time t + v − v , when they were at age v ) but that did not develop effective protection. Iteration of Eq. (1) for an initial number of susceptibles S0 yields a cumulative balance equation for susceptibles: S t = S0 +
t i=1
Bi−m −
t i=1
ηi C i −
t i=1
Vi +
t
ei .
(3)
i=1
This equation shows that the numbers of susceptibles can be reconstructed from the available data on reported births
The first unknown in the cumulative balance equation for susceptibles (Eq. (3)) is the initial number of susceptibles S0 . Its value can be estimated, together with the standard deviation of the book-keeping errors σe , from serological survey data according to the method of maximum likelihood. A condition for this maximum likelihood estimation is that we know the values of the notification ratios ηt . Appendix A shows how the likelihood function is derived from the balance equation for susceptibles (Eq. (1)). With the maximum likelihood estimate for σe we can check whether the book-keeping errors are small enough to make use of the balance equation. With the maximum likelihood estimate for S0 and Eq. (3), we can reconstruct the time series of the absolute numbers of susceptibles St . Dividing the number of susceptibles St by population size Nt gives the proportion of susceptibles. Averaging over the whole time series gives the average proportion of susceptibles, which we denote by E(S/N) (see Appendix A). The remaining unknowns in the balance equation are the notification ratios ηt . If we know the average proportion of susceptibles, and if we assume this proportion is stationary over each epidemic cycle, we know what the expected number of infections over the epidemic cycle must be (see Appendix B). The notification ratio at time t is then estimated as the ratio of the expected number of infections and the number of notified cases over the epidemic cycle. Appendix B shows how this estimator of notification ratio is derived from the balance equation for susceptibles (Eq. (1)). We take as an initial estimate for the notification ratios ηt = 1 and apply the following steps: 1. Estimate the initial number of susceptibles from the serological data, given the estimated notification ratios (using Eq. (A.4)). 2. Calculate the corresponding average proportion of susceptibles over the whole time series (using Eq. (A.5)). 3. Estimate the notification ratios, given the calculated average proportion of susceptibles (using Eq. (B.4)). We iterate these steps until the estimates have converged. The resulting estimates are substituted into Eq. (3) to reconstruct the time series of susceptibles St . The estimated notification ratios are multiplied with the observed number of cases per 4-week period to reconstruct the time series of infections ηt Ct . The convergence and accuracy of the estimation algorithm are tested by a parametric bootstrap method. The bootstrap method uses for each time interval t−1, t] a binomial distribution of notified cases with parameters ηt Ct and 1/ηt . Successive sampling from these binomial distributions yields a bootstrapped time series of case notifications. Bootstrapped
2646
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
serological data are obtained by sampling from Gaussian distributions with mean and standard error as observed in serological surveys. To each of the bootstrapped data sets the same reconstruction procedure is applied as for the original data set. The resulting 95% bootstrap confidence intervals on the estimates indicate the uncertainty due to the incomplete reporting of cases and sampling error in serological data. 2.4. Analysis of reconstructed time series Spectral analysis is used to identify periodic cycles in the reconstructed time series of measles infections ηt Ct . The power spectral density is calculated using the Lomb algorithm [9]. Here, we report results from analysing the complete time series of measles infections. If only a very small proportion of all measles infections is notified, there will be many short periods with successive zero case notifications, which might alter the power spectral densities. To check for such an artefact, we also carried out spectral analysis for a time series where 4-week periods with zero case notifications were excluded. The results for the non-zero time series were found to be almost identical to the results for the complete time series. 2.5. Estimation of the threshold level for susceptibles and critical vaccination coverage The basic reproduction ratio R0 is defined as the expected number of secondary cases produced by an infectious individual in a completely susceptible population [5,10–12]. For a homogeneously mixing population that consists of susceptibles and immunes, the expected number of cases produced per infectious individual at time t is R0 (St /Nt ), where St is the number of susceptibles and Nt is the population size. Over a long period of time, the geometric mean of the numbers of secondary cases produced per infectious individual must be one if the number of infections remains larger than zero but bounded. In mathematical terms this reads R0 (St /Nt ) = 1, where · represents the geometric mean. Rewriting this expression shows that the basic reproduction ratio is estimated by the geometric mean of the inverse proportion of susceptibles [12]: Nt R0 = . (4) St The threshold level for susceptibles sc is defined as the proportion of susceptibles where each infectious individual produces on average one secondary case. A major outbreak is only possible if the proportion of susceptibles exceeds the threshold level. In a homogeneously mixing population, this level is equal to the proportion of susceptibles in an endemic equilibrium [5]: sc =
1 . R0
(5)
The critical vaccination coverage is defined for a hypothetical vaccine that is 100% effective and administered soon after birth. It is defined as the coverage that maintains the proportion of susceptibles at the threshold level [5]: pc = 1 −
1 . R0
(6)
All three equations in this section rely on the assumption that the population is mixing homogeneously. Of course, this assumption can only be true by approximation.
3. Results Fig. 2a shows the estimated proportion of cases actually notified. The estimated proportion of cases that were notified over the years 1976–1980 is on average 0.024, which is in line with an independent estimate of 0.028 obtained by comparing nation-wide notifications with measles incidence as reported in sentinel studies (van den Hof et al., 1998, RIVM report 213676 008). The increased proportion of cases notified around 1999 coincides with the set-up of a case-register and more active searching for measles cases. The estimated proportion of cases that were notified in 1999–2000 is on average 0.067, which is in line with an independent estimate of 0.09 based on comparing questionnaires that revealed 164 cases among a study population with 15 case reports among the same population during a measles outbreak (van Isterdael et al., in preparation). Therefore, the estimated time series of notification ratios agrees with all available independent information about notification ratios. The reconstructed time series of infections, as shown in Fig. 2b, reveals a slightly different picture than the time series of case notifications (compare Fig. 1a) in that the largest outbreak occurred just after the start of vaccination in 1976–1977. Other major outbreaks occurred in 1987–1988 and 1999–2000, and in between these major outbreaks were episodes with smaller outbreaks in 1980–1984 and 1992–1996. Time series analysis reveals a significant periodicity at a frequency of 0.17 per year, corresponding to an interepidemic period of 5.9 years (the null-hypothesis that this periodicity is produced by uncorrelated white noise is rejected with P 0.05). As predicted by epidemic theory [13,14], the estimated interepidemic period is much larger for a vaccinated population than for unvaccinated populations. The book-keeping errors in the balance equation for susceptibles are estimated to have an extremely small standard deviation of σe = 3 × 10−13 . Therefore, we can almost exactly reconstruct the time series of susceptibles from births, vaccinations, and the estimated number of infections. The estimated number of susceptibles at onset of the mass vaccination programme is estimated to be 704 × 103 (95% bootstrap confidence interval (634 − 803) × 103 ); this number corresponds to a proportion 0.05 of the entire population being susceptible at onset of the mass vaccination
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
2647
tervals is of the order of 0.01 and much of this uncertainty is due to sampling error in serological survey data. The basic reproduction ratio R0 is estimated as 23.01 (95% bootstrap confidence interval 20.81–25.33). The corresponding critical vaccination coverage pc is estimated as 0.957 (95% bootstrap confidence interval 0.952–0.960). As far as we are aware of, this is the first estimate of measles critical vaccination coverage based on data obtained from a population after introduction of mass vaccination. The corresponding threshold level of susceptibles sc is estimated as 0.043 (95% bootstrap confidence interval 0.037–0.049). Fig. 2d shows the difference between the proportion of susceptibles and the threshold level. The excess proportion of susceptibles can be estimated accurately; the 95% bootstrap confidence limits almost overlap the estimate. Fig. 2d reveals the episodes during which the proportion of susceptibles exceeds the threshold level, that is, when the vaccination programme fails to provide herd immunity and the vaccinated population is susceptible to major measles outbreaks. Comparison with the reconstructed times series of infections (Fig. 2b) confirms that major outbreaks indeed occur only when the proportion of susceptibles has reached or exceeded the threshold level.
4. Discussion
Fig. 2. Reconstructed dynamics of measles in The Netherlands after introduction of mass vaccination: (a) proportion of cases actually notified; (b) number of infections in 4-weekly intervals; (c) proportion of susceptibles among the population, and the observed proportions of susceptibles with 95% confidence intervals (markers); (d) proportion of susceptibles in excess of the threshold level for major outbreaks. Dashed lines indicate the 95% bootstrap confidence limits.
programme. Fig. 2c shows that the reconstructed proportion of susceptibles is stationary around a mean value of 0.043 and that temporal fluctuations remain within a limited range of 0.04–0.05. The width of the 95% bootstrap confidence in-
This study reveals the values of hitherto inaccessible variables that characterize measles in a vaccinated population, such as proportion of susceptibles, number of measles infections, notification ratio, and critical vaccination coverage. The values are inferred from combined information on vaccination coverage, serological survey data and case notifications. The framework for combining information from the various data sources is provided by the balance equation for susceptibles (Eq. (1)). The use of such a balance equation to analyse trends in susceptibles is not entirely new; the use of an explicit balance equation dates back to 1933 [16] and the more implicit use of balance considerations dates back even further to 1906 [17]. The contribution of the present study to the rapidly expanding literature on this topic [1–3,7] is to use the balance equation to estimate the temporal changes in the absolute number of susceptibles from serological data (Appendix A) and to use the balance equation to estimate the time-varying notification ratios (Appendix B). We would like to point out that the use of the balance equation can be extended even further by keeping track of susceptibles per age-cohort [18] or per region [3]. This makes it possible to infer age-specific or region-specific notification ratios and identify the age groups or regions that are least protected. A reason for not including this refinement here is that the Dutch case-notification data is structured by age and by geographical location only from 1988 onward: the time series is too short to obtain reliable estimates of the notification ratios.
2648
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
The estimation of the time-varying notification ratios relies on the assumption that the proportion of susceptibles remains stationary over epidemic cycles. That is, we assume similarity of expected proportions of susceptibles that are in the same phase of the epidemic cycle. This is a considerable improvement over the assumption of an equilibrium. The proposed method for estimating notification ratios is at first sight more or less equivalent to the local regression approach proposed by Finkenstädt and Grenfell [2]. This approach implicitly assumes stationarity of absolute number of susceptibles over an arbitrarily chosen period. However, applying the local regression approach to a set of case notification data for subpopulations of varying sizes shows that it is the proportion of susceptibles, rather than the absolute number, that remains stationary [3]. In this sense, the assumptions underlying the local regression approach are inconsistent with its results, whereas our approach does not have such a problem. As in the local regression approach, the assumption of stationarity over epidemic cycles does not hold in all situations. For example, it does not hold when the vaccination programme has eliminated measles and an epidemic cycle is absent. The excellent agreement between the estimated notification ratios and the available independent information suggests that at least for measles in The Netherlands the assumption of stationarity holds. The estimation of the basic reproduction ratio, the critical vaccination coverage, and the threshold for susceptibles, relies on the assumption of a homogeneously mixing population. This assumption is motivated by parsimony, as homogeneous mixing is the simplest assumption that produces meaningful and useful estimates. However, heterogeneous mixing may be more realistic as it accounts for the socio-geographical clustering of unvaccinated individuals in The Netherlands [8] and for the observed distribution of immunity across ages [15]. We are currently analyzing the estimated time series of measles infections under various mixing assumptions, and the preliminary results suggest that some features of the time series are better explained by heterogeneous mixing than homogeneous mixing, although this has little effect on estimated values of the basic reproduction ratio and the critical vaccination coverage. Notwithstanding these reservations concerning the homogeneous mixing assumption, the practical usefulness of the resulting estimates should not be underestimated. The pattern of major outbreaks is consistent with the pattern of excess of susceptibles over the threshold value, that is, with the pattern of losing herd immunity according to the homogeneous mixing assumption. For practical purposes, the assumption of homogeneous mixing allows us to estimate when a vaccinated population is protected by herd immunity and when it is not. Our estimate of the basic reproduction ratio R0 for measles in a vaccinated population is 23.01. This value is higher than earlier estimates based on the average age at infection as observed before introduction of mass vaccination (range: 14–18) [5], but lower than the estimate of 29.9 based on time series of case notifications before introduction of
mass vaccination [3]. The slight discrepancies between the estimated values could result from differences between the estimation methods, differences between the observed populations, and differences between the times of observation. The fact that our estimate for a vaccinated population lies within the range of previously published values for unvaccinated populations suggests that the basic reproduction ratio for measles has not been drastically altered by vaccination. Our estimates of the proportion of measles cases that are notified are remarkably low and show fluctuations over time. Apparently, the number of cases notified by practicing physicians is a rather poor indicator of the number of actual measles infections among their patients. Therefore case-notification data alone gives a distorted picture of the impact of the measles vaccination programme on the number of measles infections. Time series reconstruction, as proposed here, is necessary to correct for the low and fluctuating notification ratios. But the statistical algorithms for time series reconstruction cannot retrieve all the information that is lost due to low reporting levels, and a higher level of reporting would certainly improve the estimates of the basic reproduction ratio and the critical vaccination coverage. An approach that has proven to increase the reporting level is to search actively for measles cases by public health services, another approach might be to increase the motivation for reporting cases by communicating the insights obtained from notification data back to the reporting physicians. The reconstruction of measles dynamics in the vaccinating Dutch population makes it possible to identify some of the causes of the large measles outbreak in 1999–2000. Firstly, the actual vaccination coverage remained near or below the critical level needed for elimination, and the proportion of susceptibles increased beyond the threshold level, such that from 1997 onwards a measles outbreak could occur. Secondly, the rate of import of measles infections is such that previous outbreaks started on average 2 years after the proportion of susceptibles increased beyond the threshold level, and therefore an outbreak could be expected around 1999, approximately 2 years after the threshold was exceeded in 1997. Thirdly, there has been a drastic increase in the proportion of infections that were notified in 1999. Together, these three causes account for the outbreak with the largest number of notified cases since introduction of vaccination. It is of course not implied that these three causes exhaust all possible relevant causes. The reconstruction of measles dynamics also reveals the epidemiology of measles in a vaccinated population in much greater detail than ever before. A striking feature is that regular epidemics are maintained, even when measles transmission is interrupted after each major outbreak. Already in 1956, Bartlett [13] proposed a mathematical explanation for regularity in such recurrent epidemic outbreaks: the epidemic period consists of a fixed duration required for replenishing the proportion of susceptibles up to the threshold level, and a variable duration for waiting until measles is reintroduced and another outbreak starts. Bartlett assumed
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
that all susceptibles are infected during a major outbreak. The reconstructed time series of proportion of susceptibles, as presented here, shows this is not true: the proportion of susceptibles varies during an epidemic cycle only within a very small range. We are currently adapting Bartlett’s mathematical argument to allow for smaller, more realistic, fluctuations in proportion of susceptibles. In conclusion, this work demonstrates that we can overcome the problems of incomplete observation of dynamics of measles infections and susceptibles by combining case-notification data and serological data. Much more information relevant to the evaluation of vaccination programmes can be extracted from the combined data sources than from each source in separation. We hope that the proposed methods of analysis will be applied to other data sets for other infectious diseases, other vaccination programmes and other populations. Reconstructing the dynamics of infectious diseases in vaccinated populations greatly increases our under− (S0 , σe |xk , sk , ηt ) =
kobs 1 k=1
2
log(2π(sk2 + tk σe2 )) +
2649
where σk is the standard deviation of proportion of susceptibles at time tk . The expected number of susceptibles at time tk follows from Eq. (3), and dividing this number by the population size Nk gives the proportion of susceptibles: tk tk tk Bi−m − i=1 ηi Ci − i=1 S0 + i=1 Vi Sk = . (A.2) Nk Nk The variance of proportion of susceptibles at time tk can be split into two factors, one factor is due to the book-keeping error and follows from Eq. (3), the other factor is the observation error sk : σk2 = sk2 + tk σe2 .
(A.3)
Substituting Eqs. (A.2) and (A.3) into Eq. (A.1) and summing over the number of serological surveys, gives the log-likelihood function for initial number of susceptibles S0 and the standard deviation of the book-keeping error σe
2 tk tk tk kobs x − S + B − /N η C − V 0 i−m i i i k k i=1 i=1 i=1 2(sk2 + tk σe2 )
k=1
standing of how vaccination influences the epidemiology of the infectious disease, it indicates when the vaccination programme fails to provide herd immunity, and thus enhances our capability to measure the actual impact of vaccination.
Acknowledgements This manuscript has profited from discussion with Bryan Grenfell at the University of Cambridge, and from comments on an earlier version by Susan van den Hof and Nico Nagelkerke at the National Institute for Public Health and the Environment.
.
(A.4) Maximizing the log-likelihood for S0 and σe gives the maximum likelihood estimates. With the maximum likelihood estimate for S0 and Eq. (3), we can reconstruct the time series of the absolute numbers of susceptibles St . The average proportion of susceptibles is then calculated as T S 1 St E , (A.5) = N T Nt t=0
where T is the length of the time series of case notifications.
Appendix B. Estimation of notification ratios Appendix A. Maximum likelihood estimation of initial number of susceptibles The initial number of susceptibles S0 and the standard deviation of the book-keeping error σe are estimated from serological data. We use the method of maximum likelihood and assume that the notification ratios ηt are given. We have data on serological surveys conducted at times tk (k = 1, . . . , kobs ; kobs denotes the number of serological surveys available). The population sizes at time of each survey are indicated as Nk . For each serological survey we have the mean proportion of susceptibles xk and standard error of the mean sk . For a Gaussian distribution of mean proportion of susceptibles, the log-likelihood function for number of susceptibles at the time of the serological survey Sk takes the following form: −(Sk , σk |xk ) =
1 (xk − Sk /Nk )2 , log(2πσk2 ) + 2 2σk2
(A.1)
The notification ratio at time t, denoted as ηt , is estimated as the ratio of expected number of infections over the number of reported cases. Both numbers are taken over the interval t−d, t+d, where d is a positive integer. We derive this ratio from the balance equation for susceptibles (Eq. (1)) by first computing St−d − St+d and then solving for ηt , assuming that the cumulative error over this interval is negligible. This yields: ηt ≈
(St−d − St+d ) +
t+d
i=t−d+1 Bi−m t+d i=t−d+1 Ci
−
t+d
i=t−d+1 Vi
.
(B.1) The length of the interval t − d, t + d is 2d. We choose values of 2d close to the length of an epidemic interval or a multiple thereof. We find such values by calculating a weight ρd that is approximately proportional to the autocorrelation function of case reports at lag 2d while
2650
J. Wallinga et al. / Vaccine 21 (2003) 2643–2650
remaining positive: ρd = q
T
Ct−d Ct+d ,
(B.2)
t=1
where T is the length of the time series of case notifications and q is a normalization constant, chosen such that all ρd sum to one. For the specific time series of case notifications studied here, the weights are highest around a length of about 6 years, which is the mean epidemic interval. We exclude lengths shorter than 4 years to ensure that the denominator in Eq. (B.1) is large enough, and lengths longer than 12 years because it is unlikely that the notification ratio remains constant over such a long interval. Thus, we consider only lengths between 4 and 12 years (that is, 26 ≤ d ≤ 78). Furthermore, we have to express the unknown term for change in number of susceptibles in Eq. (B.1), (St−d −St+d ), in terms of known variables. We rely on the approximate stationarity of proportion of susceptibles over epidemic cycles. That is, the change in number of susceptibles is approximately the change in population size times the average proportion of susceptibles in the population, E(S/N), for values of interval length 2d close to an epidemic interval or a multiple thereof: S St−d − St+d ≈ E (Nt−d − Nt+d ). (B.3) N The expected value of proportion of susceptibles E(S/N) is calculated according to Eq. (A.5). The estimator for notification ratio is obtained by substitution of Eq. (B.3) into Eq. (B.1) and taking a weighted average over the appropriate values of d: ηt ≈
d=78 d=26
ρd
t+d E(S/N)(Nt−d −Nt+d )+ t+d i=t−d+1 Bi−m − i=t−d+1 Vi . t+d i=t−d+1 Ci (B.4)
References [1] Fine PEM, Clarkson JA. Measles in England and Wales. I. An analysis of factors underlying seasonal patterns. Int J Epidemiol 1982;11:5–14.
[2] Finkenstädt BF, Grenfell BT. Time series modelling of childhood diseases: a dynamical systems approach. Appl Statist 2000;49:187– 205. [3] Bjørnstad ON, Finkenstädt BF, Grenfell BT. Dynamics of measles epidemics: estimating scaling of transmission rates using a time series SIR model. Ecol Monogr 2002;72:169–84. [4] Farrington CP, Kanaan MN, Gay NJ. Estimation of the basic reproduction number for infectious diseases from age-stratified serological survey data. Appl Statist 2001;50:251–9. [5] Anderson RM, May RM. Directly transmitted infectious diseases: control by vaccination. Science 1982;215:1053–60. [6] vanden Hof S, Berbers GAM, de Melker H, Conyn-van Spaendonck MAE. Sero-epidemiology of measles antibodies in The Netherlands: a cross-sectional study in a national sample and in communities with low vaccine coverage. Vaccine 2000;18:931–40. [7] Grenfell BT, Bjørnstad ON, Finkenstädt BF. Dynamics of measles epidemics: scaling noise, determinism, and predictability with the TSIR model. Ecol Monogr 2002;72:185–202. [8] van den Hof S, Wallinga J, Widdowson M-A, Conyn-van Spaendonck MAE. Protecting the vaccinating population in the face of a measles epidemic: assessing the impact of adjusted vaccination schedules. Epidemiol Infect 2002;128:47–57. [9] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical Recipes in C. Cambridge: Cambridge University Press; 1992. [10] Diekmann O, Heesterbeek JAP, Metz JAJ. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J Math Biol 1990;28:365–82. [11] Dietz K. The estimation of the basic reproduction number for infectious diseases. Stat Methods Med Res 1993;2:23–41. [12] Keeling MJ, Grenfell BT. Individual-based perspectives on R0 . J Theor Biol 2000;203:51–61. [13] Bartlett MS. Deterministic and stochastic models for recurrent epidemics. In: Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability. Berkeley: University of California Press; 1956. p. 81–109. [14] Griffiths DA. The effect of measles vaccination on the incidence of measles in the community. J R Stat Soc A 1973;136:441–9. [15] Wallinga J, Lévy-Bruhl D. Gay NJ, Wachmann CH. Estimation of measles reproduction ratios and prospects for elimination of measles by vaccination in some Western European countries. Epidemiol Infect 2001;127:281–95. [16] Hedrich AW. Monthly estimates of the child population ‘susceptible’ to measles, 1900–1931, Baltimore, Md. Am J Hyg 1933;17:613– 36. [17] Hamer WH. Epidemic disease in England. Lancet 1906;1:733– 9. [18] Fine PEM, Clarkson JA. Measles in England and Wales. II. The impact of the measles vaccination programme on the distribution of immunity in the population. Int J Epidemiol 1982;11:15– 25.