Hostname: page-component-cd9895bd7-dk4vv Total loading time: 0 Render date: 2024-12-27T07:56:45.276Z Has data issue: false hasContentIssue false
  • English
  • Français

Assessment of varicella vaccine effectiveness in Germany: a time-series approach

Published online by Cambridge University Press:  15 December 2010

M. HÖHLE ,
A. SIEDLER ,
H.-M. BADER ,
M. LUDWIG ,
U. HEININGER  and
R. VON KRIES
M. HÖHLE*
Affiliation:
Department for Infectious Disease Epidemiology, Robert Koch Institute, Berlin, Germany Department of Statistics, Ludwig-Maximilians-University München, Munich, Germany
A. SIEDLER
Affiliation:
Department for Infectious Disease Epidemiology, Robert Koch Institute, Berlin, Germany
H.-M. BADER
Affiliation:
Working Vaccination Committee, Ministry for Health, Schleswig-Holstein, Germany
M. LUDWIG
Affiliation:
Kassenärztliche Vereinigung Schleswig-Holstein, Germany
U. HEININGER
Affiliation:
University Children's Hospital Basel, Basel, Switzerland
R. VON KRIES
Affiliation:
Institute of Social Paediatrics and Adolescent Medicine, Ludwig-Maximilians-University München, Munich, Germany
*
*Author for correspondence: Dr M. Höhle, Department for Infectious Disease Epidemiology, Robert Koch Institute, DGZ-Ring 1, 13086 Berlin, Germany. (Email: hoehlem@rki.de)
Rights & Permissions [Opens in a new window]

Summary

A multivariate time-series regression model was developed in order to describe the 2005–2008 age-specific time-course of varicella sentinel surveillance data following the introduction of a varicella childhood vaccination programme in Germany. This ecological approach allows the assessment of vaccine effectiveness under field conditions by relating vaccine coverage in cohorts of 24-month-old children to the mean number of cases per reporting unit in the sentinel network. For the 1–2 years age group, which is directly affected by the vaccination programme, a one-dose vaccine effectiveness of 83·2% (95% CI 80·2–85·7) was estimated which corresponds to previous approaches assessing varicella vaccine effectiveness in the field in the USA.

Keywords

GP surveillance systemsinfectious disease epidemiologystatisticsvaccinesvaricella zoster
Type
Original Papers
Information
Epidemiology & Infection , Volume 139 , Issue 11 , November 2011 , pp. 1710 - 1719
Copyright
Copyright © Cambridge University Press 2010

INTRODUCTION

Varicella vaccination programmes have been implemented in several countries during the last decade [Reference Nguyen, Jumaan and Seward1–Reference Seward4]. Evaluations of the varicella vaccination programme in the USA have demonstrated an impressive effect with respect to varicella-related mortality [Reference Nguyen, Jumaan and Seward1], hospitalization [Reference Reynolds2, Reference Davis, Patel and Gebremariam3] and incidence [Reference Seward4], despite limited vaccine efficacy of one dose [Reference Seward, Marin and Vazquez5] and even two doses of varicella vaccine [Reference Kuter6, Reference Gould7]. However, the epidemiological characteristics for the spread of varicella in unvaccinated populations vary considerably [Reference Nardone8] making assessment of the effects of vaccination programmes on a national level mandatory.

Since varicella is very common in early childhood [Reference Heininger and Seward9] assessment of the incidence of temporal trends of varicella incidence at a national level is not feasible. Therefore trends are either monitored regarding incidence in regional samples [Reference Nguyen, Jumaan and Seward1, Reference Nardone8] or in sentinel surveys [Reference Siedler and Arndt10]. While trends in case ascertainment can be depicted by these approaches, quantifying vaccine effects and age-group interdependencies may require a more sophisticated approach. We therefore attempted to model the time-course of varicella following the introduction of a universal childhood vaccination programme in Germany in 2004 [Reference Siedler and Arndt10] targeted at children aged 1–2 years and sought to apply these analyses to an ecological approach in order to assess varicella vaccine effectiveness.

Measuring the success of an immunization programme under field conditions is a particular challenge which can be addressed by epidemiological means. The study of vaccine effectiveness typically investigates direct effects in individuals using different vaccination models and study designs [Reference Halloran, Longini and Struchiner11]. However, the literature on quantifying population effects of vaccination is less developed and typically occurs within the framework of a detailed dynamic model for disease transmission (see e.g. [Reference Halloran12, Reference Brisson13] for such an approach in the context of varicella vaccination). However, for data originating from routine monitoring such models are too detailed. The screening method [Reference Farrington14] is one approach operating under such monitoring conditions, but it requires knowledge about both the vaccination status of cases and the proportion of population vaccinated, which are not always available. Furthermore, the approach does not specifically address the time-varying aspect of coverage following the introduction of a vaccination programme. Our selected ecological approach provides a solution for estimating population-level effects of vaccination under field conditions by using a coarse disease transmission model based on time-series methodology which indirectly takes herd immunity effects into account. Classical time-series modelling has already received some attention in analysing vaccine effects in populations [Reference Anderson, Grenfell and May15–Reference Mullooly19]. However, previous analyses are typically restricted to univariate analyses with information on vaccination coverage entered only in binary form. Our approach continues along these lines and is novel in as much as it considers coverage percentage as an exogenous time series and, being a multivariate time-series model, handles progression in different age strata.

