1928-1931 Detroit Mumps Infection Analysis
2018/4/25
1 Introduction
Mumps is a worldwide common disease, especially for children.1 It caused by the mumps virus. Mumps is highly contagious and spreads rapidly among people living in close quarters. The virus is transmitted by respiratory droplets or direct contact with an infected person. People are infectious to each other from about seven days before the start of symptoms to about eight days after.2 Symptoms typically occur 16 to 18 days after exposure and resolve after seven to ten days. About a third of people have mild or no symptoms. Without immunization about 0.1 percent to one percent of the population are affected per year. These characteristics give mumps great space to spread before the introduction of its vaccination. This vaccination is extensively used since the start of U.S. mumps vaccination program in 1967.3
I will analyse the process of mumps transmission in Detroit, Michigan, from January 1928 to December 1931, which is earlier than the vaccination coming out. The data is gained from Project Tycho freely.4 The data is the weekly count of mumps infected people. There are two week count missing in the data, ‘1929.09.22-1929.09.28’ and ‘1930.12.21-1930.12.27’. I just fill them with the average of their neighbouring week counts. I build a SEIR model and use the partial observed markov process to estimate parameters. The SEIR model code is based on the case study: Measles in large and small towns, Aaron King5 and He et al6.
2 Data Exploratory Analysis
From the line plot of weekly cases, we can notice that the time series show a strong seasonality with 1 year period. The infection peaks always happen in spring. Summer is always the vales. And the number of infected people would continuously increase from fall to winter. This phenomenon may be caused by the University of Michigan that locates nearby Detroit. During summer vacations, most students go home and make Detroit less crowded, which brings the decreasing of infected people number. The decomposition of cases shows the trend more clearly.
3 Covariates
According to Aaron King, the SEIR model has the prior distribution related to two covariates, the population size and the birth rate in Detroit. I find them from the Michigan Department of Health & Human Services website.7 The population is increasing over 1928-1931, while the birth rate is decreasing from 1930.
Both the data of population size and birth rate are collected by year. To fit a partially observed Markov process model more precisely, the demographical information is smoothed by the spline method.
4 The Partially Observed Markov Process Model
4.1 The Process Model
Model diagram:
\(b = \text{births}\)
\(S = \text{susceptibles}\)
\(E = \text{exposed, non-infectious, incubating}\)
\(I = \text{Infectious}\)
\(R = \text{Recovered}\)
Transmission from B -> S:
\[\mu_{BS}=(1-c)B(t-\tau)+c\delta(t-t_0)\int_{t-1}^tB(t-\tau-s)ds\]
- \(B(t)=\text{birth rate}\)
- \(N(t)=\text{population size}\)
- \(c=\text{cohort effect}\)
- \(\tau=\text{entry delay}\)
Like the measles case study8, \(\mu_{BS}\) is the transmission from new birth babies to susceptible population. The cohort effect and the entry delay should be taken into account. Since the vaccination of mumps did not exist from 1928 to 1932, I consider that all new births enter the susceptibles.
Transmission from S -> E:
\[\mu_{SE}(t)=\frac{\beta(t)}{N(t)}(I+\iota)\zeta(t)\]
\[ \beta(t) = \left\{ \begin{array}{ll} \beta_0(1+a) & \textrm{during term}\\ \beta_0(1-a) & \textrm{during vacation} \end{array} \right. \]
- To get the vacation date, I just use University of Michigan 2017 fall to 2018 summer schedule.
- \(\iota=\text{imported infections}\)
- \(\zeta(t)=\text{Gamma white noise with intensity }\sigma_{SE} \text{(He et al. 2010, Bhadra et al. 2011)}\)
- \(\text{Overdispersed binomial measurement model: cases}_t|\triangle N_{IR}=z_t \sim N(\rho z_t,\rho (1-\rho)z_t+(\psi \rho z_t)^2)\)
4.2 The Measurement Model
The assumptions are \(E(cases|C)=\rho C\) where C is the true incidence and \(0<\rho<1\) is the reporting efficiency and \(Var[cases|C]=\rho (1-\rho)C+(\psi\rho C)^2\), where \(\psi\) quantifies overdispersion. Then cases can be chosen by:
\[cases|C \sim f(c|\rho,\psi,C)=\Phi(c+\frac{1}{2},\rho C,\rho(1-\rho)C+(\psi \rho C)^2)-\Phi(c-\frac{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 code is to calculate \(P(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;
}
")5 Parameters Estimating
5.1 Parameters Annotation
The pomp object is created by the following code.
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,
toEstimationScale=toEst,
fromEstimationScale=fromEst,
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")
) -> m1Then, I would explain the meaning of each parameters:
- \(R0\): This parameter corresponds with the transmission rate. The higher \(R0\) is, the higher transmission rate is.
- \(\mu\): The death rate.
- \(\gamma\): The rate of recovery.
- \(\alpha\): The proportion of mixture. The closer to 0.5 \(\alpha\) is, the more balanced the mixture is.
- \(\iota\): imported infections.
- \(\rho\): The reporting rate.
- \(\psi\): Overdispersion parameter in the reporting process.
- \(\sigma\): The rate of ending the latent stage.
- \(\sigma_{SE}\): The extra-demographic stochasticity.9
- \(cohort\): Cohort effect.
- \(amplitude\): The seasonality parameter
5.2 Local Searching
First, I run a local searching to investigate the tendency of likelihood changing according to each parameter. The initial point is the maximum likelihood estimators of parameters in an assistant large range global searching. Because the death rate varies little in 1928-1932, I set \(\mu\) to a fixed number 1%. This data is found from the Michigan Department of Health & Human Services website.10
run_level=3
switch(run_level,
{mumps_Np=100; mumps_Nmif=10; mumps_Neval=2; mumps_Nglobal=10; mumps_Nlocal=10; mumps_Nsim=50},
{mumps_Np=500; mumps_Nmif=100; mumps_Neval=10; mumps_Nglobal=20; mumps_Nlocal=20; mumps_Nsim=100},
{mumps_Np=1000; mumps_Nmif=300; mumps_Neval=20; mumps_Nglobal=100; mumps_Nlocal=40; mumps_Nsim=500}
)
mumps_rw.sd <- 0.02
mumps_rw.sd.R0 <- 0.1
mumps_rw.sd.simga <- 0.1
mumps_cooling.fraction.50 <- 0.5
init_param = c(R0=90.4585,gamma=25.0531,alpha=0.7964782,iota=0.05067204,rho=0.986742,psi=0.3063079,sigma=23.34644,sigmaSE=0.3221183,cohort=0.05082893,amplitude=0.2669821,mu=0.01,S_0=0.004984235,E_0=2.779469e-05,I_0=2.174464e-05,R_0=0.9949662)stew(file=sprintf("/Users/mayumeng/Downloads/STATS_531/local_search-%d.rda",run_level),{
t_local <- system.time({
mifs_local <- foreach(i=1:mumps_Nlocal,
.packages='pomp',
.combine=c,
.options.multicore=mcopts) %dopar% {
mif2(
m1,
start=init_param,
Np=mumps_Np,
Nmif= mumps_Nmif,
cooling.type="geometric",
cooling.fraction.50=mumps_cooling.fraction.50,
transform=TRUE,
rw.sd=rw.sd(
R0=mumps_rw.sd.R0,
gamma=mumps_rw.sd,
alpha=mumps_rw.sd,
iota=mumps_rw.sd,
rho=mumps_rw.sd,
psi=mumps_rw.sd,
sigma=mumps_rw.sd.simga,
sigmaSE=mumps_rw.sd,
cohort=mumps_rw.sd,
amplitude=mumps_rw.sd,
mu=0,
S_0=mumps_rw.sd,
E_0=mumps_rw.sd,
I_0=mumps_rw.sd,
R_0=mumps_rw.sd
)
)
}
})
},seed=900242057,kind="L'Ecuyer")
stew(file=sprintf("/Users/mayumeng/Downloads/STATS_531/lik_local_eval-%d.rda",run_level),{
t_local_eval <- system.time({
liks_local <- foreach(i=1:mumps_Nlocal,.packages='pomp',.combine=rbind, .options.multicore=mcopts) %dopar% {
evals <- replicate(mumps_Neval, logLik(pfilter(m1,params=coef(mifs_local[[i]]),Np=mumps_Np)))
logmeanexp(evals, se=TRUE)
}
})
},seed=442141592,kind="L'Ecuyer")## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -770.7 -760.0 -757.0 -758.1 -754.6 -751.2This investigation took 28 minutes for the maximization and 1 minutes for the likelihood evaluation. From the point plot below, we can see the pairwise relationship between parameters and the likelihood values. Only the largerest likelihood values are present. Then I choose the suitable range intervals for each parameters for the global searching.
From the superimposed convergence diagnostic plots for multiple Monte Carlo replications of the maximization procedure, I notice the effective sample size appears to be good, but \(\rho\), \(\sigma_{SE}\), cohort and amplitude don’t converge. Their value range need to be adjusted.
5.3 Global Searching
There are parameters value range I use for global searching. The global investigation lasts for 2 hours and 22 minutes. S_0, E_0, I_0, R_0 are fixed to the MLE in local searching.
mumps_box <- rbind(
R0=c(20,40),
gamma=c(0,70),
alpha=c(.5,1),
iota=c(-0.5,.5),
rho=c(.4,1),
psi=c(.15,.5),
sigma=c(0,40),
sigmaSE=c(0.01,.5),
cohort=c(0,1),
amplitude=c(0.01, 0.5)
)
mumps_fixed_params = c(mu=.01, S_0=.103, E_0=6.62e-05, I_0=3.38e-06, R_0=.897)stew(file=sprintf("/Users/mayumeng/Downloads/STATS_531/box_eval-%d.rda",run_level),{
t_global <- system.time({
mifs_global <- foreach(i=1:mumps_Nglobal,.packages='pomp', .combine=c, .options.multicore=mcopts) %dopar%
mif2(
m1,
start=c(apply(mumps_box,1,function(x)runif(1,x[1],x[2])),mu=.01),
Np=mumps_Np,
Nmif=mumps_Nmif,
cooling.type="geometric",
cooling.fraction.50=0.2,
transform=TRUE,
rw.sd=rw.sd(
R0=0.02,
sigma=0.02,
gamma=0.02,
alpha=0.02,
iota=0.02,
rho=0.02,
sigmaSE=0.02,
psi=0.02,
cohort=0.02,
amplitude=0.02,
S_0=0.02,
E_0=0.02,
I_0=0.02,
R_0=0.02
)
)
})
},seed=1270401374,kind="L'Ecuyer")
stew(file=sprintf("/Users/mayumeng/Downloads/STATS_531/lik_global_eval-%d.rda",run_level),{
t_global_eval <- system.time({
liks_global <- foreach(i=1:mumps_Nglobal,.packages='pomp',.combine=rbind, .options.multicore=mcopts) %dopar% {
evals <- replicate(mumps_Neval, logLik(pfilter(m1,params=coef(mifs_global[[i]]),Np=mumps_Np)))
logmeanexp(evals, se=TRUE)
}
})
},seed=442141592,kind="L'Ecuyer")## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1032.6 -772.1 -769.4 -775.1 -767.2 -765.1
- In all these iteration, the maximum value of the log-likelihood is as large as -765.1. And the log-likelihood values also show a good converge in the diagnostic plot.
- From the diagnostic plot we can see that among all the 15 parameters, \(\alpha\), \(cohort\), \(amplitude\) converge badly, while \(\psi\) converge well. To get a resonable computing time, I just fix 5 parameters. I believe fixing \(\mu\) is a good trick, but S_0, E_0, I_0, R_0’s values should be explored more carefully. To improve the MLE paramters’, we still need more iterations to make them converge.
- When it comes to the effective sample size plot, the values are unstable. The effective sample size is very small at the beginning (1928), then increasing quickly with the time. It drops a little bit around 1932. Increasing the number of particles may improve the performance. Overall, the effective sample size performs fine.
6 Simulation
The following code is to simulate a time series with a set of given parameters. Then I do the simulation with the maximum likelihood estimator in global searching.
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
")
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;
")
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;
}
")
- From the 5 simulated time series I find that this model will over estimate the case number seriously after 1931. I restrict the y-axis to [0,100] in the plot. The problem is probably caused by the lack of iteration. The three run level only get (100,500,2000) iterations. The ideal iteration level is at least (1000,5000,20000) in my mind. But the model has too many parameters to finish computing in a reasonable time.
- Then I try to use local searching MLE to simulate for the data. The results is better than global searching. But the simulated time series lack of obvious seasonality. They mainly show the first, third and fourth peaks, but ignore the second peak in the original data. It is probably caused by not converging of the parameter “amplitude”. The blue shadow in the second plot below is the area covered by several simulated time series, intuitively like the confidence interval of the simulation.
7 Conclusion
There are three conclusion that can be demonstrated in this report:
- The SEIR model is suitable for analysing the infection of mumps overall. It can reproduce the first, third and fourth peaks in the original data. The second peak is truely weaker than other three peaks. So with more iteration, the SEIR model will probably fit the original data pretty well.
- The MLE of parameters are around:
## R0 gamma alpha iota rho psi sigma
## result.37 22.27849 37.03844 0.700264 1.026439 0.8578209 0.2669429 26.37818
## sigmaSE cohort amplitude mu S_0 E_0
## result.37 0.214664 0.172121 0.005685944 0.01 0.1031331 6.619662e-05
## I_0 R_0
## result.37 3.384217e-06 0.8967973- We need to explore the S_0, E_0, I_0 and R_0 more carefully, instead of fixing them to some unsuitable values.
8 Reference
Wikipedia, Mumps, https://en.wikipedia.org/wiki/Mumps↩
Centers of Diease Control and Prevention, Mumps Vaccination Introduction, https://www.cdc.gov/mumps/vaccination.html↩
Project Tycho, https://www.tycho.pitt.edu/dataset/US.36989005/↩
King, A. (2015, June 7). Case study: Measles in large and small towns. Retrieved from http://kingaa.github.io/sbied/measles/measles.html↩
He, D., Ionides, E. L., & King, A. A. (2010). Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface, 7(43), 271-283. http://doi.org/10.1098/rsif.2009.0151↩
Michigan Department of Health & Human Services https://www.mdch.state.mi.us/osr/natality/tab4.1.asp↩
King, A. (2015, June 7). Case study: Measles in large and small towns. Retrieved from http://kingaa.github.io/sbied/measles/measles.html↩
He, D., Ionides, E. L., & King, A. A. (2010). Plug-and-play inference for disease dynamics: measles in large and small populations as a case study. Journal of the Royal Society Interface, 7(43), 271-283. http://doi.org/10.1098/rsif.2009.0151↩
Michigan Department of Health & Human Services https://www.mdch.state.mi.us/pha/osr/deaths/g21.asp↩