An Ensemble Approach to Predicting the Impact of Vaccination on Rotavirus Disease in Niger

TL;DR

This study uses an ensemble of Bayesian SIR models to estimate rotavirus burden and vaccination impact in Niger, accounting for model uncertainty and providing projections for disease reduction.

Contribution

It introduces an ensemble modeling approach combining multiple SIR variants with Bayesian inference to improve rotavirus impact predictions in sub-Saharan Africa.

Findings

Current rotavirus burden is 2.9-4.1% annually.

A 70% coverage vaccine could reduce burden by 37-43%.

Ensemble models effectively quantify uncertainty in epidemiological estimates.

Abstract

Recently developed vaccines provide a new way of controlling rotavirus in sub-Saharan Africa. Models for the transmission dynamics of rotavirus are critical both for estimating current burden from imperfect surveillance and for assessing potential effects of vaccine intervention strategies. We examine rotavirus infection in the Maradi area in southern Niger using hospital surveillance data provided by Epicentre collected over two years. Additionally, a cluster survey of households in the region allows us to estimate the proportion of children with diarrhea who consulted at a health structure. Model fit and future projections are necessarily particular to a given model; thus, where there are competing models for the underlying epidemiology an ensemble approach can account for that uncertainty. We compare our results across several variants of Susceptible-Infectious-Recovered (SIR) compartmental models to quantify the impact of modeling assumptions on our estimates. Model-specific parameters are estimated by Bayesian inference using Markov chain Monte Carlo. We then use Bayesian model averaging to generate ensemble estimates of the current dynamics, including estimates of $R_0$, the burden of infection in the region, as well as the impact of vaccination on both the short-term dynamics and the long-term reduction of rotavirus incidence under varying levels of coverage. The ensemble of models predicts that the current burden of severe rotavirus disease is 2.9 to 4.1% of the population each year and that a 2-dose vaccine schedule achieving 70% coverage could reduce burden by 37-43%.

Tables

Table 1. Table 1: For each model we provide posterior model probability (PMP), the basic reproductive number R 0 subscript 𝑅 0 R_{0} , and estimated burden. Burden corresponds to yearly cases with severe RVGE (% of population). The last row corresponds to the model-averaged (via Bayesian model averaging) versions of these estimates.

Model	Probability	$R_0$	Burden
A	0	30.7 (25.8,34.3)	9.2 (8.1,10.1)
B	0.01	13.9 (12.7,15.4)	3.5 (3.1,3.8)
C	0.92	13.4 (11.7,15.3)	3.5 (3.1,3.9)
D	0.03	11.2 (9.4,12.7)	3.6 (3.2,4.1)
E	0.04	10.3 (9.5,12.6)	3.2 (2.9,3.5)
BMA		13.4 (10.3,15.4)	3.5 (2.9,4.2)

Table 2. Table 2: Seasonal variation in birth rate in Niger, estimated from 1980-2000 using Demographic and Health Surveys. Dorélien, ( 2013 ) An amplitude of − .17 .17 -.17 for January tells us the birth rate is 17 % percent 17 17\% below the mean.

Month	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
Amplitude	-.17	.01	.03	.25	.12	.03	-.01	.09	.01	.13	-.31	-.17

Table 3. Table 3: Posterior means and 95% HPD intervals

Model	$b$	$\varphi$	$r$	$\rho$
A	0.50 (0.48,0.51)	7.4 (7.3,7.5)	1.5 (1.4,1.5)	0.039 (0.035,0.044)
B	0.41 (0.37,0.45)	7.4 (7.3,7.6)	2.6 (2.3,2.8)	0.096 (0.089,0.104)
C	0.42 (0.36,0.46)	7.4 (7.3,7.5)	2.5 (2.0,2.8)	0.097 (0.089,0.104)
D	0.30 (0.26,0.34)	7.2 (7.0,7.4)	2.7 (2.5,2.8)	0.098 (0.090,0.107)
E	0.32 (0.27,0.36)	7.1 (7.0,7.3)	2.6 (2.5,2.7)	0.109 (0.099,0.117)

Equations

$${\beta }_i(t)={\beta }_{0i}\left(1+b\text{cos}\left(\frac{2\pi t-52\varphi }{52}\right)\right),$$

$$\frac{d{M}_i}{dt}={\alpha }_{i-1}{M}_{(i-1)}-{\alpha }_i{M}_i+\mu N-\delta M_i$$

$$\frac{d{S}_i}{dt}={\alpha }_{i-1}{S}_{(i-1)}-{\alpha }_i{S}_i+\delta M_i-{\lambda }_i{S}_i+\tau R_i$$

$$\frac{d{I}_s,i}{dt}={\alpha }_{i-1}{I}_{s,(i-1)}-{\alpha }_i{I}_{s,i}+{\lambda }_{s,i}{S}_i-{\gamma }_s{I}_{s,i}$$

$$\frac{d{I}_m,i}{dt}={\alpha }_{i-1}{I}_{m,(i-1)}-{\alpha }_i{I}_{m,i}+{\lambda }_{m,i}{S}_i-{\gamma }_m{I}_{m,i}$$

$$\small C=\left(\begin{array}[]{cccccc}1&1&1&1&1&1\\ 1&1&1&1&1&1\\ 1&1&1&1&1&1\\ 3&3&3&1&1&1\\ 6&6&6&2&1&1\\ 18&18&18&6&3&1\end{array}\right)$$

$$\frac{d{M}_i}{dt}={\alpha }_{i-1}{M}_{(i-1)}-{\alpha }_i{M}_i+\mu N-\delta M_i$$

$$\frac{d{S}_{1i}}{dt}={\alpha }_{i-1}{S}_{1(i-1)}-{\alpha }_i{S}_{1i}+\delta M_i-{\lambda }_i{S}_{1i}$$

$$\frac{d{I}_{1i}}{dt}={\alpha }_{i-1}{I}_{1(i-1)}-{\alpha }_i{I}_{1i}+{\lambda }_i{S}_{1i}-{\gamma }_1{I}_{1i}$$

$$\frac{d{R}_{1i}}{dt}={\alpha }_{i-1}{R}_{1(i-1)}-{\alpha }_i{R}_{1i}+{\gamma }_1{I}_{1i}-\tau R_{1i}$$

