Vaccination coverage
time (years)
coverage rate
2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
0%20%40%60%80%100%
●
●
●
●
● ● ● ●
95% credibility interval for φt
mean and range of prior distribution
Figure 5.6: Posterior estimate of the yearly rotavirus vaccination coverage rate in the EFS from 2004 till 2013 among children aged three months. The black line denotes the posterior mean, the shaded area illustrates the 95% equi-tailed credibility intervals. Prior distribution ranges and means are marked grey. Prior distributions for the years following 2012 and later were based on the coverage estimates from 2011 (Dudareva et al., 2012).
5.3 Epidemiological impact of rotavirus routine
(a)
Risk ratio for acquiring infection
ηI
density
0.0 0.2 0.4 0.6 0.8 1.0
012345
(b)
Risk ratio for developing symptoms
ηS
density
0.0 0.2 0.4 0.6 0.8 1.0
0.00.51.01.52.02.53.03.5
(c)
Risk ratio for acquiring symptomatic infection
ηIηS
density
0.0 0.2 0.4 0.6 0.8 1.0
05101520
(d)
Mean duration of immunity loss
1 (52ηW)
density
0 2 4 6 8 10
0.00.10.20.30.4
Figure 5.7: Histogram of vaccine effectiveness parameter distribution according to the posterior sample. Posterior medians and 95% equitailed credibility intervals for each pa-rameter are indicated on the x-axis, respectively. Figures (a) and (b) display the vaccine protection against acquiring rotavirus infection (ηI) and developing symptoms (ηS). Figure (c) provides the combined risk ratio for acquiring symptomatic infection (ηIηS). Figure (d) shows the mean duration of immunity loss in years (52·ηW−1), i.e. loss of one immunity level.
rates, we assumed that the respective demographic processes would evolve stochastically subject to their past trend until 2013.
For each sampled parameter vector and demographic development both the expected and observed incidence for each week and age group was computed, where the latter was sampled from the negative binomial distribution of the observations. Finally, 95% pre-diction intervals were calculated for the expected and observed incidence based on their respective sampled distribution. The complete sample procedure for model prediction is
given by Algorithm 6.
Algorithm 6: Predictive incidence sampling from a posterior sample of model pa-rameters
Input: Θ: sample from the model parameter’s posterior distribution π(· |D) Input: Fκ, Fγ, Fm: distributions for the future development of birth, death and
migration rates.
Input: φlong: long term vaccination coverage level.
Input: K: size of the predictive sample.
Output: nh(k)·Y(k)(j)o
k=1,...,K
nX(k)(j)o
k=1,...,K: predictive samples (of size K) for the expected and reported incidence time series for each age group
• for k= 1 to K do
1. Draw a parameter vector ϑk from the posterior sample Θ and draw samples for the demographic processes subject to Fκ, Fγ, Fm.
2. Compute the expected incidences hϑk(t)·Yϑ(j)
k (t) resulting from ODE system corresponding to model WVβ subject to parameter ϑ(k) and the sampled demographic processes
3. set h(k)·Y(k)(j) =nhϑk(t)·Yϑ(j)
k (t)o
t=1,...,T
4. draw a sample for the observed incidence X(j)(t)∼NegBinhϑk(t)·Yϑ(j)
k (t), dϑ(k)
for t= 1, . . . , T and j = 1, . . . , nD 5. set X(k)(j) =nX(j)(t)o
t=1,...,T
end
5.3.2 Model validation using WFS data
In order to investigate whether the estimated vaccine effectiveness and the overall calibrated model is transferable to other settings, we used the posterior parameter distribution ob-tained using the EFS-data to compute the predictive distribution for the notified rotavirus incidence of the western federal states (WFS) from 2001 until 2013 (using Algorithm 6) and compared with the available data. To do so, we acquired the necessary data on de-mographics (Federal Bureau of Statistics, 2013), vaccine coverage (up to 28% (Dudareva et al., 2012)) and the corresponding notification data (Krause et al., 2007). The only model parameters requiring adjustment were those regarding the reporting rates in the WFS (h, q), for which we used posterior estimates from the model inference in Chapter 4
as seen in Table 4.4.
The resulting predictive distribution together with the data from the WFS is displayed in Figure 5.8. Overall, 87% of the data points are contained in the calculated 95% pointwise prediction bands. The good model fit implies that the employed model structure and parameter estimates are able to reproduce the rotavirus epidemiology also from external settings, which suggests a general applicability of our model and confirms the estimated distributions of our parameters.
020406080
0 − 4 years of age
time (years)
weekly incidence (per 100,000)
2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
95% prediction interval for reported incidence 95% prediction interval for expected incidence mean incidence prediction
reported incidence data
012345
5 − 59 years of age
time (years)
weekly incidence
2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
012345
60 − 99 years of age
time (years)
weekly incidence
2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
Figure 5.8: Model prediction for expected and observed number of weekly reported tavirus incidence including pointwise 95% prediction intervals together with observed ro-tavirus incidence for the three age groups 0-4, 5-59 and 60+ years of age in the WFS from 2001 until 2013.
5.3.3 Investigating demographic uncertainty
To assess the impact of uncertainty around the future demographic development deter-mined by migration, birth and death rates, we fitted suitable stochastic processes to the available demographic data up to 2013.
Considering future migration, we assumed that for each upcoming year the yearly migration rate per age group follows an autoregressive process
mt−m¯ =α(m)(mt−1−m) +¯ t, t i.i.d.
∼ N (0, σm2),
with parameters ¯m, −1 < α(m) < 1 and σm > 0 being calibrated by the corresponding age-specific migration rates in the years 2001-2012. As the migration rates already exhibit
heavy fluctuation even in short periods (Statistisches Bundesamt, 2009) and since the assumption of migration being Markovian appears reasonable we did not account for any higher order correlations within the linear model.
Regarding the age group specific death rates, we modelled the yearly log death rates log(γt) to follow a random walk, i.e. the increments of log(γt) were assumed to be indepen-dently normally distributed
logγt−logγt−1 =t, ti.i.d.∼ N (µγ, σ2γ),
again with meanµγand variance σγ >0 being estimated by the death rates from the years 2001-2012. The linear trend of the death rate process is thus captured in the estimated mean of the increments µγ. Simulating from this fitted model yields a time series for the future log death rates, which continue the linear trend from the past years but with an added random walk controlled by the estimated variance σγ2. More advanced mortility models can be found in, e.g., Lee and Carter (1992).
Finally, for the birth rates κt we used a model which also accounts for seasonality and fertility. As fertility yt we defined the crude number of births per week κt divided by the population aged 20-39 years. We then fitted the following linear model to the weekly fertility rates {yt}t=1,...,T from 2002 till 2012
logyt=α(y)1 +α(y)2 t+α(y)3 cos
2πt 52
+α(y)4 sin
2πt 52
+t,
where t refers to the week number and the residuals {t} follow an autoregressive process of order one, i.e.
t=βt−1+ξt, ξti.i.d.∼ N(0, σ2y).
A posterior mode estimate for (α(y), β) andσy under the assumption of flat priors was then obtained by using the R functionarima(). The log fertility rates yt for future times t > T were simulated from the fitted model enhanced with an added random walk to account for changing trends.
logyt = ˆα(y)1 + ˆα(y)2 t+ ˆα(y)3 cos
2πt 52
+ ˆα(y)4 sin
2πt 52
+
t
X
k=T+1
Wk+t, Wk i.i.d.∼ N(0, σW2 ),
where ˆα(y)i denotes the estimated coefficients and t was simulated from the fitted au-toregressive process defined by ( ˆβ,σˆy). The future weekly number of births was then recalculated by multiplying the simulated fertility rates with the momentaneous model population aged 20-39.
Note that the further the transmission model predicts into the future the more those predictions are subject to uncertainty regarding the demographic development. For a graphical representation of the demographic data and simulations from the fitted models see Figure 5.9.
−4000−200002000
Stochastic evolution of migration
years
absolute yearly migration
2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024 2026 2028
Demographic data Simulations from fitted model 95% pointwise prediction intervals
0.00120.00160.0020
Stochastic evolution of death rate
years
weekly deathrate
2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024 2026 2028
0.00050.00070.00090.0011
Stochastic evolution of fertility
years
fertility rate
2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 2024 2026 2028 2030
Figure 5.9: Demographic data on yearly migration (top), weekly death rates in the age group >80 years (middle), and weekly fertility rates (bottom) together with simulations from the corresponding time series regression models fitted to that data. For model fitting all available data (black) from 2001 until 2013 was used. For each of the three time series 200 sample paths (grey) were simulated into the future from the corresponding models.
The impact of the additional demographic uncertainty was measured by comparing the 95%-prediction bands from the original incidence sampling with model predictions assuming the stochastic demographic development instead of fixed rates remaining at the respective levels from 2013.
5.3.4 Herd immunity
A major goal of the dynamic modelling approach is to assess the indirect effects of a routine vaccination program, i.e. the additional incidence decrease in certain age groups which is not attributable to the direct protection granted by the vaccine but rather to the decreased risk of infection due to the changed prevalence.
To measure the indirect epidemiological effects it is necessary to know which direct pro-tection can be expected from the vaccination programme, i.e. by which proportion is the original incidence expected to decrease. An approximation for this expected decrease can be obtained using the estimates for the vaccine effectiveness parameter η = (ηI, ηS, ηW) regarding protection from symptomatic infection and waning. Since the vaccine induced
immunity wanes over time the expected direct effects differ for various age groups. There-fore, the expected relative incidence decrease was calculated for each agea(a= 0, . . . ,99).
The vaccine protection against symptomatic infection corresponding to compartment Vk (k = 1, . . . ,4) was calculated according to combined estimated risk ratios RR(k)I and RR(k)S for acquiring infection and developing symptoms in case of infection compared to being fully susceptible, respectively. As stated in Section 5.1.1 this relative risk was defined by
RR(k)= RR(k)I ·RR(k)S = (ηIηS)k4 .
Since the effectiveness parameters ηI, ηS are to be considered as random variables within the Bayesian framework we defined the estimator RRd(k) as the posterior mean
RRd(k)=E
(ηIηS)k4
,
where the expectation is with respect to the posterior distribution corresponding to the best fitting model scenario WVβ. Thus, the estimated compartment-specific relative risks against symptomatic infection are as follows:
RRd(4) = 0.041, RRd(3)= 0.085, RRd(2) = 0.175, RRd(1) = 0.414.
However, considering the underlying vaccination model and the most likely model sce-nario WVβ, there is no unique vaccine compartment Vk linked to each age a. Instead individuals of the same age are distributed over the different vaccine compartments Vk (k = 1, . . . ,4) while a complete loss of protection is also possible. The exact distribution is determined by the ODE system (5.3) and depends on the waning rate ηW and the time passed since complete vaccine administration at 4 months of age. We assumed that the mean time since vaccine administration for individuals of agea was approximatelya years.
Assuming that individuals would not leave the vaccine compartments in case of infection, the population distribution over the four vaccine compartments (and eventually S1) after a years coincides with the probability masses of a Poisson distribution with parameter 52aηW, since the waning rate can be interpreted as the jump intensity of a Poisson process (counting the number of jumps) and we are interested in the distribution of this process.
Thus, the ratio pa(Vk) of people ageda years being in Vk is given by
pa(V4) =PPoiss(0; 52aηW), pa(V3) = PPoiss(1; 52aηW), pa(V2) =PPoiss(2; 52aηW),
pa(V1) =PPoiss(3; 52aηW), pa(S1) =PPoiss(≥4; 52aηW),
where PPoiss(·;λ) denotes the probability mass function of the Poisson distribution with rate λ and pa(S1) is the probability of the vaccine immunity being completely vanished.
Using this age-specific compartment distribution and the compartment-specific relative risks against symptomatic infection one can calculate the overall vaccine protection for a given age a by
RR(a) =d
4
X
k=1
pa(Vk)·RRd(k)+pa(S1)·1.
For larger age groups the expected protection against symptomatic infection was averaged over the corresponding age range, e.g. for the age group 5 till 19 years of age
RRd5−19= 1 15
19
X
a=5
RR(a).d
Finally, this yields expected relative incidence decreases for each age group that was con-tained in the data as well as in the model:
RRd0 = 0.042, RRd1 = 0.075, RRd2 = 0.126, RRd3 = 0.199, RRd4 = 0.286, RRd5−19= 0.799, RRd20−39 = 0.998, RRd40−59=RRd60−79=RRd80−99= 1.
Note that these stated relative decreases are based on a complete vaccination coverage.
For incomplete coverage the expected decrease RRd(a)(φ) scales with the vaccine coverage rate φ:
RRd(a)(φ) = 1−φ·1−RR(a)d .
Herd immunity was then calculated as the relative difference, i.e. the ratio, between the model predicted and the expected direct decrease of the annual incidence. The model prediction was measured for the year 2025, 10 years after long term vaccination coverage was achieved, assuming that the annual incidence reached a steady state by then.
5.3.5 Epidemiological results
In the base case we assumed that vaccination coverage for each birth cohort would increase to a level of 90% until the end of 2015 (see Figure 5.10). Considering children aged <5 years, formerly causing 60% of all notified cases, the model predicted an rotavirus incidence decrease to 437 cases per 100,000 population for the year 2020 (95% prediction interval (PI): 283 – 701) whereas the corresponding incidence without vaccination was 2,709 (95%
PI: 1677 – 4053) cases per 100,000 yielding an 84% decrease. In the age group 5-59 years, the model predicted incidences of 107 (95% PI: 84 – 135) and 101 (95% PI: 61 – 154) with and without vaccination program, respectively, suggesting a slight increase. In the group of adults aged >60 years, the annual incidence was predicted to increase from 171 (95%
PI: 96 – 265) to 199 (95% PI: 158 – 248) cases per 100,000. The resulting overall reduction in rotavirus incidence was predicted at 35%.
The rotavirus seasonality was predicted to shift after adoption of routine rotavirus vaccination as seen in Figure 5.11. According to our model, peak incidence before vaccine introduction occurred during March with week 12 of each year whereas after introduction the maximum incidence was predicted to be reached in week 15. Thereby, the relative decrease was predicted to be much more pronounced from January to April compared to the off-season.
From analysing the effect of demographic uncertainty (Section 5.3.3), we found that the width of the weekly prediction intervals for the reported incidence increased by up to
0100200300
0 − 4 years of age
time (years)
weekly incidence (per 100,000)
Licensure of RV vaccines Recommendation of routine RV−vaccination
2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020
95% prediction interval for reported incidence 95% prediction interval for expected incidence mean incidence prediction
reported incidence data
051020
5 − 59 years of age
time (years)
weekly incidence
2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020
051020
60 − 99 years of age
time (years)
weekly incidence
2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020
Figure 5.10: Model prediction of the weekly rotavirus incidence in the EFS for the age groups 0-4, 5-59 and 60+ years of age from 2004 to 2020 with introduction of routine rotavirus vaccination in 2013. The dark grey area provides 95% prediction intervals for the expected incidence h·Y(j) incorporating only parameter uncertainty. The light grey area which provides the reported incidenceX(j)additionally includes uncertainty from seasonal fluctuations. Notified incidence data up to 2013 is given in black.
4% compared to a scenario with the demographic rates remaining at their respective levels from 2013.
The indirect effects of routine rotavirus vaccination were predicted to be most promi-nent in young children (Figure 5.12). At a vaccination coverage of 90% herd protection was predicted to prevent 14% of those childhood rotavirus cases remaining when only ac-counting for direct protection, whereas at 80% coverage herd protection prevents 11% of these cases. In the age-groups 5-59 and>60 years of age indirect effects lead to increased incidences of 14% and 18% above the expected level, respectively.