This document is derived from King, Domenech de Cellès, Magpantay, and Rohani (2015) and the data supplement archived on datadryad.org.
Produced in R version 4.3.2 using pomp version 5.6.
All codes needed to reproduce the results of the original King et al. (2015) paper are available on datadryad.org.
Download the data from the WHO Situation Report of 1 October 2014:
read_csv("https://kingaa.github.io/sbied/ebola/ebola_data.csv") -> dat
dat
# A tibble: 74 × 4
week date country cases
<dbl> <date> <chr> <dbl>
1 1 2014-01-05 Guinea 2.24
2 2 2014-01-12 Guinea 2.24
3 3 2014-01-19 Guinea 0.073
4 4 2014-01-26 Guinea 5.72
5 5 2014-02-02 Guinea 3.95
6 6 2014-02-09 Guinea 5.44
7 7 2014-02-16 Guinea 3.27
8 8 2014-02-23 Guinea 5.76
9 9 2014-03-02 Guinea 7.62
10 10 2014-03-09 Guinea 7.62
# ℹ 64 more rows
Supplementing these data are population estimates for the three countries. These are census figures from 2014.
c(Guinea=10628972,Liberia=4092310,SierraLeone=6190280) populations <-
|>
dat ggplot(aes(x=date,y=cases,group=country,color=country))+
geom_line()
Many of the early modeling efforts used variants on the simple SEIR model. Here, we’ll focus on a variant that attempts a more careful description of the duration of the latent period. Specifically, this model assumes that the amount of time an infection remains latent is LP∼Gamma(m,1mα), where m is an integer. This means that the latent period has expectation 1/α and variance 1/(mα). In this document, we’ll fix m=3.
We implement Gamma distributions using the so-called linear chain trick.
Csnippet("
rSim <- double lambda, beta;
double *E = &E1;
beta = R0 * gamma; // Transmission rate
lambda = beta * I / N; // Force of infection
int i;
// Transitions
// From class S
double transS = rbinom(S, 1.0 - exp(- lambda * dt)); // No of infections
// From class E
double transE[nstageE]; // No of transitions between classes E
for(i = 0; i < nstageE; i++){
transE[i] = rbinom(E[i], 1.0 - exp(-nstageE * alpha * dt));
}
// From class I
double transI = rbinom(I, 1.0 - exp(-gamma * dt)); // No of transitions I->R
// Balance the equations
S -= transS;
E[0] += transS - transE[0];
for(i=1; i < nstageE; i++) {
E[i] += transE[i-1] - transE[i];
}
I += transE[nstageE-1] - transI;
R += transI;
N_EI += transE[nstageE-1]; // No of transitions from E to I
N_IR += transI; // No of transitions from I to R
")
Csnippet("
rInit <- double m = N/(S_0+E_0+I_0+R_0);
double *E = &E1;
int j;
S = nearbyint(m*S_0);
for (j = 0; j < nstageE; j++) E[j] = nearbyint(m*E_0/nstageE);
I = nearbyint(m*I_0);
R = nearbyint(m*R_0);
N_EI = 0;
N_IR = 0;
")
The deterministic skeleton is a vectorfield (i.e., a system of ordinary differential equations). The following C snippet computes the components of this vectorfield as functions of the state variables and parameters.
Csnippet("
skel <- double lambda, beta;
const double *E = &E1;
double *DE = &DE1;
beta = R0 * gamma; // Transmission rate
lambda = beta * I / N; // Force of infection
int i;
// Balance the equations
DS = - lambda * S;
DE[0] = lambda * S - nstageE * alpha * E[0];
for (i=1; i < nstageE; i++)
DE[i] = nstageE * alpha * (E[i-1]-E[i]);
DI = nstageE * alpha * E[nstageE-1] - gamma * I;
DR = gamma * I;
DN_EI = nstageE * alpha * E[nstageE-1];
DN_IR = gamma * I;
")
Ct|Ht is negative binomial with E[Ct|Ht]=ρHt and Var[Ct|Ht]=ρHt(1+kρHt).
Csnippet("
dObs <- double f;
if (k > 0.0)
f = dnbinom_mu(nearbyint(cases),1.0/k,rho*N_EI,1);
else
f = dpois(nearbyint(cases),rho*N_EI,1);
lik = (give_log) ? f : exp(f);
")
Csnippet("
rObs <- if (k > 0) {
cases = rnbinom_mu(1.0/k,rho*N_EI);
} else {
cases = rpois(rho*N_EI);
}")
The following function constructs a pomp
object to hold the data for any one of the countries. It demonstrates one level of abstraction above the basic pomp
constructor.
function (country=c("Guinea", "SierraLeone", "Liberia"),
ebolaModel <-timestep = 0.1, nstageE = 3) {
match.arg(country)
ctry <- unname(populations[ctry])
pop <- as.integer(nstageE)
nstageE <-
paste0("static int nstageE = ",nstageE,";")
globs <-
|> filter(country==ctry) |> select(-country) -> dat
dat
## Create the pomp object
|>
dat select(week,cases) |>
pomp(
times="week",
t0=min(dat$week)-1,
globals=globs,
accumvars=c("N_EI","N_IR"),
statenames=c("S",sprintf("E%1d",seq_len(nstageE)),
"I","R","N_EI","N_IR"),
paramnames=c("N","R0","alpha","gamma","rho","k",
"S_0","E_0","I_0","R_0"),
dmeasure=dObs, rmeasure=rObs,
rprocess=discrete_time(step.fun=rSim, delta.t=timestep),
skeleton=vectorfield(skel),
partrans=parameter_trans(
log=c("R0","k"),logit="rho",
barycentric=c("S_0","E_0","I_0","R_0")),
rinit=rInit
po
) ->
}
ebolaModel("Guinea") -> gin
ebolaModel("SierraLeone") -> sle
ebolaModel("Liberia") -> lbr
King AA, Domenech de Cellès M, Magpantay FMG, Rohani P (2015). “Avoidable Errors in the Modelling of Outbreaks of Emerging Pathogens, with Special Reference to Ebola.” Proc R Soc Lond B, 282(1806), 20150347. https://doi.org/10.1098/rspb.2015.0347.
Licensed under the Creative Commons Attribution-NonCommercial license. Please share and remix noncommercially, mentioning its origin.