Case study: Measles in large and small towns

Data sets

• Twenty towns, including
  • 10 largest
  • 10 smaller, chosen at random
• Population sizes: 2k–3.4M
• Weekly case reports, 1950–1963
• Annual birth records and population sizes, 1944–1963

Continuous-time Markov process model

Diagram of the model:

• $B(t) = \text{birth rate, from data}$
• $N(t) = \text{population size, from data}$

• Overdispersed binomial measurement model: $\mathrm{cases}_t\,\vert\,{{\Delta}{N}}_{IR}=z_t \sim {\mathrm{Normal}\left(\rho\,z_t,\rho\,(1-\rho)\,z_t+(\psi\,\rho\,z_t)^2\right)}$

• Entry into susceptible class: $$\mu_{BS}(t) = (1-c)\,B(t-\tau)+c\,\delta(t-t_0)\,\int_{t-1}^{t}\,B(t-\tau-s)\,ds$$

• Force of infection $$\mu_{SE}(t) = \tfrac{\beta(t)}{N(t)}\,(I+\iota)\,\zeta(t)$$

• $c = \text{cohort effect}$
  
• $\tau = \text{school-entry delay}$
  
• $\beta(t) = \text{school-term transmission} = \begin{cases}\beta_0\,(1+a) &\text{during term}\\\beta_0\,(1-a) &\text{during vacation}\end{cases}$
  
• $\iota = \text{imported infections}$

• $\zeta(t) = \text{Gamma white noise with intensity}\,\sigma_{SE}$ (He et al. 2010,Bhadra et al. (2011))

The (unobserved) process model

We propose a variant of the SEIR model as an explanation for these data. This is a compartmental model that, diagrammatically, looks as follows.

model diagram

$B = \text{births}$
$S = \text{susceptibles}$
$E = \text{exposed, incubating}$
$I = \text{infectious}$
$R = \text{recovered}$

We require a simulator for this model. The following code implements a simulator.

Notable complexities include:

1. Incorporation of the known birthrate.
2. The birth-cohort effect: a specified fraction (cohort) of the cohort enter the susceptible pool aall at once.
3. Seasonality in the transmission rate: during school terms, the transmission rate is higher than it is during holidays.
4. Extra-demographic stochasticity in the form of a Gamma white-noise term acting multiplicatively on the force of infection.
5. Demographic stochasticity implmented using Euler-multinomial distributions.

