Predictive Modeling of Cholera Outbreaks in Bangladesh

Despite seasonal cholera outbreaks in Bangladesh, little is known about the relationship between environmental conditions and cholera cases. We seek to develop a predictive model for cholera outbreaks in Bangladesh based on environmental predictors. To do this, we estimate the contribution of environmental variables, such as water depth and water temperature, to cholera outbreaks in the context of a disease transmission model. We implement a method which simultaneously accounts for disease dynamics and environmental variables in a Susceptible-Infected-Recovered-Susceptible (SIRS) model. The entire system is treated as a continuous-time hidden Markov model, where the hidden Markov states are the numbers of people who are susceptible, infected, or recovered at each time point, and the observed states are the numbers of cholera cases reported. We use a Bayesian framework to fit this hidden SIRS model, implementing particle Markov chain Monte Carlo methods to sample from the posterior distribution of the environmental and transmission parameters given the observed data. We test this method using both simulation and data from Mathbaria, Bangladesh. Parameter estimates are used to make short-term predictions that capture the formation and decline of epidemic peaks. We demonstrate that our model can successfully predict an increase in the number of infected individuals in the population weeks before the observed number of cholera cases increases, which could allow for early notification of an epidemic and timely allocation of resources.


Introduction
In Bangladesh, cholera is an endemic disease that demonstrates seasonal outbreaks [Huq et al., 2005, Koelle and Pascual, 2004, Koelle et al., 2005, Longini et al., 2002. The burden of cholera is high in that country, with an estimated 352,000 cases and 3,500 to 7,000 deaths annually [International Vaccine Institute, 2012]. We seek to understand the dynamics of cholera and to develop a model that will be able to predict outbreaks several weeks in advance. If the timing and size of a seasonal epidemic could be predicted reliably, vaccines and other resources could be allocated effectively to curb the impact of the disease.
Specifically, we want to understand how the disease dynamics are related to environmental covariates. It is currently not known what triggers the seasonal cholera outbreaks in Bangladesh, but it has been shown that Vibrio cholerae, the causative bacterial agent of cholera, can be detected in the environment year round [Huq et al., 1990, Colwell andHuq, 1994]. Environmental forces are thought to contribute to the spread of cholera, evident from the many cholera disease dynamics models that incorporate the role of the aquatic environment on cholera transmission through an environmental reservoir effect [Codeço, 2001, Tien andEarn, 2010]. One hypothesis is that proliferation of V. cholerae in the environment triggers the seasonal epidemic, feedback from infected individuals drives the epidemic, and then cholera outbreaks wane, either due to an exhaustion of the susceptibles or due to the deteriorating ecological conditions for propagation of V. cholerae in the environment. We probe this hypothesis using cholera incidence data and ecological data collected from multiple thanas (administrative subdistricts with a police station) in rural Bangladesh over sixteen years. There have been three phases of data collection so far, each lasting approximately three years and being separated by gaps of a few years; the current collection phase is ongoing. For a subset of these data, Huq et al. [2005] used Poisson regression to study the association between lagged predictors from a particular water body to cholera cases in that thana. This resulted in different lags and different significant covariates across multiple water bodies and thanas. Thus, it was hard to derive a cohesive model for predicting cholera outbreaks from the environmental covariates. Also, there is no easy way to account for disease dynamics in this Poisson regression framework. We want to measure the effect of the environmental covariates while accounting for disease dynamics via mechanistic models of disease transmission. Moreover, we want to see if we can make reliable short-term predictions with our model -a task that was not attempted by Huq et al. [2005].
Mechanistic infectious disease models use scientific understanding of the transmission process to develop dynamical systems that describe the evolution of the process [Bretó et al., 2009]. Realistic models of disease transmission incorporate non-linear dynamics [He et al., 2010], which leads to difficulties with statistical inference under these models, specifically in the tractability of the likelihood. Keeling and Ross [2008] demonstrate some of these difficulties; they use an exact stochastic continuous-time, discrete-state model which evolves Markov processes using the deterministic Kolmogorov forward equations to express the probabilities of being in all possible states. However, that method only works for small populations due to computational limitations. To overcome this intractability, Finkenstädt and Grenfell [2000] develop a time-series Susceptible-Infected-Recovered (SIR) model which extends mechanistic models of disease dynamics to larger populations. A similar development is the auto-Poisson model of Held et al. [2005]. To facilitate tractability of the likelihood, both of the above approaches make simplifying assumptions that are difficult to test. Moreover, these discretetime approaches work only for evenly spaced data or require aggregating the data into evenly spaced intervals. Cauchemez and Ferguson [2008] develop a different, continuous-time, approach to analyze epidemiological time-series data, but assume the transmission parameter and number of susceptibles remain relatively constant within an observation period. Our current understanding of cholera disease dynamics leads us to think that this assumption is not appropriate for modeling endemic cholera with seasonal outbreaks.
To implement a mechanistic approach without these approximations, both maximum likelihood and Bayesian methods can be used. Maximum likelihood based statistical inference techniques use Monte Carlo methods to allow maximization of the likelihood without explicitly evaluating it [He et al., 2010, Bretó et al., 2009, Ionides et al., 2006, Bhadra et al., 2011. Ionides et al. [2006] use this methodology to study how large scale climate fluctuations influence cholera transmission in Bangladesh. Bhadra et al. [2011] use this framework to study malaria transmission in India. They are able to incorporate a rainfall covariate into their model and study how climate fluctuations influence disease incidence when one controls for disease dynamics, such as waning immunity. Under a Bayesian approach, particle filter Markov chain Monte Carlo (MCMC) methods have been developed which require only an unbiased estimate of the likelihood [Andrieu et al., 2010]. Rasmussen et al. [2011] use this particle MCMC methodology to simultaneously estimate the epidemiological parameters of a SIR model and past disease dynamics from time series data and gene genealogies. Using Google flu trends data [Ginsberg et al., 2008], Dukic et al. [2012] implement a particle filtering algorithm which sequentially estimates the odds of a pandemic. Notably, Dukic et al. [2012] concentrate on predicting influenza activity. Similarly, here we develop a model-based predictive framework for seasonal cholera epidemics in Bangladesh.
In this paper, we use sequential Monte Carlo methods in a Bayesian framework. Specifically, we develop a hidden Susceptible-Infected-Recovered-Susceptible (SIRS) model for cholera transmission in Bangladesh, incorporating environmental covariates. We use a particle MCMC method to sample from the posterior distribution of the environmental and transmission parameters given the observed data, as described by Andrieu et al. [2010]. Further, we predict future behavior of the epidemic within our Bayesian framework. Cholera transmission dynamics in our model are described by a continuous-time, rather than a discretetime, Markov process to easily incorporate data with irregular observation times. Also, the continuous-time framework allows for greater parameter interpretability and comparability to models based on deterministic differential equations. We test our Bayesian inference procedure using simulated cholera data, generated from a model with a time-varying environmental covariate. We then analyze cholera data from Mathbaria, Bangladesh, similar to the data studied by Huq et al. [2005]. Parameter estimates indicate that most of the transmission is coming from environmental sources. We test the ability of our model to make short-term predictions during different time intervals in the data observation period and find that the pattern of predictive distribution dynamics matches the pattern of changes in the reported number of cases. Moreover, we find that the predictive distribution of the hidden states, specifically the unobserved number of infected individuals, clearly pinpoints the beginning of an epidemic approximately two to three weeks in advance, making our methodology potentially useful during cholera surveillance in Bangladesh.

SIRS model with environmental predictors
We consider a compartmental model of disease transmission [May andAnderson, 1991, Keeling andRohani, 2008], where the population is divided into three disease states, or compartments: susceptible, infected, and recovered. We model a continuous process observed at discrete time points. The vector X t = (S t , I t , R t ) contains the numbers of susceptible, infected, and recovered individuals at time t, and we consider a closed population of size N such that N = S t + I t + R t for all t. Individuals move between the compartments with different rates; for cholera transmission we consider the transition rates shown in Figure 1. In this framework, a susceptible individual's rate of infection is proportional to the number of infected people and the covariates that serve as proxy for the amount of V. cholerae in the environment. Thus, the hazard rate of infection, also called the force of infection, is βI t + α(t) for each time t, where β represents the infectious contact rate between infected individuals and susceptible individuals and α(t) represents the time-varying environmental force of infection. Possible mechanisms for infectious contact include direct person-to-person transmission of cholera and consumption of water that has been contaminated by infected individuals. Infected individuals recover from infection at a rate γ, where 1/γ is the average length of the infectious period. Once the infected individual has recovered from infection, they move to the recovered compartment. Recovered individuals develop a temporary immunity to the disease after infection. They move from the recovered compartment to the susceptible compartment with rate µ, where 1/µ is the average length of immunity. Similar to Codeço [2001] and Koelle and Pascual [2004], birth and death are incorporated into the system indirectly through the waning of immunity; thus, instead of representing natural loss of immunity only, µ also represents the loss of immunity through the death of recovered individuals and birth of new susceptible individuals.
We model X t as an inhomogeneous Markov process [Taylor and Karlin, 1998] with infinitesimal rates   3   S, I, R   S+1, I, R-1  S, I-1, R+1 S-1, I+1, R From the current state (S, I, R), the system can transition to one of three new states. These new states correspond to a susceptible becoming infected, an infected recovering from infection, or a recovered individual losing immunity to infection and becoming susceptible. The parameter β is the infectious contact rate, α(t) is the time-varying environmental force of infection, γ is the recovery rate, and µ is the rate at which immunity is lost.
where X = (S, I, R) is the current state and X = (S , I , R ) is a new state. Because R t = N − S t − I t , we keep track of only susceptible and infected individuals, S t and I t .
This type of compartmental model is similar to other cholera models in the literature. The time-series SIRS model of Koelle and Pascual [2004] also includes the effects of both intrinsic factors (disease dynamics) and extrinsic factors (environment) on transmission. King et al. [2008] examine both a regular SIRS model and a two-path model to include asymptomatic infections, and use a time-varying transmission term that incorporates transmission via the environmental reservoir and direct person-to-person transmission, but does not allow for feedback from infected individuals into the environmental reservoir. The SIWR model of Tien and Earn [2010] and Eisenberg et al. [2013] allows for infections from both a water compartment (W) and direct transmission and considers the feedback created by infected individuals contaminating the water. To allow for the possibility of asymptomatic individuals, Longini et al. [2007] use a model with a compartment for asymptomatic infections; that model only considers direct transmission. Codeço [2001] uses an SIR model with no direct person-to-person transmission; infected individuals excrete directly into the environment and susceptible individuals are infected from exposure to contaminated water. Our SIRS model is not identical to any of the above models, but it borrows from them two important features: explicit modeling of disease transmission from either direct person-to-person transmission of cholera or consumption of water that has been contaminated by infected individuals and a time-varying environmental force of infection.

Hidden SIRS model
While the underlying dynamics of the disease are described by X t , these states are not directly observed. The number y t of infected individuals observed at each time point t is only a random fraction of the number of infected individuals. This fraction depends on both the number of infected individuals that are symptomatic and the fraction of symptomatic infected individuals that seek treatment and get reported (the reporting rate). Thus, y t i , the number of observed infections at time t i for observation i ∈ {0, 1, ..., n}, has a binomial distribution with size I t i , the number of infected individuals at time t i , and success probability ρ, the probability of infected individuals seeking treatment, so Given X t i , y t i is independent of the other observations and other hidden states.
We use a Bayesian framework to estimate the parameters of the hidden SIRS model, where the unobserved states X t are governed by the infinitesimal rates in Equation (1). The parameters that we want to estimate are β, γ, µ, ρ, and the k + 1 parameters that will be incorporated into α(t), the time-varying environmental force of infection. We assume α(t) = exp (α 0 + α 1 C 1 (t) + · · · + α k C k (t)), where C 1 (t), . . . , C k (t) denote the k time-varying environmental covariates.
We assume independent Poisson initial distributions for S t 0 and I t 0 , with means φ S and φ I . Thus Parameters that are constrained to be greater than zero, such as β, γ, µ, φ S , and φ I , are transformed to the log scale. A logit transformation is used for the probability ρ. We assume independent normal prior distributions on all of the transformed parameters, incorporating biological information into the priors where possible.
We are interested in the posterior distribution Pr(θ|y) ∝ Pr(y|θ)Pr(θ), where y = (y t 0 , . . . , y tn ), θ = (log(β), log(γ), log(µ), logit(ρ), α 0 , . . . , α k , log(φ S ), log(φ I )), and Here p(X t i |X t i−1 , θ) for i = 1, . . . , n are the transition probabilities of the continuous-time Markov chain (CTMC). However, this likelihood is intractable; there is no practical method to compute the finite time transition probabilities of the SIRS CTMC because the size of the state space of X t grows on the order of N 2 . For the same reason, summing over X with the forward-backward algorithm [Baum et al., 1970] is not feasible. To use Bayesian inference despite this likelihood intractability, we turn to a particle marginal Metropolis-Hastings (PMMH) algorithm.

Overview
The PMMH algorithm, introduced by Beaumont [2003] and studied in Andrieu and Roberts [2009] and Andrieu et al. [2010], constructs a Markov chain that targets the joint posterior distribution π(θ, X|y), where X is a set of auxiliary or hidden variables, and requires only an unbiased estimate of the likelihood. To construct this likelihood estimate, we use an SMC algorithm, also known as a bootstrap particle filter [Doucet et al., 2001]. The SMC algorithm sequentially estimates the likelihood using weighted particles; it requires the ability to propagate the unobserved data, X t , forward in time and the calculation of the probability of the observed data given the simulated unobserved data. For the hidden SIRS model, where ρ depends on the number of symptomatic infected individuals that seek treatment, as described in Section 3. Thus the probability of the observed data given the simulated unobserved data is given by Equation (2). To propagate the hidden variables forward in time, we first simulate initial states X t 0 = (S t 0 , I t 0 ) from Poisson distributions with means φ S and φ I . We then use properties of CTMCs to simulate the trajectories of the unobserved states.
Thus, the PMMH algorithm has two parts: an SMC algorithm, which is used to estimate the marginal likelihood of the data given a particular set of parameters, θ, and a Metropolis-Hastings step [Metropolis et al., 1953, Hastings, 1970, which uses the estimated likelihood in the acceptance ratio. At each step, a new θ * is proposed from the proposal distribution q(·|θ). An SMC algorithm is used to generate and weight K particle trajectories corresponding to the hidden state processes using the proposed parameter set θ * . A proposed X * t 0:n = X * t 0 , . . . , X * tn trajectory is sampled from the K particle trajectories based on the final particle weights of the SMC algorithm. The marginal likelihood is estimated by summing the 6 weights of the SMC algorithm, and the proposed θ * and X * t 0:n are accepted with probability equal to the familiar Metropolis-Hastings acceptance ratio.
To propagate the unobserved X t = (S t , I t , R t ) forward in time, we simulate from a cholera transmission model with a time-varying environmental force of infection. CTMCs which incorporate time-varying transition rates are inhomogeneous. The details of the discretelyobserved inhomogeneous CTMC simulations are now described.

Simulating inhomogeneous SIRS using tau-leaping
Gillespie developed two methods for exact stochastic simulation of trajectories with constant rates: the direct method [Gillespie, 1977] and the first reaction method [Gillespie, 1976]. Details of these methods are given in Appendix A. The exact algorithms work for small populations, but for large state spaces these methods require a prohibitively long computing time. This is a common problem in the chemical kinetics literature, where an approximate method called the tau-leaping algorithm originated [Gillespie, 2001, Cao et al., 2005. This method simulates CTMCs by jumping over a small amount of time τ and approximating the number of events that happen in this time using a series of Poisson distributions. As τ approaches zero, this approximation theoretically approaches the exact algorithm. The value of τ must be chosen such that the rates remain roughly constant over the period of time; this is referred to as the "leap condition".
Specifically, for our simulation, using the methods outlined in Cao et al. [2005], we define the rate functions h 1 (X t ) = (βI t + α(t)) S t , h 2 (X t ) = γI t , and h 3 (X t ) = µR t , corresponding to the infinitesimal rates of the CTMC. Then k 1 ∼ Poisson(h 1 (X t )τ ) represents the number of infections in time [t, t + τ ), k 2 ∼ Poisson(h 2 (X t )τ ) represents the number of recoveries in time [t, t + τ ), and k 3 ∼ Poisson(h 3 (X t )τ ) represents the number of people that become susceptible to infection in time [t, t+τ ). We make the assumption that the time-varying force of infection, α(t), remains constant each day. We define daily time intervals A i := [i, i + 1) for i ∈ {t 0 , t 0 + 1, . . . , t n − 1}, and α(t) = α A i for t ∈ A i . Using τ = 1 day, our rates now remain constant within each tau jump. To see if this value for τ is reasonable, we perform a simulation study; see Appendix A for details.

Metropolis-Hastings proposal for model parameters
Our implementation of the PMMH algorithm starts with a preliminary run, which consists of a burn-in run plus a secondary run, both using independent normal random walk proposal distributions for the parameters. From the secondary run, we calculate the approximate posterior covariance of the parameters and use it to construct the covariance of the multivariate normal random walk proposal distribution in the final run of the PMMH algorithm. In all runs, parameters are proposed and updated jointly.

Prediction
One of the main goals of this analysis is to be able to predict cholera outbreaks in advance using environmental predictors. To assess the predictive ability of our model, we estimate the parameters of the model using a training set of data and then predict future behavior of the epidemic process. We examine the posterior predictive distributions of cholera counts by simulating data forward in time under the time-varying SIRS model using the accepted parameter values explored by the particle MCMC algorithm and the accepted values of the hidden states S T and I T at the final observation time, t = T , of the training data. These hidden states are sampled in the PMMH algorithm by sampling the last set of particles using the last set of weights [Andrieu et al., 2010]. Under each set of parameters, we generate possible future hidden states and observed data, and we compare the posterior predictive distribution of observed cholera cases to the test data. In the analyses below, the PMMH output is always thinned to 500 iterations for prediction purposes by saving only every kth iteration, where k depends on the total number of iterations.

Simulation results
To test the PMMH algorithm on simulated infectious disease data, we generate data from a hidden SIRS model with a time-varying environmental force of infection. We then use our Bayesian framework to estimate the parameters of the simulated model and compare the posterior distributions of the parameters with the true values. To simulate endemic cholera where many people have been previously infected, we start with a population size of N = 10000 and assume independent Poisson initial distributions for S t 0 and I t 0 , with means φ S = 2100 and φ I = 15. The other parameters are set at β = 1.25 × 10 −5 , γ = 0.1, and µ = 0.0009. All rates are measured in the number of events per day. The average length of the infectious period, 1/γ, is set to be 10 days, and the average length of immunity, 1/µ, is set to be about 3 years. Parameter values are chosen such that the simulated data are similar to the data collected from Mathbaria, Bangladesh. We use the daily time intervals A i := [i, i + 1) for i ∈ {t 0 , t 0 + 1, . . . , t n − 1}, as in Section 4.2, and define α(t) = α A i for t ∈ A i where α A i = exp [α 0 + α 1 sin (2πi/365)] . The intercept α 0 and the amplitude α 1 are parameters to be estimated. The frequency of the sine function is set to mimic the annual peak seen in the environmental data collected from Bangladesh. For the simulations we set α 0 = −7 and α 1 = 3.5. Using the modified Gillespie algorithm described in Appendix A, we simulate the (S t , I t ) chain given in the left plot of Figure 2. The observed number of infections y t ∼ Binomial(I t , ρ), where ρ = 0.015 and is treated as an unknown parameter.
We simulate three years of training data because this is approximately how long the data collection phases last in our data from Bangladesh [Huq et al., 2005]. Therefore we do not attempt to estimate the loss of immunity rate µ since it is on the scale of three years. Also, there is not enough information in the data to estimate the means of the Poisson initial  Figure 2: Plots of simulated hidden states (counts of susceptible, S t , and infected, I t , individuals) and the observed data (number of observed infections = y t ∼ Binomial(I t , ρ) ) plotted over time, t. Simulation with seasonally varying α(t) generates data with seasonal epidemic peaks. The dashed vertical black line represents the first cut off between the training sets and the test data. Data before the line are used to estimate parameters, and we use those estimates to predict the data after the line. Other data cut offs are shown in Figure 3 distributions, φ S and φ I , since estimation of these parameters is only informed by the very beginning of the observed data. We set these parameters to different values and compare parameter estimation and prediction between models with parameter assumptions which differ from the truth. We also assume that we know the population size, N = 10000.
We assume normal prior distributions on all of the parameters, with means and standard deviations chosen such that the mass of each prior distribution is not centered at true value of the parameter in this simulation setting. We use relatively uninformative, diffuse priors for log(β), α 0 , and α 1 , centered at log(1.25 × 10 −4 ), −8, and 0, respectively, and with standard deviations of 5. The prior distribution for logit(ρ) is centered at logit(0.03) and has a standard deviation of 2. For log(γ), the prior is centered at log(0.1) with a relatively small standard deviation of 0.09, since this value is well studied for cholera. Thus, a priori 1/γ falls between 8.4 to 11.9 days with probability 0.95.
Using these data, the PMMH algorithm starts with a burn-in run of 10000 iterations, a secondary run of 10000 iterations, and a final run of 50000 iterations. To thin the chains, we save only every 10th iteration. We use K = 100 particles in the SMC algorithm. We compare results from models with different assumptions on the values of φ S and φ I : assumed φ S /N and φ I /N are above the true values (0.31 and 0.003), at the true values (0.21 and 0.0015), below the true values (0.11 and 0.00075), or further below the true values (0.055 and 0.000375). Marginal posterior distributions for the parameters of the SIRS model from the final runs of these PMMH algorithms are in Appendix B. The posterior distributions are similar, regardless of assumed values for φ S and φ I . Trace plots, auto-correlation plots, bivariate scatterplots, and effective sample sizes for the posterior samples under the situation in which the true values of φ S and φ I were assumed are also given in Appendix B. We report β × N and ρ × N , since in sensitivity analyses we found these to be robust to assumptions about the total population size N . From the posterior distributions, it is clear that the algorithm is providing good estimates of the true parameter values, though estimates of the parameter ρ × N are slightly different than the truth, especially when φ S and φ I are not set at the true values. To test the predictive ability of the model, we use multiple cut off times to separate our simulated data into staggered training sets and test sets. The simulated observed data are shown in the right plot of Figure 2. For each cut off time, parameters were drawn from the posterior distribution based on the training data. These parameter values were then used to simulate possible realizations of reported infections after the training data until the next cut off, 28 days later. The distributions of these predicted reported cases are shown in the top plot of Figure 3. The test data are denoted by the purple diamonds, connected by straight lines to help visualize ups and downs in the case counts. Case counts are observed once every 14 days. On each observation day, the colored bar represents the distribution of predicted counts for that day. As desired, the posterior predictive distribution shifts its mass as time progresses to follow the case counts in the test data. The plot of the predicted hidden states in the bottom row of Figure 3 also shows that our model is capturing the formation and decline of the epidemic peak well, as seen in the trajectory of the predicted fraction of infected individuals. This plot illustrates the interplay of the hidden states of the underlying compartmental model. During an epidemic, the fraction of susceptibles decreases while the fraction of infected individuals quickly increases. Afterwards, the fraction of infected individuals drops and the pool of susceptibles slowly begins to increase as both immunity is lost and more susceptible individuals are born.

Prediction results
These predictions were made under the assumption that φ S and φ I are set to the true values. To test sensitivity to these assumptions, we compare predictions made from models that assume other values; these are shown in Appendix D. Predicted distributions are similar for all values of φ S and φ I .
6 Using cholera incidence data and covariates from Mathbaria, Bangladesh Huq et al. [2005] found that water temperature (WT) and water depth (WD) in some water bodies had a significant lagged relationship with cholera incidence. Therefore, we use these covariates and cholera incidence data from Mathbaria, Bangladesh collected between April 2004 to September 2007 and again from October 2010 to July 2013. During these time periods, cholera incidence data were collected over a period of three days approximately every two weeks. Environmental data were also collected approximately every two weeks from six water bodies. To get a smooth summary of the covariates using data from all water bodies, we fit a cubic spline to the covariate values. We then slightly modify our environmental force of infection to allow for a lagged covariate effect. Let κ denote the length of the lag. We consider the daily time intervals A i := [i, i + 1) for i ∈ {t 0 , t 0 + deviation. We consider and compare results from models assuming three different lags: κ = 14, κ = 18, and κ = 21. Predictions from all three models are similar, so we report only results from the model assuming κ = 21, in order to receive the earliest warning of upcoming epidemics; see Appendix C for details and prediction comparisons. The smoothed, standardized, 21 day lagged covariates and cholera incidence data are shown in Figure 4. Since there are only about six years of data, estimating the loss of immunity rate µ is infeasible. Thus, we set µ = 0.0009 so that 1/µ is 3 years [Sack et al., 2004]. Also, the population size N , which quantifies the size catchment area for the medical center, is assumed to be 10000 for computational convenience. We do not know the true value of N , but 10000 is a reasonable estimate and is small enough that simulations run quickly. We studied sensitivity to these assumptions by setting both µ and N to different values, obtaining similar results. We also again set φ S and φ I to various values and the results were insensitive. See Appendix D for details.
In these analyses, we use relatively uninformative, diffuse normal prior distributions on the time-varying environmental covariates α 1 and α 2 , centered at 0 and with standard deviations of 5. The diffuse normal prior distributions on the transformed parameter values log(β) and α 0 are centered at log(1.25 × 10 −7 ) and −8, respectively, with standard deviations of 5. We know that the average infectious period for cholera, 1/γ, should be between 8 and 12 days. Thus, the transformed parameter log(γ) is given a normal prior distribution with mean log(0.1) and standard deviation 0.09 to give 0.95 prior probability of 1/γ falling within the interval (8, 12). We also know that ρ should be very close to zero, since only a small proportion of cholera infections are symptomatic and a smaller proportion will be treated at the health complex [Sack et al., 2003]. Thus, the transformed parameter logit(ρ) is given a normal prior distribution with mean logit(0.0008) and standard deviation equal to 2, to give 0.95 prior probability of ρ falling within the interval (1.6 × 10 −5 , 0.04).
We run the PMMH algorithm with a burn-in run of 30000 iterations, a secondary run of 20000 iterations, and a final run of 400000 iterations. We again save only every 10th iteration and use K = 100 particles in the SMC algorithm. Posterior medians and 95% Bayesian credible intervals for the parameters β × N , γ, α 0 , α 1 , α 2 , and ρ × N generated by the final run of the PMMH algorithm are given in Table 1. We report β × N and ρ × N since we found these parameter estimates to be robust to changes in the population size N during sensitivity analyses. For more details, see Appendix D. The credible intervals for α 1 and α 2 do not include zero, so both water depth and water temperature have a significant relationship with the force of infection. Decreasing water depth increases the force of infection, likely due to the higher concentration and resulting proliferation of V. cholerae in the environment; increasing water temperature increases the force of infection [Huq et al., 2005].
The basic reproductive number, R 0 , is the average number of secondary cases caused by a typical infected individual in a completely susceptible population [Diekmann et al., 1990]. We report (β × N )/γ, the part of the reproductive number that is related to the number of infected individuals in the population under our model assumptions. Our estimate of 4.35 is fairly large; it is very similar to the reproductive number of 5 (sd=3.3) estimated by Longini et al. [2007] using data from Matlab, Bangladesh. However, the 95% credible interval is wide, with the lower end being approximately 1. Moreover, posterior median values for α(t) range from 0.00003 to 0.38, while posterior median values for βI t only range from 0 to 0.03, suggesting that the epidemic peaks in our model are driven mostly by the environmental force of infection. See Appendix F for more details. However, the infectious contact rate is not zero and is not negligible compared to the environmental force of infection.

Prediction Results
For the data collected from Mathbaria, we begin prediction at multiple points around the time of the two epidemic peaks that occur in 2012 and 2013. Figure 4 shows the full cholera data with smoothed and standardized covariates. Figure 5 shows the posterior predictive distribution of observed cholera cases (top row) and hidden states from the time-varying SIRS model (bottom row). Parameters used to simulate the SIRS forward in time have been sampled using the PMMH algorithm applied to the training data, with data being cut off at different points during the 2012 and 2013 epidemic peaks. From each of these cut offs, parameter values are then used to simulate possible realizations of the test data. Predictions are run until the next cut off point, with cut off points chosen based on the length of the lag κ. Realistically, at time t we have covariate information to use for prediction only until time t + κ, where κ is the covariate lag. Since the smallest lag considered is 14 days, we make only 14 day ahead predictions where possible to mimic a realistic prediction set up. Due to the sparse sampling between epidemic peaks (June 2012 to February 2013), we use longer prediction intervals for these cut-offs than would be possible in real time data analysis in order to evaluate our model predictions.
In the top row of Figure 5, the coloring of the bars again represents the distribution of predicted cases. Between the two peaks of case counts (June 2012 to February 2013), the frequency of predicted zero counts is very high, so we conclude that the model is doing well with respect to predicting the lack of an epidemic. During the epidemics, the distribution of the counts shifts its mass away from zero. The plot in the bottom row of Figure  , as observed in the test data for those years. The predicted fraction of infected people in the population increases before an increase can be seen in the case counts, which could allow for early warning of an epidemic.
We also use a quasi-Poisson regression model similar to the one used by Huq et al. [2005] to predict the mean number of cholera cases (Appendix E). Although the quasi-Poisson model predicts reasonably well the timing of epidemic peaks, it appears to overestimate the duration of the outbreaks. The predicted means under both the quasi-Poisson and SIRS models most likely underestimate the true mean of the observed counts, with the quasi-Poisson model performing slightly better. However, the SIRS predicted fraction of infected individuals -a hidden variable in the SIRS model -provides a more detailed picture of how cholera affects a population. By providing not only accurate prediction of the time of epidemic peaks, but also the predicted fraction of the population that is infected, the SIRS model predictions could be used for efficient resource allocation to treat infected individuals. See Appendix E for additional details.

15
We use a Bayesian framework to fit a nonlinear dynamic model for cholera transmission in Bangladesh which incorporates environmental covariate effects. We demonstrate these techniques on simulated data from a hidden SIRS model with a time-varying environmental force of infection, and the results show that we are recovering well the true parameter values. We also estimate the effect of two environmental covariates on cholera case counts in Mathbaria, Bangladesh while accounting for infectious disease dynamics, and we test the predictive ability of our model. Overall, the prediction results look promising. Based on data collected, the predicted hidden states show a noticeable increase in the fraction of infected individuals weeks before the observed number of cholera cases increases, which could allow for early notification of an epidemic and timely allocation of resources. The predicted hidden states show that the fraction of infected individuals in the population decreases greatly between epidemics, supporting the hypothesis that the environmental force of infection triggers outbreaks. Estimates of βI t are low, but not negligible, compared to estimates of α(t), suggesting that most of the transmission is coming from environmental sources.
Computational efficiency is an important factor in determining the usefulness of this approach in the field. We have written an R package which implements the PMMH algorithm for our hidden SIRS model, available at https://github.com/vnminin/bayessir. The computationally expensive portions of the PMMH code are primarily written in C++ to optimize performance, using Rcpp to integrate C++ and R [Eddelbuettel andFrançois, 2011, Eddelbuettel, 2013]; however there is still room for improvement. Running 400000 iterations of the PMMH algorithm on the six years of data from Mathbaria takes 2.5 days on a 4.3 GHz i7 processor. Since we can predict three weeks into the future using a 21 day covariate lag, we do not think timing is a big limitation for using our model predictions in practice.
Plots of residuals over time, shown in Appendix C, show that we are modeling well case counts between the epidemic peaks but not the epidemic peaks themselves, either due to missing the timing of the epidemic peak or the latent states not being modeled accurately. This possible model misspecification might be fixed by including more covariates, using different lags, or modifying the SIRS model. Also, we assume a constant reporting rate, ρ, rather than using a time-varying ρ t [Finkenstädt and Grenfell, 2000]. With better quality data we might be able to allow for a reporting rate that varies over time; we will try to address these model refinements in future analyses.
In the future, we will extend this analysis to allow for variable selection over a large number of covariates. This will allow us to include many covariates at many different lags and incorporate information from all of the water bodies in a way that does not involve averaging.
In the current PMMH framework, choosing an optimal proposal distribution to explore a much larger parameter space would be difficult. We want to include a way of automatically selecting covariates or shrinking irrelevant covariate effects to zero with sparsity inducing priors. The particle Gibbs sampler, introduced by Andrieu et al. [2010], would allow for such extensions. Approximate Bayesian computation is also an option for further model development [McKinley et al., 2009]. In addition, the available data consist of observations from multiple thanas during the same time period. Future analyses will look into sharing information across space and time and accounting for correlations between thanas. Another challenging future direction involves exploring models which incorporate a feedback loop from infected individuals back into the environment to capture the effect of infected individuals excreting V. cholerae into the environment. To accomplish this, we could add a water compartment to our SIRS model that quantifies the concentration of V. cholerae in the environment, similar to the model of Tien and Earn [2010]. However, adding an additional latent state leads to identifiability problems, even with fully observed data [Eisenberg et al., 2013], so such an extension will require rigorous testing and fine tuning.

Supplementary Materials
Appendix A: Simulation Details PMMH pseudocode The following exposition of the algorithm follows closely the pseudocode of Andrieu et al. [2010] and Wilkinson [2011].
Thus we have a simple, sequential, likelihood-free algorithm which generates an unbiased estimate of the marginal likelihood, p(y|θ(0)). A X t 0:n (0) trajectory is sampled from the K trajectories (X k t 0:n , for k = 1, . . . , K) based on the final set of particle weights, W (X k tn ).

Simulating homogeneous SIRS
Gillespie's direct method [Gillespie, 1977] simulates the time to the next event and then determines which event happens at that time. The first reaction method [Gillespie, 1976] calculates the time to the next reaction for each of the possible events, and the minimum time to next reaction determines the next step of the chain.
Using the direct method, we can think of our continuous-time Markov chain (CTMC) as a chemical system with three different reactions. These reactions and their rate functions are given by the infinitesimal rates Thus the three reactions have the rate functions h 1 (X t ) = (βI t + α)S t , h 2 (X t ) = γI t , and h 3 (X t ) = µR t , corresponding to the infinitesimal rates of the CTMC. Then the time to the next reaction, τ , has an exponential distribution with rate λ = h 1 (X t ) + h 2 (X t ) + h 3 (X t ), and the kth reaction occurs with probability h k (X t )/λ, for k = {1, 2, 3}.
The first reaction method instead simulates the time τ k that the kth reaction happens for k = {1, 2, 3}, given no other reactions happen in that time. Then the time to the next reaction τ = min k (τ k ), and the reaction with the reaction time equal to τ is the event that happens.
Both the direct method and the first reaction method work only for homogeneous Markov chains. If we want to assume that the additional force of infection, α, varies over time, the associated Markov chain is inhomogeneous and we must account for the fact that the transition rate could change before the next reaction occurs. Gibson and Bruck [2000] introduce the next reaction method, an efficient exact algorithm to simulate stochastic chemical systems. They extend this next reaction method to include timedependent rates and non-Markov processes. Anderson [2007] deviates from these methods a bit, using Poisson processes to represent the reaction times, with time to next reaction given by integrated rate functions. This leads to a more efficient modified next reaction method which they extend to systems with more complicated reaction dynamics.

Simulating inhomogeneous SIRS
Using the methods described by Gibson and Bruck [2000] and Anderson [2007], to incorporate a time-varying force of infection into the SIRS model we must integrate over the rate function h 1 (X t , s) = (βI t + α(s)) S t . Thus, to find the time τ 1 that the first reaction happens, given no other reactions happen in that time, we generate u ∼ Uniform(0, 1) and solve τ 1 t h 1 (X t , s)ds = ln(1/u) for τ 1 . Since the other two reactions have no time-varying parameters, we can solve for τ 2 and τ 3 , the reaction times of the second and third reactions, using the methods of the previous section. Then we can continue, using the first reaction method to simulate the process.
We simplify this approach by assuming that the time-varying force of infection, α(t), remains constant each day. We define daily time intervals A i := [i, i+1) for i ∈ {t 0 , t 0 +1, . . . , t n −1}, and α(t) = α A i for t ∈ A i . Then we can take advantage of the memoryless property of exponentials and propagate the chain forward in daily increments. Thus, we use the direct method, but when the time to next event exceeds the right end point of the current interval A i , we restart CTMC simulation from the beginning of the interval A i+1 using α A i+1 in the waiting time distribution rate λ(α This modified Gillespie algorithm is depicted and detailed in Figure

Selecting Tau
Unchecked, tau-leaping can lead to negative population sizes in a compartment if the compartment has a low number of individuals. To avoid this, we use a simplified version of the modified tau-leaping algorithm presented by Cao et al. [2005]. If the population of a Figure A-1: Depiction of the modified Gillespie algorithm. We assume the environmental force of infection, α(t), is a step function which changes daily. Daily time intervals are denoted by Starting at time t = 1, the waiting time to the next event, τ , has an exponential distribution with rate λ(α In the depiction, τ = 0.5. The simulated waiting time plus the current time, t * = t + τ , remains in the interval A 1 , so we use t * as the next time in our CTMC and propagate X t forward at that time using Gillespie's direct method. Since we are still in the interval A 1 , we again simulate the time to the next event as an exponential random variable with rate λ(α A 1 ) = h 1 (X t * , α A 1 )+h 2 (X t * )+h 3 (X t * ). In this iteration, the waiting time plus the current time, t * + τ , exceeds the boundary of the interval A 1 , so we discard this simulated waiting time τ . Using the memoryless property of exponentials, we restart our simulation from the beginning of the interval A 2 using the new α(t) value, α A 2 . We continue in this manner until we have simulated the Markov process X t up to time t n . compartment is lower than some prespecified critical size, a single step algorithm (like the Gillespie algorithm) is used until the population gets above that critical size. If the size of the compartment is not critically low but the current value of τ still produces a negative population, we reject that simulation and try again with a smaller τ (reduced by a factor of 1/2). The subsequent value of τ is picked based on how long the current daily time-varying force of infection remains constant. We choose a value of τ that simulates what happens during the remainder of the day, until the value of the transition rate changes. This modified tau-leaping algorithm is depicted and detailed in Figure A-2.
For our simulations, we have chosen τ = 1 day; we perform a simulation study to see if this value for τ is reasonable. Using the posterior estimates of the parameters, we simulate data forward in time 5000 times using both the modified Gillespie algorithm and the modified tau-leaping algorithm. We simulate data over the entire epidemic curve to see how the comparison changes for varying values of α(t). Figure A-3 shows estimates of the median and 95% intervals for the simulated values. The Monte Carlo standard error is very small for all estimates. For the numbers of susceptible individuals, the estimates under Gillespie and tau-leaping are almost identical over the entire epidemic. For the numbers of infected, the values are very close except at the epidemic peaks. However, the differences are very small. We conclude that for our application τ = 1 day is a good compromise between computational efficiency and accuracy.

Binomial tau-leaping
Another solution to the negative population size problem is to use Binomial tau-leaping [Chatterjee et al., 2005, Tian andBurrage, 2004], which further approximates k j as a binomial random variable with mean h j (X t )τ and upper limit chosen such that k j cannot be large enough to simulate a negative population. We opt instead to use the simplified version of the modified tau-leaping algorithm.
Negative population Negative population SSA τ=1−0.2=0.8 Figure A-2: Depiction of the modified tau-leaping algorithm. We assume the environmental force of infection, α(t), is a step function which changes daily. Daily time intervals are denoted by A i := [i, i + 1) for i ∈ {t 0 , t 0 + 1, . . . , t n − 1}, so α(t) = α A i for t ∈ A i . As a default, we use τ = 1 day. Starting at time t = 1, we simulate the changes in compartment populations over the interval t ∈ [1, 2). At time t = 2, we again use τ = 1 day to simulate the changes over the interval t ∈ [2, 3). This value of τ produces a negative population so we reject that simulation and try again with a smaller τ (reduced by a factor of 1/2). The next value of τ is then calculated based on how long the current daily time-varying force of infection remains constant, so τ = 0.5. At time t = 3, the population of a compartment is lower than some prespecified critical size, so a single step algorithm (SSA), in our case the Gillespie algorithm, is used until the population gets above that critical size. Once the compartment populations are all above the critical size again, at time t = 3.2, the subsequent value of τ is again picked based on how long the current daily time-varying force of infection remains constant, so τ = 0.8.  Table B-1 and the parameter estimates from the chains that started at these initial values are given in the top third of Table B-2. Credible intervals for β × N vary slightly for different initial values; this is likely due to a heavy tail in the posterior distribution that is not yet explored in the run initializing from the second set of starting value and is most often explored in the run initializing at the third set of starting values. If we obtained larger samples from the posterior, the credible intervals would be more similar.

Appendix C: Model fit
To select a lag for the environmental covariates in the Mathbaria analysis, we compare prediction results from models assuming three different lags: κ = 14, κ = 18, and κ = 21. These are shown in Figures C-1 and C-2. The predictive distributions of the hidden states look similar across lags, so we use the 21 day lag model in order to predict an upcoming epidemic furthest in advance. With a three week lag, we would be able to make predictions three weeks in advance.   q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q qq q q

Appendix F: Routes of transmission
An important question in cholera modeling is: what is the relative contribution of different routes of transmission at different points of the epidemic? We hypothesized that environmental forces trigger the seasonal cholera epidemics and that infectious contact between susceptible and infected individuals drives the epidemics. To examine this possible dynamic, we compare the forces of infection from the environment, α(t), to that from infected individuals, β × I t , over time. Values are computed by sampling 5000 sets of parameter values from the posterior. For each set of parameters, we generate data using our hidden SIRS model. Figure F-1 shows median and 95% quantiles for α(t) vs β × I t plotted over time. The median values of α(t) are almost always higher than values of β × I t ; the only time it is not (in early 2011) is most likely caused by model misspecification, specifically in setting φ I to the wrong value for this phase. This model misspecification is similarly reflected in Figure  C-3. This supports the hypothesis that the epidemics are driven by the environmental force of infection.