Measles Outbreak Example 2

library(multigroup.vaccine)
library(epiworldR)
#> Thank you for using epiworldR! Please consider citing it in your work.
#> You can find the citation information by running
#>   citation("epiworldR")

schpops <- c(592, 334, 32, 211, 657, 796, 728, 1305, 386, 447, 625, 706, 306, 27,
  1002, 560, 27, 124, 692, 190, 607, 448, 567, 529, 1107, 457, 597, 438, 420, 690,
  806, 461, 537, 1177, 1263, 200, 741, 793, 555, 930, 1424, 409, 635, 486, 408,
  1119, 19, 736, 24, 469, 42, 526, 847, 1234, 504, 289, 37, 439, 911, 618, 432, 546)

schvacc <- c(543, 194, 29, 132, 567, 726, 676, 1151, 299, 411, 541, 619, 291, 19,
  831, 486, 26, 118, 638, 80, 524, 398, 517, 508, 1026, 417, 537, 379, 372, 600,
  740, 425, 444, 1016, 1122, 165, 653, 744, 475, 874, 1229, 345, 473, 424, 318,
  870, 18, 643, 12, 404, 42, 449, 785, 1060, 461, 143, 34, 390, 753, 547, 360, 471)

onlpops <- 1612
onlvacc <- 1400

# 0-4, 18-24 25-44, 45-64, 65 plus
agepops <- c(11772, 19219, 47701, 41451, 45347)
agecovr <- c(0.84, 0.9, 0.92, 0.95, 1)
agevacc <- round(agepops * agecovr)

pops <- c(schpops, onlpops, agepops)
initV <- c(schvacc, onlvacc, agevacc)

initGrp <- length(schpops) + 4
initI <- rep(0, length(pops))
initI[initGrp] <- 1
initR <- rep(0, length(pops))

R0 <- 15
meaninf <- 7
incontact <- rep(0.2, length(pops))
f <- (1 - incontact) * pops
cmat <- (diag(incontact) + outer((1 - incontact), f / sum(f)))

betaij <- transmissionRates(
  R0 = R0,
  meaninf = meaninf,
  reltransm = cmat
)

ode_size <- multigroup.vaccine::getFinalSizeODE(
  transmrates = betaij,
  recoveryrate = 1 / meaninf,
  popsize = pops,
  initR = initR,
  initI = initI,
  initV = initV
)

fsODE <- sum(ode_size$totalSize)

numsims <- 1000
size_dist <- getFinalSizeDist(
  n = numsims,
  transmrates = betaij,
  recoveryrate = 1 / meaninf,
  popsize = pops,
  initR = initR,
  initV = initV,
  initI = initI)

fs1 <- rowSums(size_dist)
hist(fs1, breaks = 50)