rproc <- Csnippet("
  double beta, br, seas, foi, dw, births;
  double rate[6], trans[6];
  
  // cohort effect
  if (fabs(t-floor(t)-251.0/365.0) < 0.5*dt) 
    br = cohort*birthrate/dt + (1-cohort)*birthrate;
  else 
    br = (1.0-cohort)*birthrate;

  // term-time seasonality
  t = (t-floor(t))*365.25;
  if ((t>=7&&t<=100) || (t>=115&&t<=199) || (t>=252&&t<=300) || (t>=308&&t<=356))
      seas = 1.0+amplitude*0.2411/0.7589;
    else
      seas = 1.0-amplitude;

  // transmission rate
  beta = R0*(gamma+mu)*seas;
  // expected force of infection
  foi = beta*pow(I+iota,alpha)/pop;
  // white noise (extrademographic stochasticity)
  dw = rgammawn(sigmaSE,dt);

  rate[0] = foi*dw/dt;  // stochastic force of infection
  rate[1] = mu;             // natural S death
  rate[2] = sigma;        // rate of ending of latent stage
  rate[3] = mu;             // natural E death
  rate[4] = gamma;        // recovery
  rate[5] = mu;             // natural I death

  // Poisson births
  births = rpois(br*dt);
  
  // transitions between classes
  reulermultinom(2,S,&rate[0],dt,&trans[0]);
  reulermultinom(2,E,&rate[2],dt,&trans[2]);
  reulermultinom(2,I,&rate[4],dt,&trans[4]);

  S += births   - trans[0] - trans[1];
  E += trans[0] - trans[2] - trans[3];
  I += trans[2] - trans[4] - trans[5];
  R = pop - S - E - I;
  W += (dw - dt)/sigmaSE;  // standardized i.i.d. white noise
  C += trans[4];           // true incidence
")

In the above, $C$ represents the true incidence, i.e., the number of new infections occurring over an interval. Since recognized measles infections are quarantined, we argue that most infection occurs before case recognition so that true incidence is a measure of the number of individuals progressing from the I to the R compartment in a given interval.

We complete the process model definition by specifying the distribution of initial unobserved states. The following codes assume that the fraction of the population in each of the four compartments is known.

initlz <- Csnippet("
  double m = pop/(S_0+E_0+I_0+R_0);
  S = nearbyint(m*S_0);
  E = nearbyint(m*E_0);
  I = nearbyint(m*I_0);
  R = nearbyint(m*R_0);
  W = 0;
  C = 0;
")

The measurement model

We’ll model both under-reporting and measurement error. We want $\mathbb{E}[\text{cases}|C] = \rho\,C$, where $C$ is the true incidence and $0<\rho<1$ is the reporting efficiency. We’ll also assume that $\mathrm{Var}[\text{cases}|C] = \rho\,(1-\rho)\,C + (\psi\,\rho\,C)^2$, where $\psi$ quantifies overdispersion. Note that when $\psi=0$, the variance-mean relation is that of the binomial distribution. To be specific, we’ll choose $\text{cases|C} \sim f(\cdot|\rho,\psi,C)$, where $$f(c|\rho,\psi,C) = \Phi(c+\tfrac{1}{2},\rho\,C,\rho\,(1-\rho)\,C+(\psi\,\rho\,C)^2)-\Phi(c-\tfrac{1}{2},\rho\,C,\rho\,(1-\rho)\,C+(\psi\,\rho\,C)^2),$$ where $\Phi(x,\mu,\sigma^2)$ is the c.d.f. of the normal distribution with mean $\mu$ and variance $\sigma^2$.

The following computes $\mathbb{P}[\text{cases}|C]$.

dmeas <- Csnippet("
  double m = rho*C;
  double v = m*(1.0-rho+psi*psi*m);
  double tol = 1.0e-18;
  if (cases > 0.0) {
    lik = pnorm(cases+0.5,m,sqrt(v)+tol,1,0)-pnorm(cases-0.5,m,sqrt(v)+tol,1,0)+tol;
  } else {
    lik = pnorm(cases+0.5,m,sqrt(v)+tol,1,0)+tol;
  }
")

The following codes simulate $\text{cases} | C$.

rmeas <- Csnippet("
  double m = rho*C;
  double v = m*(1.0-rho+psi*psi*m);
  double tol = 1.0e-18;
  cases = rnorm(m,sqrt(v)+tol);
  if (cases > 0.0) {
    cases = nearbyint(cases);
  } else {
    cases = 0.0;
  }
")

Constructing the pomp object

We put all the model components together with the data in a call to pomp:

dat %>% 
  pomp(t0=with(dat,2*time[1]-time[2]),
       time="time",
       rprocess=euler.sim(rproc,delta.t=1/365.25),
       initializer=initlz,
       dmeasure=dmeas,
       rmeasure=rmeas,
       covar=covar,
       tcovar="time",
       zeronames=c("C","W"),
       statenames=c("S","E","I","R","C","W"),
       paramnames=c("R0","mu","sigma","gamma","alpha","iota",
                    "rho","sigmaSE","psi","cohort","amplitude",
                    "S_0","E_0","I_0","R_0")
  ) -> m1

The following codes plot the data and covariates together.

m1 %>% as.data.frame() %>% 
  melt(id="time") %>%
  ggplot(aes(x=time,y=value))+
  geom_line()+
  facet_grid(variable~.,scales="free_y")

He et al. (2010) estimated the parameters of this model. The full set is included in the R code accompanying this document, where they are read into a data frame called mles.

mles %>% subset(town=="London") -> mle
paramnames <- c("R0","mu","sigma","gamma","alpha","iota",
                "rho","sigmaSE","psi","cohort","amplitude",
                "S_0","E_0","I_0","R_0")
mle %>% extract(paramnames) %>% unlist() -> theta
mle %>% subset(select=-c(S_0,E_0,I_0,R_0)) %>%
  knitr::kable(row.names=FALSE)

Verify that we get the same likelihood as He et al. (2010).

library(foreach)
library(doParallel)

registerDoParallel()

set.seed(998468235L,kind="L'Ecuyer")

foreach(i=1:4,
        .packages="pomp",
        .options.multicore=list(set.seed=TRUE)
) %dopar% {
  pfilter(m1,Np=10000,params=theta)
} -> pfs
logmeanexp(sapply(pfs,logLik),se=TRUE)

##                          se 
## -3800.9542945     0.3404236

Simulations at the MLE:

m1 %>% 
  simulate(params=theta,nsim=9,as.data.frame=TRUE,include.data=TRUE) %>%
  ggplot(aes(x=time,y=cases,group=sim,color=(sim=="data")))+
  guides(color=FALSE)+
  geom_line()+facet_wrap(~sim,ncol=2)

Parameter transformations

The parameters are constrained to be positive, and some of them are constrained to lie between $0$ and $1$. We can turn the likelihood maximization problem into an unconstrained maximization problem by transforming the parameters. The following Csnippets implement such a transformation and its inverse.

toEst <- Csnippet("
  Tmu = log(mu);
  Tsigma = log(sigma);
  Tgamma = log(gamma);
  Talpha = log(alpha);
  Tiota = log(iota);
  Trho = logit(rho);
  Tcohort = logit(cohort);
  Tamplitude = logit(amplitude);
  TsigmaSE = log(sigmaSE);
  Tpsi = log(psi);
  TR0 = log(R0);
  to_log_barycentric (&TS_0, &S_0, 4);
")

fromEst <- Csnippet("
  Tmu = exp(mu);
  Tsigma = exp(sigma);
  Tgamma = exp(gamma);
  Talpha = exp(alpha);
  Tiota = exp(iota);
  Trho = expit(rho);
  Tcohort = expit(cohort);
  Tamplitude = expit(amplitude);
  TsigmaSE = exp(sigmaSE);
  Tpsi = exp(psi);
  TR0 = exp(R0);
  from_log_barycentric (&TS_0, &S_0, 4);
")


pomp(m1,toEstimationScale=toEst,
     fromEstimationScale=fromEst,
     statenames=c("S","E","I","R","C","W"),
     paramnames=c("R0","mu","sigma","gamma","alpha","iota",
                  "rho","sigmaSE","psi","cohort","amplitude",
                  "S_0","E_0","I_0","R_0")) -> m1

Imported infections

$$\text{force of infection} = \mu_{SE}=\frac{\beta(t)}{N(t)}\,(I+\iota)\,\zeta(t)$$

Seasonality

Cohort effect

Birth delay

Profile likelihood for birth-cohort delay, showing 95% and 99% critical values of the log likelihood.

Reporting rate

Predicted vs observed critical community size

$R_0$

• Recall that $R_0$ is the basic reproduction number: a measure of how communicable an infection is.
• Existing estimates of $R_0$ (c. 15–20) come from two sources:
• serology surveys
• models fit to data using feature-based methods

Parameter estimates

$r = \mathrm{cor}(\log{\hat\theta},\log{N_{1950}})$

Extrademographic stochasticity

$$\mu_{SE}=\frac{\beta(t)}{N(t)}\,(I+\iota)\,\zeta(t)$$