The R package EpiILMCT provides tools for simulating from continuous-time individual level models of disease transmission, and carrying out infectious disease data analyses with the same models. The epidemic models considered are distance-based and contact network-based models within Susceptible-Infectious-Removed (SIR) or Susceptible-Infectious-Notified-Removed (SINR) compartmental frameworks.
You can install the EpiILMCT version from CRAN.
install.packages('EpiILMCT', dependencies = TRUE)You can install the development version from Github
# install.packages("devtools")
devtools::install_github("waleedalmutiry/EpiILMCT-package")Different types of undirected unweighted contact networks can be generated through the function contactnet. This function function has a 'type' option to specify the type of network to be used. It has three available options ("powerlaw", "Cauchy", and "random") for the network model, where the first two options simulate spatial contact networks in which the probability of connections between individuals are based on required XY coordinate input. The output of this function is stored as an object of class 'contactnet'. Also, an S3 method plot function is introduced which uses this object as its input to provide a fancy plot of the contact network.
library(EpiILMCT)
set.seed(22)
# to generate the XY coordinates of 50 individuals:
loc<- matrix(cbind(runif(50, 0, 10),runif(50, 0, 10)), ncol = 2, nrow = 50)net1<- contactnet(type = "powerlaw", location = loc, beta = 1.5,
nu = 0.5)
plot(net1)net2<- contactnet(type = "Cauchy", location = loc, beta = 0.5)
plot(net2)net3<- contactnet(type = "random", num.id = 50, beta = 0.08)
plot(net3) # the default options in igraph package.
plot(net3, vertex.color = "red", vertex.size = 15, edge.color = "black",
vertex.label.cex = 1, vertex.label.color = "black")
The function datagen allows the user to generate epidemics from the continuous time ILMs under the SIR or SINR compartmental frameworks. This function can generate epidemics based on different kernels through the kerneltype argument which takes one of three options: "distance" for distance-based, "network" for network-based, or "both" for distance and network-based. The appropriate kernel matrix must also be provided via the kernelmatrix argument. If "distance" is chosen as the kerneltype, the user must choose a spatial kernel ("powerlaw" or "Cauchy") through the distancekernel argument. Here is an example of generating an epidemic from the SIR continuous time network-based ILM.
library("EpiILMCT")
set.seed(91938)
# To simulate the XY coordinate of 50 individuals and their corresponding binary covariate values:
loc <- matrix(cbind(runif(50, 0, 10), runif(50, 0, 10)), ncol = 2, nrow = 50)
cov <- cbind(rep(1, 50), rbinom(50, 1, 0.5))
# To simulate the contact network:
net <- contactnet(type = "powerlaw", location = loc, beta = 1.8, nu = 1)
# To simulate the epidemic:
epi <- datagen(type = "SIR", kerneltype = "network", kernelmatrix = net, suspar = c(0.08, 0.5), delta = c(4, 2),
suscov = cov)The output of the datagen function is stored as an datagen object which takes a list of:
- type
- kerneltype
- epidat (event times)
- location (XY coordinates of individuals)
- network (contact network matrix), in the case of setting the kerneltype to distance, a zero contact network matrix will be created for the network option.
The package also contains an S3 method plot.datagen function, which illustrates disease spread through the epidemic timeline. This function can be used for either distance-based or network-based ILMs. The object of this function has to be of class datagen. The plot S3 function has a plottype argument that can be set to "history", to produce epidemic curves of infection and removal times, or set to "propagation" to produce plots of the epidemic propagation over time. The following commands are to produce the below two graphs for the generated epidemic above.
plot(epi, plottype = "propagation", time.index = seq_len(6))
plot(epi, plottype = "history")
Metropolis-Hastings MCMC is performed to estimate the joint posterior of the model parameters and latent variables (the latter if various event times are assumed unknown). This is achieved using the function epictmcmc. The below figure illustrates the structure of most arguments of this function.
The output of this function is an object of class epictmcmc. There are S3 methods: print.epictmcmc, summary.epictmcmc and plot.epictmcmc that depend on the coda package. The latter function produced the trace plots of the posterior distributions of the model parameters with the same options of the plot.mcmc function in the coda package, such as, start, thin, and density. An argument plottype needs to be specified to parameter in order to get the above trace plots. Also, as the parameter.samples in the object epictmcmc is stored as an mcmc object of the coda package, the plot.mcmc function in the coda package can be used directly using the input as the parameter.samples.
In the case of datatype is set to either "known removal" (unknown infection times only) or "unknown removal" (unknown infection and removal times), plots of the average posterior and 95% CI of the unobserved event times can also be produced through specifying the argument plottype to inf.times or rem.times.
The class epictmcmc contains the MCMC samples of the model parameters and the missing information (in case datatype is not set to known epidemic), and other useful information to be used in other functions such as the above S3 methods. For example, when datatype = known epidemic, the class epictmcmc has a list contained the following:
- "compart.framework"
- "kernel.type"
- "data.assumption"
- "parameter.samples"
- "log.likelihood"
- "acceptance.rate"
- "number.iteration"
- "number.parameter"
- "number.chains"
The following commands are to perform the MCMC for analyzing the above epidemic data set using the epictmcmc function.
suscov <- list(NULL)
suscov[[1]] <- list(c(0.01, 0.1), matrix(c("gamma", "gamma", 1, 1, 0.1, 0.1, 0.5, 1), ncol = 4, nrow = 2))
suscov[[2]] <- cov
mcmc1 <- epictmcmc(object = epi, datatype = "known epidemic", nsim = 150000, control.sus = suscov)The estimates of the model parameters can be then obtained either through using S3 summary function of the epictmcmc for the mcmc1, or using the summary function of coda package for mcmc1$parameter.samples, for example. The posterior means and 95% credible intervals of the model parameters using the former summary can be obtained via the following command:
summary(mcmc1, start = 10000)
*********************************************************
Model: SIR network-based continuous-time ILM
Method: Markov chain Monte Carlo (MCMC)
Data assumption: fully observed epidemic
number.chains : 1 chains
number.iteration : 140000 iterations
number.parameter : 2 parameters
*********************************************************
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
Alpha_s[1] 0.0851393 0.0268787 7.18362e-05 0.000279848
Alpha_s[2] 0.5090538 0.1284466 3.43287e-04 0.000936706
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
Alpha_s[1] 0.0409208 0.0656406 0.0825226 0.101275 0.144555
Alpha_s[2] 0.2820138 0.4187457 0.4997392 0.590065 0.783330
3. Empirical mean, standard deviation, and quantiles for the log likelihood,
Mean SD Naive SE Time-series SE
-55.80474998 1.01640137 0.00271644 0.00958428
2.5% 25% 50% 75% 97.5%
-58.5354 -56.1984 -55.4910 -55.0804 -54.8096
4. acceptance.rate :
Alpha_s[1] Alpha_s[2]
0.112127 0.223155 The MCMC trace plots for the model parameters can be produced using the S3 method plot.epictmcmc as follows:
plot(mcmc1, plottype = "parameter", start = 10000, thin = 10, density = FALSE)
- R - The Comprehensive R Archive Network
- Waleed Almutiry - Author and Maintainer
- Rob Deardon - Author
- Vineetha Warriyar K. V. - Contributer
This project is licensed under the GNU General Public License, version 3 - see the (http://www.r-project.org/Licenses/GPL-3) file for details.