tests/testthat/test-netcensus.R

In EpiModel: Mathematical Modeling of Infectious Disease Dynamics

context("Observed Network Layers")

# Shared fixtures ------------------------------------------------------------

# A static census: a fixed observed network with a nodal attribute
set.seed(41)
N <- 60
nw_static <- network_initialize(n = N)
nw_static <- set_vertex_attribute(nw_static, "sex", rep(c("F", "M"), N / 2))
pairs <- t(combn(N, 2))
pairs <- pairs[sample.int(nrow(pairs), 90), ]
nw_static <- add.edges(nw_static, tail = pairs[, 1], head = pairs[, 2])

# A dynamic census: edge spells over an observation window of 0 to 30. Some
# dyads start before the window (censored onsets), some are still active at
# its end, and a third of the dyads have a second, later spell, as recurring
# contacts in proximity-logger data do.
make_dynamic_census <- function(n = 60, n.edges = 150, end = 30) {
  pairs <- t(combn(n, 2))
  pairs <- pairs[sample.int(nrow(pairs), n.edges), ]
  onset <- sample(c(-Inf, 0:(end - 5)), n.edges, replace = TRUE)
  terminus <- pmax(ifelse(is.finite(onset), onset, 0), 0) +
    sample(2:12, n.edges, replace = TRUE)
  terminus[terminus >= end] <- Inf
  spells <- data.frame(onset = onset, terminus = terminus,
                       tail = pairs[, 1], head = pairs[, 2])
  again <- which(is.finite(terminus))[seq_len(n.edges %/% 3)]
  second <- spells[again, ]
  second$onset <- second$terminus + 2
  second$terminus <- pmin(second$onset + 3, Inf)
  second <- second[second$onset < end, ]
  spells <- rbind(spells, second)
  base <- network::network.initialize(n, directed = FALSE)
  base <- set_vertex_attribute(base, "sex", rep(c("F", "M"), n / 2))
  nd <- networkDynamic::networkDynamic(
    base.net = base,
    edge.spells = spells
  )
  nd %n% "net.obs.period" <- list(observations = list(c(0, end)),
                                  mode = "discrete", time.increment = 1,
                                  time.unit = "step")
  nd
}
nw_dyn <- suppressMessages(make_dynamic_census())

active_at <- function(nd, at) {
  el <- networkDynamic::get.dyads.active(nd, at = at)
  el <- cbind(pmin(el[, 1], el[, 2]), pmax(el[, 1], el[, 2]))
  el[order(el[, 1], el[, 2]), , drop = FALSE]
}

init <- init.net(i.num = 5)

# networkDynamic() prints the observation period it creates
quiet_nd <- function(...) {
  out <- NULL
  utils::capture.output(out <- suppressMessages(
    networkDynamic::networkDynamic(...)
  ))
  out
}


# Constructor ------------------------------------------------------------------

test_that("netcensus wraps a static network", {
  obs <- netcensus(nw_static)
  expect_s3_class(obs, "netcensus")
  expect_false(obs$dynamic)
  expect_null(obs$window)
  expect_false(obs$edapprox)
  expect_null(obs$target.stats)
  expect_equal(obs$summary$edges, 90)
  expect_equal(obs$summary$n, N)
  expect_equal(obs$summary$mean.degree, 2 * 90 / N)
  expect_equal(obs$summary$attributes, "sex")
  out <- capture.output(print(obs))
  expect_true(any(grepl("static \\(network\\)", out)))
  expect_true(any(grepl("Edges: 90", out)))
})

