data("hagelloch")
plot(hagelloch)
# a simplistic twinSIR model
fit <- twinSIR(~ household, data = hagelloch)
# overall total intensity
plot(fit, which = "total")
# overall epidemic proportion
epi <- plot(fit, which = "epidemic", ylim = c(0, 1))
head(epi)
# add overall endemic proportion = 1 - epidemic proportion
ende <- plot(fit, which = "endemic", add = TRUE, col = 2)
legend("topleft", legend = "endemic proportion", lty = 1, col = 2, bty = "n")
# individual intensities
tmp <- plot(fit, which = "total", aggregate = FALSE,
col = rgb(0, 0, 0, alpha = 0.1),
main = expression("Individual infection intensities " *
lambda[i](t) == Y[i](t) %.% (e[i](t) + h[i](t))))
# return value: matrix of individual intensity paths
str(tmp)
# plot intensity path only for individuals 3 and 99
matplot(x = tmp[,1], y = tmp[,1+c(3,99)], type = "S",
ylab = "Force of infection", xlab = "time",
main = expression("Paths of the infection intensities " *
lambda[3](t) * " and " * lambda[99](t)))
legend("topright", legend = paste("Individual", c(3,99)),
col = 1:2, lty = 1:2)
Run the code above in your browser using DataLab