An ensemble approach to predicting the impact of vaccination on rotavirus disease in Niger.

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 R0, 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.6-3.7% of the population each year and that a 2-dose vaccine schedule achieving 70% coverage could reduce burden by 39-42%.

Introduction

Diarrheal disease is the second leading cause of death around the world for children under 5 years of age [4]. Though there are many infectious causes of diarrheal disease in children, rotavirus is the leading cause of gastroenteritis [6], [39], [37]. In many countries, better sanitation, hygiene and access to care have reduced the burden of diarrhea [9], [20]. Despite this trend, the proportion of diarrheal hospitalizations attributable to rotavirus increased [28]. 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 [29], [23].

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 [33], [8]. In sub-Saharan Africa, the cause of diarrheal disease is often unknown due to a lack of diagnostic capacity [26]. 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 [17], [14], [5]. 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 likelihood 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 [19]. 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 [17], [7], [5], [35]. 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 [30]. 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 [30] proposed alternative models for boosting of immunity following sequential exposure to rotavirus. Bayesian model averaging (BMA) [2], [16] 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.6–3.7% of the population each year and that a 2-dose vaccine schedule achieving 70% coverage could reduce burden by 39–42%.

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 9590 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. Using the rapid test, 2921 cases tested positive for rotavirus. 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 [27] of children approximately under 5 years old 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 . 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).

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 [30]. 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 [30]. While the structure of our model is the same as [30], there are two distinctions: (1) our focus is on age-structured modeling for children under the age of five instead of across all age groups; (2) we use a Bayesian inferential approach (as opposed to the maximum likelihood approach in [30]). The latter difference is crucial because we use Bayesian methodology to obtain posterior model probabilities for each of five potential model structures given the observations. This not only enables a probabilistic comparison of the different models but also allows for Bayesian model averaging, thereby providing an ensemble-based projection of rotavirus burden as well as the impact of vaccination. 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 including infection period, immunity period, and exposed period in the SIR models are obtained from Table 2 in previous work [30]; these estimates are from data from England and Wales.

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 disproportionately contribute to the force of infection (e.g. because of higher rates of shedding [30]). 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. In all models, we assume that the transmission varies as a cosine function with a period of 1 year; the mean and amplitude of that function are parameters to be estimated in our approach. We further assume that birth rate varies seasonally based on estimates for Niger from the Demographic and Health Surveys [13]. Then the number of observed reported cases is modeled as a negative binomial with mean equal to the expected number of cases from each model. Transmission rate-related parameters, a reporting rate, a dispersion parameter of the above negative binomial, as well as the mean and amplitude of the time-varying transmission function are estimated via Markov chain Monte Carlo (MCMC). The estimated burden over time was then obtained from each fitted model (see supplementary material for details).

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 [22], for low income countries, we assume a seroconversion rate of 63%. Second, a recent study of a 3-dose vaccination strategy in Niger [18] 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.

Results

In this section we describe the results of fitting various models and calculating ensemble estimates across those models using Bayesian model averaging (BMA) [2], [16] (see supplementary material for details). Our BMA-based estimates, which provide a weighted average of estimates across five different models, formalize uncertainty in model selection. This is because the weighting is done by using posterior model probabilities which measure how well each model is supported by the data. 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.

Table 1

For each model we provide posterior model probability (PMP), the basic reproductive number , 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	PMP	R_0	Burden (severe)	Burden (any)
A	0	30.7 (25.8, 34.3)	9.2 (8.1, 10.1)	38.3 (33.6, 42.4)
B	0.06	12.8 (11.9, 13.6)	3.1 (2.8, 3.4)	25.6 (23.1, 28.1)
C	0.92	12.9 (11.6, 14.0)	3.1 (2.8, 3.5)	25.5 (22.8, 28.1)
D	0.01	10.1 (9.4, 11.9)	3.2 (2.8, 3.6)	20.8 (18.6, 23.0)
E	0.01	9.3 (8.3, 10.2)	3.0 (2.7, 3.2)	14.9 (13.8, 16.0)
BMA		13.3 (12.0, 16.7)	3.1 (2.6, 3.7)	25.5 (22.6, 28.9)

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.1% per year develop severe RVGE as estimated by Models B-E, though this estimate is significantly larger for Model A. The basic reproductive number, is the expected number of secondary cases caused by a single infection in a fully susceptible population. is found as the largest eigenvalue of the next-generation matrix [12] and significantly larger for Model A. BMA for burden and are close to those of Model C, which has the highest weight. Values of reported in the literature range from 1 to 100 [10], [31], [30] depending on the assumptions made. We estimate 3.1% of children under 5 years of age will experience severe RVGE per year; this corresponds to an estimate of 25.5% of children under 5 experiencing rotavirus infection per year. This compares to estimates of between 13% and 33% of children under 5y experiencing rotavirus infection in low and middle income countries reported in [15]. In Fig. 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.