$$\frac{d{S}_{2i}}{dt}={\alpha }_{i-1}{S}_{2(i-1)}-{\alpha }_i{S}_{2i}+\tau R_{1i}-{\lambda }_{2i}{S}_{2i}$$

$$\frac{d{I}_{2i}}{dt}={\alpha }_{i-1}{I}_{2(i-1)}-{\alpha }_i{I}_{2i}+{\lambda }_{2i}{S}_{2i}-{\gamma }_{2i}{I}_{2i}$$

$$\frac{d{R}_{2i}}{dt}={\alpha }_{i-1}{R}_{2(i-1)}-{\alpha }_i{R}_{2i}+{\gamma }_2{I}_{2i}-\tau R_{2i}$$

$$\frac{d{S}_{3i}}{dt}={\alpha }_{i-1}{S}_{3(i-1)}-{\alpha }_i{S}_{3i}+\tau R_{2i}+\tau R_{3i}-{\lambda }_{3i}{S}_{3i}$$

$$\frac{d{I}_{3i}}{dt}={\alpha }_{i-1}{I}_{3(i-1)}-{\alpha }_i{I}_{3i}+{\lambda }_{3i}{S}_{3i}-{\gamma }_2{I}_{3i}$$

$$\frac{d{R}_{3i}}{dt}={\alpha }_{i-1}{R}_{3(i-1)}-{\alpha }_i{R}_{3i}+{\gamma }_2{I}_{3i}-\tau R_{3i}$$

$$\frac{d{M}_i}{dt}={\alpha }_{i-1}{M}_{(i-1)}-{\alpha }_i{M}_i+\mu N-\delta M_i$$

$$\frac{d{S}_{1i}}{dt}={\alpha }_{i-1}{S}_{1(i-1)}-{\alpha }_i{S}_{1i}+\delta M_i-{\lambda }_i{S}_{1i}$$

$$\frac{d{E}_{1i}}{dt}={\alpha }_{i-1}{E}_{1(i-1)}-{\alpha }_i{E}_{1i}+{\lambda }_i{S}_{1i}-\xi E_{1i}$$

$$\frac{d{I}_{1i}}{dt}={\alpha }_{i-1}{I}_{1(i-1)}-{\alpha }_i{I}_{1i}+\xi E_{1i}-{\gamma }_1{I}_{1i}$$

$$\frac{d{R}_{1i}}{dt}={\alpha }_{i-1}{R}_{1(i-1)}-{\alpha }_i{R}_{1i}+{\gamma }_1{I}_{1i}-\tau R_{1i}$$

$$\frac{d{S}_{2i}}{dt}={\alpha }_{i-1}{S}_{2(i-1)}-{\alpha }_i{S}_{2i}+\tau R_{1i}-{\lambda }_i{S}_{2i}$$

$$\frac{d{E}_{2i}}{dt}={\alpha }_{i-1}{E}_{2(i-1)}-{\alpha }_i{E}_{2i}+{\lambda }_i{S}_{2i}-\xi E_{2i}$$

$$\frac{d{I}_{2i}}{dt}={\alpha }_{i-1}{I}_{2(i-1)}-{\alpha }_i{I}_{2i}+\xi E_{2i}-{\gamma }_2{I}_{2i}$$

$$\frac{d{R}_{2i}}{dt}={\alpha }_{i-1}{R}_{2(i-1)}-{\alpha }_i{R}_{2i}+{\gamma }_2{I}_{2i}-\tau R_{2i}$$

$$\frac{d{S}_{3i}}{dt}={\alpha }_{i-1}{S}_{3(i-1)}-{\alpha }_i{S}_{3i}+\tau R_{2i}-{\lambda }_i{S}_{3i}$$

$$\frac{d{E}_{3i}}{dt}={\alpha }_{i-1}{E}_{3(i-1)}-{\alpha }_i{E}_{3i}+{\lambda }_i{S}_{3i}-\xi E_{3i}$$

$$\frac{d{I}_{3i}}{dt}={\alpha }_{i-1}{I}_{3(i-1)}-{\alpha }_i{I}_{3i}+\xi E_{3i}-{\gamma }_2{I}_{3i}$$

$$\frac{d{R}_{3i}}{dt}={\alpha }_{i-1}{R}_{3(i-1)}-{\alpha }_i{R}_{3i}+{\gamma }_2{I}_{3i}-\tau R_{3i}$$

$$\frac{d{M}_i}{dt}={\alpha }_{i-1}{M}_{(i-1)}-{\alpha }_i{M}_i+\mu N-\delta M_i$$

$$\frac{d{S}_{1i}}{dt}={\alpha }_{i-1}{S}_{1(i-1)}-{\alpha }_i{S}_{1i}+\delta M_i-{\lambda }_i{S}_{1i}$$

$$\frac{d{I}_{1i}}{dt}={\alpha }_{i-1}{I}_{1(i-1)}-{\alpha }_i{I}_{1i}+{\lambda }_i{S}_{1i}-{\gamma }_1{I}_{1i}$$

$$\frac{d{S}_{2i}}{dt}={\alpha }_{i-1}{S}_{2(i-1)}-{\alpha }_i{S}_{2i}+{\gamma }_1{I}_{1i}-{\lambda }_{2i}{S}_{2i}$$

$$\frac{d{I}_{2i}}{dt}={\alpha }_{i-1}{I}_{2(i-1)}-{\alpha }_i{I}_{2i}+{\lambda }_{2i}{S}_{2i}-{\gamma }_{2i}{I}_{2i}$$

$$\frac{d{S}_{3i}}{dt}={\alpha }_{i-1}{S}_{3(i-1)}-{\alpha }_i{S}_{3i}+{\gamma }_2{I}_{2i}-{\lambda }_{3i}{S}_{3i}$$

$$\frac{d{I}_{3i}}{dt}={\alpha }_{i-1}{I}_{3(i-1)}-{\alpha }_i{I}_{3i}+{\lambda }_{3i}{S}_{3i}-{\gamma }_2{I}_{3i}$$

$$\frac{d{S}_{4i}}{dt}={\alpha }_{i-1}{S}_{4(i-1)}-{\alpha }_i{S}_{4i}+{\gamma }_2{I}_{3i}-{\lambda }_{4i}{S}_{4i}$$