Given the above, our research question can be formulated as follows: does modelling the time-course of varicella following the introduction of the varicella vaccination programme in Germany allow for quantification of age-specific effects and for an ecological assessment of vaccine effectiveness? Specifically, we wanted to quantify direct vaccine effects in the vaccinated 1–2 years age group by using seasonal and autoregressive terms in the model as proxy for adjusting for indirect vaccine effects in all regarded age groups.

METHODS

Data sources

A convenience sample of 1176 primary-care physicians, consisting of 57% paediatricians and 43% general practitioners (GPs) which accounted for about 15% of all German paediatricians and about 1% of all German GPs in private practice was the sampling frame for the sentinel network. Sentinel network physicians of both groups are distributed in German federal states in the same proportion as the total number of respective physicians in private practice. Physicians reported aggregated monthly numbers of varicella cases by age group. Clear case definitions for reporting were provided. A case of varicella is defined as a person presenting with a clinical picture compatible with varicella, i.e. the presence of skin exanthema and concomitant presentation of papules, blisters, pustules, crusts [Reference Siedler and Arndt10]. Zero-reporting and active reminders were included, but as physicians participated voluntarily without any incentives for their reports, the additional workload for them needed be kept as low as possible. Hence, vaccination status of cases was not recorded and only aggregated case numbers by age groups were reported. Furthermore, as patients are free to choose and change their physician and only a sample of physicians report to the sentinel network, the population size under surveillance can not be defined [Reference Siedler and Arndt10].

Vaccine uptake was extrapolated from health insurance data in one [Schleswig-Holstein (SH)] of the 16 German federal states. Vaccination billing for about 90% of all children in Germany is handled by the Association of Statutory Health Insurance Physicians (‘Kassenärztliche Vereinigungen’) in the different German federal states [Reference Kalies20]. Children in the health insurance system were traced from birth to age 24 months upon billing for the well-baby check visits. Varicella vaccination status by age 24 months was assessed for consecutive birth cohorts from 2003 to 2006 and this constituted the numerator for assessing vaccination coverage for the different birth cohorts.

Analyses

Data from all paediatricians and GPs participating in the sentinel network were used in the analysis. As response for the time-series modelling we chose the monthly ‘mean number of varicella cases per reporting unit’. This quantity accounts for the varying number of reporting units in the sentinel network, whereas the size of the reporting units is considered to be stable over time. The dynamics of the response variable over time is described by a two-phase model: at the first level the mean number of varicella cases at time t (measured in months) per sentinel reporting unit was modelled. This was done using a linear predictor with intercept, linear trend, seasonal terms and an additional ARMA(p,q) model with p autoregressive (AR) and q moving-average (MA) terms. At the second level the age distribution of the observed cases at time t in the <1, 1–2, 3–4, 5–9 and >9 years age groups was handled by a multinomial logistic regression model (for further details on the modelling see the Appendix). For levels 1 and 2, Akaike's Information Criterion (AIC) was used to select from a set of candidate models. A strong correlation between vaccination coverage and time (measured in months) was expected. Including both month and vaccination coverage as covariables would thus lead to strong collinearity and hence lead to stability problems if they are both included. Therefore, the following strategy was applied: fit separate models including either month or coverage as covariate and use AIC to select the better model.

We used R [21] – a free software environment for statistical computing and graphics – to perform all statistical computations.

In order to estimate vaccination coverage nationally, coverage estimates for SH were extrapolated to the whole of Germany as follows. For each of the 16 federal states the SH time-series was shifted in time such that the time of reimbursement in SH (August 2005) was moved to the corresponding month of reimbursement in the federal state. The resulting 16 time series were then averaged by weighting with the number of sentinel network units for the respective federal state in that month.

Vaccine effectiveness is conventionally defined as the reduction between the risk of unvaccinated and vaccinated individuals divided by the risk of unvaccinated individuals [Reference Halloran, Longini and Struchiner22], i.e. 1 minus the relative risk. In our model we use a before-and-after approach to estimate vaccine effectiveness under field conditions. We consider 0% coverage as baseline – representing the unvaccinated situation – and compare this with 100% coverage, corresponding to the vaccinated situation. Using this approach, vaccine effectiveness can be estimated in the unstratified situation as 1 minus the relative reduction in the expectation of the response, i.e. 1 – μ100% coverage/μ0% coverage=1 – exp(β1)100, where β1 represents the effect of coverage in the level 1 model described in the Appendix. Since coverage is only available for the 1–2 years age group, we obtained an estimate of vaccine effectiveness for this specific age group by computing the relative reduction in expectation for this age group, i.e. 1 – μ100% coverage,1–2 yr/μ0% coverage,1–2 yr (see Appendix for details on how these expectations are obtained by a combination of level 1 and level 2 models).

RESULTS

