Working with Custom Attributes and Summary Statistics in EpiModel

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():

get_core_attributes()
#> $active
#> [1] 1
#> 
#> $entrTime
#> [1] 1
#> 
#> $exitTime
#> [1] NA
#> 
#> $unique_id
#> [1] 1

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:

Node 3 departs at step 10. The vectors still have 5 elements, but active[3] = 0:

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:

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

active <- get_attr(dat, "active")

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)
}

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.

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():

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
#       <int>  <int> <chr>  <chr>
# 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():

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
#   <dbl> <int> <chr>     <int> <chr>
# 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
#   <dbl> <int> <chr>     <int> <chr>
# 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).

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_ids 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.

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:

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():

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:

# 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:

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:

d <- as.data.frame(sim)
names(d)
#>  [1] "sim"         "time"        "s.num"       "s.num.risk0" "s.num.risk1"
#>  [6] "i.num"       "i.num.risk0" "i.num.risk1" "num"         "num.risk0"  
#> [11] "num.risk1"   "si.flow"

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.

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)
}

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)
}

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
)
#> Starting simulated annealing (SAN)
#> Iteration 1 of at most 4
#> Finished simulated annealing
#> Starting maximum pseudolikelihood estimation (MPLE):
#> Obtaining the responsible dyads.
#> Evaluating the predictor and response matrix.
#> Maximizing the pseudolikelihood.
#> Finished MPLE.

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

d <- as.data.frame(sim)

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