In this vignette, we will demonstrate how to implement the Particle Marginal Metropolis-Hastings (PMMH) algorithm (which has been implemented in the pomp package) for analyzing prevalence counts from an outbreak of influenza in a British boarding school. This dataset was pulled from the pomp package, and was described and analyzed in Fintzi, Cui, Wakefield, and Minin (2016).
We begin by read and plotting in the dataset.
library(pomp) library(ggplot2) library(gtools) set.seed(52787) st.date <- as.Date("1978/1/22") en.date <- as.Date("1978/2/4") date.seq <- seq(st.date,en.date, by = "day") bbs <- data.frame(time = rep(date.seq,1), count = c(rep("Observed count", 14)), I = c(3,8,28,76,222,293,257,237,192,126,70,28,12,5))
We use the SIR model to model the British Boarding data here with a negative binomial sampling distribution. The hosts are divided into three classes, according to their status. The susceptible class (S) contains those that have not yet been infected and are thereby still susceptible to it; the infected class (I) comprises those who are currently infected and, by assumption, infectious; the removed class (R) includes those who are recovered or quarantined as a result of the infection. Individuals in R are assumed to be immune against reinfection. It is natural to formulate this model as a continuous-time Markov process. To start with, we need to specify the measurement model.
rmeas <- " cases=rnbinom_mu(exp(phi),exp(rho)*I/(1+exp(rho))); //represent the data " dmeas<-" lik=dnbinom_mu(cases,exp(phi),exp(rho)*I/(1+exp(rho)),give_log); //return the loglikelihood "
Here, we are using the $\textbf{cases}$ to refer to the data (number of reported cases) and $\textbf{I}$ to refer to the true incidence over the reporting interval. The negative binomial simulator rnbinom_mu and density function dnbinom_mu are provided by R. Notice that, in these snippets, we never declare the variables; pomp will ensure tht the state variable ($\textbf{I}$), observables ($\textbf{cases}$), parameters ($\textbf{rho}$, $\textbf{phi}$) and likelihood ($\textbf{lik}$) are defined in the contexts within which these snippets are executed.
And for the transition of the SIR model, we use a tau-leaping algorithm, one version of which is implemented in the pomp package via the $\bf{euler.sim}$ plug-in. The algorithm holds the transition rates constant over a small interval of time and simulates the numbers of transitions that occur over the interval. The functions $\textbf{reulermultinom}$ draw random deviates from such distributions. Then we need to first specify a function that advances the states from $t$ to $t+\triangle t$, and we give the trantition of the SIR model code below:
sir.step<-" double rate[2]; double dN[2]; rate[0]=exp(bet)*I; //Infection rate rate[1]=exp(mu); //recovery rate reulermultinom(1,S,&rate[0],dt,&dN[0]); //generate the number of newly infected people reulermultinom(1,I,&rate[1],dt,&dN[1]); //generate the number of newly recovered people if(!R_FINITE(S)) Rprintf(\"%lg %lg %lg %lg %lg %lg %lg %lg %lg\\n\",dN[0],rate[0],dN[1],rate[1],bet,mu,S,I,R); S+=-dN[0]; //update the number of Susceptible I+=dN[0]-dN[1]; //update the number of Infection R+=dN[1]; //update the number of Recovery "
The pomp will ensure that the undeclared state variables and parameters are defined in the context within which the snippet is executed. And the $\textbf{rate}$ and $\textbf{dN}$ arrays hold the rates and numbers of transition events, respectively. Then we could construct the pomp objects below:
sir <- pomp( data = data.frame(cases = bbs[,3], time = seq(0, 13, by = 1)), #"cases" is the dataset, "time" is the observation time times = "time", t0 = 0, #initial time point dmeasure = Csnippet(dmeas), rmeasure = Csnippet(rmeas), #return the likelihood rprocess = euler.sim(step.fun = Csnippet(sir.step), delta.t = 1/12), statenames = c("S", "I", "R"), #state space variable name paramnames = c("bet", "mu", "rho", "theta1", "theta2", "phi"), #parameters name initializer = function(params, t0, ...) { ps <-exp(params["theta1"]) / (1 + exp(params["theta1"]) + exp(params["theta2"])) #given the initial proportion of susceptible pi <- 1 / (1 + exp(params["theta1"]) + exp(params["theta2"])) #given the initial proportion of infection pr <-exp(params["theta2"]) / (1 + exp(params["theta1"]) + exp(params["theta2"])) #given the initial proportion of recovery return(setNames(as.numeric(rmultinom( 1, 763, prob = c(ps, pi, pr) )), c("S", "I", "R"))) }, params = c( #for data generation(We didn't generate data here, but pomp need this part.) bet = -5, mu = -1, rho = 3, theta1 = 5, theta2 = 1, phi = 2 ) )
To carry out Bayesian inference we need to specify a prior distribution on unknown parame- ters. The pomp constructor function provides the rprior and dprior arguments, which can be filled with functions that simulate from and evaluate the prior density, respectively. Methods based on random walk Metropolis-Hastings require evaluation of the prior density (dprior), so we specify dprior for the SIR model as follows.
sir.dprior <- function(params, ..., log) { f <- ( dgamma( exp(params[1]), shape = 0.001, rate = 1, log = TRUE ) + params[1] + #log prior for log(infection rate) "beta" dgamma( exp(params[2]), shape = 1, rate = 2, log = TRUE ) + params[2] + #log prior for log(recovery rate) "mu" dbeta(exp(params[3]) / (1 + exp(params[3])), 2, 1, log = TRUE) + params[3] - log((1 + exp(params[3])) ^ 2) + #log prior for logit(sampling probablity) "rho" log(ddirichlet(c( exp(params[4]) / (1 + exp(params[4]) + exp(params[5])), 1 / (1 + exp(params[4]) + exp(params[5])), exp(params[5]) / (1 + exp(params[4]) + exp(params[5])) ), c(900, 3, 9))) + params[4] + params[5] - 3 * log(1 + exp(params[4]) + exp(params[5])) + #log prior for logit(initial value) dgamma(exp(params[6]), 1, 0.1, log = TRUE) + params[6] #log prior for log(negative size parameter) "phi" ) if (log) f else exp(f) }
And starting value of the PMMH is specified below:
param.initial <- c( #initial value of the parameters of the PMCMC bet = -6, mu = -1, rho = 2.5, theta1 = 6, theta2 = 1, phi = 2 )
And the following runs 1 PMMH chain for 2000 tunning steps (using 500 particles) to get a suitable multivariate normal random walk proposal diagonal variance-covariance matrix value.
pmcmc1 <- pmcmc( pomp(sir, dprior = sir.dprior), #given the prior function start = param.initial, #given the initial value of the parameters Nmcmc = 10, #number of mcmc steps Np = 500, max.fail = Inf, proposal = mvn.rw.adaptive(0.1 * c( bet = 0.1, #sampling variance for beta mu = 0.1, #sampling variance for mu rho = 0.3, #sampling variance for rho theta1 = 0.2, #sampling variance for theta1 theta2 = 0.2, #sampling variance for theta2 phi = 0.5 #sampling variance for phi ), shape.start = 100, scale.start = 100, max.scaling = 2) )
And the following runs 1 PMMH chain for 100000 steps (using 500 particles) with suitable multivariate normal random walk proposal diagonal variance-covariance matrix value we get above to get the posterior distribution of the interested parameters.
start_time <- Sys.time(); #calculation of time pmcmc1 <- pmcmc( pmcmc1, #given the prior function start = param.initial, #given the initial value of the parameters Nmcmc = 10, max.fail = Inf, proposal = mvn.rw(covmat(pmcmc1)) ) end_time <- Sys.time(); run_time <- difftime(end_time, start_time, units = "hours") pomp_results <- list(time = run_time, results = pmcmc1)
Here, we use the SEIR model to model the British Boarding data here with a negative binomial sampling distribution. The hosts are divided into four classes, according to their status. The susceptible class (S) contains those that have not yet been infected and are thereby still susceptible to it; the exposure class (E) contains those that have not yet been infected but already exposed to it; the infected class (I) comprises those who are currently infected and, by assumption, infectious; the removed class (R) includes those who are recovered or quarantined as a result of the infection. Individuals in R are assumed to be immune against reinfection. It is natural to formulate this model as a continuous-time Markov process. To start with, we need to specify the measurement model.
rmeas <- " cases=rnbinom_mu(exp(phi),exp(rho)*I/(1+exp(rho))); //represent the data " dmeas<-" lik=dnbinom_mu(cases,exp(phi),exp(rho)*I/(1+exp(rho)),give_log); //return the loglikelihood "
Here, we are using the $\textbf{cases}$ to refer to the data (number of reported cases) and $\textbf{I}$ to refer to the true incidence over the reporting interval. The binomial simulator rbinom_mu and density function dbinom_mu are provided by R. Notice that, in these snippets, we never declare the variables; pomp will ensure tht the state variable ($\textbf{I}$), observables ($\textbf{cases}$), parameters ($\textbf{rho}$) and likelihood ($\textbf{lik}$) are defined in the contexts within which these snippets are executed.
And for the transition of the SEIR model, we use a tau-leaping algorithm, one version of which is implemented in the pomp package via the $\bf{euler.sim}$ plug-in. The algorithm holds the transition rates constant over a small interval of time and simulates the numbers of transitions that occur over the interval. The functions $\textbf{reulermultinom}$ draw random deviates from such distributions. Then we need to first specify a function that advances the states from $t$ to $t+\triangle t$, and we give the trantition of the SEIR model code below:
seir.step<-" double rate[3]; double dN[3]; rate[0]=exp(beta)*I; //Infection rate[1]=exp(gamma); //Rate of S to E rate[2]=exp(mu); //Rate of E to I reulermultinom(1,S,&rate[0],dt,&dN[0]); //generate the number of newly from S to E reulermultinom(1,E,&rate[1],dt,&dN[1]); //generate the number of newly from E to I reulermultinom(1,I,&rate[2],dt,&dN[2]); //generate the number of newly from I to R if(!R_FINITE(S)) Rprintf(\"%lg %lg %lg %lg %lg %lg %lg %lg %lg\\n\",dN[0],rate[0],dN[1],rate[1],beta,mu,S,I,R); S+=-dN[0]; //update the number of Susceptible E+=dN[0]-dN[1]; //update the number of E I+=dN[1]-dN[2]; //update the number of I R+=dN[2]; //update the number of R "
The pomp will ensure that the undeclared state variables and parameters are defined in the context within which the snippet is executed. And the $\textbf{rate}$ and $\textbf{dN}$ arrays hold the rates and numbers of transition events, respectively. Then we could construct the pomp objects below:
seir <- pomp( data = data.frame(cases = bbs[,3], time = seq(0, 13, by = 1)), #"cases" is the dataset, "time" is the observation time times = "time", t0 = 0, #initial time point dmeasure = Csnippet(dmeas), rmeasure = Csnippet(rmeas), #return the likelihood rprocess = euler.sim(step.fun = Csnippet(seir.step), delta.t = 1/12), #delta.t is here statenames = c("S", "E", "I", "R"), #state space variable name paramnames = c("beta", "gamma", "mu", "rho", "theta1", "theta2", "theta3", "phi"), #parameters name initializer = function(params, t0, ...) { #Initial proportion of S, E, I, R ps <- exp(params["theta1"]) / (1 + exp(params["theta1"]) + exp(params["theta2"]) + exp(params["theta3"])) pe <- exp(params["theta2"]) / (1 + exp(params["theta1"]) + exp(params["theta2"]) + exp(params["theta3"])) pi <- 1 / (1 + exp(params["theta1"]) + exp(params["theta2"]) + exp(params["theta3"])) pr <- exp(params["theta3"]) / (1 + exp(params["theta1"]) + exp(params["theta2"]) + exp(params["theta3"])) return(setNames(as.numeric(rmultinom(1, 763, prob = c(ps, pe, pi, pr))), c("S", "E", "I", "R"))) }, params = c( #for data generation(We didn't generate data here, but pomp need this part.) beta = -6, gamma = 1.5, mu = -1, rho = 2.5, theta1 = 6, theta2 = 1, theta3 = 1, phi = 2 ) )
To carry out Bayesian inference we need to specify a prior distribution on unknown parame- ters. The pomp constructor function provides the rprior and dprior arguments, which can be filled with functions that simulate from and evaluate the prior density, respectively. Methods based on random walk Metropolis-Hastings require evaluation of the prior density (dprior), so we specify dprior for the SEIR model as follows.
trans<-function(a,b,c){ a_t<-exp(a) b_t<-exp(b) c_t<-exp(c) return(a_t*b_t*c_t/((1+a_t+b_t+c_t)^4)) } seir.dprior <- function(params, ..., log) { f <- (dgamma(exp(params[1]), 0.001, 1, log = TRUE) + params[1] + #log prior for beta dgamma(exp(params[2]), 0.001, 1, log = TRUE) + params[2] + #log prior for gamma dgamma(exp(params[3]), 1, 2, log = TRUE) + params[3] + #log prior for mu dbeta(exp(params[4]) / (1 + exp(params[4])), 2, 1, log = TRUE) + params[4] - log((1 + exp(params[4])) ^ 2) + #log prior for rho log(ddirichlet(c(exp(params[5]) / (1 + exp(params[5]) + exp(params[6]) + exp(params[7])), exp(params[6]) / (1 + exp(params[5]) + exp(params[6]) + exp(params[7])), 1 / (1 + exp(params[5]) + exp(params[6]) + exp(params[7])), exp(params[7]) / (1 + exp(params[5]) + exp(params[6]) + exp(params[7]))), c(900, 6, 3, 9)) * trans(params[5], params[6], params[7])) + # log prior for p_t1 dgamma(exp(params[6]), 1, 0.1, log = TRUE) + params[6] #log prior for log(negative size parameter) "phi" ) if (log) { f } else { exp(f) } }
And starting value of the PMMH is specified below:
param.initial <- c( #initial value of the parameters of the PMCMC beta = -6, gamma = 2, mu = -1, rho = 2.5, theta1 = 6, theta2 = 1, theta3 = 1, phi = 2 )
And the following runs 1 PMMH chain for 2000 tunning steps (using 5000 particles) to get a suitable multivariate normal random walk proposal diagonal variance-covariance matrix value.
pmcmc1 <- pmcmc( pomp(seir, dprior = seir.dprior), #given the prior function start = param.initial, #given the initial value of the parameters Nmcmc = 2000, #number of mcmc steps Np = 10, max.fail = Inf, proposal = mvn.rw.adaptive(0.1 * c( beta = 0.1, #sampling variance for beta gamma = 0.1, # proposal variance for gamma mu = 0.1, #sampling variance for mu rho = 0.3, #sampling variance for rho theta1 = 0.2, #sampling variance for theta1 theta2 = 0.2, #sampling variance for theta2 theta3 = 0.2, phi = 0.5 #sampling variance for phi ), shape.start = 100, scale.start = 100, max.scaling = 1.2) )
And the following runs 1 PMMH chain for 100000 steps (using 5000 particles) with suitable multivariate normal random walk proposal diagonal variance-covariance matrix value we get above to get the posterior distribution of the interested parameters.
start_time <- Sys.time(); #calculation of time pmcmc1 <- pmcmc( pmcmc1, #given the prior function start = param.initial, #given the initial value of the parameters Nmcmc = 10, max.fail = Inf, proposal = mvn.rw(covmat(pmcmc1)) ) end_time <- Sys.time(); run_time <- difftime(end_time, start_time, units = "hours") pomp_results <- list(time = run_time, results = pmcmc1)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.