$$\frac{d{I}_{4i}}{dt}={\alpha }_{i-1}{I}_{4(i-1)}-{\alpha }_i{I}_{4i}+{\lambda }_{4i}{S}_{4i}-{\gamma }_2{I}_{4i}$$

$$\frac{d{M}_i}{dt}={\alpha }_{i-1}{M}_{(i-1)}-{\alpha }_i{M}_i+\mu N-\delta M_i$$

$$\frac{d{S}_{1i}}{dt}={\alpha }_{i-1}{S}_{1(i-1)}-{\alpha }_i{S}_{1i}+\delta M_i-{\lambda }_i{S}_{1i}$$

$$\frac{d{I}_{1i}}{dt}={\alpha }_{i-1}{I}_{1(i-1)}-{\alpha }_i{I}_{1i}+{\lambda }_i{S}_{1i}-{\gamma }_1{I}_{1i}$$

$$\frac{d{S}_{2i}}{dt}={\alpha }_{i-1}{S}_{2(i-1)}-{\alpha }_i{S}_{2i}+{\kappa }_1{\gamma }_1{I}_{1i}-{\lambda }_{2i}{S}_{2i}$$

$$\frac{d{I}_{2i}}{dt}={\alpha }_{i-1}{I}_{2(i-1)}-{\alpha }_i{I}_{2i}+{\lambda }_{2i}{S}_{2i}-{\gamma }_{2i}{I}_{2i}$$

$$\frac{d{S}_{3i}}{dt}={\alpha }_{i-1}{S}_{3(i-1)}-{\alpha }_i{S}_{3i}+{\kappa }_2{\gamma }_2{I}_{2i}-{\lambda }_{3i}{S}_{3i}$$

$$\frac{d{I}_{3i}}{dt}={\alpha }_{i-1}{I}_{3(i-1)}-{\alpha }_i{I}_{3i}+{\lambda }_{3i}{S}_{3i}-{\gamma }_2{I}_{3i}$$

$$\frac{d{S}_{4i}}{dt}={\alpha }_{i-1}{S}_{4(i-1)}-{\alpha }_i{S}_{4i}+{\kappa }_3{\gamma }_2{I}_{3i}-{\lambda }_{4i}{S}_{4i}$$

$$\frac{d{I}_{4i}}{dt}={\alpha }_{i-1}{I}_{4(i-1)}-{\alpha }_i{I}_{4i}+{\lambda }_{4i}{S}_{4i}-{\gamma }_2{I}_{4i}$$

$${\xi }_i(t)=0.24\rho {\lambda }_{i,t}{S}_{i,t}$$

$${\xi }_i(t)=\rho {\lambda }_{i,t}(0.13S_{1i,t}+0.03S_{2i,t})$$

$$\sum _{i=1}^6\sum _{t=t^{\ast }}^{t^{\ast }+52}|{\xi }_i(t)-{\xi }_i(t-52)|<ϵ=0.01$$

Full text

Keywords: Rotavirus, Vaccine, Niger, Susceptible-Infectious-Recovered models, Bayesian model averaging

1 Introduction

Diarrheal disease is the second leading cause of death around the world for children under 5 years of age (Black et al., , 2010). Though there are many infectious causes of diarrheal disease in children, rotavirus is the leading cause of gastroenteritis (Bryce et al., , 2005; Unicef et al., , 2010; Tate et al., , 2012). In many countries, better sanitation, hygiene and access to care have reduced the burden of diarrhea (Clasen et al., , 2007; Kilgore et al., , 1995). Despite this trend, the proportion of diarrheal hospitalizations attributable to rotavirus increased between 2000 and 2004 (Parashar et al., , 2006). The recent development of new prophylactic vaccines for rotavirus is a promising advance in the prevention of diarrheal disease and the reduction of overall childhood mortality (Patel et al., , 2009; Madhi et al., , 2010).

Observation of rotavirus dynamics and estimation of the burden of rotavirus disease is limited both by non-specific surveillance and under-reporting. The dynamics of rotavirus transmission must often be inferred from non-specific temporal surveillance of diarrheal disease that includes multiple causes. This is analogous to the dynamics of specific influenza strains, which are commonly inferred from non-specific time series surveillance of influenza-like illness (ILI) that includes infection by multiple influenza strains (influenza A and B), as well as additional viral infections, for example parainfluenza, coronavirus, rhinovirus (Riley et al., , 2003; Chowell et al., , 2011). In sub-Saharan Africa, the cause of diarrheal disease is often unknown due to a lack of diagnostic capacity (Mwenda et al., , 2010). Even when the cause of disease is known, an unknown fraction of cases will occur in the community and never be recorded by the health system, leading to a potentially significant level of under-reporting. Dynamic models in general, and so-called state-space models in particular, have been an important tool in the assessment of disease burden from non-specific or imperfect surveillance (Ionides et al., , 2006; Ferrari et al., , 2008; Bretó et al., , 2009). We estimate the burden of rotavirus in the Maradi region of Niger by synthesizing two sources of data. We use hospital surveillance data collected by Epicentre for the incident cases over time, including lab-confirmation to assess the likelilhood that a case of severe diarrhea is caused by rotavirus. In addition, we use a cluster survey of households conducted to estimate the proportion of diarrheal disease cases in the region seeking care. The latter is used to help estimate the reporting rate.

State-space models rely on the temporal correlation in a dynamic model to make the unobservable true state of the system, that is, the incidence of the pathogen of interest, estimable from noisy or imperfectly sampled data (Jones, , 1993). Thus, the inference about disease burden is conditional on the structure of the underlying dynamic model. For pathogens with well characterized epidemiology, such as measles and influenza, the application of state-space models to infer disease burden and transmission dynamics has become common (Ionides et al., , 2006; Cauchemez and Ferguson, , 2008; Bretó et al., , 2009; Simons et al., , 2012). The dynamics of rotavirus, which itself comprises multiple strains that result in varying levels of cross-protective immunity to other strains, has been variously described by a suite of different models (Pitzer et al., , 2012). Therefore, inference about rotavirus burden is limited both by imperfect surveillance of rotavirus infection and uncertainty about the underlying transmission dynamics. Rather than condition our analysis on any one model, we fit the observed time series to a suite of 5 different model structures and assumptions to account for uncertainty in model parameters as well as the dynamics represented in the models themselves.