test_that("netcensus wraps a dynamic network and finds its window", {
  obs <- netcensus(nw_dyn)
  expect_true(obs$dynamic)
  expect_equal(obs$window, c(0, 30))
  expect_equal(obs$summary$edges.ever, network.edgecount(nw_dyn))
  # the fixture has dyads with more than one spell
  sp <- networkDynamic::get.edge.activity(nw_dyn, as.spellList = TRUE)
  expect_gt(nrow(sp), network.edgecount(nw_dyn))
  # the window is half-open, so its end is not an observed step
  ea <- obs$summary$edges.active
  expect_equal(names(ea), as.character(0:29))
  expect_equal(obs$summary$steps, c(0, 29))
  expect_equal(unname(ea["10"]), nrow(active_at(nw_dyn, 10)))
  out <- capture.output(print(obs))
  expect_true(any(grepl("dynamic \\(networkDynamic\\)", out)))
  expect_true(any(grepl("Observation window: [0, 30), time steps 0 to 29",
                        out, fixed = TRUE)))

  # without net.obs.period, the window is the range of finite spell times
  nd <- nw_dyn
  network::delete.network.attribute(nd, "net.obs.period")
  obs2 <- netcensus(nd)
  sp <- networkDynamic::get.edge.activity(nd, as.spellList = TRUE)
  on <- sp$onset[is.finite(sp$onset)]
  term <- sp$terminus[is.finite(sp$terminus)]
  expect_equal(obs2$window, c(min(on, term - 1), max(on + 1, term)))

  # an explicit window overrides
  obs3 <- netcensus(nw_dyn, window = c(5, 20))
  expect_equal(obs3$window, c(5, 20))
  expect_equal(names(obs3$summary$edges.active), as.character(5:19))
})

test_that("netcensus finds a finite window when net.obs.period is infinite", {
  # the networkDynamic() constructor sets an infinite net.obs.period for open
  # spells (terminus Inf) and censored onsets (-Inf)
  base <- network::network.initialize(10, directed = FALSE)
  spells <- data.frame(onset = c(1, 3, -Inf), terminus = c(5, Inf, 4),
                       tail = c(1, 2, 3), head = c(2, 3, 4))
  nd <- quiet_nd(base.net = base, edge.spells = spells)
  expect_false(all(is.finite(unlist((nd %n% "net.obs.period")$observations))))
  obs <- netcensus(nd)
  expect_equal(obs$window, c(1, 5))
  expect_equal(obs$summary$steps, c(1, 4))
  expect_equal(unname(obs$summary$edges.active),
               vapply(1:4, function(t) {
                 NROW(networkDynamic::get.dyads.active(nd, at = t))
               }, numeric(1)))

  expect_error(netcensus(nd, window = c(0, Inf)), "finite times")
  expect_error(netcensus(nd, window = c(1.2, 1.8)), "no whole time step")

  # a census without any finite spell time needs an explicit window
  nd2 <- quiet_nd(
    base.net = base,
    edge.spells = data.frame(onset = -Inf, terminus = Inf, tail = 1, head = 2)
  )
  expect_error(netcensus(nd2), "Set `window`")
  expect_equal(netcensus(nd2, window = c(0, 10))$summary$steps, c(0, 9))

  # a dynamic census without edges is refused like a static one
  nd3 <- networkDynamic::as.networkDynamic(base)
  expect_true(networkDynamic::is.networkDynamic(nd3))
  expect_error(netcensus(nd3), "no edges")
})

test_that("netcensus refuses loops and parallel edges", {
  # add.edges() changes its network in place, so build each one anew
  flagged <- function(tail = c(1, 2), head = c(2, 3)) {
    nw <- network::network.initialize(5, directed = FALSE, loops = TRUE,
                                      multiple = TRUE)
    network::add.edges(nw, tail = tail, head = head)
  }
  expect_error(netcensus(flagged(c(1, 2, 4), c(2, 3, 4))), "loops")
  expect_error(netcensus(flagged(c(1, 2, 2), c(2, 3, 1))), "parallel edges")
  nw <- flagged()

  # a network only flagged as allowing them runs under tergmLite
  obs <- netcensus(nw)
  expect_false(obs$newnetwork %n% "loops")
  expect_false(obs$newnetwork %n% "multiple")
  control <- control.net(type = "SI", nsteps = 3, nsims = 1, tergmLite = TRUE,
                         resimulate.network = TRUE, verbose = FALSE)
  sim <- netsim(obs, param.net(inf.prob = 0.5), init.net(i.num = 1), control)
  expect_true(all(get_nwstats(sim, network = 1)$edges == 2))
})