Figure 1 shows the monthly vaccination coverage among the 24-month-old birth cohort in SH. Here, reimbursement started in August 2005. Also shown in Figure 1 is our extrapolated national coverage time-series; the extrapolated coverage percentage at the beginning of the sentinel network in April 2005 was 11%. The corresponding coverages were 38%, 63%, 78%, respectively, at 1, 2 and 3 years after the sentinel network started. A rough comparison (not shown) with differently collected yearly coverage rates for the 24-month-olds from the 2004 and 2005 birth cohorts in an available subset of ten federal states showed that our extrapolation provided a fair description of the coverage of these two cohorts.

Fig. 1. Vaccination coverage in the 24-month-old birth cohort in Schleswig-Holstein (SH) and the extrapolated coverage for Germany.

The mean number of monthly reporting units in the sentinel network for the available 42 months of monitoring was 679 (s.d.=41). Figure 2 shows the observed monthly time series of mean number of cases per reporting unit together with the fit of our level 1 model containing intercept, trend and periodic function, f(t), consisting of two harmonics with frequencies of 12 and 6 months, respectively. An ARMA(2,1) model captured the remaining autocorrelation of the error terms. This level 1 model was determined as the best-fitting model according to AIC in the investigated set of candidate models. The estimated coefficient for the time trend in this model is −0·023 (95% CI −0·025 to −0·020), i.e. the reduction in the mean number of cases per reporting unit is significant (P<0·001) and corresponds to a decrease in the mean number of cases per reporting unit by a factor of exp(−0·023)=0·978 per month or exp(−0·023×12)=0·762 per year.

Fig. 2. Observed and model-fitted mean number of varicella cases per reporting unit in the sentinel network.

Figure 3 shows the model-fitted values combining level 1 and level 2 models (using month as linear covariate). Also shown are the estimated time trends, μt,i, for each of the five age groups (<1, 1–2, 3–4, 5–9, >9 years) as described in the Appendix.

Fig. 3. Two-level model-fitted mean number of cases for each of the five age groups. Also shown are the trends μt,i (smooth lines) for each age group i and time-point t obtained by combining the intercept and trend components of the level 1 and level 2 models.

A decline in all five age groups was observed. For better interpretation, Table 1 shows the model-predicted yearly relative reductions for the first three full years of the sentinel network. As expected, the highest decline is observed in the 1–2 years age group, but all other age groups also experienced a decline. For example, the percentage reduction over 3 years in the <1 year age group is given by 1−0·77×0·76×0·76=56%. The corresponding reductions in the other groups were 69%, 62%, 42% and 30%. Since no vaccinations were performed in the <1 year age group, the reduction here may be entirely due to herd effect.

Table 1. Model-predicted relative reduction factors for each of the five age groups

However, the extrapolated national vaccination coverage percentage shown in Figure 1 has a Pearson correlation coefficient with month of 0·97. Using either month or coverage percentage in the model shows that the model with the covariate coverage provides the better fit (AIC of −44·66 vs. −39·73). When replacing the term β1t with β1(coverage in %) in the level 1 model, the estimated coefficient of coverage is −0·013 (95% CI −0·014 to −0·011), which means that the coverage effect is significant (P<0·001) and for each percentage point of vaccination coverage in the 24-month-old cohort the mean number of cases per reporting unit is reduced by a factor of 0·987. While this approach describes the overall impact of the vaccination programme, the assessment of vaccine effectiveness requires confining the analysis to vaccinated age groups only.

To specifically target the reduction in the 1–2 years age group, for whom the vaccine uptake could be estimated and who experience the direct vaccine effects, we combined the above findings with the level 2 model where the linear time trend is also replaced with coverage, i.e. coverage is used as single linear covariate in both level 1 and level 2 models, but the harmonic functions at both levels remain functions of time. This allows us to predict the reduction factor such as μ100%,1–2 yr/μ0%,1–2 yr for the 1–2 years age group, i.e. the reduction of a hypothetical 100% vaccination coverage in the 24-month-old birth cohort compared to a coverage baseline of 0%. This comparison mimics the conventional vaccine effectiveness assessment comparing attack rates in vaccinated and unvaccinated individuals. Table 2 gives the results for coverage percentages between 60% and 100%. This yields a vaccine effectiveness of 83·2% (95% CI 80·2–85·7) in the 1–2 years age group.

Table 2. Model-predicted reduction factors for five hypothetical coverage percentages

CI, Confidence interval.

DISCUSSION

Since the implementation of most vaccination programmes is progressive with gradual increments of vaccination coverage, the decrease of the targeted diseases due to direct and indirect vaccine effects will also be progressive. During the implementation of a universal varicella childhood vaccination programme in Germany we demonstrated that modelling of the time-course of disease allows estimation of the reduction for different age groups for comparison of direct and indirect vaccination effects. With further inclusion of vaccine uptake, an ecological estimate for vaccine effectiveness was possible and yielded almost identical results to other varicella vaccine effectiveness estimates.