While the development of several novel rotavirus vaccines is a promising advance for the control of diarrheal disease in children, the potential impact of the introduction of these vaccines at the population-scale is uncertain. The predicted impact of vaccine introduction may depend both on the efficacy of the vaccine and model structure; for example (Pitzer et al., , 2012) proposed alternative models for boosting of immunity following sequential exposure to rotavirus. Bayesian model averaging (BMA) (Bates and Granger, , 1969; Hoeting et al., , 1999) allows for the integration of predictions of multiple models, weighted by their posterior support, to generate a single ensemble estimate that accounts for uncertainty in model selection. Here, via BMA, we use the ensemble of fitted models to predict the short-term and long-term impact of vaccination on rotavirus incidence. We then estimate the predicted impact using the vaccine efficacy from two different studies. Our ensemble approach predicts that the current burden of severe rotavirus disease is 2.9 to 4.1% of the population each year and that a 2-dose vaccine schedule achieving 70% coverage could reduce burden by 37-43%.

2 Material and Methods

We use data from two sources: a time series of clinic admissions for diarrheal disease and a community based survey of health-seeking behavior. Clinic surveillance covers a collection of health centers and district hospitals from four districts in the Maradi region of Niger including Aguie, Guidan Roumdji, Madarounfa, and the city of Maradi. A total of 9,590 cases of diarrhea in children under 5 were recorded from December 23, 2009 to March 31, 2012 (118 weeks). For each patient age in months, village of origin, date of consultation were recorded. Also noted were potential symptoms including temperature, duration of diarrhea before consultation, presence of blood in the stool, presence and duration of vomiting, and level of dehydration. In each case a rapid test was administered for detecting rotavirus. 2,921 cases tested positive for rotavirus via the rapid test. A subset of 378 cases testing positive for rotavirus were also genotyped. While 32 separate strains were identified, more than two thirds of positive cases were of strains G12P[8] or G2P[4].

The distributed nature of Niger’s healthcare system is a challenge for surveillance. Roughly a third of all health centers in these districts were included. Notably absent were the many local health posts staffed by community health workers. To estimate both the fraction of cases seeking care at a health center, and the fraction seeking any level of care, a second source of data is needed. We use a community survey (Page et al., , 2011) in the region of children under 5 to get estimates of these reporting rates.

A total of 2940 children under 5 were selected for inclusion in the cluster survey from households across the four districts. Clusters were allotted according to the population of each village from census data. Sampling weights accounted for household composition and the relative populations of the districts. Among those surveyed, 1099 caregivers reported at least one episode of diarrhea during the recall period of 27 days. Respondents reported whether they sought care at a health structure. We use the reporting rate of severe diarrhea, which is defined as the presence of acute watery diarrhea and the presence of two or more of the signs of loss of consciousness, sunken eyes, and an incapacity to drink or drinking very little.

From the cluster survey we determine that an estimated total of 42.9% of caregivers who reported severe diarrhea consulted at a health center $\left(95\%\text{CI}:(33.1\%,52.7\%)\right)$. The rest either sought care at a district hospital, local health post or did not seek care at a formal health structure. This estimate is used as a proxy for the reporting rate of rotavirus. More specifically, this information is used to construct an informative prior for our Bayesian approach (as described in the supplementary material).

2.1 Model Overview

We consider a range of dynamic models for rotavirus transmission. Information linking individual-level data on the course of infection to the between-person transmission of rotavirus is lacking, leading to variation in the structure of mathematical models for rotavirus (Pitzer et al., , 2012). Using a range of different models allows us to account for the uncertainty in estimation due to model choice. The five models we consider are SIR-like compartmental models of transmission, building upon the models in Pitzer et al., (2012). We incorporate age into the model with separate compartments for ages from 0-1 month, 2-3 months, 4-5 months, 6-11 months, 12-23 months, and 24-59 months. Fixed parameters are estimated from England and Wales data as described in Pitzer et al., (2012).

Here we very briefly outline the main features of five models, Models A through E, based on the SIR framework. Details of the model and inferential procedure are described in the supplementary material. Model A tracks severe and mild rotavirus separately. Severe infections have larger force of infections than mild infections. Unlike Model A, Models B-E assume successive infections and immunity are obtained through repeated infections. Subsequent infections will have a reduced susceptibility to infection and level of infectiousness. Model C allows for an incubation period of infections as well. In Model D there is no temporary immunity during successive infections and immunity is granted after all repeated infections. Model E assumes that full immunity can be obtained during successive infections. All model parameters are estimated via Markov chain Monte Carlo (MCMC) and estimated burden over time were obtained from each model.

2.2 Vaccination

We assume vaccination imparts immunity comparable to a natural infection, and consider a strategy wherein a first dose is administered at 2 months of age and a second dose is administered at 4 months. The vaccine was assumed to confer protection comparable to protection conferred by primary infection following the first dose. The second dose confers additional protection comparable to that conferred by secondary infection. For Model A, where the risk of infection does not decrease based previous number of infections, a separate input parameter is used for the vaccine efficacy. The vaccine efficacy is set to be equal to the predicted efficacy for Models B-E (see supplementary material for details). We study the effect of the vaccine under varying levels of coverage. The short-term effect of vaccination is assessed by looking at incidence over a five year period following introduction of the vaccine. The long-term effect is measured by the yearly reduction in incident cases of Rotavirus gastroenteritis (RVGE) measured 20 years after introduction of the vaccine. Field efficacy of a multi-dose rotavirus vaccination strategy is uncertain. To reflect this uncertainty, we investigate the impact of vaccination using the value of efficacy from two different studies. First based on the results of (Lopman et al., , 2012), for low income countries, we assume a seroconversion rate of 63%. Second, a recent study of a 3-dose vaccination strategy in Niger (Isanaka et al., , 2017) estimated efficacy of 66.7% with all doses. The details of representing these two estimates of efficacy in the 5 models are presented in the supplement.

3 Results