Fig. 1

Burden estimates under the five fitted models for 4 districts in Niger. Dashed lines denote 95% confidence interval. Top: weekly reported severe cases of RVGE and model projections. Middle: model projections of all severe RVGE cases. Bottom: model projections of severe RVGE incidence by age. Lines are model projections while points represent observed cases.

Fig. 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.

Fig. 2

Model-averaged (BMA) burden estimates from the five fitted models for 4 districts in Niger. Dashed lines denote 95% confidence interval. Left: weekly reported severe cases of RVGE and model projections. Middle: model projections of all severe RVGE cases. Right: model projections of severe RVGE incidence by age. Lines are model projections while points represent observed cases.

All of the fitted models are able to successfully capture the observed age distribution of cases (Fig. 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 in transmission rate when the seasonal forcing is stronger than that estimated here. BMA shows a similar trend as the Model C, which has the highest weight.

Fig. 3

Distribution of cases across age groups observed in the data (black dots), predicted by the models (solid lines), and predicted 20 years after vaccination has been introduced at 70% coverage (dashed lines).

Projected impact of vaccination

Here we investigate the impact of vaccination based on the seroconversion rate for low socio-economic settings [22]. 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 [18]. 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 (Fig. 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 [25]. 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 (Fig. 4 ). This happens because the incidence in the inter-epidemic trough drops significantly as well. So, the reduction in cases in the low season means a larger build-up of susceptibles, and a higher peak in the high season. The total annual burden is still lower.

Fig. 4

Relative incidence of severe RVGE after vaccination has been introduced into the models assuming 70% coverage, out to five years after vaccination has been introduced. The vaccination has been introduced at 0 year.

Fig. 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 40.1% (indirect effect: 1.4%) reduction of severe RVGE over the long-term. Based on the recent vaccine efficacy study in [18], we predict 31.1% (indirect effect: 1.0%) reduction in RVGE over the long-term. Though we observe that the indirect effect increases for higher coverage levels, the indirect effect is still marginal. Details are provided in the supplementary material.

Fig. 5

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 for the total effect. The vaccination has been introduced at 0 year. Variation in reduction for a fixed (70%) level of coverage is demonstrated.

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 [30]. 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 [24], [36], [38], 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 [34] and assessing the impact of control actions for foot and mouth disease outbreaks [21].

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. Early March, when urban population density is at its maximum due to seasonal rural-urban migration, was the peak season for transmission of measles [3]. Though measles is transmitted through aerosolized droplets, the similarity in the peak seasonality suggests that higher population density may also facilitate transmission of rotavirus. However, it is important to note that there are other factors that could contribute to the seasonality that may affect rotavirus quite differently when compared to measles. Hence, the comparison to measles above may not be clear-cut. We simply point this out as an intriguing parallel to measles (and meningitis, in fact) all of which have different transmission mechanisms, but the same seasonal pattern.

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 9.3 to 12.9 in Table 1, though Model A has a much larger .

There is an observed double peak in incidence (Fig. 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 [11].

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.1% (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. Ref. [1] estimated that rotavirus vaccine could result in 2.46 million childhood deaths between 2011 and 2030. Of course, uncertainty in the seroconversion rate [22] 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 [22] we estimate that 70% coverage could result in 39–42% reduction in severe RVGE in children under 5. Ref. [18] reported a lower efficacy from a 3-dose schedule in Niger; this would lower the projected reduction of severe RVGE to 29–32%. Notably, although BMA estimates a total reduction in yearly cases using both the efficacy reported in [22], [18], 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 [32]. 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, 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.