library(pomp)
library(plyr)
library(reshape2)
library(magrittr)
options(stringsAsFactors=FALSE)
##' ## Components of the Lavine et al (2013) model.
##' the process model simulator
Csnippet("
int nrate = 6;
int nseas = 3;
double trans[nrate]; // transition numbers
double Beta, lambda, nu, pboost, leaveW, boost, wane2;
double dQ;
int k;
double surv = exp(-mu*dt);
const double *seas = &seas1;
const double *logbeta = &logbeta1;
for (k = 0, Beta = 0; k < nseas; k++) Beta += logbeta[k]*seas[k];
Beta = exp(Beta);
dQ = rgammawn(beta_sd,dt);
lambda = Beta*I/pop*dQ/dt;
nu = kappa*lambda;
pboost = (kappa == 0) ? 0.0 : ((sigma == 0) ? 1.0 : nu/(nu+sigma));
leaveW = nearbyint(W*(1-exp(-(nu+sigma)*dt)));
boost = nearbyint(pboost*leaveW);
wane2 = leaveW-boost;
// compute the transition numbers
trans[0] = rpois(birth*dt); // births are Poisson
trans[1] = rbinom(S,1-exp(-lambda*dt)); // infections
trans[2] = nearbyint(I*(1-exp(-gamma*dt))); // recoveries
trans[3] = nearbyint(R*(1-exp(-sigma*dt))); // R->W
trans[4] = boost; // W->R
trans[5] = wane2; // W->S
// balance the equations
S += trans[0]-trans[1]+trans[5];
I += trans[1]-trans[2];
R += trans[2]+trans[4]-trans[3];
W += trans[3]-trans[4]-trans[5];
CASE += trans[2]; // cases are cumulative recoveries
BOOSTS += boost; // W->R events
WANES += wane2; // W->S events
if (beta_sd > 0.0) {
Q += (dQ-dt)/beta_sd; // mean zero, variance = dt
}
S = nearbyint(surv*S);
I = nearbyint(surv*I);
R = nearbyint(surv*R);
W = nearbyint(surv*W);
") -> rproc
##' the deterministic skeleton
Csnippet("
int nrate = 6;
int nseas = 3;
double trans[nrate]; // transition numbers
double Beta, lambda, nu, pboost, leaveW, boost, wane2;
int k;
double surv = exp(-mu);
const double *seas = &seas1;
const double *logbeta = &logbeta1;
for (k = 0, Beta = 0; k < nseas; k++) Beta += logbeta[k]*seas[k];
Beta = exp(Beta);
lambda = Beta*I/pop;
nu = kappa*lambda;
pboost = nu/(nu+sigma);
leaveW = W*(1-exp(-(nu+sigma)));
boost = pboost*leaveW;
wane2 = leaveW-boost;
// compute the transition numbers
trans[0] = birth; // births
trans[1] = S*(1-exp(-lambda)); // infections
trans[2] = (I*(1-exp(-gamma))); // recoveries
trans[3] = (R*(1-exp(-sigma))); // R->W
trans[4] = boost; // W->R
trans[5] = wane2; // W->S
// balance the equations
DS = (S+trans[0]-trans[1]+trans[5])*surv;
DI = (I+trans[1]-trans[2])*surv;
DR = (R+trans[2]+trans[4]-trans[3])*surv;
DW = (W+trans[3]-trans[4]-trans[5])*surv;
DCASE = trans[2]; // cases are cumulative recoveries
DBOOSTS = boost; // W->R events
DWANES = wane2; // W->S events
") -> skel
##' the measurement model density
Csnippet("
double tol = 1.0e-20;
double mean, sd, f;
mean = obsprob*CASE;
sd = sqrt(mean*(1-obsprob)+tau*CASE*CASE)+tol;
if (reports > 0)
f = pnorm(reports+0.5,mean,sd,1,0)-pnorm(reports-0.5,mean,sd,1,0);
else
f = pnorm(0.5,mean,sd,1,0);
lik = (give_log) ? log(f+tol) : f+tol;
") -> dmeas
##' the measurement model simulator
Csnippet("
double tol = 1.0e-20;
double mean, sd, rep;
mean = obsprob*CASE;
sd = sqrt(mean*(1-obsprob)+tau*CASE*CASE)+tol;
rep = nearbyint(rnorm(mean,sd));
reports = (rep > 0) ? rep : 0;
") -> rmeas
##' transformation to the estimation scale
Csnippet("
double sum;
Tgamma = exp(gamma);
Tsigma = exp(sigma);
Tkappa = exp(kappa);
Ttau = exp(tau);
Tmu = exp(mu);
Tbeta_sd = exp(beta_sd);
Tobsprob = expit(obsprob);
TS_0 = exp(S_0);
TI_0 = exp(I_0);
TW_0 = exp(W_0);
TR_0 = exp(R_0);
sum = TS_0+TI_0+TW_0+TR_0;
TS_0 /= sum;
TI_0 /= sum;
TW_0 /= sum;
TR_0 /= sum;
") -> totrans
##' transformation from the estimation scale
Csnippet("
double sum;
Tgamma = log(gamma);
Tsigma = log(sigma);
Tkappa = log(kappa);
Ttau = log(tau);
Tmu = log(mu);
Tbeta_sd = log(beta_sd);
Tobsprob = logit(obsprob);
sum = S_0+I_0+W_0+R_0;
TS_0 = log(S_0/sum);
TI_0 = log(I_0/sum);
TW_0 = log(W_0/sum);
TR_0 = log(R_0/sum);
") -> fromtrans
##' initializer
Csnippet("
double tot;
tot = S_0+I_0+W_0+R_0;
S = nearbyint(pop*S_0/tot);
I = nearbyint(pop*I_0/tot);
W = nearbyint(pop*W_0/tot);
R = nearbyint(pop*R_0/tot);
CASE = 0;
BOOSTS = 0;
WANES = 0;
Q = 0;
") -> init
statenames <- c("S","I","W","R","CASE","BOOSTS","WANES","Q")
zeronames <- c("WANES","BOOSTS","CASE","Q")
c(logbeta1=log(5.02),logbeta2=log(3.22),logbeta3=log(4.80),
gamma=0.273,sigma=0.00113,kappa=6.56,tau=0.000755,
mu=1.28e-04,beta_sd=0.512,obsprob=0.152,
S_0=0.0493,I_0=0.000770,R_0=0.947,W_0=0.00294
) -> mle
##' ## The data, which consist of:
##' - pertussis case reports
##' - smoothed weekly birth numbers
##' - smoothed population sizes
##' all for the city of Copenhagen.
Lavine2013 %>%
mutate(week=julian(date,origin=as.Date("1899-12-31"))/7) -> dat
dat %<>% cbind(
periodic.bspline.basis(x=dat$week,nbasis=3,degree=3,
period=52,names="seas%d"))
dat %>%
subset(date<"1938-01-01",select=c("week","reports")) %>%
tail(-1) %>%
pomp(
times="week", t0=0,
covar=subset(dat,select=-c(date,reports)), tcovar="week",
paramnames=c("logbeta1","gamma","sigma","kappa","tau",
"mu","beta_sd","obsprob","S_0","I_0","R_0","W_0"),
statenames=statenames, zeronames=zeronames,
rprocess=discrete.time.sim(step.fun=rproc,delta.t=1),
skeleton=map(skel,delta.t=1),
rmeasure=rmeas, dmeasure=dmeas,
fromEstimationScale=fromtrans, toEstimationScale=totrans,
initializer=init, params=mle
) -> lavine2013
c("lavine2013")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.