Although year-to-year variation in the incidence of varicella (measured by GP visits or hospitalizations) has been reported from various countries and settings, the predicted and observed decline of varicella cases in the German sentinel network system, mainly in the 1–2 years age group, is consistent with an effect of vaccination rather than a secular trend. Given the short time-span of the monitoring period it appears unlikely that the reduction is attributable to secular changes in factors other than vaccine use [Reference Mullooly19]. Other epidemiological investigations and modelling approaches [Reference Bramley and Jones23–Reference Garnett and Grenfell26] also indicate a lack of large multi-year periodicities. However, our approach remains sensitive to unexplained trends in incidence, whereas ordinary year-to-year fluctuations are taken into account by the stochastic nature of the model.

Early detection of a herd immunity effect in children aged <1 year and the finding of similar vaccine effectiveness estimates for the ecological compared to other approaches are important findings with respect to varicella vaccination programmes.

Impressive herd immunity effects of varicella vaccination programmes have previously been shown in adults, e.g. a reduction of 74% over a 10-year period with increasing childhood vaccination coverage in the USA [Reference Marin, Meissner and Seward27]. Even more impressive reductions than in adults or adolescents have been reported in infants too young to be immunized themselves [Reference Nguyen, Jumaan and Seward1, Reference Anderson, Grenfell and May15, Reference Grijalva16]. These impressive indirect effects were achieved over a 6-year period with >90% coverage since programme implementation [Reference Quian28] and over 10 years with vaccination coverage >80% achieved after the fifth year [Reference Guris29]. Our results show, that even with coverage <80% first effects of indirect protection might occur.

Scientifically more relevant, however, might be the suggested ecological approach to estimate vaccine effectiveness based on time-series methodology. Using vaccine uptake in the model instead of time to estimate the reduction of reported varicella cases per increment in vaccine uptake appears justified because the model fit for vaccine uptake was considerably better than the fit for year of observation. A hypothetical vaccine uptake of 100% would reflect the effectiveness of the vaccine in vaccinated individuals. Vaccination status at age 24 months was used as an estimate for the 1–2 years age group, which covers children aged 12–35 months. We expected that choosing coverage at the midpoint of the age group balances between lower and higher coverage rates in children aged 12–23 and 24–35 months. The reduction of the reported varicella cases for each 1% increment in vaccination coverage was used to estimate vaccine effectiveness from an ecological approach. As the time-series model contains autoregressive terms, herd immunity effects are indirectly addressed in the modelling. A decrease in the mean number of cases (and hence the infection pressure) during previous months implies a smaller contribution from the past for the current response.

Estimates of varicella vaccine effectiveness vary depending on the setting of the assessment [Reference Bayer30], severity of the infection [Reference Seward, Marin and Vazquez5], number of vaccine doses [Reference Kuter6] and duration of follow-up [Reference Chaves31]. Most of these were based on outbreak investigations [Reference Bayer30, Reference Spackova32], a cohort study [Reference Clements33] a household contact study [Reference Seward34] or a case-control study [Reference Vazquez35]. The median effectiveness has been estimated at 72% and 71% [Reference Bayer30, Reference Spackova32] for all and 96·5% or 89% [Reference Seward, Marin and Vazquez5, Reference Bayer30], respectively, for moderate to severe infections from outbreak investigations. The respective figures were 83% and 100% in a prospective cohort study [Reference Clements33], 85% and 97% in a case-control study confined to laboratory-confirmed cases [Reference Vazquez35] and 79% and 92% in a household contact study for moderate cases [Reference Seward34]. The German varicella childhood vaccination programme was initiated with one-dose varicella vaccination and only recently switched to a two-dose schedule in 2009. The observed vaccine effectiveness of 83·2% in our ecological approach is close to the estimates from outbreaks in the USA which similarly included all levels of severity of varicella.

Vaccination uptake was only assessed in one German federal state; however, vaccination coverage measured for children at school entry is similar in SH to that in Germany on average [36]. The challenge in estimating vaccination coverage from health insurance data, which provide complete information on all vaccines given, is the generation of an appropriate denominator [Reference Kalies20]. The approach to include children with follow-up information up to age 24 months allows defining of a valid denominator. Almost 90% of children have completed all well-baby check-up visits scheduled within the first 2 years of life. If children without complete well-baby check-up visits were less completely vaccinated, this would account for an overestimation of the vaccination coverage by a maximum of 10%. However, formal assessment of this figure in SH, found an overestimation of only 1% [Reference Bader and Ludwig37]. As a sensitivity analysis incorporating uncertainty from well-baby check-ups and the extrapolation to other federal states we investigated the following: if using 90% of the coverage rate observed in Figure 1 as covariate, the resulting estimate for vaccine effectiveness would be 86·3%. As a consequence, our reported figure of 83·2% is conservative with respect to potential overestimation of coverage.

Our time-series modelling makes up an easy to implement alternative to assessing vaccine effectiveness. The similarity of our ecological estimate of varicella vaccine efficacy compared to other approaches can be explained by absence of secular trends and consistent reporting. Lack of large multi-year periodicities is suggested from other epidemiological investigations and modelling approaches [Reference Guris29, Reference Guris34–36]. Furthermore, there was no specific varicella catch-up vaccination programme in Germany. Even though vaccinations have occurred in the 3–4, 5–9 and >9 years age groups due to delay between recommended and actual vaccination age [Reference Kalies20], the effects can be considered as minor. Finally, ascertainment of cases and vaccination coverage was consistent over time due to the structure of the reporting instruments.

