inst/doc/Predict.R

In EpiILM: Spatial and Network Based Individual Level Models for Epidemics

### R code from vignette source 'Predict.Rnw'

###################################################
### code chunk number 1: Predict.Rnw:40-43
###################################################
options(continue="  ", width=60)
pdf.options(pointsize=8) #text in graph about the same as regular text
library(EpiILM, quietly=TRUE)


###################################################
### code chunk number 2: Predict.Rnw:47-54
###################################################
set.seed(101)
n <- 100
contact <- matrix(0, n, n)
for(i in 1:(n-1)) {
   contact[i, ((i+1):n)] <- rbinom((n-i), 1, 0.05)
   contact[((i+1):n), i] <- contact[i, ((i+1):n)]
}


###################################################
### code chunk number 3: Predict.Rnw:65-68
###################################################
set.seed(101)
netdat <- epidata( type = "SI", n = 100, tmax = 25,  sus.par = 0.1, 
                   contact = contact)


###################################################
### code chunk number 4: Predict.Rnw:73-77
###################################################
t_end <- max(netdat$inftime)
mcmcout_net <- epimcmc(object = netdat,  tmax = t_end,
        niter = 50000, sus.par.ini = 0.01,pro.sus.var = 0.01,
        prior.sus.dist = "uniform", prior.sus.par = c(0, 10000))


###################################################
### code chunk number 5: nettrace
###################################################
plot(mcmcout_net, partype = "parameter", start = 10001, density = FALSE)


###################################################
### code chunk number 6: Predict.Rnw:89-96
###################################################
set.seed(1001)
pred.net.15 <- pred.epi(netdat, xx = mcmcout_net, tmin = 15, 
                       burnin = 1000, criterion = "newly infectious", 
                       n.samples = 500)
pred.net.20 <- pred.epi(netdat, xx = mcmcout_net, tmin = 20, 
                       burnin = 1000, criterion = "newly infectious", 
                       n.samples = 500)


###################################################
### code chunk number 7: t15
###################################################
plot(pred.net.15, col = "red", lwd = 2, pch = 19)


###################################################
### code chunk number 8: t20
###################################################
plot(pred.net.20, col = "red", lwd = 2, pch = 19)


###################################################
### code chunk number 9: Predict.Rnw:123-131
###################################################
# simulate spatial locations
set.seed(101)
x <- runif(100, 0, 10)
y <- runif(100, 0, 10)
# simulate covariate, number of animals on each farm
A <- round(rexp(100,1/50))
SI.cov <- epidata(type = "SI", n = 100, tmax = 25, x = x, y = y,
                 Sformula = ~A, sus.par = c(0.5, 0.5), beta = 6)


###################################################
### code chunk number 10: tM1
###################################################
t_end <- max(SI.cov$inftime)
prior_par <- matrix(rep(1, 4), ncol = 2, nrow = 2)
mcmcout_M1 <- epimcmc(SI.cov, Sformula = ~A, tmax = t_end, niter = 50000, 
                      sus.par.ini = c(0.001, 0.001), beta.ini = 0.01, 
                      pro.sus.var = c(0.08, 0.4), pro.beta.var = 0.5, 
                      prior.sus.dist = c("gamma","gamma"), 
                      prior.sus.par = prior_par, 
                      prior.beta.dist = "uniform", prior.beta.par = c(0, 10000) )
summary(mcmcout_M1, start = 10001)
#  MCMC traceplot for the estimation of the model (2) parameters
plot(mcmcout_M1, partype = "parameter", start = 10001, density = FALSE)


###################################################
### code chunk number 11: Predict.Rnw:161-167
###################################################
set.seed(101)
mcmcout_M2 <- epimcmc(SI.cov, tmax = t_end, niter = 50000, sus.par.ini = 0.01, 
                      beta.ini = 0.01,  pro.sus.var = 0.1, pro.beta.var = 0.5, 
                      prior.sus.dist = "uniform",
                      prior.sus.par = c(0, 10000), prior.beta.dist  = "uniform",
                      prior.beta.par = c(0, 10000))


###################################################
### code chunk number 12: tM2
###################################################
summary(mcmcout_M2,  start = 10001)  
# MCMC traceplot for the estimation of the model (3) parameters
plot(mcmcout_M2, partype = "parameter", start = 10001, density = FALSE)


###################################################
### code chunk number 13: Predict.Rnw:179-187
###################################################
set.seed(101)
pred.model1 <- pred.epi(SI.cov, Sformula = ~A, xx = mcmcout_M1, 
                    criterion = "newly infectious",  n.samples = 500)
# convert ILM model (2) into epidata object
model2.data <- as.epidata(type = "SI", n = 100, x = x, y = y, 
                          inftime = SI.cov$inftime)
pred.model2 <- pred.epi(model2.data, xx = mcmcout_M2, 
                     criterion = "newly infectious", n.samples = 500)


###################################################
### code chunk number 14: predictM1
###################################################
plot(pred.model1, col = "red", type = "b", lwd = 2)


###################################################
### code chunk number 15: predictM2
###################################################
plot(pred.model2, col = "red", type = "b", lwd = 2)


###################################################
### code chunk number 16: Predict.Rnw:209-214
###################################################
loglike <- epilike(SI.cov, tmax = t_end, Sformula = ~A, sus.par = c(0.597, 1.071),
                                         beta = 6.606)
dic1 <- epidic(burnin = 10000, niter = 50000, LLchain = mcmcout_M1$Loglikelihood,
                                         LLpostmean = loglike)
dic1


###################################################
### code chunk number 17: Predict.Rnw:218-221
###################################################
loglike <- epilike(model2.data, tmax = t_end,  sus.par = 6.210, beta = 4.942)
dic2 <- epidic(burnin = 10000, niter = 50000, LLchain = mcmcout_M2$Loglikelihood,
                                          LLpostmean = loglike)