Diarrheal disease is the second leading cause of death around the world for children under 5 years of age [6] . Though there are many infectious causes of diarrheal disease in children, rotavirus is the leading cause of gastroenteritis [8, 43, 41] . In many countries, better sanitation, hygiene and access to care have reduced the burden of diarrhea [11, 22] . Despite this trend, the proportion of diarrheal hospitalizations attributable to rotavirus increased between 2000 and 2004 [30] . 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 [31, 25] .
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 [35, 10] . In sub-Saharan Africa, the cause of diarrheal disease is often unknown due to a lack of diagnostic capacity [28] . 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 [19, 16, 7] . 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 [21] . 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 [19, 9, 7, 38] . 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 [32] . 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 [32] proposed alternative models for boosting of immunity following sequential exposure to rotavirus. Bayesian model averaging (BMA) [4, 18] 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%.
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 [29] 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 (95%CI : (33.1%, 52.7%)). 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).
We consider a range of dynamic models for rotavirus transmission. Information linking individuallevel 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 [32] . 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 [32] .
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 [32] .
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.
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 multidose 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 [24] , for low income countries, we assume a seroconversion rate of 63%. Second, a recent study of a 3-dose vaccination strategy in Niger [20] 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.
We fit each model independently and estimate parameters. Then we calculate ensemble estimates using Bayesian model averaging (BMA) [4, 18] 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.
Our fitted models allow us to construct estimates of the burden in these four districts ( basic reproductive number R 0 is found as the largest eigenvalue of the next-generation matrix [14] 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. 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.
Here we investigate the impact of vaccination based on the seroconversion rate for low socioeconomic settings [24] . 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 [20] . 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 [27] . 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 ). 
Diarrheal disease is a major source of childhood morbidity and mortality. However, the multietiology 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 healthseeking 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 [32] . 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 [26, 40, 42] , 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 [36] and assessing the impact of control actions for foot and mouth disease outbreaks [23] .
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. [5] 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 [13] .
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. [2] estimated that rotavirus vaccine could result in 2.46 million childhood deaths between 2011-2030. Of course, uncertainty in the seroconversion rate [24] 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 [24] we estimate that 70% coverage could result in 37-43% reduction in severe RVGE in children under 5. [20] 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 [24] and [20] , 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 [33] . 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.
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.
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.
A Model Details Figure 6 : Structure of the compartmental models adapted from [32] .
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,
where t is time in weeks, β 0i is the baseline rate for age class i, and b and φ are the amplitude and offset of the seasonal variation.
We also assume the birth rate µ(t) varies with time. The mean weekly birthrate is estimated bȳ µ = 1/(5 × 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 [37, 3] 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:
Movement between age classes occurs at rates dependent on the length of the interval in weeks, α = 1 8
,
. The force of infection for age class i is given by
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 [17] , 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 [32] or unreported for older children and adults, but could potentially play a role in transmission. The contact matrix is 
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 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ρλ i S i where ρ is the reporting rate. We make the simplifying assumption for all models that ρ is constant across time and does not vary by age group.
Model B [33] 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.
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 λ i = 6 j=1 β j (t)C ij (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 λ 2i = 0.62λ i and λ 3i = 0.37λ i as in [32] . Only 13% of first infections and 3% of second infections are assumed to develop severe RVGE, based on data from a Mexico cohort study [45] . So the expected number of reported cases for age class i is given by ρλ i (0.13S 1i + 0.03S 2i ). Following [45] , we assume that the mean infectious period for the first infection is one week (γ 1 = 1) and for subsequent infections is half a week (γ 2 = 2).
Model C [12] 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:
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 (ξ = 7).
Model D [44] 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.
Model D
The force of infection is λ i = Model E
The chance of returning to the susceptible class varies by number of previous infections. Following [1] we fix κ 1 = 0.62, κ 2 = 0.65, κ 3 = 0.85. The remaining modeling assumptions are the same as for models B-D.
Denote the observed data by C = {C i,t ; t ∈ (1, ..., t obs ), i ∈ (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 ξ i (t). For Model A,
While for models B-E,
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 [39] 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.
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 (Ξ i (t); t ∈ (1, ..., t obs ), i ∈ (1, ..., 6)).
Define random variables N i,t ∼ N B(Ξ i (t), r). The likelihood is
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. 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 φ is similar across all models (Table 3) .
For reference, an estimated φ 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 [5] . 
Based on the results of [24] 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:
Where c is the coverage and s = 0.63 is the rate of seroconversion [24] 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, η s = .796 and η m = .609. 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. Figure 7 : Distribution of cases across age groups observed in the data (black dots), predicted by the models (red lines), and predicted after vaccination has been introduced at 70% coverage (blue lines).
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 ). Year Relative incidence of severe RVGE Model A Model B Model C Model D Model E Figure 8 : Relative incidence of severe RVGE after vaccination has been introduced into the models assuming 70% coverge, out to five years after vaccination has been introduced. The vaccination has been introduced at 0 year. 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 (99%CI : (28.0%, 32.7%)). Long−term reduction of severe RVGE (#cases) Figure 9 : Relative incidence of severe RVGE (Left), percent (Middle) and absolute (Right) long term reduction in cases by coverage for Bayesian model averaging from the five fitted models. Dashed lines denote 99% confidence interval. The vaccination has been introduced at 0 year. Variation in reduction for a fixed (70%) level of coverage is demonstrated.
For k = 1, ..., 5, consider M k , the kth model, with prior p(Θ k |M k ) and likelihood function L(C|Θ 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 [34] . Then the posterior model 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.