A limitation of the current modelling is that coverage percentage enters the predictor of the log-linear model in linear form. Predictions about the effect of hypothetical 100% coverage are sensitive to this assumed parametric form. An alternative could be to use more flexible functional forms based on semi-parametric penalized spline approaches [Reference Fahrmeir and Tutz38], but in our case such modelling did not indicate strong deviations from the linear form. Furthermore, the multinomial phase 2 model implies a fixed variance-covariance structure of the count data series, which – similar to binomial time series – might be subject to overdispersion [Reference Mebane and Sekhon39]. An alternative modelling approach based on a multivariate AR(1) Poisson model with population offsets can be found in Herzog et al. [Reference Herzog, Paul and Held40]. However, their model does not specifically address the sentinel network situation found in our context and is less flexible with respect to the autoregressive modelling.

It should be noted that standard epidemic models suggest that the number of new cases is proportional to the product of the proportions being infectious and susceptible. Thus, our use of infection pressure in the autoregression represents a pragmatic approximation when no information about the susceptible proportion is available. Theoretically, changes in the susceptible proportions due to vaccination could be reflected by a time-varying autoregression parameter. Furthermore, as herd effects are reflected in the unstratified level 1 model, we implicitly assume that these effects are homogeneous over all age groups.

Finally, given the short period of the German vaccination programme, our analyses of vaccine effectiveness do not take into account that immunity from vaccination might wane over time which would confound our effectiveness estimate. Re-assessment of waning immunity might become necessary, once the programme has been in operation for a longer time.

Altogether, our vaccine effectiveness estimates correspond to immediate effects and thus do not reflect long-term impact. It is well known that vaccination effects can be exaggerated during the initial period of a vaccination programme, e.g. due to ‘honeymoon effects’ [Reference McLean and Anderson41] or boosting from natural varicella cases which is expected to decline in the future [Reference Goldman42]. Hence our estimates may not be directly comparable with effectiveness based on individual-level data since they are subject to transient herd effects. However, the proposed two-phase multivariate time-series modelling represents a novel view on investigating population-level vaccine effectiveness for an age-stratified population indirectly taking current herd effects into account. With this model, the dynamics using a monthly time resolution were investigated and quantified. Such modelling is especially helpful in case of short time-series when starting monitoring immunization programmes without available pre-vaccination data.

ACKNOWLEDGEMENTS

We thank the participating sentinel physicians for their data reporting and the German Green Cross for collecting the data in the varicella sentinel. The work of the first author was financially supported by the LMUinnovativ project ‘Munich Center of Health Sciences’.

APPENDIX

Statistical methodology

Allow n t to be the number of sentinel network reporting units at time t; c t the total number of reported varicella cases at this time and (c t,<1, c t,1–2, c t,3–4, c t,5–9, c t,>9)T a column vector of length 5 containing the number of cases in each of the five age groups. The elements of this vector sum to c t. At the first level, our interest is in modelling the response y t=c t/n t, i.e. a random variable reflecting the mean number of cases per reporting unit at time t. To reflect the fact that the response is non-negative, we transform the response using the natural logarithm and hence model log(y t). Now, a time-series model is used to model trend and seasonality:

\log \lpar y_{t} \rpar \equals \beta _{\setnum{0}} \plus \beta _{\setnum{1}} t \plus f\lpar t\rpar \plus \epsiv _{t} \comma

where t=1, 2, … , 42 reflects time (in months) and εt is possibly correlated zero mean Gaussian random variables with variance σ2. Because a log-transformation is used, all effects should be interpreted in multiplicative fashion, i.e. exp(β0) corresponds to the base level and exp(β1) is the factor multiplied on this base level for each additional month. Furthermore, f(t) is a periodical function with period 12 months, e.g. f(1)=f(13). The function f(t) is obtained by combining several sine and cosine terms with different frequencies. An advantage of the above regression formulation is that additional covariates can be easily be added to the model.

The expectation of the response can be calculated as

\mu _{t} \equals E\lpar y_{t} \rpar \equals \exp \left( {\beta _{\setnum{0}} \plus \beta _{\setnum{1}} t \plus f\lpar t\rpar \plus {\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 2$}}\sigma ^{\setnum{2}} } \right).

When using, e.g. vaccination coverage instead of time as linear covariate, vaccine effectiveness is defined as 1 – μ100% coverage/μ0% coverage=1 – exp(β1)100.

Since our data are of time-series nature, we expected the observations y t to be correlated. Even after modelling trend and seasonality as explained above, we expected remaining correlation to be present within the series. If this additional correlation (i.e. autocorrelation) is not taken into account, the standard errors obtained from the modelling would be too optimistic and hence wrong results about significance will be obtained. To capture this autocorrelation, a ARMA(p,q) model with p autoregressive (AR) and q moving-average (MA) terms is fitted to the residuals, εt. This modelling approach is a well-known procedure from time-series analysis for fitting seasonal time-series data [Reference Shumway and Stoffer43].

