extfieldsofstudy:
Mathematics
Geography
Medicine
Computer Science

text:
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.


s2fieldsofstudy: Economics Medicine

text:
Maternal and Neonatal Mortality in South-West Ethiopia: Estimates and Socio-Economic Inequality

Introduction Ethiopia has achieved the fourth Millennium Development Goal by reducing under 5 mortality. Nevertheless, there are challenges in reducing maternal and neonatal mortality. The aim of this study was to estimate maternal and neonatal mortality and the socio-economic inequalities of these mortalities in rural south-west Ethiopia. Methods We visited and enumerated all households but collected data from those that reported pregnancy and birth outcomes in the last five years in 15 of the 30 rural kebeles in Bonke woreda, Gamo Gofa, south-west Ethiopia. The primary outcomes were maternal and neonatal mortality and a secondary outcome was the rate of institutional delivery. Results We found 11,762 births in 6572 households; 11,536 live and 226 stillbirths. There were 49 maternal deaths; yielding a maternal mortality ratio of 425 per 100,000 live births (95% CI:318–556). The poorest households had greater MMR compared to richest (550 vs 239 per 100,000 live births). However, the socio-economic factors examined did not have statistically significant association with maternal mortality. There were 308 neonatal deaths; resulting in a neonatal mortality ratio of 27 per 1000 live births (95% CI: 24–30). Neonatal mortality was greater in households in the poorest quartile compared to the richest; adjusted OR (AOR): 2.62 (95% CI: 1.65–4.15), headed by illiterates compared to better educated; AOR: 3.54 (95% CI: 1.11–11.30), far from road (≥6 km) compared to within 5 km; AOR: 2.40 (95% CI: 1.56–3.69), that had three or more births in five years compared to two or less; AOR: 3.22 (95% CI: 2.45–4.22). Households with maternal mortality had an increased risk of stillbirths; OR: 11.6 (95% CI: 6.00–22.7), and neonatal deaths; OR: 7.2 (95% CI: 3.6–14.3). Institutional delivery was only 3.7%. Conclusion High mortality with socio-economic inequality and low institutional delivery highlight the importance of strengthening obstetric interventions in rural south-west Ethiopia.

Introduction
In 2010, there were 287,000 maternal deaths in the world from causes of pregnancy related complications, which is 50% down from the 1990 baseline [1]. Every year, four million newborns die during the neonatal period (28 days after birth) [2], while 3.2 million pregnancies end with stillbirths [3]. Among these deaths, 99% of the neonatal deaths and stillbirths as well as 98% of maternal deaths, occur in low-and middle-income countries. Moreover, 40% of global, and 29% of African child mortality is caused by neonatal deaths [4]. The Millennium Development Goals (MDG-4 and 5) aim to reduce child mortality by two-thirds, and maternal mortality by three-quarters, between 1990 and 2015. High maternal and neonatal deaths and stillbirths often occur because of inadequate care during pregnancy and childbirth. Accordingly, these deaths are considered sensitive indicators of the quality of a healthcare system in an area [5]. High neonatal deaths and stillbirths are often related to maternal health; complications during pregnancy, childbirth, and the post-natal period increase the risk of death for both the baby and the mother [6].
Ethiopia has achieved MDG-4 by reducing two-thirds of under 5 mortality per 1000 live births from 204 in 1990 to 68 in 2012 [7]. Ethiopia's achievement has been attributed to the country's community-based health promotion and disease prevention programme through a health extension package [8]. However in terms of neonatal mortality, maternal mortality, and stillbirth reduction, challenges still remain. Ethiopia is among 10 countries that contributed to two-thirds of global neonatal deaths in 2005 [2] and stillbirths in 2009 [9]. UN agencies (WHO and UNFPA) and World Bank data also indicate a decline in the maternal mortality ratio (MMR) per 100,000 live births in Ethiopia from 950 in 1990 to 350 in 2011 [1]. Nevertheless, there are often controversies over these estimates. For example, the 2011 national Demographic and Health Survey (DHS) reported an MMR of 676 per 100,000 live births, which is close to two-folds higher compared to the UN estimate [10]. According to the 2007 National Census, over 84% of the estimated 90 million people in Ethiopia live in rural areas with limited access to quality health care [11]. In 2008, a nationwide study showed that 7% of all deliveries took place at health institutions while only 3% in facilities that could provide comprehensive essential obstetric care [12].
In places such as rural Ethiopia, where there is no birth and death registration and the majority of births and deaths take place outside of health institutions, measuring maternal and neonatal mortality is difficult [13]. In some places, methods such as demographic surveillance systems help to find and measure maternal and neonatal mortalities [14]. Unfortunately, this alternative method does not exist in the Gamo Gofa province of Ethiopia. Maternal mortality can also be estimated through lowcost innovative options of the sisterhood method, which asks adult siblings about their sisters' death related to pregnancy or childbirth during reproductive age [15]. Nevertheless, results from the sisterhood method refer to many years before the survey and may not show the current magnitude of the problem [16]. Findings from well-planned household surveys that use large samples in high fertility and high mortality areas can be useful in providing real-time data to motivate actions [17]. The aim of this study was to estimate maternal and neonatal mortality, the stillbirth rate, the institutional delivery rate, and household risk factors associated with these mortality outcomes in rural south-west Ethiopia.