We fit each model independently and estimate parameters. Then we calculate ensemble estimates using Bayesian model averaging (BMA) (Bates and Granger, , 1969; Hoeting et al., , 1999) to formalize uncertainty in model selection (see supplementary material for details). Posterior model probabilities (PMP) measures how much each model is supported by the data. Following the BMA approach, based on these probabilities, we provide a weighted average of estimates from five different models. There is significant discordance across models in the measures of model fit (Table 1). Model C, the model with incubation periods performs the best. Notably, Model A, the only model that does not allow for successive infections with decreased levels of infectiousness, performs significantly worse as measured by posterior model probability.

3.1 Pre-vaccination

Our fitted models allow us to construct estimates of the burden in these four districts (Table 1). Of children under five, an approximate 3.5% per year develop severe RVGE as estimated by Models B-D, though this estimate is lower for Model E and significantly larger for Model A. The basic reproductive number $R_{0}$ is found as the largest eigenvalue of the next-generation matrix (Diekmann et al., , 1990) and significantly larger for Model A. BMA for burden and $R_{0}$ are close to those of Model C, which has the highest weight. In Figure 1, we plot our model projections with uncertainty for reported cases of rotavirus as well as for all cases of severe RVGE. We also note that Models B-E predict a steep decline in cases in children over 1y of age following the epidemic peak; cases in infants under 1y, by contrast, are predicted to decline more slowly.

Figure 2 shows the BMA-based model projections which are close to those of Model C. However we note that BMA-based projections have wider confidence intervals because averaged projections incorporate model uncertainty.

All of the fitted models are able to successfully capture the observed age distribution of cases (Figure 3), though Models C and E predict noticeably more cases than observed for older children (2-5 years). The models vary in their ability to capture the temporal dynamics. During the second year of hospital surveillance we can see a secondary peak in the number of cases that is not captured by our fitted model, although we did find that the model dynamics can produce this double peak through an interaction of a high birth rate and seasonal variation when the seasonal forcing is stronger than that estimated here. BMA shows a similar trend as the Model C, which has the highest weight.

3.2 Projected Impact of Vaccination

Here we investigate the impact of vaccination based on the seroconversion rate for low socio-economic settings (Lopman et al., , 2012). In the supplementary material we provide the impact of vaccination using a different value of efficacy which is measured based on a 3-dose strategy (Isanaka et al., , 2017). This was qualitatively similar, but quantitatively small compared to the results in the main paper. Vaccination causes a noticeable shift in the age distribution across Models B-E (Figure 3), with a higher proportion of RVGE cases occurring in older children. This has significant benefits when considering the age-specific mortality of rotavirus is higher for children under 2 years of age (Morris et al., , 2012). The BMA-based burden shows a similar trend.

Over the short term, Models A-E predict an overall decline in total burden, but an increase in the magnitude of peak incidence (Figure 4).

Figure 5 provides short term and long term impact of vaccination which are model-averaged values from five different models. The short term trend of vaccination impacts based on BMA is similar to that of Model C. At equilibrium (long term), we can observe the reduction in severe rotavirus cases with higher levels of coverage. For a fixed (70%) level of coverage, we predict 38.9% reduction of severe RVGE $\left(99\%\text{CI}:(37.1\%,42.6\%)\right)$ over the long-term. Based on the recent vaccine efficacy study in Isanaka et al., (2017), we predict 29.6% reduction $\left(99\%\text{CI}:(28.0\%,32.7\%)\right)$ in RVGE over the long-term. Details are provided in the supplementary material.

4 Discussion

Diarrheal disease is a major source of childhood morbidity and mortality. However, the multi-etiology nature of diarrheal disease means that it is difficult, in the absence of lab confirmation, to infer total burden or project the consequences of novel interventions. We have rich but short-term data with which to understand the dynamic process; in combination with survey data on health-seeking behavior; however, we can bring additional information to bear on the observation rate to interpret the patterns from the non-specific clinic surveillance.

For rotavirus, the uncertainty inherent in imperfectly observed incidence is compounded by the lack of a generally accepted model and debate about the underlying mechanisms that drive the epidemiological response (Pitzer et al., , 2012). This motivates an ensemble approach, using a combination of different models along with quantitative surveillance to get practical measures of burden and projections about the operational impact of controls. This multi-model ensemble approach is common in geosciences (McAvaney et al., , 2001; Stainforth et al., , 2005; Tebaldi et al., , 2005), where different assumptions on complex underlying processes can produce different climate projections, which motivates a probabilistic forecast from a variety of models. A competing models approach has been adapted to epidemiological problems as well, such as choosing an optimal strategy for measles vaccination (Shea et al., , 2014) and assessing the impact of control actions for foot and mouth disease outbreaks (Lindström et al., , 2015).

Here we formally address these two sources of uncertainty, using a state-space model to address the problem of incidence from non-specific surveillance data, and comparing the inference from an ensemble of proposed models to address the uncertainty in the underlying dynamics. Our ensemble approach suggests robust support for some general patterns of rotavirus dynamics. The peak transmission is well estimated, with a maximum in early March, with little variation between models. Rainfall, which is a primary driver of seasonality in the region, peaks in August. Bharti et al., (2011) found that early March, when urban population density is at its maximum due to seasonal rural-urban migration, was the peak season for transmission of measles. Though measles is transmitted through aerosolized droplets, the similarity in the peak seasonality suggests that higher population density may also facilitate transmission of rotavirus.

We find the SEIRS structure in Model C (model with incubation period) best explains the observed data. In this model, subsequent infections have decreased levels of infectiousness and lower risk of infection compared to the initial infection. All models except for Model A, which offers the worst fit to the data, include this dynamic. The estimated basic reproductive number is fairly robust across Models B-E. In particular, point estimates for models B-E vary from 10.3 to 13.9 in Table 1, though Model A has a much larger $R_{0}$.

There is an observed double peak in incidence (Figure 1) during the second year of observation which our fitted models do not capture. However, this may be an anomaly, as the double peak is not seen strongly during the first and third years. We note that our models are capable of reproducing this behavior when the seasonal variation in transmission is stronger than the best fit estimate, via an interaction between seasonal effects and the high birth rate in the region. More complex explanations for such double peaks have been observed elsewhere. In cholera, similar to rotavirus in transmission, local ecological variations were responsible for bimodal incidence (de Magny et al., , 2011).

Our estimate of overall burden of severe RVGE is robust across Models B-E. In spite of the fact that the full epidemiological processes are unknown, we can be fairly sure that the total yearly burden among children under 5 is in the vicinity of 3.5% (Table 1). Model A predicts a 3-fold greater incidence of severe RGVE; however, this model has the weakest support and model-averaged burden is similar to Models B-E.