Estimation of the model parameters β0, β1, σ2, the parameters of f(t) and the coefficients in the seasonal ARMA model given known values of p,q was performed by standard likelihood methods. The appropriate parameters p,q of the ARMA(p,q) time-series model were selected according to AIC in all nine models with p,q being either 0, 1 or 2. We also investigated seasonal ARMA(p,q) models, but as the seasonality is modelled explicitly using f(t), no such additional complexity improved the AIC.

The age-specific dynamics in the mean number of cases per reporting unit is modelled by a two-stage approach. First, the above presented seasonal ARMA model captures the development in the mean number of visits per reporting unit over all age groups. Second, the vector of proportions (p t,<1, p t,1–2, p t,3–4, p t,5–9, p t,>9)T in age groups at a specific time-point is modelled by multinomial logistic regression using category-specific intercept, trend and period function (see e.g. [Reference Mebane and Sekhon39]). Altogether,

Level 1:

\eqalign{\tab{y_{t} \equals \tab {{c_{t} } \over {n_{t} }} \equals \exp \lpar \beta _{\setnum{0}} \plus \beta _{\setnum{1}} t \plus f\lpar t\rpar \plus \epsiv _{t} \rpar \comma }\cr \tab\epsiv_{t} \sim{\rm ARMA}\lpar p\comma q\semi \sigma ^{\setnum{2}} \rpar \comma}

Level 2:

\eqalign{\tab{\lpar c_{t\comma \lt \setnum{1}} \comma c_{t\comma \setnum{1} \minus \setnum{2}} \comma c_{t\comma \setnum{3} \minus \setnum{4}} \comma c_{t\comma \setnum{5} \minus \setnum{9}} \comma c_{t\comma \gt \setnum{9}} \rpar ^{\rm {T}}} \cr \tab {\sim {\rm {Multinom}}}\lpar c_{t} \comma \lpar p_{t\comma \lt \setnum{1}} \comma p_{t\comma \setnum{1} \minus \setnum{2}} \comma p_{t\comma \setnum{3} \minus \setnum{4}} \comma p_{t\comma \setnum{5} \minus \setnum{9}} \comma p_{t\comma \gt \setnum{9}} \rpar ^{\rm {T}} \rpar .}

When estimating parameters at level 2 we used the actual observed c t and not the level 1 predicted c t – this makes it possible to fit the two models independently from each other. However, for the later predictions we combined level 1 and level 2 models using the multiplication y t×(p t,<1, p t,1–2, p t,3–4, p t,5–9, p t,>9)T, where y t is now the predicted value from level 1 and the vector of probabilities is the prediction from level 2. With this combination a time-series model for the mean number of cases per sentinel unit in each age group is obtained. We now calculate

\eqalign{\tab{\mu_{t\comma i} \equals E\lpar y_{t\comma i} \rpar}\cr \tab\hskip14pt \equals {{{\rm exp}\left( {\gamma _{{\rm \setnum{0}\comma }i} \plus \beta _{\rm \setnum{0}} \plus \lpar \gamma _{{\rm \setnum{1}\comma }i} \plus \beta _{\rm \setnum{1}} \rpar \cdot t \plus f\lpar t\rpar \plus f_{i} \lpar t\rpar \plus {1\over 2}\sigma ^{\setnum{2}} } \right)} \over {\sum\limits_{j \in \lcub \lt \setnum{1}\comma \setnum{1} \minus \setnum{2}\comma \setnum{3} \minus \setnum{4}\comma \setnum{5} \minus \setnum{9}\comma \gt \setnum{9}\rcub }{{\rm exp}\lpar \gamma _{{\rm \setnum{0}\comma }j} \plus \gamma _{{\rm \setnum{1}\comma }j} t \plus f_{j} \lpar t\rpar \rpar } }}\comma }

for i=<1, 1–2, 3–4, 5–9, >9, with the β values and the γ values being the intercept and trend parameters in the level 1 and level 2 models, respectively, f(t) and f i(t) being the harmonic functions, and σ2 being the variance of the level 1 model error term. Note that for the reference age group (e.g. 1–2 years), the restrictions γ0,1–2=γ1,1–2=0 and f 1–2(t)=0 for all t applies, i.e. no parameters need to be estimated for this group.

A big advantage of our proposed multivariate two-phase model, compared to modelling the mean number of varicella cases separately for each age group, is that it takes into account, that the number of cases summed over the age groups equals the observed total number of cases. However, when considering the dynamics within a single age group in our model, the interpretation of the model parameters is not straightforward. Instead, we report ratios such as μ13,i/μ1,i giving the relative reduction within the first year of the sentinel network for a specific age group. Confidence intervals for these reduction factors can be obtained by parametric bootstrap exploiting the asymptotic normality of both the β values and the γ values in the computation of the respective μt,i values. When using vaccination coverage as linear covariate instead of time, the relative reduction μ100% coverage,i/μ0% coverage,i denotes the relative reduction in the expectation of the response in age group i – by comparing the 0% coverage with the 100% coverage situation. Again, this can be used to define vaccine effectiveness. In case the harmonic functions of time remain in the predictors of the level 1 and level 2 models, one would typically fix time to some arbitrary value, e.g. t=0, for the comparison.

