Use joint GLM g-computation for build_netstats target stats (#62)

[smjenness]

Summary

Adds method = c("existing", "joint") to build_netstats(). Under method = "joint", per-synthetic-node expected degree is predicted from the joint Poisson GLMs produced by build_netparams(..., method = "joint") in #61, and target statistics are aggregated via g-computation. This makes edges and every nodefactor target internally consistent by construction (sum(nodefactor_<attr>) == 2 * edges to machine precision).

Default method = "existing" is byte-identical to pre-refactor. Closes #62.

Scope note: This PR mixes surface area from #61 and #62. During review, PI pointed out that the concurrent target should also be joint-modeled — using a binomial GLM on the deg > 1 indicator with the same joint RHS as the Poisson. That fit logically belongs in build_netparams (#61 territory) but its use is in build_netstats (#62). Fixing it split across two branches would have been awkward, so the joint concurrent binomial is included here as commit cf558a6.

What method = "joint" does NOT yet cover — these target stats still use univariate marginals and will be handled in #63 (extended to include dissolution/duration per PI request):

• nodematch_age.grp, nodematch_race, nodematch_race_diffF
• absdiff_age, absdiff_sqrt.age
• nodematch(role.class)
• diss.homog, diss.byage (durations / dissolution coefficients)

Until #63 lands, running method = "joint" end-to-end through EpiModelHIV-Template would mix joint and univariate sources in the ERGM target-stat vector — coherent enough for unit validation but not meaningful for a production estimation run.

What's refactored under method = "joint"

Field	Joint formula	Model
edges	sum(pred) / 2	Poisson joint
nodefactor_race	sum(pred[race == r]) per level	Poisson joint
nodefactor_age.grp	sum(pred[age.grp == k]) per level	Poisson joint
nodefactor_deg.casl (main)	sum(pred[deg.casl == d]) per level	Poisson joint
nodefactor_deg.main (casl)	same pattern	Poisson joint
nodefactor_deg.tot (inst)	same pattern	Poisson joint
nodefactor_diag.status	sum(pred[diag.status == h]) per level	Poisson joint
concurrent (main/casl)	sum(pred_conc) over synthetic nodes	Binomial joint on deg > 1
nodematch_*, absdiff_*	rescaled with new edges/nodefactor	(univariate ratio retained — see #63)

On concurrent specifically

First pass (commit 5ef9737) retained concurrent from the univariate binomial because the Poisson-derived P(deg > 1) = 1 - exp(-λ)(1+λ) was inflated ~5× by the truncation of deg.main at 2 and deg.casl at 3 in the training data. The Poisson rejection was correct but the conclusion was wrong: the right fix is to joint-model the binomial on the concurrency indicator, same RHS as the Poisson, different family. Commit cf558a6 does that:

• fit_joint_poisson() generalized to fit_joint_glm(..., family = poisson())
• Under method = "joint" in build_netparams, fit glm(deg.<layer>.conc ~ <joint terms>, family = binomial()) and store at netparams$<layer>$joint_concurrent_model (main and casl only — inst has no concurrent target)
• Under method = "joint" in build_netstats, out$<layer>$concurrent = sum(predict(joint_concurrent_model, synth, type = "response"))

Empirical result (Atlanta, race = TRUE, N = 10k):

• main: existing 97.09 → joint 95.41 (-1.7%)
• casl: existing 1456.31 → joint 1614.05 (+10.8%)

Both in reasonable ranges — no Poisson-style inflation. Main barely moves because the univariate intercept-only fit already captures the population-average rate well; casl moves more because the AIC-selected age.grp:deg.main interaction captures real heterogeneity the univariate can't.

Validation

Backward-compat snapshot harness (extended with netstats_extra + symmetric method =):

Call	Result
compare_to_snapshot() (default)	ALL MATCH 3/3
compare_to_snapshot(method = "existing") (now forwards to both)	ALL MATCH 3/3

Internal consistency (Atlanta, race = TRUE, N = 10k):

| Layer | |2*edges − Σnf_race| (existing) | under joint |
|---|---:|---:|
| main | 537 | 0.0 |
| casl | 319 | 0.0 |
| inst | 70 | 0.0 |

Shifted target population (race.prop = c(0.35, 0.25, 0.40)):

Layer	existing edges	joint edges	Δ
main	1990.3	1808.3	−9.1%
casl	2669.9	2929.2	+9.7%
inst	382.7	356.7	−6.8%

Tests

• tests/testthat/test-joint-netstats.R — 7 blocks, 29 assertions.
• tests/testthat/test-joint-model.R — extended to 74 assertions (new: joint_concurrent_model presence + convergence + marginal recovery).
• All 103 assertions pass.

Known minor (pre-existing, not a regression)

Dissolution-coefficient formulas capture the caller's scope environment. Under joint, that scope includes the pred_deg_* / pred_conc_* vectors (~10k doubles × 5), which silently bloats a serialized netstats. Consumers only use $coef.diss, so this is cosmetic — worth a follow-up if on-disk size matters.

Test plan

• Backward-compat snapshot match 3/3 (default + explicit method = "existing")
• Internal consistency to machine precision, all 3 layers, all nodefactor attrs
• Joint ≠ existing under a shifted target population
• Concurrent joint binomial: converges, recovers marginal, no inflation
• Unit tests 103/103 total
• End-to-end estimation with EpiModelHIV-Template using method = "joint" — deferred to after [Phase 1.3] Joint modeling for nodematch + absdiff + dissolution (dyad-level target stats) #63. The netstats object is not yet coherently joint-driven: nodematch_*, absdiff_*, and dissolution terms still use univariate marginals. Running the end-to-end now would mix joint and univariate sources in the ERGM target-stat vector, which wouldn't be a meaningful validation of either. The real end-to-end run goes on the PR that lands [Phase 1.3] Joint modeling for nodematch + absdiff + dissolution (dyad-level target stats) #63.

Depends on #61 (merged). Unblocks #63, #64, #65.