Biometrics 62, 1170–1177
December 2006
DOI: 10.1111/j.1541-0420.2006.00609.x
Statistical Inference in a Stochastic Epidemic SEIR Model
with Control Intervention: Ebola as a Case Study
Phenyo E. Lekone
∗
and B¨arbel F. Finkenst¨adt
Department of Statistics, University of Warwick, Coventry CV4 7AL, U.K.
∗
email: lekonepe@mopipi.ub.bw
Summary. A stochastic discrete-time susceptible-exposed-infectious-recovered (SEIR) model for infectious
diseases is developed with the aim of estimating parameters from daily incidence and mortality time series
for an outbreak of Ebola in the Democratic Republic of Congo in 1995. The incidence time series exhibit
many low integers as well as zero counts requiring an intrinsically stochastic modeling approach. In order to
capture the stochastic nature of the transitions between the compartmental populations in such a model we
specify appropriate conditional binomial distributions. In addition, a relatively simple temporally varying
transmission rate function is introduced that allows for the effect of control interventions. We develop
Markov chain Monte Carlo methods for inference that are used to explore the posterior distribution of the
parameters. The algorithm is further extended to integrate numerically over state variables of the model,
which are unobserved. This provides a realistic stochastic model that can be used by epidemiologists to
study the dynamics of the disease and the effect of control interventions.
Key words: Control intervention; Ebola epidemics; Estimating transition rates; Latent process; Stochastic
SEIR model.
1. Introduction
Mathematical modeling has emerged as an important tool for
gaining understanding of the dynamics of the spread of in-
fectious diseases. The theoretical framework most commonly
used is based on the division of the human host population
into categories containing susceptible, infected but not yet
infectious (exposed), infectious, and recovered individuals.
These susceptible-exposed-infectious-recovered (SEIR) mod-
els are usually expressed as a system of differential equations
(see Anderson and May, 1991), where the rates of flow be-
tween compartments are determined by parameters specific
to the natural history of the disease. It is also recognized that
stochastic modeling is important (see, e.g., Bailey, 1975; Ball,
Mollison, and Scalia-Tomba, 1997; Andersson, 1999), in par-
ticular if disease incidence is small in the sense that the dis-
crete and stochastic nature of the transmission process may
not be neglected. This is the case for the observed Ebola out-
break studied in this article.
Inference for epidemic models is complicated by the fact
that first, one or several of the model variables may be unob-
served (latent), and second, data are often available at discrete
time points while the underlying true process is continuous in
time. In addition, model parameters may change over time,
for instance, if interventions to control the spread of the dis-
ease were introduced during the epidemics. In the case that
data are observed in continuous time, parameter estimates can
be obtained for complete data (Becker, 1976) and depend-
ing on the nature of incomplete data, estimators have also
been developed using martingale methods (Andersson and
Britton, 2000), exact forward-backward filtering (Fearnhead
and Meligkotsidou, 2004), or Markov chain Monte Carlo
(MCMC) methods (Gibson and Renshaw, 1998; O’Neill
and Roberts, 1999; O’Neill, 2002; Neal and Roberts, 2004;
Streftaris and Gibson, 2004). All of these studies are based
on the assumption that all event times or a subset thereof are
available to the investigator. Unfortunately, times at which
single events occur are rarely recorded. More commonly, ob-
served data sets are time series of counts of events that have
occurred during time intervals such as a day or a week.
In this study, we develop methods of inference in the case
that discrete-time epidemic data are available on an infec-
tious disease whose natural history pertains to an SEIR epi-
demic model. Likelihood-based inference for the case when
the time interval is equal to the generation period of the dis-
ease (chain-binomial model) has been considered by Bailey
(1975) and O’Neill and Roberts (1999). In the chain-binomial
model it is assumed that the generation period, that is, the
latent and infectious periods taken together, is fixed. Here we
relax this assumption by allowing both the latent and infec-
tious period to be stochastic and the data to be observed at
time points whose distances may be different from the length
of the generation period. The introduction of probability den-
sities for the transition of state variables allows us to formu-
late a probabilistic discrete-time model that, for a sufficiently
small interval length, provides a good approximation to the
underlying continuous-time process generated by the stochas-
tic SEIR model. The likelihood function of the data can then
be approximated on the basis of the transition densities.
Depending on what measurements are available, re-
searchers dealing with infectious disease data are confronted
1170
C
2006, The International Biometric Society
Statistical Inference in a Stochastic Epidemic SEIR Model
1171
with different scenarios of what are important observed and
unobserved variables. The first step is to construct a proba-
bility model for the disease to be studied and to investigate
parameter identifiability under the scenario of the available
data. The use of MCMC allows for numerical integration over
the distribution of unobserved variables. This is important
in achieving parameter estimates with standard errors that
realistically reflect increased uncertainty when variables are
unobserved. We illustrate the performance of our approach
by applying it to a data set of an outbreak of Ebola in the
Democratic Republic of Congo in 1995. Chowell et al. (2004)
considered only a part of the available data on this outbreak
for estimating the parameters of a set of deterministic SEIR
differential equations with a time-varying transmission rate to
allow for the control intervention. Their estimation approach
is based on simulating the solutions to the deterministic SEIR
equations and identifying the parameter values that minimize
the sum of squared errors between the observed and simu-
lated cumulative number of cases. The optimization process
was started from 10 different initial parameter values and the
reported parameter estimate is the one that resulted in the
smallest sum of squares of error.
In this article we introduce a fully probabilistic model and
therefore perform a rigid likelihood-based inference using all
available data while incorporating uncertainty about unob-
served variables as well as errors in the reporting process.
More generally, we provide a model that can be used by
epidemic researchers to realistically simulate the stochastic
dynamics of Ebola epidemics in order to study the effect of
control intervention and other questions of biological interest.
We describe the data and the modeling approach in the
next section and, in Section 3, address inference based on
the likelihood of the model. We start by considering the case
0 20 40 60 80 100 120 140
5
10
15
Time of symptom onset (days)
Cases
0 20 40 60 80 100 120 140
5
10
15
Time of death (days)
Cases
Figure 1. Data from the 1995 Ebola outbreak in the Democratic Republic of Congo recorded from March 1 (corresponding
to day 1 on the x-axis) to July 16. Top: daily counts of Ebola cases by date of symptom onset. Bottom: daily counts of death
from Ebola.
that all relevant variables are observed at daily time intervals.
In reality, the number of susceptible individuals who become
infected is a latent process and this, as well as some under-
reporting of variables, has to be considered by the estima-
tion algorithm. The application to the Ebola data produces
parameter estimates that characterize the natural history of
this disease and allows us to investigate the effectiveness of
the control intervention that took place.
2. Data and Model
2.1 Data
Ebola hemorrhagic fever, commonly known as Ebola, is trans-
mitted via physical contact with body fluids, secretions, tis-
sues, or semen from infected individuals (Chowell et al.,
2004). The disease is characterized by initial flu-like symp-
toms, which rapidly progress to vomiting, diarrhea, rash, and
internal and external bleeding. Once exposed, individuals go
through a latent period of approximately 6.3 days after which
they become infectious for a period that is estimated to be be-
tween 3.5 and 10.7 days (Breman et al., 1977). Most cases die
within 10 days of their initial infection, with the disease hav-
ing a mortality rate of 50–90% (World Health Organization,
2003). There are two known strains of the Ebola virus, the
Ebola-Sudan and the Ebola-Zaire, named after the countries
where they were first discovered (World Health Organization,
2003). Here we study data from the 1995 Ebola outbreak in
the Democratic Republic of Congo, which is well documented
by Khan et al. (1999). The data consist of two time series
(see Figure 1) recorded from March 1 to July 16, namely,
daily counts of Ebola cases by date of symptom onset, ac-
counting for a total of 291 cases, and daily counts of deaths
from Ebola, accounting for a total of 236 deaths. It is also
documented that the first case became ill on January 6, 1995,
1172 Biometrics, December 2006
the last case died on July 16, and a total of 316 cases were
identified resulting in a rate of 81% fatality. The Ebola virus
was confirmed as the causative agent on May 9 when tests
were carried out on specimens collected from some of the early
cases. Control measures were immediately introduced. These
included, among others, the use of protective clothing, active
surveillance, and community education (Khan et al., 1999).
The epidemic lasted for about 200 days with control mea-
sures being introduced about 130 days after the start of the
epidemic. The exact starting time and evolution of the epi-
demic prior to March 1 is unobserved. Furthermore, from the
total number of 316 identified cases, it can be deduced that
the dates of symptom onset for 25 cases and the dates of re-
moval from the infectious class for 80 cases are not reported
in the given time series.
2.2 Model
Consider a time interval (t, t + h], where h represents the
length between the time points at which measurements are
taken, here h =1day.LetB(t) denote the number of sus-
ceptible individuals who become infected, C(t) the number
of cases by date of symptom onset, and D(t) the number of
cases who are removed (die or recover) from the infectious
class during that time interval. Furthermore, let τ
∗
denote
the time point when the epidemic goes extinct, that is, the
first time point at which there are no exposed or infectious
individuals in the population. Let
B
= {B(t)}
τ∗
t=0
represent
the time series of B(t) from the beginning to the end of the
epidemic and define C and D similarly. We use a discrete-time
approximation to the stochastic continuous-time SEIR model
(see Gibson and Renshaw, 1998). Define S(t), E(t), I(t), and
R(t) as the number of susceptible, exposed, infectious, and
removed individuals in the population at time t, respectively.
Given initial conditions S(0) = s
0
, E(0) = e
0
, I(0) = a, and
the population size N, the discretized stochastic SEIR model
is specified by
S(t + h)=S(t) − B(t), (1)
E(t + h)=E(t)+B(t) − C(t), (2)
I(t + h)=I(t)+C(t) − D(t), (3)
S(t)+E(t)+I(t)+R(t)=N, (4)
where
B(t) ∼ Bin(S(t),P(t))
,C(t) ∼ Bin(E(t),p
C
),
D (t) ∼ Bin(I(t),p
R
)
(5)
are random variables with binomial Bin(n, p) distributions
with probabilities:
P (t)= 1− exp
−
β(t)
N
hI (t)
,p
C
=1− exp(−h),
p
R
=1−exp(−γh).
(6)
The parameters β(t), 1/, and 1/γ are the time-dependent
transmission rate, the mean incubation period, and the mean
infectious period, respectively. Mode and Sleeman (2000) de-
rived binomial densities as specified in (5, 6) for an SIR model
with three disease stages or compartments. The basic idea is
that transitions of individuals from the previous to the next
stage of the disease are seen as stochastic movements be-
tween the corresponding population compartments. In each
period an individual either stays or moves on to the next
compartment. Assuming that the time length that an indi-
vidual spends in a compartment is exponentially distributed
with some compartment-specific rate λ(t), then the probabil-
ity of extending the stay by a further period of length h is
exp(−λ(t)h) and the probability of leaving is therefore
1 − exp(−λ(t)h). The binomial distributions (5) result from
summation over the individual Bernoulli trials assuming that
they are independent and identical for all members of a com-
partment. Furthermore, noting that the compartment-specific
exponential rates are
β(t)
N
I(t),, and γ for the susceptible, ex-
posed, and infectious compartment leads to the probabilities
of staying in a compartment as specified in (6) (see Mode and
Sleeman, 2000). It follows that the exponential distribution
of the incubation and the infectious period is approximated
by the corresponding geometric distribution with means 1/p
C
and 1/p
R
, respectively. Conditional on all information up to
time t, the binomial random variables B(t), C(t), and D(t)
are independent. The model further assumes that the pop-
ulation size N remains constant and that individuals mix
homogeneously.
In order to account for the control intervention we assume
that the transmission parameter β(t) is constant up to the
time point when the control measures are introduced and after
that decays exponentially. This can be formulated as
β(t)=
β, t<t
∗
βe
−q(t−t
∗
)
,t≥ t
∗
,
(7)
where t
∗
is the time point at which control measures are in-
troduced, β is the initial transmission rate, and q > 0isthe
rate at which β(t) decays for t > t
∗
. Chowell et al. (2004)
consider an exponential decay of the transmission rate but
introduce a third parameter that additionally characterizes
the decay. We found that for realistic sample sizes the geom-
etry of the likelihood function does not permit identification
of this parameter as its estimator is correlated with the expo-
nent q. Our model for the transmission rate in (7) is therefore
a more parsimonious parameterization. Note that the inter-
vention does not affect γ unless the disease is curable, which
is not the case for Ebola. The basic reproduction number R
0
is defined as the average number of secondary cases gener-
ated by a primary case over his/her infectious period when
introduced into a large population of susceptible individuals
(Diekmann, Heesterbeek, and Metz, 1990). The constant R
0
thus measures the initial growth rate of the epidemic and for
the model above it can be shown that R
0
= β/γ (Hethcote,
2000). Furthermore, Chowell et al. (2004) define the time-
dependent effective reproductive number R
0
(t)=
β(t)
γ
S(t)
N
as
the number of secondary cases per infectious case at time
t. Because S(t) ≈ N, it follows that R
0
(t) ≈
β(t)
γ
is a function
proportional to the time-varying transmission rate in (7). The
time point at which R
0
(t) assumes values smaller than 1 in-
dicates when control measures have become effective in con-
trolling the epidemic.
The epidemic model specified in (1)–(6) together with the
contact rate model (7) has parameter vector Θ = {β, q,
, γ}, which we would like to estimate from knowledge of
initial conditions, population size, and from observation of
Statistical Inference in a Stochastic Epidemic SEIR Model
1173
{B, C, D} or a subset thereof. The temporal evolution
of the effective R
0
(t) is then derived from the estimated
parameters.
3. Inference
If initial conditions, the population size, and vectors
{B, C, D} are observed at time intervals of length h, then
here we say that the data are complete. Note that the time
series for {S(t), E(t), I(t)} are in this case fully determined
by applying equations (1)–(3) from given initial conditions.
Because B(t), C(t), and D(t) are conditionally independent,
the likelihood of the data can be approximated by
L(
B
,
C
,
D
|Θ) =
τ
∗
t=0
g
1
(B(t) |·)g
2
(C(t) |·)g
3
(D(t) |·), (8)
where g
1
, g
2
, and g
3
stand for the binomial transition densities
specified in (5) and (6) conditioned on Θ and on all the infor-
mation up to time t. The maximum likelihood (ML) estimator
for Θ, and subsequently for R
0
and R
0
(t), can be obtained by
maximizing (8).
As is the case in reality we furthermore assume that B
is not observed. For the purpose of this article we perform
data imputation within the framework of Bayesian MCMC
methods, which allows us to numerically integrate over the
probability distribution of the unobserved process. Because
the epidemic is observed until the end, we have that E(τ
∗
)=
0, I(τ
∗
) = 0, and the final number of susceptible individuals
is given by S(τ
∗
)=N −
τ
∗
j=0
C(t). The final size of the epi-
demic, defined as the total number of individuals who eventu-
ally contract the disease, is in this case given by m = S(τ
∗
) −
s
0
. The series {I(t)} is also known as it can be reconstructed
from observed variables via equation (3). The series {S(t),
E(t)}, however, depends on the unobserved B. At each sweep
of the Markov chain we now impute the stochastic process
B and reconstruct the values of {S(t)} and {E(t)} using (1)
and (2). The augmented likelihood of B, C, and D is L(B, C,
D |Θ) as given in (8). Multiplying the likelihood by the prior
π(Θ) gives, up to a constant of proportionality, the posterior
distribution
π(Θ,
B
|
C
,
D
) ∝ L(
B
,
C
,
D
|Θ)π(Θ) (9)
that we wish to sample from. An MCMC algorithm that sam-
ples in turn from the conditional distributions π(B |C, D,Θ)
and π(Θ |C, D, B) produces draws from the desired π(Θ,
B |C, D). The general structure of the algorithm is thus as
follows:
(1) Initialize B. This can be done, for example, by setting
B(0) = m, with all the other positions in the vector filled
with zeros. For any initial B the series {S(t), E(t)} are
reconstructed using (1) and (2), respectively.
(2) Initialize the parameter vector Θ.
(3) Update B from B |C, D, Θ and compute new series for
{S(t)} and {E(t)} from (1) and (2).
(4) Update Θ from Θ |C, D, B.
(5) Repeat steps (3) and (4) until the required sample is
obtained after the chain has converged.
3.1 Sampling from B |C, D,Θ
A natural way of proposing B would be to sample each
B(
t) from its conditional binomial distribution at each
time point. In addition, it is checked that proposals are
consistent with the observed final size and length of the
observed epidemic
τ
∗
−1
t=0
B(t)=m, E(t) ≥ 0,E(t)+I(t) >
0, ∀t<τ
∗
, and I(τ
∗
)+E(τ
∗
) = 0. Unfortunately, this is rarely
the case for such proposals and thus this scheme is not very
efficient. Instead, we explicitly condition our proposals on
the observed extinction time by using the following updat-
ing scheme. To update the current realization of B, select a
time point t
(satisfying B(t
) > 0) uniformly at random from
{0, h,2h, ..., τ
∗
− h} and set B(t
)=B(t
) − 1. Then, select
˜
t
∈{0, h, ..., τ
∗
− h} uniformly at random and set B(
˜
t)=B(
˜
t)
+ 1. Update the series {S(t), E(t)} using (1) and (2), respec-
tively, and check that all the conditions above are satisfied. If
so then the new configuration B
is taken as a candidate pro-
posal for B and we accept B
with probability min[1,
π(
B
|·)
π(
B
|·)
].
At each iteration of the Markov chain, more than one ele-
ment of B may be updated by iterating the procedure above
a number of times (about 10% of the final size). This is known
to improve mixing and convergence of the MCMC algorithm
(Neal and Roberts, 2004).
3.2 Sampling from Θ |C, D, B
We update each element of Θ using a random walk proposal
where the variance of the Gaussian perturbations is tuned
such that the overall acceptance rate is between 20% and 40%
(Roberts and Rosenthal, 2001). Preliminary analysis showed
that imputation of B introduces a negative correlation be-
tween the chains for q and and we find that convergence is
improved if we sample them jointly and independently of their
previous values using a bivariate normal proposal centered at
the mode of the bivariate conditional distribution (see, e.g.,
Chib and Greenberg, 1994). Independent gamma priors are
assigned to each of the parameters in Θ, that is, π(ζ) ∼
Γ(ν
ζ
, λ
ζ
), where ζ = β, q, γ, and , where Γ(a, b) refers
to a gamma distribution with parameters a and b, mean a/b,
and variance a/b
2
.
3.3 Application to Simulated Data
First we test and demonstrate the performance of the algo-
rithm by applying it to simulated epidemic time series from
model (1)–(6) with parameter values s
0
= 5,364,500, e
0
=1,
a =0,β = 0.2, q = 0.2, =1/5 = 0.2, γ =1/7 ≈ 0.143,
h = 1 day, and intervention time t
∗
= 130. These parameter
values are, as far as possible, tuned to the data as they are
motivated by earlier studies on Ebola and by what is known
about its natural history (Breman et al., 1977; Khan et al.,
1999; Chowell et al., 2004). The simulated epidemic time se-
ries resulted in a final size of m = 297 and terminated at
τ
∗
= 172.
We apply the MCMC algorithm described in the pre-
vious section to the simulated data, pretending that B is
unobserved, and conduct inference using three parameter sets
for prior distributions. First {(ν
ζ
, λ
ζ
); ζ = β, , γ, q} = {(2,
10), (2, 10), (2, 14), (2, 10)}, second {(20, 100), (20, 100),
(20, 140), (2, 10)}, and finally {(20, 40), (20, 40), (20, 40),
(20, 40)}. In the first two cases the mean of the distribution
is the true parameter value, with the priors being more infor-
mative in the latter. In the third case, the mean of each dis-
tribution is 0.5, which is substantially different from the true
values. We shall refer to this case as the noncentered prior
case. The prior distribution for q is chosen to be less informa-
tive in the first two cases because it is difficult to know a priori
1174 Biometrics, December 2006
Table 1
Parameter estimates for simulated data. The first row displays the true parameter setting for the simulation.
ML estimates for complete data are shown in the second row. The third, fourth, and last rows give posterior
means using a vague, informative, and noncentered prior distributions, respectively. Posterior standard
deviations are given in parentheses.
Method
ˆ
β ˆq ˆ ˆγ
ˆ
R
0
True values 0.2 0.2 0.2 0.143 1.4
MLE 0.194 (0.0083) 0.170 (0.0173) 0.201 (0.0117) 0.144 (0.00837) 1.35 (0.097)
MCMC (vague) 0.191 (0.0115) 0.166 (0.0247) 0.192 (0.0246) 0.144 (0.00863) 1.33 (0.115)
MCMC (informative) 0.192 (0.0112) 0.164 (0.0246) 0.195 (0.0222) 0.144 (0.00783) 1.33 (0.104)
MCMC (noncentered) 0.203 (0.0120) 0.193 (0.0249) 0.211 (0.0237) 0.151 (0.00870) 1.35 (0.110)
the effect of the intervention. In each case, after discarding a
burn-in period of 10,000, a sample of size 1000 taken every 100
iterations of the chain was used to obtain posterior distribu-
tions. At each iteration of the Markov chain B was updated 30
times. Convergence of the Markov chain was assessed using a
series of runs for different starting values and also inspecting
the autocorrelation function. In all three cases the Markov
chain appeared to have converged after the burn-in period.
The posterior means and standard deviations of the parame-
ters are reported in Table 1. The posterior means are well in
agreement with the ML estimates obtained from the complete
data. As can be expected, the posterior standard deviations
for β, , and q are larger demonstrating that the algorithm has
incorporated the increased level of uncertainty, as effectively
about 30% of the complete data are not available. Uncertainty
about the parameter γ is not affected by latent data, because
the component of the likelihood involving this parameter de-
pends only on {I(t)}, which is available from the observed
data. Table 1 suggests that the estimation method is not very
sensitive to the choice of prior information, as different choices
result in parameter estimates that are in agreement with the
ML estimates. The same estimation algorithm was also ap-
0.15 0.20 0.25 0.30
β
Informative vague
4 5 6 7 8 9
1/γ
8 9 10 11 12 13
1/ρ
0.13 0.14 0.15 0.16 0.17 0.18 0.19
q
Figure 2. Estimated posterior distributions for observed Ebola data with a vague prior (solid curves) and informative prior
(dotted curves).
plied imputing about 10% unknown dates on which exposed
individuals became infectious and about 25% unknown dates
on which infected cases recovered. Results not reported here
show that imputing these unknown dates together with B re-
sulted in posterior means in agreement with results presented
in Table 1.
4. Inference from Observed Ebola Data
4.1 Estimation and Results
The algorithm described above is applied to the real Ebola
data set where we impute the unobserved process B as well
as the unknown dates on which 25 cases became infectious
and on which 80 cases were removed. We assume that there
was initially one exposed individual and that the very first
case was exposed for 6 days (approximate mean incubation
period) prior to being diagnosed. Parameter constraints are
given by 0 <β<1; q > 0 and, following Chowell et al.
(2004), we choose initial conditions from 1 < 1/<21 and
3.5 < 1/γ<10.7 as well as their estimated population size
of N = 5,363,500 as in the simulations above. The choice of
prior distributions and the sampling methods are exactly the
same as in the previous section. Figure 2 shows the estimated