While uncertainty in retrospective dynamics and disease burden can be characterized using different models, additional uncertainty about the efficacy of proposed interventions limits the ability to predict future dynamics and disease burden. Atherly et al., (2012) estimated that rotavirus vaccine could result in 2.46 million childhood deaths between 2011-2030. Of course, uncertainty in the seroconversion rate (Lopman et al., , 2012) and achievable vaccination coverage means that the true benefit of these vaccines is unknown. Here, we used the ensemble prediction to project the potential impact of rotavirus vaccination in the Niger setting under two scenarios for vaccine efficacy; thus integrating both dynamic uncertainty due to different models and sensitivity to the realized effectiveness of a vaccine program. Using a vaccine efficacy derived from Lopman et al., (2012) we estimate that 70% coverage could result in 37-43% reduction in severe RVGE in children under 5. Isanaka et al., (2017) reported a lower efficacy from a 3-dose schedule in Niger; this would lower the projected reduction of severe RVGE to 28-33%. Notably, although BMA estimates a total reduction in yearly cases using both the efficacy reported in Lopman et al., (2012) and Isanaka et al., (2017), it also predicts higher peaks where more cases are observed than pre-vaccination. This short-term difference in cycle amplitude for these models is a phenomenon anticipated by Pitzer et al., (2009). Anticipation of this shift in dynamic regime caused by vaccination may be critical to the interpretation of short-term surveillance as the observation of higher peak incidence following the introduction of vaccination may be wrongly interpreted as a failure in the vaccination program.

Dynamic models are a powerful tool to interpret disease surveillance data and anticipate the potential consequences of interventions. The method we describe here addresses two main sources of uncertainty: imperfectly observed data and scientific uncertainty about epidemiological dynamics. Our methods also allow us to identify key epidemiological interpretations – transmission seasonality and the proportional impact of vaccination – that are robust to model choice, and those that are model dependent, that is, $R_{0}$ and the annual burden of severe RVGE. By assessing the fit of the observed surveillance to each model, we find that these latter measures are robust within the subset of well supported models.

5 Acknowledgments

The authors are grateful to Epicentre for providing the data sets for this research project. MF is funded by a grant from the Ecology and Evolution of Infectious Disease program of the NSF/NIH (award number 1 R01 GM105247-01).

Potential conflict of interest statement: none of the authors has conflicts of interest.

**Supplementary Material for An Ensemble Approach to Predicting the Impact of Vaccination on Rotavirus Disease in Niger

**

Joshua Goldstein, Jaewoo Park, Murali Haran, and Matthew Ferrari

We provide details below about the five different dynamic models along with information about computational methods used to perform inference for each of them. In addition, we describe the implemention of the Bayesian model averaging approach used in the manuscript.

Appendix A Model Details

The structure of the models is given in Figure 6, which we explain in detail below. Common to each of the models we describe, we assume a time-varying transmission rate with a period of one year to account for seasonality,

[TABLE]

where $t$ is time in weeks, $\beta_{0i}$ is the baseline rate for age class $i$, and $b$ and $\phi$ are the amplitude and offset of the seasonal variation.

We also assume the birth rate $\mu(t)$ varies with time. The mean weekly birthrate is estimated by $\bar{\mu}=1/(5\times 52)$. The variation in monthly birth rate is shown in Table 2. Finally, for each model we assume a negative binomial observation process with mean equal to the number of weekly reported cases and dispersion parameter $r$.

We describe in detail the dynamics of each of the five models outlined in Figure 6. Model A (Shim et al., , 2006; Atkins et al., , 2012) is an SIRS model in which severe and mild rotavirus are tracked separately. Severe infections have a longer duration and contribute more to the overall force of infection. Following infection, there is a period of temporary immunity that wanes over time. The model is age structured with age groups 0-1 month, 2-3 months, 4-5 months, 6-11 months, 1 year, and 2-5 years indexed by $i$. The differential equations describing the model dynamics are:

[TABLE]

Movement between age classes occurs at rates dependent on the length of the interval in weeks, $\boldsymbol{\alpha}=\left\{\dfrac{1}{8},\dfrac{1}{8},\dfrac{1}{8},\dfrac{1}{24},\dfrac{1}{48},\dfrac{1}{144}\right\}$. The force of infection for age class $i$ is given by $\lambda_{i}=\displaystyle\sum_{j=1}^{6}\beta_{j}(t)C_{ij}\dfrac{(I_{s}+0.5I_{m})}{N_{j}}$, assuming that relative infectiousness for mild infections is less than for severe RVGE. Here $C_{ij}$ represents the frequency of contact from age class $i$ onto class $j$ (Hethcote, , 1996), and satisfies $f_{i}C_{ij}=f_{j}C_{ji}$ where $f_{i}$ is the fraction of the population in class $i$. We make the simplifying assumption that contact between age groups is homogeneous. With the absence of data on rotavirus infections for children over 5 and adults, we also assume the population of children under 5 is closed and consider child-child transmission only. Infection with rotavirus is typically asymptomatic (Pitzer et al., , 2012) or unreported for older children and adults, but could potentially play a role in transmission. The contact matrix is

[TABLE]

.

The differences in our age groups means that the contact matrix is not symmetric, for example we assume the population from 2-5 years is 18 times larger than the population from 0-1 months.

After a period of maternal immunity $(M_{i})$, individuals can be susceptible $(S_{i})$, infected with either mild $(I_{m,i})$ or severe $(I_{s,i})$ rotavirus, or recovered $(R_{i})$. These represent the number of individuals in each class. In (1) we see how the number of children protected by maternal immunity change over time. Newborns are added to this class at rate $\mu$ and individuals leave this class when maternal immunity wanes with rate $\delta$, where the mean period of maternal immunity is assumed to be 13 weeks ($\delta=\frac{1}{13}$).

When maternal immunity wanes children are susceptible to rotavirus infection. In (2), we see that individuals enter the susceptible class when maternal immunity wanes. They become infected at a rate given by the force of infection $\lambda_{i}$. After recovery, individuals may reenter the susceptible class at rate $\tau$, where the mean period of immunity following infection is fixed at one year ($\tau={1}/{52}$).

