From 0ac844f7e4a6de7019a7c5922ad0704b9df0d92c Mon Sep 17 00:00:00 2001 From: Samuel Jenness Date: Sat, 7 Nov 2026 09:22:46 -0500 Subject: [PATCH 1/4] Balance race mixing under method = "joint", check target feasibility, seed the partner-race imputation Under method = "joint", nodematch_race was standardized on the ego side only, so the partner side was still ARTnet's and the race targets came out infeasible in mostly-White cities. Partner race is now modeled with a multinomial on ego attributes, the reported mixing matrix is built at the synthetic population and symmetrically scaled with nodefactor_race held fixed, and nodematch_race and nodematch_race_diffF are read from the balanced matrix (#9). The opt-in nodematch.race.floor covers the legacy method (#9). build_netstats() now stops when a layer's race targets cannot be realized by any network, and check_race_feasibility() is exported (#10). build_netparams() gains a seed argument for the partner-race imputation (#11). --- DESCRIPTION | 5 +- NAMESPACE | 32 ++- R/NetParams.R | 81 +++++- R/NetStats.R | 84 ++++-- R/RaceMixing.R | 264 +++++++++++++++++++ inst/validation/validate_backward_compat.R | 2 +- man/build_netparams.Rd | 16 +- man/build_netstats.Rd | 30 ++- man/check_race_feasibility.Rd | 55 ++++ tests/testthat/test-joint-dyad.R | 5 +- tests/testthat/test-race-mixing.R | 291 +++++++++++++++++++++ 11 files changed, 823 insertions(+), 42 deletions(-) create mode 100644 R/RaceMixing.R create mode 100644 man/check_race_feasibility.Rd create mode 100644 tests/testthat/test-race-mixing.R diff --git a/DESCRIPTION b/DESCRIPTION index 13c5e29..0e40118 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,7 +13,8 @@ Depends: R (>= 4.1), EpiModel (>= 2.5) Imports: - dplyr + dplyr, + nnet Suggests: ARTnetData, knitr, @@ -26,4 +27,4 @@ Remotes: github::EpiModel/ARTnetData@main Roxygen: list(markdown = TRUE) Config/testthat/edition: 3 -Config/roxygen2/version: 8.0.0 +Config/roxygen2/version: 8.1.1 diff --git a/NAMESPACE b/NAMESPACE index 8edaa0d..f78b6fe 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ export(build_epistats) export(build_netparams) export(build_netstats) +export(check_race_feasibility) export(make_race_combo) export(reweight_age_pyr) export(trim_epistats) @@ -11,18 +12,21 @@ export(update_asmr) import(EpiModel) import(dplyr) import(utils) -importFrom(stats,AIC) -importFrom(stats,binomial) -importFrom(stats,coef) -importFrom(stats,formula) -importFrom(stats,gaussian) -importFrom(stats,glm) -importFrom(stats,lm) -importFrom(stats,median) -importFrom(stats,poisson) -importFrom(stats,predict) -importFrom(stats,rbinom) -importFrom(stats,reformulate) -importFrom(stats,runif) -importFrom(stats,update) +importFrom(nnet,multinom) +importFrom(stats, + AIC, + binomial, + coef, + formula, + gaussian, + glm, + lm, + median, + poisson, + predict, + rbinom, + reformulate, + runif, + update +) importFrom(utils,head) diff --git a/R/NetParams.R b/R/NetParams.R index df507f5..322b123 100644 --- a/R/NetParams.R +++ b/R/NetParams.R @@ -117,6 +117,48 @@ fit_joint_glm <- function(d, response, main_terms, } +# Joint multinomial (baseline-category logit) for partner race on ego +# attributes: the partner-race analog of the binomial same-race fit from +# fit_joint_glm (issue #9). Same term assembly and the same AIC-selected +# age.grp:race interaction. `response` must name a factor column whose +# levels are the race categories in order; the fit is consumed by +# build_netstats(method = "joint") to form the reported race mixing matrix +# at the synthetic population before balancing it. +#' @importFrom nnet multinom +fit_joint_multinom <- function(d, response, main_terms, race, geog.lvl) { + if (isTRUE(race)) { + main_terms <- c(main_terms, "as.factor(race.cat.num)") + } + if (!is.null(geog.lvl)) { + main_terms <- c("geogYN", main_terms) + } + fit_one <- function(terms) { + fml <- reformulate(terms, response = response) + tryCatch( + suppressWarnings(nnet::multinom(fml, data = d, trace = FALSE, maxit = 500)), + error = function(e) NULL + ) + } + best <- fit_one(main_terms) + if (is.null(best)) { + stop("The joint multinomial partner-race model for ", response, + " failed to fit.", call. = FALSE) + } + selected <- character(0) + candidates <- if (isTRUE(race)) "age.grp:as.factor(race.cat.num)" else character(0) + for (ix in candidates) { + cand <- fit_one(c(main_terms, selected, ix)) + if (!is.null(cand) && cand$convergence == 0 && AIC(cand) < AIC(best)) { + best <- cand + selected <- c(selected, ix) + } + } + attr(best, "selected_interactions") <- selected + attr(best, "candidate_interactions") <- candidates + best +} + + #' Calculate Individual-Level Network Parameters #' #' @description Builds statistical models predicting mean degree, mixing, and duration of sexual @@ -169,12 +211,20 @@ fit_joint_glm <- function(d, response, main_terms, #' \item \strong{Dyad-level} (partnership data with ego attrs on RHS): #' binomial for same-age-group at `$$joint_nm_age_model`, binomial #' for same-race (when `race = TRUE`) at `$$joint_nm_race_model`, -#' and Gaussian fits for the age absolute-difference terms at -#' `$$joint_absdiff_age_model` and +#' a multinomial (`nnet::multinom`) for partner race (when `race = TRUE`) at +#' `$$joint_prace_model`, and Gaussian fits for the age +#' absolute-difference terms at `$$joint_absdiff_age_model` and #' `$$joint_absdiff_sqrtage_model`. #' } #' These models are consumed by [`build_netstats`] under `method = "joint"` to produce -#' internally-consistent g-computation target statistics. +#' internally-consistent g-computation target statistics. The partner-race multinomial +#' is what lets `build_netstats` balance race mixing to the target population instead of +#' standardizing the ego side only. +#' @param seed Optional integer. If supplied, `set.seed(seed)` is called at the start of the +#' function. The one stochastic step is the imputation of missing partner race, which +#' draws from each ego race's observed partner-race distribution (pooled across +#' partnership types), so two unseeded calls give slightly different race mixing +#' parameters. Defaults to `NULL` (no seeding), which preserves prior behavior. #' @param browser If `TRUE`, run `build_netparams` in interactive browser mode. #' #' @details @@ -219,10 +269,14 @@ build_netparams <- function(epistats, oo.nquants = 5, duration.method = c("empirical", "joint_lm"), method = c("existing", "joint"), + seed = NULL, browser = FALSE) { method <- match.arg(method) duration.method <- match.arg(duration.method) + # Optional reproducibility of the partner-race imputation (issue #11). + if (!is.null(seed)) set.seed(seed) + # Ensures that ARTnetData is installed if (system.file(package = "ARTnetData") == "") stop(missing_data_msg) @@ -800,6 +854,15 @@ build_netparams <- function(epistats, race = race, geog.lvl = geog.lvl, family = binomial() ) + # Partner race itself (not just sameness), for the balanced mixing + # matrix in build_netstats (#9). + lmain$p_race_f <- factor(lmain$p_race.cat.num, levels = race.categories) + out$main$joint_prace_model <- fit_joint_multinom( + lmain, + response = "p_race_f", + main_terms = c("age.grp", "sqrt(age.grp)", "hiv2"), + race = race, geog.lvl = geog.lvl + ) } out$main$joint_absdiff_age_model <- fit_joint_glm( lmain, @@ -1211,6 +1274,12 @@ build_netparams <- function(epistats, race = race, geog.lvl = geog.lvl, family = binomial() ) + lcasl$p_race_f <- factor(lcasl$p_race.cat.num, levels = race.categories) + out$casl$joint_prace_model <- fit_joint_multinom( + lcasl, response = "p_race_f", + main_terms = c("age.grp", "sqrt(age.grp)", "hiv2"), + race = race, geog.lvl = geog.lvl + ) } out$casl$joint_absdiff_age_model <- fit_joint_glm( lcasl, response = "ad", @@ -1526,6 +1595,12 @@ build_netparams <- function(epistats, race = race, geog.lvl = geog.lvl, family = binomial() ) + linst$p_race_f <- factor(linst$p_race.cat.num, levels = race.categories) + out$inst$joint_prace_model <- fit_joint_multinom( + linst, response = "p_race_f", + main_terms = c("age.grp", "sqrt(age.grp)", "hiv2"), + race = race, geog.lvl = geog.lvl + ) } out$inst$joint_absdiff_age_model <- fit_joint_glm( linst, response = "ad", diff --git a/R/NetStats.R b/R/NetStats.R index aa4cdb2..85b6748 100644 --- a/R/NetStats.R +++ b/R/NetStats.R @@ -250,12 +250,36 @@ #' partner-race uncertainty via `joint_nm_race_model`, then median-aggregated within #' stratum. `nodefactor_risk.grp` and `diss.homog` still use the within-ARTnet #' univariate / aggregated values — those are not consumed by the standard -#' EpiModelHIV-Template dissolution offset. +#' EpiModelHIV-Template dissolution offset. The race mixing targets are the exception +#' to the per-ego aggregation: `nodematch_race` and `nodematch_race_diffF` come from a +#' race mixing matrix built at the synthetic population from the partner-race +#' multinomial (`netparams$$joint_prace_model`) and balanced so that it describes +#' a closed population, with `nodefactor_race` held fixed. Standardizing the ego side +#' alone leaves the partner side as ARTnet's, and in a city whose race composition is far +#' from ARTnet's the two sides do not add up (Seattle: White egos report more Black +#' partners than Black egos have partnership endpoints). The balance is a symmetric +#' scaling of the averaged reports that keeps their homophily odds ratios and refits the +#' margins, so the targets are always feasible. The matrices are stored as +#' `mixmat_race` (balanced expected edge counts) and `mixmat_race_reported` (ego-side +#' reported endpoints, rows are the ego race). #' @param seed Optional integer. If supplied, `set.seed(seed)` is called at the start of the #' function for reproducibility of the stochastic sampling path (default `target_pop` and #' list forms). Defaults to `NULL` (no seeding), which preserves prior behavior. The #' data.frame form of `target_pop` bypasses attribute sampling and is deterministic #' regardless of this argument. +#' @param nodematch.race.floor Logical. If `TRUE`, each layer's `nodematch_race_diffF` is raised +#' to the random-mixing value `edges * sum(p^2)`, with `p` the endpoint shares +#' `nodefactor_race / sum(nodefactor_race)`, wherever it falls below, with a message. Any +#' network with the layer's `nodefactor_race` targets can realize random mixing, so the +#' floored targets are always feasible. This is a stopgap for `method = "existing"`, whose +#' race mixing is not balanced to the target population; under `method = "joint"` the +#' balanced targets are feasible already and the floor should not bind. The per-race +#' `nodematch_race` vector is left unchanged. Default `FALSE`. +#' @param check.feasibility Logical. If `TRUE` (default), the race targets are checked with +#' [`check_race_feasibility`] at the end and the function stops when any layer has a +#' race whose `nodefactor_race` exceeds `edges + nodematch_race_diffF`, naming the layer, +#' the race and the excess. Such targets cannot be fit by MCMLE. Set to `FALSE` to return +#' the targets anyway. #' @param browser If `TRUE`, run `build_netparams` in interactive browser mode. #' #' @details @@ -321,6 +345,8 @@ build_netstats <- function(epistats, netparams, method = c("existing", "joint"), target_pop = NULL, seed = NULL, + nodematch.race.floor = FALSE, + check.feasibility = TRUE, browser = FALSE) { method <- match.arg(method) @@ -925,12 +951,17 @@ build_netstats <- function(epistats, netparams, n_race <- length(netparams$main$nf.race) out$main$nodefactor_race <- vapply(seq_len(n_race), function(r) sum(pred_deg_main[out$attr$race == r]), numeric(1)) - # nodematch_race[r] = expected count of race-r-to-race-r edges - race_match_contrib <- edge_wt_main * dyad_main$pred_nm_race - out$main$nodematch_race <- vapply(seq_len(n_race), - function(r) sum(race_match_contrib[out$attr$race == r]) / 2, - numeric(1)) - out$main$nodematch_race_diffF <- sum(race_match_contrib) / 2 + # Race mixing from both sides (#9): the reported endpoint matrix from + # the partner-race multinomial, balanced to a closed population with + # nodefactor_race held fixed. nodematch_race[r] is then the expected + # count of race-r-to-race-r edges in the balanced matrix. + race_mix <- .joint_race_mixing( + netparams$main$joint_prace_model, synth, pred_deg_main, + out$attr$race, n_race, flattened_race_level, layer = "main") + out$main$nodematch_race <- race_mix$nodematch_race + out$main$nodematch_race_diffF <- race_mix$nodematch_race_diffF + out$main$mixmat_race <- race_mix$mixmat_race + out$main$mixmat_race_reported <- race_mix$mixmat_race_reported } # nodefactor("age.grp") + nodematch("age.grp") -------------------------- @@ -1046,11 +1077,13 @@ build_netstats <- function(epistats, netparams, n_race <- length(netparams$casl$nf.race) out$casl$nodefactor_race <- vapply(seq_len(n_race), function(r) sum(pred_deg_casl[out$attr$race == r]), numeric(1)) - race_match_contrib <- edge_wt_casl * dyad_casl$pred_nm_race - out$casl$nodematch_race <- vapply(seq_len(n_race), - function(r) sum(race_match_contrib[out$attr$race == r]) / 2, - numeric(1)) - out$casl$nodematch_race_diffF <- sum(race_match_contrib) / 2 + race_mix <- .joint_race_mixing( + netparams$casl$joint_prace_model, synth, pred_deg_casl, + out$attr$race, n_race, flattened_race_level, layer = "casl") + out$casl$nodematch_race <- race_mix$nodematch_race + out$casl$nodematch_race_diffF <- race_mix$nodematch_race_diffF + out$casl$mixmat_race <- race_mix$mixmat_race + out$casl$mixmat_race_reported <- race_mix$mixmat_race_reported } n_age <- length(netparams$casl$nf.age.grp) @@ -1160,11 +1193,13 @@ build_netstats <- function(epistats, netparams, n_race <- length(netparams$inst$nf.race) out$inst$nodefactor_race <- vapply(seq_len(n_race), function(r) sum(pred_deg_inst[out$attr$race == r]), numeric(1)) - race_match_contrib <- edge_wt_inst * dyad_inst$pred_nm_race - out$inst$nodematch_race <- vapply(seq_len(n_race), - function(r) sum(race_match_contrib[out$attr$race == r]) / 2, - numeric(1)) - out$inst$nodematch_race_diffF <- sum(race_match_contrib) / 2 + race_mix <- .joint_race_mixing( + netparams$inst$joint_prace_model, synth, pred_deg_inst, + out$attr$race, n_race, flattened_race_level, layer = "inst") + out$inst$nodematch_race <- race_mix$nodematch_race + out$inst$nodematch_race_diffF <- race_mix$nodematch_race_diffF + out$inst$mixmat_race <- race_mix$mixmat_race + out$inst$mixmat_race_reported <- race_mix$mixmat_race_reported } n_age <- length(netparams$inst$nf.age.grp) @@ -1191,6 +1226,21 @@ build_netstats <- function(epistats, netparams, function(h) sum(pred_deg_inst[out$attr$diag.status == h]), numeric(1)) } + # Race targets: the opt-in random-mixing floor (#9) and the feasibility + # check (#10), both methods. A race holding more partnership endpoints than + # edges + nodematch_race_diffF allows cannot be fit by MCMLE; better to + # stop here than on the cluster. + if (isTRUE(race)) { + if (isTRUE(nodematch.race.floor)) { + for (layer in c("main", "casl", "inst")) { + out[[layer]] <- .apply_race_floor(out[[layer]], layer) + } + } + if (isTRUE(check.feasibility)) { + .stop_if_race_infeasible(out, labels = flattened_race_level) + } + } + # The dissolution formulas were created inside this function, so their # environment holds `epistats` and, through it, the ARTnet data. Point them # at the global environment so a saved netstats object carries no data. diff --git a/R/RaceMixing.R b/R/RaceMixing.R new file mode 100644 index 0000000..ccc47cf --- /dev/null +++ b/R/RaceMixing.R @@ -0,0 +1,264 @@ +# Race-mixing helpers for build_netstats(): the balanced partner-race mixing +# matrix under method = "joint" (issue #9), the opt-in random-mixing floor +# (#9), and the feasibility check on the race targets (#10). +# +# Race levels are integer-coded 1..K in the order of epistats$race.level +# (default 1 = Black, 2 = Hispanic, 3 = White/Other). + + +# Predicted partner-race probabilities from an nnet::multinom fit, always +# returned as an n x K matrix with columns in race-level order. nnet returns +# a bare vector of P(level 2) when there are two levels, and a bare vector of +# length K when newdata has a single row. +.predict_prace <- function(model, newdata, n_levels) { + p <- predict(model, newdata = newdata, type = "probs") + if (is.null(dim(p))) { + p <- if (n_levels == 2) cbind(1 - p, p) else matrix(p, nrow = 1) + } + p <- unname(as.matrix(p)) + if (ncol(p) != n_levels || nrow(p) != nrow(newdata)) { + stop("Partner-race prediction has dimension ", nrow(p), " x ", ncol(p), + "; expected ", nrow(newdata), " x ", n_levels, ".", call. = FALSE) + } + p +} + + +# Reported endpoint matrix at the synthetic population. A[r, s] is the +# expected number of partnership endpoints held by egos of race r whose +# partner is of race s: sum over egos of race r of predicted degree times +# the predicted probability that a partner is of race s. Rows are the ego +# side and sum to nodefactor_race by construction. In a closed population +# A would be symmetric; egocentric reports standardized on the ego side +# only are not. +.reported_race_mixing <- function(pred_deg, probs, race, n_levels) { + A <- matrix(0, n_levels, n_levels) + for (r in seq_len(n_levels)) { + sel <- race == r + if (any(sel)) { + A[r, ] <- colSums(pred_deg[sel] * probs[sel, , drop = FALSE]) + } + } + A +} + + +# Balance a reported endpoint matrix so that it describes a closed +# population, holding each race's total endpoints (nodefactor_race) fixed. +# +# With S = (A + t(A)) / 2, find the diagonal scaling d such that +# B = diag(d) S diag(d) has rowSums(B) == rowSums(A). B is symmetric, so it +# is the endpoint matrix of a real network: B[r, r] / 2 same-race edges and +# B[r, s] cross-race edges. The scaling keeps the homophily odds ratios of +# S (B[r, r] B[s, s] / B[r, s]^2 == S[r, r] S[s, s] / S[r, s]^2) and refits +# the margins to the target population, so it is the log-linear transport +# of ARTnet's selective mixing to a population with a different race +# composition. Plain averaging of the two reports of each cross-race cell +# is not used because, with the margins held, it can drive a minority +# group's same-race cell negative when the majority's reports of that +# group exceed what the group's own degree can supply. +# +# Solved by the symmetric Sinkhorn iteration d <- d * sqrt(m / (S d)). +# A race with no endpoints (m == 0) gets d == 0 and an all-zero row and +# column. Returns B with the scaling factors and the iteration count as +# attributes. +.balance_race_mixing <- function(A, tol = 1e-10, maxit = 5000L) { + if (!is.matrix(A) || nrow(A) != ncol(A)) { + stop("A must be a square matrix.", call. = FALSE) + } + if (any(A < 0)) stop("A must be non-negative.", call. = FALSE) + m <- rowSums(A) + S <- (A + t(A)) / 2 + active <- m > 0 + d <- as.numeric(active) + it <- 0L + converged <- !any(active) + while (!converged && it < maxit) { + it <- it + 1L + Sd <- as.vector(S %*% d) + if (any(Sd[active] <= 0)) { + stop("Race mixing matrix cannot be balanced: a race with partnership ", + "endpoints has no partners of any race present in the population.", + call. = FALSE) + } + # rowSums(B) = d * (S d); ratio is the target margin over the current one + ratio <- rep(1, length(m)) + ratio[active] <- m[active] / (d[active] * Sd[active]) + d[active] <- d[active] * sqrt(ratio[active]) + converged <- max(abs(ratio - 1)) < tol + } + if (!converged) { + stop("Race mixing balance did not converge in ", maxit, " iterations.", + call. = FALSE) + } + B <- outer(d, d) * S + dimnames(B) <- dimnames(A) + attr(B, "scaling") <- d + attr(B, "iterations") <- it + B +} + + +# Full joint-method race-mixing step for one layer. Predicts partner race +# for every synthetic node, forms the reported endpoint matrix, balances it, +# and returns the nodematch targets plus the two matrices that +# build_netstats stores as additive fields under method = "joint". +.joint_race_mixing <- function(prace_model, synth, pred_deg, race, n_levels, + race_labels = NULL, layer = "") { + if (is.null(prace_model)) { + stop("method = 'joint' needs the partner-race model ", + "netparams$", layer, "$joint_prace_model, which build_netparams() ", + "fits under method = 'joint' from ARTnet 2.10.0. Re-run ", + "build_netparams(..., method = 'joint').", call. = FALSE) + } + probs <- .predict_prace(prace_model, synth, n_levels) + A <- .reported_race_mixing(pred_deg, probs, race, n_levels) + if (!is.null(race_labels) && length(race_labels) == n_levels) { + dimnames(A) <- list(race_labels, race_labels) + } + B <- .balance_race_mixing(A) + nodematch_race <- diag(B) / 2 + mixmat <- B + diag(mixmat) <- diag(B) / 2 + attr(mixmat, "scaling") <- NULL + attr(mixmat, "iterations") <- NULL + list( + nodematch_race = unname(nodematch_race), + nodematch_race_diffF = sum(nodematch_race), + mixmat_race = mixmat, + mixmat_race_reported = A + ) +} + + +# Opt-in random-mixing floor on nodematch_race_diffF for one layer (#9). +# With endpoint shares p = nodefactor_race / sum(nodefactor_race), random +# mixing gives edges * sum(p^2) same-race edges, which every network with +# those nodefactor targets can realize (the feasibility condition +# nodefactor[r] <= edges + nodematch at that value reduces to +# (1 - p[r])^2 >= 0). The floor raises the target to that value when it +# falls below and says so. Only nodematch_race_diffF is touched, the +# target that the EpiModelHIV template formulas consume. +.apply_race_floor <- function(layer_stats, layer_name) { + E <- layer_stats$edges + nf <- layer_stats$nodefactor_race + M <- layer_stats$nodematch_race_diffF + if (is.null(E) || is.null(nf) || is.null(M)) return(layer_stats) + p <- as.numeric(nf) / sum(nf) + floor_val <- E * sum(p^2) + if (M < floor_val) { + message(sprintf( + paste0("nodematch.race.floor: raising the %s same-race target from ", + "%.1f to %.1f edges (share %.3f -> %.3f), the random-mixing floor."), + layer_name, M, floor_val, M / E, floor_val / E)) + layer_stats$nodematch_race_diffF <- floor_val + } + layer_stats +} + + +#' Check Feasibility of Race-Mixing Target Statistics +#' +#' @description Checks, for each network layer, whether the race target +#' statistics in a `netstats` object can be realized by any network. +#' +#' @param netstats Output from [`build_netstats`]. +#' @param tol Numeric slack above zero that is still treated as feasible, +#' to absorb floating-point noise. +#' @param labels Optional character vector of race labels in race-code order, +#' for the `label` column. When `NULL`, the labels are taken from the +#' dimnames of `mixmat_race` when present (`method = "joint"`), and are the +#' integer codes otherwise. +#' +#' @details +#' With a single same-race count `M` (the `nodematch("race", diff = FALSE)` +#' target, `nodematch_race_diffF`) and `edges` `E`, every race `r` needs +#' `nodefactor_race[r] <= E + M`: a race cannot hold more partnership +#' endpoints than the same-race edges can contribute twice plus the +#' cross-race edges can contribute once. The function reports +#' `slack = nodefactor_race - edges - nodematch_race_diffF` for every layer +#' and race; a positive slack means no network satisfies the targets. +#' +#' Infeasible targets fail in `ergm` MCMLE estimation with "The MPLE does +#' not exist!" and then the simulated-density guard, usually after a +#' warning that `nodefactor.race.` is linearly dependent on +#' `edges + nodematch.race`, which is this same boundary seen from the +#' other side. Stochastic-Approximation estimation settles on the nearest +#' reachable network without complaint, so the problem can go unnoticed +#' until a model is refit with MCMLE. +#' +#' `build_netstats()` runs this check by default (`check.feasibility = TRUE`) +#' and stops on a positive slack. The condition covers the +#' `nodematch("race", diff = FALSE)` specification used by the EpiModelHIV +#' template formulas; the per-race `nodematch_race` vector is not checked. +#' +#' @return A data frame with one row per layer and race: `layer`, `race` +#' (integer code), `label`, `nodefactor`, `edges`, `nodematch_race_diffF`, +#' `slack`, and `feasible`. Zero rows when `netstats` has no race targets. +#' +#' @export +#' +#' @examples +#' netstats <- readRDS(system.file("netstats-example.rds", package = "ARTnet")) +#' check_race_feasibility(netstats) +#' +check_race_feasibility <- function(netstats, tol = 1e-6, labels = NULL) { + rows <- list() + for (l in c("main", "casl", "inst")) { + x <- netstats[[l]] + if (is.null(x$edges) || is.null(x$nodefactor_race) || + is.null(x$nodematch_race_diffF)) { + next + } + nf <- as.numeric(x$nodefactor_race) + # demog$props columns follow ARTnetData::race.dist, not the race codes, + # so they are not used for labels. + lab <- labels + if (is.null(lab) && !is.null(x$mixmat_race)) lab <- rownames(x$mixmat_race) + if (is.null(lab) || length(lab) != length(nf)) lab <- as.character(seq_along(nf)) + slack <- nf - x$edges - x$nodematch_race_diffF + rows[[l]] <- data.frame( + layer = l, + race = seq_along(nf), + label = lab, + nodefactor = nf, + edges = x$edges, + nodematch_race_diffF = x$nodematch_race_diffF, + slack = slack, + feasible = slack <= tol, + stringsAsFactors = FALSE + ) + } + if (length(rows) == 0) { + message("No race targets in netstats (race = FALSE); nothing to check.") + return(data.frame(layer = character(0), race = integer(0), + label = character(0), nodefactor = numeric(0), + edges = numeric(0), nodematch_race_diffF = numeric(0), + slack = numeric(0), feasible = logical(0), + stringsAsFactors = FALSE)) + } + out <- do.call(rbind, rows) + rownames(out) <- NULL + out +} + + +# Stop with a layer-by-race account of the excess when the race targets +# are infeasible. Called at the end of build_netstats(). +.stop_if_race_infeasible <- function(netstats, labels = NULL) { + feas <- check_race_feasibility(netstats, labels = labels) + bad <- feas[!feas$feasible, , drop = FALSE] + if (nrow(bad) == 0) return(invisible(TRUE)) + lines <- sprintf( + paste0(" %s: nodefactor_race[%d] (%s) exceeds edges + nodematch_race_diffF ", + "by %.1f endpoints (%.1f%% of nodefactor)"), + bad$layer, bad$race, bad$label, bad$slack, 100 * bad$slack / bad$nodefactor) + stop("Race-mixing targets are infeasible: no network can carry these ", + "nodefactor_race and nodematch_race_diffF values together, and MCMLE ", + "estimation will fail on them.\n", + paste(lines, collapse = "\n"), "\n", + "Options: method = 'joint' (partner-race mixing balanced to the target ", + "population), nodematch.race.floor = TRUE (raise the same-race target ", + "to random mixing), or check.feasibility = FALSE to return the targets ", + "anyway. See ?check_race_feasibility.", call. = FALSE) +} diff --git a/inst/validation/validate_backward_compat.R b/inst/validation/validate_backward_compat.R index baea81e..183e0a8 100644 --- a/inst/validation/validate_backward_compat.R +++ b/inst/validation/validate_backward_compat.R @@ -86,7 +86,7 @@ PARAM_SETS <- list( # fields are added (e.g., $joint_model). .strip_additive <- function(netparams) { additive <- c("joint_model", "joint_concurrent_model", - "joint_nm_age_model", "joint_nm_race_model", + "joint_nm_age_model", "joint_nm_race_model", "joint_prace_model", "joint_absdiff_age_model", "joint_absdiff_sqrtage_model") for (layer in c("main", "casl", "inst", "all")) { if (is.null(netparams[[layer]])) next diff --git a/man/build_netparams.Rd b/man/build_netparams.Rd index 14cc261..bad1c95 100644 --- a/man/build_netparams.Rd +++ b/man/build_netparams.Rd @@ -10,6 +10,7 @@ build_netparams( oo.nquants = 5, duration.method = c("empirical", "joint_lm"), method = c("existing", "joint"), + seed = NULL, browser = FALSE ) } @@ -64,12 +65,21 @@ concurrency indicator at \verb{$$joint_concurrent_model}. \item \strong{Dyad-level} (partnership data with ego attrs on RHS): binomial for same-age-group at \verb{$$joint_nm_age_model}, binomial for same-race (when \code{race = TRUE}) at \verb{$$joint_nm_race_model}, -and Gaussian fits for the age absolute-difference terms at -\verb{$$joint_absdiff_age_model} and +a multinomial (\code{nnet::multinom}) for partner race (when \code{race = TRUE}) at +\verb{$$joint_prace_model}, and Gaussian fits for the age +absolute-difference terms at \verb{$$joint_absdiff_age_model} and \verb{$$joint_absdiff_sqrtage_model}. } These models are consumed by \code{\link{build_netstats}} under \code{method = "joint"} to produce -internally-consistent g-computation target statistics.} +internally-consistent g-computation target statistics. The partner-race multinomial +is what lets \code{build_netstats} balance race mixing to the target population instead of +standardizing the ego side only.} + +\item{seed}{Optional integer. If supplied, \code{set.seed(seed)} is called at the start of the +function. The one stochastic step is the imputation of missing partner race, which +draws from each ego race's observed partner-race distribution (pooled across +partnership types), so two unseeded calls give slightly different race mixing +parameters. Defaults to \code{NULL} (no seeding), which preserves prior behavior.} \item{browser}{If \code{TRUE}, run \code{build_netparams} in interactive browser mode.} } diff --git a/man/build_netstats.Rd b/man/build_netstats.Rd index 8bf008f..e197be8 100644 --- a/man/build_netstats.Rd +++ b/man/build_netstats.Rd @@ -16,6 +16,8 @@ build_netstats( method = c("existing", "joint"), target_pop = NULL, seed = NULL, + nodematch.race.floor = FALSE, + check.feasibility = TRUE, browser = FALSE ) } @@ -68,7 +70,18 @@ per-ego predicted log(duration) from \code{joint_duration_model}, marginalized o partner-race uncertainty via \code{joint_nm_race_model}, then median-aggregated within stratum. \code{nodefactor_risk.grp} and \code{diss.homog} still use the within-ARTnet univariate / aggregated values — those are not consumed by the standard -EpiModelHIV-Template dissolution offset.} +EpiModelHIV-Template dissolution offset. The race mixing targets are the exception +to the per-ego aggregation: \code{nodematch_race} and \code{nodematch_race_diffF} come from a +race mixing matrix built at the synthetic population from the partner-race +multinomial (\verb{netparams$$joint_prace_model}) and balanced so that it describes +a closed population, with \code{nodefactor_race} held fixed. Standardizing the ego side +alone leaves the partner side as ARTnet's, and in a city whose race composition is far +from ARTnet's the two sides do not add up (Seattle: White egos report more Black +partners than Black egos have partnership endpoints). The balance is a symmetric +scaling of the averaged reports that keeps their homophily odds ratios and refits the +margins, so the targets are always feasible. The matrices are stored as +\code{mixmat_race} (balanced expected edge counts) and \code{mixmat_race_reported} (ego-side +reported endpoints, rows are the ego race).} \item{target_pop}{Optional specification of the synthetic target population. Defaults to NULL, which uses the legacy patchwork of reference sources (NCHS age pyramid + @@ -111,6 +124,21 @@ list forms). Defaults to \code{NULL} (no seeding), which preserves prior behavio data.frame form of \code{target_pop} bypasses attribute sampling and is deterministic regardless of this argument.} +\item{nodematch.race.floor}{Logical. If \code{TRUE}, each layer's \code{nodematch_race_diffF} is raised +to the random-mixing value \code{edges * sum(p^2)}, with \code{p} the endpoint shares +\code{nodefactor_race / sum(nodefactor_race)}, wherever it falls below, with a message. Any +network with the layer's \code{nodefactor_race} targets can realize random mixing, so the +floored targets are always feasible. This is a stopgap for \code{method = "existing"}, whose +race mixing is not balanced to the target population; under \code{method = "joint"} the +balanced targets are feasible already and the floor should not bind. The per-race +\code{nodematch_race} vector is left unchanged. Default \code{FALSE}.} + +\item{check.feasibility}{Logical. If \code{TRUE} (default), the race targets are checked with +\code{\link{check_race_feasibility}} at the end and the function stops when any layer has a +race whose \code{nodefactor_race} exceeds \code{edges + nodematch_race_diffF}, naming the layer, +the race and the excess. Such targets cannot be fit by MCMLE. Set to \code{FALSE} to return +the targets anyway.} + \item{browser}{If \code{TRUE}, run \code{build_netparams} in interactive browser mode.} } \description{ diff --git a/man/check_race_feasibility.Rd b/man/check_race_feasibility.Rd new file mode 100644 index 0000000..dabd251 --- /dev/null +++ b/man/check_race_feasibility.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RaceMixing.R +\name{check_race_feasibility} +\alias{check_race_feasibility} +\title{Check Feasibility of Race-Mixing Target Statistics} +\usage{ +check_race_feasibility(netstats, tol = 1e-06, labels = NULL) +} +\arguments{ +\item{netstats}{Output from \code{\link{build_netstats}}.} + +\item{tol}{Numeric slack above zero that is still treated as feasible, +to absorb floating-point noise.} + +\item{labels}{Optional character vector of race labels in race-code order, +for the \code{label} column. When \code{NULL}, the labels are taken from the +dimnames of \code{mixmat_race} when present (\code{method = "joint"}), and are the +integer codes otherwise.} +} +\value{ +A data frame with one row per layer and race: \code{layer}, \code{race} +(integer code), \code{label}, \code{nodefactor}, \code{edges}, \code{nodematch_race_diffF}, +\code{slack}, and \code{feasible}. Zero rows when \code{netstats} has no race targets. +} +\description{ +Checks, for each network layer, whether the race target +statistics in a \code{netstats} object can be realized by any network. +} +\details{ +With a single same-race count \code{M} (the \code{nodematch("race", diff = FALSE)} +target, \code{nodematch_race_diffF}) and \code{edges} \code{E}, every race \code{r} needs +\code{nodefactor_race[r] <= E + M}: a race cannot hold more partnership +endpoints than the same-race edges can contribute twice plus the +cross-race edges can contribute once. The function reports +\code{slack = nodefactor_race - edges - nodematch_race_diffF} for every layer +and race; a positive slack means no network satisfies the targets. + +Infeasible targets fail in \code{ergm} MCMLE estimation with "The MPLE does +not exist!" and then the simulated-density guard, usually after a +warning that \verb{nodefactor.race.} is linearly dependent on +\code{edges + nodematch.race}, which is this same boundary seen from the +other side. Stochastic-Approximation estimation settles on the nearest +reachable network without complaint, so the problem can go unnoticed +until a model is refit with MCMLE. + +\code{build_netstats()} runs this check by default (\code{check.feasibility = TRUE}) +and stops on a positive slack. The condition covers the +\code{nodematch("race", diff = FALSE)} specification used by the EpiModelHIV +template formulas; the per-race \code{nodematch_race} vector is not checked. +} +\examples{ +netstats <- readRDS(system.file("netstats-example.rds", package = "ARTnet")) +check_race_feasibility(netstats) + +} diff --git a/tests/testthat/test-joint-dyad.R b/tests/testthat/test-joint-dyad.R index 1a8f977..a417204 100644 --- a/tests/testthat/test-joint-dyad.R +++ b/tests/testthat/test-joint-dyad.R @@ -37,6 +37,7 @@ test_that("method='joint' produces all dyad models per layer", { for (layer in c("main", "casl", "inst")) { expect_s3_class(np[[layer]]$joint_nm_age_model, "glm") expect_s3_class(np[[layer]]$joint_nm_race_model, "glm") + expect_s3_class(np[[layer]]$joint_prace_model, "multinom") expect_s3_class(np[[layer]]$joint_absdiff_age_model, "glm") expect_s3_class(np[[layer]]$joint_absdiff_sqrtage_model, "glm") } @@ -48,6 +49,7 @@ test_that("method='existing' does not produce dyad models", { for (layer in c("main", "casl", "inst")) { expect_null(np[[layer]]$joint_nm_age_model) expect_null(np[[layer]]$joint_nm_race_model) + expect_null(np[[layer]]$joint_prace_model) expect_null(np[[layer]]$joint_absdiff_age_model) expect_null(np[[layer]]$joint_absdiff_sqrtage_model) } @@ -134,9 +136,10 @@ test_that("race = FALSE skips nm_race model gracefully", { ) set.seed(20260419L) np <- build_netparams(epistats, smooth.main.dur = TRUE, method = "joint") - # nm_age, absdiff models should be present; nm_race should not + # nm_age, absdiff models should be present; nm_race and prace should not expect_s3_class(np$main$joint_nm_age_model, "glm") expect_null(np$main$joint_nm_race_model) + expect_null(np$main$joint_prace_model) expect_s3_class(np$main$joint_absdiff_age_model, "glm") set.seed(20260419L) diff --git a/tests/testthat/test-race-mixing.R b/tests/testthat/test-race-mixing.R new file mode 100644 index 0000000..c182d78 --- /dev/null +++ b/tests/testthat/test-race-mixing.R @@ -0,0 +1,291 @@ +# Tests for the race-mixing helpers (issues #9, #10, #11): the balanced +# partner-race mixing matrix under method = "joint", the opt-in random-mixing +# floor, the feasibility check, and the build_netparams seed. + +skip_without_artnetdata <- function() { + testthat::skip_if(system.file(package = "ARTnetData") == "", + "ARTnetData not installed") +} + +# Seattle main-layer reported endpoint matrix from the planning prototype +# (joint, 10,000 nodes, default race.prop). Rows are the ego race (Black, +# Hispanic, White/Other), columns the partner race. White egos report more +# Black partners (231) than Black egos have endpoints in total (174). +seattle_main <- matrix(c(91.7, 7.3, 74.5, + 17.4, 75.7, 177.4, + 231.1, 277.4, 3181.7), + nrow = 3, byrow = TRUE) + +slack_of <- function(B) { + nf <- rowSums(B) + nf - sum(nf) / 2 - sum(diag(B)) / 2 +} + +# ---- .balance_race_mixing --------------------------------------------------- + +test_that("a symmetric matrix is a fixed point of the balance", { + A <- matrix(c(10, 2, 1, + 2, 20, 3, + 1, 3, 30), nrow = 3) + B <- .balance_race_mixing(A) + expect_equal(unclass(B)[, ], A, tolerance = 1e-9) + expect_equal(attr(B, "scaling"), rep(1, 3), tolerance = 1e-9) +}) + +test_that("the balance holds the margins, is symmetric, and is feasible on Seattle", { + A <- seattle_main + expect_gt(max(abs(A - t(A))), 100) # the reports do not balance + B <- .balance_race_mixing(A) + expect_equal(rowSums(B), rowSums(A), tolerance = 1e-8) + expect_equal(unclass(B)[, ], t(unclass(B)[, ]), tolerance = 1e-9) + expect_true(all(B >= 0)) + expect_true(all(diag(B) > 0)) + expect_true(all(slack_of(B) < 0)) + # the homophily odds ratios of the symmetrized input are preserved + S <- (A + t(A)) / 2 + or <- function(M, r, s) M[r, r] * M[s, s] / M[r, s]^2 + for (r in 1:2) for (s in (r + 1):3) { + expect_equal(or(B, r, s), or(S, r, s), tolerance = 1e-8) + } +}) + +test_that("plain averaging with the margins held goes negative where the balance does not", { + # Seattle casual from the prototype: averaging drives the Hispanic + # diagonal negative; the symmetric scaling keeps every cell positive. + A <- matrix(c(191.3, 51.5, 175.3, + 46.6, 181.6, 249.7, + 443.8, 961.6, 4131.8), nrow = 3, byrow = TRUE) + expect_gt(max(slack_of(A)), 0) # the unbalanced reports are infeasible as targets + S <- (A + t(A)) / 2 + off <- S + diag(off) <- 0 + avg_diag <- rowSums(A) - rowSums(off) + expect_lt(min(avg_diag), 0) + B <- .balance_race_mixing(A) + expect_true(all(diag(B) > 0)) + expect_true(all(slack_of(B) < 0)) +}) + +test_that("a race with no endpoints gets an empty row and column", { + A <- seattle_main + A[2, ] <- 0 + B <- .balance_race_mixing(A) + expect_equal(B[2, ], rep(0, 3)) + expect_equal(B[, 2], rep(0, 3)) + expect_equal(rowSums(B), rowSums(A), tolerance = 1e-8) + expect_equal(attr(B, "scaling")[2], 0) +}) + +test_that("the balance rejects bad input", { + expect_error(.balance_race_mixing(matrix(1, 2, 3)), "square") + expect_error(.balance_race_mixing(matrix(c(1, -1, 0, 1), 2)), "non-negative") +}) + +# ---- .reported_race_mixing / .predict_prace -------------------------------- + +test_that("the reported matrix sums predicted degree by ego race and partner race", { + pred_deg <- c(1, 2, 3) + probs <- rbind(c(0.5, 0.5), c(0.2, 0.8), c(0.9, 0.1)) + race <- c(1L, 2L, 2L) + A <- .reported_race_mixing(pred_deg, probs, race, 2) + expect_equal(A[1, ], c(0.5, 0.5)) + expect_equal(A[2, ], c(2 * 0.2 + 3 * 0.9, 2 * 0.8 + 3 * 0.1)) + expect_equal(rowSums(A), c(1, 5)) +}) + +test_that(".predict_prace returns an n x K matrix for two and three levels", { + set.seed(1) + d <- data.frame(x = rnorm(300)) + d$y3 <- factor(sample(1:3, 300, TRUE), levels = 1:3) + d$y2 <- factor(sample(1:2, 300, TRUE), levels = 1:2) + m3 <- nnet::multinom(y3 ~ x, data = d, trace = FALSE) + m2 <- nnet::multinom(y2 ~ x, data = d, trace = FALSE) + nd <- data.frame(x = c(-1, 0, 1)) + p3 <- .predict_prace(m3, nd, 3) + p2 <- .predict_prace(m2, nd, 2) + expect_equal(dim(p3), c(3, 3)) + expect_equal(dim(p2), c(3, 2)) + expect_equal(rowSums(p3), rep(1, 3), tolerance = 1e-9) + expect_equal(rowSums(p2), rep(1, 3), tolerance = 1e-9) + # single-row newdata + expect_equal(dim(.predict_prace(m3, nd[1, , drop = FALSE], 3)), c(1, 3)) + expect_equal(dim(.predict_prace(m2, nd[1, , drop = FALSE], 2)), c(1, 2)) +}) + +# ---- check_race_feasibility / floor ---------------------------------------- + +fake_netstats <- function(nf, edges, M) { + list(race = TRUE, + demog = list(props = data.frame(Black = 0.1, Hispanic = 0.1, White.Other = 0.8)), + main = list(edges = edges, nodefactor_race = nf, nodematch_race_diffF = M), + casl = list(edges = edges, nodefactor_race = nf, nodematch_race_diffF = M), + inst = list(edges = edges, nodefactor_race = nf, nodematch_race_diffF = M)) +} + +test_that("check_race_feasibility reports slack per layer and race", { + ns <- fake_netstats(nf = c(100, 100, 800), edges = 500, M = 400) + res <- check_race_feasibility(ns) + expect_s3_class(res, "data.frame") + expect_equal(nrow(res), 9) + expect_named(res, c("layer", "race", "label", "nodefactor", "edges", + "nodematch_race_diffF", "slack", "feasible")) + expect_equal(res$slack[res$layer == "main"], c(100, 100, 800) - 500 - 400) + expect_true(all(res$feasible)) + # labels: integer codes by default (demog$props columns are not in code + # order), explicit when supplied, from mixmat_race dimnames when present + expect_equal(unique(res$label), c("1", "2", "3")) + lab <- c("Black", "Hispanic", "White.Other") + expect_equal(unique(check_race_feasibility(ns, labels = lab)$label), lab) + ns$main$mixmat_race <- matrix(0, 3, 3, dimnames = list(lab, lab)) + res2 <- check_race_feasibility(ns) + expect_equal(res2$label[res2$layer == "main"], lab) + expect_equal(res2$label[res2$layer == "casl"], c("1", "2", "3")) + # M too small for the large race: nodefactor 800 > 500 + 250 + bad <- check_race_feasibility(fake_netstats(c(100, 100, 800), 500, 250)) + expect_false(all(bad$feasible)) + expect_true(all(bad$race[!bad$feasible] == 3)) + expect_equal(bad$slack[bad$race == 3], rep(50, 3)) +}) + +test_that("check_race_feasibility is empty without race targets", { + expect_message(res <- check_race_feasibility(list(race = FALSE, main = list(edges = 1))), + "nothing to check") + expect_equal(nrow(res), 0) +}) + +test_that("the random-mixing floor binds only below random mixing", { + layer <- list(edges = 500, nodefactor_race = c(100, 100, 800), + nodematch_race_diffF = 250) + p <- c(0.1, 0.1, 0.8) + rnd <- 500 * sum(p^2) + expect_message(out <- .apply_race_floor(layer, "main"), "random-mixing floor") + expect_equal(out$nodematch_race_diffF, rnd) + # feasible after flooring + expect_true(all(c(100, 100, 800) - 500 - rnd <= 0)) + above <- layer + above$nodematch_race_diffF <- 400 + expect_silent(out2 <- .apply_race_floor(above, "main")) + expect_equal(out2$nodematch_race_diffF, 400) + # layers without race targets pass through + expect_identical(.apply_race_floor(list(edges = 1), "inst"), list(edges = 1)) +}) + +# ---- end to end on Seattle (ARTnetData) ------------------------------------- + +seattle_env <- new.env(parent = emptyenv()) + +seattle_epistats <- function() { + if (is.null(seattle_env$epi)) { + set.seed(20260420L) + seattle_env$epi <- build_epistats( + geog.lvl = "city", geog.cat = "Seattle", + init.hiv.prev = c(0.33, 0.137, 0.084), + race = TRUE, time.unit = 7 + ) + } + seattle_env$epi +} + +seattle_netparams <- function(method) { + key <- paste0("np_", method) + if (is.null(seattle_env[[key]])) { + epi <- seattle_epistats() + dur <- if (method == "joint") "joint_lm" else "empirical" + seattle_env[[key]] <- build_netparams(epi, smooth.main.dur = TRUE, + method = method, duration.method = dur, + seed = 20260420L) + } + seattle_env[[key]] +} + +test_that("build_netstats stops on Seattle's infeasible legacy race targets", { + skip_without_artnetdata() + epi <- seattle_epistats() + np <- seattle_netparams("existing") + expect_error( + build_netstats(epi, np, network.size = 5000, expect.mort = 0.00025, + seed = 20260420L), + regexp = "infeasible.*nodefactor_race\\[3\\] \\(White.Other\\)" + ) + ns <- build_netstats(epi, np, network.size = 5000, expect.mort = 0.00025, + seed = 20260420L, check.feasibility = FALSE) + feas <- check_race_feasibility(ns) + expect_false(all(feas$feasible)) + expect_true(all(feas$feasible[feas$race != 3])) + expect_true(all(!feas$feasible[feas$race == 3])) +}) + +test_that("the floor makes Seattle's legacy race targets feasible", { + skip_without_artnetdata() + epi <- seattle_epistats() + np <- seattle_netparams("existing") + expect_message( + ns <- build_netstats(epi, np, network.size = 5000, expect.mort = 0.00025, + seed = 20260420L, nodematch.race.floor = TRUE), + "random-mixing floor" + ) + feas <- check_race_feasibility(ns) + expect_true(all(feas$feasible)) + for (layer in c("main", "casl", "inst")) { + p <- ns[[layer]]$nodefactor_race / sum(ns[[layer]]$nodefactor_race) + expect_equal(ns[[layer]]$nodematch_race_diffF, + ns[[layer]]$edges * sum(p^2), tolerance = 1e-9, + info = layer) + } +}) + +test_that("the joint method balances Seattle's race mixing and passes the check", { + skip_without_artnetdata() + epi <- seattle_epistats() + np <- seattle_netparams("joint") + ns <- build_netstats(epi, np, network.size = 5000, expect.mort = 0.00025, + method = "joint", seed = 20260420L) + feas <- check_race_feasibility(ns) + expect_true(all(feas$feasible)) + for (layer in c("main", "casl", "inst")) { + L <- ns[[layer]] + expect_s3_class(np[[layer]]$joint_prace_model, "multinom") + expect_equal(sum(L$nodematch_race), L$nodematch_race_diffF, tolerance = 1e-9) + expect_true(all(L$nodematch_race <= L$nodefactor_race / 2 + 1e-9)) + expect_true(all(L$nodematch_race > 0)) + # balanced matrix: symmetric edge counts; each race's endpoints + # (twice the diagonal plus the off-diagonal row) equal nodefactor_race, + # and the distinct cells sum to edges + M <- L$mixmat_race + expect_equal(dim(M), c(3, 3)) + expect_equal(M, t(M), tolerance = 1e-9) + expect_equal(unname(2 * diag(M) + (rowSums(M) - diag(M))), L$nodefactor_race, + tolerance = 1e-8, info = layer) + expect_equal(sum(M[upper.tri(M, diag = TRUE)]), L$edges, tolerance = 1e-8) + # reported matrix: rows sum to nodefactor_race, and it is not symmetric + R <- L$mixmat_race_reported + expect_equal(unname(rowSums(R)), L$nodefactor_race, tolerance = 1e-8) + expect_gt(max(abs(R - t(R))), 1) + expect_equal(rownames(M), c("Black", "Hispanic", "White.Other")) + # the checker picks the labels up from the matrix + expect_equal(feas$label[feas$layer == layer], rownames(M)) + } + # same-race shares are assortative relative to random mixing in every layer + for (layer in c("main", "casl", "inst")) { + L <- ns[[layer]] + p <- L$nodefactor_race / sum(L$nodefactor_race) + expect_gt(L$nodematch_race_diffF, L$edges * sum(p^2)) + } +}) + +test_that("build_netparams(seed) makes the partner-race imputation reproducible", { + skip_without_artnetdata() + epi <- seattle_epistats() + a <- build_netparams(epi, smooth.main.dur = TRUE, seed = 1L) + b <- build_netparams(epi, smooth.main.dur = TRUE, seed = 1L) + c2 <- build_netparams(epi, smooth.main.dur = TRUE, seed = 2L) + for (layer in c("main", "casl", "inst")) { + expect_identical(a[[layer]]$nm.race, b[[layer]]$nm.race) + expect_identical(a[[layer]]$nm.race_diffF, b[[layer]]$nm.race_diffF) + expect_identical(a[[layer]]$md.main, c2[[layer]]$md.main) + } + moved <- vapply(c("main", "casl", "inst"), + function(l) !isTRUE(all.equal(a[[l]]$nm.race, c2[[l]]$nm.race)), + logical(1)) + expect_true(any(moved)) +}) From 9b8bb21face77b7adcf52d1dd44051bd011eb4f1 Mon Sep 17 00:00:00 2001 From: Samuel Jenness Date: Sat, 7 Nov 2026 09:41:54 -0500 Subject: [PATCH 2/4] Fit every scenario city with MCMLE in the method comparison fit_comparison_ergms() builds each scenario's targets under three arms (existing, existing with the random-mixing floor, joint), runs check_race_feasibility() on each, fits the three template ERGMs with MCMLE for the feasible ones, and records convergence and static netdx recovery. render_fit_report() writes method_comparison_fits.md. The comparison report is regenerated for the balanced race targets (#12). --- inst/validation/README.md | 18 + inst/validation/method_comparison.R | 421 +++++++++++++++++++++- inst/validation/method_comparison.md | 64 ++-- inst/validation/method_comparison_fits.md | 284 +++++++++++++++ 4 files changed, 741 insertions(+), 46 deletions(-) create mode 100644 inst/validation/method_comparison_fits.md diff --git a/inst/validation/README.md b/inst/validation/README.md index ab77d46..2756efd 100644 --- a/inst/validation/README.md +++ b/inst/validation/README.md @@ -70,3 +70,21 @@ Two reasons: Unit tests for individual joint-GLM behaviors (convergence, marginal recovery, coefficient sanity) should still live in `tests/testthat/` as normal. + +## MCMLE fits by scenario city (issue #12) + +`method_comparison.R` also fits the three template ERGMs for every scenario +city with MCMLE, under three arms (`existing`, `existing_floor`, `joint`), +after running `check_race_feasibility()` on each target set: + +```r +devtools::load_all() +source(system.file("validation/method_comparison.R", package = "ARTnet")) +fits <- fit_comparison_ergms(ncores = 4) # ~15-20 min +render_fit_report(fits) # writes inst/validation/method_comparison_fits.md +``` + +Stochastic-Approximation settles on the nearest reachable network when the +targets are infeasible, so the earlier Atlanta-only SA fits could not have +caught the Seattle race-mixing problem (issue #9). MCMLE fails on it, which +is what the cluster runs use. diff --git a/inst/validation/method_comparison.R b/inst/validation/method_comparison.R index a99c55f..be79fea 100644 --- a/inst/validation/method_comparison.R +++ b/inst/validation/method_comparison.R @@ -103,24 +103,35 @@ COMPARISON_SCENARIOS <- lapply( do.call(rbind, rows[!vapply(rows, is.null, logical(1))]) } +# Build the netstats object for one scenario under one method. The seed is +# reset before each build_* call so the three stochastic steps (HIV model, +# partner-race imputation, attribute sampling) are reproducible and the +# same across methods. The feasibility check is off so that infeasible +# target sets (Seattle under method = "existing") come back and can be +# reported rather than stopping the run; fit_comparison_ergms() runs the +# check itself. `floor = TRUE` turns on the random-mixing floor (#9). +.build_scenario_netstats <- function(scenario, method, floor = FALSE) { + dur_method <- if (method == "joint") "joint_lm" else "empirical" + set.seed(.COMPARISON_SEED) + epistats <- do.call(ARTnet::build_epistats, scenario$epistats) + np_args <- c(list(epistats = epistats), scenario$netparams, + list(method = method, duration.method = dur_method, + seed = .COMPARISON_SEED)) + netparams <- do.call(ARTnet::build_netparams, np_args) + ns_args <- c(list(epistats = epistats, netparams = netparams), + scenario$netstats, + list(method = method, seed = .COMPARISON_SEED, + nodematch.race.floor = floor, + check.feasibility = FALSE)) + suppressMessages(do.call(ARTnet::build_netstats, ns_args)) +} + # Run one scenario through both methods, return both extracted stat tables. .run_one_scenario <- function(scenario) { - do_run <- function(np_method, dur_method) { - set.seed(.COMPARISON_SEED) - epistats <- do.call(ARTnet::build_epistats, scenario$epistats) - set.seed(.COMPARISON_SEED) - np_args <- c(list(epistats = epistats), scenario$netparams, - list(method = np_method, duration.method = dur_method)) - netparams <- do.call(ARTnet::build_netparams, np_args) - set.seed(.COMPARISON_SEED) - ns_args <- c(list(epistats = epistats, netparams = netparams), - scenario$netstats, list(method = np_method)) - do.call(ARTnet::build_netstats, ns_args) - } list( - existing = .extract_target_stats(do_run("existing", "empirical"), + existing = .extract_target_stats(.build_scenario_netstats(scenario, "existing"), scenario$name, "existing"), - joint = .extract_target_stats(do_run("joint", "joint_lm"), + joint = .extract_target_stats(.build_scenario_netstats(scenario, "joint"), scenario$name, "joint") ) } @@ -300,3 +311,385 @@ render_comparison_report <- function(comparison, out("```") invisible(file) } + + +# ---- ERGM fits with MCMLE (issue #12) --------------------------------------- +# +# Fits the three EpiModelHIV-Template ERGMs (formulas and target mappings +# as pinned in inst/validation/epimodelhiv_template_ref/) for every +# scenario x arm with MCMLE, the production estimation method, after running +# check_race_feasibility() on each target set. Stochastic-Approximation +# settles on the nearest reachable network when targets are infeasible, so +# it could not have caught the Seattle problem (#9); MCMLE fails on it. +# +# Arms: +# existing legacy targets +# existing_floor legacy targets with nodematch.race.floor = TRUE +# joint joint g-computation with balanced partner-race mixing +# +# Infeasible target sets are recorded and not fit. Each fit records whether +# netest() returned, whether MCMLE converged, whether ergm warned about a +# linear dependence (nodefactor.race = edges + nodematch.race, the +# feasibility boundary seen from the other side), the fit time, and from a +# static netdx the largest |Z| and |% diff| over the formation targets. +# +# Usage: +# fits <- fit_comparison_ergms(ncores = 4) # ~15-20 min on a laptop +# render_fit_report(fits) # writes method_comparison_fits.md +# +# Requires EpiModel (and through it ergm) in addition to ARTnetData. + +.ERGM_CONTROL <- function() { + ergm::control.ergm( + main.method = "MCMLE", + MCMLE.maxit = 500, + SAN.maxit = 3, + SAN.nsteps.times = 4, + MCMC.samplesize = 1e4, + MCMC.interval = 5e3, + # 0, not 1: ergm starts a one-worker PSOCK cluster at parallel = 1, and + # concurrent jobs under mclapply then collide on the cluster port. + parallel = 0 + ) +} + +# Formation formula, target.stats vector and dissolution coefs for one +# layer, verbatim from the template (tests/workflows/workflow-standard.R). +.template_layer_spec <- function(netstats, layer) { + ns <- netstats[[layer]] + if (layer == "main") { + list( + formation = ~ edges + + nodematch("age.grp", diff = TRUE) + + nodefactor("age.grp", levels = -1) + + nodematch("race", diff = FALSE) + + nodefactor("race", levels = -1) + + nodefactor("deg.casl", levels = -1) + + concurrent + + degrange(from = 3) + + nodematch("role.class", diff = TRUE, levels = c(1, 2)), + target.stats = unname(c( + ns$edges, ns$nodematch_age.grp, ns$nodefactor_age.grp[-1], + ns$nodematch_race_diffF, ns$nodefactor_race[-1], + ns$nodefactor_deg.casl[-1], ns$concurrent, 0, 0, 0)), + coef.diss = ns$diss.byage + ) + } else if (layer == "casl") { + list( + formation = ~ edges + + nodematch("age.grp", diff = TRUE) + + nodefactor("age.grp", levels = -5) + + nodematch("race", diff = FALSE) + + nodefactor("race", levels = -1) + + nodefactor("deg.main", levels = -3) + + concurrent + + degrange(from = 4) + + nodematch("role.class", diff = TRUE, levels = c(1, 2)), + target.stats = unname(c( + ns$edges, ns$nodematch_age.grp, ns$nodefactor_age.grp[-5], + ns$nodematch_race_diffF, ns$nodefactor_race[-1], + ns$nodefactor_deg.main[-3], ns$concurrent, 0, 0, 0)), + coef.diss = ns$diss.byage + ) + } else { + list( + formation = ~ edges + + nodematch("age.grp", diff = FALSE) + + nodefactor("age.grp", levels = -1) + + nodematch("race", diff = FALSE) + + nodefactor("race", levels = -1) + + nodefactor("risk.grp", levels = -5) + + nodefactor("deg.tot", levels = -1) + + nodematch("role.class", diff = TRUE, levels = c(1, 2)), + target.stats = unname(c( + ns$edges, sum(ns$nodematch_age.grp), ns$nodefactor_age.grp[-1], + ns$nodematch_race_diffF, ns$nodefactor_race[-1], + ns$nodefactor_risk.grp[-5], ns$nodefactor_deg.tot[-1], 0, 0)), + coef.diss = EpiModel::dissolution_coefs(~ offset(edges), 1) + ) + } +} + +.empty_fit_row <- function() { + data.frame(fit_ok = NA, converged = NA, lin_dep_warning = NA, + fit_minutes = NA_real_, netdx_max_abs_z = NA_real_, + netdx_max_abs_pct = NA_real_, error = NA_character_, + stringsAsFactors = FALSE) +} + +# Fit one layer with MCMLE and run a static netdx on the formation targets. +.fit_one_layer <- function(netstats, layer, nsims_dx = 1000) { + spec <- .template_layer_spec(netstats, layer) + nw <- EpiModel::network_initialize(netstats$demog$num) + nw <- EpiModel::set_vertex_attribute(nw, names(netstats$attr), netstats$attr) + warns <- character(0) + t0 <- Sys.time() + fit <- withCallingHandlers( + tryCatch( + EpiModel::netest(nw, formation = spec$formation, + target.stats = spec$target.stats, + coef.diss = spec$coef.diss, + set.control.ergm = .ERGM_CONTROL(), + verbose = FALSE), + error = function(e) { + structure(list(message = conditionMessage(e)), class = "fit_error") + } + ), + warning = function(w) { + warns <<- c(warns, conditionMessage(w)) + invokeRestart("muffleWarning") + }, + message = function(m) invokeRestart("muffleMessage") + ) + row <- .empty_fit_row() + row$fit_minutes <- round(as.numeric(difftime(Sys.time(), t0, units = "mins")), 2) + row$lin_dep_warning <- any(grepl("linear", warns, ignore.case = TRUE)) + if (inherits(fit, "fit_error")) { + row$fit_ok <- FALSE + row$converged <- FALSE + row$error <- gsub("\\s+", " ", fit$message) + return(row) + } + row$fit_ok <- TRUE + row$converged <- !any(grepl("did not converge", warns, ignore.case = TRUE)) + dx <- tryCatch( + withCallingHandlers( + EpiModel::netdx(fit, dynamic = FALSE, nsims = nsims_dx, + nwstats.formula = "formation", + set.control.ergm = ergm::control.simulate.formula(MCMC.burnin = 1e5), + verbose = FALSE), + message = function(m) invokeRestart("muffleMessage"), + warning = function(w) invokeRestart("muffleWarning") + ), + error = function(e) NULL + ) + if (!is.null(dx)) { + st <- dx$stats.table.formation + z <- st[, grep("^Z", colnames(st))[1]] + pct <- st[, grep("^Pct", colnames(st))[1]] + ok <- is.finite(z) & is.finite(pct) + row$netdx_max_abs_z <- round(max(abs(z[ok])), 2) + row$netdx_max_abs_pct <- round(max(abs(pct[ok])), 2) + } + row +} + +# Build netstats for one scenario x arm, check feasibility, fit what is +# feasible. Returns a list(rows = data.frame, mixmat = list or NULL). +.run_one_fit_job <- function(scenario, arm, nsims_dx = 1000) { + method <- if (arm == "joint") "joint" else "existing" + ns <- .build_scenario_netstats(scenario, method, floor = arm == "existing_floor") + feas <- ARTnet::check_race_feasibility(ns) + rows <- lapply(c("main", "casl", "inst"), function(l) { + f <- feas[feas$layer == l, , drop = FALSE] + worst <- which.max(f$slack) + base <- data.frame( + scenario = scenario$name, arm = arm, layer = l, + feasible = all(f$feasible), + worst_slack = round(max(f$slack), 1), + worst_race = f$race[worst], + edges = round(ns[[l]]$edges, 1), + nodematch_race_diffF = round(ns[[l]]$nodematch_race_diffF, 1), + stringsAsFactors = FALSE + ) + if (!all(f$feasible)) { + r <- .empty_fit_row() + r$fit_ok <- FALSE + r$converged <- FALSE + r$error <- "skipped: infeasible race targets" + } else { + r <- .fit_one_layer(ns, l, nsims_dx) + } + cbind(base, r) + }) + mixmat <- NULL + if (method == "joint") { + mixmat <- lapply(c(main = "main", casl = "casl", inst = "inst"), function(l) { + list(reported = ns[[l]]$mixmat_race_reported, balanced = ns[[l]]$mixmat_race) + }) + } + list(rows = do.call(rbind, rows), mixmat = mixmat) +} + + +#' Fit the template ERGMs with MCMLE for every scenario and arm. +#' +#' @param scenarios List of scenario specs (defaults to COMPARISON_SCENARIOS). +#' @param arms Character vector of arms to run; see the header comment. +#' @param ncores Number of cores for parallel::mclapply over scenario x arm +#' jobs. Each ergm fit itself runs on one core. +#' @param nsims_dx Number of static netdx simulations per fitted layer. +#' @param save_to Optional path for saveRDS() of the result. +#' @return A list with `fits` (one row per scenario x arm x layer), `mixmats` +#' (the joint arm's reported and balanced race mixing matrices by scenario +#' and layer), and run metadata. +fit_comparison_ergms <- function(scenarios = COMPARISON_SCENARIOS, + arms = c("existing", "existing_floor", "joint"), + ncores = 1, + nsims_dx = 1000, + save_to = NULL) { + .require_artnetdata() + if (!requireNamespace("EpiModel", quietly = TRUE)) { + stop("EpiModel is required to fit the ERGMs.") + } + jobs <- expand.grid(scenario = seq_along(scenarios), arm = arms, + stringsAsFactors = FALSE) + message("Fitting ", nrow(jobs), " scenario x arm jobs (3 layers each) with MCMLE ", + "on ", ncores, " core(s)...") + run_job <- function(i) { + s <- scenarios[[jobs$scenario[i]]] + arm <- jobs$arm[i] + tryCatch( + .run_one_fit_job(s, arm, nsims_dx), + error = function(e) { + rows <- data.frame(scenario = s$name, arm = arm, + layer = c("main", "casl", "inst"), + feasible = NA, worst_slack = NA_real_, worst_race = NA_integer_, + edges = NA_real_, nodematch_race_diffF = NA_real_, + stringsAsFactors = FALSE) + r <- .empty_fit_row()[rep(1, 3), ] + r$fit_ok <- FALSE + r$error <- paste("job error:", gsub("\\s+", " ", conditionMessage(e))) + list(rows = cbind(rows, r), mixmat = NULL) + } + ) + } + res <- if (ncores > 1) { + parallel::mclapply(seq_len(nrow(jobs)), run_job, mc.cores = ncores) + } else { + lapply(seq_len(nrow(jobs)), function(i) { + message(" ", scenarios[[jobs$scenario[i]]]$name, " / ", jobs$arm[i]) + run_job(i) + }) + } + fits <- do.call(rbind, lapply(res, `[[`, "rows")) + rownames(fits) <- NULL + mixmats <- list() + for (i in seq_along(res)) { + if (!is.null(res[[i]]$mixmat)) { + mixmats[[scenarios[[jobs$scenario[i]]]$name]] <- res[[i]]$mixmat + } + } + out <- list( + fits = fits, + mixmats = mixmats, + network.size = scenarios[[1]]$netstats$network.size, + seed = .COMPARISON_SEED, + artnet_version = as.character(packageVersion("ARTnet")), + ergm_version = as.character(packageVersion("ergm")), + date = as.character(Sys.Date()) + ) + if (!is.null(save_to)) saveRDS(out, save_to) + out +} + + +.fmt_matrix_md <- function(M, digits = 1) { + if (is.null(M)) return("_not available_") + labs <- rownames(M) + if (is.null(labs)) labs <- as.character(seq_len(nrow(M))) + lines <- c(paste0("| | ", paste(labs, collapse = " | "), " |"), + paste0("|---|", paste(rep("---:", ncol(M)), collapse = "|"), "|")) + for (i in seq_len(nrow(M))) { + lines <- c(lines, paste0("| ", labs[i], " | ", + paste(formatC(M[i, ], format = "f", digits = digits), + collapse = " | "), " |")) + } + paste(lines, collapse = "\n") +} + +#' Render a markdown report of the ERGM fit results. +render_fit_report <- function(fits_result, + file = "inst/validation/method_comparison_fits.md") { + fits <- fits_result$fits + con <- file(file, "w") + on.exit(close(con)) + out <- function(...) cat(..., "\n", sep = "", file = con) + + out("# Method comparison: MCMLE fits of the template ERGMs by scenario city") + out() + out("Generated by `fit_comparison_ergms()` in `inst/validation/method_comparison.R` on ", + fits_result$date, ". Addresses issue #12.") + out() + out("ARTnet version: ", fits_result$artnet_version, ". ergm version: ", + fits_result$ergm_version, ". Seed: ", fits_result$seed, + ". Network size: ", fits_result$network.size, ".") + out() + out("Every scenario x arm target set is first run through `check_race_feasibility()`. ", + "Feasible sets are fit for all three layers with the template formulas and ", + "`control.ergm(main.method = \"MCMLE\", MCMLE.maxit = 500, SAN.maxit = 3, ", + "SAN.nsteps.times = 4, MCMC.samplesize = 1e4, MCMC.interval = 5e3)`, then checked ", + "with a static `netdx` on the formation targets. Infeasible sets are skipped.") + out() + out("Arms: `existing` is the legacy per-attribute method; `existing_floor` is the same ", + "with `nodematch.race.floor = TRUE`; `joint` is the joint g-computation with ", + "partner-race mixing balanced to the target population.") + out() + out("## Summary") + out() + n_fit <- sum(fits$fit_ok, na.rm = TRUE) + n_conv <- sum(fits$converged, na.rm = TRUE) + n_skip <- sum(!is.na(fits$feasible) & !fits$feasible) + out("- Target sets (scenario x arm x layer): ", nrow(fits)) + out("- Skipped as infeasible: ", n_skip) + out("- Fits returned: ", n_fit, "; MCMLE converged: ", n_conv) + by_arm <- split(fits, fits$arm) + for (a in names(by_arm)) { + f <- by_arm[[a]] + out("- `", a, "`: ", sum(!is.na(f$feasible) & !f$feasible), " infeasible, ", + sum(f$converged, na.rm = TRUE), " of ", nrow(f), " converged, max netdx |Z| ", + if (all(is.na(f$netdx_max_abs_z))) "n/a" else sprintf("%.2f", max(f$netdx_max_abs_z, na.rm = TRUE)), + ", max |% diff| ", + if (all(is.na(f$netdx_max_abs_pct))) "n/a" else sprintf("%.2f%%", max(f$netdx_max_abs_pct, na.rm = TRUE))) + } + out() + out("## Fits by scenario") + out() + tf <- function(x) ifelse(is.na(x), "n/a", ifelse(x, "yes", "no")) + num <- function(x, d = 1) ifelse(is.na(x), "n/a", formatC(x, format = "f", digits = d)) + for (s in unique(fits$scenario)) { + out("### ", s) + out() + out("| Arm | Layer | Feasible | Worst slack | Edges | Same-race target | Fit | Converged | ", + "Lin. dep. warning | Minutes | netdx max abs Z | netdx max abs % | Note |") + out("|---|---|---|---:|---:|---:|---|---|---|---:|---:|---:|---|") + sub <- fits[fits$scenario == s, , drop = FALSE] + for (i in seq_len(nrow(sub))) { + r <- sub[i, ] + out(sprintf("| %s | %s | %s | %s | %s | %s | %s | %s | %s | %s | %s | %s | %s |", + r$arm, r$layer, tf(r$feasible), num(r$worst_slack), num(r$edges), + num(r$nodematch_race_diffF), tf(r$fit_ok), tf(r$converged), + tf(r$lin_dep_warning), num(r$fit_minutes, 2), num(r$netdx_max_abs_z, 2), + num(r$netdx_max_abs_pct, 2), + if (is.na(r$error)) "" else r$error)) + } + out() + mm <- fits_result$mixmats[[s]] + if (!is.null(mm)) { + out("Joint arm race mixing matrices (rows = ego race). Reported: expected partnership ", + "endpoints by ego race and partner race at the synthetic population, standardized ", + "on the ego side only; row sums are `nodefactor_race`. Balanced: expected edge ", + "counts after symmetric scaling with `nodefactor_race` held fixed; the diagonal is ", + "`nodematch_race`.") + out() + for (l in c("main", "casl", "inst")) { + out("**", l, ", reported endpoints**") + out() + out(.fmt_matrix_md(mm[[l]]$reported)) + out() + out("**", l, ", balanced edges**") + out() + out(.fmt_matrix_md(mm[[l]]$balanced)) + out() + } + } + } + out("## Reproducibility") + out() + out("```r") + out("source(system.file(\"validation/method_comparison.R\", package = \"ARTnet\"))") + out("fits <- fit_comparison_ergms(ncores = 4)") + out("render_fit_report(fits)") + out("```") + invisible(file) +} diff --git a/inst/validation/method_comparison.md b/inst/validation/method_comparison.md index 305c719..d516a37 100644 --- a/inst/validation/method_comparison.md +++ b/inst/validation/method_comparison.md @@ -1,8 +1,8 @@ # Method comparison: marginal vs joint g-computation -Generated by `inst/validation/method_comparison.R` on 2026-06-08. Phase 1.5 of the joint g-comp refactor; closes part of issue #65. +Generated by `inst/validation/method_comparison.R` on 2026-11-07. Phase 1.5 of the joint g-comp refactor; closes part of issue #65. -ARTnet version: 2.9.0. Seed: 20260420. Network size: 5000. +ARTnet version: 2.10.0. Seed: 20260420. Network size: 5000. ## Scenarios @@ -18,12 +18,12 @@ Derived from the scenarios actually present in the comparison object, so this ta ## High-level summary - Total target-stat cells across scenarios: 384 -- Cells where |joint vs existing %diff| > 5%: **232** +- Cells where |joint vs existing %diff| > 5%: **253** -- `atlanta_default`: 96 cells, 63 materially shifted (>5%) -- `boston_default`: 96 cells, 64 materially shifted (>5%) -- `chicago_default`: 96 cells, 61 materially shifted (>5%) -- `seattle_default`: 96 cells, 44 materially shifted (>5%) +- `atlanta_default`: 96 cells, 68 materially shifted (>5%) +- `boston_default`: 96 cells, 66 materially shifted (>5%) +- `chicago_default`: 96 cells, 66 materially shifted (>5%) +- `seattle_default`: 96 cells, 53 materially shifted (>5%) ## Per-scenario top shifts (by |% diff|) @@ -31,7 +31,11 @@ Derived from the scenarios actually present in the comparison object, so this ta | Layer | Stat | Level | Existing | Joint | %Δ | |---|---|---:|---:|---:|---:| +| inst | nodematch_race | 2 | 3.64 | 0.67 | -81.7% | +| casl | nodematch_race | 1 | 368.50 | 579.90 | +57.4% | +| casl | nodematch_race | 2 | 21.54 | 10.29 | -52.2% | | inst | nodematch_age.grp | 5 | 7.90 | 3.86 | -51.1% | +| main | nodematch_race | 2 | 22.14 | 11.04 | -50.1% | | main | dissolution_duration | 6 | 934.45 | 491.48 | -47.4% | | inst | nodefactor_diag.status | 2 | 150.36 | 89.69 | -40.3% | | main | dissolution_duration | 4 | 539.20 | 323.74 | -40.0% | @@ -39,42 +43,41 @@ Derived from the scenarios actually present in the comparison object, so this ta | main | dissolution_duration | 5 | 682.59 | 428.11 | -37.3% | | casl | dissolution_duration | 2 | 50.44 | 68.93 | +36.7% | | inst | nodefactor_age.grp | 1 | 46.94 | 63.36 | +35.0% | +| inst | nodematch_race | 1 | 38.16 | 51.16 | +34.1% | | inst | nodefactor_age.grp | 5 | 85.81 | 57.25 | -33.3% | -| inst | nodematch_age.grp | 4 | 11.22 | 7.99 | -28.8% | -| casl | dissolution_duration | 4 | 113.16 | 81.76 | -27.7% | -| casl | dissolution_duration | 3 | 73.17 | 92.29 | +26.1% | -| casl | nodefactor_race | 2 | 119.99 | 149.17 | +24.3% | -| casl | dissolution_duration | 6 | 150.10 | 186.04 | +23.9% | -| inst | nodefactor_age.grp | 4 | 85.36 | 65.21 | -23.6% | +| main | nodematch_race | 1 | 202.90 | 266.02 | +31.1% | -_...and 48 more cells over threshold._ +_...and 53 more cells over threshold._ ### boston_default | Layer | Stat | Level | Existing | Joint | %Δ | |---|---|---:|---:|---:|---:| | main | dissolution_duration | 6 | 788.80 | 394.43 | -50.0% | +| casl | nodematch_race | 2 | 88.01 | 126.75 | +44.0% | | inst | nodematch_age.grp | 5 | 9.66 | 5.41 | -44.0% | | main | dissolution_duration | 4 | 455.19 | 260.07 | -42.9% | | main | dissolution_duration | 5 | 576.22 | 343.82 | -40.3% | | casl | nodefactor_deg.main | 3 | 22.83 | 31.88 | +39.6% | | inst | nodefactor_diag.status | 2 | 138.96 | 88.21 | -36.5% | +| casl | nodematch_race | 1 | 146.86 | 199.14 | +35.6% | +| inst | nodematch_race | 1 | 13.53 | 18.28 | +35.1% | | inst | nodefactor_age.grp | 1 | 63.38 | 82.87 | +30.8% | | main | dissolution_duration | 3 | 214.19 | 157.61 | -26.4% | | inst | nodefactor_age.grp | 5 | 115.85 | 86.20 | -25.6% | | casl | absdiff_age | — | 15284.68 | 18886.97 | +23.6% | +| inst | nodematch_race | 2 | 12.79 | 15.76 | +23.1% | | casl | nodefactor_race | 2 | 546.68 | 669.29 | +22.4% | -| main | nodematch_race_diffF | — | 679.50 | 530.72 | -21.9% | -| casl | dissolution_duration | 5 | 145.38 | 113.85 | -21.7% | -| casl | absdiff_sqrt.age | — | 1245.73 | 1512.13 | +21.4% | -| main | concurrent | — | 90.09 | 108.99 | +21.0% | -_...and 49 more cells over threshold._ +_...and 51 more cells over threshold._ ### chicago_default | Layer | Stat | Level | Existing | Joint | %Δ | |---|---|---:|---:|---:|---:| +| main | nodematch_race | 2 | 88.62 | 169.59 | +91.4% | +| casl | nodematch_race | 2 | 165.67 | 293.12 | +76.9% | +| main | nodematch_race | 1 | 80.44 | 128.75 | +60.0% | | main | dissolution_duration | 6 | 1034.43 | 492.95 | -52.3% | | main | dissolution_duration | 4 | 596.86 | 324.98 | -45.6% | | inst | nodematch_age.grp | 5 | 11.55 | 6.32 | -45.3% | @@ -84,36 +87,33 @@ _...and 49 more cells over threshold._ | inst | nodefactor_age.grp | 1 | 77.34 | 103.64 | +34.0% | | casl | nodefactor_deg.main | 3 | 25.95 | 34.03 | +31.2% | | casl | dissolution_duration | 1 | 92.48 | 64.08 | -30.7% | +| casl | nodematch_race | 1 | 256.25 | 333.58 | +30.2% | | main | dissolution_duration | 3 | 280.77 | 197.07 | -29.8% | | casl | dissolution_duration | 4 | 98.86 | 69.92 | -29.3% | -| inst | nodefactor_age.grp | 5 | 140.35 | 103.65 | -26.1% | -| inst | nodematch_race_diffF | — | 212.84 | 157.75 | -25.9% | -| main | nodematch_race_diffF | — | 592.50 | 449.49 | -24.1% | -| main | concurrent | — | 82.87 | 99.96 | +20.6% | -_...and 46 more cells over threshold._ +_...and 51 more cells over threshold._ ### seattle_default | Layer | Stat | Level | Existing | Joint | %Δ | |---|---|---:|---:|---:|---:| +| inst | nodematch_race | 2 | 4.51 | 1.43 | -68.2% | | main | dissolution_duration | 6 | 1380.99 | 471.87 | -65.8% | | main | dissolution_duration | 4 | 796.76 | 311.16 | -60.9% | | main | dissolution_duration | 5 | 1008.72 | 411.30 | -59.2% | +| casl | nodematch_race | 2 | 34.94 | 15.52 | -55.6% | | main | dissolution_duration | 3 | 374.72 | 188.65 | -49.7% | +| main | nodematch_race | 2 | 25.61 | 14.32 | -44.1% | | casl | dissolution_duration | 5 | 153.51 | 86.81 | -43.5% | +| inst | nodematch_race | 1 | 4.12 | 2.39 | -41.9% | +| main | nodematch_race | 1 | 20.79 | 12.11 | -41.7% | | inst | nodematch_age.grp | 5 | 7.70 | 4.66 | -39.4% | | casl | dissolution_duration | 1 | 101.86 | 63.71 | -37.5% | | casl | dissolution_duration | 4 | 108.89 | 69.18 | -36.5% | +| casl | nodematch_race | 1 | 50.13 | 32.39 | -35.4% | | main | dissolution_duration | 2 | 107.02 | 74.70 | -30.2% | -| inst | nodefactor_age.grp | 1 | 56.39 | 71.43 | +26.7% | -| inst | nodefactor_diag.status | 2 | 88.57 | 65.13 | -26.5% | -| casl | dissolution_duration | 6 | 144.43 | 107.50 | -25.6% | -| casl | nodefactor_deg.main | 3 | 24.27 | 30.42 | +25.4% | -| main | dissolution_duration | 1 | 307.95 | 233.27 | -24.2% | -| casl | dissolution_duration | 3 | 70.42 | 54.00 | -23.3% | - -_...and 29 more cells over threshold._ + +_...and 38 more cells over threshold._ ## Edges / concurrent comparison (all scenarios) diff --git a/inst/validation/method_comparison_fits.md b/inst/validation/method_comparison_fits.md new file mode 100644 index 0000000..291e3ac --- /dev/null +++ b/inst/validation/method_comparison_fits.md @@ -0,0 +1,284 @@ +# Method comparison: MCMLE fits of the template ERGMs by scenario city + +Generated by `fit_comparison_ergms()` in `inst/validation/method_comparison.R` on 2026-11-07. Addresses issue #12. + +ARTnet version: 2.10.0. ergm version: 4.12.0. Seed: 20260420. Network size: 5000. + +Every scenario x arm target set is first run through `check_race_feasibility()`. Feasible sets are fit for all three layers with the template formulas and `control.ergm(main.method = "MCMLE", MCMLE.maxit = 500, SAN.maxit = 3, SAN.nsteps.times = 4, MCMC.samplesize = 1e4, MCMC.interval = 5e3)`, then checked with a static `netdx` on the formation targets. Infeasible sets are skipped. + +Arms: `existing` is the legacy per-attribute method; `existing_floor` is the same with `nodematch.race.floor = TRUE`; `joint` is the joint g-computation with partner-race mixing balanced to the target population. + +## Summary + +- Target sets (scenario x arm x layer): 36 +- Skipped as infeasible: 3 +- Fits returned: 33; MCMLE converged: 33 +- `existing`: 3 infeasible, 9 of 12 converged, max netdx |Z| 3.67, max |% diff| 6.69% +- `existing_floor`: 0 infeasible, 12 of 12 converged, max netdx |Z| 4.45, max |% diff| 6.69% +- `joint`: 0 infeasible, 12 of 12 converged, max netdx |Z| 4.40, max |% diff| 6.30% + +## Fits by scenario + +### atlanta_default + +| Arm | Layer | Feasible | Worst slack | Edges | Same-race target | Fit | Converged | Lin. dep. warning | Minutes | netdx max abs Z | netdx max abs % | Note | +|---|---|---|---:|---:|---:|---|---|---|---:|---:|---:|---| +| existing | main | yes | -843.3 | 995.1 | 759.5 | yes | yes | no | 1.37 | 2.09 | 2.07 | | +| existing | casl | yes | -627.2 | 1335.0 | 852.1 | yes | yes | no | 0.25 | 1.62 | 2.18 | | +| existing | inst | yes | -150.7 | 191.4 | 133.3 | yes | yes | no | 1.42 | 2.87 | 6.69 | | +| existing_floor | main | yes | -843.3 | 995.1 | 759.5 | yes | yes | no | 1.37 | 2.09 | 2.07 | | +| existing_floor | casl | yes | -627.2 | 1335.0 | 852.1 | yes | yes | no | 0.25 | 1.62 | 2.18 | | +| existing_floor | inst | yes | -150.7 | 191.4 | 133.3 | yes | yes | no | 1.42 | 2.87 | 6.69 | | +| joint | main | yes | -569.6 | 844.8 | 620.3 | yes | yes | no | 0.23 | 1.81 | 2.08 | | +| joint | casl | yes | -906.5 | 1437.4 | 990.5 | yes | yes | no | 0.24 | 1.98 | 2.68 | | +| joint | inst | yes | -109.7 | 171.3 | 109.8 | yes | yes | no | 2.08 | 2.74 | 5.38 | | + +Joint arm race mixing matrices (rows = ego race). Reported: expected partnership endpoints by ego race and partner race at the synthetic population, standardized on the ego side only; row sums are `nodefactor_race`. Balanced: expected edge counts after symmetric scaling with `nodefactor_race` held fixed; the diagonal is `nodematch_race`. + +**main, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 504.9 | 17.9 | 177.5 | +| Hispanic | 13.6 | 24.0 | 56.2 | +| White.Other | 124.8 | 62.8 | 707.9 | + +**main, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 266.0 | 15.5 | 152.8 | +| Hispanic | 15.5 | 11.0 | 56.2 | +| White.Other | 152.8 | 56.2 | 343.2 | + +**casl, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 1050.6 | 92.3 | 378.5 | +| Hispanic | 27.7 | 43.6 | 77.9 | +| White.Other | 228.9 | 170.8 | 804.6 | + +**casl, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 579.9 | 43.3 | 318.3 | +| Hispanic | 43.3 | 10.3 | 85.2 | +| White.Other | 318.3 | 85.2 | 400.3 | + +**inst, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 77.7 | 16.7 | 59.7 | +| Hispanic | 2.7 | 4.6 | 9.7 | +| White.Other | 22.3 | 27.1 | 122.1 | + +**inst, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 51.2 | 6.0 | 45.8 | +| Hispanic | 6.0 | 0.7 | 9.7 | +| White.Other | 45.8 | 9.7 | 58.0 | + +### boston_default + +| Arm | Layer | Feasible | Worst slack | Edges | Same-race target | Fit | Converged | Lin. dep. warning | Minutes | netdx max abs Z | netdx max abs % | Note | +|---|---|---|---:|---:|---:|---|---|---|---:|---:|---:|---| +| existing | main | yes | -491.1 | 923.4 | 679.5 | yes | yes | no | 0.24 | 2.13 | 2.70 | | +| existing | casl | yes | -639.0 | 1418.9 | 907.7 | yes | yes | no | 0.21 | 1.68 | 2.14 | | +| existing | inst | yes | -91.4 | 249.5 | 143.1 | yes | yes | no | 1.37 | 3.05 | 5.15 | | +| existing_floor | main | yes | -491.1 | 923.4 | 679.5 | yes | yes | no | 0.23 | 2.13 | 2.70 | | +| existing_floor | casl | yes | -639.0 | 1418.9 | 907.7 | yes | yes | no | 0.21 | 1.68 | 2.14 | | +| existing_floor | inst | yes | -91.4 | 249.5 | 143.1 | yes | yes | no | 1.38 | 3.05 | 5.15 | | +| joint | main | yes | -325.8 | 839.7 | 566.4 | yes | yes | no | 0.23 | 1.64 | 1.35 | | +| joint | casl | yes | -724.8 | 1587.8 | 949.5 | yes | yes | no | 0.20 | 3.19 | 3.66 | | +| joint | inst | yes | -79.9 | 243.8 | 142.7 | yes | yes | no | 1.61 | 2.43 | 2.90 | | + +Joint arm race mixing matrices (rows = ego race). Reported: expected partnership endpoints by ego race and partner race at the synthetic population, standardized on the ego side only; row sums are `nodefactor_race`. Balanced: expected edge counts after symmetric scaling with `nodefactor_race` held fixed; the diagonal is `nodematch_race`. + +**main, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 150.1 | 13.4 | 85.3 | +| Hispanic | 26.2 | 132.3 | 191.7 | +| White.Other | 88.5 | 122.5 | 869.3 | + +**main, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 72.2 | 21.3 | 83.2 | +| Hispanic | 21.3 | 80.1 | 168.7 | +| White.Other | 83.2 | 168.7 | 414.2 | + +**casl, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 346.6 | 88.2 | 259.0 | +| Hispanic | 54.9 | 279.9 | 334.6 | +| White.Other | 158.8 | 390.4 | 1263.3 | + +**casl, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 199.1 | 73.0 | 222.5 | +| Hispanic | 73.0 | 126.8 | 342.8 | +| White.Other | 222.5 | 342.8 | 623.6 | + +**inst, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 36.0 | 11.7 | 33.8 | +| Hispanic | 12.3 | 32.6 | 54.7 | +| White.Other | 31.4 | 59.3 | 215.9 | + +**inst, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 18.3 | 11.9 | 33.0 | +| Hispanic | 11.9 | 15.8 | 56.2 | +| White.Other | 33.0 | 56.2 | 108.7 | + +### chicago_default + +| Arm | Layer | Feasible | Worst slack | Edges | Same-race target | Fit | Converged | Lin. dep. warning | Minutes | netdx max abs Z | netdx max abs % | Note | +|---|---|---|---:|---:|---:|---|---|---|---:|---:|---:|---| +| existing | main | yes | -741.5 | 1008.3 | 592.5 | yes | yes | no | 1.05 | 3.13 | 4.11 | | +| existing | casl | yes | -1355.5 | 1726.5 | 1059.9 | yes | yes | no | 0.24 | 2.35 | 1.34 | | +| existing | inst | yes | -259.3 | 319.8 | 212.8 | yes | yes | no | 2.29 | 3.67 | 3.21 | | +| existing_floor | main | yes | -741.5 | 1008.3 | 592.5 | yes | yes | no | 1.03 | 3.13 | 4.11 | | +| existing_floor | casl | yes | -1355.5 | 1726.5 | 1059.9 | yes | yes | no | 0.25 | 2.35 | 1.34 | | +| existing_floor | inst | yes | -259.3 | 319.8 | 212.8 | yes | yes | no | 2.29 | 3.67 | 3.21 | | +| joint | main | yes | -637.0 | 907.2 | 572.3 | yes | yes | no | 0.22 | 2.15 | 2.08 | | +| joint | casl | yes | -1405.1 | 1857.6 | 1003.5 | yes | yes | no | 0.24 | 1.54 | 2.73 | | +| joint | inst | yes | -178.2 | 301.0 | 152.0 | yes | yes | no | 1.65 | 3.40 | 4.57 | | + +Joint arm race mixing matrices (rows = ego race). Reported: expected partnership endpoints by ego race and partner race at the synthetic population, standardized on the ego side only; row sums are `nodefactor_race`. Balanced: expected edge counts after symmetric scaling with `nodefactor_race` held fixed; the diagonal is `nodematch_race`. + +**main, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 248.3 | 30.4 | 112.4 | +| Hispanic | 43.3 | 293.4 | 244.0 | +| White.Other | 81.1 | 150.9 | 610.5 | + +**main, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 128.7 | 40.3 | 93.4 | +| Hispanic | 40.3 | 169.6 | 201.1 | +| White.Other | 93.4 | 201.1 | 274.0 | + +**casl, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 472.8 | 170.3 | 459.5 | +| Hispanic | 69.4 | 516.7 | 570.6 | +| White.Other | 95.7 | 343.0 | 1017.3 | + +**casl, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 333.6 | 151.7 | 283.8 | +| Hispanic | 151.7 | 293.1 | 418.8 | +| White.Other | 283.8 | 418.8 | 376.8 | + +**inst, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 73.3 | 15.1 | 54.6 | +| Hispanic | 30.0 | 50.1 | 104.3 | +| White.Other | 36.2 | 43.2 | 195.4 | + +**inst, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 37.2 | 27.9 | 40.6 | +| Hispanic | 27.9 | 37.9 | 80.5 | +| White.Other | 40.6 | 80.5 | 76.9 | + +### seattle_default + +| Arm | Layer | Feasible | Worst slack | Edges | Same-race target | Fit | Converged | Lin. dep. warning | Minutes | netdx max abs Z | netdx max abs % | Note | +|---|---|---|---:|---:|---:|---|---|---|---:|---:|---:|---| +| existing | main | no | 115.8 | 1071.4 | 726.4 | no | no | n/a | n/a | n/a | n/a | skipped: infeasible race targets | +| existing | casl | no | 155.4 | 1575.6 | 996.2 | no | no | n/a | n/a | n/a | n/a | skipped: infeasible race targets | +| existing | inst | no | 40.3 | 231.9 | 142.5 | no | no | n/a | n/a | n/a | n/a | skipped: infeasible race targets | +| existing_floor | main | yes | -16.2 | 1071.4 | 858.3 | yes | yes | no | 1.39 | 3.26 | 2.19 | | +| existing_floor | casl | yes | -44.4 | 1575.6 | 1196.0 | yes | yes | no | 0.24 | 2.17 | 2.89 | | +| existing_floor | inst | yes | -2.3 | 231.9 | 185.1 | yes | yes | no | 2.03 | 4.45 | 5.63 | | +| joint | main | yes | -56.8 | 1034.8 | 871.4 | yes | yes | no | 0.91 | 4.40 | 5.22 | | +| joint | casl | yes | -106.1 | 1607.1 | 1258.1 | yes | yes | no | 0.21 | 2.09 | 2.08 | | +| joint | inst | yes | -9.1 | 221.4 | 176.8 | yes | yes | no | 1.63 | 3.78 | 6.30 | | + +Joint arm race mixing matrices (rows = ego race). Reported: expected partnership endpoints by ego race and partner race at the synthetic population, standardized on the ego side only; row sums are `nodefactor_race`. Balanced: expected edge counts after symmetric scaling with `nodefactor_race` held fixed; the diagonal is `nodematch_race`. + +**main, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 45.4 | 3.6 | 36.2 | +| Hispanic | 8.9 | 37.5 | 88.6 | +| White.Other | 115.4 | 139.0 | 1595.0 | + +**main, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 12.1 | 4.0 | 57.0 | +| Hispanic | 4.0 | 14.3 | 102.4 | +| White.Other | 57.0 | 102.4 | 845.0 | + +**casl, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 98.8 | 23.6 | 88.0 | +| Hispanic | 18.9 | 87.0 | 138.7 | +| White.Other | 219.9 | 487.8 | 2051.5 | + +**casl, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 32.4 | 10.3 | 135.4 | +| Hispanic | 10.3 | 15.5 | 203.3 | +| White.Other | 135.4 | 203.3 | 1210.2 | + +**inst, reported endpoints** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 5.5 | 4.2 | 11.9 | +| Hispanic | 1.9 | 11.2 | 19.1 | +| White.Other | 18.1 | 81.4 | 289.6 | + +**inst, balanced edges** + +| | Black | Hispanic | White.Other | +|---|---:|---:|---:| +| Black | 2.4 | 1.4 | 15.3 | +| Hispanic | 1.4 | 1.4 | 27.8 | +| White.Other | 15.3 | 27.8 | 173.0 | + +## Reproducibility + +```r +source(system.file("validation/method_comparison.R", package = "ARTnet")) +fits <- fit_comparison_ergms(ncores = 4) +render_fit_report(fits) +``` From 4d373989b6f987ae9cceb8023cd68b567b732ac5 Mon Sep 17 00:00:00 2001 From: Samuel Jenness Date: Sat, 7 Nov 2026 09:41:54 -0500 Subject: [PATCH 3/4] ARTnet 2.10.0 --- DESCRIPTION | 4 ++-- NEWS.md | 9 +++++++++ 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0e40118..f768a8c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ARTnet -Version: 2.9.3 -Date: 2026-10-18 +Version: 2.10.0 +Date: 2026-11-07 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 bda67bd..f337377 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,12 @@ +## ARTnet v2.10.0 + +* `build_netstats(method = "joint")` now balances race mixing to the target population: partner race is modeled with a multinomial (`netparams$$joint_prace_model`), the reported mixing matrix is built at the synthetic population and symmetrically scaled with `nodefactor_race` held fixed, and `nodematch_race` and `nodematch_race_diffF` are read from the balanced matrix. The matrices are stored as `mixmat_race` and `mixmat_race_reported` (#9). +* `build_netstats()` gains `nodematch.race.floor` (default `FALSE`), which raises each layer's `nodematch_race_diffF` to the random-mixing value where it falls below, the opt-in stopgap for `method = "existing"` (#9). +* `build_netstats()` now stops when a layer's race targets are infeasible (`nodefactor_race[r] > edges + nodematch_race_diffF`), naming the layer, race and excess; `check.feasibility = FALSE` returns the targets anyway. The check is exported as `check_race_feasibility()` (#10). +* `build_netparams()` gains `seed`, which makes the partner-race imputation reproducible (#11). +* `inst/validation/method_comparison.R` gains `fit_comparison_ergms()`, which fits the three template ERGMs with MCMLE for every scenario city and method after checking feasibility, and `render_fit_report()` (#12). +* `nnet` is a new import. + ## ARTnet v2.9.3 * `trim_epistats()` now also trims `hiv.mod`, which is fit when `init.hiv.prev` is `NULL`. From d09a4b63819f38325a6147d5a49656a9602b5028 Mon Sep 17 00:00:00 2001 From: Samuel Jenness Date: Sat, 7 Nov 2026 10:01:55 -0500 Subject: [PATCH 4/4] Add duration.cap to build_netstats() and warn when joint_lm drops an edited by-age duration Under duration.method = "joint_lm", build_netstats(method = "joint") recomputes the by-age durations at the synthetic population and used the netparams values only as a fallback, so a value edited in netparams$$durs..byage$mean.dur.adj after build_netparams() was dropped without a word. build_netstats() now takes duration.cap, applied with pmin() to every main and casual duration passed to dissolution_coefs() under both methods, and build_netparams() records the by-age durations it produced under joint_lm so build_netstats() can warn when the column has been changed. --- NEWS.md | 2 + R/NetParams.R | 21 ++++++- R/NetStats.R | 99 +++++++++++++++++++++++++++--- man/build_netparams.Rd | 7 ++- man/build_netstats.Rd | 13 ++++ tests/testthat/test-duration-cap.R | 96 +++++++++++++++++++++++++++++ 6 files changed, 226 insertions(+), 12 deletions(-) create mode 100644 tests/testthat/test-duration-cap.R diff --git a/NEWS.md b/NEWS.md index f337377..f37b1f3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,8 @@ * `build_netstats()` now stops when a layer's race targets are infeasible (`nodefactor_race[r] > edges + nodematch_race_diffF`), naming the layer, race and excess; `check.feasibility = FALSE` returns the targets anyway. The check is exported as `check_race_feasibility()` (#10). * `build_netparams()` gains `seed`, which makes the partner-race imputation reproducible (#11). * `inst/validation/method_comparison.R` gains `fit_comparison_ergms()`, which fits the three template ERGMs with MCMLE for every scenario city and method after checking feasibility, and `render_fit_report()` (#12). +* `build_netstats()` gains `duration.cap`, which caps every main and casual duration passed to `dissolution_coefs()` under both methods, with a message naming the strata it binds on. +* Under `duration.method = "joint_lm"`, `build_netstats(method = "joint")` now warns when `netparams$$durs..byage$mean.dur.adj` has been edited after `build_netparams()`, since the by-age durations are recomputed at the synthetic population and the edit was silently dropped before. * `nnet` is a new import. ## ARTnet v2.9.3 diff --git a/R/NetParams.R b/R/NetParams.R index 322b123..8b4aa02 100644 --- a/R/NetParams.R +++ b/R/NetParams.R @@ -199,7 +199,12 @@ fit_joint_multinom <- function(d, response, main_terms, race, geog.lvl) { #' output can't be faithfully honored by a geometric TERGM without a non-standard #' age-dependent dissolution term, making it a diagnostic rather than a #' production target. The one-off layer has no durations, so this flag only -#' affects main and casual-layer output. +#' affects main and casual-layer output. Note that under `"joint_lm"`, +#' [`build_netstats`] with `method = "joint"` recomputes the by-age durations at the +#' synthetic population and uses `durs..byage$mean.dur.adj` only as a fallback +#' for empty strata, so editing that column after `build_netparams()` has no effect +#' there (a warning says so). To cap durations, use the `duration.cap` argument of +#' [`build_netstats`], which applies under both methods. #' @param method Character. Either `"existing"` (default) or `"joint"`. `"existing"` reproduces #' the pre-refactor behavior byte-for-byte: a separate univariate Poisson/binomial/linear #' fit for each ERGM target statistic. `"joint"` leaves all of those outputs intact **and** @@ -810,6 +815,15 @@ build_netparams <- function(epistats, out$main$durs.main.byage <- rbind(out$main$durs.main.byage, df) } + # Under joint_lm, build_netstats() recomputes these by-age durations at the + # synthetic population and uses the values here only as a fallback, so an + # edit made to them afterwards would be lost. Record what was produced so + # build_netstats() can warn when the column has been changed. + if (duration.method == "joint_lm") { + attr(out$main$durs.main.byage, "mean.dur.adj.built") <- + out$main$durs.main.byage$mean.dur.adj + } + ## joint g-computation models (additive outputs; see issues #61/#62/#63) ---- if (method == "joint") { @@ -1237,6 +1251,11 @@ build_netparams <- function(epistats, out$casl$durs.casl.byage <- rbind(out$casl$durs.casl.byage, df) } + if (duration.method == "joint_lm") { + attr(out$casl$durs.casl.byage, "mean.dur.adj.built") <- + out$casl$durs.casl.byage$mean.dur.adj + } + ## joint g-computation models (additive outputs; see issues #61/#62/#63) ---- if (method == "joint") { diff --git a/R/NetStats.R b/R/NetStats.R index 85b6748..4cb70fc 100644 --- a/R/NetStats.R +++ b/R/NetStats.R @@ -170,6 +170,47 @@ } +# Apply the optional build_netstats(duration.cap) to a vector of durations +# (in time units) right before dissolution_coefs(), with a message naming +# the layer and the strata it binds on. NULL cap returns the input as is. +.cap_durations <- function(x, cap, layer, what) { + if (is.null(cap)) return(x) + over <- which(!is.na(x) & x > cap) + if (length(over) > 0) { + message(sprintf( + "duration.cap: capping %d of %d %s %s duration%s at %s (%s: %s).", + length(over), length(x), layer, what, if (length(x) > 1) "s" else "", + format(cap), + if (what == "by-age") paste0("stratum ", paste(over, collapse = ", ")) else "value", + paste(sprintf("%.1f", x[over]), collapse = ", "))) + x[over] <- cap + } + x +} + + +# Under duration.method = "joint_lm", build_netstats() recomputes the by-age +# durations at the synthetic population and uses the netparams values only +# as a fallback for empty strata, so an edit to +# netparams$$durs..byage$mean.dur.adj made after +# build_netparams() has no effect. build_netparams() records the values it +# produced as an attribute; warn when the column no longer matches. +.warn_if_byage_edited <- function(byage, layer) { + built <- attr(byage, "mean.dur.adj.built") + if (is.null(built)) return(invisible(FALSE)) + if (!isTRUE(all.equal(as.numeric(byage$mean.dur.adj), as.numeric(built)))) { + warning("netparams$", layer, "$durs.", layer, ".byage$mean.dur.adj differs from ", + "what build_netparams() produced. Under duration.method = 'joint_lm', ", + "build_netstats() recomputes the by-age durations at the synthetic ", + "population and uses the netparams values only as a fallback for empty ", + "strata, so this edit is ignored. To cap durations, use the ", + "duration.cap argument of build_netstats().", call. = FALSE) + return(invisible(TRUE)) + } + invisible(FALSE) +} + + #' Calculate Network Target Statistics #' #' @description Calculates the final target statistics for the network models by applying @@ -180,6 +221,17 @@ #' @param netparams Output from [`build_netparams`]. #' @param network.size Size of the starting network. #' @param expect.mort Expected average mortality level to pass into [`dissolution_coefs`] function. +#' @param duration.cap Optional single number, in the same time units as the durations +#' (`time.unit` days). When supplied, every main and casual partnership duration passed +#' to [`dissolution_coefs`] (the homogeneous and the by-age values, under both methods) +#' is capped at this value with `pmin()`, and a message names the layer and strata it +#' binds on. [`dissolution_coefs`] rejects a duration at or above the departure ceiling, +#' about `1 / (2 * expect.mort)`, so a cap just below that keeps every stratum fittable +#' when a small ARTnet sub-sample produces a very long stratum duration. Default `NULL` +#' (no cap). This replaces editing `netparams$$durs..byage$mean.dur.adj` +#' by hand, which `method = "joint"` with `duration.method = "joint_lm"` does not honor +#' (the by-age durations are recomputed at the synthetic population; a warning is +#' raised if the column has been changed). #' @param age.pyramid Numerical vector of length equal to the length of the age range specified in #' [`build_epistats`], containing probability distribution of each year of age, summing to #' one. If `NULL`, then a default distribution based on 2020 NCHS data is used. @@ -347,9 +399,18 @@ build_netstats <- function(epistats, netparams, seed = NULL, nodematch.race.floor = FALSE, check.feasibility = TRUE, + duration.cap = NULL, browser = FALSE) { method <- match.arg(method) + if (!is.null(duration.cap)) { + if (!is.numeric(duration.cap) || length(duration.cap) != 1 || + is.na(duration.cap) || duration.cap < 1) { + stop("`duration.cap` must be a single number >= 1, in the same time ", + "units as the partnership durations.", call. = FALSE) + } + } + # Optional reproducibility. The default sampling path draws the synthetic # population stochastically (sample / runif / rbinom / apportion_lr), so two # calls with identical arguments otherwise differ at the Monte Carlo noise @@ -871,6 +932,14 @@ build_netstats <- function(epistats, netparams, sex_cess_extra_row = sex.cess.mod, fallback_mean_dur_adj = netparams$casl$durs.casl.byage$mean.dur.adj ) + # The override above drops any edit made to the netparams by-age + # durations after build_netparams(); say so rather than stay silent. + if (!is.null(synth_dur_main_byage)) { + .warn_if_byage_edited(netparams$main$durs.main.byage, "main") + } + if (!is.null(synth_dur_casl_byage)) { + .warn_if_byage_edited(netparams$casl$durs.casl.byage, "casl") + } } @@ -995,13 +1064,18 @@ build_netstats <- function(epistats, netparams, # ARTnet aggregation from build_netparams; it is not consumed by # EpiModelHIV-Template's tergm offset, and its synth analog can be added # later without changing the byage interface that matters for production. - out$main$diss.homog <- dissolution_coefs(dissolution = ~offset(edges), - duration = netparams$main$durs.main.homog$mean.dur.adj, - d.rate = expect.mort) + out$main$diss.homog <- dissolution_coefs( + dissolution = ~offset(edges), + duration = .cap_durations(netparams$main$durs.main.homog$mean.dur.adj, + duration.cap, "main", "homogeneous"), + d.rate = expect.mort + ) out$main$diss.byage <- dissolution_coefs( dissolution = ~offset(edges) + offset(nodematch("age.grp", diff = TRUE)), - duration = if (!is.null(synth_dur_main_byage)) synth_dur_main_byage - else netparams$main$durs.main.byage$mean.dur.adj, + duration = .cap_durations( + if (!is.null(synth_dur_main_byage)) synth_dur_main_byage + else netparams$main$durs.main.byage$mean.dur.adj, + duration.cap, "main", "by-age"), d.rate = expect.mort ) @@ -1109,13 +1183,18 @@ build_netstats <- function(epistats, netparams, } # Dissolution (see note on diss.byage / diss.homog at the main layer block) - out$casl$diss.homog <- dissolution_coefs(dissolution = ~offset(edges), - duration = netparams$casl$durs.casl.homog$mean.dur.adj, - d.rate = expect.mort) + out$casl$diss.homog <- dissolution_coefs( + dissolution = ~offset(edges), + duration = .cap_durations(netparams$casl$durs.casl.homog$mean.dur.adj, + duration.cap, "casl", "homogeneous"), + d.rate = expect.mort + ) out$casl$diss.byage <- dissolution_coefs( dissolution = ~offset(edges) + offset(nodematch("age.grp", diff = TRUE)), - duration = if (!is.null(synth_dur_casl_byage)) synth_dur_casl_byage - else netparams$casl$durs.casl.byage$mean.dur.adj, + duration = .cap_durations( + if (!is.null(synth_dur_casl_byage)) synth_dur_casl_byage + else netparams$casl$durs.casl.byage$mean.dur.adj, + duration.cap, "casl", "by-age"), d.rate = expect.mort ) diff --git a/man/build_netparams.Rd b/man/build_netparams.Rd index bad1c95..1485028 100644 --- a/man/build_netparams.Rd +++ b/man/build_netparams.Rd @@ -52,7 +52,12 @@ underlying mean full-partnership duration was explored and removed: its output can't be faithfully honored by a geometric TERGM without a non-standard age-dependent dissolution term, making it a diagnostic rather than a production target. The one-off layer has no durations, so this flag only -affects main and casual-layer output.} +affects main and casual-layer output. Note that under \code{"joint_lm"}, +\code{\link{build_netstats}} with \code{method = "joint"} recomputes the by-age durations at the +synthetic population and uses \verb{durs..byage$mean.dur.adj} only as a fallback +for empty strata, so editing that column after \code{build_netparams()} has no effect +there (a warning says so). To cap durations, use the \code{duration.cap} argument of +\code{\link{build_netstats}}, which applies under both methods.} \item{method}{Character. Either \code{"existing"} (default) or \code{"joint"}. \code{"existing"} reproduces the pre-refactor behavior byte-for-byte: a separate univariate Poisson/binomial/linear diff --git a/man/build_netstats.Rd b/man/build_netstats.Rd index e197be8..96ca6c8 100644 --- a/man/build_netstats.Rd +++ b/man/build_netstats.Rd @@ -18,6 +18,7 @@ build_netstats( seed = NULL, nodematch.race.floor = FALSE, check.feasibility = TRUE, + duration.cap = NULL, browser = FALSE ) } @@ -139,6 +140,18 @@ race whose \code{nodefactor_race} exceeds \code{edges + nodematch_race_diffF}, n the race and the excess. Such targets cannot be fit by MCMLE. Set to \code{FALSE} to return the targets anyway.} +\item{duration.cap}{Optional single number, in the same time units as the durations +(\code{time.unit} days). When supplied, every main and casual partnership duration passed +to \code{\link[EpiModel:dissolution_coefs]{dissolution_coefs}} (the homogeneous and the by-age values, under both methods) +is capped at this value with \code{pmin()}, and a message names the layer and strata it +binds on. \code{\link[EpiModel:dissolution_coefs]{dissolution_coefs}} rejects a duration at or above the departure ceiling, +about \code{1 / (2 * expect.mort)}, so a cap just below that keeps every stratum fittable +when a small ARTnet sub-sample produces a very long stratum duration. Default \code{NULL} +(no cap). This replaces editing \verb{netparams$$durs..byage$mean.dur.adj} +by hand, which \code{method = "joint"} with \code{duration.method = "joint_lm"} does not honor +(the by-age durations are recomputed at the synthetic population; a warning is +raised if the column has been changed).} + \item{browser}{If \code{TRUE}, run \code{build_netparams} in interactive browser mode.} } \description{ diff --git a/tests/testthat/test-duration-cap.R b/tests/testthat/test-duration-cap.R new file mode 100644 index 0000000..7450bb6 --- /dev/null +++ b/tests/testthat/test-duration-cap.R @@ -0,0 +1,96 @@ +# Tests for build_netstats(duration.cap) and the warning that a by-age +# duration edited in netparams is ignored under duration.method = "joint_lm". + +skip_without_artnetdata <- function() { + testthat::skip_if(system.file(package = "ARTnetData") == "", + "ARTnetData not installed") +} + +cap_env <- new.env(parent = emptyenv()) + +get_inputs <- function() { + if (is.null(cap_env$epi)) { + set.seed(20260419L) + cap_env$epi <- build_epistats( + geog.lvl = "city", geog.cat = "Atlanta", + init.hiv.prev = c(0.33, 0.137, 0.084), + race = TRUE, time.unit = 7 + ) + cap_env$np_e <- build_netparams(cap_env$epi, smooth.main.dur = TRUE, seed = 1L) + cap_env$np_j <- build_netparams(cap_env$epi, smooth.main.dur = TRUE, seed = 1L, + method = "joint", duration.method = "joint_lm") + } + cap_env +} + +build <- function(np, method, ...) { + build_netstats(get_inputs()$epi, np, network.size = 2000, expect.mort = 0.00025, + method = method, seed = 1L, check.feasibility = FALSE, ...) +} + +test_that("duration.cap caps the by-age and homogeneous durations under both methods", { + skip_without_artnetdata() + inp <- get_inputs() + cap <- 100 + for (method in c("existing", "joint")) { + np <- if (method == "joint") inp$np_j else inp$np_e + ns0 <- suppressWarnings(build(np, method)) + expect_true(any(ns0$main$diss.byage$duration > cap)) # the cap binds + expect_message(ns1 <- suppressWarnings(build(np, method, duration.cap = cap)), + "duration.cap: capping") + for (layer in c("main", "casl")) { + for (diss in c("diss.byage", "diss.homog")) { + expect_true(all(ns1[[layer]][[diss]]$duration <= cap + 1e-9), + info = paste(method, layer, diss)) + expect_equal(ns1[[layer]][[diss]]$duration, + pmin(ns0[[layer]][[diss]]$duration, cap), + info = paste(method, layer, diss)) + } + } + # a cap that does not bind changes nothing + ns2 <- suppressWarnings(build(np, method, duration.cap = 1e6)) + for (layer in c("main", "casl")) { + expect_equal(ns2[[layer]]$diss.byage$duration, ns0[[layer]]$diss.byage$duration) + expect_equal(ns2[[layer]]$diss.byage$coef.adj, ns0[[layer]]$diss.byage$coef.adj) + expect_equal(ns2[[layer]]$diss.homog$duration, ns0[[layer]]$diss.homog$duration) + } + } +}) + +test_that("duration.cap is validated", { + skip_without_artnetdata() + inp <- get_inputs() + expect_error(build(inp$np_e, "existing", duration.cap = 0), "duration.cap") + expect_error(build(inp$np_e, "existing", duration.cap = "a"), "duration.cap") + expect_error(build(inp$np_e, "existing", duration.cap = c(10, 20)), "duration.cap") +}) + +test_that("build_netparams records the joint_lm by-age durations it produced", { + skip_without_artnetdata() + inp <- get_inputs() + for (layer in c("main", "casl")) { + df <- inp$np_j[[layer]][[paste0("durs.", layer, ".byage")]] + expect_equal(attr(df, "mean.dur.adj.built"), df$mean.dur.adj) + expect_null(attr(inp$np_e[[layer]][[paste0("durs.", layer, ".byage")]], + "mean.dur.adj.built")) + } +}) + +test_that("an edited by-age duration is honored under empirical and warned about under joint_lm", { + skip_without_artnetdata() + inp <- get_inputs() + # legacy: the edit flows through to the dissolution coefficients + np_e2 <- inp$np_e + np_e2$main$durs.main.byage$mean.dur.adj[1] <- 50 + expect_no_warning(ns_e <- build(np_e2, "existing")) + expect_equal(ns_e$main$diss.byage$duration[1], 50) + # joint_lm: the synth aggregation replaces the column, so the edit is lost; + # build_netstats says so + ns_j0 <- expect_no_warning(build(inp$np_j, "joint")) + np_j2 <- inp$np_j + np_j2$main$durs.main.byage$mean.dur.adj[1] <- 50 + expect_warning(ns_j <- build(np_j2, "joint"), + "durs.main.byage\\$mean.dur.adj differs from what build_netparams") + expect_equal(ns_j$main$diss.byage$duration, ns_j0$main$diss.byage$duration) + expect_false(isTRUE(all.equal(ns_j$main$diss.byage$duration[1], 50))) +})