test_that("netcensus drops temporally extended vertex attributes", {
  nd <- nw_dyn
  nd <- networkDynamic::activate.vertex.attribute(nd, "status", "s",
                                                  onset = 0, terminus = Inf)
  expect_true("status.active" %in% list.vertex.attributes(nd))
  expect_message(obs <- netcensus(nd), "status.active")
  expect_false("status.active" %in% list.vertex.attributes(obs$newnetwork))
  expect_equal(obs$summary$attributes, "sex")
})

test_that("trim_netest converts a static census and keeps a dynamic one", {
  obs <- netcensus(nw_static)
  tr <- trim_netest(obs)
  expect_s3_class(tr, "netcensus")
  expect_s3_class(tr$newnetwork, "networkLite")
  expect_equal(tr$summary, obs$summary)
  expect_equal(unclass(as.edgelist(tr$newnetwork))[, 1:2],
               unclass(as.edgelist(nw_static))[, 1:2])
  dyn <- suppressMessages(netcensus(nw_dyn))
  expect_identical(trim_netest(dyn), dyn)

  # the trimmed static census simulates in either storage mode
  param <- param.net(inf.prob = 0.3, act.rate = 1)
  for (tergmLite in c(TRUE, FALSE)) {
    control <- control.net(type = "SI", nsteps = 5, nsims = 1,
                           tergmLite = tergmLite,
                           resimulate.network = tergmLite, verbose = FALSE)
    sim <- netsim(tr, param, init, control)
    expect_true(all(get_nwstats(sim, network = 1)$edges == 90))
  }
})

test_that("netcensus validates its inputs", {
  expect_error(netcensus(list()), "class `network`")
  expect_error(netcensus(network::network.initialize(10, directed = TRUE)),
               "undirected")
  expect_error(netcensus(network_initialize(10)), "no edges")
  expect_error(netcensus(nw_static, window = c(0, 10)), "networkDynamic")
  expect_error(netcensus(nw_dyn, window = c(10, 0)), "start < end")
  expect_error(netcensus(nw_dyn, window = 1), "start < end")
})

test_that("netdx refuses a netcensus object", {
  expect_error(netdx(netcensus(nw_static), nsims = 1, nsteps = 5),
               "does not apply")
})

test_that("make_stats_table accepts NULL targets", {
  stats <- list(matrix(c(1, 2, 3), ncol = 1, dimnames = list(NULL, "edges")))
  tab <- make_stats_table(stats, NULL)
  expect_true(is.na(tab$Target))
  expect_equal(tab[["Sim Mean"]], 2)
})


# Simulation: static census -----------------------------------------------------

test_that("netsim runs on a static census, tergmLite and networkDynamic", {
  obs <- netcensus(nw_static)
  param <- param.net(inf.prob = 0.3, act.rate = 1)
  for (tergmLite in c(TRUE, FALSE)) {
    control <- control.net(type = "SI", nsteps = 15, nsims = 1,
                           tergmLite = tergmLite,
                           resimulate.network = tergmLite, verbose = FALSE,
                           save.transmat = TRUE, save.run = TRUE)
    sim <- netsim(obs, param, init, control)
    expect_s3_class(sim, "netsim")
    expect_s3_class(sim$nwparam[[1]], "netcensus")
    test_net(sim)

    # the edge set is the observed one throughout
    ns <- get_nwstats(sim, network = 1)
    expect_equal(nrow(ns), 15)
    expect_true(all(ns$edges == 90))
    el <- if (tergmLite) sim$run[[1]]$el[[1]] else {
      active_at(get_network(sim), 15)
    }
    expect_equal(unclass(el)[, 1:2], unclass(as.edgelist(nw_static))[, 1:2],
                 ignore_attr = TRUE)

    # every transmission crossed an observed edge
    tm <- get_transmat(sim)
    obs_el <- as.edgelist(nw_static)
    keys <- paste(obs_el[, 1], obs_el[, 2])
    expect_true(all(paste(pmin(tm$sus, tm$inf), pmax(tm$sus, tm$inf)) %in% keys))

    # the nodal attribute came through
    expect_equal(sim$run[[1]]$attr$sex, rep(c("F", "M"), N / 2))

    # print and plot show a target of NA for the observed layer
    expect_output(print(sim), "netcensus")
    plot(sim)
    plot(sim, type = "formation")
  }
})