Equation (3) models the change in total infections with severe rotavirus. We assume the proportion of infections with severe rotavirus is lower than mild by setting $\lambda_{s,i}=0.24\lambda_{i}$. Individuals leave the infected with severe rotavirus for a mean period of one week ($\gamma_{s}=1$) following which they are considered to be recovered. Similarly, (4) tracks the total infections with mild rotavirus, with $\lambda_{m,i}=0.76\lambda_{i}$ and a mean infectious period of just half a week ($\gamma_{m}=2$).

Only a fraction of infections with rotavirus develop RVGE (fixed at 24%), and we assume only severe cases are reported, so the expected number of reported cases for age class $i$ is given by $0.24\rho\lambda_{i}S_{i}$ where $\rho$ is the reporting rate. We make the simplifying assumption for all models that $\rho$ is constant across time and does not vary by age group.

Model B (Pitzer et al., , 2009) is an SIRS model allowing for successive infections in which a second, third or subsequent infection will have a reduced susceptibility to infection and level of infectiousness. This represents partial immunity granted through repeated infections. Only a fraction of individuals in the a first or second infectious class are assumed to develop severe RVGE. The model dynamics are described by as follows.

[TABLE]

Here in addition to an initial period of maternal immunity, individuals can be in the susceptible, infected, or recovered classes for their first $(S_{1},I_{1},R_{1})$, second $(S_{2},I_{2},R_{2})$, or third and subsequent $(S_{3},I_{3},R_{3})$ infections. The force of infection for age class $i$ is given by $\lambda_{i}=\displaystyle\sum_{j=1}^{6}\beta_{j}(t)C_{ij}\dfrac{(I_{1j}+0.5I_{2j}+0.2I_{3j})}{N_{j}}$, assuming that relative infectiousness decreases for subsequent infections. We assume the relative risk of infection decreases for subsequent infections, setting $\lambda_{2i}=0.62\lambda_{i}$ and $\lambda_{3i}=0.37\lambda_{i}$ as in Pitzer et al., (2012). Only 13% of first infections and 3% of second infections are assumed to develop severe RVGE, based on data from a Mexico cohort study (Velázquez et al., , 1996). So the expected number of reported cases for age class $i$ is given by $\rho\lambda_{i}(0.13S_{1i}+0.03S_{2i})$. Following Velázquez et al., (1996), we assume that the mean infectious period for the first infection is one week ($\gamma_{1}=1$) and for subsequent infections is half a week ($\gamma_{2}=2$).

Model C (de Blasio et al., , 2010) is an SEIRS model, similar to Model B but allowing for an additional exposed or incubation period. Individuals in the exposed class are infected but not yet infectious. The dynamic equations are given by:

[TABLE]

The modeling assumptions are the same as Model B but for the addition of an exposed class for the first, second, or subsequent infections $(E_{1},E_{2},E_{3})$. We assume a mean exposed period of 1 day ($\xi=7$).

Model D (Van Effelterre et al., , 2010) is an SIS model which also allows for successive infections with different levels of infectiousness, but assumes there is no period of temporary immunity following infection. After four infections individuals are assumed to be fully immune to infection. The dynamics are described as follows.

[TABLE]

The force of infection is $\lambda_{i}=\displaystyle\sum_{j=1}^{6}\dfrac{\beta_{j}(t)C_{ij}(I_{1}+0.5I_{2}+0.2I_{3}+0.2I_{4})}{N}$, assuming that relative infectiousness decreases for subsequent infections. We also assume the relative risk of infection decreases for subsequent infections, setting $\lambda_{2i}=0.62\lambda_{i}$ and $\lambda_{3i}=\lambda_{4i}=0.37\lambda_{i}$. Again, we assume only 13% of first infections and 3% of second infections are assumed to develop severe RVGE. So the expected number of reported cases in age group $i$ is given by $\rho\lambda(0.13S_{1i}+0.03S_{2i})$.

Finally, Model E (Atchison et al., , 2010) is an SIR-SIS hybrid wherein following infection, individuals have a chance to either return to the susceptible class or gain full immunity. The equations for the dynamics are as follows.

[TABLE]

The chance of returning to the susceptible class varies by number of previous infections. Following Atchison et al., (2010) we fix $\kappa_{1}=0.62$, $\kappa_{2}=0.65$, $\kappa_{3}=0.85$. The remaining modeling assumptions are the same as for models B-D.

A.1 Computational Details

Denote the observed data by $C=\{C_{i,t};t\in(1,...,t_{obs}),i\in(1,...,6)\}$ where $C_{i,t}$ is the number of reported cases in age group $i$ during week $t$. Cases were observed over $t_{obs}=118$ weeks. Denote the number of cases in age group $i$ during week $t$ predicted by our models by $\xi_{i}(t)$. For Model A,

[TABLE]

While for models B-E,

[TABLE]

For each model, the periodic solution to the system of ODEs specified above determines the number of reported cases in age group $i$ during a given week. Model dynamics are integrated forward using the deSolve (Soetaert et al., , 2010) package in R until a periodic solution is reached. Solutions have a period of one year; that is, starting from arbitrary initial conditions, we run the dynamics forward until our expected number of cases is identical from one 52 week period to the next, to within a small tolerance; i.e.

[TABLE]

In practice, numerical integration for 20 years was enough to ensure the periodic solution was reached. After reaching a periodic solution, the models are integrated forward an additional 118 weeks to get the expected number of reported cases $\left(\Xi_{i}(t);t\in(1,...,t_{obs}),i\in(1,...,6)\right)$.

Define random variables $N_{i,t}\sim NB(\Xi_{i}(t),r)$. The likelihood is

[TABLE]

The number of observed reported cases is modeled as a Negative Binomial with mean equal to the expected number of cases and dispersion parameter $r$.

Inference for our model parameters is done via Markov chain Monte Carlo (MCMC) for models A-E. At each step of the Markov chain, new parameters $\Theta^{\prime}$ are proposed and the model dynamics are integrated forward until the periodic solution $\Xi_{i}(t;\Theta^{\prime})$ is reached in order to calculate $\mathcal{L}(C|\Theta^{\prime})$. The parameters estimated by MCMC are $\Theta=\left(b,\phi,r,\rho,\beta_{0i};i\in(1,...,6)\right)$, including seasonal amplitude $b$, seasonal phase $\phi$, the dispersion $r$ of the Negative Binomial observation process, the reporting rate $\rho$, and the baseline transmission rate for age class $\beta_{0i}$.

