In fintzij/BDAepimodel: Fit Stochastic Epidemic Models in R via Bayesian Data Augmentation
library(BDAepimodel)
library(coda)
library(Rcpp)
library(pomp)
library(gtools)
library(ggplot2)
Analyzing an epidemic with SIRS dynamics
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 simulated prevalence counts. We fit SIRS models to binomially distributed prevalence counts.
We simulated an epidemic with SIRS dynamics in a population of 200 individuals, 1% of whom were initially infected, with the rest being suscpetible. The per–contact infectivity rate was $\beta = 0.0009$, the average infectious period duration was $1/\mu = 14$ days, and the average time until loss of immunity was $1/\gamma = 150$ days, which together correspond to a basic reproduction number of $R_0 = \beta N / \mu = 2.52$. Binomially sampled prevalence was observed at one week intervals over a period of one year with sampling probability $\rho=0.8$.
The data simulation code is in "SIR_SEIR_SIRS" vignettes. And the data is shown below.
data<-c(1,4,6,6,4,9,12,23,29,32,29,48,32,22,15,13,8,7,3,2,1,0,0,1,3,1,0,1,3,3,3,1,5,3,5,3,2,1,0,6,2,3,6,3,6,5,5,8,15,17,14,8,8)
Using Approximate Tau-leaping Algorithm:
We use the SIRS model to model the dataset. 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 lose immunity after a few days of the recovery. 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=rbinom(I,exp(rho)/(1+exp(rho)));//represent the data
"
dmeas<-"
lik=dbinom(cases,I,exp(rho)/(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 and density function dbinom 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 SIRS model, we use a approximating 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 SIRS model code below:
sirs.step<-"
double rate[3];
double dN[3];
rate[0]=exp(beta)*I; //Infection rate
rate[1]=exp(mu); //Recovery rate
rate[2]=exp(gamma); //Loss of immunity 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 infected people
reulermultinom(1,R,&rate[2],dt,&dN[2]); //generate the number of people who newly lose immunity
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]+dN[2]; //update the number of Susceptible
I+=dN[0]-dN[1]; //update the number of Infection
R+=dN[1]-dN[2]; //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:
sirs <- pomp(
data = data.frame(cases = data, time = seq(1, 365, by = 7)), #"cases" is the dataset, "time" is the observation time
times = "time",
t0 = 1, #initial time point
dmeasure = Csnippet(dmeas),
rmeasure = Csnippet(rmeas), #return the likelihood
rprocess = euler.sim(step.fun = Csnippet(sirs.step), delta.t = 1/3), #delta.t is here
statenames = c("S", "I", "R"), #state space variable name
paramnames = c("beta", "mu", "gamma", "rho", "theta1", "theta2"), #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, 200, prob = c(ps, pi, pr))), c("S", "I", "R")))
},
# not used, but pomp needs this
params = c(
beta = log(0.0001),
mu = log(1 / 14),
gamma = log(1 / 140),
rho = 1.5,
theta1 = log(49),
theta2 = log(0.00000001)
)
)
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 SIRS model as follows.
sirs.dprior <- function(params, ..., log) {
f <- (
dgamma(exp(params[1]), 0.1, 100, log = TRUE) + params[1] + #log prior for log(infection rate) "beta"
dgamma(exp(params[2]), 1.8, 14, log = TRUE) + params[2] + #log prior for log(recovery rate) "mu"
dgamma(exp(params[3]), 0.0625, 10, log = TRUE) + params[3] + #log prior for log(loss of susceptible) "gamma"
dbeta(exp(params[4]) / (1 + exp(params[4])), 5, 1, log = TRUE) +
params[4] - log((1 + exp(params[4])) ^ 2) + #log prior for logit(sampling probablity) "rho"
log(ddirichlet(c(exp(params[5]) / (1 + exp(params[5]) + exp(params[6])),
1 / (1 + exp(params[5]) + exp(params[6])),
exp(params[6]) / (1 + exp(params[5]) + exp(params[6]))),
c(9, 0.15, 0.001))) +
params[5] + params[6] - 3 * log(1 + exp(params[5]) + exp(params[6])) #log prior for logit(initial value)
)
if (log) {
f
} else {
exp(f)
}
}
And starting value of the PMMH is specified below:
param.initial <- c(
beta = log(abs(rnorm(1, 0.00077, 1e-4))),
mu = log(abs(rnorm(1, 1/18, 1e-4))),
gamma = log(abs(rnorm(1, 1/140, 1e-4))),
rho = 1.5,
theta1 = 3,
theta2 = -5
)
And the following runs 1 PMMH chain for 5000 tunning steps (using 500 particles) to get a suitable multivariate normal random walk proposal diagonal variance-covariance matrix value.
pmcmc1 <- pmcmc(
pomp(sirs, dprior = sirs.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(
beta = 0.1, #sampling variance for beta
mu = 0.1, #sampling variance for mu
gamma = 0.1, #sampling variance for gamma
rho = 0.1, #sampling variance for rho
theta1 = 0.1,#sampling variance for theta1
theta2 = 0.1 #sampling variance for theta2
),
scale.start = 100,
shape.start = 100)
)
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(); #calculation of time
run_time <- difftime(end_time, start_time, units = "hours") #calculation of time
pomp_results <- list(time = run_time, results = pmcmc1)
Using Exact Gillespie Algorithm:
We use the SIRS model to model the dataset. 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 lose immunity after a few days of the recovery. 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=rbinom(I,exp(rho)/(1+exp(rho)));//represent the data
"
dmeas<-"
lik=dbinom(cases,I,exp(rho)/(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 and density function dbinom 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 SIRS model, we use an exact gillespie algorithm, one version of which is implemented in the pomp package via the $\bf{gillespie.sim}$ plug-in. We give the trantition of the SIR model code below:
# define the rate function, stoichiometry matrix, and rate-event dependency matrix
SIRS_stoich <- cbind(infection = c(-1, 1, 0), # infections yield S-1, I+1
recovery = c(0, -1, 1), # recoveries yield I-1, R+1
susceptibility = c(1, 0, -1)) # susceptibility events yield R-1, S+1
SIRS_depmat <- cbind(infection = c(1, 1, 0), # infectivity rate updated by changes to S and I
recovery = c(0, 1, 0), # recovery rate updated by changes to I
susceptibility = c(0, 0, 1)) # susceptibility rate updated by changes to R
# SIRS rate function
# j the number of the elementary event (1 = infection, 2 = recovery, 3 = loss of immunity)
# x named numeric vector with the value of the state of the process at time t
# t time
# params named numeric vector containing the parameters
# returns single numerical value with the rate of the elementary event
SIRS_rates <- function(j, x, t, params, ...) {
switch(j,
exp(params["beta"]) * x["S"] * x["I"], # infection
exp(params["mu"]) * x["I"], # recovery
exp(params["gamma"]) * x["R"] # loss of immunity
)
}
# instatiate the gillespie stepper function
SIRS_sim <- gillespie.sim(rate.fun = SIRS_rates,
v = SIRS_stoich,
d = SIRS_depmat)
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:
# initialize the pomp object
sirs <- pomp(
data = data.frame(time = seq(1, 365, by = 7), cases = data), #"cases" is the dataset, "time" is the observation time
times = "time",
t0 = 1, # initial time point
dmeasure = Csnippet(dmeas), # evaluates the density of the measurement process
rmeasure = Csnippet(rmeas), # simulates from the measurement process
rprocess = SIRS_sim, # simulates from the latent process
statenames = c("S", "I", "R"), #state space variable name
paramnames = c("beta", "mu", "gamma", "rho", "theta1", "theta2"), #parameters name
initializer = function(params, t0, ...) {
ps <-exp(params["theta1"]) / (1 + exp(params["theta1"]) + exp(params["theta2"])) # initial prob of susceptible
pi <- 1 / (1 + exp(params["theta1"]) + exp(params["theta2"])) # initial prob of infection
pr <-exp(params["theta2"]) / (1 + exp(params["theta1"]) + exp(params["theta2"])) # initial prob of recovery
return(setNames(as.numeric(rmultinom(1, 200, prob = c(ps, pi, pr))), c("S", "I", "R")))
},
# not used, but pomp needs this
params = c(
beta = log(0.0001),
mu = log(1 / 14),
gamma = log(1 / 150),
rho = log(0.85),
theta1 = log(49),
theta2 = log(0.00000001)
)
)
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 SIRS model as follows.
sirs.dprior <- function(params, ..., log) {
f <- (
dgamma(exp(params[1]), 0.1, 100, log = TRUE) + params[1] + #log prior for log(infection rate) "beta"
dgamma(exp(params[2]), 1.8, 14, log = TRUE) + params[2] + #log prior for log(recovery rate) "mu"
dgamma(exp(params[3]), 0.0625, 10, log = TRUE) + params[3] + #log prior for log(loss of susceptible) "gamma"
dbeta(exp(params[4]) / (1 + exp(params[4])), 5, 1, log = TRUE) +
params[4] - log((1 + exp(params[4])) ^ 2) + #log prior for logit(sampling probablity) "rho"
log(ddirichlet(c(exp(params[5]) / (1 + exp(params[5]) + exp(params[6])),
1 / (1 + exp(params[5]) + exp(params[6])),
exp(params[6]) / (1 + exp(params[5]) + exp(params[6]))),
c(90, 1.5, 0.01))) +
params[5] + params[6] - 3 * log(1 + exp(params[5]) + exp(params[6])) #log prior for logit(initial value)
)
if (log) {
f
} else {
exp(f)
}
}
And starting value of the PMMH is specified below:
param.initial <- c(
beta = log(abs(rnorm(1, 0.00077, 1e-4))),
mu = log(abs(rnorm(1, 1/18, 1e-4))),
gamma = log(abs(rnorm(1, 1/140, 1e-4))),
rho = 1.5,
theta1 = 3,
theta2 = -5
)
And the following runs 1 PMMH chain for 5000 tunning steps (using 500 particles) to get a suitable multivariate normal random walk proposal diagonal variance-covariance matrix value.
pmcmc1 <- pmcmc(
pomp(sirs, dprior = sirs.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(
beta = 0.1,
#sampling variance for beta
mu = 0.1,
#sampling variance for mu
gamma = 0.1,
#sampling variance for gamma
rho = 0.1,
#sampling variance for rho
theta1 = 0.1,
#sampling variance for theta1
theta2 = 0.1 #sampling variance for theta2
),
scale.start = 100,
shape.start = 100,
max.scaling = 1.1)
)
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(); #calculation of time
run_time <- difftime(end_time, start_time, units = "hours") #calculation of time
pomp_results <- list(time = run_time, results = pmcmc1)
fintzij/BDAepimodel documentation built on Sept. 20, 2020, 1:44 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(BDAepimodel) library(coda) library(Rcpp) library(pomp) library(gtools) library(ggplot2)
Analyzing an epidemic with SIRS dynamics
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 simulated prevalence counts. We fit SIRS models to binomially distributed prevalence counts.
We simulated an epidemic with SIRS dynamics in a population of 200 individuals, 1% of whom were initially infected, with the rest being suscpetible. The per–contact infectivity rate was $\beta = 0.0009$, the average infectious period duration was $1/\mu = 14$ days, and the average time until loss of immunity was $1/\gamma = 150$ days, which together correspond to a basic reproduction number of $R_0 = \beta N / \mu = 2.52$. Binomially sampled prevalence was observed at one week intervals over a period of one year with sampling probability $\rho=0.8$.
The data simulation code is in "SIR_SEIR_SIRS" vignettes. And the data is shown below.
data<-c(1,4,6,6,4,9,12,23,29,32,29,48,32,22,15,13,8,7,3,2,1,0,0,1,3,1,0,1,3,3,3,1,5,3,5,3,2,1,0,6,2,3,6,3,6,5,5,8,15,17,14,8,8)
Using Approximate Tau-leaping Algorithm:
We use the SIRS model to model the dataset. 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 lose immunity after a few days of the recovery. 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=rbinom(I,exp(rho)/(1+exp(rho)));//represent the data " dmeas<-" lik=dbinom(cases,I,exp(rho)/(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 and density function dbinom 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 SIRS model, we use a approximating 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 SIRS model code below:
sirs.step<-" double rate[3]; double dN[3]; rate[0]=exp(beta)*I; //Infection rate rate[1]=exp(mu); //Recovery rate rate[2]=exp(gamma); //Loss of immunity 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 infected people reulermultinom(1,R,&rate[2],dt,&dN[2]); //generate the number of people who newly lose immunity 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]+dN[2]; //update the number of Susceptible I+=dN[0]-dN[1]; //update the number of Infection R+=dN[1]-dN[2]; //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:
sirs <- pomp( data = data.frame(cases = data, time = seq(1, 365, by = 7)), #"cases" is the dataset, "time" is the observation time times = "time", t0 = 1, #initial time point dmeasure = Csnippet(dmeas), rmeasure = Csnippet(rmeas), #return the likelihood rprocess = euler.sim(step.fun = Csnippet(sirs.step), delta.t = 1/3), #delta.t is here statenames = c("S", "I", "R"), #state space variable name paramnames = c("beta", "mu", "gamma", "rho", "theta1", "theta2"), #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, 200, prob = c(ps, pi, pr))), c("S", "I", "R"))) }, # not used, but pomp needs this params = c( beta = log(0.0001), mu = log(1 / 14), gamma = log(1 / 140), rho = 1.5, theta1 = log(49), theta2 = log(0.00000001) ) )
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 SIRS model as follows.
sirs.dprior <- function(params, ..., log) { f <- ( dgamma(exp(params[1]), 0.1, 100, log = TRUE) + params[1] + #log prior for log(infection rate) "beta" dgamma(exp(params[2]), 1.8, 14, log = TRUE) + params[2] + #log prior for log(recovery rate) "mu" dgamma(exp(params[3]), 0.0625, 10, log = TRUE) + params[3] + #log prior for log(loss of susceptible) "gamma" dbeta(exp(params[4]) / (1 + exp(params[4])), 5, 1, log = TRUE) + params[4] - log((1 + exp(params[4])) ^ 2) + #log prior for logit(sampling probablity) "rho" log(ddirichlet(c(exp(params[5]) / (1 + exp(params[5]) + exp(params[6])), 1 / (1 + exp(params[5]) + exp(params[6])), exp(params[6]) / (1 + exp(params[5]) + exp(params[6]))), c(9, 0.15, 0.001))) + params[5] + params[6] - 3 * log(1 + exp(params[5]) + exp(params[6])) #log prior for logit(initial value) ) if (log) { f } else { exp(f) } }
And starting value of the PMMH is specified below:
param.initial <- c( beta = log(abs(rnorm(1, 0.00077, 1e-4))), mu = log(abs(rnorm(1, 1/18, 1e-4))), gamma = log(abs(rnorm(1, 1/140, 1e-4))), rho = 1.5, theta1 = 3, theta2 = -5 )
And the following runs 1 PMMH chain for 5000 tunning steps (using 500 particles) to get a suitable multivariate normal random walk proposal diagonal variance-covariance matrix value.
pmcmc1 <- pmcmc( pomp(sirs, dprior = sirs.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( beta = 0.1, #sampling variance for beta mu = 0.1, #sampling variance for mu gamma = 0.1, #sampling variance for gamma rho = 0.1, #sampling variance for rho theta1 = 0.1,#sampling variance for theta1 theta2 = 0.1 #sampling variance for theta2 ), scale.start = 100, shape.start = 100) )
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(); #calculation of time run_time <- difftime(end_time, start_time, units = "hours") #calculation of time pomp_results <- list(time = run_time, results = pmcmc1)
Using Exact Gillespie Algorithm:
We use the SIRS model to model the dataset. 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 lose immunity after a few days of the recovery. 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=rbinom(I,exp(rho)/(1+exp(rho)));//represent the data " dmeas<-" lik=dbinom(cases,I,exp(rho)/(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 and density function dbinom 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 SIRS model, we use an exact gillespie algorithm, one version of which is implemented in the pomp package via the $\bf{gillespie.sim}$ plug-in. We give the trantition of the SIR model code below:
# define the rate function, stoichiometry matrix, and rate-event dependency matrix SIRS_stoich <- cbind(infection = c(-1, 1, 0), # infections yield S-1, I+1 recovery = c(0, -1, 1), # recoveries yield I-1, R+1 susceptibility = c(1, 0, -1)) # susceptibility events yield R-1, S+1 SIRS_depmat <- cbind(infection = c(1, 1, 0), # infectivity rate updated by changes to S and I recovery = c(0, 1, 0), # recovery rate updated by changes to I susceptibility = c(0, 0, 1)) # susceptibility rate updated by changes to R # SIRS rate function # j the number of the elementary event (1 = infection, 2 = recovery, 3 = loss of immunity) # x named numeric vector with the value of the state of the process at time t # t time # params named numeric vector containing the parameters # returns single numerical value with the rate of the elementary event SIRS_rates <- function(j, x, t, params, ...) { switch(j, exp(params["beta"]) * x["S"] * x["I"], # infection exp(params["mu"]) * x["I"], # recovery exp(params["gamma"]) * x["R"] # loss of immunity ) } # instatiate the gillespie stepper function SIRS_sim <- gillespie.sim(rate.fun = SIRS_rates, v = SIRS_stoich, d = SIRS_depmat)
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:
# initialize the pomp object sirs <- pomp( data = data.frame(time = seq(1, 365, by = 7), cases = data), #"cases" is the dataset, "time" is the observation time times = "time", t0 = 1, # initial time point dmeasure = Csnippet(dmeas), # evaluates the density of the measurement process rmeasure = Csnippet(rmeas), # simulates from the measurement process rprocess = SIRS_sim, # simulates from the latent process statenames = c("S", "I", "R"), #state space variable name paramnames = c("beta", "mu", "gamma", "rho", "theta1", "theta2"), #parameters name initializer = function(params, t0, ...) { ps <-exp(params["theta1"]) / (1 + exp(params["theta1"]) + exp(params["theta2"])) # initial prob of susceptible pi <- 1 / (1 + exp(params["theta1"]) + exp(params["theta2"])) # initial prob of infection pr <-exp(params["theta2"]) / (1 + exp(params["theta1"]) + exp(params["theta2"])) # initial prob of recovery return(setNames(as.numeric(rmultinom(1, 200, prob = c(ps, pi, pr))), c("S", "I", "R"))) }, # not used, but pomp needs this params = c( beta = log(0.0001), mu = log(1 / 14), gamma = log(1 / 150), rho = log(0.85), theta1 = log(49), theta2 = log(0.00000001) ) )
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 SIRS model as follows.
sirs.dprior <- function(params, ..., log) { f <- ( dgamma(exp(params[1]), 0.1, 100, log = TRUE) + params[1] + #log prior for log(infection rate) "beta" dgamma(exp(params[2]), 1.8, 14, log = TRUE) + params[2] + #log prior for log(recovery rate) "mu" dgamma(exp(params[3]), 0.0625, 10, log = TRUE) + params[3] + #log prior for log(loss of susceptible) "gamma" dbeta(exp(params[4]) / (1 + exp(params[4])), 5, 1, log = TRUE) + params[4] - log((1 + exp(params[4])) ^ 2) + #log prior for logit(sampling probablity) "rho" log(ddirichlet(c(exp(params[5]) / (1 + exp(params[5]) + exp(params[6])), 1 / (1 + exp(params[5]) + exp(params[6])), exp(params[6]) / (1 + exp(params[5]) + exp(params[6]))), c(90, 1.5, 0.01))) + params[5] + params[6] - 3 * log(1 + exp(params[5]) + exp(params[6])) #log prior for logit(initial value) ) if (log) { f } else { exp(f) } }
And starting value of the PMMH is specified below:
param.initial <- c( beta = log(abs(rnorm(1, 0.00077, 1e-4))), mu = log(abs(rnorm(1, 1/18, 1e-4))), gamma = log(abs(rnorm(1, 1/140, 1e-4))), rho = 1.5, theta1 = 3, theta2 = -5 )
And the following runs 1 PMMH chain for 5000 tunning steps (using 500 particles) to get a suitable multivariate normal random walk proposal diagonal variance-covariance matrix value.
pmcmc1 <- pmcmc( pomp(sirs, dprior = sirs.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( beta = 0.1, #sampling variance for beta mu = 0.1, #sampling variance for mu gamma = 0.1, #sampling variance for gamma rho = 0.1, #sampling variance for rho theta1 = 0.1, #sampling variance for theta1 theta2 = 0.1 #sampling variance for theta2 ), scale.start = 100, shape.start = 100, max.scaling = 1.1) )
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(); #calculation of time run_time <- difftime(end_time, start_time, units = "hours") #calculation of time pomp_results <- list(time = run_time, results = pmcmc1)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.