# Simulation: dynamic census ----------------------------------------------------

test_that("netsim reads the active edges of a dynamic census each step", {
  skip_on_cran()
  obs <- netcensus(nw_dyn)
  param <- param.net(inf.prob = 0.3, act.rate = 1)
  for (tergmLite in c(TRUE, FALSE)) {
    control <- control.net(type = "SI", nsteps = 25, nsims = 1,
                           tergmLite = tergmLite,
                           resimulate.network = tergmLite, verbose = FALSE,
                           save.transmat = TRUE, save.run = TRUE,
                           cumulative.edgelist = TRUE,
                           save.cumulative.edgelist = TRUE,
                           truncate.el.cuml = Inf)
    sim <- netsim(obs, param, init, control)
    test_net(sim)

    # the recorded edge count at each step is the census count at that time
    ns <- get_nwstats(sim, network = 1)
    census_counts <- vapply(1:25, function(t) nrow(active_at(nw_dyn, t)),
                            numeric(1))
    expect_equal(ns$edges, census_counts)

    # under tergmLite the final edgelist is the census at the last step
    if (tergmLite) {
      el <- sim$run[[1]]$el[[1]]
      expect_equal(unclass(el)[, 1:2], active_at(nw_dyn, 25),
                   ignore_attr = TRUE)
      expect_equal(attr(el, "n"), N)
    } else {
      # without tergmLite the observed spells are untouched
      expect_identical(networkDynamic::get.edge.activity(get_network(sim)),
                       networkDynamic::get.edge.activity(nw_dyn))
    }

    # every transmission crossed an edge active at that time
    tm <- get_transmat(sim)
    ok <- vapply(seq_len(nrow(tm)), function(i) {
      el <- active_at(nw_dyn, tm$at[i])
      any(el[, 1] == min(tm$sus[i], tm$inf[i]) &
            el[, 2] == max(tm$sus[i], tm$inf[i]))
    }, logical(1))
    expect_true(all(ok))

    # the cumulative edgelist follows the census: an edge that dissolved in
    # the census has a stop time
    cel <- sim$cumulative.edgelist[[1]]
    expect_true(any(!is.na(cel$stop)))
    expect_true(all(cel$start >= 0))
  }
})

test_that("a dynamic census works next to estimated and clique layers", {
  skip_on_cran()
  nw <- network_initialize(N)
  nw <- set_vertex_attribute(nw, "sex", rep(c("F", "M"), N / 2))
  nw <- set_vertex_attribute(nw, "hh_id", rep(1:20, each = 3))
  est <- netest(nw, ~edges, target.stats = 30,
                coef.diss = dissolution_coefs(~offset(edges), 10),
                verbose = FALSE)
  hh <- netclique(nw, group.attr = "hh_id")
  obs <- netcensus(nw_dyn)
  param <- param.net(inf.prob = multilayer(0.1, 0.3, 0.2), act.rate = 1)
  for (tergmLite in c(TRUE, FALSE)) {
    control <- control.net(type = "SI", nsteps = 20, nsims = 1,
                           tergmLite = tergmLite, resimulate.network = TRUE,
                           verbose = FALSE, save.transmat = TRUE,
                           save.run = TRUE)
    sim <- netsim(list(est, obs, hh), param, init, control)
    test_net(sim)
    expect_s3_class(sim$nwparam[[2]], "netcensus")
    expect_s3_class(sim$nwparam[[3]], "netclique")
    # only the estimated layer has a coefficient
    expect_length(sim$nwparam[[1]]$coef.form, 1)
    expect_null(sim$nwparam[[2]]$coef.form)
    tm <- get_transmat(sim)
    expect_true(all(tm$network %in% 1:3))
    expect_true(all(tm$transProb[tm$network == 2] == 0.3))
    ns <- get_nwstats(sim, network = 2)
    expect_equal(ns$edges,
                 vapply(1:20, function(t) nrow(active_at(nw_dyn, t)),
                        numeric(1)))
  }
})

