Title: | Spatial and Network Based Individual Level Models for Epidemics |
---|---|
Description: | Provides tools for simulating from discrete-time individual level models for infectious disease data analysis. This epidemic model class contains spatial and contact-network based models with two disease types: Susceptible-Infectious (SI) and Susceptible-Infectious-Removed (SIR). |
Authors: | Vineetha Warriyar. K. V., Waleed Almutiry, and Rob Deardon |
Maintainer: | Waleed Almutiry <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.5.2 |
Built: | 2024-11-04 04:34:41 UTC |
Source: | https://github.com/waleedalmutiry/epiilm |
The R package EpiILM is provided for simulating from, and carrying out Bayesian MCMC-based statistical inference for spatial and/or network-based individual-level modelling framework. The package allows for the incorporation of individual-level susceptibility and transmissibility covariates in models, and provides various methods of summarizing epidemic data sets.
The R package EpiILM can be used to carry out simulation of epidemics, estimate the basic reproduction number, plot various epidemic summary graphics, calculate the log-likelihood, carry out Bayesian inference using Metropolis-Hastings MCMC, and implement posterior predictive checks and model selection for a given data set and model. The key functions for this package are detailed in the value section. One of the important functions epimcmc
depends heavily on the MCMC
from the adaptMCMC package for performing the MCMC analysis. This function implements the robust adaptive Metropolis sampler of Vihola (2012) for tuning the covariance matrix of the (normal) jump distribution adaptively to achieve the desired acceptance rate. The package has other features for making predictions or forecasting for a specific model via the pred.epi
function. The main functions, including for epidemic simulation (epidata
) and likelihood calculation (epilike
) are coded in Fortran in order to achieve the goal of agile implementation.
Key functions for this package:
\link{epidata} |
Simulates epidemics for the specified model type and parameters. |
\link{epilike} |
Calculates the log-likelihood for the specified model and data set. |
\link{epimcmc} |
Runs an MCMC algorithm for the estimation of specified model parameters. |
\link{pred.epi} |
Computes posterior predictions for a specified epidemic model. |
Vineetha Warriyar. K. V., Waleed Almutiry, and Rob Deardon
Maintainer: Waleed Almutiry <[email protected]>
Deardon, R., Brooks, S. P., Grenfell, B. T., Keeling, M. J., Tildesley, M. J., Savill, N. J., Shaw, D. J., and Woolhouse, M. E. (2010). Inference for individual level models of infectious diseases in large populations. Statistica Sinica, 20, 239-261.
Vihola, M. (2012) Robust adaptive Metropolis algorithm with coerced acceptance rate. Statistics and Computing, 22(5), 997-1008. doi:10.1007/s11222-011-9269-5.
## Not run: demo(EpiILM.spatial) demo(EpiILM.network) ## End(Not run)
## Not run: demo(EpiILM.spatial) demo(EpiILM.network) ## End(Not run)
This function allows the user to generate objects of class "epidata"
. The output of this function provides information required to be used in the other functions in the package.
as.epidata (type, n, x = NULL, y = NULL, inftime, infperiod = NULL, contact = NULL)
as.epidata (type, n, x = NULL, y = NULL, inftime, infperiod = NULL, contact = NULL)
type |
Type of compartment framework, with the choice of "SI" for Susceptible-Infectious diseases and "SIR" for Susceptible-Infectious-Removed. |
n |
Population size |
x |
X coordinates of individuals. |
y |
Y coordinates of individuals. |
inftime |
Times at which individuals are infected to initialize epidemic simulation. |
infperiod |
Length of infectious period for each individual. |
contact |
A contact network matrix or an array of contact network matrices. |
An object of class epidata
is returned containing the following:
Type of compartment framework, with the choice of "SI" for Susceptible-Infectious diseases and "SIR" for Susceptible-Infectious-Removed.
The XY-coordinates of individuals if distance-based ILM is used, otherwise is NULL.
Contact network matrix if network-based ILM is used, otherwise is NULL.
The infection times of individuals.
The removal times of individuals when type
= “SIR”.
# generate 100 X-Y coordinates for a distance-based ILM x <- runif(100, 0, 10) y <- runif(100, 0, 10) # suppose we know the infection times for a spatial SI model based on above X-Y coordinates inftime <- rpois(100, 8) # Now we can convert above information to an epidata object with out <- as.epidata(type = "SI", n = 100, x = x, y = y, inftime = inftime ) out
# generate 100 X-Y coordinates for a distance-based ILM x <- runif(100, 0, 10) y <- runif(100, 0, 10) # suppose we know the infection times for a spatial SI model based on above X-Y coordinates inftime <- rpois(100, 8) # Now we can convert above information to an epidata object with out <- as.epidata(type = "SI", n = 100, x = x, y = y, inftime = inftime ) out
Gives a Monte Carlo estimate of the basic reproduction number for a specified SIR model and data set
epiBR0 (x = NULL, y = NULL, contact = NULL, sus.par, trans.par = NULL, beta, spark = NULL, infperiod, Sformula = NULL, Tformula = NULL, tmax, niter)
epiBR0 (x = NULL, y = NULL, contact = NULL, sus.par, trans.par = NULL, beta, spark = NULL, infperiod, Sformula = NULL, Tformula = NULL, tmax, niter)
x |
X coordinates of individuals |
y |
Y coordinates of individuals |
contact |
Contact network(s) |
sus.par |
Susceptibility parameter(>0) |
trans.par |
Transmissibility parameter(>0) |
beta |
Spatial parameter(s) (>0) or network parameter (s) (>0) if contact is used |
spark |
Sparks parameter (>=0), representing infections unexplained by other parts of the model or infections coming in from outside the observed population, default value is zero |
infperiod |
Length of infectious period for each individual |
Sformula |
An object of class formula. See formula Individual-level covariate information associated with susceptibility can be passed through this argument. An expression of the form |
Tformula |
An object of class formula. See formula Individual-level covariate information associated with transmissibility can be passed through this argument. An expression of the form |
tmax |
The last time point of simulation |
niter |
Number of epidemic simulations to calculate basic reproduction number |
A list is returned with the following components:
BasicR0 |
The basic reproduction number value |
simulated_BR0 |
Number of infections per simulation |
# generate 100 X-Y coordinates for a distance-based ILM x <- runif(100, 0, 10) y <- runif(100, 0, 10) # Suppose we know the length of infectious period for each individual. Also, assume # susceptibility parameter = 1.5 and spatial parameter = 5 for this SIR model infperiod <- rep(3, 100) # For a 1000 iteration with a last observed time point 15, we can estimate the basic # reproduction number using Monte Carlo simulation out <- epiBR0(x = x, y = y, sus.par = 1.5, beta = 5, infperiod= infperiod, tmax = 15, niter = 1000) out$BasicR0
# generate 100 X-Y coordinates for a distance-based ILM x <- runif(100, 0, 10) y <- runif(100, 0, 10) # Suppose we know the length of infectious period for each individual. Also, assume # susceptibility parameter = 1.5 and spatial parameter = 5 for this SIR model infperiod <- rep(3, 100) # For a 1000 iteration with a last observed time point 15, we can estimate the basic # reproduction number using Monte Carlo simulation out <- epiBR0(x = x, y = y, sus.par = 1.5, beta = 5, infperiod= infperiod, tmax = 15, niter = 1000) out$BasicR0
This function allows the user to simulate epidemics under different models and scenarios
epidata (type, n, tmin = NULL, tmax, sus.par, trans.par = NULL, beta = NULL, spark = NULL, Sformula = NULL, Tformula = NULL, x = NULL, y = NULL, inftime = NULL, infperiod = NULL, contact = NULL)
epidata (type, n, tmin = NULL, tmax, sus.par, trans.par = NULL, beta = NULL, spark = NULL, Sformula = NULL, Tformula = NULL, x = NULL, y = NULL, inftime = NULL, infperiod = NULL, contact = NULL)
type |
Type of compartment framework, with the choice of "SI" for Susceptible-Infectious diseases and "SIR" for Susceptible-Infectious-Removed. |
n |
Population size |
tmin |
The time point at which simulation begins, default value is one. |
tmax |
The last time point of simulation. |
sus.par |
Susceptibility parameter (>0). |
trans.par |
Transmissibility parameter (>0). |
beta |
Spatial parameter(s) (>0) or network parameter (s) (>0) if contact network is used. |
spark |
Sparks parameter (>=0), representing infections unexplained by other parts of the model (eg. infections coming in from outside the observed population), default value is zero. |
Sformula |
An object of class formula. See formula. Individual-level covariate information associated with susceptibility can be passed through this argument. An expression of the form |
Tformula |
An object of class formula. See formula. Individual-level covariate information associated with transmissibility can be passed through this argument. An expression of the form |
x |
X coordinates of individuals. |
y |
Y coordinates of individuals. |
inftime |
Times at which individuals are infected to initialize epidemic simulation. |
infperiod |
Length of infectious period for each individual. |
contact |
A contact network matrix or an array of contact network matrices. |
We consider following two individual level models:
Spatial model:
Network model:
where is the probability that susceptible individual i is infected at time point t, becoming infectious at time t+1;
is a susceptibility function which accommodates potential risk factors associated with susceptible individual i contracting the disease;
is a transmissibility function which accommodates potential risk factors associated with infectious individual j;
is a sparks term which represents infections originating from outside the population being observed or some other unobserved infection mechanism.
The susceptibility function can incorporate any individual-level covariates of interest and is treated as a linear function of the covariates, i.e.,
, where
denote
covariates associated with susceptible individual $i$, along with susceptibility parameters
. If the model does not contain any susceptibility covariates then
is used. In a similar way, the transmissibility function can incorporate any individual-level covariates of interest associated with infectious individual.
is also treated as a linear function of the covariates, but without the intercept term, i.e.,
, where
denote the
covariates associated with infectious individual j, along with transmissibility parameters
. If the model does not contain any transmissibility covariates then
is used.
An object of class epidata
is returned containing the following:
Type of compartment framework, with the choice of "SI" for Susceptible-Infectious diseases and "SIR" for Susceptible-Infectious-Removed
The XY-coordinates of individuals.
Contact network matrix.
The infection times of individuals.
The removal times of individuals when type
= “SIR”.
Deardon, R., Brooks, S. P., Grenfell, B. T., Keeling, M. J., Tildesley, M. J., Savill, N. J., Shaw, D. J., and Woolhouse, M. E. (2010). Inference for individual level models of infectious diseases in large populations. Statistica Sinica, 20, 239-261.
Deardon, R., Fang, X., and Kwong, G.P.S. (2014). Statistical modelling of spatio-temporal infectious disease transmission in analyzing and modeling Spatial and temporal dynamics of infectious diseases, (Ed: D. Chen, B. Moulin, J. Wu), John Wiley & Sons. Chapter 11.
plot.epidata
, epimcmc
, epilike
, pred.epi
.
## Example 1: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI",n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) # Plots of epidemic progression (optional) plot(out1, plottype = "spatial") plot(out1, plottype = "curve", curvetype = "newinfect") ## Example 2: spatial SIR model # generate infectious period(=3) for 100 individuals lambda <- rep(3, 100) out2 <- epidata(type = "SIR", n = 100, tmax = 15, sus.par = 0.3, beta = 5.0, infperiod = lambda, x = x, y = y) plot(out2, plottype = "spatial") plot(out2, plottype = "curve", curvetype = "newinfect") ## Example 3: SI network model contact1 <- matrix(rbinom(10000, 1, 0.1), nrow = 100, ncol = 100) contact2 <- matrix(rbinom(10000, 1, 0.1), nrow = 100, ncol = 100) diag(contact1[,] ) <- 0 diag(contact2[,] ) <- 0 contact <- array(c(contact1, contact2), dim = c(100, 100, 2)) out3 <- epidata(type = "SI", n = 100, tmax = 15, sus.par = 0.3, beta = c(3.0, 5.0), contact = contact) plot(out3, plottype = "curve", curvetype = "complete") plot(out3, plottype = "curve", curvetype = "susceptible") plot(out3, plottype = "curve", curvetype = "newinfect") plot(out3, plottype = "curve", curvetype = "totalinfect")
## Example 1: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI",n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) # Plots of epidemic progression (optional) plot(out1, plottype = "spatial") plot(out1, plottype = "curve", curvetype = "newinfect") ## Example 2: spatial SIR model # generate infectious period(=3) for 100 individuals lambda <- rep(3, 100) out2 <- epidata(type = "SIR", n = 100, tmax = 15, sus.par = 0.3, beta = 5.0, infperiod = lambda, x = x, y = y) plot(out2, plottype = "spatial") plot(out2, plottype = "curve", curvetype = "newinfect") ## Example 3: SI network model contact1 <- matrix(rbinom(10000, 1, 0.1), nrow = 100, ncol = 100) contact2 <- matrix(rbinom(10000, 1, 0.1), nrow = 100, ncol = 100) diag(contact1[,] ) <- 0 diag(contact2[,] ) <- 0 contact <- array(c(contact1, contact2), dim = c(100, 100, 2)) out3 <- epidata(type = "SI", n = 100, tmax = 15, sus.par = 0.3, beta = c(3.0, 5.0), contact = contact) plot(out3, plottype = "curve", curvetype = "complete") plot(out3, plottype = "curve", curvetype = "susceptible") plot(out3, plottype = "curve", curvetype = "newinfect") plot(out3, plottype = "curve", curvetype = "totalinfect")
Computes the Deviance Information Criterion for individual level models
epidic (burnin, niter, LLchain, LLpostmean)
epidic (burnin, niter, LLchain, LLpostmean)
burnin |
Burnin period for MCMC |
niter |
Number of MCMC iterations |
LLchain |
Loglikelihood values from the MCMC output |
LLpostmean |
Loglikelihood value of the model with posterior mean of estimates |
Spiegelhalter, D., Best, N., Carlin, B., Van der Linde, A. (2002). Bayesian Measures of Model Complexity and Fit. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 64(4), 583-639.
## Example 1: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) unif_range <- matrix(c(0, 0, 10000, 10000), nrow = 2, ncol = 2) # estimate parameters mcmcout <- epimcmc(out1, tmax = 15, niter = 1500, Sformula = ~covariate, sus.par.ini = c(0.003, 0.01), beta.ini =0.01, pro.sus.var = c(0.1, 0.1),pro.beta.var = 0.5, prior.sus.par = unif_range, prior.sus.dist = c("uniform","uniform"), prior.beta.dist = "uniform", prior.beta.par = c(0, 10000), adapt = TRUE, acc.rate = 0.5 ) # store the estimates sus.parameters = c(mean(unlist(mcmcout$Estimates[1])), mean(unlist(mcmcout$Estimates[2]))) beta.par = mean(unlist(mcmcout$Estimates[3])) # likelihood value loglike <- epilike(out1, tmax = 15, Sformula = ~covariate, sus.par = sus.parameters, beta = beta.par) # deviance information criterion calculation for the above epidemic dic <- epidic(burnin = 500, niter = 1500, LLchain = mcmcout$Loglikelihood, LLpostmean = loglike) dic
## Example 1: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) unif_range <- matrix(c(0, 0, 10000, 10000), nrow = 2, ncol = 2) # estimate parameters mcmcout <- epimcmc(out1, tmax = 15, niter = 1500, Sformula = ~covariate, sus.par.ini = c(0.003, 0.01), beta.ini =0.01, pro.sus.var = c(0.1, 0.1),pro.beta.var = 0.5, prior.sus.par = unif_range, prior.sus.dist = c("uniform","uniform"), prior.beta.dist = "uniform", prior.beta.par = c(0, 10000), adapt = TRUE, acc.rate = 0.5 ) # store the estimates sus.parameters = c(mean(unlist(mcmcout$Estimates[1])), mean(unlist(mcmcout$Estimates[2]))) beta.par = mean(unlist(mcmcout$Estimates[3])) # likelihood value loglike <- epilike(out1, tmax = 15, Sformula = ~covariate, sus.par = sus.parameters, beta = beta.par) # deviance information criterion calculation for the above epidemic dic <- epidic(burnin = 500, niter = 1500, LLchain = mcmcout$Loglikelihood, LLpostmean = loglike) dic
Calculates the log likelihood for the specified individual level model and data set
epilike (object, tmin = NULL, tmax, sus.par, trans.par = NULL, beta = NULL, spark = NULL, Sformula = NULL, Tformula = NULL)
epilike (object, tmin = NULL, tmax, sus.par, trans.par = NULL, beta = NULL, spark = NULL, Sformula = NULL, Tformula = NULL)
object |
An object of class |
tmin |
The first time point at which data is observed, default value is one. |
tmax |
The last time point at which data is observed. |
sus.par |
Susceptibility parameter(>0). |
trans.par |
Transmissibility parameter(>0). |
beta |
Spatial parameter(s) (>0) or network parameter (s) (>0) if contact network is used. |
spark |
Sparks parameter(>=0), representing infections unexplained by other parts of the model or infections coming in from outside the observed population, default value is zero. |
Sformula |
An object of class formula. See formula. Individual-level covariate information associated with susceptibility can be passed through this argument. An expression of the form |
Tformula |
An object of class formula. See formula. Individual-level covariate information associated with transmissibility can be passed through this argument. An expression of the form |
Returns the value of the log-likelihood function.
Deardon R, Brooks, S. P., Grenfell, B. T., Keeling, M. J., Tildesley, M. J., Savill, N. J., Shaw, D. J., Woolhouse, M. E. (2010). Inference for individual level models of infectious diseases in large populations. Statistica Sinica, 20, 239-261.
## Example 1: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) epilike(out1, tmax = 15, sus.par = c(0.1, 0.3), beta = 5, Sformula = ~covariate) ## Example 2: spatial SIR model # generate infectious period (=3) for 100 individuals lambda <- rep(3, 100) out2 <- epidata(type = "SIR", n = 100, tmax = 15, sus.par =0.3, beta = 5.0, infperiod = lambda, x = x, y = y) epilike(out2, tmax = 15, sus.par = 0.3, beta = 5.0)
## Example 1: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) epilike(out1, tmax = 15, sus.par = c(0.1, 0.3), beta = 5, Sformula = ~covariate) ## Example 2: spatial SIR model # generate infectious period (=3) for 100 individuals lambda <- rep(3, 100) out2 <- epidata(type = "SIR", n = 100, tmax = 15, sus.par =0.3, beta = 5.0, infperiod = lambda, x = x, y = y) epilike(out2, tmax = 15, sus.par = 0.3, beta = 5.0)
Runs an MCMC algorithm for the estimation of specified model parameters
epimcmc (object, tmin = NULL, tmax, niter, sus.par.ini, trans.par.ini = NULL, beta.ini = NULL, spark.ini = NULL, Sformula = NULL, Tformula = NULL, pro.sus.var, pro.trans.var = NULL, pro.beta.var = NULL, pro.spark.var = NULL, prior.sus.dist, prior.trans.dist = NULL, prior.beta.dist = NULL, prior.spark.dist = NULL, prior.sus.par, prior.trans.par, prior.beta.par = NULL, prior.spark.par = NULL, adapt = FALSE, acc.rate = NULL)
epimcmc (object, tmin = NULL, tmax, niter, sus.par.ini, trans.par.ini = NULL, beta.ini = NULL, spark.ini = NULL, Sformula = NULL, Tformula = NULL, pro.sus.var, pro.trans.var = NULL, pro.beta.var = NULL, pro.spark.var = NULL, prior.sus.dist, prior.trans.dist = NULL, prior.beta.dist = NULL, prior.spark.dist = NULL, prior.sus.par, prior.trans.par, prior.beta.par = NULL, prior.spark.par = NULL, adapt = FALSE, acc.rate = NULL)
object |
An object of class |
tmin |
The first time point at which the infection occurs, default value is one. |
tmax |
The last time point at which data is observed. |
niter |
Number of MCMC iterations. |
sus.par.ini |
Initial value(s) of the susceptibility parameter(s) (>0). |
trans.par.ini |
Initial value(s) of the transmissibility parameter(s) (>0). |
beta.ini |
Initial value(s) of the spatial parameter(s) (>0) or the network parameter(s) (>0) if contact network is used. |
spark.ini |
Initial value of the spark parameter (>=0). |
Sformula |
An object of class formula. See formula Individual-level covariate information associated with susceptibility can be passed through this argument. An expression of the form |
Tformula |
An object of class formula. See formula Individual-level covariate information associated with transmissibility can be passed through this argument. An expression of the form |
pro.sus.var |
Proposal density variance(s) for susceptibility parameter(s). If a zero value is assigned to the proposal variance of any parameter, the parameter is considered fixed to its |
pro.trans.var |
Proposal density variance(s) for transmissibility parameter(s). If a zero value is assigned to the proposal variance of any parameter, the parameter is considered fixed to its |
pro.beta.var |
Proposal density variance(s) for beta parameter(s). If a zero value is assigned to the proposal variance of any parameter, the parameter is considered fixed to its |
pro.spark.var |
Proposal density variance for the spark parameter. |
prior.sus.dist |
Select the prior distribution(s) for the susceptibility parameter(s) with the choice of "halfnormal" for positive half normal distribution, "gamma" for gamma distribution and "uniform" for uniform distribution |
prior.trans.dist |
Select the prior distribution(s) for the transmissibility parameter(s) with the choice of "halfnormal" for positive half normal distribution, "gamma" for gamma distribution and "uniform" for uniform distribution |
prior.beta.dist |
Select the prior distribution(s) for the beta parameter(s) with the choice of "halfnormal" for half normal distribution, "gamma" for gamma distribution and "uniform" for uniform distribution |
prior.spark.dist |
Select the prior distribution for the spark parameter with the choice of "halfnormal" for half normal distribution, "gamma" for gamma distribution and "uniform" for uniform distribution |
prior.sus.par |
A vector (matrix) of the prior distribution parameters for updating the susceptibility parameter(s). |
prior.trans.par |
A vector (matrix) of the prior distribution parameters for updating the transmissibility parameter(s). |
prior.beta.par |
A vector (matrix) of the prior distribution parameters for updating the kernel parameter(s). |
prior.spark.par |
A vector of the prior distribution parameters for updating the spark parameter. |
adapt |
To enable the adaptive MCMC method in the |
acc.rate |
To set an acceptance rate. This option will be ignored if |
Independent Gaussian random walks are used as the Metropolis-Hastings MCMC proposal for all parameters. The epimcmc
function depends on the MCMC
function from the adaptMCMC package.
Returns an object of class epimcmc
that contains:
the compartmental framework model used in the analysis.
the used kernel.type
in the function (distance-based or network-based).
the MCMC output of the updated model parameters.
the loglikelihood of the updated model parameters.
the MCMC output of all the model parameters (including fixed parameters).
the number of parameters in the susceptibility function.
the number of parameters in the transmissibility function.
the number of parameters in the kernel function.
Rob Deardon, Xuan Fang, and Grace P. S. Kwong (2015). Statistical modelling of spatio-temporal infectious disease tranmission in Analyzing and Modeling Spatial and Temporal Dynamics of Infectious Diseases, (Ed: D. Chen, B. Moulin, J. Wu), John Wiley & Sons.. Chapter 11.
summary.epimcmc
, plot.epimcmc
, epidata
, epilike
, pred.epi
.
## Example 1: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) alphapar1 <- matrix(c(1, 1, 1, 1), ncol = 2, nrow = 2) betapar1 <- c(10, 2) epi <- epimcmc(object = out1, tmin = 1, tmax = 15, niter = 1000, sus.par.ini = c(1, 1), beta.ini = 1, Sformula = ~covariate, pro.sus.var = c(0.5, 0.3), pro.beta.var = 0.1, prior.sus.dist = c("gamma", "gamma"), prior.beta.dist = "gamma", prior.sus.par = alphapar1, prior.beta.par = betapar1, adapt = TRUE, acc.rate = 0.5) epi ## Example 2: spatial SIR model lambda <- rep(3, 100) out2 <- epidata(type = "SIR", n = 100, tmax = 15, sus.par = 0.3, beta = 5.0, infperiod = lambda, x = x, y = y) alphapar2 <- c(1, 1) betapar2 <- c(1, 1) epi2 <- epimcmc(object = out2, tmin = 1, tmax = 15, niter = 1000, sus.par.ini = 1, beta.ini = 1, Sformula = NULL, pro.sus.var = 0.3, pro.beta.var = 0.1, prior.sus.dist = "gamma", prior.beta.dist = "gamma", prior.sus.par = alphapar2, prior.beta.par = betapar2, adapt = FALSE, acc.rate = NULL) epi2
## Example 1: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) alphapar1 <- matrix(c(1, 1, 1, 1), ncol = 2, nrow = 2) betapar1 <- c(10, 2) epi <- epimcmc(object = out1, tmin = 1, tmax = 15, niter = 1000, sus.par.ini = c(1, 1), beta.ini = 1, Sformula = ~covariate, pro.sus.var = c(0.5, 0.3), pro.beta.var = 0.1, prior.sus.dist = c("gamma", "gamma"), prior.beta.dist = "gamma", prior.sus.par = alphapar1, prior.beta.par = betapar1, adapt = TRUE, acc.rate = 0.5) epi ## Example 2: spatial SIR model lambda <- rep(3, 100) out2 <- epidata(type = "SIR", n = 100, tmax = 15, sus.par = 0.3, beta = 5.0, infperiod = lambda, x = x, y = y) alphapar2 <- c(1, 1) betapar2 <- c(1, 1) epi2 <- epimcmc(object = out2, tmin = 1, tmax = 15, niter = 1000, sus.par.ini = 1, beta.ini = 1, Sformula = NULL, pro.sus.var = 0.3, pro.beta.var = 0.1, prior.sus.dist = "gamma", prior.beta.dist = "gamma", prior.sus.par = alphapar2, prior.beta.par = betapar2, adapt = FALSE, acc.rate = NULL) epi2
Produces various graphs summarizing epidemic of class epidata
.
## S3 method for class 'epidata' plot(x, plottype, curvetype = NULL, time_id = NULL, tmin = NULL, timepoints = NULL, ...)
## S3 method for class 'epidata' plot(x, plottype, curvetype = NULL, time_id = NULL, tmin = NULL, timepoints = NULL, ...)
x |
An object of class |
plottype |
Provide two types of plots. When |
curvetype |
It has four options: "complete", "susceptible", "newinfect", and "totalinfect". See details for more information. |
time_id |
Specify time points at which the spatial square is plotted. |
tmin |
Initial time point at which infection occurs, default value is one. |
timepoints |
Specify time points at which the curve is plotted |
... |
........ |
The argument plottype
has two options. When plottype
=“spatial” spatial plots are produced for the epidemic progression over time, and when plottype
=“curve”, the argument curvetype
has to be specified to one of the four available options: "complete" for plotting the number of susceptible, infected and removed individuals at each time point, "susceptible" for plotting the number of susceptible individuals at each time point, "newinfect" for plotting the number of newly infected individuals at each time point, and "totalinfect" for plotting the cumulative number of infected individuals at each time point.
plot
epidata
, plot.epimcmc
, plot.pred.epi
.
## Example : spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI",n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) # Plots of epidemic progression plot(out1, plottype = "spatial") plot(out1, plottype = "curve", curvetype = "newinfect")
## Example : spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI",n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) # Plots of epidemic progression plot(out1, plottype = "spatial") plot(out1, plottype = "curve", curvetype = "newinfect")
epimcmc
objectplot.epimcmc
is an S3 method that plots the output of an S3 object of class epimcmc
.
## S3 method for class 'epimcmc' plot(x, partype, start = 1, end = NULL, thin = 1, ...)
## S3 method for class 'epimcmc' plot(x, partype, start = 1, end = NULL, thin = 1, ...)
x |
An S3 object of class |
partype |
Determines which of two options to plot the output of the |
start , end , thin
|
options for creating |
... |
additional arguments that are passed to the generic |
plot.
epimcmc
, summary.epimcmc
, mcmc
, plot.mcmc
.
## Example : spatial SI model # generate 100 individuals set.seed(59991) x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) alphapar1 <- matrix(c(1, 1, 1, 1), ncol = 2, nrow = 2) betapar1 <- c(10, 2) epi <- epimcmc(object = out1, tmin = 1, tmax = 15, niter = 1000, sus.par.ini = c(0.1, 0.1), beta.ini = 5, Sformula = ~covariate, pro.sus.var = c(0.2, 0.3), pro.beta.var = 0.8, prior.sus.dist = c("gamma", "gamma"), prior.beta.dist = "gamma", prior.sus.par = alphapar1, prior.beta.par = betapar1, adapt = TRUE, acc.rate = 0.5) # plot estimates plot(epi, partype = "parameter", start = 100)
## Example : spatial SI model # generate 100 individuals set.seed(59991) x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) alphapar1 <- matrix(c(1, 1, 1, 1), ncol = 2, nrow = 2) betapar1 <- c(10, 2) epi <- epimcmc(object = out1, tmin = 1, tmax = 15, niter = 1000, sus.par.ini = c(0.1, 0.1), beta.ini = 5, Sformula = ~covariate, pro.sus.var = c(0.2, 0.3), pro.beta.var = 0.8, prior.sus.dist = c("gamma", "gamma"), prior.beta.dist = "gamma", prior.sus.par = alphapar1, prior.beta.par = betapar1, adapt = TRUE, acc.rate = 0.5) # plot estimates plot(epi, partype = "parameter", start = 100)
Produces various graphs for the output of the posterior predictive check of class pred.epi
.
## S3 method for class 'pred.epi' plot(x, ...)
## S3 method for class 'pred.epi' plot(x, ...)
x |
An object of class |
... |
additional arguments to be passed to |
plot
pred.epi
, plot.epidata
, plot.epimcmc
.
Computing the posterior predictive check based on different summary statistics.
pred.epi (object, xx, criterion , n.samples, burnin = NULL, tmin = NULL, Sformula = NULL, Tformula = NULL, showProgressBar = interactive())
pred.epi (object, xx, criterion , n.samples, burnin = NULL, tmin = NULL, Sformula = NULL, Tformula = NULL, showProgressBar = interactive())
object |
An object of class |
xx |
An object of class |
criterion |
The (multivariate) statistical criteria used in the posterior predictive check. It has three options: “newly infectious” which is a multivariate statistics represented by the number of newly infectious individuals over time, “epidemic length” represents the length of epidemic, and “peak time” represents the time of the peak of epidemic. |
n.samples |
The number of epidemics that needs to be simulated in the posterior predictive check procedure. |
burnin |
A scalar value which represents the number of samples needs to be discarded from the MCMC output. |
tmin |
The first time point at which the infection occurs, default value is one. |
Sformula |
An object of class formula. See formula. Individual-level covariate information associated with susceptibility can be passed through this argument. An expression of the form |
Tformula |
An object of class formula. See formula. Individual-level covariate information associated with transmissibility can be passed through this argument. An expression of the form |
showProgressBar |
logical. If TRUE a progress bar is shown. |
An object of class pred.epi
that contains the following:
The compartmental framework model used in the analysis.
The (multivariate) statistical criteria used in the posterior predictive check.
The output of the evaluated criterion
on the simulated epidemics.
The output of the evaluated criterion
on the observed epidemics.
The last time point at which data is observed.
The number of simulated epidemics used in the posterior predictive check procedure.
Deardon R, Brooks, S. P., Grenfell, B. T., Keeling, M. J., Tildesley, M. J., Savill, N. J., Shaw, D. J., Woolhouse, M. E. (2010). Inference for individual level models of infectious diseases in large populations. Statistica Sinica, 20, 239-261.
epimcmc
, epidata
, epilike
, plot.pred.epi
.
## Example 1: spatial SI model # generate 100 individuals set.seed(59991) x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- cbind(runif(100, 0, 2), rbinom(100, 1, 0.5)) out <- epidata(type = "SI",n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3, 0.01), beta = 5.0, x = x, y = y) alphapar2 <- matrix(c(1, 1, 1, 1, 1, 1), ncol = 2, nrow = 3) betapar2 <- c(1, 1) epi<-epimcmc(object = out, tmin = 1, tmax = 15, niter = 500, sus.par.ini = c(1, 1, 1), beta.ini = 1, Sformula = ~covariate, pro.sus.var = c(0.5, 0.3, 0.2), pro.beta.var = 0.1, prior.sus.dist = c("gamma", "gamma", "gamma"), prior.beta.dist = "gamma", prior.sus.par = alphapar2, prior.beta.par = betapar2, adapt = TRUE, acc.rate = 0.5) epipred1 <- pred.epi (object = out, xx = epi, criterion = "newly infectious", n.samples = 100, burnin = 200, tmin = 1, Sformula = ~covariate) plot(epipred1, col = "red", type = "b", lwd = 2) epipred2 <- pred.epi (object = out, xx = epi, criterion = "peak time", n.samples = 100, burnin = 200, tmin = 1, Sformula = ~covariate) plot(epipred2, col = "dark gray")
## Example 1: spatial SI model # generate 100 individuals set.seed(59991) x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- cbind(runif(100, 0, 2), rbinom(100, 1, 0.5)) out <- epidata(type = "SI",n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3, 0.01), beta = 5.0, x = x, y = y) alphapar2 <- matrix(c(1, 1, 1, 1, 1, 1), ncol = 2, nrow = 3) betapar2 <- c(1, 1) epi<-epimcmc(object = out, tmin = 1, tmax = 15, niter = 500, sus.par.ini = c(1, 1, 1), beta.ini = 1, Sformula = ~covariate, pro.sus.var = c(0.5, 0.3, 0.2), pro.beta.var = 0.1, prior.sus.dist = c("gamma", "gamma", "gamma"), prior.beta.dist = "gamma", prior.sus.par = alphapar2, prior.beta.par = betapar2, adapt = TRUE, acc.rate = 0.5) epipred1 <- pred.epi (object = out, xx = epi, criterion = "newly infectious", n.samples = 100, burnin = 200, tmin = 1, Sformula = ~covariate) plot(epipred1, col = "red", type = "b", lwd = 2) epipred2 <- pred.epi (object = out, xx = epi, criterion = "peak time", n.samples = 100, burnin = 200, tmin = 1, Sformula = ~covariate) plot(epipred2, col = "dark gray")
epimcmc
objectsSummarize a epimcmc
object and return an object of class summary.epimcmc
.
## S3 method for class 'epimcmc' summary(object, ...)
## S3 method for class 'epimcmc' summary(object, ...)
object |
an S3 object of class |
... |
potential further arguments (require by generic). |
## Example: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) alphapar1 <- matrix(c(1, 1, 1, 1), ncol = 2, nrow = 2) betapar1 <- c(1, 1) epi <- epimcmc(object = out1, tmin = 1, tmax = 15, niter = 1000, sus.par.ini = c(1, 1), beta.ini = 1, Sformula = ~covariate, pro.sus.var = c(0.5, 0.3), pro.beta.var = 0.1, prior.sus.dist = c("gamma", "gamma"), prior.beta.dist = "gamma", prior.sus.par = alphapar1, prior.beta.par = betapar1, adapt = TRUE, acc.rate = 0.5) # summary of mcmc output summary(epi)
## Example: spatial SI model # generate 100 individuals x <- runif(100, 0, 10) y <- runif(100, 0, 10) covariate <- runif(100, 0, 2) out1 <- epidata(type = "SI", n = 100, Sformula = ~covariate, tmax = 15, sus.par = c(0.1, 0.3), beta = 5.0, x = x, y = y) alphapar1 <- matrix(c(1, 1, 1, 1), ncol = 2, nrow = 2) betapar1 <- c(1, 1) epi <- epimcmc(object = out1, tmin = 1, tmax = 15, niter = 1000, sus.par.ini = c(1, 1), beta.ini = 1, Sformula = ~covariate, pro.sus.var = c(0.5, 0.3), pro.beta.var = 0.1, prior.sus.dist = c("gamma", "gamma"), prior.beta.dist = "gamma", prior.sus.par = alphapar1, prior.beta.par = betapar1, adapt = TRUE, acc.rate = 0.5) # summary of mcmc output summary(epi)
Data extracted from Hughes et al. (1997). Data obtained from a field experiment as the spatial dynamics of tomato spotted wilt virus (tswv).
data(tswv)
data(tswv)
A data frame with following variables
X coordinate
Y coordinate
Infection times
Times at which individuals are removed
Hughes, G., McRoberts,N., Madden, L.V., Nelson, S. C. (1997). Validating mathematical models of plant disease progress in space and time. IMA Journal of Mathematics Applied in Medicine and Biology, 14, 85-112.
data("tswv") x <- tswv$x y <- tswv$y inftime <- tswv$inftime removaltime <- tswv$removaltime infperiod <- rep(3, length(x)) # change to epilate object epidat.tswv <- as.epidata(type = "SIR", n = 520, x = x, y = y, inftime = inftime, infperiod = infperiod) # plot plot(epidat.tswv, plottype = "spatial", tmin = 2)
data("tswv") x <- tswv$x y <- tswv$y inftime <- tswv$inftime removaltime <- tswv$removaltime infperiod <- rep(3, length(x)) # change to epilate object epidat.tswv <- as.epidata(type = "SIR", n = 520, x = x, y = y, inftime = inftime, infperiod = infperiod) # plot plot(epidat.tswv, plottype = "spatial", tmin = 2)