pompExamples: Additional POMP Examples

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")