Ethics statement
The Ethical Review Committee for the Health Research of Southern Nations Nationalities and Peoples' Regional State (SNNPRS) Health Bureau in Ethiopia, and the Regional Committee for Health Research Ethics of North Norway (REK Nord) approved the study. We obtained informed verbal consent from all respondents and the response was recorded on the questionnaire as ''accepted'' or ''declined'' to participate. Almost all approached households were willing to be interviewed, and written consent was not considered because a large number of the respondents were illiterates. The study involved only interview and the ethics committee approved the verbal consent procedure. Additionally, minors were not included in this study.

Definitions
Verbal autopsy for maternal deaths. A method for finding out the medical causes of death and ascertaining the factors that may have contributed to the death in women who died outside of a medical facility. It consists of interviewing people (family members, neighbours, traditional birth attendants) who had knowledge about the events leading to the death [18].
Neonatal mortality. A death within 28 days of an alive born baby.
Maternal mortality. A death of a woman while pregnant, in labour, or within 42 days of the termination of pregnancy, irrespective of the duration and site of the pregnancy, from any cause related to or aggravated by the pregnancy or its management, but not from accidental or incidental causes (ICD-10) [19].
Stillbirth. A birth of a dead fetus after 28 weeks of gestation. We did not use the baby-weight criteria of classifying stillbirths, as it was not possible to measure weight in the rural area.
Household. A person or a group of people living in a room or rooms and sharing common things together. In cases of polygamy (more than one wife for a man, we considered each wife as a separate household as they culturally have separate houses. Maternal mortality ratio (MMR). Is the number of maternal deaths in a population during a given time period per 100,000 live births during the same period.
Neonatal mortality ratio (NMR). Is the number of newborn deaths (within 0-28 days) in a population per 1,000 live births in the same population.
Stillbirth rate. Is the number of births of dead fetuses after 28 weeks of gestation per 1,000 births.

Study area
We conducted this study in Bonke, one of the 15 woredas in the Gamo Gofa zone in south-west Ethiopia. In 2010, the woreda had a population of 173,240 people [11]. A kebele is the lowest administrative structure with 5,000 to 7,000 residents in the Ethiopian government system. Bonke has 31 kebeles and one of these, the administrative centre, has a town status with a population of 6,347 people in 2007. Table 1 shows the profile of the 15 kebeles included in this study. Bonke is 618 km from Addis Ababa, and 68 km from Arba Minch (zonal capital) where the nearest hospital is situated. Nevertheless, over half of the remote areas of Bonke are more than 100 km (20 hours walking distance) away from the hospital in Arba Minch. An estimated three-fourths of the population also live in villages far from the motorable road ($6 km). The only road to the woreda is the road from Arba Minch to Kamba, which crosses parts of Bonke. Overflowing rivers during the rainy season often interrupt the road. So, often people have to carry critical patients or use transport animals such as horses and mules to go to the hospital.
Health care is provided by a health centre in the town (Geresse) as well as three other rural health centres. There are no medical doctors working in the woreda; a few health officers (people with a bachelor's degree in medical training), nurses, and midwives staff the health centres. In Bonke, there is no access to lifesaving comprehensive essential obstetric care that can provide caesarean sections, blood transfusions, and effective care to sick and low birth-weight newborns. This study was part of an implementation project to reduce maternal mortality in Gamo Gofa. The project trains health officers in emergency obstetric services, community health workers in identifying and referring high-risk mothers, in addition to equipping health centres and hospitals with essential instruments.

Study design and period
The study was a cross-sectional household survey with a fiveyear recall of events prior to the data collection. We collected the data in February 2011 from all households that had births and pregnancy outcomes between January 2006 and December 2010. We purposely selected January 2006 as the starting reference period for the recall because it was the immediate period after the 2005 National Election, an event well known to all respondents.

Sampling
We based our sample size calculation on the assumptions of a crude national birth rate of 35 per 1,000 population, and a neonatal mortality ratio of 35 per 1,000 live births [20]. We aimed to detect a minimum of 350 neonatal deaths to make empirical estimates and assess the household risk factors associated with neonatal deaths. To find 350 neonatal deaths we needed to find at least 10,000 live births in five years (an average of 2,000 per year) within the population. With a fertility rate of 35 births per 1,000 people, 2,000 live births per year could be obtained from an estimated population of 57,143. Assuming a constant birth rate over the five years, we projected the population of 57,143 in 2006 to be 67,244 in 2010 (half the rural population in Bonke).
We used OpenEpi, open source calculator (www.openepi.com), and calculated the minimum sample needed based on information that three-fourths of the study households resided far from the motorable road ($6 km), thereby expecting a neonatal mortality prevalence twice that among households far from the road. We used a statistical power of 80%, a 95% confidence interval, and an assumption that 4% of households far from motorable road could experience neonatal mortality to calculate the number of households needed for the study. This provided 5,187 households with expected births in the five years before the survey. On average, we expected two births per household over five years yielding 10,374 births. The neonatal mortality rate from the estimated number of households was also assumed to provide enough power to detect other risk factors (wealth, education, nonspaced births).
We also assumed that a number of maternal deaths among the estimated 10,374 births would give an optimum MMR estimate, and a similar assumption was applied for stillbirths. Taking into account a potential 10% of non-responders, we decided to study 50% of the rural population in the Bonke woreda and we randomly selected 15 of the 30 rural kebeles in Bonke. Data collectors visited all of the households in the selected 15 kebeles asking about any pregnancy and birth outcomes (abortion, alive and stillbirths, neonatal and maternal deaths) in the households over the previous five years. Enumerators noted the number of households in each kebele, collecting data from the households that had pregnancy and birth outcomes during the stated time period. As projected from the Ethiopian 2007 census [11], the selected kebeles had a population of 78,181 people in 2010.

Variables and questions
The primary outcome variables were maternal and neonatal mortality and stillbirth rates while a secondary outcome variable was rate of skilled delivery service utilization. We also used the following household predictor variables: household wealth (assets index), educational level of heads of households, the number of births in a household over five years, and household distance to a motorable road.
Survey questions included: where each delivery took place (home, health post, health centre, and hospital) and who attended the birth (family member, community health extension worker, skilled health professional-doctor, midwife). We also asked about the place where the maternal death occurred (at home or within a health facility) as well as what had happened to the fetus (stillbirth, neonatal death, or alive at the time of data collection).
Respondents were asked about whether the newborn was alive or dead at the end of the 4th week after delivery. If the response was ''dead'', then we asked about the timing of death (in weeks) in relation to the birth. However, we did not investigate the causes of deaths for neonatal deaths and stillbirths, assuming it would be difficult for rural respondents to answer it properly. In the households that had deaths of women -in reproductive age (15-49 years), we used questions modified from the WHO manual for verbal autopsy for maternal death to investigate the causes of deaths [18].
The questions included: whether the mother was pregnant, in the process of giving birth, or in postnatal period after birth, what main medical condition or symptom was associated with her death, what assistance she received, and from whom she received help. A nurse decided on pre-coded choices of the major causes of maternal deaths (bleeding, prolonged labour, fever and convulsions, including the option of ''others'') based on quick algorithmic analysis of information provided by the respondents. Sensitive questions related to abortion deaths were placed at the end of the interview to minimize the intentional hiding of information. We collected information on the estimated walking distance (in hours) from each house to the nearest health centre, the nearest motorable road, and the nearest hospital. Based on the local experience of one hour of walking time per 5 km for an average person, we converted the walking distance into kilometers.

Data collection
We recruited 15 natives from the respective study villages who had completed the 12 th grade for data collection. The purpose of selecting data collectors from their respective kebele of data collection was to reduce the potential recall bias by the respondents. Data collectors are aware of many vital events in the villages they collected data by living and participating in social events such as birth celebrations, mourning rituals, and burials at the time of the deaths. Five diploma graduates who had a thorough knowledge of the culture and language of the area supervised the data collection.
The data collectors were trained for two days on pre-testing field interviews, translating the questions from ''Amharic'' (the official Ethiopian state language) to ''Gamotho'' (the language of the ethnic ''Gamo'' community) and how to introduce the simplified verbal autopsy questions. Depending whoever was present at home during the visit, the respondent was the father or mother for a recently deceased newborn. In cases of death of a married woman, we interviewed a husband while in the absence of a husband, an adult relative or an adult child of the deceased was interviewed. For those who were unmarried, we asked parents or siblings. If the respondents were not present at home during the first visit, the data collectors re-visited the next day in the early morning. Less than 1% of households were missed after two visits.

The wealth-index creation
For the wealth index, we selected 10 variables of household assets with the highest standard deviation (.0.20), as recommended by Seema Vyas and colleagues [21]. The types of asset variables and their standard deviations are presented in Table S1 in File S1. We transformed the categorical variables into dichotomous (0-1) indicators: 0 for indicators of poor wealth and 1 for indicators of good wealth. We examined the dichotomous variables by using the principal component analysis (PCA) to produce a factor score for each household with households being assigned a rank according to the factor score. Because of the low number of maternal deaths in the socioeconomic classes for calculation, we divided households into four equal categories (quartiles), rather than the widely used five classes. Each category was comprised of 25% of the households studied. Table S2 in File S1 shows the mean score, standard deviations, communalities, and correlations of the variables to the first (main) component. The total variance explained by the first component was 20.58%, with an eigenvalue of 2.06.

Data analysis
We used two units of analysis (household and birth). By using births as the unit of analysis, we presented descriptive tabulations of outcomes in the form of rates and ratios. By using the household as a unit of analysis, and applying logistic regression, we present household risk factors associated with the mortality outcomes (Tables 2, 3, and 4). We used SPSS 16 (Statistical Package for Social Sciences) for the data entry and analysis [22]. Data are freely available from the corresponding author on request. Table 1 describes households and population of the study area. Data collectors enumerated all households (11,920) in the selected 15 kebeles but collected data from 6,572 households that had pregnancy and birth outcomes in the last five years before the data collection. In the 6,572 households that had pregnancy and birth outcomes, there were 40,357 persons, an average of 6.1 persons per household.

Demographic description
Of the 6,572 household heads, 3,842 (58.5%) were not able to read and write (illiterate), 2,446 (37.2%) had an elementary education (grade 1-8), and 279 (4.2%) had completed 9th grade or more. Regarding the occupations of the head of the households, 6,289 (95.7%) engaged in farming, 204 (3.1%) in farming mixed with small trade, while 79 (1.2%) were salary employed ($ 25-50 per month) ( Table 2). Table 2 presents the backgrounds of households with pregnancy and birth outcomes and description of births. We found 11,762 births in 6572 households over 5 years (3.2% annual crude birth rate), of which 11,536 were live births. Of the 11,762 total births, 10,861 (92.3%) took place at home, 470 (4%) at health posts, 338 (2.9%) at health centres and 93 (0.8%) at hospitals. A health post is a two-room building which is the lowest ranked health facility and staffed by two community agents known as health extension workers whom received one year of general health training.

Description of births
Health centres are medium level institutions with health officers and nurses, and some have midwives. The overall institutional delivery rate (health centre and hospital combined) was 3.7% (431/11762).

Health centre births
Births that occurred within a 10 km distance of health centres, were five-times more likely to use health centres for delivery, 4.75% (274/5766), compared to households $11 km away, 1.07% (64/5996); OR: 4.62 (95% CI: 3.51-6.09) ( Table 2). A 10 km distance is the Ethiopian government plan to achieve health centre physical access to the population.

Household risk for mortality outcomes
Maternal mortality. Table 4 shows the MMR differences across different risk factors. The MMR was increased in households in the poorest quartile compared to the richest (550 vs 239 per 100,000 live births); OR: 2.29 (95% CI: 0.91-6.44). However, socio-economic factors examined (wealth, distance from road, education, and non-spaced births) did not have statistically significant association with maternal mortality because of the relative rarity of maternal deaths in terms of absolute numbers.
Neonatal mortality. Table 5 describes the household risk factors associated with neonatal mortality. Neonatal mortality was greater among the poorest quartile households compared to the richest; adjusted OR (AOR): 2.62 (95% CI: 1.65-4.15). However, the highest risk was in the wealthy class; AOR: 3.57 (95% CI:2.37-5.38). The poorest were in the second highest at risk groups. The likelihood of neonatal mortality was also increased among households with illiterate heads compared to where the heads had a higher education (9th grade or more); AOR: 3.54 (95% CI: 1.11-11.30), in households far from a motorable road ($6 km) compared to those within 5 km of a road; AOR: 2.40 (95% CI: 1. 56-3.69), and greater among households that had three or more births in five years compared to those that had two or less births; AOR: 3.22 (95% CI: 2. 45-4.22).
Stillbirths. Table 6 shows the factors associated with stillbirth which were greater in households in the poorest quartile compared to the richest; AOR: 3.13 (95% CI: 1.66-5.87). However, similar to neonatal mortality, the highest risk of stillbirth was in the households in the rich category compared to the richest; AOR: 6.40 (95% CI: 3.69-11.11). Stillbirth was greater among households far from motorable road ($6 km) compared to those at a distance of 5 km to the road; AOR: 3.40 (95% CI: 1.91-6.06), and greater in households that had three or more births in five years compared to two or fewer births; adjusted OR: 4.55 (95% CI: 3.35-6.16). Table 7 presents the concentration of maternal and newborn mortality in certain households. Of the 49 households that had maternal deaths, nearly half (46.9%) also experienced either a stillbirth or neonatal death; 12 (24.5%) had stillbirths, and 11 (22.4%) had neonatal deaths. The likelihood of neonatal death in households that had maternal deaths was seven times higher compared to households that had no maternal mortality; OR: 7.2 (95% CI: 3.6-14.3). The similar likelihood of having a stillbirth in households that had maternal mortality was 11 times greater compared to households without maternal death; OR: 11.6 (95% CI: 6.0-22.7).

Discussion
In this household study we found a maternal mortality ratio of 425 per 100,000 live births, a neonatal mortality ratio of 27 per 1,000 live births, and a stillbirth rate of 19 per 1,000 births. The risk of neonatal mortality was associated with the wealth status of the households, literacy status of the head, non-spaced births, and the distance from motorable road of households. The risk of stillbirth was also associated with the wealth, distance of the household from motorable road, and non-spaced births. The maternal mortality was also high among the poorest households compared to the richest and households that had maternal mortality also experienced a clustering of neonatal mortality or stillbirths. Moreover, the institutional delivery rate was unacceptably low.
To the best of our knowledge, this study is the first to describe three mortality estimates and associated household-risk factors using a large sample with high response rate in Gamo Gofa. In fact, there is a limited amount of evidence of maternal and neonatal mortality and stillbirths using community-based data from southern Ethiopia [14]. We used data collectors who had experience and were sensitive to the cultural taboos of the respective villages where they collected the data. Being the residents in their villages, they participated in all social events, including involvements in vital events such as celebrating births, caring for the sick, and funerals for the dead. The experience and deep knowledge of the area enabled them with the skill to handle sensitive questions and recall many of the deaths that occurred in their villages. We also used an experienced nurse as a field technical supervisor to help in classifying deaths by using the verbal autopsy method.
Because of lack of previous reports from the area, we compare our findings with community-based studies from other provinces and the national level estimates for Ethiopia. As a result, we do not know whether our findings were under-reported or represent the reality of the area. The maternal mortality ratio of 425 per 100,000 live births was similar to the findings of community-based studies decades ago in other parts of Ethiopia: the MMR per 100,000 live births was 402 in Jimma in 1990 [23], and 440 in Butajira in 1996 [24]. However, our estimate is higher than the UN and World Bank's estimate for Ethiopia of 350 per 100,000 live births in 2010 [1]. If we adjust our maternal mortality estimate upward by a factor of 1.6, it yields an MMR of 680 per 100,000 live births. The 2011 DHS reported MMR of 676 per 100,000 live births for Ethiopia [25] which is similar to our upwardly adjusted estimate. The suggestion as well as the factor of adjustment were conducted according to the recommendation by Santon and colleagues for the correction of the potential under-reporting of demographic studies [26]. The under-reporting of maternal deaths is a well-recognized global problem [27][28][29], and our study may not have escaped the challenge. However, given a general downward trend of Ethiopia's MMR estimate by the UN inter-  agency [1], and Hogan et al [30], our direct (unadjusted) estimate of an MMR of 425 per 100,000 may be close to reality. Our neonatal mortality estimate of 27 per 1,000 live births is lower than the DHS national estimate of 37 [31], and the UN Inter-agency estimate of 31 for Ethiopia in 2011 [20]. However, since mathematically modeled UN estimations often have inconsistencies with survey findings, and we have limited locally available estimates, we cannot judge whether our finding is due to under-reporting. Regarding stillbirths, Cousens et al estimated Ethiopia's stillbirth rate as 26 per 1,000 births in 2009 [9], which is still higher than our finding of 19 per 1,000 births. In the ethnic ''Gamo'' culture, a birth of a dead fetus and the early death of a newborn are not publicized and not publicly mourned. Only close relatives and family members are informed and other people are told the incident as ''something went wrong'' in local terms. This can make it difficult for an outsider to distinguish between a neonatal death and a stillbirth. Accordingly, though we tried to minimize under-reporting by carefully choosing the data collectors, some of the neonatal deaths and stillbirths may have been missed in this study. Global estimates suggest that 75% of neonatal deaths take place in the first week of life [4]. In our finding, firstweek deaths were approximately half of the neonatal deaths (46.4%) which may also suggest an under-reporting of early neonatal deaths and stillbirths.
The present study showed that indicators of the socio-economic status (SES) of households (wealth, education of head of household, the distance of households from a motorable road), and a factor related to reproductive health (non-spaced births) were associated with stillbirths and neonatal mortality outcomes. We examined the effect of household wealth (poor-rich differences) on mortality outcomes and found a significant variation. The poorest households were more likely to have mortality outcomes compared to the richest. Nevertheless, when compared to the richest quartile, the greatest risk of both neonatal mortality and stillbirth was to the rich households, not to the poorest.
Households in the poorest category were the second most at-risk group for the mortality outcomes. This finding needs further investigation with a focus on how the true wealth status can be determined at the particular stage of socio-economic development in the rural Ethiopia [21]. Our opinion is that the asset variables selected to indicate wealth status may not have correctly reflected the actual wealth status of all households. For example, farming land, cattle, and crops may not be owned by young people who returned to rural residence after certain years of education in urban areas. In addition, educated residents such as government employees who usually do not have these rural assets may have been wrongly classified as poor. Our finding of the association of education with neonatal mortality could also support the idea that the educated, but classified as poor and poorest in terms of rural wealth, may have been relatively better-off compared to those in the rich category because of the knowledge advantage.
On the other hand, households in the richest wealth category had the lowest risk of mortality possibly because of their economic access to health services and expected better living conditions. In general, the findings that the richest households had advantages while households in the three other categories experience greater risk could be due to the situation of wealth where few households have greater possession of assets and the other majority is homogeneously poor. The finding agrees with a previous analysis of Ethiopian rural asset data that classified the majority (up to 60%) of households as having low SES, thereby suggesting a homogeneity of most households in asset ownership [21].
Economists measure economic status indicators through information from income or expenditure, which is difficult to gather in low-income countries, and asset-based wealth is an alternative proxy in less developed areas [32]. Practically speaking, the wealth index is equally valid to income or expenditure data for health surveys in Africa [33], and the effect of household wealth on health outcome is well known [34][35][36]. Less clear, however, is how wealth causes mortality difference in areas where the overall access to health services and service utilization is very low to all people such as those in our study area. For instance, in our results, there is the poorest-richest difference in maternal mortality and stillbirths. Even so, very few households (including the richest) utilize health institutions for delivery service. This indicates that household wealth contributes to maternal mortality in mechanisms other than those that create economic access to health service. These mechanisms may include, e.g., improved, clean housing and better nutrition [37]. Unfortunately, this is beyond the aim and limitations of our study.
Regarding the association of wealth with neonatal mortality, the poorest-richest difference may be due to economic access to antibiotics and other drugs from rural private venders, where a wealthy family often has a better access in addition to better   nutrition and improved housing since access to antibiotics plays an important role in newborn survival [38]. In addition, the association of neonatal mortality with the education of the heads of a households and the distance to a motorable road further suggests the importance of these variables as tools to access health interventions. The education of parents has a positive correlation with better health of children through better knowledge of solutions and critical decisions during crisis, in addition to the opportunity education creates for job and economic access [39][40][41][42]. However, in the rural area where this study took place, few people had achieved better jobs despite higher levels of education. The knowledge advantage may have played a role related to access and utilization of treatments in households with educated heads.
Our study also demonstrated an association of household distance from a motorable road with neonatal death. Travel distance is clearly an important factor once the decision is made to seek medical care during critical conditions and distance to a motorable road has dual effects: 1) as a disincentive to seek health care, and 2) as a barrier to reach the relevant facility [43]. As such, people in households closer to a motorable-road are more likely to seek health care and save lives. We also found an increased likelihood of neonatal mortality and stillbirths among households where there was a maternal death. This is in agreement with a finding in a WHO multi-country maternal and newborn health survey, which showed a seven fold greater early neonatal mortality, in which mothers had died or developed nearmis (dangerous illness) compared to mothers without these events [44]. The clustering of maternal and neonatal mortality, as well as stillbirths in certain households, illustrates how impoverished households are trapped in many adverse outcomes.
The findings of a low utilization of facilities for skilled delivery, compounded with high mortality rates, call for educating mothers and other family members the importance of seeking skilled delivery service (institutional deliveries). Furthermore, improving the quality of the existing poorly equipped health institutions in the area [45] might help to increase the willingness and trust of families to utilize these institutions. The Ethiopian health extension workers have the opportunity to educate women and family, as they have close contact during the antenatal care and the routine household visits. In summary, our findings highlighted that households living in the villages far from road access, poor SES, and having illiterate household heads and non-spaced births, had greater risks of mortality. Therefore, a targeted follow-up of pregnancies in these households could help to achieve reduced mortality outcomes. Interventions such as family planning education and the availability of FP technology choices to support women in poor households may help prevent deaths caused due to the risk of non-spaced births.
In order to obtain ongoing data to help monitor progress, the community-based registration of pregnancy and birth outcomes (abortions, stillbirths, livebirths, neonatal and maternal deaths) may be the ideal option. In addition to providing sustainable data for evaluation of effectiveness of policy and programmes, ongoing registration can provide evidence for rights-based advocacy for improvement of health services, and many other benefits outside of the health sector [46]. For prospective registration to occur, Ethiopia must utilize the privilege of the available two health extension workers (HEWs) responsible for each kebele (average 500 households per HEW). Registration-based information may have benefits of reduced risk of recall bias and the cost of surveys by using the already available community health workers for active data collection.
In the following, we address some of the limitations of the study. First, recall bias and under-reporting are widely recognized problems in studies that ask respondents about past events. The intensity of the bias depends on the time interval between the event and the sensitivity of the event to memory [3]. We tried to reduce recall bias in two ways: 1) by selecting data collectors from the respective villages of the data collection that the data collectors helped the respondents to recall the events through their in depth knowledge of social events that happened in their villages, and 2) by choosing a memorable and short time reference period for the event to be recalled. However, there might have been some deaths that were missed due to recall bias in the current data. Second, as it is well-known in survey studies, we cannot show the temporality (time sequence) of the occurrence of exposures and outcomes [47]. For example, having more than three births in the last five years in households was associated with both neonatal mortality and stillbirth. Nevertheless, we cannot assure whether neonatal mortality and stillbirths led to more births or whether more births led to a greater risk.
Third, we used reported information from family members on mortality outcomes and mid-level expert decisions of the cause-ofdeath classifications using a simplified verbal autopsy technique, which also has the potential for misclassification (misdiagnosis). The existence of misclassification in using the verbal autopsy was reported from a well-designed prospective study of maternal deaths in Guinea-Bissau in Africa as 30% of maternal deaths were left unclassified [48]. The confirmatory diagnostic method used to ascertain the cause of maternal death is an autopsy test which has been used in a hospital setting in Mozambique in Africa [49]. However, such a modern technology cannot be applied in a rural community such as Bonke. Fourth, we were not able to show yearly changes in mortality. We aimed to describe at least aggregated measures by using a single reference period in a community where the date of the event is not easy to identify. From a practical point, it is difficult and inappropriate to expect a specific time-related response from a largely illiterate rural society where there is no vital registration system.

Conclusions
Mortality rates are still lagging behind the MDG targets for Ethiopia. There also exist socio-economic inequalities in maternal and neonatal mortalities, as well as stillbirths in the area. The socio-economic inequality in mortality and the low utilization of existing institutions for delivery care highlight the importance of quality emergency obstetric care service. The services need to target the poorest households where mortalities cluster and disproportionately high. It is important to address barriers to accessing institutional delivery services in a way that is acceptable for rural women.

Supporting Information
File S1 Supporting information file (two tables informing variables included in Principal Component Analysis (PCA). Table S1. Variables included in Principal Component Analysis (PCA) for wealth index creation.