feat(public-health): NHEFS g-methods certification + full R EValue parity (#21)

Squash-merge of #21 (qijiang3). Adds NHEFS dataset + What If g-methods parity, full R EValue package parity (HR rare/common, OLS measure, exact RD grid, CI clamp), competing risks (Aalen-Johansen/Fine-Gray), epi study-design power, summary-data meta-analysis, and epi guides. Registry 1029->1031.

Pre-merge verification: new-feature suites green locally (108 passed); the 8 CI 'failures' verified as environment-sensitive pins / pre-existing budget drift, none a regression introduced by this PR (see follow-up commits for tolerance hardening + the required evalue CHANGELOG/MIGRATION correctness note).

Author: qijiang3

.gitignore

@@ -215,6 +215,11 @@ MORNING_REPORT_2026-04-29.md
 .cov_decomp/*
 !.cov_decomp/DECOMP_CAMPAIGN.md
 
+# NHEFS / What If parity dumps: regenerable byte-for-byte from the bundled
+# sp.datasets.nhefs() by running the Python side of each paired script, so we
+# do not commit ~3.4 MB of duplicated copies of the already-bundled CSV.
+tests/orig_parity/data/*_nhefs_*.csv
+
 # Tier-D/Tier-B validation campaign working dir (regenerable artifacts ignored; tracker + worklist kept)
 .tierd_campaign/*
 !.tierd_campaign/CAMPAIGN.md

README.md

@@ -138,7 +138,7 @@ StatsPAI's focus is **causal inference**. The grid below summarizes method-famil
 
 **Legend**: B = broad API coverage within this comparison table; Y = implemented entry points; P = partial, scattered, or single-algorithm support; N = no first-class entry point. These are API-breadth labels, not validation tiers.
 
-**StatsPAI at a glance**: 1,020 registered functions in the live agent registry · 81 submodules · 269k LOC (core) + 114k LOC (tests). All four numbers are reproducible from the canonical generator (`python scripts/registry_stats.py`); the per-module table in [`docs/stats.md`](docs/stats.md) is regenerated from the same script. For the API-breadth matrix (23 method families) and cross-ecosystem line-count comparison, see [`docs/stats.md`](docs/stats.md).
+**StatsPAI at a glance**: 1,031 registered functions in the live agent registry · 81 submodules · 271k LOC (core) + 134k LOC (tests). All four numbers are reproducible from the canonical generator (`python scripts/registry_stats.py`); the per-module table in [`docs/stats.md`](docs/stats.md) is regenerated from the same script. For the API-breadth matrix (23 method families) and cross-ecosystem line-count comparison, see [`docs/stats.md`](docs/stats.md).
 
 **Validation tiers matter**: `stability="stable"` means the public API is SemVer-stable; it does not by itself mean R/Stata/paper parity. Use `sp.list_functions(validation_status="certified")` for cross-language or published-reference evidence, and inspect `sp.describe_function(name)["limitations"]` before production use. See [`docs/guides/stability.md`](docs/guides/stability.md).
 

README_CN.md

@@ -46,7 +46,7 @@ StatsPAI 聚焦**因果推断**。下表描述方法家族层面的 API 覆盖
 
 **图例**：B = 本表范围内 API 覆盖较宽；Y = 有已实现入口；P = 部分、分散或单算法支持；N = 无一等入口。这些只是 API 广度标签，不是 validation tier。
 
-**StatsPAI 一句话概览**：live agent registry 中有 1,020 个注册函数 · 81 个子模块 · 269k 行核心代码 + 114k 行测试。这四个数字都可以由唯一的生成器 (`python scripts/registry_stats.py`) 现场复算；[`docs/stats.md`](docs/stats.md) 中的按模块拆分表也由同一个脚本回写。23 个方法家族的 API 广度矩阵以及跨生态行数对比，详见 [`docs/stats.md`](docs/stats.md)。这些覆盖数字描述 API 广度，不等同于每个函数都有 R/Stata 数值验证；生产使用请查看 `validation_status`。
+**StatsPAI 一句话概览**：live agent registry 中有 1,031 个注册函数 · 81 个子模块 · 271k 行核心代码 + 134k 行测试。这四个数字都可以由唯一的生成器 (`python scripts/registry_stats.py`) 现场复算；[`docs/stats.md`](docs/stats.md) 中的按模块拆分表也由同一个脚本回写。23 个方法家族的 API 广度矩阵以及跨生态行数对比，详见 [`docs/stats.md`](docs/stats.md)。这些覆盖数字描述 API 广度，不等同于每个函数都有 R/Stata 数值验证；生产使用请查看 `validation_status`。
 
 **📦 v1.16.0（2026-05-29）— 正确性修复与跨语言对齐扩展**
 

docs/guides/competing_risks.md

@@ -0,0 +1,135 @@
+# Competing risks — cumulative incidence & Fine-Gray
+
+> **When a subject can fail from more than one cause, Kaplan-Meier lies.**
+> Treating "death from other causes" as ordinary censoring makes the
+> naïve `1 - KM` curve over-state the risk of the cause you care about,
+> because it implicitly assumes the censored subjects could still have
+> had the event. StatsPAI ships the two standard competing-risks tools:
+> the **Aalen-Johansen cumulative incidence function** (`sp.cuminc`) for
+> description and Gray's test, and the **Fine-Gray subdistribution
+> hazards model** (`sp.finegray`) for regression
+> [@aalen1978nonparametric; @gray1988class; @fine1999proportional].
+
+**Event coding.** Throughout, the `event` column is an integer:
+`0` = right-censored, and `1, 2, ...` = the competing causes. This
+matches R's `cmprsk` and `survival` conventions.
+
+Scope note: these estimators are validated internally (CIF
+self-consistency, analytic-vs-bootstrap variance agreement, Fine-Gray
+recovery on simulated data, Gray's-test calibration) but are **not yet
+parity-certified** against R's `cmprsk` / `survival`. Validate your
+headline number before publication and read
+`sp.describe_function("cuminc")["limitations"]`.
+
+---
+
+## 1. Why not Kaplan-Meier?
+
+In a competing-risks setting the quantity you almost always want is the
+**cumulative incidence function** (CIF) — the probability of failing from
+cause *k* by time *t*, accounting for the fact that a competing event
+removes a subject from ever experiencing cause *k*. The Aalen-Johansen
+estimator weights each cause-specific increment by the overall (all-cause)
+survival:
+
+$$
+\hat F_k(t) = \sum_{t_i \le t} \hat S(t_{i-1}) \, \frac{d_{ki}}{n_i}
+$$
+
+so that the CIFs of all causes plus the overall survival sum to one at
+every time — a property `1 - KM_k` does not have.
+
+---
+
+## 2. Cumulative incidence — `sp.cuminc`
+
+```python
+import statspai as sp
+
+ci = sp.cuminc(df, duration="time", event="status")
+print(ci.summary())
+print(ci.cif_table.head())   # group / cause / time / cif / se / ci_lower / ci_upper
+ci.plot(cause=1)             # step CIF curve with the other causes overlaid
+```
+
+Each cause's CIF comes with a delta-method standard error and a
+confidence band (Marubini-Valsecchi / Klein-Moeschberger variance). To
+read off the cumulative incidence at a specific horizon:
+
+```python
+ci.cif_at(time=5.0, cause=1)   # CIF and CI at t = 5
+```
+
+### Comparing groups — Gray's test
+
+Pass `group=` to estimate per-group CIFs and get Gray's (1988) K-sample
+test for equality of cumulative incidence, reported per cause:
+
+```python
+ci = sp.cuminc(df, duration="time", event="status", group="arm")
+print(ci.gray_test[1])    # {'statistic', 'df', 'p_value'} for cause 1
+ci.plot(cause=1)          # one CIF curve per arm
+```
+
+Gray's test targets the **subdistribution** hazard, so it answers the
+clinically relevant question "do the groups differ in cumulative
+incidence of cause 1?" rather than the cause-specific-hazard question.
+
+---
+
+## 3. Regression — `sp.finegray`
+
+The Fine-Gray model puts covariates on the subdistribution hazard, so its
+coefficients exponentiate to **subdistribution hazard ratios (sHR)** that
+map monotonically to the cumulative incidence — a covariate with sHR > 1
+*increases* the cause-of-interest CIF. (Cause-specific Cox coefficients do
+not have this property, which is why they are easy to misread in a
+competing-risks setting.)
+
+```python
+import statspai as sp
+
+fg = sp.finegray(df, duration="time", event="status",
+                 x=["treatment", "age", "stage"], cause=1)
+print(fg.summary())
+fg.tidy()        # term / coef / shr / std_err / z / p_value / shr_lower / shr_upper
+```
+
+Subjects who fail from a competing cause are kept in the risk set with
+time-decaying inverse-probability-of-censoring weights
+`Ĝ(t)/Ĝ(T_i)` (Fine & Gray 1999); the weighted partial likelihood is
+maximised by Newton-Raphson.
+
+---
+
+## 4. Cause-specific vs. subdistribution — which to report?
+
+Both are legitimate and answer **different** questions:
+
+| Question | Tool |
+| --- | --- |
+| "What is the probability of failing from cause 1 by time *t*?" | `sp.cuminc` (CIF) |
+| "Does treatment change the **cumulative incidence** of cause 1?" | `sp.finegray` (sHR) or Gray's test |
+| "Does treatment change the **rate** of cause 1 among those still event-free?" | cause-specific `sp.cox` (censor competing events) |
+
+A common, defensible reporting choice is to present the CIF curves
+descriptively and **both** a cause-specific Cox model and a Fine-Gray
+model, since they decompose the effect on the rate vs. the effect on the
+cumulative incidence.
+
+---
+
+## 5. Limitations (read before you publish)
+
+- **Standard errors for `sp.finegray` are model-based** (inverse
+  information). A fully robust sandwich variance that accounts for
+  estimating the censoring distribution Ĝ is not yet implemented; for
+  small samples or heavy censoring, validate the SEs against R's
+  `cmprsk::crr`.
+- **No time-varying covariates** in `sp.finegray` yet.
+- **Parity is not yet certified** against `cmprsk` / `survival`.
+
+## Where to next
+
+- [Survival analysis for public health](survival_ph.md)
+- [G-methods for time-varying confounding](g_methods_ph.md)

docs/guides/epi_measures.md

@@ -0,0 +1,171 @@
+# Epidemiological measures — the `sp.epi` family
+
+> **The primitives epidemiologists reach for first.** Before any model,
+> public-health analysis starts with measures of association from
+> contingency tables, rate standardization, and screening-test accuracy.
+> StatsPAI ships these as top-level functions modelled after R's
+> `epiR`, `epitools`, and `fmsb`, with the standard small-sample
+> conventions (Woolf and exact CIs, Haldane–Anscombe continuity
+> correction, exact Poisson rate intervals) built in
+> [@rothman2008modern].
+
+All functions return a small typed result with `.estimate`, `.ci`, and
+(where defined) `.p_value`, so they compose with the rest of StatsPAI's
+reporting layer.
+
+---
+
+## 1. Measures of association (2×2 tables)
+
+The 2×2 layout throughout is `(a, b, c, d)` =
+`(exposed cases, exposed non-cases, unexposed cases, unexposed non-cases)`.
+
+```python
+import statspai as sp
+
+# Cohort study
+rr = sp.relative_risk(50, 950, 10, 990)
+print(f"RR = {rr.estimate:.2f}  95% CI {rr.ci}")
+print(rr.risk_exposed, rr.risk_unexposed)
+
+rd = sp.risk_difference(50, 950, 10, 990)
+print(f"RD = {rd.estimate:.4f}  95% CI {rd.ci}")
+
+nnt = sp.number_needed_to_treat(50, 950, 10, 990)
+print(f"NNT = {nnt.estimate:.1f}")
+
+# Case-control study (odds ratio)
+orr = sp.odds_ratio(50, 20, 30, 40, method="woolf")
+print(f"OR = {orr.estimate:.2f}  95% CI {orr.ci}")
+
+# Population attributable fraction (Levin)
+ar = sp.attributable_risk(50, 950, 10, 990)
+print(ar.estimate)
+```
+
+**Person-time / rates.** When the denominator is person-time rather than
+people, use the incidence-rate ratio with an exact Poisson interval:
+
+```python
+irr = sp.incidence_rate_ratio(
+    events_exposed=30, pt_exposed=1000,
+    events_unexposed=15, pt_unexposed=1200,
+)
+print(f"IRR = {irr.estimate:.2f}  95% CI {irr.ci}")
+```
+
+**Cross-sectional prevalence.** For prevalent (not incident) outcomes the
+prevalence ratio is preferred over the OR; `sp.prevalence_ratio` uses the
+modified-Poisson / log-binomial approach [@zou2004modified].
+
+**Zero cells** are handled automatically with the Haldane–Anscombe 0.5
+continuity correction on the log-scale standard error, so a single empty
+cell does not blow up the interval.
+
+---
+
+## 2. Stratified analysis — Mantel–Haenszel
+
+When a categorical variable confounds the exposure–outcome relationship,
+pool stratum-specific tables with Mantel–Haenszel and test whether the
+effect is homogeneous across strata [@mantel1959statistical]:
+
+```python
+import numpy as np
+import statspai as sp
+
+# K = 2 strata; each is a 2x2 [[a, b], [c, d]]
+tables = np.array([
+    [[10, 5], [8, 12]],
+    [[20, 15], [6, 9]],
+])
+
+mh = sp.mantel_haenszel(tables, measure="OR")   # or measure="RR"
+print(f"MH pooled OR = {mh.estimate:.3f}  95% CI {mh.ci}")
+print(f"stratum-specific: {mh.stratum_estimates}")
+print(f"Breslow–Day homogeneity p = {mh.homogeneity_p:.3f}")
+```
+
+A small homogeneity p-value warns that a single pooled estimate hides
+effect-measure modification — report stratum-specific estimates instead.
+`sp.breslow_day_test` exposes the test directly (with the Tarone
+correction by default).
+
+---
+
+## 3. Rate standardization
+
+Crude rates are not comparable across populations with different age
+structures. Standardize directly (when you have stratum-specific rates in
+the study population) or indirectly (SMR, when study strata are small):
+
+```python
+import statspai as sp
+
+# Direct standardization to a reference population structure
+std = sp.direct_standardize(
+    events=[50, 60, 40],            # events per age band
+    population=[1000, 2000, 1500],  # population/person-time per band
+    standard_weights=[0.3, 0.45, 0.25],
+)
+print(f"Age-standardized rate = {std.rate:.5f}  95% CI {std.ci}")
+
+# Indirect standardization → standardized mortality/morbidity ratio (SMR)
+smr = sp.indirect_standardize(
+    observed=120,
+    events_reference=[50, 60, 40],
+    population_reference=[1000, 2000, 1500],
+    population_study=[800, 1500, 1000],
+)
+print(f"SMR = {smr.estimate:.3f}  95% CI {smr.ci}")
+```
+
+---
+
+## 4. Screening and diagnostic-test accuracy
+
+```python
+import statspai as sp
+
+# From a confusion matrix
+acc = sp.sensitivity_specificity(tp=80, fn=20, fp=10, tn=90)
+print(acc.summary())   # sensitivity, specificity, PPV, NPV, LR+/LR- with Wilson CIs
+
+# From scores + labels
+roc = sp.roc_curve(y_true, scores)
+print(roc.auc)
+
+# Inter-rater agreement
+kappa = sp.cohen_kappa(rater_a, rater_b, weights="linear")
+print(kappa.estimate)
+```
+
+---
+
+## 5. Weighing the evidence — Bradford–Hill
+
+`sp.bradford_hill` turns the nine Bradford–Hill viewpoints into a
+structured, transparent scoring object — useful for causal-narrative
+sections and for forcing each viewpoint to be argued explicitly rather
+than asserted:
+
+```python
+import statspai as sp
+
+bh = sp.bradford_hill(
+    strength=1.0, temporality=1.0, consistency=0.5,
+    biological_gradient=0.8, plausibility=0.7,
+)
+print(bh.summary())
+```
+
+It is a structuring aid, **not** a causal test — the score does not
+establish causation, it documents how each viewpoint was judged.
+
+---
+
+## Where to next
+
+- [Public health & epidemiology overview](public_health.md)
+- [G-methods for time-varying confounding](g_methods_ph.md)
+- [Survival analysis](survival_ph.md)