Duration methods: empirical + joint_lm (#63 phase 3)#71
Open
Duration methods: empirical + joint_lm (#63 phase 3)#71
Conversation
… phase 3)
New `duration.method` argument to `build_netparams()`. Default
`"empirical"` preserves byte-identical behavior. Two alternatives:
- `"weibull_strat"`: per-stratum Weibull AFT fits via survival::survreg
with proper right-censoring (ongoing partnerships = censored, completed
= observed events). Uses more of the data than empirical (which only
uses ongoing). Attaches the Weibull shape parameter k as a
`"weibull_shape"` attribute on `durs.<layer>.byage` as a diagnostic
for the constant-hazard assumption (k == 1 means geometric; the
current ARTnet data gives k ~= 0.5 in every stratum, a clean
rejection of the geometric assumption).
- `"joint_lm"`: log-linear regression `lm(log(duration.time) ~
joint ego + partner + matching terms)` on ongoing partnerships.
Stratum medians computed from model predictions. Fitted model stored
at `netparams$<layer>$joint_duration_model` for use by future
build_netstats refactors that want per-dyad predictions.
All three methods produce `durs.<layer>.{homog,byage}` data.frames with
identical shapes and column names so the geometric transformation
(`mean.dur.adj = 1/(1 - 2^(-1/median))`) and downstream
`dissolution_coefs()` pipeline are unchanged. TERGM parameterization
identical regardless of method.
Implementation is a minimal surgical override at the end of each
layer's empirical duration block in NetParams.R. When duration.method
!= "empirical", the stratum median / mean values are replaced with
model-based estimates and rates.*.adj / mean.dur.adj are recomputed.
Smoothing (smooth.main.dur) and sex.cess.mod logic apply uniformly
after the override. Per-stratum fits fall back to empirical when the
new method fails (sparse data, convergence issues).
New dependency: survival (base-R recommended package) added to
Suggests. The weibull_strat path guards with a requireNamespace
check and errors clearly if survival isn't installed.
Findings from Atlanta + race = TRUE (empirical vs weibull_strat vs
joint_lm on mean.dur.adj):
Main layer match.grp.1: 73 | 379 | 76
Main layer match.grp.5: 934 | 1,482,404 | 446
Casl layer nonmatch: 106 | 359 | 93
Casl layer match.grp.5: 150 | 1,157 | 129
Weibull shape (main): 0.48, 0.72, 0.60, 0.83, 0.35, 0.46 -- all well
below 1, consistent with decreasing hazard. The extreme Weibull
mean.dur.adj values in high-age-matched strata reflect extrapolation
under k << 1 on heavily censored data; the roxygen doc flags this
explicitly and recommends joint_lm as the more production-safe
non-default option.
New tests: tests/testthat/test-duration-method.R (8 blocks, 112
assertions). Covers: default is empirical; weibull_strat produces valid
durs shape + shape diagnostic attribute; joint_lm stores lm model;
output shape invariant across all three methods (for dissolution_coefs
compatibility); non-duration netparams fields preserved across
methods; weibull shapes are all < 1 (substantive finding locked in);
joint_lm composes cleanly with method = "joint" in build_netstats.
Backward-compat snapshot: default and explicit method = "existing"
match 3/3 on all parameter sets.
Full test suite: 563 assertions pass in 10.3s.
R CMD check: 0 errors / 0 warnings / 0 relevant notes
(one "unable to verify current time" environment note, unrelated).
Part of #63.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Replaces the naive survreg-based Weibull fit with a length-biased MLE. ARTnet is a cross-sectional prevalence sample, not an incidence cohort: duration.time for ongoing partnerships is the backward recurrence time (age at inspection in a renewal process), not the full partnership duration. Under stationarity and Weibull(k, lambda) durations, the density of observed ages is: f_A(a) = S(a; k, lambda) / mu(k, lambda) where S is the Weibull survival function and mu = lambda * Gamma(1 + 1/k) is the Weibull mean. Maximizing sum log f_A(t_i) over ongoing obs via optim() gives unbiased estimates of (k, lambda) of the underlying full-duration distribution. The naive survreg treatment of the same data interprets length-biased elapsed durations as incidence-cohort survival times and produces catastrophically biased estimates (scale parameter off by 3+ orders of magnitude on ARTnet, median extrapolations to 15-20 year main partnerships in the oldest matched age groups, etc.). Completed partnerships (ongoing2 == 0) are deliberately not used here: they are subject to a different truncation (past-12-month recall window) that would require separate modeling. Restricting to ongoing matches the convention of empirical and joint_lm. Additionally, under weibull_strat, mean.dur.adj is now set directly to the Weibull analytical mean (wt * lambda * Gamma(1 + 1/k)) rather than being derived from median.dur via the geometric-distribution formula. When k is far from 1, these diverge by orders of magnitude; the geometric formula is inappropriate and was masking the underlying fit quality. Empirical comparison on Atlanta + race = TRUE, mean.dur.adj: MAIN empirical weibull_lb joint_lm nonmatch 208.5 160.0 235.5 matched.1 72.6 72.2 75.7 matched.5 934.4 1368.4 445.9 (was 1,482,404 under naive) CASL empirical weibull_lb joint_lm nonmatch 105.8 64.4 93.3 matched.5 150.1 96.6 129.1 (was 1,157,200 under naive) All three methods now produce values in the same order of magnitude. Weibull shape parameters are also more informative: MAIN k per stratum: 0.68 1.03 1.96 2.62 2.84 4.96 (overall 0.63) CASL k per stratum: 0.62 0.68 0.55 0.59 0.53 0.59 (overall 0.57) Casual partnerships show consistent decreasing hazard (k < 1) across all strata. Main shows a more interesting pattern: non-matched and youngest-matched near k = 1 (geometric-like), older-matched strata with k > 1 (increasing hazard -- stable mature relationships with break points). Issue #72 opened for the broader length-bias and 5-partnership truncation issue affecting formation stats. Test: "all shapes < 1" assertion replaced with a generous sanity range (0.1 <= k <= 20) per stratum and a narrower band around the overall pooled k (0.3-1.2) that would catch a regression back to naive-fit behavior. Full suite: 570/570 pass. R CMD check: 0/0/0 after adding stats::optim to @importFrom. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Per PI + Steve Goodreau review on PR #71: the Weibull-fit mean full-partnership duration can't be faithfully reproduced by a geometric-dissolution TERGM. Under Weibull with k != 1, the mean full duration and the mean cross-sectional age of extant ties diverge; TERGM with `offset(edges)` can only match the latter (via the inspection-paradox identity at exponential). Passing the length-bias-corrected Weibull mean to `dissolution_coefs()` produces a sim whose individual-partnership duration matches the Weibull but whose cross-sectional age structure is mis-matched to reality by a factor of (1 + CV^2)/2 — i.e., ~2x under k = 0.5. Adding non-geometric dissolution to tergm is a real upstream project (C-level ergm term that reads dynamic-network toggle history to get partnership age), not in ARTnet's scope. So weibull_strat as a production target was structurally wrong given the framework we're stuck with. Keeping it as a diagnostic-only option would have been defensible, but a diagnostic nobody acts on is cruft in production code. Dropped entirely. The remaining duration.method menu is: - "empirical" (default) — unadjusted median age of extant ties, geometric transformation to mean.dur.adj. - "joint_lm" — covariate-adjusted median age of extant ties (same target, regression-based), same geometric transformation. Steve's reframe makes the menu read cleanly: both methods target the same observable cross-sectional quantity that geometric TERGM can honor. joint_lm is the covariate-adjusted version of empirical — the direct duration analog of the marginal-vs-joint fix the rest of #63 applies to formation stats. Removed: - fit_weibull_dur() helper and its optim-based length-biased MLE. - weibull_strat branches in compute_alt_durs() and in the main/casl block override logic (use_direct_mean paths, weibull_shape attrs, survival package guard). - Four weibull-specific test blocks in test-duration-method.R. - stats::optim from @importFrom (no longer used). Opened issue #73 to track the remaining joint_lm gap: the stratum medians are currently computed by predicting the fitted log-duration lm at ARTnet observations, not at the synthetic target population constructed in build_netstats. This is "confounding-corrected within ARTnet" but not fully g-computed. Same pattern as #68 closing the loop on #61. Issue references point here and at the existing roxygen for forward discoverability. Full suite: 508/508 pass. R CMD check: 0/0/0. Backward-compat snapshot: 3/3 on both default and explicit method = "existing". Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
smjenness
changed the title
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds
duration.method = c("empirical", "joint_lm")tobuild_netparams()for partnership-duration estimation feeding the TERGM dissolution offset. Default"empirical"is byte-identical to pre-refactor behavior.Closes phase 3 of issue #63. Partially addresses #63; duration work is complete but one synthetic-population gap is tracked as a follow-up (issue #73).
Method menu
"empirical"(default)median(duration.time[ongoing2 == 1])per(age-match x index.age.grp)stratum"joint_lm"lm(log(duration.time) ~ index.age.grp + sqrt(index.age.grp) + hiv2 + same.age.grp + race + same.race [+ geogYN])on ongoing partnerships, per-stratum medians of fitted valuesBoth methods flow through the same geometric transformation (
mean.dur.adj = 1 / (1 - 2^(-1/median))) thatbuild_netstatspasses todissolution_coefs(), so the TERGM offset structure is identical across methods.joint_lmis the direct analog of the marginal-vs-joint fix thatmethod = "joint"applies to formation stats in #61 / #62 / #69 — same ARTnet data, same confounding-correction logic, just applied to the duration model instead of the degree / nodematch / absdiff models. The fittedlmis stored atnetparams$<layer>$joint_duration_modelfor follow-up consumption (see #73).Why no Weibull option
The PR originally introduced a
weibull_stratmethod that fit a length-bias-corrected Weibull MLE to the observed elapsed durations and used the fitted analytical mean asmean.dur.adj. Review from @sgoodreau on #72 + internal discussion surfaced that this was structurally incompatible with TERGM's geometric dissolution:offset(edges)) dissolution, mean simulated full duration equals mean simulated age of extant ties at cross-section (inspection paradox at exponential, CV = 1).k != 1, real mean full duration and real mean age of extant ties differ by a factor of(1 + CV^2) / 2— ~2x fork = 0.5.dissolution_coefs(), the simulated network gets individual-partnership durations that match reality but a cross-sectional age distribution mis-matched by 2x.empiricalandjoint_lmtarget), cross-sectional structure matches reality but individual partnerships are too short.The Weibull shape parameter
kis real and informative (ARTnet data rejects the constant-hazard assumption fairly cleanly) but fixing it would require a non-geometric dissolution term in tergm — a real upstream project, not in ARTnet's scope. Aweibull_stratkept as diagnostic-only is cruft in production code, so it was removed. The discussion is preserved in the commit history of this branch.Commit
40288d2removes:fit_weibull_dur()helper and itsoptim-based length-biased MLEweibull_stratbranches incompute_alt_durs()and in the main / casl block override logic (use_direct_meanpaths,weibull_shapeattributes,survivalpackage guard)stats::optimfrom@importFrom(no longer used)Empirical comparison (Atlanta + race = TRUE)
mean.dur.adjby stratum, passed todissolution_coefs():joint_lmshrinks the stratum medians toward the covariate-adjusted expectation — the biggest effect is on older main-matched strata where empirical, with small per-stratum samples, produced larger-than-warranted estimates driven by a few long ongoing partnerships.Validation
Backward-compat snapshot harness (default and explicit
method = "existing"):Full testthat suite: 508 / 508 pass, clean tabular output, 9.5s total.
R CMD check: 0 errors / 0 warnings / 0 notes.
GHA CI: green on
ubuntu-latestR release withARTnetDatainstalled viaEPIMODEL_PAT.Follow-up gap: #73
joint_lmcurrently aggregates stratum medians by predicting the fittedlmat ARTnet observations, not at the synthetic target population constructed downstream inbuild_netstats. This is confounding-corrected within ARTnet but not fully g-computed — the same half-step gap we had between #61 (joint models fit) and #68 (consumed inbuild_netstatsfor formation stats).Closing that gap is tracked as #73. The implementation path: move the aggregation out of
build_netparamsand intobuild_netstatsundermethod = "joint", where the synthetic population exists. The dyad-level infrastructure from #69 (pairing synthetic nodes for nodematch / absdiff) is reusable.The gap doesn't affect this PR's usability under default
method = "existing"— byte-identical output is preserved. Undermethod = "joint"+duration.method = "joint_lm", the current behavior is a strict improvement over empirical (confounding correction within ARTnet's own joint distribution), just not the full g-computation one would want for truly different target populations.Test plan
method = "existing"joint_lmstores a validlmobject atnetparams$<layer>$joint_duration_modeldissolution_coefs()needsduration.methoddoes not mutate non-duration netparams fieldsjoint_lmcomposes cleanly withmethod = "joint"inbuild_netstatsmethod = "joint"+duration.method = "joint_lm"— moves to the PR that lands Duration g-computation: predict joint_lm per-dyad on synthetic population in build_netstats #73 (first PR where the netstats object is coherently joint-driven including the synthetic-population per-dyad duration prediction).Depends on #69 (merged). Part of #63. Follow-up: #73.