The following document presents a set of notes from "EpiModel: An R Package for Mathematical Modeling of Infectious Disease over Networks".
Another interesting example is provided here
The process is as follows:
Initialize and set the network calling network::network.initialize
,
Estimate the network parameters calling netest
passing network effects, this will be used to simulate the dynamic network, and initialize the simulation calling init.net
(initial infected, recovered, etc.)
Define the set of parameters that will be used during the simulation using param.net
Define the control functions for the simulation (the modules) using control.net
Run your simulation by calling netsim
and passing your estimates from netest
, the output from param.net
and from control.net
Internal work of netsim
from EpiModel
# Loop through steps for (at in max(2, control$start):control$nsteps) { # Fetching the ordering in which modules happen morder <- control$module.order if (is.null(morder)) { lim.bi.mods <- control$bi.mods[-which(control$bi.mods %in% c("initialize.FUN", "verbose.FUN"))] morder <- c(control$user.mods, lim.bi.mods) } # Applying the models (here is where the simulation actually happens) for (i in seq_along(morder)) { dat <- do.call(control[[morder[i]]], list(dat, at)) } # Printing results on screen if (!is.null(control[["verbose.FUN"]])) { do.call(control[["verbose.FUN"]], list(dat, type = "progress", s, at)) } }
In this section, I review some of the examples presented in the document. The subsections are titled as the subsection in the original document. We first need to load the R packages that we will be used.
library(EpiModel) library(sna) library(networkDynamic)
The following lines of code introduce a new module that sets the mortality rate as a function of age (which is not included by default)
aging <- function(dat, at) { if (at == 2) { # Initializing age n <- sum(dat$attr$active == 1) dat$attr$age <- sample( seq(from = 18, to = 69 + 11 / 12, by = 1 / 12), n, replace = TRUE) } else { # Increasing age dat$attr$age <- dat$attr$age + 1 / 12 } # Saving some stats: Average age if (at == 2) { dat$epi$meanAge <- c(NA, mean(dat$attr$age, na.rm = TRUE)) } else { dat$epi$meanAge[at] <- mean(dat$attr$age, na.rm = TRUE) } return(dat) }
This example shows how to create a new mortality extension, which is to replace the default function passed to control.net
as the argument deaths.FUN
. Currently the default is deaths.net
.
ages <- 18:69 death.rates <- 1/(70*12 - ages*12) dfunc <- function(dat, at) { # Finding active nodes idsElig <- which(dat$attr$active == 1) nElig <- length(idsElig) nDeaths <- 0 # If there are active nodes if (nElig > 0) { # Fetching their ages (from attr) ages <- dat$attr$age[idsElig] # The maximun age is fetched from params max.age <- dat$param$max.age death.rates <- pmin(1, 1 / (max.age * 12 - ages * 12)) vecDeaths <- which(rbinom(nElig, 1, death.rates) == 1) idsDeaths <- idsElig[vecDeaths] nDeaths <- length(idsDeaths) # If there are deads, then: # - Set them as inactive, # - Set the exit time, # - And deactivate the vertices (and edges) in the network if (nDeaths > 0) { dat$attr$active[idsDeaths] <- 0 dat$attr$exitTime[idsDeaths] <- at dat$nw <- deactivate.vertices( dat$nw, onset = at, terminus = Inf, v = idsDeaths, deactivate.edges = TRUE ) } } # Adding summary stats if (at == 2) { # In the first run, the statistic d.flow is created dat$epi$d.flow <- c(0, nDeaths) } else { # Then it is just filled dat$epi$d.flow[at] <- nDeaths } return(dat) }
The following module replaces the birth.FUN
(which by default is births.net
):
bfunc <- function(dat, at) { # Getting some parameters growth.rate <- dat$param$growth.rate exptPopSize <- dat$epi$num[1] * (1 + growth.rate * at) n <- network.size(dat$nw) numNeeded <- exptPopSize - sum(dat$attr$active == 1) # If new births are needed, then simulate them using a poisson if (numNeeded > 0) { nBirths <- rpois(1, numNeeded) } else { nBirths <- 0 } # If there were new birdths, then add them to the network and # activate them if (nBirths > 0) { dat$nw <- add.vertices(dat$nw, nv = nBirths) newNodes <- (n + 1):(n + nBirths) dat$nw <- activate.vertices( dat$nw, onset = at, terminus = Inf, v = newNodes ) } # Initializing vertices atributes (if any) if (nBirths > 0) { dat$attr$active <- c(dat$attr$active, rep(1, nBirths)) dat$attr$status <- c(dat$attr$status, rep("s", nBirths)) dat$attr$infTime <- c(dat$attr$infTime, rep(NA, nBirths)) dat$attr$age <- c(dat$attr$age, rep(18, nBirths)) } # Saving stats if (at == 2) { dat$epi$b.flow <- c(0, nBirths) } else { dat$epi$b.flow[at] <- nBirths } return(dat) }
The following new module computes and stores the mean degree through the simulation, notice that in the case of networkDynamic
, the time starts at 0.
dfun <- function(dat, at) { if (at == 2) { dat$epi$meandeg <- c(0, mean(degree(network.extract(dat$nw, at = at - 1L)))) } else { dat$epi$meandeg[at] <- mean(degree(network.extract(dat$nw, at = at - 1L))) } dat }
# Step 1: Bernoulli network nw <- network::network.initialize(500, directed = FALSE) # Step 2 est3 <- netest( nw, formation = ~edges, target.stats = 150, coef.diss = dissolution_coefs(~offset(edges), 60, mean(death.rates)) ) # Step 3: Passing the parameters param <- param.net(inf.prob = 0.15, growth.rate = 0.01/12, max.age = 70) init <- init.net(i.num = 50) # Step 4: Passing the new functions control <- control.net( type = "SI", nsims = 1, nsteps = 100, deaths.FUN = dfunc, births.FUN = bfunc, aging.FUN = aging, depend = TRUE, verbose = FALSE, # This function is new in -control.net- deg.FUN = dfun ) # Step 5: Running the simulation sim3 <- netsim(est3, param, init, control)
# If our simulation worked, then we should be able to replicate the # mean degree degs <- sapply(1:100, function(i) { mean(degree(network.extract(sim3$network$sim1, at = i-1))) }) cor( sim3$epi$meandeg[-1,1], degs[-1] ) # This should be one!
In this example the authors describe the a small set of modifications that can be done to the function infection.net
so that there's an intermediate step in the process, exposed, the initial stage of the infection.
The code of infection.net
follows:
infection.net <- function (dat, at) { # Vector of attributes active <- dat$attr$active status <- dat$attr$status # Parameters of the network modes <- dat$param$modes mode <- idmode(dat$nw) inf.prob <- dat$param$inf.prob inf.prob.m2 <- dat$param$inf.prob.m2 act.rate <- dat$param$act.rate # The actual network nw <- dat$nw tea.status <- dat$control$tea.status # Which ids are susceptible and infected? idsSus <- which(active == 1 & status == "s") idsInf <- which(active == 1 & status == "i") # How many are active? nActive <- sum(active == 1) nElig <- length(idsInf) nInf <- nInfM2 <- totInf <- 0 # If there are elegible, then... if (nElig > 0 && nElig < nActive) { # Finding edges in which there are inf and suscept # this returns a data frame with ids of infect and suscept. del <- discord_edgelist(dat, idsInf, idsSus, at) # If there are any if (!(is.null(del))) { # How long have the infected been infected # This is relevant as the probability of infecting is a function # of this time del$infDur <- at - dat$attr$infTime[del$inf] del$infDur[del$infDur == 0] <- 1 linf.prob <- length(inf.prob) # Specifying infection probabilities from mode 1 to 2 (if not specified) if (is.null(inf.prob.m2)) { del$transProb <- ifelse(del$infDur <= linf.prob, inf.prob[del$infDur], inf.prob[linf.prob]) } # Otherwise, specify these according to the type of edge else { del$transProb <- ifelse(del$sus <= nw %n% "bipartite", ifelse(del$infDur <= linf.prob, inf.prob[del$infDur], inf.prob[linf.prob]), ifelse(del$infDur <= linf.prob, inf.prob.m2[del$infDur], inf.prob.m2[linf.prob])) } if (!is.null(dat$param$inter.eff) && at >= dat$param$inter.start) { del$transProb <- del$transProb * (1 - dat$param$inter.eff) } lact.rate <- length(act.rate) del$actRate <- ifelse(del$infDur <= lact.rate, act.rate[del$infDur], act.rate[lact.rate]) del$finalProb <- 1 - (1 - del$transProb)^del$actRate transmit <- rbinom(nrow(del), 1, del$finalProb) del <- del[which(transmit == 1), ] idsNewInf <- unique(del$sus) nInf <- sum(mode[idsNewInf] == 1) nInfM2 <- sum(mode[idsNewInf] == 2) totInf <- nInf + nInfM2 if (totInf > 0) { if (tea.status == TRUE) { nw <- activate.vertex.attribute(nw, prefix = "testatus", value = "i", onset = at, terminus = Inf, v = idsNewInf) } dat$attr$status[idsNewInf] <- "i" dat$attr$infTime[idsNewInf] <- at form <- get_nwparam(dat)$formation fterms <- get_formula_terms(form) if ("status" %in% fterms) { nw <- set.vertex.attribute(nw, "status", dat$attr$status) } } if (any(names(nw$gal) %in% "vertex.pid")) { del$sus <- get.vertex.pid(nw, del$sus) del$inf <- get.vertex.pid(nw, del$inf) } } } if (totInf > 0) { del <- del[!duplicated(del$sus), ] if (at == 2) { dat$stats$transmat <- del } else { dat$stats$transmat <- rbind(dat$stats$transmat, del) } } if (at == 2) { dat$epi$si.flow <- c(0, nInf) if (modes == 2) { dat$epi$si.flow.m2 <- c(0, nInfM2) } } else { dat$epi$si.flow[at] <- nInf if (modes == 2) { dat$epi$si.flow.m2[at] <- nInfM2 } } dat$nw <- nw return(dat) }
contagion <- function(dat, at) { }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.