MCMC samples are obtained from the posterior distribution

[TABLE]

where we take priors $p(\beta_{0i})=N(20,5)$, $p(b)=\text{Unif}(0,1)$, $p(\phi)=\text{Unif}(2,2\pi+2)$, $p(r)=\text{Gamma}(0.001,0.001)$, and $p(\rho)=N(0.117,0.06)$. The prior of our reporting rate $\rho$ is centered at 11.7%, determined from the estimated reporting rate from the cluster survey (42.9%) and the estimated proportion of the population under 5 in the four districts that is covered by hospital surveillance (27.3%, from 2009 census data). In practice, we find that our estimates are robust to the choice of standard deviation of $p(\rho)$.

Table 3 provides parameter estimates from five different models. Estimates of the strength of transmission are similar for models B-D, higher for Model E and significantly lower for Model A. The same holds true for the reporting rate (Model A’s estimate of the reporting rate is dramatically lower, and does not agree with estimates from the cluster survey, evidence that it is performing poorly). Notably, the estimated phase of the transmission $\phi$ is similar across all models (Table 3). For reference, an estimated $\phi$ of $7.4$ corresponds to a peak transmission in early March. This is quite close to the period of peak night time brightness in Maradi as measured by satellite imagery (Bharti et al., , 2011).

A.2 Dynamics Accounting for Vaccination

Based on the results of (Lopman et al., , 2012) we assume that $63\%$ of vaccinated individuals are successfully seroconvert after a single dose. Our models with vaccination allow for the red transitions in Figure 6. For example, Model B allows for transitions directly from $M_{i=1}$ and $S_{1,i=1}$ to $R_{1,i=2}$ on the first dose, and from $R_{1,i=2}$ and $S_{2,i=2}$ to $R_{2,i=3}$ on the second dose. The dynamics equations will be modified by the following terms:

[TABLE]

Where $c$ is the coverage and $s=0.63$ is the rate of seroconversion (Lopman et al., , 2012) for low socio-economic settings. This means that an individual who is vaccinated with a single dose has a lower risk of infection, comparable to the effect of recovering from a natural infection. Vaccination with a second dose further reduces risk of infection.

In Model A, the risk of infection does not decline with the previous number of infections. Therefore, an additional vaccinated state $V_{i}$ is added to the model for age group $i$. Two additional input parameters are required for the vaccine efficacy against severe and mild RVGE. We assume the vaccination happens at $2$ months, but the vaccine efficacy is equal to the efficacy predicted under models B-E for the two dose strategy, $\eta_{s}=.796$ and $\eta_{m}=.609$.

[TABLE]

Given our vaccination strategy for models B-E, the vaccine efficacy for severe RVGE after two doses is 79.6%, in line with efficacy studies of rotavirus vaccines. This is calculated by multiplying the proportion of individuals who are successfully immunized twice, once, or zero times by the expected reduction in RVGE incidence for each case.

[TABLE]

We assume following Velázquez et al., (1996) that 47% of first infections and 25% of second infections and 32% of third infections are assumed to develop any RVGE (mild RVGE is unreported). Therefore the vaccine efficacy for all RVGE is

[TABLE]

In practice, first model parameters $\Theta$ are estimated via MCMC for the models without vaccination. Using the fitted model, the dynamics are then integrated forward at the posterior mean of $\Theta$ until the periodic solution has been reached. Then, the dynamics are modified to allow for transitions between compartments by vaccination.

A.2.1 Projections Based on Vaccine Efficacy from a Recent Study

Recently Isanaka et al., (2017) estimated that 3 doses of vaccine had 66.7% efficacy against severe RVGE among children in Niger. Though we do not explicitly account for 3 doses of vaccine, we can calculate the effective seroconversion rate for our model above that would yield this observed efficacy after a complete sequence of doses. Thus, we set $\eta_{s}=.667$ and use (5) to calculate the effective seroconversion rate as 49%. Then the vaccine efficacy for all RVGE is $\eta_{m}=.515$. We estimate the predicted impact of vaccination using different $\eta_{s},\eta_{m}$ and $s$ values with the same dynamic equations.

Although we used the two dose strategy, by using different value of the efficacy, our study can account for uncertainty in the seroconversion rate. Figures 7-9 are matched to Figures 3-5 in the main paper. Because of the lower seroconversion rate, the projected results were qualitatively similar, quantitatively smaller. Vaccination causes a shift in the age distribution across models (Figure 7), with a higher proportion of RVGE cases occurring for older children.

Over the short term, Models A-E predict an overall decline in total burden, but an increase in the magnitude of peak incidence (Figure 8).

Figure 9 indicates that the short term trend of vaccination impacts based on BMA is similar to that of Model C. BMA predicts 29.6% of long term reduction $\left(99\%\text{CI}:(28.0\%,32.7\%)\right)$.

Appendix B Bayesian Model Averaging

For $k=1,...,5$, consider $M_{k}$, the $k$th model, with prior $p(\Theta_{k}|M_{k})$ and likelihood function $\mathcal{L}(C|\Theta_{k},M_{k})$. Note that we take the uniform model prior for $p(M_{l})$ and model evidence $P(C|M_{k})$ is approximated via Bayesian information criterion (BIC) as in Raftery, (1996). Then the posterior model probability (PMP) for $M_{k}$ given the observed data $C$ is

[TABLE]

where

[TABLE]

is the model evidence for $M_{k}$ which measures how well each model is supported by the observed data. Then the BMA estimate of the burden is

[TABLE]

A summary of our implementation of BMA is as follows: (1) We construct a separate MCMC algorithm for each of the models A-E. (2) For each model, the burden estimate $\xi_{k}(t)$ is evaluated for the MCMC samples of the posterior distribution of that model. (3) The expected burden for model $k$, $E[\xi_{k}(t)|C,M_{k}]$, is estimated through the sample mean of the $\xi_{k}(t)$s obtained from Step (2). (4) We take the weighted average of the burden across all models, with the weights equal to the posterior model probabilities, $p(M_{k}|C)$, obtained above.