## load the 'imdepi' data and a model fit
data("imdepi")
data("imdepifit")
## calculate individual and type-specific reproduction numbers
R0s <- R0(imdepifit)
tapply(R0s, imdepi$events@data[names(R0s), "type"], summary)
## untrimmed R0 for a specific event
R0(imdepifit, newevents=marks(imdepi)[1,], trimmed=FALSE)
### compute a Monte Carlo confidence interval
## use a simpler model with constant 'siaf' for speed
simplefit <- update(imdepifit, epidemic=~type, siaf=NULL, subset = NULL)
## we'd like to compute the mean R0's by event type
meanR0ByType <- function (newcoef) {
R0events <- R0(simplefit, newcoef=newcoef)
tapply(R0events, imdepi$events@data[names(R0events),"type"], mean)
}
(meansMLE <- meanR0ByType(newcoef=NULL))
## sample B times from asymptotic multivariate normal of the MLE
B <- 5 # CAVE: toy example! In practice this has to be much larger
set.seed(123)
parsamples <- MASS::mvrnorm(B, mu=coef(simplefit), Sigma=vcov(simplefit))
## for each sample compute the 'meanR0ByType'
meansMC <- apply(parsamples, 1, meanR0ByType)
## get the quantiles and print the result
cisMC <- apply(cbind(meansMLE, meansMC), 1, quantile, probs=c(0.025,0.975))
print(rbind(MLE=meansMLE, cisMC))
### R0 for a simple epidemic model (homogeneous individuals)
mepi1 <- update(simplefit, epidemic = ~1, data = imdepi,
subset = type == "B", model = TRUE, verbose = FALSE)
## using the default spatial and temporal ranges of interaction
(R0B <- simpleR0(mepi1)) # eps.s=200, eps.t=30
stopifnot(identical(R0B, R0(mepi1, trimmed = FALSE)[[1]]))
## assuming less interaction
simpleR0(mepi1, eps.s = 50, eps.t = 15)
Run the code above in your browser using DataLab