infect <- function(dat, at) { ## Uncomment this to run environment interactively # browser() ## Attributes ## active <- get_attr(dat, "active") status <- get_attr(dat, "status") infTime <- get_attr(dat, "infTime") ## Parameters ## inf.prob <- get_param(dat, "inf.prob") act.rate <- get_param(dat, "act.rate") ## Find infected nodes ## idsInf <- which(active == 1 & status == "i") nActive <- sum(active == 1) nElig <- length(idsInf) ## Initialize default incidence at 0 ## nInf <- 0 ## If any infected nodes, proceed with transmission ## if (nElig > 0 && nElig < nActive) { ## Look up discordant edgelist ## del <- discord_edgelist(dat, at) ## If any discordant pairs, proceed ## if (!(is.null(del))) { # Set parameters on discordant edgelist data frame del$transProb <- inf.prob del$actRate <- act.rate del$finalProb <- 1 - (1 - del$transProb)^del$actRate # Stochastic transmission process transmit <- rbinom(nrow(del), 1, del$finalProb) # Keep rows where transmission occurred del <- del[which(transmit == 1), ] # Look up new ids if any transmissions occurred idsNewInf <- unique(del$sus) nInf <- length(idsNewInf) # Set new attributes and transmission matrix if (nInf > 0) { status[idsNewInf] <- "e" infTime[idsNewInf] <- at dat <- set_attr(dat, "status", status) dat <- set_attr(dat, "infTime", infTime) dat <- set_transmat(dat, del, at) } } } ## Save summary statistic for S->E flow dat <- set_epi(dat, "se.flow", at, nInf) return(dat) } progress <- function(dat, at) { ## Uncomment this to function environment interactively # browser() ## Attributes ## active <- get_attr(dat, "active") status <- get_attr(dat, "status") ## Parameters ## ei.rate <- get_param(dat, "ei.rate") ir.rate <- get_param(dat, "ir.rate") ## E to I progression process ## nInf <- 0 idsEligInf <- which(active == 1 & status == "e") nEligInf <- length(idsEligInf) if (nEligInf > 0) { vecInf <- which(rbinom(nEligInf, 1, ei.rate) == 1) if (length(vecInf) > 0) { idsInf <- idsEligInf[vecInf] nInf <- length(idsInf) status[idsInf] <- "i" } } ## I to R progression process ## nRec <- 0 idsEligRec <- which(active == 1 & status == "i") nEligRec <- length(idsEligRec) if (nEligRec > 0) { vecRec <- which(rbinom(nEligRec, 1, ir.rate) == 1) if (length(vecRec) > 0) { idsRec <- idsEligRec[vecRec] nRec <- length(idsRec) status[idsRec] <- "r" } } ## Write out updated status attribute ## dat <- set_attr(dat, "status", status) ## Save summary statistics ## dat <- set_epi(dat, "ei.flow", at, nInf) dat <- set_epi(dat, "ir.flow", at, nRec) dat <- set_epi(dat, "e.num", at, sum(active == 1 & status == "e")) dat <- set_epi(dat, "r.num", at, sum(active == 1 & status == "r")) return(dat) } # ---- Aging ---- # # A time step in this model is one DAY, so one step of aging advances every # node by 1/365 of a year. Age is kept as a continuous attribute rather than # as an age group because two other parts of the model read it directly: the # age-specific mortality lookup in dfunc below, and the absdiff("age") term # in the formation formula. # # meanAge is a tracker, not part of the process. Aging is the only force # pushing mean age up; arrivals entering at age 0 and age-skewed mortality # both push it down. Watching meanAge is how we see which dominates over a # long run, and confirms the age structure is behaving rather than drifting # in some unintended way. aging <- function(dat, at) { age <- get_attr(dat, "age") age <- age + 1/365 dat <- set_attr(dat, "age", age) dat <- set_epi(dat, "meanAge", at, mean(age, na.rm = TRUE)) return(dat) } dfunc <- function(dat, at) { ## Attributes active <- get_attr(dat, "active") exitTime <- get_attr(dat, "exitTime") age <- get_attr(dat, "age") status <- get_attr(dat, "status") ## Parameters dep.rates <- get_param(dat, "departure.rates") dep.dis.mult <- get_param(dat, "departure.disease.mult") ## Query alive idsElig <- which(active == 1) nElig <- length(idsElig) ## Initialize trackers nDepts <- 0 idsDepts <- NULL if (nElig > 0) { ## Calculate age-specific departure rates for each eligible node ## ## Everyone older than 85 gets the final mortality rate # # dep.rates is a plain vector of 86 daily rates, one per year of age, and # R vectors are 1-indexed. So a continuous age in years has to become a # 1-based position: ceiling() maps an age in (0, 1] to index 1, an age in # (1, 2] to index 2, and so on. pmin(..., 86) puts everyone past 85 on # the final rate rather than running off the end of the vector. # # There is a trap worth knowing here. ceiling(0) is 0, and R SILENTLY # DROPS index 0 from a vector subscript rather than raising an error, so # a node aged exactly 0 would shorten departure_rates_of_elig by one # element. rbinom() would then recycle the short vector and assign every # subsequent node the wrong mortality rate, with no warning at all. # # That never happens here, because of module ORDER rather than anything # in this function. Arrivals enter at age 0, but aging() has already run # this step and moved them to 1/365, so dfunc() never sees an age of 0. # The order holds for a slightly subtle reason: aging.FUN is not one of # control.net's own arguments, so it is a USER module, and EpiModel runs # all user modules before all built-in ones (departures.FUN is built-in). # Shuffling the arguments in the control.net call therefore cannot break # it, but setting the module.order control explicitly could. whole_ages_of_elig <- pmin(ceiling(age[idsElig]), 86) departure_rates_of_elig <- dep.rates[whole_ages_of_elig] ## Multiply departure rates for diseased persons idsElig.inf <- which(status[idsElig] == "i") departure_rates_of_elig[idsElig.inf] <- departure_rates_of_elig[idsElig.inf] * dep.dis.mult ## Simulate departure process vecDepts <- which(rbinom(nElig, 1, departure_rates_of_elig) == 1) idsDepts <- idsElig[vecDepts] nDepts <- length(idsDepts) ## Update nodal attributes if (nDepts > 0) { active[idsDepts] <- 0 exitTime[idsDepts] <- at } } ## Set updated attributes dat <- set_attr(dat, "active", active) dat <- set_attr(dat, "exitTime", exitTime) ## Summary statistics ## dat <- set_epi(dat, "total.deaths", at, nDepts) # covid deaths covid.deaths <- length(intersect(idsDepts, which(status == "i"))) dat <- set_epi(dat, "covid.deaths", at, covid.deaths) return(dat) } afunc <- function(dat, at) { ## Parameters ## # Note this reads the population size at the PREVIOUS step. The current # step's "num" has not been written yet: the prevalence module, which # records it, is a built-in that runs near the end of the step, after # arrivals. So at - 1 is the most recent count available here. n <- get_epi(dat, "num", at - 1) a.rate <- get_param(dat, "arrival.rate") ## Process ## nArrivalsExp <- n * a.rate nArrivals <- rpois(1, nArrivalsExp) # Update attributes. The ORDER of these calls is load-bearing. # # append_core_attr must come first. It extends EpiModel's own bookkeeping # attributes (active, entrTime, exitTime, and the unique node ids) and sets # the new population size. Only once the population has grown can # append_attr add the model-specific attributes for those same new nodes. # # After that, there must be one append_attr call for EVERY attribute the # model tracks. Omit one and its vector stays short while the others grow, # so every attribute vector is misaligned against the node ids from that # point on, with no error raised. This is exactly what the demography lab # exercises when it asks you to add an age.grp attribute: adding it to the # network is not enough, the arrivals module has to append it too. if (nArrivals > 0) { dat <- append_core_attr(dat, at = at, n.new = nArrivals) dat <- append_attr(dat, "status", "s", nArrivals) dat <- append_attr(dat, "infTime", NA, nArrivals) dat <- append_attr(dat, "age", 0, nArrivals) } ## Summary statistics ## dat <- set_epi(dat, "a.flow", at, nArrivals) return(dat) }