test_that("a census layer refuses vital dynamics and node changes", {
  obs <- netcensus(nw_dyn)
  param <- param.net(inf.prob = 0.3, a.rate = 0.01, ds.rate = 0.01,
                     di.rate = 0.01)
  control <- control.net(type = "SI", nsteps = 10, nsims = 1,
                         tergmLite = TRUE, resimulate.network = TRUE,
                         verbose = FALSE)
  expect_error(netsim(obs, param, init, control), "fixed node set")

  # a custom module that adds nodes is stopped in arrive_nodes
  control <- control.net(type = "SI", nsteps = 5, nsims = 1, tergmLite = TRUE,
                         resimulate.network = TRUE, verbose = FALSE)
  dat <- crosscheck.net(obs, param.net(inf.prob = 0.3), init, control)
  dat <- initialize.net(obs, param.net(inf.prob = 0.3), init, control, 1)
  expect_error(arrive_nodes(dat, 1), "fixed node set")
  expect_error(depart_nodes(dat, 1), "fixed node set")
  expect_s3_class(arrive_nodes(dat, 0), "netsim_dat")
})

test_that("duration tracking and runs past the window are caught", {
  obs <- netcensus(nw_dyn)
  param <- param.net(inf.prob = 0.3)
  control <- control.net(type = "SI", nsteps = 10, nsims = 1,
                         tergmLite = TRUE, resimulate.network = TRUE,
                         tergmLite.track.duration = TRUE, verbose = FALSE)
  expect_error(netsim(obs, param, init, control), "track.duration")

  control <- control.net(type = "SI", nsteps = 40, nsims = 1,
                         tergmLite = TRUE, resimulate.network = TRUE,
                         verbose = FALSE)
  expect_warning(sim <- netsim(obs, param, init, control),
                 "observation window")
  ns <- get_nwstats(sim, network = 1)
  # once the last finite spell has ended, the edge set is frozen
  sp <- networkDynamic::get.edge.activity(nw_dyn, as.spellList = TRUE)
  last <- ceiling(max(sp$terminus[is.finite(sp$terminus)]))
  expect_lt(last, 40)
  expect_true(all(ns$edges[last:40] == ns$edges[last]))
  expect_equal(ns$edges[40], nrow(active_at(nw_dyn, 40)))

  # a static census has no window and no warning
  control <- control.net(type = "SI", nsteps = 40, nsims = 1,
                         tergmLite = TRUE, resimulate.network = TRUE,
                         verbose = FALSE)
  expect_silent(sim <- netsim(netcensus(nw_static), param, init, control))
})

test_that("a census with vertex spells and censored onsets runs without tergmLite", {
  skip_on_cran()
  # vertex activity spells that start after 0 and end before the last edge,
  # as in networkDynamicData::concurrencyComparisonNets; netsim ignores them
  nd <- nw_dyn
  nd <- networkDynamic::activate.vertices(nd, onset = 2, terminus = 28)
  obs <- netcensus(nd)
  param <- param.net(inf.prob = 0.3)
  for (resim in c(FALSE, TRUE)) {
    control <- control.net(type = "SI", nsteps = 25, nsims = 1,
                           tergmLite = FALSE, resimulate.network = resim,
                           verbose = FALSE, save.nwstats = TRUE)
    sim <- netsim(obs, param, init, control)
    ns <- get_nwstats(sim, network = 1)
    expect_equal(ns$edges,
                 vapply(1:25, function(t) nrow(active_at(nd, t)), numeric(1)))
    test_net(sim)
  }
})

test_that("a census layer after the first one keeps its vertex attributes", {
  nw <- network_initialize(N)
  est <- netest(nw, ~edges, target.stats = 30,
                coef.diss = dissolution_coefs(~offset(edges), 10),
                verbose = FALSE)
  obs <- netcensus(nw_static)
  param <- param.net(inf.prob = 0.3, act.rate = 1)
  for (tergmLite in c(TRUE, FALSE)) {
    control <- control.net(type = "SI", nsteps = 3, nsims = 1,
                           tergmLite = tergmLite, resimulate.network = TRUE,
                           verbose = FALSE, save.run = TRUE, epi.by = "sex")
    sim <- netsim(list(est, obs), param, init, control)
    expect_equal(sim$run[[1]]$attr$sex, rep(c("F", "M"), N / 2))
    expect_true(all(c("i.num.sexF", "i.num.sexM") %in% names(sim$epi)))
  }

  # an attribute on both layers must agree
  nw2 <- set_vertex_attribute(nw, "sex", rep(c("M", "F"), N / 2))
  est2 <- netest(nw2, ~edges, target.stats = 30,
                 coef.diss = dissolution_coefs(~offset(edges), 10),
                 verbose = FALSE)
  control <- control.net(type = "SI", nsteps = 3, nsims = 1, tergmLite = TRUE,
                         resimulate.network = TRUE, verbose = FALSE)
  expect_error(netsim(list(est2, obs), param, init, control),
               "differs from the one on census layer 2")
})