DECLARATION OF INTEREST

The Varicella Sentinel is part of the work of the Measles and Varicella Working Group (Arbeitsgemeinschaft Masern und Varizellen), which is a joint initiative of the Robert Koch Institute and the vaccine manufacturers, GlaxoSmithKline and Sanofi Pasteur MSD, together with the German Green Cross. The manufacturers finance the work of the German Green Cross, which is responsible for the realization of the project. The Robert Koch Institute supplied scientific guidance and does not receive any financial support from the vaccine manufacturers.

References

REFERENCES

1.Nguyen, HQ, Jumaan, AO, Seward, JF. Decline in mortality due to varicella after implementation of varicella vaccination in the United States. New England Journal of Medicine 2005; 352: 450–458.CrossRefGoogle ScholarPubMed
2.Reynolds, MA, et al. Epidemiology of varicella hospitalizations in the United States, 1995–2005. Journal of Infectious Diseases 2008; 197 (Suppl. 2): S120–126.CrossRefGoogle ScholarPubMed
3.Davis, MM, Patel, MS, Gebremariam, A. Decline in varicella-related hospitalizations and expenditures for children and adults after introduction of varicella vaccine in the United States. Pediatrics 2004; 114: 786–792.CrossRefGoogle ScholarPubMed
4.Seward, JF, et al. Varicella disease after introduction of varicella vaccine in the United States, 1995–2000. Journal of the American Medical Association 2002; 287: 606–611.CrossRefGoogle ScholarPubMed
5.Seward, JF, Marin, M, Vazquez, M. Varicella vaccine effectiveness in the US vaccination program: a review. Journal of Infectious Diseases 2008; 197 (Suppl. 2): S82–89.CrossRefGoogle ScholarPubMed
6.Kuter, B, et al. Ten year follow-up of healthy children who received one or two injections of varicella vaccine. Pediatric Infectious Disease Journal 2004; 23: 132–137.CrossRefGoogle ScholarPubMed
7.Gould, PL, et al. An outbreak of varicella in elementary school children with two-dose varicella vaccine recipients – Arkansas, 2006. Pediatric Infectious Disease Journal 2009; 28: 678–681.CrossRefGoogle ScholarPubMed
8.Nardone, A, et al. The comparative sero-epidemiology of varicella zoster virus in 11 countries in the European region. Vaccine 2007; 25: 7866–7872.CrossRefGoogle ScholarPubMed
9.Heininger, U, Seward, JF. Varicella. Lancet 2006; 368: 1365–1376.CrossRefGoogle ScholarPubMed
10.Siedler, A, Arndt, U. The impact of the routine varicella vaccination programme on varicella epidemiology in Germany. Eurosurveillance 2010; 15(13): pii=19530.CrossRefGoogle ScholarPubMed
11.Halloran, ME, Longini, IM, Struchiner, CJ. Design and Analysis of Vaccine Studies. New York: Springer, 2010.CrossRefGoogle Scholar
12.Halloran, ME, et al. Theoretical epidemiologic and morbidity effects of routine varicella immunization of preschool children in the United States. American Journal of Epidemiology 1994; 140: 81–104.CrossRefGoogle ScholarPubMed
13.Brisson, M, et al. Analysis of varicella vaccine breakthrough rates: implications for the effectiveness of immunisation programmes. Vaccine 2000; 18: 2775–2778.CrossRefGoogle ScholarPubMed
14.Farrington, CP. Estimation of vaccine effectiveness using the screening method. International Journal of Epidemiology 1993; 22: 742–746.CrossRefGoogle ScholarPubMed
15.Anderson, RM, Grenfell, BT, May, RM. Oscillatory fluctuations in the incidence of infectious disease and the impact of vaccination: time series analysis. Journal of Hygiene (London) 1984; 93: 587–608.CrossRefGoogle ScholarPubMed
16.Grijalva, CG, et al. Decline in pneumonia admissions after routine childhood immunisation with pneumococcal conjugate vaccine in the USA: a time-series analysis. Lancet 2007; 369: 1179–1186.CrossRefGoogle Scholar
17.Antunes, JL, et al. Effectiveness of influenza vaccination and its impact on health inequalities. International Journal of Epidemiology 2007; 36: 1319–1326.CrossRefGoogle ScholarPubMed
18.Girard, DZ. Intervention times series analysis of pertussis vaccination in England and Wales. Health Policy 2000; 54: 13–25.CrossRefGoogle ScholarPubMed
19.Mullooly, JP, et al. Evaluation of the impact of an HMO's varicella vaccination program on incidence of varicella. Vaccine 2004; 22: 1480–1485.CrossRefGoogle ScholarPubMed
20.Kalies, H, et al. Vaccination coverage in children can be estimated from health insurance data. BMC Public Health 2008; 8: 82.CrossRefGoogle ScholarPubMed
21.R Development Core Team.R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing, 2009.Google Scholar
22.Halloran, ME, Longini, IM Jr., Struchiner, CJ. Design and interpretation of vaccine field studies. Epidemiological Review 1999; 21: 73–88.CrossRefGoogle ScholarPubMed
23.Bramley, JC, Jones, IG. Epidemiology of chickenpox in Scotland: 1981 to 1998. Communicable Disease and Public Health 2000; 3: 282–287.Google ScholarPubMed
24.Brisson, M, Edmunds, WJ, Gay, NJ. Varicella vaccination: impact of vaccine efficacy on the epidemiology of VZV. Journal of Medical Virology 2003; 70 (Suppl. 1): S31–37.CrossRefGoogle ScholarPubMed
25.Ferguson, NM, Anderson, RM, Garnett, GP. Mass vaccination to control chickenpox: the influence of zoster. Proceedings of the National Academy of Sciences USA 1996; 93: 7231–7235.CrossRefGoogle ScholarPubMed
26.Garnett, GP, Grenfell, BT. The epidemiology of varicella-zoster virus infections: a mathematical model. Epidemiology and Infection 1992; 108: 495–511.CrossRefGoogle ScholarPubMed
27.Marin, M, Meissner, HC, Seward, JF. Varicella prevention in the United States: a review of successes and challenges. Pediatrics 2008; 122: E744–E751.CrossRefGoogle ScholarPubMed
28.Quian, J, et al. Impact of universal varicella vaccination on 1-year-olds in Uruguay: 1997–2005. Archives of Diseases in Childhood 2008; 93: 845–850.CrossRefGoogle ScholarPubMed
29.Guris, D, et al. Changing varicella epidemiology in active surveillance sites – United States, 1995–2005. Journal of Infectious Diseases 2008; 197 (Suppl. 2): S71–75.CrossRefGoogle ScholarPubMed
30.Bayer, O, et al. Metaanalysis of vaccine effectiveness in varicella outbreaks. Vaccine 2007; 25: 6655–6660.CrossRefGoogle ScholarPubMed
31.Chaves, SS, et al. Loss of vaccine-induced immunity to varicella over time. New England Journal of Medicine 2007; 356: 1121–1129.CrossRefGoogle ScholarPubMed
32.Spackova, M, et al. Comparative varicella vaccine effectiveness during outbreaks in day-care centres. Vaccine 2010; 28: 686–691.CrossRefGoogle ScholarPubMed
33.Clements, DA, et al. Postlicensure study of varicella vaccine effectiveness in a day-care setting. Pediatric Infectious Disease Journal 1999; 18: 1047–1050.CrossRefGoogle Scholar
34.Seward, JF, et al. Contagiousness of varicella in vaccinated cases: a household contact study. Journal of the American Medical Association 2004; 292: 704–708.CrossRefGoogle Scholar
35.Vazquez, M, et al. Effectiveness over time of varicella vaccine. Journal of the American Medical Association 2004; 291: 851–855.CrossRefGoogle ScholarPubMed
36.Anon.Vaccine coverage at school health examinations in Germany, 2007 [in German]. Epidemiologisches Bulletin 2009; 16: 143–145.Google Scholar
37.Bader, HM, Ludwig, M. Vaccine coverage from association of SHI physicians billing: regional determination of immunization from birth until the age of 24 months 2003–2007 in the state of Schleswig-Holstein/Northern Germany. In: Abstracts of the 1st National Vaccination Conference, 5–7 March 2009, pp. 22–23 (http://www.nationale-impfkonferenz.de/media/Abstractband_5_03_09_.pdf), 2009.Google Scholar
38.Fahrmeir, L, Tutz, G. Multivariate Statistical Modelling Based on Generalized Linear Models, 2nd edn. New York: Springer, 2001.CrossRefGoogle Scholar
39.Mebane, WR, Sekhon, JS. Robust estimation and outlier detection for overdispersed multinomial models of count data. American Journal of Political Science 2004; 48: 392–411.CrossRefGoogle Scholar
40.Herzog, SA, Paul, M, Held, L. Heterogeneity in vaccination rates explains the size and occurence of measles epidemics in German surveillance data. Epidemiology and Infection. Published online: 12 July 2010. doi:10.1017/S0950268810001664.Google Scholar
41.McLean, AR, Anderson, RM. Measles in developing countries. Part II. The predicted impact of mass vaccination. Epidemiology and Infection 1988; 100: 419–442.CrossRefGoogle ScholarPubMed
42.Goldman, GS. The case against varicella vaccination. International Journal of Toxicology 2006; 25: 313–317.CrossRefGoogle ScholarPubMed
43.Shumway, RH, Stoffer, DS. Time Series Analysis and its Applications. New York: Springer, 2000.CrossRefGoogle Scholar
Figure 0

Fig. 1. Vaccination coverage in the 24-month-old birth cohort in Schleswig-Holstein (SH) and the extrapolated coverage for Germany.

Figure 1

Fig. 2. Observed and model-fitted mean number of varicella cases per reporting unit in the sentinel network.

Figure 2

Fig. 3. Two-level model-fitted mean number of cases for each of the five age groups. Also shown are the trends μt,i (smooth lines) for each age group i and time-point t obtained by combining the intercept and trend components of the level 1 and level 2 models.

Figure 3

Table 1. Model-predicted relative reduction factors for each of the five age groups

Figure 4

Table 2. Model-predicted reduction factors for five hypothetical coverage percentages