## Generate a data frame containing a hypothetic population with 100 individuals
set.seed(1234)
n <- 100
pos <- matrix(rnorm(n*2), ncol=2, dimnames=list(NULL, c("x", "y")))
pop <- data.frame(id=1:n, x=pos[,1], y=pos[,2],
gender=sample(0:1, n, replace=TRUE),
I0col=c(rep(1,3),rep(0,n-3)), # 3 initially infectious
start=rep(0,n), stop=rep(Inf,n))
## Simulate an SIR epidemic in this population
set.seed(123)
infPeriods <- setNames(c(1:3/10, rexp(n-3, rate=1)), 1:n)
epi <- simEpidata(
cbind(start,stop) ~ cox(gender), data = pop,
id.col = "id", I0.col = "I0col", coords.cols = c("x","y"),
beta = c(-2), h0 = -1, alpha = c(B1=0.1), f = list(B1=function(u) u<=1),
infPeriod = function(ids) infPeriods[ids],
##remPeriod = function(ids) rexp(length(ids), rate=0.1), end = 30 # -> SIRS
)
## extract event times by id
head(summary(epi)$byID)
## Plot the numbers of susceptible, infectious and removed individuals
plot(epi)
## load the 1861 Hagelloch measles epidemic
data("hagelloch")
summary(hagelloch)
plot(hagelloch)
## fit a simplistic twinSIR model
fit <- twinSIR(~ household, data = hagelloch)
## simulate a new epidemic from the above model
## with simulation period = observation period, re-using observed infPeriods
sim1 <- simulate(fit, data = hagelloch)
plot(sim1)
## check if we find similar parameters in the simulated epidemic
fitsim1 <- update(fit, data = sim1)
cbind(base = coef(fit), new = coef(fitsim1))
if (surveillance.options("allExamples")) {
## simulate only 10 days, using random infPeriods ~ Exp(0.1)
sim2 <- simulate(fit, data = hagelloch, seed = 2, end = 10,
infPeriod = function(ids) rexp(length(ids), rate = 0.1))
plot(sim2)
## simulate from a different model with manually specified parameters
set.seed(321)
simepi <- simEpidata(~ cox(AGE), data = hagelloch,
beta = c(0.1), h0 = -4, alpha = c(household = 0.05),
f = list(household = function(u) u == 0),
infPeriod = function(ids) rexp(length(ids), rate=1/8))
plot(simepi)
intensityplot(simepi)
## see if we correctly estimate the parameters
fitsimepi <- twinSIR(~ cox(AGE) + household, data = simepi)
cbind(true = c(0.05, -4, 0.1), est = coef(fitsimepi), confint(fitsimepi))
}
Run the code above in your browser using DataLab