diff --git a/DESCRIPTION b/DESCRIPTION index eac6171..0ba7828 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ARTnet -Version: 2.9.1 -Date: 2026-06-08 +Version: 2.9.2 +Date: 2026-08-21 Title: Epidemic and Network Model Parameterization with the ARTnet Dataset Description: Epidemic and Network Model Parameterization with the ARTnet Dataset. Maintainer: Samuel Jenness diff --git a/NEWS.md b/NEWS.md index 4e025d7..07bd818 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,28 @@ +## ARTnet v2.9.2 + +Two fixes to the sexual-cessation configuration of `build_netstats()`, where +`age.limits[2]` is set above `age.sexual.cessation` so the population keeps +aging past the point where it stops forming partnerships. + +* `active.sex` is now zeroed on the age group the breaks declare rather than on + the largest age group realized in the sampled population. When the + post-cessation group came out empty, the old behavior marked the oldest + *sexually active* group inactive and zeroed its degree on all three layers, + with no error and no warning. This affected both the sampling path and the + `target_pop` data frame path, where a projected population whose oldest + member sits below the cessation age would have hit the same problem. +* `young.prop` must now be strictly between 0 and 1 when a cessation band + exists, and is rejected with a message pointing at a usable value. Either + endpoint empties one side of the age split by construction, and an age group + with no members carries no `ergm` term, so the failure used to surface much + later as a `target.stats` length mismatch in `netest()`. A post-cessation + group that comes out empty for other reasons (a `network.size` too small for + the seeded share, or a `target_pop` with no members above the cessation age) + now warns for the same reason. + +Neither change affects the default configuration, where the population upper +limit and the cessation age coincide and there is no post-cessation band. + ## ARTnet v2.9.0 * `ARTnetData` is now a suggested package and example `EpiStats` and `Netstats` diff --git a/R/NetStats.R b/R/NetStats.R index db06bea..92ccc4d 100644 --- a/R/NetStats.R +++ b/R/NetStats.R @@ -191,6 +191,10 @@ #' @param young.prop The proportion of the population that should be below the age of sexual cessation. #' Default is NULL (meaning no re-weighting of the `age.pyramid` parameter is performed). #' This parameter is only used if the age of sexual cessation is less than the upper age bound. +#' Must be strictly between 0 and 1: either endpoint empties one side of the split, and an +#' age group with no members cannot carry an `ergm` term. To start a model with effectively +#' nobody above the cessation age, use a value just below 1 (for example 0.999) rather than +#' 1 itself, and let the post-cessation group fill over the burn-in. #' @param target_pop Optional specification of the synthetic target population. Defaults to NULL, #' which uses the legacy patchwork of reference sources (NCHS age pyramid + #' `ARTnetData::race.dist` + ARTnet's own degree / role / risk-quintile distributions). @@ -373,6 +377,30 @@ build_netstats <- function(epistats, netparams, time.unit <- epistats$time.unit + # `young.prop` splits the age pyramid at the sexual cessation age, giving + # `young.prop` of the sampling mass to the sexually active ages and the + # remainder to the post-cessation ages. The two endpoints are degenerate: + # each empties one side of the split entirely, and an age group with no + # members cannot carry an `ergm` term, because levels are derived from the + # values actually present on the network. `netest()` would then fail on a + # target.stats length mismatch, a long way from the cause. Reject them here. + if (sex.cess.mod == TRUE && !is.null(young.prop)) { + if (young.prop >= 1) { + stop("`young.prop` must be less than 1 when the age of sexual cessation ", + "is below the upper age limit. At 1 no node is placed in the ", + "post-cessation age group, so that group's `nodefactor` and ", + "`nodematch` terms would be missing from any model fitted on the ", + "resulting population. Use a value just below 1 (for example ", + "0.999, which seeds about 0.1% of `network.size` above the ", + "cessation age) and let the group fill over the model burn-in.") + } + if (young.prop <= 0) { + stop("`young.prop` must be greater than 0. At 0 the entire population ", + "is placed above the age of sexual cessation, leaving no sexually ", + "active nodes to fit a network model to.") + } + } + # Demographic Initialization ---------------------------------------------- @@ -510,6 +538,14 @@ build_netstats <- function(epistats, netparams, age.breaks <- out$demog$age.breaks <- epistats$age.breaks nquants <- length(netparams$inst$nf.risk.grp) + # The post-cessation age group, as declared by the age breaks rather than as + # realized in the sampled population. `build_epistats` puts + # `age.sexual.cessation` into `age.breaks`, so under `sex.cess.mod` the last + # group is exactly the post-cessation band. Reading the group off the data + # instead (`max(attr_age.grp)`) silently marks the oldest *sexually active* + # group inactive whenever the post-cessation band happens to be empty. + top.age.grp <- length(age.breaks) - 1L + # List-form distribution overrides for the sampling path. NULL means # "use the existing default source" (current behavior). .ov <- if (.tp$form == "list") .tp$overrides else list() @@ -531,7 +567,7 @@ build_netstats <- function(epistats, netparams, } else { as_active <- rep(1L, num) if (sex.cess.mod == TRUE) { - as_active[attr_age.grp == max(attr_age.grp, na.rm = TRUE)] <- 0L + as_active[attr_age.grp %in% top.age.grp] <- 0L } as_active } @@ -598,7 +634,7 @@ build_netstats <- function(epistats, netparams, # sexually active attribute attr_active.sex <- rep(1L, num) if (sex.cess.mod == TRUE) { - attr_active.sex[attr_age.grp == max(attr_age.grp)] <- 0L + attr_active.sex[attr_age.grp %in% top.age.grp] <- 0L } # race attribute @@ -640,6 +676,21 @@ build_netstats <- function(epistats, netparams, attr_role.class <- apportion_lr(num, 0:2, .dist_role.class, shuffled = TRUE) } + # A post-cessation band can still come out empty for reasons other than + # `young.prop`: a `network.size` too small for the seeded share, or a + # `target_pop` data frame with no members above the cessation age. The + # netstats object is self-consistent either way (the band's target statistics + # are zero), but no model with an `age.grp` term can be fitted on the + # population it describes, so say so here rather than leaving it to `netest`. + if (sex.cess.mod == TRUE && !any(attr_age.grp %in% top.age.grp)) { + warning("No node was placed in the post-cessation age group (age.grp ", + top.age.grp, ", ages ", age.breaks[top.age.grp], " and over). ", + "`nodefactor` and `nodematch` terms on `age.grp` will be missing ", + "that level, so fitting against these target statistics will fail ", + "on a length mismatch. Raise `network.size`, lower `young.prop`, ", + "or include post-cessation members in `target_pop`.") + } + # Common attr assignments (both paths) ----------------------------------- out$attr$age <- attr_age out$attr$sqrt.age <- attr_sqrt.age diff --git a/man/build_netstats.Rd b/man/build_netstats.Rd index f6c1e61..8bf008f 100644 --- a/man/build_netstats.Rd +++ b/man/build_netstats.Rd @@ -41,7 +41,11 @@ Hispanic, and White/Other).} \item{young.prop}{The proportion of the population that should be below the age of sexual cessation. Default is NULL (meaning no re-weighting of the \code{age.pyramid} parameter is performed). -This parameter is only used if the age of sexual cessation is less than the upper age bound.} +This parameter is only used if the age of sexual cessation is less than the upper age bound. +Must be strictly between 0 and 1: either endpoint empties one side of the split, and an +age group with no members cannot carry an \code{ergm} term. To start a model with effectively +nobody above the cessation age, use a value just below 1 (for example 0.999) rather than +1 itself, and let the post-cessation group fill over the burn-in.} \item{method}{Character. Either \code{"existing"} (default) or \code{"joint"}. \code{"existing"} reproduces the pre-refactor behavior byte-for-byte: target statistics for edges, nodefactor, and diff --git a/tests/testthat/test-sex-cessation.R b/tests/testthat/test-sex-cessation.R new file mode 100644 index 0000000..3b69fb4 --- /dev/null +++ b/tests/testthat/test-sex-cessation.R @@ -0,0 +1,182 @@ +# Tests for the sexual-cessation configuration of build_netstats(), where +# `age.limits[2]` sits above `age.sexual.cessation` so the population keeps +# aging past the point where it stops forming partnerships. +# +# Two defects motivated these: +# 1. `active.sex` was zeroed on `max(attr_age.grp)` as realized rather than on +# the age group the breaks declare. Whenever the post-cessation group came +# out empty, that marked the oldest *sexually active* group inactive and +# zeroed its degree, silently and without any error. +# 2. `young.prop = 1` empties the post-cessation group by construction, which +# triggered (1) and then failed much later in `netest()` on a target.stats +# length mismatch, because an age group with no members carries no `ergm` +# term. + +skip_without_artnetdata <- function() { + testthat::skip_if(system.file(package = "ARTnetData") == "", + "ARTnetData not installed") +} + +have_data <- system.file(package = "ARTnetData") != "" +if (have_data) { + # Sexual cessation at 65, population to 100: six age groups, the sixth being + # the post-cessation band. + set.seed(20260821L) + .ep_cess <- build_epistats(geog.lvl = "city", geog.cat = "Atlanta", + init.hiv.prev = c(0.33, 0.137, 0.084), + race = TRUE, time.unit = 7, + age.limits = c(15, 100), + age.sexual.cessation = 65) + .np_cess <- build_netparams(.ep_cess, smooth.main.dur = TRUE) + + # The default configuration, where the population limit and the cessation age + # coincide and there is no post-cessation band at all. + set.seed(20260821L) + .ep_plain <- build_epistats(geog.lvl = "city", geog.cat = "Atlanta", + init.hiv.prev = c(0.33, 0.137, 0.084), + race = TRUE, time.unit = 7) + .np_plain <- build_netparams(.ep_plain, smooth.main.dur = TRUE) +} + +mk_pop <- function(n, age_lo, age_hi) { + data.frame( + age = stats::runif(n, age_lo, age_hi), + race = sample(1:3, n, replace = TRUE, prob = c(0.5, 0.05, 0.45)), + deg.casl = sample(0:3, n, replace = TRUE), + deg.main = sample(0:2, n, replace = TRUE), + role.class = sample(0:2, n, replace = TRUE), + risk.grp = sample(1:5, n, replace = TRUE) + ) +} + + +# ---- The configuration itself ------------------------------------------------ + +test_that("sexual cessation adds a post-cessation age group and real mortality", { + skip_without_artnetdata() + expect_true(.ep_cess$sex.cess.mod) + expect_equal(.ep_cess$age.breaks, c(15, 25, 35, 45, 55, 65, 100)) + expect_equal(.ep_cess$age.grps, 6) + + set.seed(101L) + ns <- build_netstats(.ep_cess, .np_cess, expect.mort = 0.000478213, + network.size = 3000) + asmr <- ns$demog$asmr + rates <- as.matrix(asmr[asmr$age >= 65 & asmr$age <= 99, -1]) + # Under the default 15-to-65 configuration mortality is forced to 1 at 65. + # Here the cut has to move to the upper age limit instead. + expect_true(all(rates > 0 & rates < 1)) + expect_true(all(asmr[asmr$age == 100, -1] == 1)) + expect_equal(min(asmr$age[apply(asmr[, -1] >= 1, 1, any)]), 100) +}) + + +# ---- active.sex keys on the declared band, not the realized maximum ---------- + +test_that("active.sex is 0 exactly for nodes at or above the cessation age", { + skip_without_artnetdata() + set.seed(102L) + ns <- build_netstats(.ep_cess, .np_cess, expect.mort = 0.000478213, + network.size = 3000) + expect_gt(sum(ns$attr$age >= 65), 0) + expect_equal(ns$attr$active.sex == 0, ns$attr$age >= 65) + expect_equal(ns$attr$active.sex == 0, ns$attr$age.grp == 6) + # Inactive nodes cannot carry degree on any layer. + inactive <- ns$attr$active.sex == 0 + expect_true(all(ns$attr$deg.main[inactive] == 0)) + expect_true(all(ns$attr$deg.casl[inactive] == 0)) + expect_true(all(ns$attr$deg.tot[inactive] == 0)) +}) + +test_that("an empty post-cessation band does not retire the 55-64 group", { + skip_without_artnetdata() + # This is the regression test for defect (1). Every member of this target + # population is sexually active, so the post-cessation group is empty and the + # realized maximum age group is 5. Reading the band off the data marked all + # 55-to-64 year olds inactive and zeroed their degree. + set.seed(103L) + df <- mk_pop(2000, 20, 64) + expect_warning( + ns <- build_netstats(.ep_cess, .np_cess, expect.mort = 0.000478213, + target_pop = df), + "post-cessation age group" + ) + expect_equal(sum(ns$attr$age.grp == 6), 0) + expect_gt(sum(ns$attr$age.grp == 5), 0) + expect_true(all(ns$attr$active.sex == 1)) + expect_true(any(ns$attr$deg.main > 0)) +}) + +test_that("a populated band is unaffected by the fix", { + skip_without_artnetdata() + set.seed(104L) + df <- mk_pop(2000, 20, 90) + ns <- build_netstats(.ep_cess, .np_cess, expect.mort = 0.000478213, + target_pop = df) + expect_gt(sum(ns$attr$age.grp == 6), 0) + expect_equal(ns$attr$active.sex == 0, ns$attr$age >= 65) +}) + + +# ---- young.prop endpoints ----------------------------------------------------- + +test_that("young.prop = 1 is rejected with an actionable message", { + skip_without_artnetdata() + expect_error( + build_netstats(.ep_cess, .np_cess, expect.mort = 0.000478213, + network.size = 3000, young.prop = 1), + "must be less than 1" + ) + expect_error( + build_netstats(.ep_cess, .np_cess, expect.mort = 0.000478213, + network.size = 3000, young.prop = 1), + "0.999" + ) +}) + +test_that("young.prop = 0 is rejected", { + skip_without_artnetdata() + expect_error( + build_netstats(.ep_cess, .np_cess, expect.mort = 0.000478213, + network.size = 3000, young.prop = 0), + "must be greater than 0" + ) +}) + +test_that("young.prop just below 1 seeds the band and nothing else changes", { + skip_without_artnetdata() + set.seed(105L) + ns <- build_netstats(.ep_cess, .np_cess, expect.mort = 0.000478213, + network.size = 25000, young.prop = 0.999) + n_retired <- sum(ns$attr$age.grp == 6) + expect_gt(n_retired, 0) + expect_lt(n_retired / ns$demog$num, 0.01) + expect_equal(ns$attr$active.sex == 0, ns$attr$age >= 65) + # The band's target statistics are zero on every layer, which is what makes + # `ergm` pin those terms off and keep retired nodes tie-free. + for (layer in c("main", "casl", "inst")) { + nf <- ns[[layer]]$nodefactor_age.grp + expect_length(nf, 6) + expect_equal(unname(nf[6]), 0) + expect_gt(sum(nf[1:5]), 0) + } +}) + + +# ---- The default configuration is untouched ---------------------------------- + +test_that("young.prop is inert without a sexual-cessation band", { + skip_without_artnetdata() + # With `age.limits[2]` equal to the cessation age there is no split to + # reweight, so no value of young.prop should be rejected and none should + # change the result. + set.seed(106L) + a <- build_netstats(.ep_plain, .np_plain, expect.mort = 0.000478213, + network.size = 3000, young.prop = 1) + set.seed(106L) + b <- build_netstats(.ep_plain, .np_plain, expect.mort = 0.000478213, + network.size = 3000) + expect_equal(a$attr, b$attr) + expect_true(all(a$attr$active.sex == 1)) + expect_equal(length(.ep_plain$age.breaks) - 1L, 5L) +})