inst/doc/attributes-and-summary-statistics.R

In EpiModel: Mathematical Modeling of Infectious Disease Dynamics

## ---- echo = FALSE, include = FALSE-------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, echo = FALSE, include = FALSE-------------------------------------
library(EpiModel)

## ----accessing-attr, eval = FALSE---------------------------------------------
#  active <- EpiModel::get_attr(dat, "active")

## ----modyfing-attr, eval = FALSE----------------------------------------------
#  aging_module <- function(dat, at) {
#  
#    # Extract current attribute
#    age <- EpiModel::get_attr(dat, "age")
#  
#    # Aging process
#    new_age <- age + 1
#  
#    # Output updated attributes
#    dat <- EpiModel::set_attr(dat, "age", new_age)
#  
#    return(dat)
#  }

## ----recordin_attr_hist, eval = FALSE-----------------------------------------
#  viral_load_logger_module <- function(dat, at) {
#  
#    # Run every 10 time steps
#    if (at %% 10 == 0) {
#  
#      # Attributes
#      status <- EpiModel::get_attr(dat, "status")
#      viral_load <- EpiModel::get_attr(dat, "viral_load")
#  
#      infected <- which(status == "infected")
#  
#      dat <- EpiModel::record_attr_history(
#        dat, at,
#        "viral_load",
#        infected,
#        viral_load[infected]
#      )
#    }
#  
#    # Output
#    return(dat)
#  }

## ----access_attr_hist, eval = FALSE-------------------------------------------
#  sim <- netsim(est, param, init, control)
#  attr_history <- EpiModel::get_attr_history(sim)

## ----get_attr_ex, eval = FALSE------------------------------------------------
#  get_attr_history(sim)
#  
#  # $viral_load
#  #    sim step attribute    uids values
#  # 1  1   10   viral_load   1001   2000
#  # 2  1   10   viral_load   1002   1878
#  # 3  1   20   viral_load   1001   1500
#  # 4  1   20   viral_load   1002    300
#  # 5  2   10   viral_load    401   2500
#  # 6  2   10   viral_load    402   1378
#  # 7  2   20   viral_load    401   1200
#  # 8  2   20   viral_load    402    100
#  # ...
#  #
#  # $status
#  #    sim step attribute      uids     values
#  # 1  1   22   status  1001   infected
#  # 2  1   64   status  1002   infected
#  # 3  1  110   status  1001  recovered
#  # 4  1  220   status  1002  recovered
#  # 5  2    7   status   401   infected
#  # 6  2   15   status   402   infected
#  # 7  2   20   status   401  recovered
#  # 8  2  120   status   402  recovered
#  # ...

## ----epi-in-modules, eval = FALSE---------------------------------------------
#  aging_track_module <- function(dat, at) {
#  
#    # Attributes
#    age <- EpiModel::get_attr(dat, "age")
#  
#    # Aging process
#    new_age <- age + 1
#  
#    # Calculate summary statistics
#    mean_age <- mean(new_age)
#    prev_mean_age <- EpiModel::get_epi(dat, at - 1, "mean_age")
#    age_change <- mean_age - prev_mean_age
#  
#    # Update nodal attributes
#    dat <- EpiModel::set_attr(dat, "age", new_age)
#  
#    # Update epidemic trackers
#    dat <- set_epi(dat, "mean_age", at, mean_age)
#    dat <- set_epi(dat, "age_change", at, age_change)
#  
#    return(dat)
#  }

## ----tracker-example----------------------------------------------------------
epi_s_num <- function(dat) {
  needed_attributes <- c("status")
  output <- with(get_attr_list(dat, needed_attributes), {
    sum(status == "s", na.rm = TRUE)
  })
  return(output)
}

## ----tracker-commented--------------------------------------------------------
epi_prop_infected <- function(dat) {
  # we need two attributes for our calculation: `status` and `active`
  needed_attributes <- c("status", "active")
  # we use `with` to simplify code
  output <- with(EpiModel::get_attr_list(dat, needed_attributes), {
    pop <- active == 1    # we only look at active nodes
    cond <- status == "i" # which are infected

    # how many are `infected` among the `active`
    sum(cond & pop, na.rm = TRUE) / sum(pop, na.rm = TRUE)
  })
  return(output)
}

## ----tracker-list-------------------------------------------------------------
# Create the `tracker.list` list
some.trackers <- list(
  prop_infected = epi_prop_infected,
  s_num         = epi_s_num
)

control <- EpiModel::control.net(
  type = "SI",
  nsims = 1,
  nsteps = 50,
  verbose = FALSE,
  .tracker.list = some.trackers
)

param <- EpiModel::param.net(
  inf.prob = 0.3,
  act.rate = 0.1
)

nw <- network_initialize(n = 50)
nw <- set_vertex_attribute(nw, "race", rbinom(50, 1, 0.5))
est <- EpiModel::netest(
  nw,
  formation = ~edges,
  target.stats = 25,
  coef.diss = dissolution_coefs(~offset(edges), 10, 0),
  verbose = FALSE
)

init <- EpiModel::init.net(i.num = 10)
sim <- EpiModel::netsim(est, param, init, control)

d <- as.data.frame(sim)

knitr::kable(tail(d, n = 15))

## ----record-raw, eval = FALSE-------------------------------------------------
#  introspect_module <- function(dat, at) {
#    # Attributes
#    age <- get_attr(dat, "age")
#  
#    if (mean(age, na.rm = TRUE) > 50) {
#      obj <- data.frame(
#          age = age,
#          status = EpiModel::get_attr(dat, "status")
#      )
#      dat <- EpiModel::record_raw_object(dat, at, "old pop", obj)
#    }
#  
#    return(dat)
#  }

## ----record-raw-access, eval = FALSE------------------------------------------
#  sim <- netsim(est, param, init, control)
#  sim[["raw.records"]][[1]]