test_that("a dynamic census with empty steps runs without tergmLite", {
  # contact data often have steps without contacts; the edgelist of such a
  # step must still be a matrix for the cumulative edgelist
  base <- network::network.initialize(20, directed = FALSE)
  spells <- data.frame(onset = 2, terminus = 6, tail = 1:5, head = 6:10)
  nd <- quiet_nd(base.net = base, edge.spells = spells)
  obs <- netcensus(nd, window = c(0, 10))
  for (tergmLite in c(TRUE, FALSE)) {
    control <- control.net(type = "SI", nsteps = 8, nsims = 1,
                           tergmLite = tergmLite, resimulate.network = tergmLite,
                           verbose = FALSE, cumulative.edgelist = TRUE,
                           save.cumulative.edgelist = TRUE,
                           truncate.el.cuml = Inf)
    sim <- netsim(obs, param.net(inf.prob = 0.5), init.net(i.num = 2), control)
    cel <- sim$cumulative.edgelist[[1]]
    expect_equal(nrow(cel), 5)
    expect_true(all(cel$start == 2 & cel$stop == 5))
    expect_equal(get_nwstats(sim, network = 1)$edges,
                 c(0, 5, 5, 5, 5, 0, 0, 0))
  }
})

test_that("netsim warns when the simulation ends before a census window", {
  base <- network::network.initialize(20, directed = FALSE)
  spells <- data.frame(onset = 100, terminus = 130, tail = 1:5, head = 6:10)
  nd <- quiet_nd(base.net = base, edge.spells = spells)
  obs <- netcensus(nd, window = c(100, 130))
  control <- control.net(type = "SI", nsteps = 20, nsims = 1, tergmLite = TRUE,
                         resimulate.network = TRUE, verbose = FALSE)
  expect_warning(netsim(obs, param.net(inf.prob = 0.5), init.net(i.num = 2),
                        control),
                 "ends at step 20, before the observation window")
})

test_that("make_restart_point refuses a dynamic census layer", {
  obs <- netcensus(nw_dyn)
  control <- control.net(type = "SI", nsteps = 5, nsims = 1,
                         tergmLite = TRUE, resimulate.network = TRUE,
                         verbose = FALSE, save.run = TRUE)
  sim <- netsim(obs, param.net(inf.prob = 0.3), init, control)
  expect_error(make_restart_point(sim, time_attrs = c()),
               "does not support dynamic `netcensus` layers")
  # a static census has no clock to follow
  sim <- netsim(netcensus(nw_static), param.net(inf.prob = 0.3), init, control)
  expect_s3_class(make_restart_point(sim, time_attrs = c()), "netsim")
})

test_that("a census layer survives a restart", {
  skip_on_cran()
  obs <- netcensus(nw_dyn)
  param <- param.net(inf.prob = 0.3)
  control <- control.net(type = "SI", nsteps = 10, nsims = 1,
                         tergmLite = TRUE, resimulate.network = TRUE,
                         verbose = FALSE, save.run = TRUE)
  sim <- netsim(obs, param, init, control)
  control$start <- 11
  control$nsteps <- 20
  sim2 <- netsim(sim, param, init, control)
  expect_equal(nrow(sim2$epi$num), 20)
  ns <- get_nwstats(sim2, network = 1)
  expect_equal(ns$edges[11:20],
               vapply(11:20, function(t) nrow(active_at(nw_dyn, t)),
                      numeric(1)))
  test_net(sim2)
})