Balance race mixing to the target population, check feasibility, seed the imputation, fit every city with MCMLE (#9 to #12)

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 <samuel.m.jenness@emory.edu>
@@ -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

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)

NEWS.md

@@ -1,3 +1,14 @@
+## ARTnet v2.10.0
+
+* `build_netstats(method = "joint")` now balances race mixing to the target population: partner race is modeled with a multinomial (`netparams$<layer>$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).
+* `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$<layer>$durs.<layer>.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
 
 * `trim_epistats()` now also trims `hiv.mod`, which is fit when `init.hiv.prev` is `NULL`.

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
@@ -157,7 +199,12 @@ fit_joint_glm <- function(d, response, main_terms,
 #'        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.<layer>.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**
@@ -169,12 +216,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 `$<layer>$joint_nm_age_model`, binomial
 #'            for same-race (when `race = TRUE`) at `$<layer>$joint_nm_race_model`,
-#'            and Gaussian fits for the age absolute-difference terms at
-#'            `$<layer>$joint_absdiff_age_model` and
+#'            a multinomial (`nnet::multinom`) for partner race (when `race = TRUE`) at
+#'            `$<layer>$joint_prace_model`, and Gaussian fits for the age
+#'            absolute-difference terms at `$<layer>$joint_absdiff_age_model` and
 #'            `$<layer>$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 +274,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)
 
@@ -756,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") {
@@ -800,6 +868,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,
@@ -1174,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") {
@@ -1211,6 +1293,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 +1614,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",