--- title: "Working with Custom Attributes and Summary Statistics in EpiModel" author: "EpiModel v`r packageVersion('EpiModel')`" date: "`r Sys.Date()`" vignette: > %\VignetteIndexEntry{Working with Custom Attributes and Summary Statistics} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console --- ```{r, echo = FALSE, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup, echo = FALSE, include = FALSE} library(EpiModel) ``` ## Introduction This vignette covers how to work with nodal attributes, attribute histories, and epidemic summary statistics in EpiModel network models with custom extension modules. It assumes familiarity with setting up and running network models with `netsim()` and with the extension API. See the [Network Modeling for Epidemics](https://epimodel.github.io/sismid/) (NME) course materials and the [EpiModel Gallery](https://epimodel.github.io/EpiModel-Gallery/) for background. For working with network objects and edgelists, see the companion vignette *Working with Network Objects in EpiModel*. For parameter management, see *Working with Model Parameters in EpiModel*. In network simulations with `netsim`, all simulation state is stored in the *main list object* (referred to as `dat`). This vignette covers three types of data on `dat`: - **Current Nodal Attributes:** Per-node characteristics at the current time step. - **Historical Nodal Attributes:** Selective recordings of attribute values over time. - **Epidemic Trackers:** Population-level summary statistics computed at each time step. ## Current Nodal Attributes *Attributes* are characteristics of the nodes (e.g., persons) in the model at the current time step. They are stored as vectors in `dat$run$attr`, all of the same length as the number of nodes. ### Core Attributes Every node has four core attributes, managed automatically by EpiModel: - `active`: 1 if the node is currently in the population, 0 if departed. - `entrTime`: the time step when the node entered the population. - `exitTime`: the time step when the node departed (`NA` if still active). - `unique_id`: a globally unique integer identifier (never reused). You can retrieve the list of core attributes and their default types with `get_core_attributes()`: ```{r} get_core_attributes() ``` In addition, the built-in epidemic modules in core EpiModel use the `status` attribute (disease status, with values like `"s"`, `"i"`, `"r"`) and `infTime` (time of infection). More complex disease-specific EpiModel implementations (e.g., for HIV or COVID) may handle disease status differently. Custom attributes of any name and scalar value can be added, such as `age`, `race`, or `viral_load`. ### Positional Indexing vs Unique IDs Each node can be referenced in two ways: - **By position:** Think of it like a row number in a spreadsheet. `get_attr(dat, "active", posit_ids = 3)` accesses the third node's value directly. This is the standard way to look up node information and is very fast. - **By `unique_id`:** A globally unique integer attribute assigned at creation and never reused. Slower to look up, but allows referencing nodes that have already departed. Consider a simple example. A model starts with 5 nodes: | Position (*row number*) | `unique_id` | `active` | |:-----------------------:|:-----------:|:--------:| | 1 | 1 | 1 | | 2 | 2 | 1 | | 3 | 3 | 1 | | 4 | 4 | 1 | | 5 | 5 | 1 | Node 3 departs at step 10. The vectors still have 5 elements, but `active[3] = 0`: | Position (*row number*) | `unique_id` | `active` | |:-----------------------:|:-----------:|:--------:| | 1 | 1 | 1 | | 2 | 2 | 1 | | 3 | 3 | **0** | | 4 | 4 | 1 | | 5 | 5 | 1 | A new node arrives at step 15. In many simulations, departed nodes are dropped from the attribute vectors, so the new node takes position 3 (the freed slot) but receives a new `unique_id = 6`: | Position (*row number*) | `unique_id` | `active` | |:-----------------------:|:-----------:|:--------:| | 1 | 1 | 1 | | 2 | 2 | 1 | | 3 | **6** | 1 | | 4 | 4 | 1 | | 5 | 5 | 1 | Positional indices are used by current edgelists and attribute accessors. Unique IDs are used by cumulative edgelists and attribute histories because they can refer to nodes no longer in the model. See `help("unique_id-tools", package = "EpiModel")` for conversion functions. ### Working With Attributes #### Accessing Attributes Use `get_attr()` to extract an attribute vector: ```{r accessing-attr, eval = FALSE} active <- get_attr(dat, "active") ``` This returns a vector of length equal to the current number of nodes. Requesting a nonexistent attribute throws an error. #### Modifying Attributes Below is a minimal module that increments the age of all nodes by 1: ```{r modyfing-attr, eval = FALSE} aging_module <- function(dat, at) { # Extract current attribute age <- get_attr(dat, "age") # Aging process new_age <- age + 1 # Output updated attributes dat <- set_attr(dat, "age", new_age) return(dat) } ``` Key points about `set_attr()`: - It returns a modified copy of `dat`, which must be assigned back. (This does not cause performance issues due to R's shallow copy semantics since R 3.1.) - The replacement vector must have the same length as the number of nodes, or the function throws an error. #### Summary The recommended workflow for attribute manipulation in modules: 1. Extract attributes into local vectors with `get_attr()`. 2. Modify the local vectors with your dynamic process. 3. Update `dat` with the revised vectors using `set_attr()`. 4. Return `dat` at the end of the function. These functions have additional arguments described in `help("get_attr", package = "EpiModel")`. ## Historical Nodal Attributes Current attributes reflect the state of each node at the present time step. Sometimes you need to track what happened to specific nodes over time. Because saving the full history of all attributes would consume too much memory, EpiModel provides two complementary mechanisms: **automatic tracking** via `control.net` and **manual recording** inside custom modules. Both write to the same `attr.history` storage, accessible after the simulation with `get_attr_history()`. ### Automatic Attribute Tracking The simplest way to record attribute histories is to declare the attributes you want tracked in `control.net`. EpiModel will then automatically record changes at each time step using delta compression---only nodes whose values actually changed (including `NA` transitions), plus newly arrived and departed nodes, are stored. The `active` attribute is always tracked when any tracked attributes are specified. ```{r tracked_attrs_setup, eval = FALSE} control <- control.net( type = "SIS", nsteps = 500, nsims = 1, tracked.attributes = c("status"), tracked.attributes.once = c("race"), verbose = FALSE ) ``` Attributes listed in `tracked.attributes` are diff-tracked: at each step, only nodes whose values changed are recorded. Attributes listed in `tracked.attributes.once` are recorded only in the initial snapshot and for newly arriving nodes---they are not compared at each step. Use `tracked.attributes.once` for attributes that are fixed at entry and never change during the simulation (e.g., `"race"`, `"birth.year"`). After the simulation completes, you can reconstruct the tracked attribute state at any time step with `get_attr_at()`: ```{r get_attr_at_ex, eval = FALSE} sim <- netsim(est, param, init, control) # Reconstruct attributes at time step 200 attrs_200 <- get_attr_at(sim, at = 200) attrs_200 # # A tibble: 500 x 4 # unique_id active status race # # 1 1 1 s W # 2 2 1 i W # 3 3 1 s W # # ... with 497 more rows ``` The raw delta-compressed history is also available via `get_attr_history()`: ```{r get_attr_history_auto, eval = FALSE} attr_history <- get_attr_history(sim) names(attr_history) # [1] "active" "status" "race" attr_history$status # # A tibble: 11,812 x 5 # sim time attribute uids values # # 1 1 1 status 1 s # 2 1 1 status 2 s # 3 1 1 status 3 s # # ... with 11,809 more rows # After the initial snapshot of the 500 nodes, only the changes are recorded attr_history$status[501:503, ] # # A tibble: 3 x 5 # sim time attribute uids values # # 1 1 2 status 97 i # 2 1 2 status 123 i # 3 1 3 status 13 i ``` Each element is a `tibble` with columns `sim` (simulation number), `time` (time step), `attribute` (the label), `uids` (unique IDs), and `values` (the recorded values). The initial time step contains the full snapshot; subsequent steps contain only the changes. ### Manual Attribute Recording For finer-grained control---recording specific subsets of nodes, at specific intervals, or attributes not listed in `tracked.attributes`---you can call `record_attr_history()` directly inside a custom module. This writes to the same `attr.history` storage as automatic tracking. The following module records the viral load of infected nodes every 10 time steps. It assumes the model defines a `status` attribute (with values `"s"` for susceptible, `"i"` for infected) and a `viral_load` attribute (continuous for infected nodes, `NA` otherwise). ```{r recording_attr_hist, eval = FALSE} viral_load_logger_module <- function(dat, at) { # Run every 10 time steps if (at %% 10 == 0) { # Attributes status <- get_attr(dat, "status") viral_load <- get_attr(dat, "viral_load") infected <- which(status == "i") dat <- record_attr_history( dat, "viral_load", viral_load[infected], posit_ids = infected ) } return(dat) } ``` `record_attr_history()` takes the following arguments: 1. `dat`: the main list object. 2. `item`: a label for the attribute (here `"viral_load"`). 3. `value`: the values to record, one per node in `posit_ids` or `unique_ids`, or a single value recycled for all of them. 4. `at`: the time step to associate with the recording. Defaults to the current time step if omitted, and must be passed by name. 5. `posit_ids` or `unique_ids`: a vector of node IDs to which the values apply, passed by name. When `posit_ids` are provided, they are converted internally to `unique_id`s so the history remains valid even after nodes depart. Use attribute tracking and histories sparingly---storage can grow large in open-population models with many nodes over many time steps. Note that automatic tracking with `tracked.attributes` also records all nodes in the initial snapshot, but subsequent steps only store the changes, making it more memory-efficient for attributes that change infrequently relative to the population size. ## Epidemic Trackers *Epidemic trackers* are population-level summary statistics computed and stored at each time step---historical summary statistics about the population, not individual nodes. For example, the built-in tracker `num` records the total population size at each step, while `s.num`, `i.num`, and `r.num` record compartment sizes. ### Working With Epidemic Trackers in Modules Inside a module, access and modify epidemic trackers with `get_epi()` and `set_epi()`. Here is an extended aging module that also tracks mean age and its change over time: ```{r epi-in-modules, eval = FALSE} aging_track_module <- function(dat, at) { # Attributes age <- get_attr(dat, "age") # Aging process new_age <- age + 1 # Calculate summary statistics mean_age <- mean(new_age) # the tracker does not exist before the module first sets it prev_mean_age <- get_epi(dat, "mean_age", at - 1, override.null.error = TRUE) age_change <- if (is.null(prev_mean_age)) NA else mean_age - prev_mean_age # Update nodal attributes dat <- set_attr(dat, "age", new_age) # Update epidemic trackers dat <- set_epi(dat, "mean_age", mean_age) dat <- set_epi(dat, "age_change", age_change) return(dat) } ``` Key points: - `set_epi()` returns a modified `dat` that must be assigned back. The value must be a scalar. It is recorded at the current time step; to record it at another one, pass `at` by name (`set_epi(dat, "mean_age", mean_age, at = 10)`). - `get_epi()` retrieves the tracker value at a specific time step. Requesting a nonexistent tracker throws an error, unless `override.null.error = TRUE`, which returns `NULL` instead; the module above uses it on its first call, before `mean_age` exists. ### Accessing Epidemic Trackers After a Simulation Epidemic trackers are the main quantities of interest after a simulation. Access them by calling `as.data.frame()` on a `netsim` object, or use `plot()` and `summary()`: ```{r, eval = FALSE} sim <- netsim(est, param, init, control) # Raw per-simulation data as.data.frame(sim) # Means across simulations as.data.frame(sim, out = "mean") # Plot and summary plot(sim) summary(sim, at = 50) ``` ### Derived Statistics with `mutate_epi()` After a simulation, you can compute derived epidemic statistics without re-running the model. `mutate_epi()` works like `dplyr::mutate()` but operates on the `epi` list: ```{r, eval = FALSE} # Add incidence rate and prevalence sim <- mutate_epi(sim, prev = i.num / num) sim <- mutate_epi(sim, ir = (si.flow / s.num) * 100) # These new variables appear in the data.frame output as.data.frame(sim) ``` ### Stratified Prevalence with `epi.by` The `epi.by` parameter in `control.net()` automatically generates stratified compartment counts by a nodal attribute. For example, to track infection by a `"risk"` attribute: ```{r tracker-epiby, results = "hide", message = FALSE} set.seed(0) nw <- network_initialize(n = 100) nw <- set_vertex_attribute(nw, "risk", rep(0:1, each = 50)) est <- netest( nw, formation = ~edges + nodefactor("risk"), target.stats = c(40, 20), coef.diss = dissolution_coefs(~offset(edges), 20), # duration = 20 verbose = FALSE ) param <- param.net(inf.prob = 0.3, act.rate = 1) init <- init.net(i.num = 5) control <- control.net( type = "SI", nsims = 1, nsteps = 50, epi.by = "risk", verbose = FALSE ) sim <- netsim(est, param, init, control) ``` With `epi.by = "risk"`, the output includes stratified columns automatically: ```{r} d <- as.data.frame(sim) names(d) ``` The naming pattern is `{tracker}.{attribute}{value}` (e.g., `i.num.risk0`, `i.num.risk1`, `num.risk0`, `num.risk1`). ### Custom Epidemic Trackers For summary statistics beyond built-in compartment counts, you can define custom tracker functions and pass them to `control.net()`. This avoids writing a full module for simple calculations. #### Writing Tracker Functions A tracker function must take **exactly one argument** (`dat`) and return a **scalar** value (numeric, logical, or character). It is called automatically at the end of each time step, after all modules have run. ```{r 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) } ``` Here is a tracker that computes the proportion of the population infected: ```{r tracker-commented} epi_prop_infected <- function(dat) { needed_attributes <- c("status", "active") output <- with(get_attr_list(dat, needed_attributes), { pop <- active == 1 infected <- status == "i" sum(infected & pop, na.rm = TRUE) / sum(pop, na.rm = TRUE) }) return(output) } ``` We recommend this pattern for tracker functions: 1. Define `needed_attributes` as a vector of attribute names. 2. Use `with(get_attr_list(dat, needed_attributes), { ... })` to work in a clean environment. 3. The last expression inside `with` becomes the output. It **must** be a scalar. 4. Return the output. #### Using Custom Tracker Functions Pass tracker functions as a named list via `.tracker.list` in `control.net()`. The names become column names in the output `data.frame`: ```{r tracker-list} some.trackers <- list( prop_infected = epi_prop_infected, s_num = epi_s_num ) control <- control.net( type = "SI", nsims = 1, nsteps = 50, verbose = FALSE, .tracker.list = some.trackers ) param <- 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 <- netest( nw, formation = ~edges, target.stats = 25, coef.diss = dissolution_coefs(~offset(edges), 10, 0), # duration = 10, closed population verbose = FALSE ) init <- init.net(i.num = 10) sim <- netsim(est, param, init, control) d <- as.data.frame(sim) knitr::kable(tail(d, n = 15)) ``` Note that in the `.tracker.list` we pass the function objects (`epi_prop_infected`, `epi_s_num`) without parentheses---we are storing the functions, not calling them. Trackers are **always** run at the end of a simulation step. The values they output reflect the state of the simulation *after* all modules have completed their work for that step. ## Recording Arbitrary Objects (Advanced Debugging) When debugging complex research models, you may want to inspect the state of an object during simulation without stopping it. `record_raw_object()` saves any R object at a given time step for later inspection. ```{r record-raw, eval = FALSE} introspect_module <- function(dat, at) { age <- get_attr(dat, "age") if (mean(age, na.rm = TRUE) > 50) { obj <- data.frame( age = age, status = get_attr(dat, "status") ) dat <- record_raw_object(dat, "old pop", obj) } return(dat) } ``` After simulation, access recorded objects from the `netsim` output: ```{r record-raw-access, eval = FALSE} sim <- netsim(est, param, init, control) # Access raw records for simulation 1, as a list of records with the # elements `time`, `label`, and `object` raw <- sim[["raw.records"]][[1]]$as_list() # Access the "old pop" data.frame recorded at step 75 old_pop <- Filter(function(r) r$label == "old pop" && r$time == 75, raw) old_pop[[1]]$object ``` **Use sparingly:** Raw records are not truncated and can consume substantial memory, especially for large attribute vectors recorded at many time steps. This feature is intended for debugging, not routine output collection.