# Inter-rater Agreement & Reliability Guide

Quantifying how well two or more raters (or a rater and a reference, or repeated
measurements) **agree**. The coefficient is easy to compute; the two ways these
analyses fail review are (1) treating **clustered** measurements as independent
(pseudoreplication) and (2) confusing **agreement** with **reliability**.

---

## When to Use

- **Cohen's kappa** — 2 raters, categorical (nominal) labels.
- **Weighted kappa** — 2 raters, **ordinal** labels (linear or quadratic weights; disagreement
  by one category counts less than by three).
- **Fleiss' kappa** — ≥3 raters, categorical.
- **Krippendorff's alpha** — any number of raters, any measurement level, tolerates missing data.
- **ICC (intraclass correlation)** — **continuous** measurements; report the model + type (below).
- **Bland–Altman** — two continuous methods/raters: bias (mean difference) + 95% limits of agreement.
- NOT for: a single 2×2 vs a reference standard (that is diagnostic accuracy — see
  `table-types/diagnostic_accuracy.md`); not for a multi-reader AI-vs-human comparison with reader +
  case variance (that is an MRMC reader study — see `table-types/reader_study.md`).

---

## Pseudoreplication comes first (this, not the coefficient, is the issue)

If each **subject contributes more than one measurement** — several lesions, aneurysms, nodules,
slices, or time-points per patient — the rows are **not independent**. Computing agreement on the
**pooled** rows (or on all pairwise distances) uses an inflated *n*, narrows the CI, and gives an
**anti-conservative** p-value. This is the single most common reliability-study error a reviewer
catches (it is flagged by the self-review probe **O18** in `observational_confounding.md`).

Two correct paths — pick one and state it:

1. **Aggregate to the independent unit first**, then compute agreement per subject. This is the
   simplest defensible analysis when a per-subject summary is meaningful (e.g. mean measurement,
   majority label, or one index lesion per subject).
2. **Model the clustering** — a mixed-effects / variance-components ICC with a **subject random
   effect** (or a GEE with an exchangeable working correlation), so the within-subject correlation
   is estimated rather than ignored.

A pooled-pairwise test can *flip* on correction: e.g. Mann–Whitney p = 0.02 on 448 pooled
pairwise distances became p = 0.59 at the per-aneurysm level (n = 112). **Report the unit of
analysis explicitly**, and when subjects have multiple measurements report a per-subject
sensitivity analysis.

### Produce the pseudoreplication-safe version

```python
import pandas as pd
import pingouin as pg   # ICC with model/type + CI

# long format: one row per (subject, measurement); rater columns rater1..raterK
df = pd.read_csv("ratings.csv")

# 1) DETECT clustering: more rows than independent subjects
n_rows, n_subjects = len(df), df["subject_id"].nunique()
if n_rows > n_subjects:
    print(f"CLUSTERED: {n_rows} measurements from {n_subjects} subjects "
          f"({n_rows / n_subjects:.1f} per subject) — do NOT pool as independent.")

# 2a) PER-SUBJECT AGGREGATION (continuous): mean per subject, then ICC on subject means
per_subj = df.groupby("subject_id")[["rater1", "rater2"]].mean().reset_index()
long = per_subj.melt(id_vars="subject_id", var_name="rater", value_name="score")
icc = pg.intraclass_corr(data=long, targets="subject_id", raters="rater", nan_policy="omit")
print(icc[["Type", "ICC", "CI95%"]])   # report Type (e.g. ICC2/ICC2k) + CI

# 2b) OR MODEL THE CLUSTERING (keep every measurement, subject random effect)
import statsmodels.formula.api as smf
df_long = df.melt(id_vars="subject_id", value_vars=["rater1", "rater2"],
                  var_name="rater", value_name="score")
m = smf.mixedlm("score ~ 1", data=df_long, groups=df_long["subject_id"])
res = m.fit()
var_between = float(res.cov_re.iloc[0, 0]); var_resid = float(res.scale)
icc_clustered = var_between / (var_between + var_resid)
print(f"variance-components ICC (subject random effect) = {icc_clustered:.3f}")
```

---

## ICC: state the model and the type (they are not interchangeable)

- **Model**: one-way random (raters differ per subject), two-way random (same raters, generalise to
  a rater population), two-way mixed (same raters, these raters only).
- **Type**: **agreement** vs **consistency** (agreement penalises systematic rater bias; consistency
  does not), and **single** vs **average** measurement (average-of-k is higher — only report it if
  the clinical use averages k raters).
- Report as e.g. **ICC(2,1) = 0.82 (95% CI 0.74–0.88), two-way random, absolute agreement, single
  rater**. An ICC with no model/type is not interpretable.

---

## Agreement is not reliability

- **Agreement** = do raters give the *same value* (absolute; Bland–Altman bias, absolute-agreement ICC).
- **Reliability** = can raters *rank/discriminate subjects* consistently (relative; consistency ICC,
  Pearson/Spearman). A method can be highly reliable yet have poor agreement (a constant offset).
  State which one the clinical claim needs, and use the matching coefficient.

---

## Reporting

- The coefficient **with a 95% CI** (bootstrap or analytic), the model/type (for ICC), and the
  **unit of analysis** (per-subject vs per-lesion, and the clustering handling).
- The interpretation band used (e.g. Landis–Koch), but do not over-interpret a point estimate whose
  CI spans two bands.
- For continuous methods: Bland–Altman **bias + 95% limits of agreement**, not just a correlation.

---

## Common failures (flag at review)

- **Pooled/pairwise agreement on clustered data** (pseudoreplication) — the headline coefficient's
  CI is too narrow; re-run per-subject or with a subject random effect (probe O18).
- **ICC reported with no model/type** — uninterpretable; the same data yields different ICCs.
- **Reliability coefficient used to claim agreement** (or vice versa) — a high consistency ICC does
  not establish that the two methods are interchangeable.
- **Correlation (r) reported as agreement** for two methods — r ignores a constant/proportional bias;
  Bland–Altman is required.
- **Kappa on ordinal labels unweighted** — treats a one-category disagreement as a full disagreement.

---

## Anti-Hallucination

- Never hand-type a coefficient or CI — compute it from the ratings CSV with a seeded script.
- Do not quote an ICC without the model/type actually estimated by the code.
- If subjects have multiple measurements, the per-subject sensitivity analysis is **mandatory** —
  do not report only the pooled number.

# Burden of Disease, Decomposition & Forecasting

Methodology reference for the **value-add analytic layers** that turn a plain descriptive
result into a top-journal burden/estimate paper: population-attributable burden, rate
decomposition, trend-break analysis, and forecasting. These are the reusable moves a
high-output epidemiology group bolts onto a standing data platform (GBD/IHME, WHO Mortality
Database, or a national registry): the disease is swapped, the shell stays fixed, and **one**
value-add layer supplies the novelty. Report the estimate against **GATHER** (`/check-reporting`)
and, for any individual-level cohort component, STROBE/RECORD.

Load this guide before generating code for: a burden-of-disease estimate, an attributable-risk
analysis, a temporal-trend / joinpoint analysis, a decomposition, or a forecast.

---

## 0. The value-add-layer playbook (pick ONE per paper)

A descriptive number (an age-standardized rate, a prevalence) is rarely publishable alone. The
lever is a single added analytic layer answering a *why / how-much / where-to* question:

| Layer | Question it answers | Method | When the clinical need supports it |
|---|---|---|---|
| **Decomposition** | *Why* did the rate change? | Das Gupta (rate → aging / population growth / epidemiological change) | a rising/falling count where "is it just ageing?" is the reviewer's first question |
| **Attribution** | *How much* is modifiable? | Comparative risk assessment / population-attributable fraction (PAF) vs a TMREL | an exposure with an established dose–response and a policy angle |
| **Trend-break** | Did a *shock* bend the trend? | Joinpoint / AAPC, pre- vs post-period | a datable policy, guideline, or event (COVID, a screening-guideline change, a reimbursement shift) |
| **Forecast** | *Where is it going*? | Age-period-cohort projection (BAPC / APC) | serial or repeated cross-sectional data with a stable trend |
| **Life-expectancy decomposition** | Which ages/causes drove ΔLE? | Arriaga (age × cause partition) | a mortality/life-table framing |

**For a single-institution health-checkup cohort** (no GBD platform): you cannot out-scale a
burden study, but three of these layers port directly onto your existing follow-up —
(1) **trend-break**: reslice your accrued follow-up around a datable guideline/scanner-era break
rather than collecting new data; (2) **forecast**: project a serial imaging trajectory (CAC,
emphysema, body-composition) forward; (3) **global framing by contextualization**: place your
individual-level effect next to the *published* GBD burden number for that disease (free from the
GHDx results tool) — this borrows the global frame without claiming to re-estimate it, so it does
**not** trigger GATHER. The one thing you cannot borrow from the ecological template is its
"robustness = uncertainty-interval propagation, no confounding control" shortcut: your data is
individual-level, so the DAG / E-value / negative-control toolkit (see `propensity_score.md`,
`survival.md`, `../style/`) still applies.

---

## 1. Uncertainty intervals (UI), not confidence intervals

Model-based estimates carry a **95% uncertainty interval** = the 2.5th–97.5th percentile of the
posterior/Monte-Carlo draws (GBD uses 250–500 draws), propagated through every modeling step.

- Report a **UI**, not a CI; a UI crossing the null means "insufficient evidence for direction,"
  not "non-significant test." Never translate a UI into *P*-value language.
- Propagate draws end-to-end: sample each input's distribution, run the full pipeline per draw,
  summarise percentiles of the output — do not add variances analytically at the end.

```r
# draw-based UI for any derived quantity f(x)
draws <- replicate(500, f(rinput()))                  # rinput() samples the input distribution
est   <- quantile(draws, c(0.5, 0.025, 0.975))        # median + 95% UI
```

---

## 2. Attributable burden — comparative risk assessment / PAF

Partition a disease burden into the fraction attributable to a modifiable exposure, against a
**theoretical minimum-risk exposure level (TMREL)**.

- **PAF** = Σ over exposure levels of `P(exposure) · (RR − 1)` / `[1 + Σ P(exposure)·(RR − 1)]`,
  evaluated relative to the TMREL. Attributable burden = PAF × total burden.
- **Relative risks** come from a meta-analysis of prospective cohorts; carry the RR curve's
  uncertainty into the draws. Reclassify the dose–response shape honestly (monotonic vs J-shaped);
  a strictly harmful exposure has TMREL = 0.
- **Mediation to avoid double-counting**: when risks act through measured mediators (a dietary
  factor → blood pressure / LDL / glucose → IHD), estimate the mediation matrix so combined PAFs
  do not exceed 100%.
- Report against **GATHER** items 8–14 (bias correction, comparability, analysis detail,
  uncertainty). R: `epitools`, `graphPAF`, or a manual draw-based implementation.

## 3. Das Gupta decomposition — why did the count change?

Partition a change in an aggregate count/rate into **population growth**, **population ageing**,
and **epidemiological (rate) change** — the three components a reviewer of any "rising burden"
claim asks to see separated.

```r
# DemoDecomp implements Das Gupta / Horiuchi decomposition of a rate function
# install.packages("DemoDecomp")
library(DemoDecomp)
# func: maps a parameter vector (age-specific rates, age structure, total pop) -> summary measure
delta <- horiuchi(func = burden_fn, pars1 = year1_vec, pars2 = year2_vec, N = 20)
# sum(delta) == observed change; group delta by its structural / rate / size components
```

Report each component with a UI; the headline is usually "X% of the increase was population
growth/ageing, not a real rise in risk."

## 4. Trend-break — joinpoint / AAPC, pre- vs post-shock

Detect where a log-linear trend changes slope and summarise it as an **average annual percentage
change (AAPC)**; compare a pre-shock and a during/post-shock window as a natural-experiment layer.

```r
library(segmented)
m  <- lm(log(rate) ~ year, data = d)                  # log-linear trend
jp <- segmented(m, seg.Z = ~year)                     # estimate joinpoint(s)
# APC per segment = 100*(exp(slope) - 1); AAPC = segment-length-weighted average of APCs
# pre/post split: fit two windows (e.g. 2010-2019 vs 2020-2023) and compare AAPCs, do not pool
```

The pre/post AAPC split is itself a robustness layer (two independent windows compared, not one
pooled trend). For a checkup cohort, the "shock" is a datable guideline/reimbursement/scanner-era
change; reslice existing follow-up around it.

## 5. Forecasting — age-period-cohort projection

Project a rate forward with a **Bayesian age-period-cohort (BAPC)** model (integrated nested
Laplace approximation), which yields predictive intervals natively.

```r
# install.packages("BAPC"); needs INLA (r-inla.org)
library(BAPC); library(INLA)
proj <- BAPC(APCList(cases, population, gf = 5),
             predict = list(npredict = 27, retro = TRUE))   # e.g. forecast to 2050 + back-test
```

- **Back-test**: withhold the last k years, forecast them, report calibration — do not ship a
  projection with no retrospective validation.
- A single predictor (SDI alone) or a single model is an ensemble caveat: state that no ensemble
  / no formal back-testing was done if it was not, rather than implying robustness.
- Alternatives: `demography`/`StMoMo` (Lee–Carter family), or a simple age-period GLM with a
  clearly stated extrapolation assumption.

## 6. Life-expectancy decomposition (Arriaga)

Partition a change in life expectancy into contributions by **age group** and, with cause-deleted
life tables, by **cause of death** — the standard framing for a mortality life-table paper.

```r
# DemoDecomp::stepwise_replacement or a direct Arriaga implementation on two life tables
library(DemoDecomp)
contrib <- stepwise_replacement(func = e0_from_mx, pars1 = mx_year1, pars2 = mx_year2)
# contrib sums to Δe0; aggregate by age band (and by cause using cause-specific mx)
```

---

## Reporting & reproducibility

- **GATHER** (`/check-reporting` `GATHER.md`) governs the estimate: define indicator/population/
  period (item 1), data sources and their biases (5–6), analysis detail with formulae (11),
  model evaluation and sensitivity (12–13), uncertainty methods and what was/was not accounted
  for (14), code access (15), and a machine-readable results file with a UI (17–18).
- Seed every RNG; emit the draw-level results, not only summaries. State the platform/version
  (GBD round, WHO Mortality Database vintage) and cite it via its data DOI (GHDx for GBD).
- Keep the causal claim honest: burden/attribution/decomposition/forecast are **descriptive or
  associational** unless an explicit causal design (a natural experiment with identification, an
  MR — see `mendelian_randomization.md`) is in place. Pair "adjustment covariate" with a
  non-causal disclaimer; do not call an adjustment covariate an instrument.

# Prediction-Model Calibration Guide

For a clinical prediction model (a logistic risk score or a Cox/survival model at a fixed
horizon), **discrimination (AUC / C-index) is not enough** — a model that ranks well can still
output probabilities that are systematically too high or too low. The ways calibration fails
review are (1) reporting **apparent (in-sample)** calibration — a slope of exactly 1.00 is the
fingerprint — with no internal-validation correction, (2) leaning on the **Hosmer–Lemeshow**
test (deprecated), and (3) omitting calibration entirely for a model meant to guide care. This
guide produces the corrected estimand; it is the produce-side of probe **S7**.

---

## When to Use

- Any model that outputs a **risk / probability** used for a decision (surveillance intensity,
  treatment eligibility, triage) — logistic or survival-at-a-horizon.
- Reported **alongside** discrimination and, for a utility claim, decision-curve net benefit —
  never discrimination alone.
- NOT a substitute for external validation: internal (bootstrap/CV) calibration corrects
  optimism but does not establish transportability.

---

## Apparent calibration is optimistic — correct it (this is the S7 issue)

Fit a logistic model by maximum likelihood and its **apparent** calibration slope on the same
data is **exactly 1.00** and its calibration-in-the-large intercept **exactly 0** — by
construction, not because the model is well-calibrated. A slope printed as `1.00` (or metrics
with no `bootstrap` / `cross-valid` / `optimism` / `held-out` token nearby) is presumptively
in-sample. Produce the **bootstrap optimism-corrected** slope instead (Harrell/Steyerberg).

```python
import numpy as np, statsmodels.api as sm
X = ...            # design matrix (add_constant), y = 0/1 outcome
def cal_slope(y, lp):
    m = sm.GLM(y, sm.add_constant(lp), family=sm.families.Binomial()).fit()
    return m.params[1], m.params[0]          # slope, calibration-in-the-large intercept

full = sm.GLM(y, X, family=sm.families.Binomial()).fit()
app_slope, app_int = cal_slope(y, X @ full.params)     # apparent: slope ~1.00, intercept ~0

rng = np.random.default_rng(42); n = len(y); opt = []
for _ in range(500):                                    # bootstrap optimism (Harrell)
    idx = rng.integers(0, n, n)
    bm = sm.GLM(y[idx], X[idx], family=sm.families.Binomial()).fit()
    s_boot, _ = cal_slope(y[idx], X[idx] @ bm.params)   # boot model on boot data (apparent)
    s_orig, _ = cal_slope(y,      X      @ bm.params)   # boot model on original data (test)
    opt.append(s_boot - s_orig)
corrected_slope = app_slope - float(np.mean(opt))       # < 1.00 when the model overfits
print(f"apparent slope {app_slope:.3f} -> optimism-corrected {corrected_slope:.3f}")
```

A corrected slope **< 1** means predictions are too extreme (overfit) and should be shrunk
(a penalized/uniform-shrinkage refit); a slope **> 1** means they are too moderate.

---

## The four levels of calibration (report weak calibration at least)

Van Calster's hierarchy — report at minimum **weak calibration** (intercept + slope):

- **Mean** (calibration-in-the-large): mean predicted = observed event rate (the intercept).
- **Weak**: intercept ≈ 0 **and** slope ≈ 1.
- **Moderate**: a **flexible calibration curve** (loess / spline of observed on predicted),
  not decile bins — the plot most reviewers now expect.
- **Strong**: correct per-covariate (rarely achievable; not required).

```python
import matplotlib.pyplot as plt
from sklearn.calibration import calibration_curve
phat = full.predict(X)
frac_pos, mean_pred = calibration_curve(y, phat, n_bins=10, strategy="quantile")
plt.plot([0, 1], [0, 1], "--"); plt.plot(mean_pred, frac_pos, "o-")   # add a loess curve for moderate
plt.xlabel("Predicted probability"); plt.ylabel("Observed frequency")
```

---

## Brier score; do NOT rely on Hosmer–Lemeshow

- **Brier score** = mean squared error of the probabilities; report the **scaled Brier**
  (1 − Brier / Brier_null) so it is interpretable against the event rate.
- **Hosmer–Lemeshow is deprecated** (Van Calster 2016; Austin & Steyerberg): its p-value depends
  on an arbitrary number of bins, it is underpowered in small samples and rejects trivially in
  large ones, and it gives no direction or magnitude. Report the **calibration slope + intercept
  + a flexible calibration plot** instead of an H–L p-value.

```python
from sklearn.metrics import brier_score_loss
brier = brier_score_loss(y, phat); scaled = 1 - brier / (y.mean() * (1 - y.mean()))
```

## Survival models

For a Cox/survival model, calibrate the **predicted vs observed risk at a fixed horizon**
(e.g. 3-year): group by predicted-risk decile and compare to a Kaplan–Meier / pseudo-value
estimate at that time, or use `rms::calibrate` / `pec` in R with bootstrap optimism correction.
State the horizon; a model can be well-calibrated at 1 year and not at 5.

---

## Reporting

- Calibration **intercept and slope** with the **internal-validation method named**
  (bootstrap/CV), not the apparent slope of 1.00; the flexible calibration plot.
- Scaled Brier; the horizon (survival); the cohort each metric was computed on (development vs
  held-out vs external) stated explicitly.
- Discrimination **and** calibration together; add decision-curve net benefit for a utility claim
  (see `table-standards/table-types/incremental_value.md` and the `make-figures` decision-curve
  exemplar).

---

## Common failures (flag at review)

- **Apparent calibration slope of exactly 1.00** (and intercept 0) with no bootstrap/CV/external
  token — in-sample fit presented as calibration (S7).
- **Hosmer–Lemeshow p-value** offered as the calibration evidence (deprecated).
- **Discrimination reported without calibration** for a model meant to guide care (S7 → MAJOR).
- **Decile-bin calibration only**, no flexible curve; or a survival calibration with no stated
  horizon.

---

## Anti-Hallucination

- Never hand-type a calibration slope/intercept, Brier, or CI — compute it from predictions with
  a seeded script.
- Do not report a calibration slope of 1.00 as evidence of good calibration — it is the
  apparent-fit artifact; report the optimism-corrected value the bootstrap produced.
- Name the validation source of every calibration number (development / bootstrap-corrected /
  external); do not present development-sample calibration as validated performance.

# Diagnostic Accuracy & Reader-Study Guide

Estimating how well an index test (or an AI model / reader) separates disease from
no-disease against a reference standard. The point estimates are easy to compute; the
ways these analyses fail review are (1) reporting AUC / sensitivity / specificity **without
CIs or on the wrong analysis unit**, (2) a **confidence-weighted / rating** score whose
novelty is never tested against the simpler **unweighted** baseline, and (3) comparing two
AUCs with a **fixed-reader** test when the claim is meant to **generalise to readers**.

---

## When to Use

- **Single index test** vs a reference standard → sensitivity, specificity, PPV, NPV, accuracy
  (each with a 95% CI), and AUC (with a DeLong or bootstrap CI).
- **Two tests / models on the same cases** → a **paired** AUC comparison (DeLong `roc.test`,
  `paired = TRUE`) — never two independent CIs eyeballed for overlap.
- **Multi-reader multi-case (MRMC)** reader study (AI-vs-reader, modality comparison) → an
  MRMC method (Obuchowski–Rockette / DBM) that carries **reader + case** variance; see
  `table-standards/table-types/reader_study.md` and `make-figures` `exemplar_plots/mrmc_roc.md`.
- NOT for: pooling accuracy across studies (that is a DTA meta-analysis — bivariate / HSROC via
  `mada`; see the `/meta-analysis` DTA path); not for rater **agreement** (that is
  `analysis_guides/agreement_reliability.md`).

---

## Every metric with a CI, on a stated analysis unit (the floor)

A bare AUC / sensitivity / specificity is not reportable. Compute the CI, and state the unit
(**per-patient vs per-lesion** — multiple lesions per patient are clustered, exactly as in
`agreement_reliability.md`; a per-lesion count inflates *n*).

```python
import numpy as np, pandas as pd
from sklearn.metrics import roc_auc_score
from statsmodels.stats.proportion import proportion_confint

df = pd.read_csv("reader_calls.csv")   # truth (0/1), call (0/1), confidence (1..K), stratum

# sensitivity / specificity / PPV / NPV at the operating point, each with a Wilson CI
tp = int(((df.truth == 1) & (df.call == 1)).sum()); fn = int(((df.truth == 1) & (df.call == 0)).sum())
tn = int(((df.truth == 0) & (df.call == 0)).sum()); fp = int(((df.truth == 0) & (df.call == 1)).sum())
for name, num, den in [("sensitivity", tp, tp+fn), ("specificity", tn, tn+fp),
                       ("PPV", tp, tp+fp), ("NPV", tn, tn+fn)]:
    lo, hi = proportion_confint(num, den, method="wilson")
    print(f"{name} = {num/den:.3f} (95% CI {lo:.3f}-{hi:.3f}), n={den}")
```

PPV and NPV are **prevalence-dependent** — do not transport them to a population with a
different base rate; report sensitivity/specificity (prevalence-invariant) plus the study
prevalence, and recompute predictive values at the target prevalence with Bayes' rule.

---

## The confidence-weighted trap comes first (this, not the AUC, is the issue)

When the novelty is a **confidence-weighted / rating-collapsed** score used as the ROC
predictor, you must (a) confirm the (call × confidence) encoding is **strictly monotone** — a
folded encoding silently collides the most-confident-positive with a negative call — and
(b) report the **unweighted binary-call AUC** beside the weighted one, so the weighting earns
its place. This is the produce-side of self-review / peer-review probe **D9**.

```python
K = int(df.confidence.max())

# intended strictly-monotone composite: negative calls rank below positive calls, and within
# a call higher confidence pushes further from the decision boundary.
def composite(call, conf, K):
    return np.where(call == 1, K + conf, K + 1 - conf)   # neg: 1..K, pos: K+1..2K

# (1) MONOTONIC-ENCODING CHECK — every distinct (call, confidence) must map to a distinct,
#     correctly-ordered score. A collision is the folded-score bug (e.g. real/1 == ai/1).
combos = df[["call", "confidence"]].drop_duplicates().copy()
combos["score"] = composite(combos.call.values, combos.confidence.values, K)
if combos.score.nunique() != len(combos):
    raise ValueError("ENCODING COLLISION — the confidence weighting is not strictly monotone "
                     "(folded-score bug); fix the encoding before computing AUC.")

# (2) UNWEIGHTED BASELINE beside the weighted primary
df["score"] = composite(df.call.values, df.confidence.values, K)
print(f"AUC (confidence-weighted) = {roc_auc_score(df.truth, df.score):.3f}")
print(f"AUC (unweighted binary call) = {roc_auc_score(df.truth, df.call):.3f}")
```

If the weighted AUC materially exceeds the unweighted one, the weighting must be justified;
if it does not, report the simpler estimator as primary. A weighted score that changed a
hypothesis's direction versus its unweighted baseline is a **post-lock change to disclose**,
not a silent primary.

---

## AUC confidence intervals and comparing two AUCs

DeLong is the field-standard analytic AUC CI and paired comparison; `pROC` (R) is canonical.

```r
library(pROC)
d  <- read.csv("reader_calls.csv")
rw <- roc(d$truth, d$score, quiet = TRUE)
ci.auc(rw)                                   # DeLong 95% CI for the weighted AUC
ru <- roc(d$truth, d$call,  quiet = TRUE)
roc.test(rw, ru, method = "delong", paired = TRUE)   # paired, SAME cases (fixed-reader)
```

A portable Python bootstrap CI (use when a DeLong implementation is unavailable):

```python
def auc_ci_bootstrap(y, s, n_boot=2000, seed=42):
    rng = np.random.default_rng(seed); y = np.asarray(y); s = np.asarray(s); n = len(y)
    boots = [roc_auc_score(y[i], s[i]) for i in (rng.integers(0, n, n) for _ in range(n_boot))
             if len(np.unique(y[i])) == 2]
    return tuple(np.percentile(boots, [2.5, 97.5]))
```

**Fixed-reader vs MRMC.** A paired DeLong test treats the readers/cases as fixed. If the claim
is that the result **generalises to readers** (a reader sample), a fixed-reader CI understates
uncertainty — use an MRMC method (Obuchowski–Rockette / DBM) that carries reader + case
variance, report per-reader and reader-averaged AUC, and state the unit (per-patient vs
per-lesion) and the superiority/non-inferiority margin. This is probe D9's MRMC lead.

---

## Per-stratum admissibility: produce the table the claim is checked against (D10)

A blanket "no subgroup met AUC ≥ 0.75 with lower bound ≥ 0.70" is falsified by a single tabled
stratum that meets it. Produce the per-stratum AUC + CI and test each row against the rule
rather than asserting a global negative.

```python
rule = lambda auc, lo: (auc >= 0.75) and (lo >= 0.70)
for name, g in df.groupby("stratum"):
    if g.truth.nunique() < 2:                       # single-class stratum: not estimable
        print(f"{name}: not estimable (one class)"); continue
    auc = roc_auc_score(g.truth, g.score); lo, hi = auc_ci_bootstrap(g.truth, g.score)
    print(f"{name}: AUC {auc:.3f} (95% CI {lo:.3f}-{hi:.3f}) "
          f"{'QUALIFIES' if rule(auc, lo) else 'below'}")
```

With *k* strata some crossings are expected — frame qualifying strata as hypothesis-generating
(note the multiplicity), but do **not** deny a row the paper's own table satisfies.

---

## One scale per comparison (D11)

Two values sharing a comparison column (or a row-wise "A vs B") must be computed under the
**same normalisation / definition** — e.g. both volume errors as a standard relative error, not
one relative and one range-normalised. If they are not identical, recompute on a common scale
before any superiority/comparability claim, or footnote "not on the same scale".

---

## Reporting

- Sensitivity, specificity, PPV, NPV, accuracy, and AUC each **with a 95% CI**; the operating
  point and how it was chosen (Youden vs a prespecified clinical threshold — a threshold picked
  on the same data is optimistic and needs a held-out or cross-validated estimate).
- The **analysis unit** (per-patient vs per-lesion) and clustering handling; study prevalence.
- For a weighted-score primary: the **unweighted-baseline AUC** beside it (D9).
- AUC alone is insufficient for a clinical claim — pair discrimination with **calibration** and a
  **decision-curve / net-benefit** pass at the relevant threshold (see the incremental-value
  table type and the `make-figures` decision-curve exemplar).

---

## Common failures (flag at review)

- **AUC / sensitivity / specificity with no CI**, or a per-lesion count reported as if per-patient.
- **Confidence-weighted AUC with no unweighted baseline** and no monotonic-encoding check (D9) —
  the weighting may have created the result or hidden a folded-score bug.
- **Two AUCs compared by CI overlap** instead of a paired DeLong test; a **fixed-reader** CI used
  for a claim that generalises to readers (needs MRMC).
- **"No stratum met the rule" contradicted by the per-stratum table** (D10).
- **Mixed-normalisation head-to-head** in one column (D11).
- **PPV/NPV transported** across a prevalence change.

---

## Anti-Hallucination

- Never hand-type an AUC, sensitivity/specificity, or CI — compute it from the calls CSV with a
  seeded script; carry each estimate together with its CI (never a bare AUC).
- Do not report a weighted-score AUC without the unweighted baseline the code produced.
- Do not quote a DeLong comparison p-value that a paired test on the same cases did not produce.

# Health Economic Evaluation Analysis Guide

Comparing the **costs and consequences** of two or more interventions to inform a coverage,
adoption, or treatment decision. The headline — an **incremental cost-effectiveness ratio (ICER)**
— is arithmetically trivial; what fails review is **the structural choices behind it** (perspective,
time horizon, discounting, the effectiveness source, the cost basis, the model, and the propagation
of uncertainty). So the analysis is mostly the costing, the model, and the sensitivity suite, not the
ratio. This is the analysis-side companion to review probes **HE1–HE8** in
`health_economic_evaluation.md` and the **CHEERS 2022** reporting checklist.

---

## When to Use

- A comparative decision question (adopt / reimburse / treat) where both **cost** and **health
  outcome** differ between options: cost-effectiveness (CEA, natural units — life-years, events
  avoided), **cost-utility** (CUA, QALYs — the default for reimbursement), cost-benefit (CBA,
  monetised outcomes), cost-minimisation (only when outcomes are demonstrably equivalent), or
  **budget-impact** (affordability, distinct from cost-effectiveness).
- **Trial-based** (within-RCT patient-level costs and outcomes) or **decision-model-based** (decision
  tree for short horizons; **Markov / state-transition** or **discrete-event simulation** when timing
  and recurrence matter; lifetime horizons).
- NOT for: a costing/burden-of-illness description with no comparator (not an economic *evaluation*);
  asserting "cost-effective" from a point ICER with no uncertainty analysis; cost-minimisation when a
  non-inferiority outcome claim has not actually been established.

## Core quantities

- **Incremental cost** ΔC and **incremental effect** ΔE between an intervention and its comparator;
  the **ICER = ΔC / ΔE** (e.g. cost per QALY gained). With ≥3 options, rank by cost, remove
  **dominated** (more costly, less effective) and **extended-dominated** strategies, then compute
  ICERs sequentially along the efficient frontier.
- **Net benefit** at a willingness-to-pay threshold λ: **INMB = λ·ΔE − ΔC** (monetary) or
  **INHB = ΔE − ΔC/λ** (health). Net benefit is linear and avoids the ICER's quadrant ambiguity, so
  it is preferred for regression and for probabilistic summaries.
- **QALYs** = time × **utility** (preference-based, 0=dead, 1=full health) from a named instrument
  (EQ-5D-3L/5L, SF-6D, HUI) and a stated **value set/tariff** for the relevant country.

## The structural choices (state and justify each)

- **Perspective** — healthcare-system/payer vs **societal**; determines which costs count (societal
  adds productivity and informal-care costs). Apply it consistently.
- **Time horizon** — long enough to capture all differential costs and effects; **lifetime** for
  chronic disease or interventions with lasting effects. Extrapolation beyond trial data must be
  modelled explicitly (e.g. parametric survival extrapolation) and tested.
- **Discounting** — apply the jurisdiction's reference-case rate (commonly ~3% or 3.5%) to **both
  costs and outcomes**; sensitivity at alternative rates.
- **Costing** — report **resource quantities and unit costs separately**; state **currency, price
  year**, and inflation/currency conversion. Match cost categories to the perspective.
- **Model** — justify structure against natural history; state cycle length and half-cycle
  correction (Markov); validate (internal/face/external/predictive); test **structural** uncertainty
  via scenario analysis.

## Uncertainty (the analytic core)

- **Deterministic** — one-way and **tornado** diagrams to find the drivers; multi-way / scenario
  analyses for methodological and structural choices (discount rate, time horizon, alternative
  model structures).
- **Probabilistic (PSA)** — assign each uncertain parameter a **distribution** (beta for
  probabilities and utilities; gamma or log-normal for costs; Dirichlet for transition-probability
  sets), propagate by **Monte Carlo**, and report the **cost-effectiveness plane** (the cloud of ΔC,
  ΔE draws) and the **cost-effectiveness acceptability curve (CEAC)** — P(cost-effective) across a
  range of λ. Report results as net benefit at the relevant threshold with its uncertainty, not a
  bare point ICER.
- **Value of information** (EVPI/EVPPI) is an optional extension quantifying the expected cost of
  current decision uncertainty / the priority parameters for further research.

## Reporting & tools

- Report to **CHEERS 2022** (28 items): perspective, horizon, discount rate, currency/price year,
  model rationale, study parameters with distributions, disaggregated costs/outcomes, the ICER,
  and the uncertainty analysis (plane + CEAC). State the **willingness-to-pay threshold** and make
  the "cost-effective" conclusion conditional on it and on the CEAC probability.
- Tools: R — `heemod` / `dampack` / `hesim` (state-transition + PSA + CEAC), `BCEA` (Bayesian
  cost-effectiveness, CEAC/EVPI), `survival`/`flexsurv` (survival extrapolation); also TreeAge, or
  spreadsheet models with a documented PSA. Make the model or the health-economic analysis plan
  available for scrutiny.
- Companion review probes: **HE1–HE8** (`peer-review`/`self-review`
  `references/domain-probes/health_economic_evaluation.md`).

# Mediation Analysis Guide

Decomposing an exposure–outcome association into a path **through** a mediator (indirect
effect) and a path **not through** it (direct effect). The estimate is easy; the
*identification* is hard — that is where mediation analyses fail review.

---

## When to Use

- A pre-specified hypothesis that exposure X affects outcome Y **partly via** mediator M
- M is measured, plausibly on the causal path, and (ideally) measured **after X and before Y**
- NOT for: a single-timepoint design used to assert a causal X→M→Y chain (see the design caveat
  below — that is the most common reason mediation papers are rejected); not a substitute for a
  longitudinal or experimental test of the mechanism

---

## Identification comes first (this, not the bootstrap, is the issue)

A bootstrapped indirect-effect CI quantifies the **sampling variability** of the a×b product. It
says nothing about whether the indirect effect is **identified**. Causal mediation requires, beyond
no exposure–outcome and no exposure–mediator confounding:

- **No unmeasured mediator–outcome confounding (sequential ignorability).** This must hold *even in
  a trial that randomizes X*, because M is not randomized. In observational data it essentially
  never holds unmodelled — so a **sensitivity analysis is mandatory**, not optional.
- **No mediator–outcome confounder affected by the exposure** (an exposure-induced M–Y confounder
  breaks the standard product/difference method; use a method that allows it, e.g. an interventional
  or G-computation estimator).
- **Correct temporal order.** X precedes M precedes Y. A **cross-sectional** design measures all
  three together and cannot establish this — report the indirect effect as *consistent with*
  mediation, state plainly that the design cannot order X/M/Y, and reserve the causal claim for a
  two-wave / longitudinal design (review probe O13 in `observational_confounding.md`).

## Method

- **Estimator**: bootstrapped product-of-coefficients (a×b) is standard for continuous M and Y
  (R `mediation`, `CMAverse`, or Hayes PROCESS; Python `pingouin`/`statsmodels`). For a binary
  outcome use the **counterfactual / natural-effects** decomposition (e.g. `CMAverse`,
  `regmedint`), not the naive product on the odds-ratio scale (non-collapsibility distorts it).
- **Bootstrap** the indirect effect (≥ 2000, ideally 5000 resamples) and report the **bias-corrected
  percentile CI** — not the Sobel test (Sobel assumes normality of a×b and is underpowered).
- **Exposure–mediator interaction**: with the counterfactual approach, report the natural direct and
  indirect effects allowing X×M interaction rather than assuming it away.

## Reporting (AGReMA)

Report against **AGReMA** (A Guideline for Reporting Mediation Analyses). Minimum:

- The full mediation model (every path), the estimator, and the confounder set for **each** of the
  X→Y, X→M, and M→Y relationships (they differ — name them separately).
- Total, direct, and indirect effects, each with a bootstrap CI.
- **Proportion mediated only with uncertainty, and only when the total effect is well-estimated.**
  Proportion mediated = indirect / total is **unstable when the total effect is small or near-null**
  — it can exceed 100% or flip sign. Do not headline a proportion-mediated when the total-effect CI
  is wide or crosses the null.
- The **sensitivity analysis** for unmeasured mediator–outcome confounding: an **E-value for the
  indirect effect** (or a ρ-/correlation-based mediation sensitivity, e.g. `mediation::medsens`),
  stating how strong unmeasured M–Y confounding would have to be to null the indirect effect.
- The temporal structure of X, M, Y (and, if cross-sectional, the explicit caveat).

## Common failures (flag at review)

- Causal "X affects Y through M" from cross-sectional data with no temporal caveat → reframe to
  association-level + add the sensitivity analysis (O13).
- A significant bootstrap CI treated as evidence of identification.
- No sensitivity analysis for unmeasured M–Y confounding.
- Proportion-mediated quoted as a stable figure when the total effect is small.
- Binary-outcome mediation via the naive OR product instead of a counterfactual decomposition.

## Anti-Hallucination

- Never fabricate path coefficients, indirect effects, or bootstrap CIs — all from executed code.
- Do not assert a causal direction the design cannot support; state the identification assumptions.
- Generate references via `/search-lit` (AGReMA, VanderWeele).

# Mendelian Randomization Analysis Guide

Using germline genetic variants as **instrumental variables (IVs)** for an exposure to
strengthen causal inference when an RCT is infeasible or unethical. The point estimate (a ratio
of GWAS associations) is easy; what fails review is **instrument validity, pleiotropy, and the
assumptions** — so the analysis is mostly the sensitivity suite, not the headline number. This is
the analysis-side companion to review probes **MR1–MR8** in `mendelian_randomization.md`.

---

## When to Use

- A causal question where an RCT is infeasible/unethical AND a **strong genetic instrument** for
  a **modifiable, well-defined, heritable** exposure exists.
- two-sample summary-data MR (most common), one-sample MR, multivariable MR (MVMR), drug-target /
  cis-MR, non-linear MR (NLMR).
- NOT for: a non-modifiable / ill-defined / composite exposure; a confirmatory **clinical-effect-
  size** claim — MR estimates a **lifelong genetic-proxy effect direction**, not the magnitude of
  a drug/behaviour change started in adulthood.

## The three IV assumptions (state and evidence each)

- **Relevance** — the instrument is strongly associated with the exposure → report the
  **F-statistic** and variance explained (R²).
- **Independence / exchangeability** — the instrument is independent of confounders → ancestry-
  matched GWAS + principal components; a PhenoScanner-style confounder scan.
- **Exclusion restriction** — the instrument affects the outcome only through the exposure (no
  **horizontal pleiotropy**). Untestable directly → addressed via the sensitivity suite below.

## Instrument construction

- Select at genome-wide significance (P < 5 × 10⁻⁸), **LD-clump** (e.g. r² < 0.001, 10 Mb window;
  report a stricter-r² sensitivity), handle **palindromic/ambiguous** SNPs, and **harmonize effect
  alleles** across the exposure and outcome GWAS (a harmonization error flips the sign).
- **Instrument strength**: per-SNP F = (β_GX / se_GX)²; overall F = [R²(N − k − 1)] / [(1 − R²)k];
  a mean F well above ~10. For **MVMR** report the **conditional F**.
- **Winner's curse**: select instruments in a GWAS **independent** of the one used for the β_GX
  estimate; selecting and estimating in the same sample inflates the instrument–exposure effect.
- **cis / drug-target**: variants within ±(gene window) of the target; correlated cis variants need
  a **GLS-corrected IVW** with the LD matrix (naive IVW assumes independent instruments).

## Estimation + sensitivity suite (pre-specify ALL, not IVW alone)

- **IVW** (random-effects) as the primary estimator.
- **MR-Egger** (its **intercept** tests directional pleiotropy), **weighted median** (valid if ≤50 %
  of weight is invalid), **weighted mode**, **MR-PRESSO** (outlier detection/correction).
- **Heterogeneity**: Cochran's Q, I²_GX; **leave-one-out** and **single-SNP** diagnostics.
- **Concordance across methods is the robustness claim** — not the IVW point estimate alone.

```r
library(TwoSampleMR)   # or the MendelianRandomization package
dat <- harmonise_data(exposure_dat, outcome_dat)          # align effect alleles
res <- mr(dat, method_list = c("mr_ivw", "mr_egger_regression",
                               "mr_weighted_median", "mr_weighted_mode"))
mr_pleiotropy_test(dat)    # MR-Egger intercept (directional pleiotropy)
mr_heterogeneity(dat)      # Cochran's Q
# MRPRESSO::mr_presso(...) for outlier detection/correction
# F-stat: with R2 and N -> F = (R2*(N-k-1))/((1-R2)*k)
```

## Direction, overlap, ancestry

- **Reverse causation**: Steiger filtering (the instrument should explain more variance in the
  exposure than the outcome) and/or **bidirectional MR**.
- **Sample overlap**: two-sample MR assumes independent samples. Report the overlap fraction; overlap
  biases toward the confounded observational estimate and, with weak instruments, inflates type-1
  error → use non-overlapping samples or an overlap-aware correction, or state the bias direction.
- **Ancestry**: instruments and LD reference must match the GWAS ancestry; control population
  stratification.

## MVMR, drug-target, non-linear

- **MVMR** for measured pleiotropy / mediation — estimate conditional (direct) effects; report the
  **conditional F** for instrument strength.
- **Drug-target / cis-MR**: report **colocalization** (`coloc`) to separate a shared causal variant
  from LD confounding, a **positive control** (a known on-target effect), and a **phenome-wide
  adverse-effect scan** for any safety claim.
- **Non-linear MR**: the residual and doubly-ranked methods can produce **artefactual shapes** — a
  quoted inflection/threshold is not validated. Require **negative + positive control outcomes**,
  **sensitivity excluding extreme strata**, biological plausibility, and triangulation before
  trusting an NLMR shape (the genetic analogue of data-driven threshold mining).

## Reporting (STROBE-MR)

Report against **STROBE-MR**: data sources + sample sizes, instrument-selection criteria and
strength, the three assumptions and how each was addressed, harmonization, the full sensitivity
suite, sample overlap, and ancestry. Interpret the estimate as a **lifelong genetic-proxy effect
direction**, not a clinical-intervention magnitude. Review-side probes: MR1–MR8 in
`mendelian_randomization.md`.

# Missing Data Guide

Strategies for handling missing data in medical research analyses.
This is a preprocessing step, not an independent analysis type.

---

## Missing Data Mechanisms

| Mechanism | Definition | Example | Implication |
|-----------|-----------|---------|-------------|
| **MCAR** | Missing completely at random | Specimen lost in transit | Complete case analysis OK |
| **MAR** | Missingness depends on observed data | Older patients more likely to drop out | **Multiple imputation needed** |
| **MNAR** | Missingness depends on unobserved value | Sicker patients miss follow-up | Sensitivity analysis required |

MNAR cannot be verified statistically — always consider as a possibility.

---

## Decision by Missing Rate

| Missing % | Action |
|-----------|--------|
| < 5% | Complete case analysis generally acceptable |
| 5-20% | **Multiple imputation recommended** |
| 20-40% | MI + sensitivity analysis mandatory |
| > 40% | Consider excluding the variable from analysis |

---

## Multiple Imputation via Chained Equations (MICE)

### Procedure
1. Specify imputation model for each variable with missing data
2. Generate m imputed datasets (m = 5 minimum; m = 20 if missing > 20%)
3. Analyze each dataset independently (same analysis model)
4. Pool results using **Rubin's rules**: estimate = mean of m estimates; variance = within + between variance

### Imputation Model Rules
- Include ALL variables from the analysis model in the imputation model
- Include the outcome variable (in imputation model only)
- Match variable type to method:
  - Continuous → predictive mean matching (PMM) or regression
  - Binary → logistic regression
  - Multinomial → polytomous regression
  - Ordinal → proportional odds

### Methods NOT Recommended
- **Mean imputation**: underestimates variance — never use
- **LOCF (Last Observation Carried Forward)**: biased in most settings — avoid
- **Single imputation**: does not account for imputation uncertainty

---

## Reporting Template

"Missing data ranged from X% (variable A) to Y% (variable B). Little's MCAR test suggested data were not MCAR (P = 0.003). Multiple imputation using chained equations (MICE) was performed with m = 20 imputed datasets. Results were pooled using Rubin's rules. Sensitivity analyses using complete case analysis yielded consistent results (Supplementary Table X)."

---

## When to Trigger in analyze-stats

During Phase 1 (Data Assessment), if any analysis variable has > 5% missing:
1. Report missing counts and percentages
2. Suggest MICE before proceeding
3. Generate imputation code as a preprocessing step
4. Run the primary analysis on imputed data
5. Run complete case analysis as sensitivity analysis

---

## Python Implementation

```python
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
import numpy as np
import pandas as pd

m = 20  # number of imputations
results = []
for i in range(m):
    imputer = IterativeImputer(max_iter=10, random_state=i, sample_posterior=True)
    df_imputed = pd.DataFrame(imputer.fit_transform(df), columns=df.columns)
    # Run analysis on each imputed dataset
    # ... append results
# Pool using Rubin's rules
pooled_estimate = np.mean(results, axis=0)
```

## R Implementation

```r
library(mice)
imp <- mice(df, m = 20, method = 'pmm', seed = 42)
fit <- with(imp, lm(y ~ x1 + x2))
pooled <- pool(fit)
summary(pooled)
```

---

## Common Reviewer Flags

1. Not reporting missing data counts per variable
2. Using listwise deletion without justification
3. Mean imputation or LOCF without acknowledging limitations
4. Not performing sensitivity analysis (complete case vs MI comparison)
5. Not stating the assumed missing mechanism (MCAR/MAR/MNAR)

# Multiple Testing & High-Dimensional Screening Guide

Testing **many exposures/predictors against one (or a few) outcomes at once** — an
exposome-/environment-/metabolome-/proteome-/nutrient-wide association scan (ExWAS / EWAS /
MWAS), or any "we screened N candidate predictors and report the significant ones" pass. The
test is easy; what fails review is the **multiplicity bookkeeping, the replication, and the
selective reporting**. This is the analysis-side companion to review probe **O17** in
`observational_confounding.md`.

---

## When to Use

- A scan that tests many exposures/biomarkers/features against an outcome simultaneously.
- Any procedure that **selects** "significant" predictors out of a large candidate set.
- NOT for: a single pre-specified exposure (use `regression.md`); a search for a cut-point or
  non-linear shape **within one** exposure (a different multiplicity — review probe O12); a
  confirmatory model with a handful of pre-registered covariates.

---

## Match the correction to the claim (FWER vs FDR)

The two control different things — pick by what the manuscript will claim, and **report the
method together with `m`, the number of tests** (the denominator).

- **FWER — family-wise error rate** (Bonferroni, Holm, or a permutation-based study-wide
  threshold). Controls the probability of **any** false positive. Use for a **confirmatory
  single-hit** claim. Bonferroni threshold = α / m (e.g. a proteome scan of 1463 proteins at a
  GWAS-style 5 × 10⁻⁸ gives ≈ 3.4 × 10⁻¹¹). A **permutation-based** "exposome-wide / metabolome-
  wide significance level" estimates the threshold that controls FWER **accounting for the
  correlation** among exposures, so it is less conservative than raw Bonferroni.
- **FDR — false discovery rate** (Benjamini–Hochberg q-value; Benjamini–Yekutieli under
  arbitrary dependence). Controls the **expected proportion of false positives among the
  declared hits**. Use for **discovery/screening**, and then **frame the result as hypothesis-
  generating, not confirmatory**.
- **Honest denominator.** Apply the correction to the **whole tested set** — count every
  exposure tested a priori, including ones too sparse to model well. Shrinking `m` to the
  favoured hits, or quoting raw p < 0.05 across hundreds of tests, is the core error.

## Replication is the real safeguard (not the correction alone)

Agnostic scans carry a **high false-discovery proportion even after FDR**. The load-bearing
control is an **independent replication**:

- a **split-half** discovery/replication within the cohort (keep the case proportion matched
  across splits), a **second cohort**, or a **second survey cycle** (e.g. NHANES cross-cycle);
- require **directional concordance** and **report the replication rate** (e.g. "110/164
  exposures replicated", "59% survived correction and were directionally concordant");
- a held-out **validation** set is used **once**, for the final model — never reused as a second
  discovery pass.

A single-cohort FDR-significant scan with no replication is **exploratory**; its top hits are
candidates, not findings.

## Correlated exposures (two consequences)

Exposome/omics exposures are heavily intercorrelated, so:

1. **Raw Bonferroni is over-conservative** — the effective number of independent tests is
   smaller than `m`. Use a permutation threshold, or an **effective-number-of-tests** adjustment
   (eigenvalue decomposition of the exposure correlation matrix; Li & Ji / `poolr::meff()`).
2. **A univariate hit can be a marker for a correlated true cause** — the single-exposure model
   ignores co-exposure confounding / mixtures. Before any causal reading, address the structure
   with **clustering / dimension reduction** or a **multi-exposure / mixture model**
   (elastic-net, WQS regression, quantile g-computation, or BKMR).

## Code patterns

```r
# FDR / FWER adjustment over a vector of p-values from the scan
padj_fdr <- p.adjust(pvals, method = "BH")     # discovery / screening
padj_fwer <- p.adjust(pvals, method = "holm")  # or "bonferroni" — confirmatory
q <- qvalue::qvalue(pvals)$qvalues             # Storey q-values

# effective number of independent tests under correlated exposures
m_eff <- poolr::meff(R = cor(exposure_matrix), method = "liji")  # then alpha / m_eff

# permutation-based study-wide significance threshold (correlation-aware)
# shuffle the outcome B times, refit the scan, keep the min p per permutation;
# the (alpha)-quantile of that min-p null is the study-wide threshold.
```

For mixture / multi-exposure modelling use `glmnet` (elastic-net), `gWQS` (WQS), `qgcomp`
(quantile g-computation), or `bkmr` (Bayesian kernel machine regression).

## Reporting checklist (what reviewers check)

- **Number of exposures tested (`m`)** stated; correction applied to the whole set.
- **Correction method named and matched to the claim** — FWER for confirmatory, FDR for
  discovery (and then framed as hypothesis-generating).
- **Replication design + rate** reported (split-half / second cohort / cross-cycle, directional
  concordance).
- **Full results in a supplement** — every tested exposure with its effect size and p/q, not
  only the winners; pre-registration of the exposure list is ideal (guards against
  selective-reporting / HARKing). Nominal (uncorrected) hits may be shown only if **labelled as
  not surviving correction**.
- **Effect sizes + the resolution floor** — with the large N these scans run on, trivial effects
  clear any threshold; report magnitudes and note that a permutation with `k` permutations
  cannot resolve p below ≈ 1/k. A bare "number of significant exposures" overstates the finding.
- **Complex-survey data** (NHANES / KNHANES / CHNS): combine **design-based standard errors**
  (weights + strata + PSU — see `survey_weighted.md`) **with** the multiplicity correction, not
  one or the other.

# Network Meta-Analysis (NMA) Analysis Guide

Comparing three or more interventions at once by combining **direct** (head-to-head) and
**indirect** evidence across a network, usually with a treatment **ranking**. The pooling is
standard meta-analysis; what fails review is the **transitivity assumption, consistency between
direct and indirect evidence, network geometry, ranking interpretation, and network-level
certainty**. This is the analysis-side companion to review probes **NM1–NM8** in
`network_meta_analysis.md`; for the pairwise machinery (search, screening, the random-effects
model, study-count thresholds) use the SR/MA workflow.

---

## When to Use

- Synthesis of ≥3 interventions with a connected evidence network (direct + indirect), or a
  **component NMA** decomposing multicomponent interventions.
- NOT for: a two-treatment pairwise meta-analysis (standard random-effects MA); a DTA meta-analysis
  (bivariate/HSROC).

## Transitivity (assess before pooling)

- An NMA is valid only if **effect modifiers are distributed similarly across the comparisons**
  (transitivity). **Assess and report it**: compare the distribution of plausible effect modifiers
  (baseline severity, age, dose, follow-up, year) across comparisons (box plots / a table), and/or
  network meta-regression. Transitivity is a clinical/epidemiological judgment — no test replaces it.

## Consistency / incoherence (the statistical footprint)

- Test **globally** (design-by-treatment interaction model) AND **locally** (node-splitting /
  loop-specific / back-calculation). A **star network (no closed loops)** cannot be checked — say so;
  the result then rests entirely on transitivity.
- Investigate the **source** of any inconsistency (often a single trial/loop), don't just report a
  global p value.

```r
library(netmeta)              # frequentist
nm <- netmeta(TE, seTE, treat1, treat2, studlab, data = d, sm = "OR", common = FALSE)
netsplit(nm)                  # local incoherence (direct vs indirect per comparison)
decomp.design(nm)             # global design-by-treatment / Q decomposition
netheat(nm)                   # net heat plot for inconsistency hotspots
netrank(nm, small.values = "bad")   # P-scores (frequentist ranking)
funnel(nm, order = ...)       # comparison-adjusted funnel plot
# Bayesian alternative: BUGSnet / gemtc / multinma (node-split, SUCRA, model fit via DIC)
```

## Network geometry, heterogeneity, ranking

- Present a **network plot** (node size ∝ sample size / #studies; edge thickness ∝ #trials) so
  thin (single-trial) edges and dominant nodes are visible; sparse/poorly-connected networks give
  fragile estimates.
- Report **between-study heterogeneity** (global τ²) and examine the usual **common-heterogeneity**
  assumption; high heterogeneity undermines pooling and transitivity.
- **Ranking (SUCRA / P-score / rankograms) is not a test of superiority** and is unstable when
  estimates are imprecise. Report ranking **with** the effect estimates (a **league table**), their
  intervals, and certainty — never "X is best" from SUCRA alone.

## Certainty (CINeMA / GRADE-NMA)

- Rate certainty **per network estimate** (within-study bias, reporting bias, indirectness,
  imprecision, heterogeneity, incoherence) via **CINeMA** or GRADE-NMA; **downgrade indirect-only**
  comparisons. Scale conclusions to certainty — a low-certainty estimate cannot support "superior."

## Component NMA

- A **component NMA** assumes **additivity** of component effects (often no component interaction) —
  state and, where possible, check it (an interaction/full-interaction model). Define the reference
  and the estimand (relative effects via a league table) explicitly.

## Reporting (PRISMA-NMA)

Report against **PRISMA-NMA**: the network plot, the transitivity assessment, the consistency
evaluation (global + local), the ranking with its caveats, comparison-adjusted publication-bias
assessment, and estimate-level certainty. Risk of bias uses **RoB-NMA**; pairwise items via the SR/MA
workflow. Review-side probes: NM1–NM8 in `network_meta_analysis.md`.

# NHIS Claims-Based Disease Definition Guide

Methodological patterns for defining exposures, outcomes, and comorbidities
using ICD-10 codes in Korean NHIS claims data. This guide covers the validated
algorithm patterns -- NOT specific code lists for individual diseases.

---

## When to Use

- Defining study populations, exposures, or outcomes in NHIS cohort studies
- Building claims-based algorithms for disease ascertainment
- Validating disease definitions using multiple data sources within NHIS
- NOT for: survey data (KNHANES/NHANES), clinical registry data

---

## Core Principle: Claims-Based Algorithms

A single ICD-10 claim is insufficient for disease ascertainment due to rule-out
diagnoses, coding errors, and provisional codes. Validated algorithms combine
multiple data elements to improve positive predictive value.

### Algorithm Components

| Component | Source in NHIS | Purpose |
|-----------|---------------|---------|
| **ICD-10 diagnostic codes** | Inpatient + outpatient claims | Primary disease identification |
| **Medication prescriptions** | Pharmacy dispensing records | Confirms active treatment |
| **Procedure codes** | Claims procedure fields | Confirms diagnostic workup or treatment |
| **Health examination results** | National health exam data | Objective measurement (lab values, vitals) |
| **Visit frequency** | Claim count within time window | Distinguishes incident from rule-out |

---

## Validated Algorithm Patterns

### Pattern 1: N-Claim Rule (Most Common)

Require N or more claims with the target ICD-10 code within a defined time
window to establish a diagnosis.

```
Disease = (ICD-10 code X appears >= N times within T months)
```

| Stringency | Claims | Window | PPV | Use case |
|------------|--------|--------|-----|----------|
| Lenient | >= 1 claim | Any | Lower | Screening, sensitivity-focused |
| Standard | >= 2 claims | 12 months | Moderate | Most cohort studies |
| Strict | >= 3 claims | 12 months | Higher | Sensitivity analysis |

**Reporting template**:
"[DISEASE] was identified by at least [N] claims with ICD-10 codes [CODE
RANGE] within [T] months during the observation period."

### Pattern 2: Claim + Medication (Recommended for Chronic Diseases)

Combine diagnostic codes with disease-specific medication prescriptions.

```
Disease = (ICD-10 code X >= 1 time) AND (related medication dispensed >= 1 time)
```

**Reporting template**:
"[DISEASE] was defined as at least one claim with ICD-10 codes [CODE RANGE]
combined with a prescription for [MEDICATION CLASS] (ATC codes [RANGE])."

**Common medication validation pairs**:

| Disease category | ICD-10 pattern | Medication confirmation |
|-----------------|----------------|----------------------|
| Hypertension | I10-I15 | Antihypertensive agents (ATC C02-C09) |
| Diabetes mellitus | E10-E14 | Antidiabetic agents (ATC A10) |
| Hyperlipidemia | E78 | Lipid-modifying agents (ATC C10) |
| Coronary artery disease | I20-I25 | Antiplatelet + statin combination |
| Depressive disorders | F32-F33 | Antidepressants (ATC N06A) |
| Asthma | J45 | Inhaled corticosteroids / bronchodilators |

### Pattern 3: Look-Back Period for Incident Cases

To identify **incident** (new-onset) cases, require a disease-free washout
period before the index date.

```
Incident case = (ICD-10 code X after index date)
              AND (NO claims with code X during look-back period)
```

| Look-back | Strength | Trade-off |
|-----------|----------|-----------|
| 1 year | Minimum acceptable | May misclassify prevalent cases as incident |
| 3 years | Recommended | Better incident case identification |
| 5+ years | Strongest | Reduces sample size due to data availability |

**Reporting template**:
"To ensure the identification of incident cases, we applied a minimum [N]-year
look-back period and excluded individuals with any diagnostic claims for
[DISEASE] during this period."

### Pattern 4: Outcome Validation with Time Windows

For complications or sequelae, require temporal proximity to the index event.

```
Complication = (ICD-10 code Y within T1-T2 days after event X)
```

**Examples of time-window patterns**:

| Complication type | Window | Rationale |
|------------------|--------|-----------|
| Post-procedural complication | 0-30 days | Immediate perioperative period |
| Delayed complication | 30-90 days | Sub-acute period |
| Chronic sequela | 90 days - 1 year | Excludes acute phase, captures chronicity |
| Hospital admission for condition | Within 30 days of onset | Links admission to the triggering event |

### Pattern 5: Hierarchical Disease Classification

When exposure or outcome categories overlap, use a hierarchical scheme.

```
Level 1 (broadest): Any code in ICD-10 block (e.g., E10-E14 for all diabetes)
Level 2: Specific subcategory (e.g., E11 for type 2 diabetes only)
Level 3 (strictest): Subcategory + medication + lab confirmation
```

Sensitivity analyses should test multiple levels to assess robustness.

---

## ICD-10-KM vs ICD-10 Differences

Korean NHIS uses ICD-10-KM (Korean Modification), which is based on ICD-10
but includes Korea-specific extensions:

| Feature | ICD-10 (WHO) | ICD-10-KM |
|---------|-------------|-----------|
| Base structure | Same alphanumeric codes | Same base + Korean extensions |
| Additional codes | -- | U-codes for Korean-specific conditions |
| Version alignment | WHO updates | HIRA publishes annual Korean versions |
| Medication mapping | ATC codes | ATC codes (same international standard) |

For most major disease categories (cardiovascular, metabolic, respiratory,
neurological), the ICD-10-KM codes are identical to international ICD-10.
Country-specific extensions primarily affect rare or culture-bound conditions.

---

## Charlson Comorbidity Index in NHIS

The CCI is computed from claims using the Quan adaptation for ICD-10:

| CCI Category | Standard grouping |
|-------------|-------------------|
| 0 | No comorbidities |
| 1 | Single comorbidity point |
| >= 2 | Multiple comorbidities |

Always specify which CCI adaptation was used (Quan 2005 is standard for
ICD-10 claims data).

---

## Covariate Definitions from NHIS Data Sources

NHIS integrates three data sources, each providing different covariate types:

### Source 1: Claims Database
- Comorbidities (via ICD-10 codes)
- Medication history (via dispensing records)
- Healthcare utilization (visit counts, hospitalization days)

### Source 2: National Health Examination
- BMI, blood pressure, fasting glucose, GFR
- Liver enzymes (AST, ALT, gamma-GTP)
- Hemoglobin, total cholesterol
- Urinalysis

### Source 3: Health Interview / Questionnaire
- Smoking status (never / former / current)
- Alcohol consumption (frequency categories)
- Physical activity (MET-based sufficiency)
- Region of residence, household income level

**Baseline alignment**: Covariates must be anchored to the most recent data
available **prior to** the index date. Do not use post-index measurements.

---

## Sensitivity Analysis Patterns

### Stricter Disease Definition
Run the primary analysis with the standard algorithm, then repeat with a
more stringent definition (e.g., 2-claim → 3-claim, or adding medication
confirmation).

### Negative Control Outcomes
Select outcomes with no plausible biological link to the exposure (e.g.,
neuromuscular junction disorders, self-harm) to detect residual confounding
or coding artifacts.

### Negative Control Exposures
Use an exposure with no expected effect on the outcome as a falsification test.

### Multiple Look-Back Periods
Test 1-year, 3-year, and 5-year look-back periods and compare incidence rates
to assess the impact of prevalent case misclassification.

---

## Reporting Checklist

When reporting NHIS claims-based definitions in a manuscript:

- [ ] ICD-10 code ranges specified (main text or supplementary table)
- [ ] Number of claims required (1 vs 2+ vs with medication)
- [ ] Time window for claims (within 12 months, any time, etc.)
- [ ] Look-back period for incident case definition
- [ ] Medication confirmation specified with ATC codes (if used)
- [ ] Data source specified (claims vs exam vs questionnaire)
- [ ] CCI adaptation cited (Quan 2005 or other)
- [ ] Sensitivity analysis with alternative definition included

---

## Common Reviewer Flags

1. Single-claim definition without justification (low PPV concern)
2. No look-back period for "incident" cases
3. Mixing ICD-10 and ICD-9 codes without noting transition dates
4. Not specifying whether inpatient, outpatient, or both claims were used
5. CCI version not cited
6. Covariate timing not anchored to index date
7. Missing sensitivity analysis with stricter/looser definitions
8. Rule-out diagnoses not addressed (diagnostic codes without treatment)

---

## References (Methodological)

- Quan H et al. Coding algorithms for defining comorbidities in ICD-9-CM
  and ICD-10 administrative data. Med Care. 2005;43(11):1130-9.
- Choi JY et al. Validation of administrative data for the identification
  of chronic diseases in South Korea. J Korean Med Sci. (methodology reference)
- NHIS data user guide: nhiss.nhis.or.kr (Korean National Health Insurance
  Sharing Service)

# Polygenic Risk Score (PRS / PGS) Analysis Guide

Building and validating a genome-wide weighted allele sum as a predictor or risk-stratifier. The
score is easy to compute; what fails review is **ancestry transferability, base/target leakage,
incremental value over established clinical risk, calibration, and the discrimination-vs-utility
gap**. This is the analysis-side companion to review probes **PG1–PG8** in
`polygenic_risk_score.md`.

---

## When to Use

- Developing, validating, or applying a PRS/PGS as a predictor or risk-stratifier.
- NOT for: the instrumental-variable use of genetics (causal inference) — see
  `mendelian_randomization.md`; a generic non-genetic prediction model — see `regression.md` /
  the clinical-prediction-model probes.

## Construction & data hygiene

- **Base (discovery GWAS)** and **target (validation)** samples must be **independent** — overlap
  inflates performance. Reuse a **PGS Catalog** score (cite the PGS ID) or document the weights,
  variant set, genome build, allele/strand alignment, and imputation/missing-genotype handling.
- **Tuning** (P+T threshold, LDpred/PRS-CS shrinkage or proportion-of-causal-variants, the chosen
  quantile cut) is selected on a **separate tuning set** and the final score evaluated **out-of-
  sample** — never tune and report performance in the same data (overfitting / winner's curse).
- Tools: `PRSice-2`, `LDpred2` (`bigsnpr`), `PRS-CS`/`PRS-CSx`, `SBayesR`; for cross-ancestry,
  `PRS-CSx` / `BridgePRS`.

## Ancestry portability (the central issue)

- PRS transfers poorly across ancestries (LD, allele-frequency, and GxE differences). **Report
  performance separately per target ancestry**; differences *within* a broad ancestry group can be
  as large as across continental groups, so a single per-group number can mislead.
- Prefer **ancestry-matched or multi-ancestry discovery** GWAS; control population stratification
  with principal components. Do not extend a European-derived score to other ancestries without
  per-ancestry validation (equity harm).

## Effect-size, discrimination, calibration

- Standardize the PRS; report **OR/HR per SD** (with CI) and **quantile stratification** (decile/
  percentile, reference group stated) with **absolute risk** by stratum for a clinical claim.
- **Discrimination**: C-statistic / AUC with CI. **Calibration is separate** — calibration plot
  (observed vs expected by stratum), slope/intercept, in the **target population** (and re-checked
  across ancestries/cohorts where baseline incidence differs). A well-discriminating score can be
  badly miscalibrated for absolute risk.

## Incremental value over established clinical risk (the clinical crux)

- A clinical claim needs the PRS **on top of the guideline clinical model** (SCORE2 / QRISK3 /
  Pooled Cohort Equations / Tyrer-Cuzick / FRAX), not PRS-alone AUC: report **ΔC-statistic (with
  CI)**, **NRI / IDI**, and ideally **net benefit (decision curve)** vs the clinical model alone.

```r
# nested models: clinical vs clinical + PRS
m_clin <- coxph(Surv(time, event) ~ score2, data = d)
m_both <- coxph(Surv(time, event) ~ score2 + scale(prs), data = d)
# Delta C-index (Uno's C) with CI; NRI/IDI via survIDINRI / nricens; decision curve via dcurves
```

## Screening / stratification utility ≠ discrimination

- A **population-screening** claim needs the **detection rate at a fixed false-positive rate** (or
  the likelihood ratio for a PRS quantile) and the **number needed to screen / absolute risk
  difference** — AUC and HR-per-SD can look favourable while the detection-rate-at-acceptable-FPR
  is poor.

## Design caveat

- Prefer **prospective/incident** validation; a cross-sectional **case–control prevalent-disease**
  association overstates predictive utility (PRS predicts incident disease less well than it
  associates with prevalent disease).

## Reporting (PGS-RS / TRIPOD+AI)

Report against the **PGS Reporting Standards** and TRIPOD+AI: development and validation samples
with **ancestry composition**, score-construction provenance (PGS Catalog ID), and the full
performance set (discrimination, **calibration**, **incremental value**, and screening operating
characteristics where a screening claim is made). Scale the conclusion to the evidence — a
quantile relative-risk gradient is not, by itself, demonstrated clinical actionability. Review-side
probes: PG1–PG8 in `polygenic_risk_score.md`.

# Propensity Score Analysis Guide

Methods for estimating causal treatment effects in observational studies by balancing
confounders between treatment groups.

---

## When to Use

- Observational study comparing treatment vs control (or two interventions)
- Multiple confounders to adjust for
- Goal: estimate causal effect analogous to an RCT
- NOT for: randomized trials (already balanced), single-arm studies

---

## Estimands — Choose Before Analysis

| Estimand | Definition | Target population | Method |
|----------|-----------|-------------------|--------|
| **ATE** | Average Treatment Effect | Entire study population | IPTW |
| **ATT** | Effect on Treated | Treatment group only | PSM, ATT weighting |
| **ATO** | Effect on Overlap population | PS 0.2-0.8 region | Overlap weighting |

Comparing PSM and IPTW results directly is inappropriate — they estimate different estimands.

---

## Step-by-Step Workflow

### Step 1: PS Estimation
- Model: logistic regression (standard)
- Dependent variable: treatment assignment (binary)
- Covariates: all variables that affect the outcome (confounders)
- Do NOT include: instrumental variables (affect treatment but not outcome)
- Individual PS model coefficients have no clinical meaning — only PS distribution matters

### Step 2: Apply PS Method

**Option A — PS Matching (PSM)**
- Nearest-neighbor matching with caliper = 0.2 x SD(logit PS)
- 1:1 matching is standard; 1:N or full matching available
- Estimand: ATT (typically)
- Drawback: unmatched subjects are excluded → sample size reduction

**Option B — IPTW (Inverse Probability of Treatment Weighting)**
- Weights: treated = 1/PS, control = 1/(1-PS) for ATE
- **Always use stabilized weights** to prevent extreme values
- Stabilized: treated = P(T=1)/PS, control = P(T=0)/(1-PS)
- All subjects included (no exclusion)
- Flag extreme weights > 10

**Option C — SIPTW (Stabilized Inverse Probability of Treatment Weighting)**
- Weights: treated = P(T=1)/PS, control = P(T=0)/(1-PS)
- Mathematically equivalent to stabilized IPTW (Option B with stabilized=True)
- Named explicitly as SIPTW in emulated target trial literature
- Maintains entire cohort sample size (no exclusion)
- Allows appropriate variance estimation of main effect
- Estimand: ATE
- Increasingly used in large-scale NHIS cohort studies
- Report effective sample size (ESS) alongside raw N

**Option D — Overlap Weighting (Recommended for most cases)**
- Weights: treated = (1-PS), control = PS
- Naturally down-weights subjects at PS extremes
- No extreme weight problem (advantage over IPTW)
- Estimand: ATO
- Increasingly recommended in recent guidelines (JAMA 2020, AJE 2024)

### Step 3: Assess Balance
- **SMD < 0.10** for all covariates (Austin, 2009)
- SMD < 0.25 is acceptable but suboptimal
- Use SMD, NOT p-values (SMD is sample-size independent)
- **Love plot**: pre/post-matching SMD comparison (mandatory figure)
- Variance ratios: should be within 0.5-2.0
- PS distribution overlap: histogram comparing treated vs control

### Step 4: Outcome Analysis
- **After PSM**: paired analysis (paired t-test, conditional logistic, stratified Cox)
- **After IPTW/OW**: weighted regression using survey methods
- Always use robust/sandwich standard errors for weighted analyses

### Step 5: Sensitivity Analysis
- **E-value**: quantifies how strong unmeasured confounding would need to be to explain away the result
- Report E-value for the point estimate and lower CI bound

---

## Balance Table (Required Output)

| Variable | Before matching | | After matching | |
|----------|------|------|------|------|
| | Treated | Control | SMD | Treated | Control | SMD |
| Age, mean (SD) | 65.2 (12.1) | 58.7 (14.3) | 0.49 | 62.1 (11.8) | 61.8 (12.0) | 0.03 |

---

## Reporting Templates

**PSM**: "Propensity scores were estimated using logistic regression with the following covariates: [list]. PS matching was performed using 1:1 nearest-neighbor matching with a caliper of 0.2 SD of the logit PS. After matching, all SMDs were below 0.10 (Figure X). In the matched cohort (n = X pairs), ..."

**IPTW/OW**: "Inverse probability of treatment weighting (or overlap weighting) was applied using stabilized weights. Covariate balance was assessed using SMDs (all < 0.10; Figure X). The weighted analysis showed ..."

**SIPTW**: "Stabilized inverse probability of treatment weighting was used to balance covariate distributions between the [exposed] and [unexposed] groups. This approach maintains the sample size of the entire cohort and allows for appropriate estimation of the variance of the main effect. Covariate balance was assessed using SMDs (all < 0.10; Figure X)."

---

## Common Reviewer Flags

1. SMD not reported (using p-values instead)
2. Caliper width not specified
3. Estimand (ATE/ATT/ATO) not stated
4. No sensitivity analysis for unmeasured confounding
5. Number of unmatched subjects not reported (for PSM)
6. Extreme weights not assessed (for IPTW)
7. Individual PS model coefficients interpreted clinically

---

## Python Packages
- `sklearn.linear_model.LogisticRegression` — PS estimation
- `causalinference` — matching (limited)
- Manual implementation for IPTW/OW (see template)
- `statsmodels` — weighted regression

## R Packages
- `MatchIt` — PS matching
- `WeightIt` — IPTW, overlap weighting
- `cobalt` — balance assessment, Love plot
- `survey` — weighted outcome analysis
- `tableone` — baseline table with SMD
- `EValue` — sensitivity analysis

# Regression Analysis Guide

Covers both logistic regression (binary outcome) and multiple linear regression (continuous outcome).

---

## Logistic Regression

### When to Use
Binary outcome variable (0/1, event/no event) with one or more predictors.
Used for risk factor identification and prediction model building.

### Assumptions
1. Binary outcome
2. Linear relationship between continuous predictors and log-odds (Box-Tidwell test)
3. Independent observations (if repeated → GEE or mixed logistic)
4. No multicollinearity: VIF < 5
5. Sufficient sample: EPV >= 10 (minimum), >= 20 (recommended)

### Variable Selection
- **Clinical rationale first** — avoid purely data-driven stepwise selection
- Include: variables significant at P < 0.10 in univariable analysis + known confounders
- STROBE/TRIPOD guideline compliance required

### Model Assessment

**Calibration:**
- Hosmer-Lemeshow test: P > 0.05 = adequate fit (overpowered with large N)
- Calibration plot: predicted probability vs observed frequency

**Discrimination:**
- C-statistic (= AUC): 0.7-0.8 acceptable, 0.8-0.9 excellent, > 0.9 outstanding

**Multicollinearity:**
- VIF > 5 → remove or combine variables

### Required Outputs
1. OR table: univariable AND multivariable OR (95% CI), P-value per variable
2. C-statistic with 95% CI
3. Hosmer-Lemeshow test result
4. VIF table (supplementary)
5. Box-Tidwell results for continuous predictors (supplementary)

### OR Table Format
| Variable | Univariable OR (95% CI) | P | Multivariable OR (95% CI) | P |
|----------|------------------------|---|--------------------------|---|
| Age (per 10 yr) | 1.45 (1.12-1.88) | 0.005 | 1.32 (1.01-1.73) | 0.042 |

### Reporting Template
"Multivariable logistic regression was performed to identify independent predictors of [outcome]. Variables with P < 0.10 in univariable analysis and clinically relevant confounders were included. The model showed good discrimination (C-statistic = 0.82, 95% CI 0.78-0.86) and calibration (Hosmer-Lemeshow P = 0.45). [Variable] was independently associated with [outcome] (adjusted OR = 2.15, 95% CI 1.43-3.24; P < 0.001)."

### Pitfalls
- OR != RR: when event rate > 10%, OR overestimates RR
- EPV < 10 → overfitting risk. Consider penalized regression (LASSO/Ridge)
- Specify reference category for categorical variables
- Specify units for continuous variables (per 1 year vs per 10 years)
- Complete separation → use Firth's penalized likelihood

---

## Multiple Linear Regression

### When to Use
Continuous outcome variable with one or more predictors.
Used for identifying determinants and estimating adjusted effects.

### Assumptions (LINE + No Multicollinearity)
1. **L**inearity: residuals vs fitted plot
2. **I**ndependence: no repeated measures (if repeated → LMM/GEE)
3. **N**ormality of residuals: Q-Q plot, Shapiro-Wilk on residuals
4. **E**qual variance (homoscedasticity): residuals vs fitted, Scale-Location plot
5. No multicollinearity: VIF < 5
6. No influential outliers: Cook's distance < 4/n

### Assumption Violations → Alternatives
| Violation | Fix |
|-----------|-----|
| Non-linearity | Log transform, polynomial terms, GAM |
| Heteroscedasticity | Robust SE, WLS |
| Non-normal residuals | Transform outcome, bootstrap CI |
| Multicollinearity | Remove variable, combine, Ridge/LASSO |

### Model Evaluation
- R² (coefficient of determination): proportion of variance explained
- Adjusted R²: penalized for number of predictors — use for model comparison
- In medical research, R² = 0.2-0.4 can be meaningful (high biological variability)

### Diagnostic Plots (4-panel, always generate)
1. Residuals vs Fitted → linearity + homoscedasticity
2. Q-Q plot → residual normality
3. Scale-Location → homoscedasticity
4. Residuals vs Leverage → influential outliers (Cook's distance)

### Required Outputs
1. Coefficient table: β (95% CI), P-value per variable
2. R² and adjusted R²
3. VIF table
4. 4-panel diagnostic plot
5. Standardized coefficients (optional, for effect size comparison)

### Coefficient Table Format
| Variable | β (95% CI) | P |
|----------|-----------|---|
| Age (per year) | 0.45 (0.32 to 0.58) | < 0.001 |
| Model R² | 0.35 | |
| Adjusted R² | 0.33 | |

### Reporting Template
"Multiple linear regression was performed with [outcome] as the dependent variable. The model explained X% of the variance (adjusted R² = 0.XX). After adjusting for [covariates], [variable] was significantly associated with [outcome] (β = X.XX, 95% CI X.XX to X.XX; P = exact). Model assumptions were verified using diagnostic plots."

### Pitfalls
- Always report β units (per 1 year, per 10 kg/m², etc.)
- Standardized β useful for comparing effect sizes but unstandardized is standard in papers
- EPV for continuous outcome: N >= 10-20 per predictor
- Always present diagnostic plots (at minimum in supplementary)

# Repeated Measures / Mixed Models Guide

Analysis methods for longitudinal data where the same subjects are measured at multiple time points.

---

## When to Use

- Same subjects measured at 2+ time points
- Key research question: does change over time differ between groups? (Time x Group interaction)
- Examples: treatment response over weeks, serial imaging measurements, before/after/follow-up

---

## Method Selection

| Condition | Recommended method |
|-----------|-------------------|
| No missing data + few time points + sphericity met | RM ANOVA |
| Missing data (MAR) + continuous outcome | **LMM (preferred)** |
| Missing data + binary/count outcome | GEE (or GLMM) |
| Individual trajectory estimation needed | LMM (random slope) |
| Population-averaged effect only | GEE |
| Non-normal outcome | GLMM or GEE |

**Default recommendation: LMM** — handles missing data, does not require sphericity, allows unequal time spacing.

---

## Data Format

All methods (LMM, GEE) require **long format**:

| id | time | group | outcome |
|----|------|-------|---------|
| 1 | 0 | A | 45.2 |
| 1 | 1 | A | 42.8 |
| 1 | 2 | A | 38.1 |

Convert from wide format before analysis.

---

## 1. RM ANOVA

### Sphericity Check (Mandatory)
- **Mauchly's test**: H0 = sphericity holds
  - P >= 0.05 → use standard F-test
  - P < 0.05 → apply correction

### Corrections for Sphericity Violation
1. **Greenhouse-Geisser (G-G)**: conservative, recommended
2. **Huynh-Feldt**: slightly liberal
3. **Multivariate test** (Pillai's trace): no sphericity assumption needed

### Limitations
- Complete cases only — any missing time point drops the entire subject
- Difficult to add covariates
- Equal time spacing required

### Key Results to Report
- Mauchly's test result (W statistic, P-value, epsilon)
- Correction method used (if sphericity violated)
- Within-subject effect: Time (F, df, P)
- Between-subject effect: Group (F, df, P)
- **Interaction: Time x Group** (most important)

---

## 2. Linear Mixed Model (LMM)

### Structure
```
Y = Xβ (fixed effects) + Zb (random effects) + ε
```

### Random Effects Selection
- **Random intercept only**: subjects differ in baseline (default, start here)
- **Random intercept + slope**: subjects differ in both baseline and rate of change
- Decision: visualize individual trajectories (spaghetti plot). If slopes vary → add random slope.
- If model does not converge with random slope → simplify to random intercept only

### Covariance Structure Selection
| Structure | Property | When to use |
|-----------|----------|-------------|
| **CS** (Compound Symmetry) | Equal correlation between all time pairs | Equivalent to RM ANOVA |
| **AR(1)** | Correlation decays with time lag | Equally spaced measurements |
| **UN** (Unstructured) | Free correlation for each pair | Few time points only (many parameters) |

Compare structures using **AIC/BIC** (lower = better).

### Advantages over RM ANOVA
- Handles missing data under MAR assumption
- No sphericity requirement
- Allows unequal time intervals
- Easy covariate adjustment
- Estimates individual trajectories

---

## 3. GEE (Generalized Estimating Equations)

### Key Properties
- Estimates **population-averaged (marginal) effects** — not individual-level
- Specify working correlation: exchangeable, AR(1), unstructured
- Works for non-normal outcomes (binary, count)
- Requires **MCAR** assumption for missing data (or use weighted GEE for MAR)
- Needs sufficient clusters: subjects >= 30-40

### LMM vs GEE Decision
| Question | LMM | GEE |
|----------|-----|-----|
| "How does each patient change?" | Yes | No |
| "How does the group average change?" | Yes | Yes |
| Binary/count outcome without GLMM | No | Yes |
| MAR missing data tolerance | Yes | No (MCAR only) |

---

## Required Outputs

1. **Spaghetti plot**: individual trajectories by group
2. **Model summary table**: fixed effects (β, 95% CI, P), random effects variance
3. **Time x Group interaction** result prominently reported
4. Model fit: AIC/BIC for covariance structure comparison
5. Missing data description (n per time point, mechanism assumed)

---

## Reporting Templates

**RM ANOVA**: "Repeated-measures ANOVA was performed with Greenhouse-Geisser correction for violation of sphericity (Mauchly's test P < 0.001, ε = 0.42). There was a significant time × group interaction (F(X, Y) = Z, P = exact)."

**LMM**: "A linear mixed-effects model with random intercepts for subjects and [CS/AR(1)] correlation structure was fitted. The time × group interaction was significant (β = -2.34, 95% CI -3.87 to -0.81; P = 0.003), indicating that the rate of change in [outcome] differed between groups."

**GEE**: "GEE with exchangeable correlation structure was used to estimate population-averaged effects. The time × group interaction was ..."

---

## Common Reviewer Flags

1. Sphericity test not reported (for RM ANOVA)
2. Covariance structure selection rationale not stated (for LMM)
3. Time × Group interaction not interpreted
4. Missing data handling not described
5. Using RM ANOVA with substantial missing data (should use LMM)
6. Not reporting individual trajectories (spaghetti plot)

---

## Python Packages
- `pingouin` — RM ANOVA with G-G correction
- `statsmodels.formula.api.mixedlm` — LMM
- `statsmodels.genmod.generalized_estimating_equations.GEE` — GEE

## R Packages
- `lme4` + `lmerTest` — LMM (standard)
- `nlme` — LMM with correlation structures (CS, AR1, UN)
- `geepack` — GEE
- Base R `aov()` with `Error()` — RM ANOVA

# Survey-Weighted Analysis Guide

Methods for analyzing complex survey data (KNHANES, NHANES, KCHS, and similar
nationally representative health surveys) that use stratified, multistage
probability sampling designs.

---

## When to Use

- Data from a national health survey with sampling weights (e.g., KNHANES, NHANES, KCHS)
- Goal: produce nationally representative prevalence, means, or associations
- Cross-national comparisons using parallel survey datasets
- NOT for: simple random samples, census data, or claims-based cohorts (NHIS, JMDC)

Claims-based databases (NHIS, JMDC) are NOT surveys -- they do not have sampling
weights or complex sampling design. Use standard regression for these.

---

## Key Concepts

### Complex Survey Design Elements

| Element | Description | Survey variable |
|---------|-------------|-----------------|
| **Stratification** | Groups the population into non-overlapping strata before sampling | `strata` |
| **Clustering (PSU)** | Primary sampling units within strata (e.g., districts, census blocks) | `cluster` / `PSU` |
| **Sampling weights** | Inverse probability of selection, adjusted for non-response and post-stratification | `weight` |

Ignoring these elements produces biased standard errors, incorrect p-values, and
non-representative point estimates.

### Dataset-Specific Design Variables

| Dataset | Strata variable | Cluster/PSU variable | Weight variable | Notes |
|---------|----------------|---------------------|-----------------|-------|
| **KNHANES** | `kstrata` | `psu` | `wt_itvex` (interview+exam) or `wt_ntr` (nutrition) | Years may be non-consecutive (PHQ-9 only in certain cycles) |
| **NHANES** | `SDMVSTRA` | `SDMVPSU` | `WTMECXYR` (exam) or `WTINTXYR` (interview) | 2-year cycles; combine cycles with adjusted weights |
| **KCHS** | varies by year | varies by year | `wt` | Annual community survey; single-stage cluster design |

### Weight Selection Rules

- **Interview-only variables**: use interview weight
- **Exam/lab variables**: use exam weight (smaller denominator)
- **Nutrition variables (KNHANES)**: use nutrition weight
- **Multi-cycle NHANES**: divide weight by number of cycles combined (e.g., 4 cycles: weight/4)
- **Single-cycle analysis**: use the cycle-specific weight as-is

### Subpopulation (domain) analysis — never row-delete

A restricted analysis (adults only, one sex, a disease subgroup) must keep the **full design** and select the domain, **not** filter the data frame and refit. Row-deletion discards the strata/PSU structure of the dropped units and gives wrong standard errors and design degrees of freedom.

- R `survey`: `subset(design, age >= 18)` on the **design object** (or `svyby`), not `svydesign(data = df[df$age>=18, ])`.
- Stata: `svy, subpop(if age>=18):` — never `keep if age>=18` before `svy:`.
- Python `samplics` / R is preferred; statsmodels has no native domain estimator.

### Reporting & common errors (these invalidate the inference, flag at review)

- **Model-based SEs on weighted points.** Applying the weight but computing SEs without strata + PSU (or replicate weights) understates uncertainty. Always declare `strata` + `id`/`cluster`.
- **Weighted total ≠ sample size.** Report the **unweighted n** as the analytic sample; the weighted figure is a *population* estimate, not "n".
- **Design effect / effective n.** Report DEFF or the effective sample size where precision is load-bearing; a large DEFF means far fewer independent observations than rows.
- **Unweighted vs weighted divergence.** If they differ materially, that signals weight-dependent selection — discuss it, do not hide it.
- **Data-driven thresholds.** A restricted-cubic-spline "non-linear/saturation" claim needs a pre-specified non-linearity test (LRT/Wald vs the linear model), and any "inflection point / threshold" from a recursive breakpoint search must carry a **confidence interval** and a stability check — a searched cutoff is exploratory, not a validated target (review-side probe O12 in `observational_confounding.md`).

---

## Step-by-Step Workflow

### Step 1: Declare Survey Design

Always declare the design before any analysis. This ensures correct variance estimation.

**Python (statsmodels)**:
```python
# statsmodels does not have a native survey design object.
# Use linearmodels or manual weight application.
# For publication-quality survey analysis, R is strongly recommended.
```

**R (survey package)**:
```r
library(survey)

# KNHANES
design_kr <- svydesign(
  id = ~psu,
  strata = ~kstrata,
  weights = ~wt_itvex,
  data = df_kr,
  nest = TRUE
)

# NHANES (2-year cycle)
design_us <- svydesign(
  id = ~SDMVPSU,
  strata = ~SDMVSTRA,
  weights = ~WTMECXYR,
  data = df_us,
  nest = TRUE
)
```

**SAS**:
```sas
/* KNHANES */
PROC SURVEYLOGISTIC DATA=kr;
  STRATA kstrata;
  CLUSTER psu;
  WEIGHT wt_itvex;
  MODEL outcome(event='1') = exposure covariates;
RUN;

/* NHANES */
PROC SURVEYLOGISTIC DATA=us;
  STRATA SDMVSTRA;
  CLUSTER SDMVPSU;
  WEIGHT WTMECXYR;
  MODEL outcome(event='1') = exposure covariates;
RUN;
```

### Step 2: Weighted Descriptive Statistics

**R**:
```r
# Weighted means
svymean(~continuous_var, design, na.rm = TRUE)

# Weighted proportions
svymean(~factor(categorical_var), design, na.rm = TRUE)

# Weighted Table 1 by group
library(tableone)
svyCreateTableOne(
  vars = c("age", "sex", "bmi", "income"),
  strata = "exposure_group",
  data = design,
  test = TRUE
)
```

**SAS**:
```sas
PROC SURVEYMEANS DATA=dataset;
  STRATA strata_var;
  CLUSTER cluster_var;
  WEIGHT weight_var;
  VAR continuous_var1 continuous_var2;
RUN;

PROC SURVEYFREQ DATA=dataset;
  STRATA strata_var;
  CLUSTER cluster_var;
  WEIGHT weight_var;
  TABLES group * categorical_var / CHISQ;
RUN;
```

### Step 3: Weighted Regression — Sequential Model Building

The standard cross-national analysis pattern uses sequential model building:

| Model | Covariates | Purpose |
|-------|-----------|---------|
| Model 1 | Age, sex | Minimal adjustment |
| Model 2 | Model 1 + income, education, smoking, alcohol, BMI, comorbidities | Full adjustment |

**R (survey-weighted logistic regression)**:
```r
# Model 1: age + sex
model1 <- svyglm(
  outcome ~ exposure + age + sex,
  design = design,
  family = quasibinomial()
)

# Model 2: full adjustment
model2 <- svyglm(
  outcome ~ exposure + age + sex + income + education +
            smoking + alcohol + bmi + cvd_history,
  design = design,
  family = quasibinomial()
)

# Extract weighted OR (wOR) with 95% CI
extract_wor <- function(model, var) {
  coef_val <- coef(model)[var]
  se_val <- summary(model)$coefficients[var, "Std. Error"]
  or <- exp(coef_val)
  ci_lo <- exp(coef_val - 1.96 * se_val)
  ci_hi <- exp(coef_val + 1.96 * se_val)
  p_val <- summary(model)$coefficients[var, "Pr(>|t|)"]
  data.frame(wOR = or, CI_lower = ci_lo, CI_upper = ci_hi, P = p_val)
}
```

**SAS (PROC SURVEYLOGISTIC)**:
```sas
/* Model 2: full adjustment */
PROC SURVEYLOGISTIC DATA=dataset;
  STRATA strata_var;
  CLUSTER cluster_var;
  WEIGHT weight_var;
  CLASS sex(ref='Male') income(ref='High') smoking(ref='Never') / PARAM=REF;
  MODEL outcome(event='1') = exposure age sex income education smoking alcohol bmi;
  ODDSRATIO exposure / CL=WALD;
RUN;
```

### Step 4: Subgroup / Stratified Analyses

```r
# Stratified by sex
svyglm(
  outcome ~ exposure + age + income + education + smoking + alcohol + bmi,
  design = subset(design, sex == "Male"),
  family = quasibinomial()
)
# Repeat for Female
# Note: exclude the stratification variable from covariates
```

Reporting pattern: "Weighted odds ratios are adjusted for all covariates except
for the stratification variable."

### Step 5: Advanced Analyses (Optional)

**Restricted cubic spline (dose-response)**:
```r
library(rms)
design_rms <- svydesign(id = ~psu, strata = ~kstrata,
                         weights = ~wt_itvex, data = df)
model_rcs <- svyglm(
  outcome ~ rcs(continuous_exposure, 3) + age + sex + covariates,
  design = design_rms,
  family = quasibinomial()
)
```

**Weighted quantile sum (WQS) regression**:
```r
library(gWQS)
# WQS for composite exposure (e.g., LE8 components)
result_wqs <- gwqs(
  outcome ~ wqs + age + sex + covariates,
  mix_name = c("comp1", "comp2", "comp3", "comp4"),
  data = df,
  q = 4,           # quartiles
  b = 500,         # bootstrap iterations
  seed = 42,
  family = "binomial",
  weights = df$weight_var
)
```

---

## Weighted SMD (Standardized Mean Difference)

For balance assessment in weighted analyses:

```r
library(survey)
library(tableone)

# SMD in survey design
tab <- svyCreateTableOne(
  vars = covariates,
  strata = "treatment",
  data = design,
  smd = TRUE
)
print(tab, smd = TRUE)
```

**Manual calculation (Python)**:
```python
def weighted_smd(x, treatment, weights, is_binary=False):
    """Calculate weighted standardized mean difference."""
    t_mask = treatment == 1
    c_mask = treatment == 0
    
    w1, w0 = weights[t_mask], weights[c_mask]
    x1, x0 = x[t_mask], x[c_mask]
    
    wm1 = np.average(x1, weights=w1)
    wm0 = np.average(x0, weights=w0)
    
    if is_binary:
        denom = np.sqrt((wm1 * (1 - wm1) + wm0 * (1 - wm0)) / 2)
    else:
        wv1 = np.average((x1 - wm1) ** 2, weights=w1)
        wv0 = np.average((x0 - wm0) ** 2, weights=w0)
        denom = np.sqrt((wv1 + wv0) / 2)
    
    return (wm1 - wm0) / denom if denom > 0 else 0.0
```

---

## Cross-National Analysis Pattern

When comparing two countries using parallel surveys:

1. **Never pool** raw data across countries into a single regression
2. Analyze each country **independently** with country-specific survey design
3. Present results **side-by-side** in the same table
4. Compare effect magnitudes narratively (not via interaction terms)

### Standard Output Table Format

| Variable | Korea (KNHANES) | | US (NHANES) | |
|----------|------|------|------|------|
| | Model 1 wOR (95% CI) | Model 2 wOR (95% CI) | Model 1 wOR (95% CI) | Model 2 wOR (95% CI) |
| Exposure | 1.42 (1.21-1.67) | 1.35 (1.14-1.59) | 1.28 (1.10-1.49) | 1.22 (1.04-1.43) |

---

## Reporting Templates

**Methods**:
"To account for the complex survey designs and ensure nationally representative
estimates, appropriate stratification, clustering, and sampling weights were
applied. Weighted logistic regression models were used to estimate weighted odds
ratios (wORs) and 95% confidence intervals (CIs). Model 1 adjusted for age and
sex, while Model 2 further included [covariates]. All analyses were performed
using SAS version 9.4 (SAS Institute Inc.) [and R version X.X.X (R Foundation
for Statistical Computing)]. A two-sided P value of less than 0.05 was
considered statistically significant."

**Results**:
"In the weighted analysis of [N] participants from [DATASET], [EXPOSURE] was
significantly associated with [OUTCOME] (wOR [X.XX]; 95% CI [X.XX-X.XX]) after
adjusting for [covariates] (Model 2)."

---

## Common Reviewer Flags

1. Survey weights not applied (unweighted analysis of survey data)
2. Strata/cluster variables not specified (incorrect SE estimation)
3. Wrong weight variable used (interview weight for lab variables)
4. Multi-cycle NHANES weights not adjusted (divided by number of cycles)
5. Data pooled across countries instead of analyzed separately
6. Weighted proportions not reported (using raw counts instead)
7. Subgroup analysis includes the stratification variable as covariate
8. Missing data handling not stated (complete case vs imputation)

---

## Python vs R Recommendation

| Task | Recommended | Reason |
|------|-------------|--------|
| Survey-weighted regression | **R (survey)** | Native support, correct variance estimation |
| Survey-weighted Table 1 | **R (tableone)** | `svyCreateTableOne()` handles design |
| WQS regression | **R (gWQS)** | Only available in R |
| Dose-response (RCS) | **R (rms + survey)** | Integrated with survey design |
| Quick descriptives | Python (statsmodels) | Adequate for simple weighted means |

For publication-quality survey analysis, **R is strongly recommended** over Python.
Python's statsmodels supports basic weighted regression but lacks full survey
design support (no strata/cluster specification for variance estimation).

---

## R Packages

- `survey` -- core survey design and analysis
- `tableone` -- survey-weighted baseline tables with SMD
- `rms` -- restricted cubic splines with survey design
- `gWQS` -- weighted quantile sum regression
- `srvyr` -- tidyverse-compatible survey analysis wrapper

## SAS Procedures

- `PROC SURVEYMEANS` -- weighted means and proportions
- `PROC SURVEYFREQ` -- weighted frequency tables and chi-square
- `PROC SURVEYLOGISTIC` -- weighted logistic regression (wOR)
- `PROC SURVEYREG` -- weighted linear regression
- `PROC SURVEYPHREG` -- weighted Cox proportional hazards

# Survival / Time-to-Event Guide

Estimating time-to-event outcomes and prognostic effects. The estimator is easy to
call; the ways these analyses fail review are (1) ignoring **competing risks** so a naive
1−KM **overestimates** the cumulative incidence, (2) reporting a **single time-averaged
hazard ratio** when the proportional-hazards assumption is violated, and (3) **estimand
drift** — quoting a subdistribution hazard for an etiologic claim or a cause-specific hazard
for an absolute-risk claim. This guide produces the right estimand; the operational caveats
(EPV gate, cluster-robust CIs, interval-censoring) live in the SKILL.md `### Survival
Analysis` section.

---

## When to Use

- **Single event type, right-censored** → Kaplan–Meier (unadjusted) + Cox proportional-hazards
  (adjusted HR with 95% CI); the log-rank test for a KM group comparison.
- **≥2 competing event types** (recurrence + competing death, or cause-specific mortality) →
  cumulative incidence functions + **cause-specific Cox** or **Fine–Gray** — see below.
- **Prognostic model discrimination** → a C-index variant matched to the censoring + a
  time-dependent AUC at a clinical horizon (S6).
- NOT for: events detected only at scheduled visits → interval-censored methods (SKILL.md);
  recurrent events per subject → a cluster-robust or frailty model (SKILL.md cluster rule);
  rater agreement → `analysis_guides/agreement_reliability.md`.

---

## Competing risks come first (this, not the HR, is the issue)

If a subject can experience an event that **precludes** the event of interest (death before
recurrence), treating the competing event as ordinary censoring is *informative* censoring:
the naive **1−KM overestimates** the cumulative incidence of the event of interest. Produce the
**cumulative incidence function (CIF)** with an Aalen–Johansen / Fine–Gray estimator instead.
This is the produce-side of probe **S3**.

```python
# lifelines: cumulative incidence for a competing-risks event of interest
from lifelines import AalenJohansenFitter, KaplanMeierFitter
import pandas as pd
df = pd.read_csv("survival.csv")   # time, event_type (0=censored, 1=interest, 2=competing)

ajf = AalenJohansenFitter()
ajf.fit(df["time"], df["event_type"], event_of_interest=1)
print(ajf.cumulative_density_.tail(1))          # correct CIF for cause 1

kmf = KaplanMeierFitter()                        # NAIVE (cause 2 censored) — overestimates
kmf.fit(df["time"], (df["event_type"] == 1).astype(int))
print("naive 1-KM:", float(1 - kmf.survival_function_.iloc[-1, 0]))
# the naive value will exceed the Aalen-Johansen CIF whenever the competing event is common.
```

**Which model for which question (S8 estimand):**

- **Cause-specific hazard** (a Cox model that censors the competing event) answers an
  **etiologic** question — "does X change the rate of recurrence among those still at risk?"
- **Fine–Gray subdistribution hazard** (`cmprsk::crr` / `survival::finegray` in R) answers a
  **prognostic / absolute-risk** question — "does X change the cumulative *incidence*?" Quote an
  sHR for an incidence claim and a cause-specific HR for an etiologic claim — not the reverse.

```r
library(survival)                       # Fine-Gray via a weighted Cox on the finegray split
fg  <- finegray(Surv(time, factor(event_type)) ~ ., data = d, etype = 1)
crr <- coxph(Surv(fgstart, fgstop, fgstatus) ~ x, weight = fgwt, data = fg)   # subdistribution HR
```

---

## Proportional hazards, and RMST when it fails

A single Cox HR is a time-average; if PH is violated it averages a changing effect. Test PH,
and when it fails report a **restricted mean survival time (RMST) difference** at a fixed
horizon (an estimand that stays interpretable under non-PH) rather than one HR.

```python
from lifelines import CoxPHFitter
cph = CoxPHFitter().fit(df, "time", "event")
cph.check_assumptions(df, p_value_threshold=0.05)   # Schoenfeld; significant -> PH violated
```

```r
library(survRM2)
rmst2(time, status, arm, tau = 3)   # RMST difference at tau=3 years, valid under non-PH
```

Do not report a single time-averaged HR alongside a significant Schoenfeld test without also
giving a piecewise/time-stratified HR or an RMST difference (SKILL.md PH-violation rule).

---

## Follow-up and discrimination (S6)

- **Reverse Kaplan–Meier median follow-up** — the honest "how long were people followed"
  (median event time answers a different question). Compute it by **swapping the event
  indicator** (censored observations become the "events"); report it per cohort and per
  outcome, with the censoring date.

```python
kmf.fit(df["time"], 1 - df["event"])            # swap: censored -> event
print("reverse-KM median follow-up:", kmf.median_survival_time_)
```

- **C-index variant** — Harrell's C is biased under heavy or non-random censoring; report
  **Uno's IPCW C** (`survC1::Est.Cval` / `timeROC`) and a **time-dependent AUC at a clinical
  horizon** (2-/3-year) beside it, and state which variant and horizon (S6).

---

## Estimand provenance (S8)

State the survival estimand explicitly and hold it consistent across Abstract / Methods /
Results — event-free survival vs cause-specific cumulative incidence vs all-cause mortality,
and subject vs population level — with the evaluation horizon fixed in advance. Do not
re-designate the primary endpoint, model, or horizon after seeing results, and make every
derived statistic (an E-value, an sHR-vs-cause-specific contrast) trace to the *declared
primary* estimand. The self-review skill automates the registration ↔ manuscript and E-value
arithmetic checks (Phase 2.5f); see also `estimand-provenance-lock`.

---

## Reporting

- KM curves **with a number-at-risk table**; median survival with 95% CI (or the reason it is
  not reached).
- Cox HR (95% CI) with the **PH check stated**; for competing risks, the CIF and whether the
  HR is cause-specific or subdistribution.
- Events and person-time per group; the **reverse-KM median follow-up**; EPV for the model.
- The estimand and horizon, stated once and consistent everywhere.

---

## Common failures (flag at review)

- **Competing risks ignored** — naive 1−KM (or a cause-specific model presented as absolute
  risk) overestimates incidence; report a CIF and name cause-specific vs subdistribution (S3/S8).
- **A single time-averaged HR under a violated PH assumption** — needs a piecewise HR or RMST.
- **Harrell's C under heavy censoring** with no Uno/IPCW variant and no horizon (S6).
- **Median *survival* reported as "follow-up"** instead of the reverse-KM follow-up.
- **Estimand drift** — a primary endpoint/model/horizon re-designated post-hoc, or a derived
  statistic quoted off a non-primary estimate (S8).

---

## Anti-Hallucination

- Never hand-type an HR, median, CIF, or CI — compute it from the time-to-event CSV with a
  seeded script, and carry each estimate together with its CI.
- Do not report a Cox HR without the proportional-hazards check the code actually ran.
- Under competing risks, do not quote a 1−KM cumulative incidence — report the Aalen–Johansen /
  Fine–Gray CIF the code produced.

# Statistical Test Selection Guide

Decision tree for selecting the appropriate statistical test based on data structure.
Reference: Petrie & Sabin flowchart, Kirkwood Summary Guide.

---

## Step 1: Outcome Variable Type

| Outcome type | Next step |
|---|---|
| Continuous (measurement, score) | Step 2A |
| Binary / Categorical | Step 2B |
| Time-to-event (survival) | → Survival analysis type |
| Agreement / Reliability | → Inter-rater Agreement type |
| Diagnostic accuracy (Se/Sp/AUC) | → Diagnostic Accuracy type |

---

## Step 2A: Continuous Outcome

| Groups | Pairing | Normal? | Test | analyze-stats type |
|--------|---------|---------|------|--------------------|
| 1 | - | Yes | One-sample t-test | Group Comparison |
| 1 | - | No | Wilcoxon signed-rank (one-sample) | Group Comparison |
| 2 | Independent | Yes | Independent t-test | Group Comparison |
| 2 | Independent | No | Mann-Whitney U | Group Comparison |
| 2 | Paired | Yes | Paired t-test | Group Comparison |
| 2 | Paired | No | Wilcoxon signed-rank | Group Comparison |
| 3+ | Independent | Yes | One-way ANOVA + post-hoc | Group Comparison |
| 3+ | Independent | No | Kruskal-Wallis + post-hoc | Group Comparison |
| 3+ | Repeated | Yes | RM ANOVA | **Repeated Measures** |
| 3+ | Repeated | No | Friedman test | **Repeated Measures** |
| - | Correlation | Yes | Pearson r | Correlation |
| - | Correlation | No | Spearman rho | Correlation |
| - | Regression | - | Multiple linear regression | **Linear Regression** |

---

## Step 2B: Categorical Outcome

| Groups | Pairing | Condition | Test | analyze-stats type |
|--------|---------|-----------|------|--------------------|
| 1 | - | - | Binomial / z-test for proportion | Group Comparison |
| 2 | Independent | Expected >= 5 | Chi-squared | Group Comparison |
| 2 | Independent | Expected < 5 | Fisher's exact | Group Comparison |
| 2 | Paired | - | McNemar's test | Group Comparison |
| 2+ | Independent | Ordered | Chi-squared trend | Group Comparison |
| 3+ | Independent | - | Chi-squared | Group Comparison |
| 3+ | Paired | - | Cochran's Q | Group Comparison |
| - | Regression | Binary outcome | Logistic regression | **Logistic Regression** |

---

## Step 3: Confounder Control

| Situation | Method | analyze-stats type |
|-----------|--------|--------------------|
| Continuous outcome + multivariable | Multiple linear regression | **Linear Regression** |
| Binary outcome + multivariable | Logistic regression | **Logistic Regression** |
| Survival outcome + multivariable | Cox proportional hazards | Survival |
| Observational study + treatment comparison | Propensity score | **Propensity Score** |
| Repeated measures + missing data | LMM / GEE | **Repeated Measures** |

---

## Normality Assessment

| Method | When | Criterion |
|--------|------|-----------|
| Shapiro-Wilk | n < 50 | p >= 0.05 → normal |
| Kolmogorov-Smirnov | n >= 50 | p >= 0.05 → normal |
| Q-Q plot | Always (visual) | Points on diagonal |
| Skewness/Kurtosis | Supplementary | |skew| < 2, |kurt| < 7 |

Practical rule: n > 30 → t-test is generally robust (CLT), unless extreme skew or outliers.

---

## Common Reviewer Flags

- Using independent test on paired data
- Using chi-squared when expected cell count < 5 (should use Fisher's exact)
- Not reporting assumption check results
- Missing multiple comparison correction for 3+ groups
- Not specifying the test selection rationale in Methods

# Publication Table Standards — Knowledge Base

> Reference document for medical journal table formatting.
> Source: YouTube tutorials, journal author guidelines, AMA Manual of Style, tool documentation.
> Last updated: 2026-04-11

---

## 1. Universal Rules (All Medical Journals)

1. **No vertical lines** — horizontal rules only (top, below header, bottom)
2. **Editable Word tables** — use Word Insert > Table. Never submit images, Excel, or PDF
3. **Sequential numbering** — Table 1, 2, 3... in order of first citation in text
4. **Define all abbreviations** — in footnotes, independently for each table
5. **Explanatory content in footnotes, not headings**
6. **Display item limit** — most journals: tables + figures combined 4-5
7. **Exact P values** — never just "significant/not significant"
8. **Variability measures required** — always state mean (SD) or median (IQR)
9. **Self-contained titles** — table title alone must convey content without reading text
10. **No duplication** — tables and figures must not repeat the same data

---

## 2. Journal-Specific Differences

### Footnote Markers

| Journal Family | Marker System | Order |
|---|---|---|
| AMA (Radiology, JAMA, AJR, Rad:AI) | Superscript lowercase letters | a, b, c, d... |
| NEJM | Symbols | *, †, ‡, §, ‖, ¶, **, ††... |
| Lancet | Symbols | *, †, ‡, §, ¶, ‖, **, ††... |
| European Radiology (Springer) | Superscript lowercase letters + asterisks for significance | a, b, c... and *, **, *** |

### P Value Formatting

| Journal | Case | Leading Zero | Examples |
|---|---|---|---|
| Radiology / JAMA / AJR / Rad:AI | Uppercase italic *P* | No | *P* = .03, *P* < .001 |
| NEJM | Uppercase *P* | No | *P* = .04, *P* < .001 |
| Lancet | Lowercase p | Yes | p = 0.03, p < 0.0001 |
| European Radiology | Lowercase p (Springer) | Varies | p = 0.03 |

### P Value Decimal Places

| Value Range | JAMA/Radiology | NEJM | Lancet |
|---|---|---|---|
| > .01 | 2 decimals (.04) | 2 decimals (.04) | 2 sig figs (0.04) |
| .01-.001 | 3 decimals (.003) | 3 decimals (.003) | Up to 4 decimals (0.0023) |
| < .001 | *P* < .001 | *P* < .001 | p < 0.0001 |

### Horizontal Lines

| Journal | Rule |
|---|---|
| AMA journals | Top, below header, bottom. Additional sparingly |
| NEJM | Minimize even internal horizontal lines |
| Lancet | Exactly 3: top, below header, bottom |

### 95% CI Format

| Journal | Format |
|---|---|
| Radiology | (XX, XX) — comma separator |
| JAMA | XX to XX — "to" separator |
| Lancet | XX to XX or XX-XX |
| NEJM | XX to XX |

### Table Title Format

| Journal | Format |
|---|---|
| AMA journals | **Table N.** Title in regular weight (bold "Table N." only) |
| NEJM | Table N. Title |
| Lancet | Table N: Title |

---

## 3. AMA Manual of Style (11th Ed) — Table Rules

### Structure
- **Column headers**: Bold
- **Row headers (stub)**: Regular weight, sentence case
- **Decking** (nested column headers): Maximum 3 levels
- **Alignment**: Text = left-aligned, Numbers = center or decimal-aligned
- **Units**: In column header parentheses, not repeated in cells
- **Maximum size**: ~40 rows x 6-8 columns, fit within 1 page

### Footnote System (AMA 11th)
- **Table footnotes**: Superscript lowercase letters (a, b, c...)
- **NOT symbols** — symbols (*, †, ‡) are for bottom-of-page text footnotes only
- **Assignment order**: Top-to-bottom, left-to-right through the table
- **Whole-table footnote**: Place superscript "a" at end of table title
- **Reference numbers + footnotes**: Reference first, comma, then footnote letter (e.g., 5,b)

### Footnote Placement Order
```
Below the table, in this sequence:

1. General note (no marker, applies to entire table)
   "Data are presented as median (IQR) unless otherwise noted."

2. Abbreviation note
   "ASA = American Society of Anesthesiologists, BMI = body mass index,
    CI = confidence interval."
   → Listed in order of appearance in table (left→right, top→bottom)
   → Some journals (JAMA) require alphabetical order — check guide

3. Specific notes (superscript letter markers)
   "a Excludes patients with missing follow-up data."
   "b Adjusted for age and sex."

4. Probability notes (asterisk markers, if used)
   "* P < .05; ** P < .01."
```

### P Value Rules
- Uppercase italic: *P*
- No leading zero: .05, not 0.05
- Thin space around operators: *P* = .03
- Only *P*, *α*, *β* omit leading zero; all other statistics keep it

---

## 4. Footnote System Specification

### When to Use Symbols vs Letters vs Numbers

| System | Markers | Use Case | Journals |
|---|---|---|---|
| Lowercase letters | a, b, c... | **Table footnotes (standard)** | AMA journals, EUR, most |
| Symbols | *, †, ‡, §, ‖, ¶ | Table footnotes (traditional) | NEJM, Lancet |
| Asterisks only | *, **, *** | **Probability notes only** | Springer journals |
| Uppercase letters | A, B, C... | Table footnotes | JCI |
| Numbers | 1, 2, 3... | Non-numeric tables only | Rarely used (confusion with references) |

### Caption vs Footnote — What Goes Where

| Content | Caption (above table) | Footnote (below table) |
|---|---|---|
| Table title | Yes | No |
| Study period/setting | Yes (brief) | No |
| Abbreviation definitions | No | **Yes** (first footnote) |
| Statistical test descriptions | No | **Yes** |
| P value thresholds | No | **Yes** (probability note) |
| Data source | No | **Yes** (source note) |
| Missing data explanation | No | **Yes** (specific note) |
| Sample size (N) | Column header or caption | Footnote (supplementary) |

### Statistical Test Footnote Patterns

**Pattern A — Individual markers (preferred when tests vary by row):**
```
Age, y                    65.3 (12.1)    62.1 (11.8)    .04ᵃ
Sex, male                 53 (52)        48 (47)        .41ᵇ

ᵃ Wilcoxon rank-sum test.
ᵇ Fisher exact test.
```

**Pattern B — General note (preferred when same test for all rows of a type):**
```
P values were calculated using the Wilcoxon rank-sum test for continuous
variables and the Fisher exact test for categorical variables.
```

---

## 5. Common Mistakes Checklist

### Structure
- [ ] Binary variables: show only one level (Male 53%, not Male 53% / Female 47%)
- [ ] Remove derivable columns (Total column when groups are shown)
- [ ] Keep tables under 50 rows (excess → supplementary)
- [ ] No sub-part numbering (Table 3A/3B → Table 3, Table 4)
- [ ] Do not repeat all table data in Results text

### Formatting
- [ ] Consistent decimal places within each column
- [ ] Units in column headers, not in cells
- [ ] No vertical lines, minimal horizontal lines
- [ ] Numbers right-aligned or center-aligned (not left)
- [ ] Specify variability: "Mean (SD)" or "Median (IQR)" in header or footnote

### Statistics
- [ ] RCTs: P values in Table 1 are usually unnecessary (randomization)
- [ ] Always name the statistical test (in footnote or general note)
- [ ] Report effect sizes per clinically meaningful unit (OR per 10-year, not per 1-year)
- [ ] CI notation consistent throughout (parentheses vs "to" — match journal style)
- [ ] Never use "NS" — report exact P values

### Submission
- [ ] Editable Word table (not image/screenshot)
- [ ] No color as sole information carrier (use bold/italic/symbols instead)
- [ ] All abbreviations defined in footnotes
- [ ] Each table's footnotes are self-contained (no "see Table 1 footnote")

---

## 6. Tool Recommendations

### Primary Pipeline: R {gtsummary}

**Why gtsummary:**
- Built-in journal themes: JAMA, Lancet, NEJM
- Auto-selects statistical tests (Wilcoxon, Fisher, Chi-square)
- Auto-generates footnotes (test names, summary statistics)
- Outputs to Word (via flextable), LaTeX (via huxtable), HTML (via gt)

**Core API:**
```r
library(gtsummary)

# Set journal theme FIRST
theme_gtsummary_journal("jama")  # or "lancet", "nejm"
theme_gtsummary_compact()

# Table 1
tbl <- df %>%
  tbl_summary(
    by = group,
    type = list(age ~ "continuous2"),
    statistic = list(
      all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})"),
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = list(all_continuous() ~ 1),
    missing = "ifany"
  ) %>%
  add_p() %>%
  add_overall() %>%
  add_stat_label() %>%
  bold_labels() %>%
  modify_footnote(all_stat_cols() ~ "Mean (SD); Median (Q1, Q3); n (%)")

# Export
tbl %>% as_flex_table() %>% flextable::save_as_docx(path = "table1.docx")
tbl %>% as_hux_table() %>% huxtable::to_latex() %>% writeLines("table1.tex")
```

**Regression table:**
```r
model <- glm(outcome ~ age + sex + bmi, data = df, family = binomial)

tbl_regression(model, exponentiate = TRUE) %>%
  add_global_p() %>%
  bold_p(t = 0.05) %>%
  bold_labels()
```

### Supporting Tools

| Tool | Role | When to Use |
|---|---|---|
| gt | Fine-grained styling | Cell-level colors, custom footnote marks, heatmaps |
| flextable | Word output engine | Final DOCX formatting, autofit, borders |
| huxtable | LaTeX output engine | LaTeX code extraction |
| tableone (R/Python) | Quick exploratory Table 1 | Early drafts, data exploration |
| python-docx | Python Word tables | When pipeline is Python-only |
| great_tables | Python gt port | Growing but immature for medical use |

### Python Limitation
Python has no equivalent to gtsummary's journal themes or automatic statistical test selection. For publication tables, R is strongly recommended. Python-only pipelines should use tableone for generation + python-docx for Word formatting, accepting manual style adjustments.

---

## 7. Format-Specific Notes

### Word (DOCX) Submission
- Use Word's Insert > Table function (never tab-separated)
- Font: Times New Roman 10-12pt (table body can be 9-10pt)
- Single spacing within cells
- Remove all vertical borders; keep 3 horizontal lines
- Bold headers, regular body text
- Each table on a separate page, after References

### LaTeX Submission
- Use `booktabs` package (`\toprule`, `\midrule`, `\bottomrule`)
- Never use `\hline` or `|` vertical separators
- `\siunitx` for decimal alignment
- `\multirow` / `\multicolumn` for merged cells
- Generate via gtsummary → huxtable → to_latex() for consistency

### HTML (Review/Proofing)
- gt produces the highest quality HTML tables
- Useful for co-author review before final Word/LaTeX export
- Can include interactive features (sorting, filtering) for supplementary materials

# Publication Table Tools — Comparison

## R Packages

| Feature | gtsummary | gt | flextable | tableone | kableExtra | huxtable |
|---|---|---|---|---|---|---|
| **Primary Use** | Table 1 + regression | Fine styling | Word output | Quick Table 1 | LaTeX/PDF | Multi-format |
| **Auto Statistics** | Yes (test selection) | No | No | Yes (p, SMD) | No | huxreg() only |
| **Journal Themes** | JAMA, Lancet, NEJM | No | Custom | No | No | No |
| **Footnotes** | Via gt engine | Native (excellent) | Native | Manual | Native | Native |
| **Word Output** | Via flextable | Limited | **Best** | Via knitr | Limited | Yes |
| **LaTeX Output** | Via huxtable | Native | Limited | Via knitr | **Best** | Native |
| **HTML Output** | Via gt | **Best** | Yes | Via knitr | Yes | Yes |
| **Maintenance** | Active (MSKCC) | Active (Posit) | Active | Maintenance | Maintenance | Maintenance |

## Python Libraries

| Feature | great_tables | tableone (py) | python-docx | pandas Styler |
|---|---|---|---|---|
| **Primary Use** | gt port | Quick Table 1 | Word tables | HTML/LaTeX |
| **Auto Statistics** | No | Yes (p, SMD) | No | No |
| **Journal Themes** | No | No | No | No |
| **Footnotes** | tab_footnote() | No | Manual | No |
| **Word Output** | No | No | **Native** | No |
| **LaTeX Output** | No | Via export | No | to_latex() |
| **Maturity** | Growing | Stable | Stable | Built-in |

## Recommended Pipelines

### Best: R gtsummary → Word
```r
theme_gtsummary_journal("jama")
theme_gtsummary_compact()
tbl %>% as_flex_table() %>% flextable::save_as_docx(path = "table.docx")
```

### Best: R gtsummary → LaTeX
```r
tbl %>% as_hux_table() %>% huxtable::to_latex()
# OR
tbl %>% as_gt() %>% gt::gtsave("table.tex")
```

### Python-only (acceptable, not ideal)
```python
# Generation
from tableone import TableOne
table1 = TableOne(df, columns=cols, categorical=cats, groupby='group', pval=True)
# Export to CSV, then format with python-docx
```

## Key Insight
> **Python lacks gtsummary-equivalent journal themes and auto-statistics.**
> For publication tables in medical journals, R is strongly preferred.
> Python pipelines require manual style adjustments that R handles automatically.

## Tool Selection Decision Tree

```
Need Table 1 or regression table?
  → gtsummary (+ journal theme)

Need pixel-perfect custom styling?
  → gt (HTML) or flextable (Word)

Submitting to Word-only journal?
  → gtsummary → as_flex_table() → save_as_docx()

Submitting LaTeX?
  → gtsummary → as_hux_table() → to_latex()
  → OR kableExtra for direct LaTeX control

Python-only constraint?
  → tableone (generation) + python-docx (formatting)
  → Accept manual style work

Quick exploratory table?
  → tableone (R or Python) — fastest setup
```

"""
Template: Inter-rater Agreement Analysis
Calculates Cohen's/Fleiss' kappa, ICC, and Krippendorff's alpha with bootstrap CIs.
Generates Bland-Altman plots for continuous measurements.

Usage:
    Modify the CONFIGURATION section below, then run:
        python agreement_analysis.py

Input:  CSV with rows=items, columns=raters (or long format)
Output: agreement_table.csv, bland_altman.pdf/.png, summary text
"""

# === REPRODUCIBILITY HEADER ===
import sys
import os
import datetime
import numpy as np
import pandas as pd
from scipy import stats

np.random.seed(42)
print(f"Date: {datetime.date.today()}")
print(f"Python: {sys.version}")
print(f"numpy: {np.__version__}, pandas: {pd.__version__}, scipy: {stats.scipy.__version__}")

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

STYLE_PATH = os.path.expanduser(
    "~/.claude/skills/analyze-stats/references/style/figure_style.mplstyle"
)
if os.path.exists(STYLE_PATH):
    plt.style.use(STYLE_PATH)

print()

# === CONFIGURATION (modify for your study) ===
INPUT_FILE = "data.csv"           # Path to input data
OUTPUT_DIR = "."                  # Output directory
RATER_COLS = []                   # Column names for each rater (wide format)
DATA_TYPE = "auto"                # "categorical", "ordinal", "continuous", or "auto"
BOOTSTRAP_N = 1000                # Number of bootstrap iterations for CIs
ALPHA = 0.05                      # Significance level
# ==============================================


def load_data(filepath: str) -> pd.DataFrame:
    """Load CSV or Excel file."""
    if filepath.endswith((".xlsx", ".xls")):
        return pd.read_excel(filepath)
    return pd.read_csv(filepath)


def detect_data_type(df: pd.DataFrame, cols: list) -> str:
    """Auto-detect whether ratings are categorical, ordinal, or continuous."""
    combined = pd.concat([df[c] for c in cols], ignore_index=True).dropna()
    n_unique = combined.nunique()

    if combined.dtype == object or combined.dtype.name == "category":
        return "categorical"
    if n_unique <= 10 and all(combined == combined.astype(int)):
        return "ordinal"
    return "continuous"


def interpret_kappa(kappa: float) -> str:
    """Interpret kappa using Landis & Koch (1977) guidelines."""
    if kappa < 0:
        return "poor"
    if kappa < 0.20:
        return "slight"
    if kappa < 0.40:
        return "fair"
    if kappa < 0.60:
        return "moderate"
    if kappa < 0.80:
        return "substantial"
    return "almost perfect"


def interpret_icc(icc: float) -> str:
    """Interpret ICC using Cicchetti (1994) guidelines."""
    if icc < 0.40:
        return "poor"
    if icc < 0.60:
        return "fair"
    if icc < 0.75:
        return "good"
    return "excellent"


def cohens_kappa(r1: np.ndarray, r2: np.ndarray) -> float:
    """Cohen's kappa for 2 raters, categorical data."""
    from sklearn.metrics import cohen_kappa_score
    return cohen_kappa_score(r1, r2)


def cohens_weighted_kappa(r1: np.ndarray, r2: np.ndarray,
                          weights: str = "linear") -> float:
    """Weighted Cohen's kappa for ordinal data."""
    from sklearn.metrics import cohen_kappa_score
    return cohen_kappa_score(r1, r2, weights=weights)


def fleiss_kappa(rating_matrix: np.ndarray) -> float:
    """Fleiss' kappa for multiple raters, categorical data.

    Args:
        rating_matrix: n_items x n_categories matrix of category counts per item
    """
    n_items, n_categories = rating_matrix.shape
    n_raters = int(rating_matrix.sum(axis=1)[0])

    p_j = rating_matrix.sum(axis=0) / (n_items * n_raters)
    P_i = (rating_matrix**2).sum(axis=1) - n_raters
    P_i = P_i / (n_raters * (n_raters - 1))

    P_bar = P_i.mean()
    P_e = (p_j**2).sum()

    if P_e == 1:
        return 1.0
    return (P_bar - P_e) / (1 - P_e)


def compute_icc(df: pd.DataFrame, cols: list,
                model: str = "two-way", type_: str = "agreement") -> float:
    """Compute Intraclass Correlation Coefficient.

    Args:
        model: "one-way", "two-way" (random), or "two-way-mixed"
        type_: "consistency" or "agreement" (only for two-way)
    Returns: ICC value
    """
    ratings = df[cols].values
    n_subjects, n_raters = ratings.shape

    grand_mean = ratings.mean()

    ss_total = np.sum((ratings - grand_mean) ** 2)
    row_means = ratings.mean(axis=1)
    col_means = ratings.mean(axis=0)

    ss_rows = n_raters * np.sum((row_means - grand_mean) ** 2)
    ss_cols = n_subjects * np.sum((col_means - grand_mean) ** 2)
    ss_error = ss_total - ss_rows - ss_cols

    ms_rows = ss_rows / (n_subjects - 1)
    ms_cols = ss_cols / (n_raters - 1) if n_raters > 1 else 0
    ms_error = ss_error / ((n_subjects - 1) * (n_raters - 1)) if (n_subjects > 1 and n_raters > 1) else 0

    if model == "one-way":
        ms_within = (ss_total - ss_rows) / (n_subjects * (n_raters - 1))
        icc = (ms_rows - ms_within) / (ms_rows + (n_raters - 1) * ms_within)
    elif model in ("two-way", "two-way-mixed"):
        if type_ == "agreement":
            icc = (ms_rows - ms_error) / (
                ms_rows + (n_raters - 1) * ms_error +
                n_raters * (ms_cols - ms_error) / n_subjects
            )
        else:  # consistency
            icc = (ms_rows - ms_error) / (ms_rows + (n_raters - 1) * ms_error)
    else:
        raise ValueError(f"Unknown model: {model}")

    return icc


def bootstrap_ci(data: pd.DataFrame, cols: list, metric_func,
                 n_bootstrap: int = 1000, alpha: float = 0.05,
                 **kwargs) -> tuple:
    """Bootstrap confidence interval for any agreement metric."""
    n = len(data)
    boot_values = []

    for _ in range(n_bootstrap):
        idx = np.random.choice(n, size=n, replace=True)
        boot_data = data.iloc[idx]
        try:
            val = metric_func(boot_data, cols, **kwargs)
            if np.isfinite(val):
                boot_values.append(val)
        except Exception:
            continue

    if len(boot_values) < 10:
        return (np.nan, np.nan)

    boot_values = np.array(boot_values)
    ci_lo = np.percentile(boot_values, 100 * alpha / 2)
    ci_hi = np.percentile(boot_values, 100 * (1 - alpha / 2))
    return (ci_lo, ci_hi)


def bland_altman_plot(r1: np.ndarray, r2: np.ndarray,
                      name1: str, name2: str, output_dir: str) -> dict:
    """Generate Bland-Altman plot for two continuous raters."""
    mean_vals = (r1 + r2) / 2
    diff_vals = r1 - r2
    mean_diff = diff_vals.mean()
    std_diff = diff_vals.std(ddof=1)

    loa_upper = mean_diff + 1.96 * std_diff
    loa_lower = mean_diff - 1.96 * std_diff

    fig, ax = plt.subplots(figsize=(3.5, 3.5))
    ax.scatter(mean_vals, diff_vals, s=15, alpha=0.6, color="#0072B2", edgecolors="none")
    ax.axhline(mean_diff, color="#D55E00", linewidth=1, label=f"Mean: {mean_diff:.2f}")
    ax.axhline(loa_upper, color="#D55E00", linewidth=0.8, linestyle="--",
               label=f"+1.96 SD: {loa_upper:.2f}")
    ax.axhline(loa_lower, color="#D55E00", linewidth=0.8, linestyle="--",
               label=f"-1.96 SD: {loa_lower:.2f}")
    ax.set_xlabel(f"Mean of {name1} and {name2}")
    ax.set_ylabel(f"Difference ({name1} - {name2})")
    ax.legend(fontsize=7, loc="upper right")

    fig.tight_layout()
    pdf_path = os.path.join(output_dir, "bland_altman.pdf")
    png_path = os.path.join(output_dir, "bland_altman.png")
    fig.savefig(pdf_path, format="pdf", bbox_inches="tight")
    fig.savefig(png_path, format="png", dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: {pdf_path}")
    print(f"Saved: {png_path}")

    return {
        "mean_diff": mean_diff,
        "sd_diff": std_diff,
        "loa_lower": loa_lower,
        "loa_upper": loa_upper,
    }


def _icc_for_bootstrap(data, cols, model="two-way", type_="agreement"):
    """Wrapper for ICC to use with bootstrap_ci."""
    return compute_icc(data, cols, model=model, type_=type_)


def _kappa_for_bootstrap_2raters(data, cols):
    """Wrapper for Cohen's kappa to use with bootstrap_ci."""
    r1 = data[cols[0]].values
    r2 = data[cols[1]].values
    mask = ~(pd.isna(r1) | pd.isna(r2))
    if mask.sum() < 2:
        return np.nan
    return cohens_kappa(r1[mask], r2[mask])


def analyze_categorical(df: pd.DataFrame, cols: list, output_dir: str) -> list:
    """Run agreement analysis for categorical data."""
    results = []
    n_raters = len(cols)
    n_items = len(df)

    if n_raters == 2:
        r1 = df[cols[0]].values
        r2 = df[cols[1]].values
        mask = ~(pd.isna(r1) | pd.isna(r2))
        r1_clean, r2_clean = r1[mask], r2[mask]

        kappa = cohens_kappa(r1_clean, r2_clean)
        ci_lo, ci_hi = bootstrap_ci(df, cols, _kappa_for_bootstrap_2raters,
                                     n_bootstrap=BOOTSTRAP_N, alpha=ALPHA)

        results.append({
            "Metric": "Cohen's kappa",
            "Value": f"{kappa:.3f}",
            "95% CI": f"({ci_lo:.3f}-{ci_hi:.3f})",
            "Interpretation": interpret_kappa(kappa),
            "n_items": int(mask.sum()),
            "n_raters": 2,
        })

        pct_agree = np.mean(r1_clean == r2_clean)
        results.append({
            "Metric": "Percent agreement",
            "Value": f"{pct_agree:.3f}",
            "95% CI": "",
            "Interpretation": "",
            "n_items": int(mask.sum()),
            "n_raters": 2,
        })
    else:
        categories = sorted(set(pd.concat([df[c] for c in cols]).dropna().unique()))
        rating_counts = np.zeros((n_items, len(categories)))
        for i in range(n_items):
            for c in cols:
                val = df[c].iloc[i]
                if pd.notna(val) and val in categories:
                    rating_counts[i, categories.index(val)] += 1

        fk = fleiss_kappa(rating_counts)
        results.append({
            "Metric": "Fleiss' kappa",
            "Value": f"{fk:.3f}",
            "95% CI": "",
            "Interpretation": interpret_kappa(fk),
            "n_items": n_items,
            "n_raters": n_raters,
        })

    return results


def analyze_continuous(df: pd.DataFrame, cols: list, output_dir: str) -> list:
    """Run agreement analysis for continuous data."""
    results = []
    n_raters = len(cols)

    for type_ in ["agreement", "consistency"]:
        icc = compute_icc(df, cols, model="two-way", type_=type_)
        ci_lo, ci_hi = bootstrap_ci(df, cols, _icc_for_bootstrap,
                                     n_bootstrap=BOOTSTRAP_N, alpha=ALPHA,
                                     model="two-way", type_=type_)
        label = f"ICC (two-way, {type_})"
        results.append({
            "Metric": label,
            "Value": f"{icc:.3f}",
            "95% CI": f"({ci_lo:.3f}-{ci_hi:.3f})",
            "Interpretation": interpret_icc(icc),
            "n_items": len(df),
            "n_raters": n_raters,
        })

    if n_raters == 2:
        r1 = df[cols[0]].values.astype(float)
        r2 = df[cols[1]].values.astype(float)
        mask = ~(np.isnan(r1) | np.isnan(r2))
        ba = bland_altman_plot(r1[mask], r2[mask], cols[0], cols[1], output_dir)
        results.append({
            "Metric": "Bland-Altman mean diff",
            "Value": f"{ba['mean_diff']:.3f}",
            "95% CI": f"LoA: ({ba['loa_lower']:.3f}-{ba['loa_upper']:.3f})",
            "Interpretation": "",
            "n_items": int(mask.sum()),
            "n_raters": 2,
        })

    return results


def analyze_ordinal(df: pd.DataFrame, cols: list, output_dir: str) -> list:
    """Run agreement analysis for ordinal data (weighted kappa + ICC)."""
    results = []

    if len(cols) == 2:
        r1 = df[cols[0]].values
        r2 = df[cols[1]].values
        mask = ~(pd.isna(r1) | pd.isna(r2))
        r1_clean, r2_clean = r1[mask], r2[mask]

        for weight in ["linear", "quadratic"]:
            wk = cohens_weighted_kappa(r1_clean, r2_clean, weights=weight)
            results.append({
                "Metric": f"Weighted kappa ({weight})",
                "Value": f"{wk:.3f}",
                "95% CI": "",
                "Interpretation": interpret_kappa(wk),
                "n_items": int(mask.sum()),
                "n_raters": 2,
            })

    results.extend(analyze_continuous(df, cols, output_dir))
    return results


def save_results(results: list, output_dir: str) -> None:
    """Save agreement metrics as CSV and print formatted output."""
    df = pd.DataFrame(results)
    csv_path = os.path.join(output_dir, "agreement_table.csv")
    df.to_csv(csv_path, index=False)
    print(f"\nSaved: {csv_path}")
    print("\n--- Agreement Analysis Results ---\n")
    print(df.to_markdown(index=False))


def print_results_text(results: list, data_type: str) -> None:
    """Print manuscript-ready text for Results section."""
    print("\n--- Results Text (copy-paste ready) ---\n")

    for r in results:
        metric = r["Metric"]
        value = r["Value"]
        ci = r.get("95% CI", "")
        interp = r.get("Interpretation", "")
        n_items = r.get("n_items", "")
        n_raters = r.get("n_raters", "")

        text = f"{metric} was {value}"
        if ci:
            text += f" (95% CI: {ci})"
        if interp:
            text += f", indicating {interp} agreement"
        text += f" (n = {n_items} items, {n_raters} raters)."
        print(text)
    print()


# === MAIN ===
if __name__ == "__main__":
    print("=" * 60)
    print("Inter-rater Agreement Analysis")
    print("=" * 60)

    df = load_data(INPUT_FILE)
    print(f"\nLoaded: {INPUT_FILE} ({df.shape[0]} rows, {df.shape[1]} columns)")

    if not RATER_COLS:
        RATER_COLS = [c for c in df.columns if c.lower().startswith("rater")]
        if not RATER_COLS:
            RATER_COLS = list(df.columns)
        print(f"Using rater columns: {RATER_COLS}")

    if DATA_TYPE == "auto":
        DATA_TYPE = detect_data_type(df, RATER_COLS)
    print(f"Data type: {DATA_TYPE}")
    print(f"Items: {len(df)}, Raters: {len(RATER_COLS)}")

    for col in RATER_COLS:
        n_miss = df[col].isna().sum()
        if n_miss > 0:
            print(f"  Missing in {col}: {n_miss} ({100*n_miss/len(df):.1f}%)")

    if DATA_TYPE == "categorical":
        results = analyze_categorical(df, RATER_COLS, OUTPUT_DIR)
    elif DATA_TYPE == "ordinal":
        results = analyze_ordinal(df, RATER_COLS, OUTPUT_DIR)
    elif DATA_TYPE == "continuous":
        results = analyze_continuous(df, RATER_COLS, OUTPUT_DIR)
    else:
        raise ValueError(f"Unknown data type: {DATA_TYPE}")

    save_results(results, OUTPUT_DIR)
    print_results_text(results, DATA_TYPE)

"""
Template: Diagnostic Accuracy Analysis
Calculates sensitivity, specificity, PPV, NPV, accuracy, AUC with 95% CIs.
Generates ROC curve and optional model comparison.

Usage:
    Modify the CONFIGURATION section below, then run:
        python diagnostic_accuracy.py

Input:  CSV with ground truth and predicted scores/labels
Output: diagnostic_accuracy_table.csv, roc_curve.pdf/.png, summary text
"""

# === REPRODUCIBILITY HEADER ===
import sys
import os
import datetime
import numpy as np
import pandas as pd
import scipy
from scipy import stats

import sklearn

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

STYLE_PATH = os.path.join(os.path.dirname(os.path.dirname(__file__)), "style", "figure_style.mplstyle")
if os.path.exists(STYLE_PATH):
    plt.style.use(STYLE_PATH)

# === CONFIGURATION (modify for your study) ===
INPUT_FILE = "data.csv"           # Path to input data
OUTPUT_DIR = "."                  # Output directory
TRUTH_COL = "ground_truth"       # Column: binary ground truth (0/1)
SCORE_COLS = ["model_score"]     # Column(s): predicted probability/score (for ROC)
PRED_COLS = ["model_pred"]       # Column(s): binary predictions (0/1) at chosen threshold
MODEL_NAMES = ["Model"]          # Display names for each model
THRESHOLD = None                  # Fixed threshold (None = use Youden's optimal)
COMPARE_MODELS = False            # True to run DeLong test between models
POSITIVE_LABEL = 1                # Value representing positive class
# ==============================================


def wilson_ci(p: float, n: int, alpha: float = 0.05) -> tuple:
    """Wilson score confidence interval for a proportion."""
    if n == 0:
        return (np.nan, np.nan)  # No denominator: undefined, not zero performance.
    z = stats.norm.ppf(1 - alpha / 2)
    denominator = 1 + z**2 / n
    center = (p + z**2 / (2 * n)) / denominator
    spread = z * np.sqrt((p * (1 - p) + z**2 / (4 * n)) / n) / denominator
    return (max(0.0, center - spread), min(1.0, center + spread))


def delong_auc_variance(y_true: np.ndarray, y_score: np.ndarray) -> float:
    """Estimate AUC variance using the DeLong method."""
    pos = y_score[y_true == 1]
    neg = y_score[y_true == 0]
    m = len(pos)
    n = len(neg)

    v_pos = np.array([np.mean(neg < p) + 0.5 * np.mean(neg == p) for p in pos])
    v_neg = np.array([np.mean(pos > nv) + 0.5 * np.mean(pos == nv) for nv in neg])

    var_auc = (np.var(v_pos, ddof=1) / m) + (np.var(v_neg, ddof=1) / n)
    return var_auc


def delong_ci(y_true: np.ndarray, y_score: np.ndarray,
              alpha: float = 0.05) -> tuple:
    """AUC with DeLong 95% CI."""
    from sklearn.metrics import roc_auc_score

    auc = roc_auc_score(y_true, y_score)
    var = delong_auc_variance(y_true, y_score)
    se = np.sqrt(var)
    z = stats.norm.ppf(1 - alpha / 2)
    ci_low = max(0.0, auc - z * se)
    ci_high = min(1.0, auc + z * se)
    return auc, ci_low, ci_high


def delong_test(y_true: np.ndarray, y_score1: np.ndarray,
                y_score2: np.ndarray) -> tuple:
    """DeLong test for comparing two AUCs on the same dataset."""
    from sklearn.metrics import roc_auc_score

    auc1 = roc_auc_score(y_true, y_score1)
    auc2 = roc_auc_score(y_true, y_score2)

    var1 = delong_auc_variance(y_true, y_score1)
    var2 = delong_auc_variance(y_true, y_score2)

    pos_mask = y_true == 1
    neg_mask = y_true == 0
    m = pos_mask.sum()
    n = neg_mask.sum()

    v1_pos = np.array([np.mean(y_score1[neg_mask] < p) +
                        0.5 * np.mean(y_score1[neg_mask] == p) for p in y_score1[pos_mask]])
    v2_pos = np.array([np.mean(y_score2[neg_mask] < p) +
                        0.5 * np.mean(y_score2[neg_mask] == p) for p in y_score2[pos_mask]])
    v1_neg = np.array([np.mean(y_score1[pos_mask] > nv) +
                        0.5 * np.mean(y_score1[pos_mask] == nv) for nv in y_score1[neg_mask]])
    v2_neg = np.array([np.mean(y_score2[pos_mask] > nv) +
                        0.5 * np.mean(y_score2[pos_mask] == nv) for nv in y_score2[neg_mask]])

    cov = np.cov(v1_pos, v2_pos)[0, 1] / m + np.cov(v1_neg, v2_neg)[0, 1] / n

    z = (auc1 - auc2) / np.sqrt(var1 + var2 - 2 * cov)
    p = 2 * stats.norm.sf(abs(z))
    return z, p


def youdens_threshold(y_true: np.ndarray, y_score: np.ndarray) -> float:
    """Find optimal threshold using Youden's J statistic."""
    from sklearn.metrics import roc_curve

    fpr, tpr, thresholds = roc_curve(y_true, y_score)
    j = tpr - fpr
    optimal_idx = np.argmax(j)
    return thresholds[optimal_idx]


def compute_metrics(y_true: np.ndarray, y_pred: np.ndarray,
                    y_score: np.ndarray = None) -> dict:
    """Compute diagnostic accuracy metrics with Wilson CIs."""
    tp = np.sum((y_pred == 1) & (y_true == 1))
    fp = np.sum((y_pred == 1) & (y_true == 0))
    tn = np.sum((y_pred == 0) & (y_true == 0))
    fn = np.sum((y_pred == 0) & (y_true == 1))
    n = len(y_true)

    sens = tp / (tp + fn) if (tp + fn) > 0 else np.nan
    spec = tn / (tn + fp) if (tn + fp) > 0 else np.nan
    ppv = tp / (tp + fp) if (tp + fp) > 0 else np.nan
    npv = tn / (tn + fn) if (tn + fn) > 0 else np.nan
    acc = (tp + tn) / n if n > 0 else np.nan

    metrics = {
        "Sensitivity": (sens, *wilson_ci(sens, tp + fn)),
        "Specificity": (spec, *wilson_ci(spec, tn + fp)),
        "PPV": (ppv, *wilson_ci(ppv, tp + fp)),
        "NPV": (npv, *wilson_ci(npv, tn + fn)),
        "Accuracy": (acc, *wilson_ci(acc, n)),
    }

    if y_score is not None:
        auc, auc_lo, auc_hi = delong_ci(y_true, y_score)
        metrics["AUC"] = (auc, auc_lo, auc_hi)

    metrics["_counts"] = {"TP": int(tp), "FP": int(fp), "TN": int(tn), "FN": int(fn)}
    return metrics


def plot_roc(y_true: np.ndarray, score_dict: dict, output_dir: str) -> None:
    """Generate ROC curve figure with AUC in legend."""
    from sklearn.metrics import roc_curve

    fig, ax = plt.subplots(figsize=(3.5, 3.5))
    colors = ["#0072B2", "#D55E00", "#009E73", "#CC79A7", "#F0E442"]

    for i, (name, y_score) in enumerate(score_dict.items()):
        fpr, tpr, _ = roc_curve(y_true, y_score)
        auc, ci_lo, ci_hi = delong_ci(y_true, y_score)
        label = f"{name}: AUC = {auc:.3f} ({ci_lo:.3f}-{ci_hi:.3f})"
        ax.plot(fpr, tpr, color=colors[i % len(colors)], linewidth=1.5, label=label)

    ax.plot([0, 1], [0, 1], color="gray", linestyle="--", linewidth=0.8)
    ax.set_xlabel("1 - Specificity (FPR)")
    ax.set_ylabel("Sensitivity (TPR)")
    ax.set_xlim([-0.02, 1.02])
    ax.set_ylim([-0.02, 1.02])
    ax.set_aspect("equal")
    ax.legend(loc="lower right", fontsize=7)

    fig.tight_layout()
    pdf_path = os.path.join(output_dir, "roc_curve.pdf")
    png_path = os.path.join(output_dir, "roc_curve.png")
    fig.savefig(pdf_path, format="pdf", bbox_inches="tight")
    fig.savefig(png_path, format="png", dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: {pdf_path}")
    print(f"Saved: {png_path}")


def plot_confusion_matrix(y_true: np.ndarray, pred_dict: dict,
                          model_names: list, output_dir: str) -> None:
    """Generate side-by-side confusion matrices using matplotlib."""
    n_models = len(pred_dict)
    fig, axes = plt.subplots(1, n_models, figsize=(3.5 * n_models, 3.5))
    if n_models == 1:
        axes = [axes]

    for ax, (name, y_pred) in zip(axes, pred_dict.items()):
        from sklearn.metrics import confusion_matrix as cm_func
        cm = cm_func(y_true, y_pred, labels=[0, 1])
        cm_pct = cm.astype(float) / cm.sum() * 100

        im = ax.imshow(cm, interpolation="nearest", cmap=plt.cm.Blues)
        ax.set_title(name, fontsize=10)
        ax.set_xlabel("Predicted")
        ax.set_ylabel("Actual")
        ax.set_xticks([0, 1])
        ax.set_yticks([0, 1])
        ax.set_xticklabels(["Neg", "Pos"])
        ax.set_yticklabels(["Neg", "Pos"])

        # Annotate cells with count and percentage
        thresh = cm.max() / 2.0
        for i in range(2):
            for j in range(2):
                ax.text(j, i, f"{cm[i, j]}\n({cm_pct[i, j]:.1f}%)",
                        ha="center", va="center", fontsize=9,
                        color="white" if cm[i, j] > thresh else "black")

    fig.tight_layout()
    pdf_path = os.path.join(output_dir, "confusion_matrix.pdf")
    png_path = os.path.join(output_dir, "confusion_matrix.png")
    fig.savefig(pdf_path, format="pdf", bbox_inches="tight",
                metadata={"CreationDate": None, "ModDate": None})
    fig.savefig(png_path, format="png", dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: {pdf_path}")
    print(f"Saved: {png_path}")


def plot_calibration(y_true: np.ndarray, score_dict: dict,
                     output_dir: str) -> None:
    """Generate calibration curves with Brier scores."""
    from sklearn.calibration import calibration_curve
    from sklearn.metrics import brier_score_loss

    fig, ax = plt.subplots(figsize=(3.5, 3.5))
    colors = ["#0072B2", "#D55E00", "#009E73", "#CC79A7", "#F0E442"]

    ax.plot([0, 1], [0, 1], color="gray", linestyle="--", linewidth=0.8,
            label="Perfect calibration")

    for i, (name, y_score) in enumerate(score_dict.items()):
        brier = brier_score_loss(y_true, y_score)
        fraction_pos, mean_predicted = calibration_curve(
            y_true, y_score, n_bins=10, strategy="uniform"
        )
        ax.plot(mean_predicted, fraction_pos, marker="o", markersize=4,
                color=colors[i % len(colors)], linewidth=1.5,
                label=f"{name} (Brier = {brier:.3f})")

    ax.set_xlabel("Mean predicted probability")
    ax.set_ylabel("Fraction of positives")
    ax.set_xlim([-0.02, 1.02])
    ax.set_ylim([-0.02, 1.02])
    ax.legend(loc="lower right", fontsize=7)

    fig.tight_layout()
    pdf_path = os.path.join(output_dir, "calibration_plot.pdf")
    png_path = os.path.join(output_dir, "calibration_plot.png")
    fig.savefig(pdf_path, format="pdf", bbox_inches="tight")
    fig.savefig(png_path, format="png", dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: {pdf_path}")
    print(f"Saved: {png_path}")


def save_performance_table(results: dict, output_dir: str) -> None:
    """Save performance metrics as CSV and print markdown."""
    rows = []
    for model_name, metrics in results.items():
        row = {"Model": model_name}
        for metric_name, vals in metrics.items():
            if metric_name.startswith("_"):
                continue
            val, ci_lo, ci_hi = vals
            row[metric_name] = f"{val:.3f} ({ci_lo:.3f}-{ci_hi:.3f})"
        counts = metrics.get("_counts", {})
        for k, v in counts.items():
            row[k] = v
        rows.append(row)

    df = pd.DataFrame(rows)
    csv_path = os.path.join(output_dir, "diagnostic_accuracy_table.csv")
    df.to_csv(csv_path, index=False)
    print(f"\nSaved: {csv_path}")
    print("\n--- Diagnostic Accuracy ---\n")
    print(df.to_markdown(index=False))


def print_results_text(results: dict) -> None:
    """Print manuscript-ready results text."""
    print("\n--- Results Text (copy-paste ready) ---\n")
    for model_name, metrics in results.items():
        parts = []
        for metric_name in ["AUC", "Sensitivity", "Specificity", "PPV", "NPV", "Accuracy"]:
            if metric_name in metrics:
                val, ci_lo, ci_hi = metrics[metric_name]
                parts.append(f"{metric_name} of {val:.3f} (95% CI: {ci_lo:.3f}-{ci_hi:.3f})")

        counts = metrics.get("_counts", {})
        n = sum(counts.values())
        n_pos = counts.get("TP", 0) + counts.get("FN", 0)
        n_neg = counts.get("TN", 0) + counts.get("FP", 0)

        print(f"{model_name} was evaluated on {n} cases "
              f"({n_pos} positive, {n_neg} negative). "
              f"The model achieved {', '.join(parts[:-1])}, and {parts[-1]}.")
        print()


# === MAIN ===
if __name__ == "__main__":
    np.random.seed(42)
    print(f"Date: {datetime.date.today()}")
    print(f"Python: {sys.version}")
    print(f"numpy: {np.__version__}, pandas: {pd.__version__}, scipy: {scipy.__version__}")
    print(f"sklearn: {sklearn.__version__}")
    print("=" * 60)
    print("Diagnostic Accuracy Analysis")
    print("=" * 60)

    df = pd.read_csv(INPUT_FILE)
    print(f"\nLoaded: {INPUT_FILE} ({df.shape[0]} rows, {df.shape[1]} columns)")

    y_true = df[TRUTH_COL].values

    # Prevalence
    prev = y_true.mean()
    print(f"Prevalence: {int(y_true.sum())}/{len(y_true)} ({100*prev:.1f}%)")

    all_results = {}
    score_dict = {}
    pred_dict = {}

    for i, (score_col, pred_col, name) in enumerate(
        zip(SCORE_COLS, PRED_COLS, MODEL_NAMES)
    ):
        print(f"\n--- {name} ---")
        y_score = df[score_col].values if score_col in df.columns else None
        if y_score is not None:
            score_dict[name] = y_score

        # Determine threshold
        if THRESHOLD is not None:
            thresh = THRESHOLD
        elif y_score is not None:
            thresh = youdens_threshold(y_true, y_score)
            print(f"Youden's optimal threshold: {thresh:.4f}")
            print(f"  WARNING: Youden's threshold optimized on evaluation data.")
            print(f"  For publication, use cross-validated thresholds or pre-specified cutoffs.")
        else:
            thresh = 0.5

        # Get predictions
        if pred_col in df.columns:
            y_pred = df[pred_col].values
        elif y_score is not None:
            y_pred = (y_score >= thresh).astype(int)
        else:
            raise ValueError(f"Neither prediction column '{pred_col}' nor "
                             f"score column '{score_col}' found.")

        pred_dict[name] = y_pred
        metrics = compute_metrics(y_true, y_pred, y_score)
        all_results[name] = metrics

    # ROC curve
    if score_dict:
        plot_roc(y_true, score_dict, OUTPUT_DIR)

    # Confusion matrix
    if pred_dict:
        plot_confusion_matrix(y_true, pred_dict, MODEL_NAMES, OUTPUT_DIR)

    # Calibration plot
    if score_dict:
        plot_calibration(y_true, score_dict, OUTPUT_DIR)

    # Model comparison (DeLong test)
    if COMPARE_MODELS and len(SCORE_COLS) >= 2:
        print("\n--- Model Comparison (DeLong Test) ---\n")
        for i in range(len(SCORE_COLS)):
            for j in range(i + 1, len(SCORE_COLS)):
                s1 = df[SCORE_COLS[i]].values
                s2 = df[SCORE_COLS[j]].values
                z, p = delong_test(y_true, s1, s2)
                from sklearn.metrics import roc_auc_score
                auc1 = roc_auc_score(y_true, s1)
                auc2 = roc_auc_score(y_true, s2)
                print(f"{MODEL_NAMES[i]} (AUC={auc1:.3f}) vs "
                      f"{MODEL_NAMES[j]} (AUC={auc2:.3f}): "
                      f"z = {z:.3f}, p = {p:.3f}")

    # Save outputs
    save_performance_table(all_results, OUTPUT_DIR)
    print_results_text(all_results)

"""
forest_plot.py — Publication-ready Forest Plot for Meta-analysis
================================================================
Generates a forest plot from a CSV file of study-level effect sizes.

Input CSV columns (required):
    study_label     : str  — Study author + year (e.g., "Kim 2022")
    effect_size     : float — Effect estimate (OR, RR, HR, MD, SMD, AUC)
    ci_lower        : float — Lower bound of 95% CI
    ci_upper        : float — Upper bound of 95% CI

Input CSV columns (optional):
    n_total         : int   — Total sample size (shown in table)
    weight          : float — Study weight % (determines box size)
    subgroup        : str   — Subgroup label (adds subgroup header rows)
    events_treat    : int   — Events in treatment/index group
    events_control  : int   — Events in control/comparator group

Pooled estimate (required as separate dict or CSV row with label="Pooled"):
    pooled_effect   : float
    pooled_ci_lower : float
    pooled_ci_upper : float
    i_squared       : float — I² heterogeneity (%)
    tau_squared     : float — τ² between-study variance
    q_p_value       : float — Cochran Q test p-value

Usage:
    python forest_plot.py --input studies.csv --pooled 0.82 0.65 1.04 \
        --i2 45.3 --tau2 0.021 --q-p 0.08 \
        --effect-label "OR" --null-value 1.0 \
        --output forest_plot --favor-left "Favors Treatment" --favor-right "Favors Control"

Outputs:
    forest_plot.pdf  — Vector format (submission quality)
    forest_plot.png  — Raster 300 DPI (web/preview)
"""

import argparse
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.lines import Line2D
from datetime import datetime

# ── Reproducibility ───────────────────────────────────────────────────────────
SCRIPT_VERSION = "1.0.0"
SEED = 42
np.random.seed(SEED)

# ── Wong (2011) colorblind-safe colors ────────────────────────────────────────
COLORS = {
    "study_ci": "#333333",      # dark gray for CI lines
    "study_box": "#0072B2",     # blue for individual study boxes
    "pooled_diamond": "#D55E00",# vermillion for pooled diamond
    "subgroup_header": "#E8E8E8",
    "grid": "#DDDDDD",
    "null_line": "#888888",
}

# ── Figure style ──────────────────────────────────────────────────────────────
plt.rcParams.update({
    "font.family": "Arial",
    "font.size": 9,
    "axes.labelsize": 9,
    "xtick.labelsize": 8,
    "ytick.labelsize": 8,
    "figure.dpi": 150,
})


def load_data(filepath: str) -> pd.DataFrame:
    """Load study data from CSV."""
    df = pd.read_csv(filepath)
    required = ["study_label", "effect_size", "ci_lower", "ci_upper"]
    missing = [c for c in required if c not in df.columns]
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    return df


def compute_box_size(weights: pd.Series, min_size=0.08, max_size=0.35) -> np.ndarray:
    """Scale box heights proportional to study weight."""
    if weights is None or weights.isna().all():
        return np.full(len(weights), (min_size + max_size) / 2)
    w = weights.fillna(weights.mean())
    w_norm = (w - w.min()) / (w.max() - w.min() + 1e-9)
    return min_size + w_norm * (max_size - min_size)


def make_forest_plot(
    df: pd.DataFrame,
    pooled_effect: float,
    pooled_ci_lower: float,
    pooled_ci_upper: float,
    i_squared: float,
    tau_squared: float,
    q_p_value: float,
    effect_label: str = "OR",
    null_value: float = 1.0,
    log_scale: bool = True,
    favor_left: str = "Favors Treatment",
    favor_right: str = "Favors Control",
    output_path: str = "forest_plot",
) -> None:
    """Generate and save the forest plot."""

    n_studies = len(df)
    has_subgroups = "subgroup" in df.columns and df["subgroup"].notna().any()
    has_weights = "weight" in df.columns

    # ── Layout calculations ───────────────────────────────────────────────────
    row_height = 0.45  # inches per row
    header_height = 1.2
    footer_height = 0.9
    n_rows = n_studies + (1 if has_subgroups else 0) + 2  # +subgroup headers + pooled + blank
    fig_height = max(4.5, header_height + n_rows * row_height + footer_height)
    fig_width = 7.0  # double column

    fig, ax = plt.subplots(figsize=(fig_width, fig_height))
    ax.set_axis_off()

    # ── Column positions (normalized 0–1 in figure space) ────────────────────
    col_label = 0.02
    col_n = 0.38
    col_events = 0.44
    col_weight = 0.51
    col_plot_left = 0.57
    col_plot_right = 0.87
    col_effect = 0.89

    # ── All effect values for x-axis range ───────────────────────────────────
    all_effects = list(df["effect_size"]) + [pooled_effect,
                                              pooled_ci_lower, pooled_ci_upper]
    if log_scale:
        all_log = [np.log(x) for x in all_effects if x > 0]
        x_min_log = min(all_log) - 0.5
        x_max_log = max(all_log) + 0.5
    else:
        x_min_log = min(all_effects) - abs(min(all_effects)) * 0.2
        x_max_log = max(all_effects) + abs(max(all_effects)) * 0.2

    def to_plot_x(value):
        """Map effect value to normalized plot x position."""
        if log_scale and value > 0:
            log_val = np.log(value)
        else:
            log_val = value
        frac = (log_val - x_min_log) / (x_max_log - x_min_log)
        return col_plot_left + frac * (col_plot_right - col_plot_left)

    # ── Headers ───────────────────────────────────────────────────────────────
    fig.text(col_label, 0.96, "Study", fontsize=9, fontweight="bold",
             transform=fig.transFigure, va="top")
    if "n_total" in df.columns:
        fig.text(col_n, 0.96, "N", fontsize=9, fontweight="bold",
                 transform=fig.transFigure, va="top", ha="center")
    if has_weights:
        fig.text(col_weight, 0.96, "Weight\n(%)", fontsize=8, fontweight="bold",
                 transform=fig.transFigure, va="top", ha="center")
    fig.text((col_plot_left + col_plot_right) / 2, 0.96,
             f"{effect_label} (95% CI)", fontsize=9, fontweight="bold",
             transform=fig.transFigure, va="top", ha="center")

    # ── Y positions for each row ───────────────────────────────────────────────
    y_start = 0.90
    y_step = row_height / fig_height
    box_sizes = compute_box_size(df.get("weight"))

    current_subgroup = None
    y = y_start

    rows = []
    if has_subgroups:
        for sg, group in df.groupby("subgroup", sort=False):
            rows.append(("subgroup_header", sg, None))
            for _, row in group.iterrows():
                rows.append(("study", row, None))
    else:
        for _, row in df.iterrows():
            rows.append(("study", row, None))

    rows.append(("blank", None, None))
    rows.append(("pooled", None, None))

    box_idx = 0
    for row_type, data, _ in rows:
        y -= y_step

        if row_type == "subgroup_header":
            fig.text(col_label, y, data, fontsize=9, fontweight="bold",
                     transform=fig.transFigure, va="center",
                     color="#333333", style="italic")

        elif row_type == "study":
            row = data
            # Study label
            fig.text(col_label, y, row["study_label"], fontsize=8,
                     transform=fig.transFigure, va="center")
            # N
            if "n_total" in row.index and pd.notna(row["n_total"]):
                fig.text(col_n, y, f"{int(row['n_total']):,}", fontsize=8,
                         transform=fig.transFigure, va="center", ha="center")
            # Weight
            if "weight" in row.index and pd.notna(row.get("weight")):
                fig.text(col_weight, y, f"{row['weight']:.1f}", fontsize=8,
                         transform=fig.transFigure, va="center", ha="center")
            # CI line
            x_lo = to_plot_x(row["ci_lower"])
            x_hi = to_plot_x(row["ci_upper"])
            x_mid = to_plot_x(row["effect_size"])
            bh = box_sizes[box_idx] * y_step

            # Draw CI line
            fig.add_artist(Line2D([x_lo, x_hi], [y, y],
                                   transform=fig.transFigure,
                                   color=COLORS["study_ci"], linewidth=0.8, zorder=2))
            # Draw study box
            rect = patches.FancyBboxPatch(
                (x_mid - bh * 0.5 * (fig_height / fig_width), y - bh * 0.5),
                bh * (fig_height / fig_width), bh,
                boxstyle="square,pad=0",
                transform=fig.transFigure,
                facecolor=COLORS["study_box"],
                edgecolor=COLORS["study_box"],
                zorder=3,
            )
            fig.add_artist(rect)

            # Effect value text
            val_str = f"{row['effect_size']:.2f} ({row['ci_lower']:.2f}–{row['ci_upper']:.2f})"
            fig.text(col_effect, y, val_str, fontsize=7.5,
                     transform=fig.transFigure, va="center")

            box_idx += 1

        elif row_type == "pooled":
            # Separator line
            y_sep = y + y_step * 0.3
            ax.axhline(y=0, xmin=col_label, xmax=0.98, color="#888888",
                        linewidth=0.5, transform=fig.transFigure, zorder=1)

            # Diamond
            x_lo = to_plot_x(pooled_ci_lower)
            x_hi = to_plot_x(pooled_ci_upper)
            x_mid = to_plot_x(pooled_effect)
            diamond_h = y_step * 0.4
            diamond = plt.Polygon(
                [[x_lo, y], [x_mid, y + diamond_h / 2],
                 [x_hi, y], [x_mid, y - diamond_h / 2]],
                transform=fig.transFigure,
                facecolor=COLORS["pooled_diamond"],
                edgecolor=COLORS["pooled_diamond"],
                zorder=4,
            )
            fig.add_artist(diamond)

            # Pooled estimate text
            val_str = f"{pooled_effect:.2f} ({pooled_ci_lower:.2f}–{pooled_ci_upper:.2f})"
            fig.text(col_label, y, "Pooled estimate", fontsize=8, fontweight="bold",
                     transform=fig.transFigure, va="center")
            fig.text(col_effect, y, val_str, fontsize=7.5, fontweight="bold",
                     transform=fig.transFigure, va="center")

    # ── Null line ─────────────────────────────────────────────────────────────
    x_null = to_plot_x(null_value)
    fig.add_artist(Line2D([x_null, x_null], [y - y_step, y_start + y_step],
                           transform=fig.transFigure,
                           color=COLORS["null_line"], linewidth=0.8,
                           linestyle="--", zorder=1))

    # ── X-axis ticks ─────────────────────────────────────────────────────────
    if log_scale:
        tick_values_raw = [0.25, 0.5, 1.0, 2.0, 4.0]
    else:
        span = x_max_log - x_min_log
        tick_values_raw = np.linspace(x_min_log, x_max_log, 5)

    y_axis = y - y_step * 1.2
    for tv in tick_values_raw:
        try:
            tx = to_plot_x(tv)
        except Exception:
            continue
        if col_plot_left <= tx <= col_plot_right:
            fig.text(tx, y_axis, f"{tv}", fontsize=7.5,
                     transform=fig.transFigure, va="top", ha="center")
            fig.add_artist(Line2D([tx, tx], [y_axis + 0.01, y - y_step * 0.8],
                                   transform=fig.transFigure,
                                   color="#888888", linewidth=0.5))

    # ── Favor labels ─────────────────────────────────────────────────────────
    y_favor = y_axis - 0.025
    fig.text((col_plot_left + x_null) / 2, y_favor, f"← {favor_left}",
             fontsize=7.5, transform=fig.transFigure, va="top", ha="center",
             color="#555555")
    fig.text((x_null + col_plot_right) / 2, y_favor, f"{favor_right} →",
             fontsize=7.5, transform=fig.transFigure, va="top", ha="center",
             color="#555555")

    # ── Heterogeneity footer ──────────────────────────────────────────────────
    y_footer = y_favor - 0.035
    het_str = (
        f"Heterogeneity: I² = {i_squared:.1f}%, τ² = {tau_squared:.3f}, "
        f"Q-test P = {q_p_value:.3f}  |  "
        f"N studies = {n_studies}"
    )
    fig.text(col_label, y_footer, het_str, fontsize=7.5,
             transform=fig.transFigure, va="top", color="#444444")

    # ── Reproducibility footer ────────────────────────────────────────────────
    rep_str = (
        f"Generated: {datetime.now().strftime('%Y-%m-%d')} | "
        f"Script v{SCRIPT_VERSION} | seed={SEED}"
    )
    fig.text(0.99, 0.01, rep_str, fontsize=6, transform=fig.transFigure,
             va="bottom", ha="right", color="#999999")

    # ── Save ──────────────────────────────────────────────────────────────────
    for ext in ["pdf", "png"]:
        outfile = f"{output_path}.{ext}"
        dpi = 300 if ext == "png" else None
        plt.savefig(outfile, dpi=dpi, bbox_inches="tight",
                    facecolor="white", edgecolor="none")
        print(f"Saved: {outfile}")

    plt.close()


def main():
    parser = argparse.ArgumentParser(description="Generate publication-ready forest plot")
    parser.add_argument("--input", required=True, help="Path to CSV file")
    parser.add_argument("--pooled", nargs=3, type=float, metavar=("EFFECT", "CI_LO", "CI_HI"),
                        required=True, help="Pooled effect size and 95% CI")
    parser.add_argument("--i2", type=float, default=0.0, help="I² (%)")
    parser.add_argument("--tau2", type=float, default=0.0, help="τ²")
    parser.add_argument("--q-p", type=float, default=1.0, dest="q_p",
                        help="Cochran Q p-value")
    parser.add_argument("--effect-label", default="OR", help="Label for effect measure")
    parser.add_argument("--null-value", type=float, default=1.0,
                        help="Null value (1.0 for OR/RR; 0 for MD)")
    parser.add_argument("--log-scale", action="store_true", default=True,
                        help="Use log scale for OR/RR (default: True)")
    parser.add_argument("--no-log-scale", action="store_false", dest="log_scale")
    parser.add_argument("--favor-left", default="Favors Treatment")
    parser.add_argument("--favor-right", default="Favors Control")
    parser.add_argument("--output", default="forest_plot", help="Output file path (no extension)")

    args = parser.parse_args()

    df = load_data(args.input)
    pooled_effect, pooled_ci_lo, pooled_ci_hi = args.pooled

    print(f"\n── Forest Plot Generation ──────────────────")
    print(f"Studies: {len(df)}")
    print(f"Pooled {args.effect_label}: {pooled_effect:.2f} "
          f"(95% CI: {pooled_ci_lo:.2f}–{pooled_ci_hi:.2f})")
    print(f"I² = {args.i2:.1f}%, τ² = {args.tau2:.3f}, Q P = {args.q_p:.3f}")

    make_forest_plot(
        df=df,
        pooled_effect=pooled_effect,
        pooled_ci_lower=pooled_ci_lo,
        pooled_ci_upper=pooled_ci_hi,
        i_squared=args.i2,
        tau_squared=args.tau2,
        q_p_value=args.q_p,
        effect_label=args.effect_label,
        null_value=args.null_value,
        log_scale=args.log_scale,
        favor_left=args.favor_left,
        favor_right=args.favor_right,
        output_path=args.output,
    )


# ── Example CSV format ────────────────────────────────────────────────────────
EXAMPLE_CSV = """study_label,effect_size,ci_lower,ci_upper,n_total,weight,subgroup
Kim 2019,0.72,0.51,1.02,234,12.4,Single-center
Park 2020,0.61,0.45,0.83,418,18.2,Single-center
Lee 2021,0.88,0.62,1.24,156,10.1,Single-center
Chen 2022,0.65,0.50,0.85,512,19.8,Multi-center
Wang 2022,0.71,0.53,0.95,345,15.6,Multi-center
Smith 2023,0.58,0.41,0.81,289,14.7,Multi-center
Jones 2023,0.77,0.55,1.08,198,9.2,Multi-center
"""

# ── Run example ───────────────────────────────────────────────────────────────
if __name__ == "__main__":
    if len(sys.argv) == 1:
        # Demo mode: generate example
        import io
        print("Running in demo mode with example data...")
        df = pd.read_csv(io.StringIO(EXAMPLE_CSV))
        make_forest_plot(
            df=df,
            pooled_effect=0.70,
            pooled_ci_lower=0.59,
            pooled_ci_upper=0.83,
            i_squared=23.4,
            tau_squared=0.008,
            q_p_value=0.24,
            effect_label="OR",
            null_value=1.0,
            log_scale=True,
            favor_left="Favors Treatment",
            favor_right="Favors Control",
            output_path="forest_plot_demo",
        )
    else:
        main()

"""
likert_summary.py — Publication-ready Likert Scale Survey Analysis
===================================================================
Analyzes Likert-scale survey data for medical education research papers.

Input:
    CSV with respondent rows and item columns.
    Likert responses should be numeric (1–5 or 1–7).
    Optional group column for subgroup comparisons.

Outputs:
    - Console: formatted descriptive statistics and test results
    - CSV: summary tables (table1_likert.csv)
    - PNG/PDF: diverging stacked bar chart (300 DPI)

Usage:
    python likert_summary.py --input survey_data.csv \
        --items Q1 Q2 Q3 Q4 Q5 \
        --labels "Strongly Disagree" "Disagree" "Neutral" "Agree" "Strongly Agree" \
        --group group_column \
        --scale 5 \
        --output likert_analysis

Dependencies:
    pip install pandas numpy matplotlib scipy pingouin
"""

import argparse
import sys
import warnings
from datetime import datetime

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from scipy import stats

try:
    import pingouin as pg
    PINGOUIN_AVAILABLE = True
except ImportError:
    PINGOUIN_AVAILABLE = False
    warnings.warn("pingouin not installed. Cronbach's alpha will be skipped. "
                  "Install with: pip install pingouin")

# ── Reproducibility ───────────────────────────────────────────────────────────
SCRIPT_VERSION = "1.0.0"
SEED = 42
np.random.seed(SEED)
print(f"Script: likert_summary.py v{SCRIPT_VERSION} | Date: {datetime.now().strftime('%Y-%m-%d')}")

# ── Style ─────────────────────────────────────────────────────────────────────
plt.rcParams.update({
    "font.family": "Arial",
    "font.size": 9,
    "figure.dpi": 150,
})

# Diverging palette: negative → neutral → positive (colorblind-safe)
DIVERGE_COLORS_5 = ["#D55E00", "#E69F00", "#CCCCCC", "#56B4E9", "#0072B2"]
DIVERGE_COLORS_7 = ["#D55E00", "#E69F00", "#F0E442", "#CCCCCC",
                     "#56B4E9", "#009E73", "#0072B2"]


def load_data(filepath: str, item_cols: list) -> pd.DataFrame:
    df = pd.read_csv(filepath)
    missing = [c for c in item_cols if c not in df.columns]
    if missing:
        raise ValueError(f"Item columns not found in CSV: {missing}")
    return df


def descriptive_stats(df: pd.DataFrame, items: list, labels: list) -> pd.DataFrame:
    """Compute descriptive statistics for each Likert item."""
    rows = []
    for item in items:
        col = df[item].dropna()
        rows.append({
            "Item": item,
            "N": len(col),
            "Mean": round(col.mean(), 2),
            "Median": col.median(),
            "SD": round(col.std(), 2),
            "Q1": col.quantile(0.25),
            "Q3": col.quantile(0.75),
            "Min": col.min(),
            "Max": col.max(),
        })
    return pd.DataFrame(rows)


def frequency_table(df: pd.DataFrame, items: list, labels: list, scale: int) -> pd.DataFrame:
    """Compute frequency distribution for each item."""
    scale_values = list(range(1, scale + 1))
    item_labels = labels if labels else [str(v) for v in scale_values]

    rows = []
    for item in items:
        col = df[item].dropna()
        row = {"Item": item, "N": len(col)}
        for val, label in zip(scale_values, item_labels):
            count = (col == val).sum()
            pct = count / len(col) * 100
            row[f"{label} (n)"] = count
            row[f"{label} (%)"] = round(pct, 1)
        rows.append(row)
    return pd.DataFrame(rows)


def apply_reverse_coding(df: pd.DataFrame, items: list, reverse_items: list,
                         scale: int) -> pd.DataFrame:
    """Recode reverse-worded items as (min+max) - x BEFORE reliability/scoring.
    Returns a copy; leaves non-reverse items untouched. min is assumed 1."""
    if not reverse_items:
        return df
    out = df.copy()
    flip_const = scale + 1  # (min=1) + (max=scale)
    for it in reverse_items:
        if it not in items:
            raise ValueError(f"--reverse-items '{it}' is not in --items")
        out[it] = flip_const - out[it]
        print(f"  Reverse-coded: {it} -> ({flip_const} - {it})")
    return out


def item_rest_correlations(df: pd.DataFrame, items: list) -> dict:
    """Corrected item-total (item-rest) correlation per item. A negative value
    means the item moves opposite the rest of the scale — the classic signature
    of a reverse-worded item that was never recoded."""
    sub = df[items].dropna()
    out = {}
    if len(sub) < 2 or len(items) < 2:
        return {it: None for it in items}
    for it in items:
        rest = sub[[c for c in items if c != it]].sum(axis=1)
        if sub[it].std() == 0 or rest.std() == 0:
            out[it] = None
        else:
            out[it] = round(float(sub[it].corr(rest)), 3)
    return out


def cronbach_alpha(df: pd.DataFrame, items: list) -> float:
    """Compute Cronbach's alpha for a set of items."""
    if not PINGOUIN_AVAILABLE:
        return None
    try:
        result = pg.cronbach_alpha(df[items].dropna())
        return round(result[0], 3)
    except Exception as e:
        print(f"  Warning: Could not compute Cronbach's alpha: {e}")
        return None


def compare_groups(df: pd.DataFrame, items: list, group_col: str) -> pd.DataFrame:
    """Mann-Whitney U test for group comparison on each item."""
    groups = df[group_col].dropna().unique()
    if len(groups) != 2:
        print(f"  Warning: Group comparison requires exactly 2 groups; found {len(groups)}")
        return pd.DataFrame()

    g1, g2 = groups
    rows = []
    for item in items:
        d1 = df.loc[df[group_col] == g1, item].dropna()
        d2 = df.loc[df[group_col] == g2, item].dropna()

        stat, p = stats.mannwhitneyu(d1, d2, alternative="two-sided")
        # Rank-biserial r (effect size)
        n1, n2 = len(d1), len(d2)
        r = 1 - (2 * stat) / (n1 * n2)

        rows.append({
            "Item": item,
            f"{g1} Median [IQR]": f"{d1.median():.1f} [{d1.quantile(0.25):.1f}–{d1.quantile(0.75):.1f}]",
            f"{g2} Median [IQR]": f"{d2.median():.1f} [{d2.quantile(0.25):.1f}–{d2.quantile(0.75):.1f}]",
            "U statistic": round(stat, 1),
            "P value": round(p, 3),
            "r (effect size)": round(r, 3),
            "Interpretation": "small" if abs(r) < 0.3 else ("medium" if abs(r) < 0.5 else "large"),
        })
    return pd.DataFrame(rows)


def prepost_comparison(df: pd.DataFrame, pre_items: list, post_items: list) -> pd.DataFrame:
    """Wilcoxon signed-rank test for pre-post comparison."""
    rows = []
    for pre, post in zip(pre_items, post_items):
        paired = df[[pre, post]].dropna()
        d_pre = paired[pre]
        d_post = paired[post]

        stat, p = stats.wilcoxon(d_pre, d_post, alternative="two-sided")
        n = len(paired)
        z = stats.norm.ppf(p / 2)  # approximation
        r = abs(z) / np.sqrt(n)

        rows.append({
            "Item": f"{pre} → {post}",
            "N (paired)": n,
            "Pre Median [IQR]": f"{d_pre.median():.1f} [{d_pre.quantile(0.25):.1f}–{d_pre.quantile(0.75):.1f}]",
            "Post Median [IQR]": f"{d_post.median():.1f} [{d_post.quantile(0.25):.1f}–{d_post.quantile(0.75):.1f}]",
            "W statistic": round(stat, 1),
            "P value": round(p, 3),
            "r (effect size)": round(r, 3),
        })
    return pd.DataFrame(rows)


def diverging_bar_chart(
    df: pd.DataFrame,
    items: list,
    labels: list,
    scale: int,
    group_col: str = None,
    output_path: str = "likert_chart",
    item_names: dict = None,
) -> None:
    """Generate diverging stacked bar chart."""
    colors = DIVERGE_COLORS_5 if scale == 5 else DIVERGE_COLORS_7
    scale_values = list(range(1, scale + 1))
    neutral_idx = scale // 2  # 0-indexed neutral position

    plot_df = df[items].copy()

    # Compute proportions for each item
    prop_data = []
    for item in items:
        col = plot_df[item].dropna()
        n = len(col)
        props = [(col == v).sum() / n for v in scale_values]
        prop_data.append(props)

    prop_arr = np.array(prop_data)  # shape: (n_items, scale)

    fig, ax = plt.subplots(figsize=(7.0, max(3.0, len(items) * 0.55 + 1.5)))

    y_pos = np.arange(len(items))
    height = 0.6

    # Diverging: negative on left, positive on right
    neg_props = prop_arr[:, :neutral_idx]     # columns before neutral
    neutral_prop = prop_arr[:, neutral_idx]   # neutral column
    pos_props = prop_arr[:, neutral_idx + 1:] # columns after neutral

    # Starting x for left bars (negative, goes left from center)
    left_start = -(neg_props.sum(axis=1) + neutral_prop / 2)

    # Draw negative bars (left to right, so right-most color first)
    x_left = left_start.copy()
    for i in range(neg_props.shape[1] - 1, -1, -1):
        ax.barh(y_pos, neg_props[:, i], left=x_left, height=height,
                color=colors[i], edgecolor="white", linewidth=0.5)
        x_left += neg_props[:, i]

    # Draw neutral bar
    x_neutral = -(neutral_prop / 2)
    ax.barh(y_pos, neutral_prop, left=x_neutral, height=height,
            color=colors[neutral_idx], edgecolor="white", linewidth=0.5)

    # Draw positive bars
    x_right = neutral_prop / 2
    for i, j in enumerate(range(neutral_idx + 1, scale)):
        ax.barh(y_pos, pos_props[:, i], left=x_right, height=height,
                color=colors[j], edgecolor="white", linewidth=0.5)
        x_right += pos_props[:, i]

    # Percentage labels on bars > 10%
    # (omitted for brevity — add if needed)

    # Axes formatting
    ax.axvline(0, color="black", linewidth=0.8, zorder=5)
    ax.set_yticks(y_pos)
    y_labels = [item_names.get(item, item) if item_names else item for item in items]
    ax.set_yticklabels(y_labels, fontsize=8)
    ax.set_xlabel("Proportion of respondents", fontsize=9)

    # X-axis: convert to percentage
    xticks = ax.get_xticks()
    ax.set_xticklabels([f"{abs(x)*100:.0f}%" for x in xticks], fontsize=8)

    # Legend
    legend_patches = [mpatches.Patch(color=colors[i], label=labels[i])
                      for i in range(scale)]
    ax.legend(handles=legend_patches, loc="lower center",
              bbox_to_anchor=(0.5, -0.25), ncol=scale, fontsize=7.5,
              frameon=False)

    ax.set_xlim(-1.0, 1.0)
    ax.grid(axis="x", color="#DDDDDD", linewidth=0.5)
    ax.spines[["top", "right"]].set_visible(False)

    plt.tight_layout()

    for ext in ["pdf", "png"]:
        outfile = f"{output_path}_diverging.{ext}"
        dpi = 300 if ext == "png" else None
        plt.savefig(outfile, dpi=dpi, bbox_inches="tight",
                    facecolor="white", edgecolor="none")
        print(f"Saved: {outfile}")
    plt.close()


def main():
    parser = argparse.ArgumentParser(description="Likert scale analysis")
    parser.add_argument("--input", required=True)
    parser.add_argument("--items", nargs="+", required=True)
    parser.add_argument("--labels", nargs="+",
                        default=["Strongly Disagree", "Disagree", "Neutral",
                                  "Agree", "Strongly Agree"])
    parser.add_argument("--scale", type=int, default=5, choices=[5, 7])
    parser.add_argument("--reverse-items", nargs="+", default=None,
                        help="reverse-worded items to recode as (scale+1)-x "
                             "BEFORE alpha/scoring (e.g. a 'I review critically' item)")
    parser.add_argument("--group", default=None)
    parser.add_argument("--pre-items", nargs="+", default=None)
    parser.add_argument("--post-items", nargs="+", default=None)
    parser.add_argument("--output", default="likert_analysis")

    args = parser.parse_args()

    df = load_data(args.input, args.items)
    print(f"\nLoaded: {len(df)} respondents, {len(args.items)} items")

    # ── Reverse-coding guard (run BEFORE any scoring/reliability) ──────────────
    # A reverse-worded item must be recoded or it sinks Cronbach's alpha (often
    # negative). See ~/.claude/rules/survey-scale-reliability.md.
    if args.reverse_items:
        print("\n── Reverse-coding applied ──────────────────")
        df = apply_reverse_coding(df, args.items, args.reverse_items, args.scale)

    # ── Descriptive statistics ────────────────────────────────────────────────
    print("\n── Descriptive Statistics ──────────────────")
    desc = descriptive_stats(df, args.items, args.labels)
    print(desc.to_string(index=False))

    # ── Frequency table ───────────────────────────────────────────────────────
    print("\n── Frequency Distribution ──────────────────")
    freq = frequency_table(df, args.items, args.labels, args.scale)
    print(freq.to_string(index=False))

    # ── Internal consistency + reverse-coding guard ───────────────────────────
    if len(args.items) >= 2:
        rests = item_rest_correlations(df, args.items)
        suspects = [it for it, r in rests.items() if r is not None and r < 0]
        if suspects:
            print("\n── Reverse-coding check ────────────────────")
            for it in args.items:
                mark = "  <-- reverse-code suspect (negative item-rest)" if it in suspects else ""
                print(f"   item-rest r  {it:<16} {rests[it]}{mark}")
            print(f"   ⚠ Suspect items: {', '.join(suspects)}. If reverse-worded, "
                  "rerun with --reverse-items before reporting alpha.")

    if PINGOUIN_AVAILABLE:
        alpha = cronbach_alpha(df, args.items)
        if alpha is not None:
            print(f"\n── Cronbach's Alpha: {alpha}")
            if alpha < 0:
                # A negative alpha is almost never a real measurement phenomenon;
                # it is a reverse-coding bug. Do NOT defend it as multidimensional.
                print("   ⚠ NEGATIVE alpha — almost always a reverse-coding bug, "
                      "not a 'multidimensional' construct.")
                print("     Recode reverse-worded items (--reverse-items) and rerun "
                      "BEFORE interpreting. See survey-scale-reliability.md.")
            else:
                interp = ("excellent" if alpha >= 0.9 else
                          "good" if alpha >= 0.8 else
                          "acceptable" if alpha >= 0.7 else
                          "questionable" if alpha >= 0.6 else "poor")
                print(f"   Interpretation: {interp}")

    # ── Group comparison ──────────────────────────────────────────────────────
    if args.group and args.group in df.columns:
        print(f"\n── Group Comparison ({args.group}) — Mann-Whitney U ──")
        grp = compare_groups(df, args.items, args.group)
        if not grp.empty:
            print(grp.to_string(index=False))
            grp.to_csv(f"{args.output}_group_comparison.csv", index=False)
            print(f"Saved: {args.output}_group_comparison.csv")

    # ── Pre-post comparison ───────────────────────────────────────────────────
    if args.pre_items and args.post_items:
        print("\n── Pre-Post Comparison — Wilcoxon Signed-Rank ──")
        pp = prepost_comparison(df, args.pre_items, args.post_items)
        print(pp.to_string(index=False))
        pp.to_csv(f"{args.output}_prepost.csv", index=False)
        print(f"Saved: {args.output}_prepost.csv")

    # ── Save tables ───────────────────────────────────────────────────────────
    desc.to_csv(f"{args.output}_descriptive.csv", index=False)
    freq.to_csv(f"{args.output}_frequency.csv", index=False)
    print(f"\nSaved: {args.output}_descriptive.csv")
    print(f"Saved: {args.output}_frequency.csv")

    # ── Diverging bar chart ───────────────────────────────────────────────────
    diverging_bar_chart(df, args.items, args.labels, args.scale,
                         output_path=args.output)

    print("\n── Session Info ─────────────────────────────")
    print(f"Python: {sys.version.split()[0]}")
    for pkg in ["pandas", "numpy", "scipy", "pingouin"]:
        try:
            import importlib
            m = importlib.import_module(pkg)
            print(f"  {pkg}: {m.__version__}")
        except Exception:
            print(f"  {pkg}: not installed")


if __name__ == "__main__":
    if len(sys.argv) == 1:
        print("No arguments provided. Run with --help for usage.")
        print("\nExample:")
        print("  python likert_summary.py --input survey.csv "
              "--items Q1 Q2 Q3 Q4 Q5 --scale 5 --group specialty")
    else:
        main()

"""
Template: Propensity Score Analysis
Supports PS matching, IPTW, and overlap weighting for observational studies.
Generates balance tables, Love plots, and weighted outcome analyses.

Usage:
    Modify the CONFIGURATION section below, then run:
        python propensity_score.py

Input:  CSV with treatment indicator, outcome, and covariate columns
Output: balance table CSV, Love plot PDF/PNG, PS distribution plot, outcome results
"""

# === REPRODUCIBILITY HEADER ===
import sys
import os
import datetime
import numpy as np
import pandas as pd
from scipy import stats

np.random.seed(42)
print(f"Date: {datetime.date.today()}")
print(f"Python: {sys.version}")
print(f"numpy: {np.__version__}, pandas: {pd.__version__}, scipy: {stats.scipy.__version__}")

try:
    from sklearn.linear_model import LogisticRegression
    from sklearn.neighbors import NearestNeighbors
    import sklearn
    print(f"sklearn: {sklearn.__version__}")
except ImportError:
    print("Error: scikit-learn not installed. Install with: pip install scikit-learn")
    sys.exit(1)

try:
    import statsmodels.api as sm
    print(f"statsmodels: {sm.__version__}")
except ImportError:
    print("Error: statsmodels not installed. Install with: pip install statsmodels")
    sys.exit(1)

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

STYLE_PATH = os.path.join(os.path.dirname(os.path.dirname(__file__)), "style", "figure_style.mplstyle")
if os.path.exists(STYLE_PATH):
    plt.style.use(STYLE_PATH)


# === CONFIGURATION ===
CONFIG = {
    # Data
    "data_path": "data.csv",
    "output_dir": ".",

    # Variables
    "treatment": "treatment",       # Binary treatment indicator (0/1)
    "outcome": "outcome",           # Outcome variable
    "outcome_type": "continuous",    # "continuous" or "binary"
    "covariates": ["age", "sex", "bmi", "comorbidity_score"],
    "categorical_covariates": ["sex"],

    # PS method: "matching", "iptw", "siptw", "overlap"
    "ps_method": "matching",

    # Matching options
    "caliper_sd_multiplier": 0.2,   # caliper = 0.2 * SD(logit PS)
    "matching_ratio": 1,            # 1:1 matching

    # Balance threshold
    "smd_threshold": 0.10,

    # IPTW options
    "stabilized_weights": True,
    "weight_truncation": 10.0,      # truncate weights > this value
}


# === HELPER FUNCTIONS ===

def estimate_ps(df, treatment_col, covariates):
    """Estimate propensity scores using logistic regression."""
    X = df[covariates].values
    y = df[treatment_col].values
    model = LogisticRegression(max_iter=1000, random_state=42)
    model.fit(X, y)
    ps = model.predict_proba(X)[:, 1]
    return ps


def calculate_smd(x1, x0, is_binary=False):
    """Calculate standardized mean difference."""
    if is_binary:
        p1, p0 = x1.mean(), x0.mean()
        denom = np.sqrt((p1 * (1 - p1) + p0 * (1 - p0)) / 2)
        if denom == 0:
            return 0.0
        return (p1 - p0) / denom
    else:
        denom = np.sqrt((x1.var() + x0.var()) / 2)
        if denom == 0:
            return 0.0
        return (x1.mean() - x0.mean()) / denom


def balance_table(df, treatment_col, covariates, categorical_covs, weights=None):
    """Generate balance table with SMD before/after adjustment."""
    treated = df[treatment_col] == 1
    results = []

    for var in covariates:
        is_cat = var in categorical_covs
        x1 = df.loc[treated, var]
        x0 = df.loc[~treated, var]

        if weights is not None:
            w1 = weights[treated]
            w0 = weights[~treated]
            # Weighted means
            wm1 = np.average(x1, weights=w1)
            wm0 = np.average(x0, weights=w0)
            # Weighted SMD (approximate)
            if is_cat:
                denom = np.sqrt((wm1 * (1 - wm1) + wm0 * (1 - wm0)) / 2)
            else:
                wv1 = np.average((x1 - wm1) ** 2, weights=w1)
                wv0 = np.average((x0 - wm0) ** 2, weights=w0)
                denom = np.sqrt((wv1 + wv0) / 2)
            smd_adj = (wm1 - wm0) / denom if denom > 0 else 0.0
        else:
            smd_adj = None

        smd_raw = calculate_smd(x1, x0, is_binary=is_cat)

        row = {
            "Variable": var,
            "Treated_mean": x1.mean(),
            "Treated_sd": x1.std(),
            "Control_mean": x0.mean(),
            "Control_sd": x0.std(),
            "SMD_before": abs(smd_raw),
        }
        if smd_adj is not None:
            row["SMD_after"] = abs(smd_adj)
        results.append(row)

    return pd.DataFrame(results)


def plot_love(bal_df, smd_threshold, output_dir, title="Love Plot"):
    """Generate Love plot comparing SMD before and after adjustment."""
    fig, ax = plt.subplots(figsize=(8, max(3, len(bal_df) * 0.5)))

    y_pos = range(len(bal_df))
    ax.scatter(bal_df["SMD_before"], y_pos, marker="o", color="gray",
               s=60, label="Before", zorder=3)
    if "SMD_after" in bal_df.columns:
        ax.scatter(bal_df["SMD_after"], y_pos, marker="s", color="navy",
                   s=60, label="After", zorder=4)

    ax.axvline(x=smd_threshold, color="red", linestyle="--", alpha=0.7,
               label=f"Threshold ({smd_threshold})")
    ax.set_yticks(y_pos)
    ax.set_yticklabels(bal_df["Variable"])
    ax.set_xlabel("Absolute Standardized Mean Difference")
    ax.set_title(title)
    ax.legend(loc="best")
    ax.set_xlim(left=0)

    plt.tight_layout()
    for ext in ["pdf", "png"]:
        fig.savefig(os.path.join(output_dir, f"love_plot.{ext}"),
                    dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: love_plot.pdf/.png")


def plot_ps_distribution(ps, treatment, output_dir):
    """Plot PS distribution by treatment group."""
    fig, ax = plt.subplots(figsize=(8, 5))

    ax.hist(ps[treatment == 1], bins=30, alpha=0.5, density=True,
            color="steelblue", label="Treated")
    ax.hist(ps[treatment == 0], bins=30, alpha=0.5, density=True,
            color="coral", label="Control")
    ax.set_xlabel("Propensity Score")
    ax.set_ylabel("Density")
    ax.set_title("Propensity Score Distribution")
    ax.legend()

    plt.tight_layout()
    for ext in ["pdf", "png"]:
        fig.savefig(os.path.join(output_dir, f"ps_distribution.{ext}"),
                    dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: ps_distribution.pdf/.png")


def ps_matching(df, ps, treatment_col, caliper_sd_mult=0.2, ratio=1):
    """Perform 1:N nearest-neighbor PS matching with caliper."""
    logit_ps = np.log(ps / (1 - ps))
    caliper = caliper_sd_mult * logit_ps.std()

    treated_idx = df.index[df[treatment_col] == 1].values
    control_idx = df.index[df[treatment_col] == 0].values

    treated_logit = logit_ps[treated_idx].reshape(-1, 1)
    control_logit = logit_ps[control_idx].reshape(-1, 1)

    nn = NearestNeighbors(n_neighbors=ratio, metric="euclidean")
    nn.fit(control_logit)
    distances, indices = nn.kneighbors(treated_logit)

    matched_treated = []
    matched_control = []
    used_controls = set()

    for i, (dist_arr, idx_arr) in enumerate(zip(distances, indices)):
        for d, j in zip(dist_arr, idx_arr):
            ctrl_orig_idx = control_idx[j]
            if d <= caliper and ctrl_orig_idx not in used_controls:
                matched_treated.append(treated_idx[i])
                matched_control.append(ctrl_orig_idx)
                used_controls.add(ctrl_orig_idx)
                break

    matched_indices = matched_treated + matched_control
    n_unmatched = len(treated_idx) - len(matched_treated)

    print(f"\nPS Matching Results:")
    print(f"  Caliper: {caliper:.4f} (= {caliper_sd_mult} x SD(logit PS))")
    print(f"  Matched pairs: {len(matched_treated)}")
    print(f"  Unmatched treated: {n_unmatched}")
    print(f"  Unmatched controls: {len(control_idx) - len(matched_control)}")

    return df.loc[matched_indices].copy(), matched_indices


def iptw_weights(ps, treatment, stabilized=True, truncation=10.0):
    """Calculate IPTW weights (ATE)."""
    if stabilized:
        p_treat = treatment.mean()
        w = np.where(treatment == 1, p_treat / ps, (1 - p_treat) / (1 - ps))
    else:
        w = np.where(treatment == 1, 1 / ps, 1 / (1 - ps))

    # Truncation
    n_truncated = (w > truncation).sum()
    if n_truncated > 0:
        print(f"  Truncated {n_truncated} weights > {truncation}")
        w = np.clip(w, None, truncation)

    print(f"\nIPTW Weights Summary:")
    print(f"  Mean: {w.mean():.2f}, SD: {w.std():.2f}")
    print(f"  Min: {w.min():.2f}, Max: {w.max():.2f}")
    print(f"  Stabilized: {stabilized}")

    return w


def siptw_weights(ps, treatment, truncation=10.0):
    """Calculate Stabilized Inverse Probability of Treatment Weights (SIPTW).

    SIPTW maintains the sample size of the entire cohort and allows for
    appropriate estimation of the variance of the main effect. Increasingly
    used in emulated target trial frameworks (Yon DK group pattern).

    Weights:
        Treated:   P(T=1) / PS
        Control:   P(T=0) / (1 - PS)

    This is equivalent to stabilized IPTW but explicitly named SIPTW in some
    literature to distinguish from unstabilized IPTW.
    """
    p_treat = treatment.mean()
    w = np.where(treatment == 1, p_treat / ps, (1 - p_treat) / (1 - ps))

    # Truncation
    n_truncated = (w > truncation).sum()
    if n_truncated > 0:
        print(f"  Truncated {n_truncated} weights > {truncation}")
        w = np.clip(w, None, truncation)

    # Effective sample size
    ess_treated = (w[treatment == 1].sum()) ** 2 / (w[treatment == 1] ** 2).sum()
    ess_control = (w[treatment == 0].sum()) ** 2 / (w[treatment == 0] ** 2).sum()

    print(f"\nSIPTW Weights Summary:")
    print(f"  Mean: {w.mean():.2f}, SD: {w.std():.2f}")
    print(f"  Min: {w.min():.2f}, Max: {w.max():.2f}")
    print(f"  Effective sample size (treated): {ess_treated:.0f} / {(treatment == 1).sum()}")
    print(f"  Effective sample size (control): {ess_control:.0f} / {(treatment == 0).sum()}")

    return w


def overlap_weights(ps, treatment):
    """Calculate overlap weights (ATO)."""
    w = np.where(treatment == 1, 1 - ps, ps)

    print(f"\nOverlap Weights Summary:")
    print(f"  Mean: {w.mean():.3f}, SD: {w.std():.3f}")
    print(f"  Min: {w.min():.3f}, Max: {w.max():.3f}")

    return w


def weighted_outcome_analysis(df, treatment_col, outcome_col, outcome_type, weights):
    """Perform weighted outcome analysis."""
    import statsmodels.api as sm

    X = sm.add_constant(df[[treatment_col]])
    y = df[outcome_col]

    if outcome_type == "binary":
        model = sm.GLM(y, X, family=sm.families.Binomial(), freq_weights=weights)
    else:
        model = sm.WLS(y, X, weights=weights)

    result = model.fit()
    print(f"\n--- Weighted Outcome Analysis ---")
    print(result.summary2())

    # Extract treatment effect
    coef = result.params[treatment_col]
    ci = result.conf_int().loc[treatment_col]
    p_val = result.pvalues[treatment_col]

    if outcome_type == "binary":
        or_val = np.exp(coef)
        or_ci = np.exp(ci)
        print(f"\nTreatment effect (OR): {or_val:.2f} (95% CI: {or_ci[0]:.2f}-{or_ci[1]:.2f}), P = {p_val:.3f}")
    else:
        print(f"\nTreatment effect (β): {coef:.3f} (95% CI: {ci[0]:.3f}-{ci[1]:.3f}), P = {p_val:.3f}")

    return result


# === MAIN ANALYSIS ===

def main():
    config = CONFIG
    df = pd.read_csv(config["data_path"])
    output_dir = config["output_dir"]
    print(f"Data loaded: {df.shape[0]} rows x {df.shape[1]} columns")

    treatment_col = config["treatment"]
    outcome_col = config["outcome"]
    covariates = config["covariates"]

    # Encode categorical variables
    for cat_var in config.get("categorical_covariates", []):
        if cat_var in df.columns and df[cat_var].dtype == "object":
            df[cat_var] = pd.Categorical(df[cat_var]).codes

    # Drop missing
    analysis_vars = [treatment_col, outcome_col] + covariates
    n_before = len(df)
    df = df.dropna(subset=analysis_vars)
    n_after = len(df)
    if n_before != n_after:
        print(f"Excluded {n_before - n_after} rows with missing data ({100*(n_before-n_after)/n_before:.1f}%)")

    treatment = df[treatment_col].values
    n_treated = treatment.sum()
    n_control = len(treatment) - n_treated

    print(f"\n{'='*60}")
    print(f"PROPENSITY SCORE ANALYSIS")
    print(f"{'='*60}")
    print(f"Method: {config['ps_method'].upper()}")
    print(f"Treated: {n_treated}, Control: {n_control}")
    print(f"Covariates: {len(covariates)}")

    # Step 1: Estimate PS
    print(f"\n--- Step 1: PS Estimation ---")
    ps = estimate_ps(df, treatment_col, covariates)
    df["ps"] = ps
    plot_ps_distribution(ps, treatment, output_dir)

    # Pre-adjustment balance
    print(f"\n--- Pre-adjustment Balance ---")
    bal_before = balance_table(df, treatment_col, covariates,
                               config.get("categorical_covariates", []))
    print(bal_before[["Variable", "SMD_before"]].to_string(index=False))
    n_imbalanced = (bal_before["SMD_before"] > config["smd_threshold"]).sum()
    print(f"Variables with SMD > {config['smd_threshold']}: {n_imbalanced}/{len(covariates)}")

    # Step 2: Apply PS method
    weights = None
    if config["ps_method"] == "matching":
        print(f"\n--- Step 2: PS Matching ---")
        df_matched, _ = ps_matching(
            df, ps, treatment_col,
            caliper_sd_mult=config["caliper_sd_multiplier"],
            ratio=config["matching_ratio"]
        )
        # Balance after matching
        bal_after = balance_table(df_matched, treatment_col, covariates,
                                  config.get("categorical_covariates", []))
        bal_combined = bal_before.copy()
        bal_combined["SMD_after"] = bal_after["SMD_before"].values

    elif config["ps_method"] == "iptw":
        print(f"\n--- Step 2: IPTW ---")
        weights = iptw_weights(ps, treatment,
                                stabilized=config["stabilized_weights"],
                                truncation=config["weight_truncation"])
        df["weights"] = weights
        bal_combined = balance_table(df, treatment_col, covariates,
                                     config.get("categorical_covariates", []),
                                     weights=weights)

    elif config["ps_method"] == "siptw":
        print(f"\n--- Step 2: SIPTW (Stabilized Inverse Probability of Treatment Weighting) ---")
        weights = siptw_weights(ps, treatment,
                                truncation=config["weight_truncation"])
        df["weights"] = weights
        bal_combined = balance_table(df, treatment_col, covariates,
                                     config.get("categorical_covariates", []),
                                     weights=weights)

    elif config["ps_method"] == "overlap":
        print(f"\n--- Step 2: Overlap Weighting ---")
        weights = overlap_weights(ps, treatment)
        df["weights"] = weights
        bal_combined = balance_table(df, treatment_col, covariates,
                                     config.get("categorical_covariates", []),
                                     weights=weights)

    # Step 3: Balance assessment
    print(f"\n--- Step 3: Post-adjustment Balance ---")
    print(bal_combined[["Variable", "SMD_before", "SMD_after"]].to_string(index=False))
    n_imbalanced_after = (bal_combined["SMD_after"] > config["smd_threshold"]).sum()
    print(f"Variables with SMD > {config['smd_threshold']} after adjustment: "
          f"{n_imbalanced_after}/{len(covariates)}")

    if n_imbalanced_after > 0:
        print("⚠ WARNING: Some covariates remain imbalanced. "
              "Consider adding interaction terms to PS model or switching method.")

    # Love plot
    plot_love(bal_combined, config["smd_threshold"], output_dir)

    # Save balance table
    bal_combined.to_csv(os.path.join(output_dir, "balance_table.csv"), index=False)
    print(f"Saved: balance_table.csv")

    # Step 4: Outcome analysis
    print(f"\n--- Step 4: Outcome Analysis ---")
    if config["ps_method"] == "matching":
        # Simple comparison in matched data
        t_outcome = df_matched.loc[df_matched[treatment_col] == 1, outcome_col]
        c_outcome = df_matched.loc[df_matched[treatment_col] == 0, outcome_col]
        if config["outcome_type"] == "continuous":
            stat, p = stats.ttest_ind(t_outcome, c_outcome)
            diff = t_outcome.mean() - c_outcome.mean()
            print(f"Mean difference: {diff:.3f}")
            print(f"Treated: {t_outcome.mean():.3f} ± {t_outcome.std():.3f}")
            print(f"Control: {c_outcome.mean():.3f} ± {c_outcome.std():.3f}")
            print(f"t = {stat:.3f}, P = {p:.3f}")
        else:
            # For binary outcome in matched data
            tab = pd.crosstab(df_matched[treatment_col], df_matched[outcome_col])
            print(tab)
    else:
        # Weighted analysis for IPTW/OW
        weighted_outcome_analysis(df, treatment_col, outcome_col,
                                   config["outcome_type"], weights)

    print(f"\n{'='*60}")
    print("Propensity score analysis complete.")


if __name__ == "__main__":
    main()

"""
Template: Regression Analysis (Logistic + Linear)
Performs logistic regression (binary outcome) or multiple linear regression (continuous outcome).
Generates OR/coefficient tables, model diagnostics, and publication-ready figures.

Usage:
    Modify the CONFIGURATION section below, then run:
        python regression.py

Input:  CSV with outcome and predictor variables
Output: coefficient/OR table CSV, diagnostic plots PDF/PNG, summary text
"""

# === REPRODUCIBILITY HEADER ===
import sys
import os
import datetime
import numpy as np
import pandas as pd
from scipy import stats

np.random.seed(42)
print(f"Date: {datetime.date.today()}")
print(f"Python: {sys.version}")
print(f"numpy: {np.__version__}, pandas: {pd.__version__}, scipy: {stats.scipy.__version__}")

try:
    import statsmodels.api as sm
    from statsmodels.stats.outliers_influence import variance_inflation_factor
    print(f"statsmodels: {sm.__version__}")
except ImportError:
    print("Error: statsmodels not installed. Install with: pip install statsmodels")
    sys.exit(1)

try:
    from sklearn.metrics import roc_auc_score, brier_score_loss
    import sklearn
    print(f"sklearn: {sklearn.__version__}")
except ImportError:
    print("Warning: scikit-learn not installed. Some metrics unavailable.")

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

STYLE_PATH = os.path.join(os.path.dirname(os.path.dirname(__file__)), "style", "figure_style.mplstyle")
if os.path.exists(STYLE_PATH):
    plt.style.use(STYLE_PATH)


# === CONFIGURATION ===
CONFIG = {
    # Data
    "data_path": "data.csv",
    "output_dir": ".",

    # Regression type: "logistic" or "linear"
    "regression_type": "logistic",

    # Variables
    "outcome": "event",
    "predictors": ["age", "sex", "bmi", "smoking"],
    "categorical_vars": ["sex", "smoking"],

    # Options
    "run_univariable": True,  # Run univariable analysis before multivariable
    "vif_threshold": 5.0,
    "epv_minimum": 10,
}


# === HELPER FUNCTIONS ===

def calculate_vif(X):
    """Calculate VIF for each predictor."""
    vif_data = pd.DataFrame()
    vif_data["Variable"] = X.columns
    vif_data["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    return vif_data


def logistic_or_table(model, var_names):
    """Generate OR table from logistic regression results."""
    or_vals = np.exp(model.params)
    ci = np.exp(model.conf_int())
    table = pd.DataFrame({
        "Variable": var_names,
        "OR": or_vals,
        "CI_lower": ci[0],
        "CI_upper": ci[1],
        "P_value": model.pvalues
    })
    table["OR_CI"] = table.apply(
        lambda r: f"{r['OR']:.2f} ({r['CI_lower']:.2f}-{r['CI_upper']:.2f})", axis=1
    )
    return table


def linear_coef_table(model, var_names):
    """Generate coefficient table from linear regression results."""
    ci = model.conf_int()
    table = pd.DataFrame({
        "Variable": var_names,
        "Coefficient": model.params,
        "CI_lower": ci[0],
        "CI_upper": ci[1],
        "P_value": model.pvalues
    })
    table["Coef_CI"] = table.apply(
        lambda r: f"{r['Coefficient']:.3f} ({r['CI_lower']:.3f} to {r['CI_upper']:.3f})", axis=1
    )
    return table


def hosmer_lemeshow_test(y_true, y_pred, n_groups=10):
    """Hosmer-Lemeshow goodness-of-fit test."""
    data = pd.DataFrame({"y": y_true, "p": y_pred})
    data["group"] = pd.qcut(data["p"], n_groups, duplicates="drop")
    grouped = data.groupby("group").agg(
        obs=("y", "sum"),
        exp=("p", "sum"),
        n=("y", "count"),
        mean_p=("p", "mean")
    )
    chi2 = ((grouped["obs"] - grouped["exp"]) ** 2 /
            (grouped["exp"] * (1 - grouped["mean_p"]))).sum()
    df = len(grouped) - 2
    p_value = 1 - stats.chi2.cdf(chi2, df)
    return chi2, df, p_value


def plot_diagnostic_4panel(model, outcome_name, output_dir):
    """Generate 4-panel diagnostic plot for linear regression."""
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))

    fitted = model.fittedvalues
    residuals = model.resid
    std_resid = model.get_influence().resid_studentized_internal
    leverage = model.get_influence().hat_matrix_diag
    cooks_d = model.get_influence().cooks_distance[0]

    # 1. Residuals vs Fitted
    axes[0, 0].scatter(fitted, residuals, alpha=0.5, s=20)
    axes[0, 0].axhline(y=0, color="red", linestyle="--")
    axes[0, 0].set_xlabel("Fitted values")
    axes[0, 0].set_ylabel("Residuals")
    axes[0, 0].set_title("Residuals vs Fitted")

    # 2. Q-Q plot
    stats.probplot(std_resid, plot=axes[0, 1])
    axes[0, 1].set_title("Normal Q-Q")

    # 3. Scale-Location
    axes[1, 0].scatter(fitted, np.sqrt(np.abs(std_resid)), alpha=0.5, s=20)
    axes[1, 0].set_xlabel("Fitted values")
    axes[1, 0].set_ylabel("√|Standardized residuals|")
    axes[1, 0].set_title("Scale-Location")

    # 4. Residuals vs Leverage
    axes[1, 1].scatter(leverage, std_resid, alpha=0.5, s=20)
    axes[1, 1].axhline(y=0, color="red", linestyle="--")
    # Cook's distance contours
    n = len(fitted)
    threshold = 4 / n
    high_cook = cooks_d > threshold
    if high_cook.any():
        axes[1, 1].scatter(leverage[high_cook], std_resid[high_cook],
                           color="red", s=40, zorder=5, label=f"Cook's D > {threshold:.3f}")
        axes[1, 1].legend()
    axes[1, 1].set_xlabel("Leverage")
    axes[1, 1].set_ylabel("Standardized residuals")
    axes[1, 1].set_title("Residuals vs Leverage")

    plt.suptitle(f"Diagnostic Plots: {outcome_name}", fontsize=14, y=1.02)
    plt.tight_layout()

    for ext in ["pdf", "png"]:
        fig.savefig(os.path.join(output_dir, f"diagnostic_plots.{ext}"),
                    dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: diagnostic_plots.pdf/.png")

    # Report influential observations
    n_influential = high_cook.sum()
    if n_influential > 0:
        print(f"\nWarning: {n_influential} observation(s) with Cook's D > {threshold:.3f}")


def plot_forest_or(or_table, output_dir, title="Multivariable Logistic Regression"):
    """Forest plot for odds ratios."""
    # Exclude intercept
    plot_data = or_table[or_table["Variable"] != "const"].copy()
    plot_data = plot_data.iloc[::-1]  # reverse for top-down display

    fig, ax = plt.subplots(figsize=(8, max(3, len(plot_data) * 0.6)))

    y_pos = range(len(plot_data))
    ax.errorbar(
        plot_data["OR"], y_pos,
        xerr=[plot_data["OR"] - plot_data["CI_lower"],
              plot_data["CI_upper"] - plot_data["OR"]],
        fmt="o", color="navy", capsize=4, markersize=6
    )
    ax.axvline(x=1, color="red", linestyle="--", alpha=0.7)
    ax.set_yticks(y_pos)
    ax.set_yticklabels(plot_data["Variable"])
    ax.set_xlabel("Odds Ratio (95% CI)")
    ax.set_title(title)
    ax.set_xscale("log")

    plt.tight_layout()
    for ext in ["pdf", "png"]:
        fig.savefig(os.path.join(output_dir, f"forest_plot_or.{ext}"),
                    dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: forest_plot_or.pdf/.png")


# === MAIN ANALYSIS ===

def run_logistic(df, config):
    """Run logistic regression analysis."""
    outcome = config["outcome"]
    predictors = config["predictors"]
    output_dir = config["output_dir"]

    y = df[outcome]
    n_events = y.sum()
    n_total = len(y)
    epv = n_events / len(predictors)

    print(f"\n{'='*60}")
    print(f"LOGISTIC REGRESSION")
    print(f"{'='*60}")
    print(f"Outcome: {outcome}")
    print(f"N = {n_total}, Events = {n_events} ({100*n_events/n_total:.1f}%)")
    print(f"Predictors: {len(predictors)}")
    print(f"EPV = {epv:.1f} (minimum recommended: {config['epv_minimum']})")
    if epv < config["epv_minimum"]:
        print(f"⚠ WARNING: EPV < {config['epv_minimum']}. Model may be unstable. "
              "Consider reducing predictors or using penalized regression.")

    # --- Univariable analysis ---
    if config["run_univariable"]:
        print(f"\n--- Univariable Analysis ---")
        uni_results = []
        for var in predictors:
            X_uni = sm.add_constant(df[[var]])
            try:
                model_uni = sm.Logit(y, X_uni).fit(disp=0)
                or_val = np.exp(model_uni.params[var])
                ci = np.exp(model_uni.conf_int().loc[var])
                p_val = model_uni.pvalues[var]
                uni_results.append({
                    "Variable": var,
                    "Uni_OR": or_val,
                    "Uni_CI_lower": ci[0],
                    "Uni_CI_upper": ci[1],
                    "Uni_P": p_val,
                    "Uni_OR_CI": f"{or_val:.2f} ({ci[0]:.2f}-{ci[1]:.2f})"
                })
            except Exception as e:
                print(f"  {var}: failed ({e})")
                uni_results.append({"Variable": var, "Uni_OR": np.nan})
        uni_df = pd.DataFrame(uni_results)
        print(uni_df[["Variable", "Uni_OR_CI", "Uni_P"]].to_string(index=False))

    # --- Multivariable analysis ---
    print(f"\n--- Multivariable Analysis ---")
    X = sm.add_constant(df[predictors])
    model = sm.Logit(y, X).fit(disp=0)
    print(model.summary2())

    # OR table
    var_names = ["const"] + predictors
    multi_table = logistic_or_table(model, var_names)

    # VIF (exclude intercept)
    vif_df = calculate_vif(df[predictors])
    print(f"\n--- VIF ---")
    print(vif_df.to_string(index=False))
    high_vif = vif_df[vif_df["VIF"] > config["vif_threshold"]]
    if len(high_vif) > 0:
        print(f"⚠ WARNING: Variables with VIF > {config['vif_threshold']}: "
              f"{', '.join(high_vif['Variable'])}")

    # C-statistic
    y_pred = model.predict(X)
    c_stat = roc_auc_score(y, y_pred)
    # Bootstrap CI for C-statistic
    n_boot = 1000
    c_boots = []
    for i in range(n_boot):
        rng = np.random.RandomState(i)
        idx = rng.choice(len(y), len(y), replace=True)
        try:
            c_boots.append(roc_auc_score(y.iloc[idx], y_pred.iloc[idx]))
        except ValueError:
            continue
    c_ci = np.percentile(c_boots, [2.5, 97.5])
    print(f"\nC-statistic (AUC) = {c_stat:.3f} (95% CI: {c_ci[0]:.3f}-{c_ci[1]:.3f})")

    # Hosmer-Lemeshow
    hl_chi2, hl_df, hl_p = hosmer_lemeshow_test(y, y_pred)
    print(f"Hosmer-Lemeshow: chi2 = {hl_chi2:.2f}, df = {hl_df}, P = {hl_p:.3f}")

    # Brier score
    brier = brier_score_loss(y, y_pred)
    print(f"Brier score = {brier:.4f}")

    # Merge univariable + multivariable
    if config["run_univariable"]:
        combined = uni_df.merge(multi_table[multi_table["Variable"] != "const"],
                                on="Variable", how="outer")
        combined.to_csv(os.path.join(output_dir, "logistic_regression_table.csv"), index=False)
    else:
        multi_table.to_csv(os.path.join(output_dir, "logistic_regression_table.csv"), index=False)

    # Forest plot
    plot_forest_or(multi_table, output_dir)

    # Results text
    print(f"\n--- Manuscript Text ---")
    print(f"The logistic regression model demonstrated a C-statistic of {c_stat:.3f} "
          f"(95% CI, {c_ci[0]:.3f}-{c_ci[1]:.3f}) and adequate calibration "
          f"(Hosmer-Lemeshow P = {hl_p:.2f}).")

    return model


def run_linear(df, config):
    """Run multiple linear regression analysis."""
    outcome = config["outcome"]
    predictors = config["predictors"]
    output_dir = config["output_dir"]

    y = df[outcome]
    n_total = len(y)

    print(f"\n{'='*60}")
    print(f"MULTIPLE LINEAR REGRESSION")
    print(f"{'='*60}")
    print(f"Outcome: {outcome}")
    print(f"N = {n_total}")
    print(f"Predictors: {len(predictors)}")
    print(f"N per predictor: {n_total / len(predictors):.0f} (recommended >= 10-20)")

    # --- Model fitting ---
    X = sm.add_constant(df[predictors])
    model = sm.OLS(y, X).fit()
    print(model.summary2())

    # Coefficient table
    var_names = ["const"] + predictors
    coef_table = linear_coef_table(model, var_names)
    coef_table["R_squared"] = ""
    coef_table.loc[0, "R_squared"] = f"R²={model.rsquared:.3f}, Adj.R²={model.rsquared_adj:.3f}"
    coef_table.to_csv(os.path.join(output_dir, "linear_regression_table.csv"), index=False)
    print(f"\nR² = {model.rsquared:.3f}")
    print(f"Adjusted R² = {model.rsquared_adj:.3f}")

    # VIF
    vif_df = calculate_vif(df[predictors])
    print(f"\n--- VIF ---")
    print(vif_df.to_string(index=False))
    high_vif = vif_df[vif_df["VIF"] > config["vif_threshold"]]
    if len(high_vif) > 0:
        print(f"⚠ WARNING: Variables with VIF > {config['vif_threshold']}: "
              f"{', '.join(high_vif['Variable'])}")

    # Diagnostic plots
    plot_diagnostic_4panel(model, outcome, output_dir)

    # Normality of residuals
    if n_total < 50:
        stat, p = stats.shapiro(model.resid)
        print(f"\nShapiro-Wilk test on residuals: W = {stat:.4f}, P = {p:.3f}")
    else:
        stat, p = stats.kstest(model.resid, "norm",
                               args=(model.resid.mean(), model.resid.std()))
        print(f"\nKolmogorov-Smirnov test on residuals: D = {stat:.4f}, P = {p:.3f}")

    # Results text
    print(f"\n--- Manuscript Text ---")
    print(f"Multiple linear regression was performed with {outcome} as the dependent variable. "
          f"The model explained {model.rsquared_adj*100:.1f}% of the variance "
          f"(adjusted R² = {model.rsquared_adj:.2f}).")

    return model


# === ENTRY POINT ===

if __name__ == "__main__":
    # Load data
    df = pd.read_csv(CONFIG["data_path"])
    print(f"Data loaded: {df.shape[0]} rows x {df.shape[1]} columns")

    # Encode categorical variables if needed
    for cat_var in CONFIG.get("categorical_vars", []):
        if cat_var in df.columns and df[cat_var].dtype == "object":
            df[cat_var] = pd.Categorical(df[cat_var]).codes

    # Drop rows with missing values in analysis variables
    analysis_vars = [CONFIG["outcome"]] + CONFIG["predictors"]
    n_before = len(df)
    df_complete = df[analysis_vars].dropna()
    n_after = len(df_complete)
    if n_before != n_after:
        print(f"Missing data: {n_before - n_after} rows excluded ({100*(n_before-n_after)/n_before:.1f}%)")
        if (n_before - n_after) / n_before > 0.05:
            print("⚠ Consider multiple imputation (> 5% missing). "
                  "See analysis_guides/missing_data.md")

    # Run appropriate regression
    if CONFIG["regression_type"] == "logistic":
        model = run_logistic(df_complete, CONFIG)
    elif CONFIG["regression_type"] == "linear":
        model = run_linear(df_complete, CONFIG)
    else:
        print(f"Error: Unknown regression type '{CONFIG['regression_type']}'")
        sys.exit(1)

    print(f"\n{'='*60}")
    print("Analysis complete.")

"""
Template: Repeated Measures / Mixed Models Analysis
Supports RM ANOVA, Linear Mixed Models (LMM), and GEE for longitudinal data.
Generates spaghetti plots, model summaries, and Time x Group interaction results.

Usage:
    Modify the CONFIGURATION section below, then run:
        python repeated_measures.py

Input:  CSV in wide or long format with repeated measurements
Output: model summary, spaghetti plot PDF/PNG, group mean trajectory plot
"""

# === REPRODUCIBILITY HEADER ===
import sys
import os
import datetime
import numpy as np
import pandas as pd
from scipy import stats

np.random.seed(42)
print(f"Date: {datetime.date.today()}")
print(f"Python: {sys.version}")
print(f"numpy: {np.__version__}, pandas: {pd.__version__}, scipy: {stats.scipy.__version__}")

try:
    import statsmodels.api as sm
    import statsmodels.formula.api as smf
    print(f"statsmodels: {sm.__version__}")
except ImportError:
    print("Error: statsmodels not installed. Install with: pip install statsmodels")
    sys.exit(1)

try:
    import pingouin as pg
    print(f"pingouin: {pg.__version__}")
except ImportError:
    print("Warning: pingouin not installed. RM ANOVA unavailable. Install: pip install pingouin")
    pg = None

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

STYLE_PATH = os.path.join(os.path.dirname(os.path.dirname(__file__)), "style", "figure_style.mplstyle")
if os.path.exists(STYLE_PATH):
    plt.style.use(STYLE_PATH)


# === CONFIGURATION ===
CONFIG = {
    # Data
    "data_path": "data.csv",
    "output_dir": ".",
    "data_format": "wide",  # "wide" or "long"

    # Variables (for wide format)
    "id_col": "subject_id",
    "group_col": "group",                # between-subject factor (optional, None if single group)
    "time_columns": ["t0", "t1", "t2", "t3"],  # column names for each time point
    "outcome_name": "score",             # name for the outcome in long format

    # Variables (for long format) -- used if data_format == "long"
    "time_col": "time",
    "outcome_col": "score",

    # Analysis method: "rm_anova", "lmm", "gee"
    "method": "lmm",

    # LMM options
    "random_effects": "intercept",  # "intercept" or "intercept_slope"

    # Covariates (optional, for LMM/GEE only)
    "covariates": [],
}


# === HELPER FUNCTIONS ===

def wide_to_long(df, id_col, group_col, time_columns, outcome_name):
    """Convert wide format to long format."""
    if group_col and group_col in df.columns:
        id_vars = [id_col, group_col]
    else:
        id_vars = [id_col]

    df_long = df.melt(
        id_vars=id_vars,
        value_vars=time_columns,
        var_name="time_label",
        value_name=outcome_name
    )

    # Create numeric time variable
    time_map = {col: i for i, col in enumerate(time_columns)}
    df_long["time"] = df_long["time_label"].map(time_map)

    return df_long


def plot_spaghetti(df_long, id_col, time_col, outcome_col, group_col, output_dir):
    """Generate spaghetti plot showing individual trajectories."""
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # Left: individual trajectories
    ax = axes[0]
    if group_col and group_col in df_long.columns:
        groups = df_long[group_col].unique()
        colors = plt.cm.tab10(np.linspace(0, 1, len(groups)))
        for grp, color in zip(groups, colors):
            grp_data = df_long[df_long[group_col] == grp]
            for subj_id in grp_data[id_col].unique():
                subj = grp_data[grp_data[id_col] == subj_id]
                ax.plot(subj[time_col], subj[outcome_col], alpha=0.2, color=color, linewidth=0.5)
            # Add group label
            ax.plot([], [], color=color, alpha=0.5, label=f"{grp}")
        ax.legend(title=group_col)
    else:
        for subj_id in df_long[id_col].unique():
            subj = df_long[df_long[id_col] == subj_id]
            ax.plot(subj[time_col], subj[outcome_col], alpha=0.15, color="steelblue", linewidth=0.5)

    ax.set_xlabel("Time")
    ax.set_ylabel(outcome_col)
    ax.set_title("Individual Trajectories")

    # Right: group means with error bars
    ax = axes[1]
    if group_col and group_col in df_long.columns:
        for grp, color in zip(groups, colors):
            grp_data = df_long[df_long[group_col] == grp]
            means = grp_data.groupby(time_col)[outcome_col].mean()
            sems = grp_data.groupby(time_col)[outcome_col].sem()
            ax.errorbar(means.index, means.values, yerr=1.96 * sems.values,
                        marker="o", capsize=4, label=f"{grp}", color=color)
        ax.legend(title=group_col)
    else:
        means = df_long.groupby(time_col)[outcome_col].mean()
        sems = df_long.groupby(time_col)[outcome_col].sem()
        ax.errorbar(means.index, means.values, yerr=1.96 * sems.values,
                    marker="o", capsize=4, color="steelblue")

    ax.set_xlabel("Time")
    ax.set_ylabel(f"{outcome_col} (mean ± 95% CI)")
    ax.set_title("Group Mean Trajectories")

    plt.tight_layout()
    for ext in ["pdf", "png"]:
        fig.savefig(os.path.join(output_dir, f"trajectories.{ext}"),
                    dpi=300, bbox_inches="tight")
    plt.close(fig)
    print(f"Saved: trajectories.pdf/.png")


def run_rm_anova(df_long, id_col, time_col, outcome_col, group_col):
    """Run repeated measures ANOVA with Greenhouse-Geisser correction."""
    if pg is None:
        print("Error: pingouin package required for RM ANOVA. Install: pip install pingouin")
        return None

    print(f"\n{'='*60}")
    print("REPEATED MEASURES ANOVA")
    print(f"{'='*60}")

    # Missing data check
    n_subjects = df_long[id_col].nunique()
    n_timepoints = df_long[time_col].nunique()
    n_expected = n_subjects * n_timepoints
    n_actual = len(df_long.dropna(subset=[outcome_col]))
    n_missing = n_expected - n_actual
    if n_missing > 0:
        pct_missing = 100 * n_missing / n_expected
        print(f"⚠ Missing observations: {n_missing}/{n_expected} ({pct_missing:.1f}%)")
        if pct_missing > 5:
            print("⚠ > 5% missing. RM ANOVA uses complete cases only. Consider LMM instead.")
        # Drop incomplete subjects
        complete_subjects = df_long.groupby(id_col)[outcome_col].count()
        complete_ids = complete_subjects[complete_subjects == n_timepoints].index
        df_long = df_long[df_long[id_col].isin(complete_ids)]
        print(f"  Complete subjects: {len(complete_ids)}/{n_subjects}")

    if group_col and group_col in df_long.columns:
        # Mixed ANOVA (between + within)
        aov = pg.mixed_anova(
            data=df_long, dv=outcome_col, within=time_col,
            between=group_col, subject=id_col, correction=True
        )
    else:
        # One-way RM ANOVA
        aov = pg.rm_anova(
            data=df_long, dv=outcome_col, within=time_col,
            subject=id_col, correction=True
        )

    print(f"\nRM ANOVA Results:")
    print(aov.to_string())

    # Sphericity check (Mauchly's test) via pingouin
    spher = pg.sphericity(
        data=df_long, dv=outcome_col, within=time_col,
        subject=id_col
    )
    if isinstance(spher, tuple):
        spher_result, W, chi2, dof, pval = spher
    else:
        spher_result = spher
        pval = None

    print(f"\nMauchly's sphericity test: {'Passed' if spher_result else 'Violated'}")
    if pval is not None:
        print(f"  P = {pval:.3f}")
    if not spher_result:
        print("  → Greenhouse-Geisser correction applied (see eps column)")

    # Post-hoc pairwise comparisons
    if group_col and group_col in df_long.columns:
        print(f"\n--- Post-hoc: Time within each group ---")
        for grp in df_long[group_col].unique():
            grp_data = df_long[df_long[group_col] == grp]
            pw = pg.pairwise_tests(
                data=grp_data, dv=outcome_col, within=time_col,
                subject=id_col, padjust="bonf"
            )
            print(f"\nGroup: {grp}")
            print(pw[["Contrast", "A", "B", "T", "p-unc", "p-corr"]].to_string(index=False))

    return aov


def run_lmm(df_long, id_col, time_col, outcome_col, group_col, random_effects, covariates):
    """Run Linear Mixed Model."""
    print(f"\n{'='*60}")
    print("LINEAR MIXED MODEL")
    print(f"{'='*60}")

    # Build formula
    fixed_parts = [time_col]
    if group_col and group_col in df_long.columns:
        fixed_parts.append(group_col)
        fixed_parts.append(f"{time_col}:{group_col}")
    fixed_parts.extend(covariates)

    formula = f"{outcome_col} ~ " + " + ".join(fixed_parts)
    print(f"Formula: {formula}")

    # Random effects
    if random_effects == "intercept_slope":
        re_formula = f"~{time_col}"
        print(f"Random effects: intercept + slope")
    else:
        re_formula = "~1"
        print(f"Random effects: intercept only")

    # Fit model
    try:
        model = smf.mixedlm(formula, data=df_long, groups=df_long[id_col],
                             re_formula=re_formula)
        result = model.fit(reml=True)
        print(f"\n{result.summary()}")

        # Extract key results
        print(f"\n--- Fixed Effects ---")
        for name, coef in result.fe_params.items():
            ci = result.conf_int().loc[name]
            p = result.pvalues[name]
            print(f"  {name}: β = {coef:.3f} (95% CI: {ci[0]:.3f} to {ci[1]:.3f}), P = {p:.3f}")

        # Random effects variance
        print(f"\n--- Random Effects ---")
        print(f"  Group variance: {result.cov_re.iloc[0, 0]:.4f}")
        print(f"  Residual variance: {result.scale:.4f}")
        print(f"  ICC (approx): {result.cov_re.iloc[0, 0] / (result.cov_re.iloc[0, 0] + result.scale):.3f}")

        # Model fit
        print(f"\nAIC: {result.aic:.1f}")
        print(f"BIC: {result.bic:.1f}")
        print(f"Log-likelihood: {result.llf:.1f}")

        # Save results
        fe_table = pd.DataFrame({
            "Variable": result.fe_params.index,
            "Coefficient": result.fe_params.values,
            "CI_lower": result.conf_int()[0].values,
            "CI_upper": result.conf_int()[1].values,
            "P_value": result.pvalues.values
        })

        return result, fe_table

    except Exception as e:
        if random_effects == "intercept_slope":
            print(f"⚠ Model with random slope failed to converge: {e}")
            print("  Falling back to random intercept only...")
            return run_lmm(df_long, id_col, time_col, outcome_col,
                           group_col, "intercept", covariates)
        else:
            print(f"Error: LMM failed: {e}")
            return None, None


def run_gee(df_long, id_col, time_col, outcome_col, group_col, covariates):
    """Run Generalized Estimating Equations."""
    from statsmodels.genmod.generalized_estimating_equations import GEE
    from statsmodels.genmod.cov_struct import Exchangeable

    print(f"\n{'='*60}")
    print("GEE (Generalized Estimating Equations)")
    print(f"{'='*60}")

    # Build formula
    fixed_parts = [time_col]
    if group_col and group_col in df_long.columns:
        fixed_parts.append(group_col)
        fixed_parts.append(f"{time_col}:{group_col}")
    fixed_parts.extend(covariates)

    formula = f"{outcome_col} ~ " + " + ".join(fixed_parts)
    print(f"Formula: {formula}")

    # Sort by id and time for proper correlation structure
    df_sorted = df_long.sort_values([id_col, time_col]).reset_index(drop=True)

    # Check cluster size
    n_clusters = df_sorted[id_col].nunique()
    print(f"Number of clusters (subjects): {n_clusters}")
    if n_clusters < 30:
        print("⚠ WARNING: GEE requires >= 30-40 clusters for reliable sandwich SE.")

    # Fit with exchangeable correlation
    try:
        model = GEE.from_formula(
            formula, groups=id_col, data=df_sorted,
            cov_struct=Exchangeable(), family=sm.families.Gaussian()
        )
        result = model.fit()
        print(f"\n{result.summary()}")

        print(f"\n--- Population-averaged Effects ---")
        for name in result.params.index:
            coef = result.params[name]
            ci = result.conf_int().loc[name]
            p = result.pvalues[name]
            print(f"  {name}: β = {coef:.3f} (95% CI: {ci[0]:.3f} to {ci[1]:.3f}), P = {p:.3f}")

        return result

    except Exception as e:
        print(f"Error: GEE failed: {e}")
        return None


# === MAIN ANALYSIS ===

def main():
    config = CONFIG
    output_dir = config["output_dir"]

    # Load data
    df = pd.read_csv(config["data_path"])
    print(f"Data loaded: {df.shape[0]} rows x {df.shape[1]} columns")

    id_col = config["id_col"]
    group_col = config.get("group_col")
    outcome_col = config.get("outcome_col", config.get("outcome_name", "score"))
    time_col = config.get("time_col", "time")

    # Convert to long format if needed
    if config["data_format"] == "wide":
        print("Converting wide → long format...")
        df_long = wide_to_long(
            df, id_col, group_col,
            config["time_columns"], config["outcome_name"]
        )
        outcome_col = config["outcome_name"]
        time_col = "time"
        print(f"Long format: {df_long.shape[0]} rows")
    else:
        df_long = df.copy()

    # Summary
    n_subjects = df_long[id_col].nunique()
    n_timepoints = df_long[time_col].nunique()
    print(f"\nSubjects: {n_subjects}")
    print(f"Time points: {n_timepoints}")
    if group_col and group_col in df_long.columns:
        print(f"Groups: {df_long[group_col].value_counts().to_dict()}")

    # Missing data report
    n_expected = n_subjects * n_timepoints
    n_observed = df_long[outcome_col].notna().sum()
    n_missing = n_expected - n_observed
    print(f"Observations: {n_observed}/{n_expected} ({100*n_missing/n_expected:.1f}% missing)")

    # Descriptive statistics per time point
    print(f"\n--- Descriptive Statistics by Time ---")
    desc = df_long.groupby(time_col)[outcome_col].agg(["count", "mean", "std"])
    print(desc.to_string())

    if group_col and group_col in df_long.columns:
        print(f"\n--- By Group and Time ---")
        desc_grp = df_long.groupby([group_col, time_col])[outcome_col].agg(["count", "mean", "std"])
        print(desc_grp.to_string())

    # Spaghetti plot
    plot_spaghetti(df_long, id_col, time_col, outcome_col, group_col, output_dir)

    # Run analysis
    method = config["method"]
    if method == "rm_anova":
        result = run_rm_anova(df_long, id_col, time_col, outcome_col, group_col)
    elif method == "lmm":
        result, fe_table = run_lmm(
            df_long, id_col, time_col, outcome_col, group_col,
            config["random_effects"], config.get("covariates", [])
        )
        if fe_table is not None:
            fe_table.to_csv(os.path.join(output_dir, "lmm_results.csv"), index=False)
            print(f"Saved: lmm_results.csv")
    elif method == "gee":
        result = run_gee(
            df_long, id_col, time_col, outcome_col, group_col,
            config.get("covariates", [])
        )
    else:
        print(f"Error: Unknown method '{method}'")
        sys.exit(1)

    print(f"\n{'='*60}")
    print("Repeated measures analysis complete.")


if __name__ == "__main__":
    main()

"""
Template: Survey-Weighted Analysis for National Health Surveys
Supports KNHANES, NHANES, KCHS and similar complex survey data.
Produces weighted descriptives, wOR tables, and subgroup analyses.

NOTE: For publication-quality survey analysis, R (survey package) is strongly
recommended. This Python template handles basic weighted analysis but cannot
fully account for strata/cluster in variance estimation. Use the companion
R code blocks in analysis_guides/survey_weighted.md for complex designs.

Usage:
    Modify the CONFIGURATION section below, then run:
        python survey_weighted_analysis.py

Input:  CSV with survey design variables (weight, strata, cluster) + analysis variables
Output: Weighted Table 1, wOR table, subgroup results
"""

# === REPRODUCIBILITY HEADER ===
import sys
import os
import datetime
import numpy as np
import pandas as pd
from scipy import stats

np.random.seed(42)
print(f"Date: {datetime.date.today()}")
print(f"Python: {sys.version}")
print(f"numpy: {np.__version__}, pandas: {pd.__version__}, scipy: {stats.scipy.__version__}")

try:
    import statsmodels.api as sm
    print(f"statsmodels: {sm.__version__}")
except ImportError:
    print("Error: statsmodels not installed. Install with: pip install statsmodels")
    sys.exit(1)

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

STYLE_PATH = os.path.join(os.path.dirname(os.path.dirname(__file__)), "style", "figure_style.mplstyle")
if os.path.exists(STYLE_PATH):
    plt.style.use(STYLE_PATH)


# === CONFIGURATION ===
CONFIG = {
    # Data
    "data_path": "data.csv",
    "output_dir": ".",

    # Survey design variables
    "weight": "wt_itvex",       # Sampling weight column
    "strata": "kstrata",        # Stratification variable (for R code generation)
    "cluster": "psu",           # Cluster/PSU variable (for R code generation)
    "dataset_name": "KNHANES",  # "KNHANES", "NHANES", "KCHS"

    # Analysis variables
    "outcome": "diabetes",           # Binary outcome (0/1)
    "exposure": "depression",        # Primary exposure (binary or categorical)
    "covariates_model1": ["age", "sex"],
    "covariates_model2": ["age", "sex", "income", "education",
                          "smoking", "alcohol", "bmi"],
    "categorical_vars": ["sex", "income", "education", "smoking"],

    # Subgroup stratification variables
    "subgroup_vars": ["sex", "age_group", "income", "obesity"],

    # Output
    "effect_measure": "wOR",  # "wOR" for logistic, "beta" for linear
}


# === HELPER FUNCTIONS ===

def weighted_mean(x, w):
    """Calculate weighted mean."""
    return np.average(x, weights=w)


def weighted_std(x, w):
    """Calculate weighted standard deviation."""
    wm = weighted_mean(x, w)
    return np.sqrt(np.average((x - wm) ** 2, weights=w))


def weighted_proportion(x, w, level=1):
    """Calculate weighted proportion for a binary/categorical variable."""
    mask = x == level
    return np.average(mask, weights=w)


def weighted_smd(x, treatment, weights, is_binary=False):
    """Calculate weighted standardized mean difference."""
    t_mask = treatment == 1
    c_mask = treatment == 0

    w1, w0 = weights[t_mask], weights[c_mask]
    x1, x0 = x[t_mask], x[c_mask]

    wm1 = np.average(x1, weights=w1)
    wm0 = np.average(x0, weights=w0)

    if is_binary:
        denom = np.sqrt((wm1 * (1 - wm1) + wm0 * (1 - wm0)) / 2)
    else:
        wv1 = np.average((x1 - wm1) ** 2, weights=w1)
        wv0 = np.average((x0 - wm0) ** 2, weights=w0)
        denom = np.sqrt((wv1 + wv0) / 2)

    return (wm1 - wm0) / denom if denom > 0 else 0.0


def weighted_table1(df, group_col, continuous_vars, categorical_vars, weight_col):
    """Generate weighted Table 1 with group comparison."""
    groups = sorted(df[group_col].unique())
    results = []

    for var in continuous_vars:
        row = {"Variable": var, "Type": "continuous"}
        for g in groups:
            mask = df[group_col] == g
            wm = weighted_mean(df.loc[mask, var].values, df.loc[mask, weight_col].values)
            ws = weighted_std(df.loc[mask, var].values, df.loc[mask, weight_col].values)
            row[f"Group_{g}"] = f"{wm:.1f} ({ws:.1f})"
        # Weighted SMD
        smd = weighted_smd(
            df[var].values, df[group_col].values,
            df[weight_col].values, is_binary=False
        )
        row["SMD"] = f"{abs(smd):.3f}"
        results.append(row)

    for var in categorical_vars:
        levels = sorted(df[var].unique())
        for level in levels:
            row = {"Variable": f"  {var} = {level}", "Type": "categorical"}
            binary = (df[var] == level).astype(int)
            for g in groups:
                mask = df[group_col] == g
                wp = weighted_proportion(binary[mask].values, df.loc[mask, weight_col].values)
                n = mask.sum()
                row[f"Group_{g}"] = f"{int(n * wp)} ({wp*100:.1f}%)"
            smd = weighted_smd(
                binary.values, df[group_col].values,
                df[weight_col].values, is_binary=True
            )
            row["SMD"] = f"{abs(smd):.3f}"
            results.append(row)

    return pd.DataFrame(results)


def weighted_logistic(df, outcome_col, exposure_col, covariates, weight_col):
    """Run weighted logistic regression and return wOR with 95% CI.

    NOTE: This uses frequency weights which approximate survey weights for
    point estimates but do NOT correctly estimate variance for complex designs.
    For publication, generate and run the R code from survey_weighted.md.
    """
    formula_vars = [exposure_col] + covariates
    X = pd.get_dummies(df[formula_vars], drop_first=True, dtype=float)
    X = sm.add_constant(X)
    y = df[outcome_col]
    w = df[weight_col]

    model = sm.GLM(y, X, family=sm.families.Binomial(), freq_weights=w)
    result = model.fit()

    # Extract exposure effect
    # Find the exposure column(s) in the dummy-encoded X
    exp_cols = [c for c in X.columns if c.startswith(exposure_col)]
    if not exp_cols:
        exp_cols = [exposure_col]

    output_rows = []
    for col in exp_cols:
        coef = result.params[col]
        se = result.bse[col]
        p = result.pvalues[col]
        wor = np.exp(coef)
        ci_lo = np.exp(coef - 1.96 * se)
        ci_hi = np.exp(coef + 1.96 * se)
        output_rows.append({
            "Variable": col,
            "wOR": wor,
            "CI_lower": ci_lo,
            "CI_upper": ci_hi,
            "P": p,
            "formatted": f"{wor:.2f} ({ci_lo:.2f}-{ci_hi:.2f})",
        })

    return result, pd.DataFrame(output_rows)


def generate_r_code(config):
    """Generate publication-ready R code for the same analysis."""
    dataset = config["dataset_name"]
    weight = config["weight"]
    strata = config["strata"]
    cluster = config["cluster"]
    outcome = config["outcome"]
    exposure = config["exposure"]
    covs_m1 = " + ".join(config["covariates_model1"])
    covs_m2 = " + ".join(config["covariates_model2"])
    subgroups = config["subgroup_vars"]

    r_code = f"""# === R Code: Survey-Weighted Analysis ({dataset}) ===
# Requires: survey, tableone
# install.packages(c("survey", "tableone"))

library(survey)
library(tableone)

df <- read.csv("{config['data_path']}")

# Step 1: Declare survey design
design <- svydesign(
  id = ~{cluster},
  strata = ~{strata},
  weights = ~{weight},
  data = df,
  nest = TRUE
)

# Step 2: Weighted Table 1
tab1 <- svyCreateTableOne(
  vars = c({', '.join([f'"{v}"' for v in config['covariates_model2']])}),
  strata = "{exposure}",
  data = design,
  test = TRUE,
  smd = TRUE
)
print(tab1, smd = TRUE)

# Step 3: Model 1 (age + sex)
model1 <- svyglm(
  {outcome} ~ {exposure} + {covs_m1},
  design = design,
  family = quasibinomial()
)
exp(cbind(wOR = coef(model1), confint(model1)))

# Step 4: Model 2 (full adjustment)
model2 <- svyglm(
  {outcome} ~ {exposure} + {covs_m2},
  design = design,
  family = quasibinomial()
)
exp(cbind(wOR = coef(model2), confint(model2)))

# Step 5: Subgroup analyses
"""
    for sg in subgroups:
        r_code += f"""
# Subgroup: {sg}
for (level in unique(df${sg})) {{
  sub_design <- subset(design, {sg} == level)
  sub_model <- svyglm(
    {outcome} ~ {exposure} + {' + '.join([v for v in config['covariates_model2'] if v != sg])},
    design = sub_design,
    family = quasibinomial()
  )
  cat("\\n{sg} =", level, "\\n")
  print(exp(cbind(wOR = coef(sub_model), confint(sub_model)))["{exposure}", ])
}}
"""
    return r_code


# === MAIN ANALYSIS ===

def main():
    config = CONFIG
    df = pd.read_csv(config["data_path"])
    output_dir = config["output_dir"]
    weight_col = config["weight"]

    print(f"Data loaded: {df.shape[0]} rows x {df.shape[1]} columns")
    print(f"Dataset: {config['dataset_name']}")
    print(f"Weight column: {weight_col}")

    # Check weight column exists
    if weight_col not in df.columns:
        print(f"ERROR: Weight column '{weight_col}' not found in data.")
        print(f"Available columns: {list(df.columns)}")
        sys.exit(1)

    # Drop missing
    analysis_vars = ([config["outcome"], config["exposure"], weight_col] +
                     config["covariates_model2"])
    n_before = len(df)
    df = df.dropna(subset=[v for v in analysis_vars if v in df.columns])
    n_after = len(df)
    if n_before != n_after:
        print(f"Excluded {n_before - n_after} rows with missing data "
              f"({100*(n_before-n_after)/n_before:.1f}%)")

    outcome_col = config["outcome"]
    exposure_col = config["exposure"]

    # Weighted sample size
    total_weight = df[weight_col].sum()
    print(f"Unweighted N: {len(df):,}")
    print(f"Weighted N: {total_weight:,.0f}")

    print(f"\n{'='*60}")
    print(f"SURVEY-WEIGHTED ANALYSIS")
    print(f"{'='*60}")

    # --- Weighted Table 1 ---
    print(f"\n--- Weighted Table 1 ---")
    continuous_vars = [v for v in config["covariates_model2"]
                       if v not in config["categorical_vars"]]
    categorical_vars = [v for v in config["covariates_model2"]
                        if v in config["categorical_vars"]]

    tab1 = weighted_table1(
        df, exposure_col, continuous_vars, categorical_vars, weight_col
    )
    print(tab1.to_string(index=False))
    tab1.to_csv(os.path.join(output_dir, "weighted_table1.csv"), index=False)
    print("Saved: weighted_table1.csv")

    # --- Model 1: Minimal adjustment ---
    print(f"\n--- Model 1: Adjusted for {', '.join(config['covariates_model1'])} ---")
    result1, wor1 = weighted_logistic(
        df, outcome_col, exposure_col,
        config["covariates_model1"], weight_col
    )
    print(wor1[["Variable", "formatted", "P"]].to_string(index=False))

    # --- Model 2: Full adjustment ---
    print(f"\n--- Model 2: Adjusted for {', '.join(config['covariates_model2'])} ---")
    result2, wor2 = weighted_logistic(
        df, outcome_col, exposure_col,
        config["covariates_model2"], weight_col
    )
    print(wor2[["Variable", "formatted", "P"]].to_string(index=False))

    # Combine wOR results
    wor_combined = pd.DataFrame({
        "Exposure": wor1["Variable"],
        "Model1_wOR": wor1["formatted"],
        "Model1_P": wor1["P"].map(lambda x: f"{x:.3f}" if x >= 0.001 else "<0.001"),
        "Model2_wOR": wor2["formatted"],
        "Model2_P": wor2["P"].map(lambda x: f"{x:.3f}" if x >= 0.001 else "<0.001"),
    })
    wor_combined.to_csv(os.path.join(output_dir, "wor_results.csv"), index=False)
    print("\nSaved: wor_results.csv")

    # --- Subgroup Analyses ---
    print(f"\n--- Subgroup Analyses ---")
    subgroup_results = []

    for sg_var in config["subgroup_vars"]:
        if sg_var not in df.columns:
            print(f"  Skipping {sg_var} (not in data)")
            continue

        covs_no_sg = [v for v in config["covariates_model2"] if v != sg_var]
        for level in sorted(df[sg_var].unique()):
            subset = df[df[sg_var] == level]
            if len(subset) < 30:
                continue
            try:
                _, wor_sg = weighted_logistic(
                    subset, outcome_col, exposure_col,
                    covs_no_sg, weight_col
                )
                for _, row in wor_sg.iterrows():
                    subgroup_results.append({
                        "Subgroup": sg_var,
                        "Level": level,
                        "wOR": row["formatted"],
                        "P": f"{row['P']:.3f}" if row["P"] >= 0.001 else "<0.001",
                    })
            except Exception as e:
                print(f"  {sg_var}={level}: analysis failed ({e})")

    if subgroup_results:
        sg_df = pd.DataFrame(subgroup_results)
        print(sg_df.to_string(index=False))
        sg_df.to_csv(os.path.join(output_dir, "subgroup_results.csv"), index=False)
        print("Saved: subgroup_results.csv")

    # --- Generate R code ---
    print(f"\n--- R Code (for publication-quality analysis) ---")
    r_code = generate_r_code(config)
    r_path = os.path.join(output_dir, "survey_analysis.R")
    with open(r_path, "w") as f:
        f.write(r_code)
    print(f"Saved: {r_path}")
    print("NOTE: Run the R code for correct variance estimation with strata/cluster.")

    # --- Summary ---
    print(f"\n{'='*60}")
    print("Survey-weighted analysis complete.")
    print(f"  Table 1: weighted_table1.csv")
    print(f"  wOR results: wor_results.csv")
    if subgroup_results:
        print(f"  Subgroup results: subgroup_results.csv")
    print(f"  R code: survey_analysis.R")
    print(f"\n⚠ Python results use frequency weights only.")
    print(f"  For publication, run survey_analysis.R with full design specification.")


if __name__ == "__main__":
    main()

"""
survival_analysis.py — Publication-ready Survival Analysis
===========================================================
Kaplan-Meier curves, log-rank test, Cox proportional hazards model.

Input CSV columns:
    time        : float — Time to event or censoring (in days/months/years)
    event       : int   — Event indicator (1 = event occurred, 0 = censored)
    group       : str   — Group variable for KM comparison (optional)
    [covariates]: float/int — Covariates for Cox model

Outputs:
    - KM plot with at-risk table (PNG/PDF, 300 DPI)
    - Cox model results table (CSV)
    - Console: complete results summary

Usage:
    python survival_analysis.py --input data.csv \
        --time time_months --event event \
        --group treatment_group \
        --covariates age sex bmi comorbidity \
        --time-unit "Months" \
        --output survival_analysis

Dependencies:
    pip install lifelines pandas numpy matplotlib scipy
"""

import argparse
import sys
from datetime import datetime

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

try:
    from lifelines import KaplanMeierFitter, CoxPHFitter
    from lifelines.statistics import logrank_test, multivariate_logrank_test
    from lifelines.utils import restricted_mean_survival_time, median_survival_times
    LIFELINES_AVAILABLE = True
except ImportError:
    LIFELINES_AVAILABLE = False
    print("ERROR: lifelines not installed. Install with: pip install lifelines")
    sys.exit(1)

# ── Reproducibility ───────────────────────────────────────────────────────────
SCRIPT_VERSION = "1.0.0"
SEED = 42
np.random.seed(SEED)


def _fmt_time(x):
    """Format a survival time; 'NR' (not reached) for inf/NaN."""
    if x is None or not np.isfinite(x):
        return "NR"
    return f"{x:.1f}"


def median_with_ci(kmf, unit):
    """Median survival with its 95% CI from a fitted KaplanMeierFitter.

    A median point estimate reported without its CI is incomplete (analyze-stats
    Survival reporting rule). Returns 'NR' when the median is not reached.
    """
    med = kmf.median_survival_time_
    try:
        ci = median_survival_times(kmf.confidence_interval_)
        lo, hi = ci.iloc[0, 0], ci.iloc[0, 1]
    except Exception:
        lo = hi = float("nan")
    return f"{_fmt_time(med)} (95% CI {_fmt_time(lo)}–{_fmt_time(hi)}) {unit}"
print(f"survival_analysis.py v{SCRIPT_VERSION} | {datetime.now().strftime('%Y-%m-%d')}")

# ── Style ─────────────────────────────────────────────────────────────────────
plt.rcParams.update({
    "font.family": "Arial",
    "font.size": 9,
    "axes.labelsize": 9,
    "xtick.labelsize": 8,
    "ytick.labelsize": 8,
})

WONG_COLORS = ["#0072B2", "#D55E00", "#009E73", "#E69F00", "#CC79A7", "#56B4E9"]
LINE_STYLES = ["-", "--", "-.", ":", (0, (3, 1, 1, 1))]


def load_data(filepath, time_col, event_col):
    df = pd.read_csv(filepath)
    for col in [time_col, event_col]:
        if col not in df.columns:
            raise ValueError(f"Column '{col}' not found in {filepath}")
    df = df.dropna(subset=[time_col, event_col])
    return df


def km_analysis(df, time_col, event_col, group_col=None,
                time_unit="Months", output_path="survival"):
    """Kaplan-Meier analysis with at-risk table."""
    fig = plt.figure(figsize=(7.0, 5.0))

    if group_col and group_col in df.columns:
        groups = sorted(df[group_col].dropna().unique())
        n_groups = len(groups)
        gs = gridspec.GridSpec(2, 1, height_ratios=[3.5, 1.5], hspace=0.05)
    else:
        groups = None
        n_groups = 1
        gs = gridspec.GridSpec(2, 1, height_ratios=[3.5, 1.5], hspace=0.05)

    ax_km = fig.add_subplot(gs[0])
    ax_risk = fig.add_subplot(gs[1])

    fitters = {}

    if groups:
        for i, g in enumerate(groups):
            mask = df[group_col] == g
            d = df[mask]
            kmf = KaplanMeierFitter()
            kmf.fit(d[time_col], d[event_col], label=str(g))
            fitters[g] = kmf

            kmf.plot_survival_function(
                ax=ax_km,
                color=WONG_COLORS[i % len(WONG_COLORS)],
                linestyle=LINE_STYLES[i % len(LINE_STYLES)],
                linewidth=1.5,
                ci_show=True,
                ci_alpha=0.15,
            )

            # Median survival with 95% CI
            print(f"  {g}: Median survival = {median_with_ci(kmf, time_unit)}")

        # Log-rank test
        if n_groups == 2:
            g1, g2 = groups
            lr = logrank_test(
                df.loc[df[group_col] == g1, time_col],
                df.loc[df[group_col] == g2, time_col],
                df.loc[df[group_col] == g1, event_col],
                df.loc[df[group_col] == g2, event_col],
            )
            p_str = f"P {'< .001' if lr.p_value < 0.001 else f'= {lr.p_value:.3f}'}"
            ax_km.text(0.98, 0.5, f"Log-rank {p_str}",
                       transform=ax_km.transAxes, fontsize=8,
                       va="center", ha="right")
            print(f"\n  Log-rank test: χ² = {lr.test_statistic:.3f}, P = {lr.p_value:.4f}")
        else:
            mlr = multivariate_logrank_test(
                df[time_col], df[group_col], df[event_col]
            )
            p_str = f"P {'< .001' if mlr.p_value < 0.001 else f'= {mlr.p_value:.3f}'}"
            ax_km.text(0.98, 0.5, f"Log-rank {p_str}",
                       transform=ax_km.transAxes, fontsize=8,
                       va="center", ha="right")
    else:
        # Single-group KM
        kmf = KaplanMeierFitter()
        kmf.fit(df[time_col], df[event_col], label="All patients")
        fitters["All"] = kmf
        kmf.plot_survival_function(ax=ax_km, color=WONG_COLORS[0],
                                    linewidth=1.5, ci_show=True, ci_alpha=0.2)
        print(f"  Median survival: {median_with_ci(kmf, time_unit)}")

    # KM axis formatting
    ax_km.set_ylim(-0.02, 1.05)
    ax_km.set_xlim(left=0)
    ax_km.set_ylabel("Survival probability", fontsize=9)
    ax_km.set_xlabel("")
    ax_km.spines[["top", "right"]].set_visible(False)
    ax_km.legend(fontsize=8, frameon=False)
    ax_km.grid(axis="y", color="#EEEEEE", linewidth=0.5)

    # At-risk table
    ax_risk.set_axis_off()
    time_points = np.linspace(0, df[time_col].max() * 0.95, 6)
    row_offset = 1.0

    for i, (label, kmf) in enumerate(fitters.items()):
        at_risk = [kmf.event_table["at_risk"].iloc[
            max(0, (kmf.event_table.index <= t).sum() - 1)
        ] for t in time_points]
        ax_risk.text(-0.01, row_offset - i * 0.35, str(label),
                     transform=ax_risk.transAxes, fontsize=7.5,
                     va="top", ha="right",
                     color=WONG_COLORS[i % len(WONG_COLORS)])
        for j, (t, n) in enumerate(zip(time_points, at_risk)):
            ax_risk.text(j / (len(time_points) - 1), row_offset - i * 0.35,
                          str(int(n)), transform=ax_risk.transAxes,
                          fontsize=7.5, va="top", ha="center")

    ax_km.set_xticks(time_points)
    ax_km.set_xlabel(time_unit, fontsize=9)

    plt.suptitle("Kaplan-Meier Survival Curves", fontsize=10, y=1.01, fontweight="bold")
    plt.tight_layout()

    for ext in ["pdf", "png"]:
        outfile = f"{output_path}_km.{ext}"
        dpi = 300 if ext == "png" else None
        plt.savefig(outfile, dpi=dpi, bbox_inches="tight",
                    facecolor="white", edgecolor="none")
        print(f"Saved: {outfile}")
    plt.close()


def cox_analysis(df, time_col, event_col, covariates, output_path="survival",
                 cluster_col=None):
    """Cox proportional hazards model.

    cluster_col: id column for nested observation units (e.g. multiple lesions /
    eyes / repeated episodes per subject). When set, lifelines computes a robust
    (cluster-sandwich) variance so the HR CIs reflect within-subject correlation
    rather than treating correlated rows as independent.
    """
    keep = [time_col, event_col] + covariates + ([cluster_col] if cluster_col else [])
    model_df = df[keep].dropna()
    print(f"\n── Cox PH Model (N = {len(model_df)}) ───────────────────────")

    # Events-per-variable (EPV) gate — mirror of the logistic EPV rule. A Wald CI
    # from a sparse-event model is not stable; warn and rely on the penalizer.
    n_events = int(model_df[event_col].sum())
    epv = n_events / max(len(covariates), 1)
    print(f"EPV = {epv:.1f} ({n_events} events / {len(covariates)} covariates; "
          f"minimum recommended: 10)")
    if epv < 10:
        print("⚠ WARNING: EPV < 10 — Cox estimates may be unstable. The penalized "
              "fit (penalizer=0.1) shrinks coefficients; consider Firth/penalized "
              "Cox or profile-likelihood CIs and interpret Wald CIs with caution.")

    cph = CoxPHFitter(penalizer=0.1)
    fit_kw = {"duration_col": time_col, "event_col": event_col}
    if cluster_col:
        fit_kw["cluster_col"] = cluster_col  # robust (cluster-sandwich) SE
        print(f"Robust cluster-sandwich SE on '{cluster_col}' (nested units).")
    cph.fit(model_df, **fit_kw)
    cph.print_summary()

    # Check PH assumption
    print("\n── Proportional Hazards Assumption (Schoenfeld residuals) ───")
    try:
        cph.check_assumptions(model_df, p_value_threshold=0.05, show_plots=False)
    except Exception as e:
        print(f"  PH assumption check: {e}")

    # Save results table
    summary = cph.summary.copy()
    summary.columns = [c.replace(" ", "_") for c in summary.columns]
    outfile = f"{output_path}_cox_results.csv"
    summary.to_csv(outfile)
    print(f"\nSaved: {outfile}")

    # Concordance index
    print(f"\nConcordance index (C-statistic): {cph.concordance_index_:.3f}")

    return cph


def rmst_analysis(df, time_col, event_col, group_col, t_star, output_path="survival"):
    """Restricted Mean Survival Time analysis."""
    if group_col not in df.columns:
        return

    groups = df[group_col].dropna().unique()
    print(f"\n── Restricted Mean Survival Time (t* = {t_star}) ────────────")
    for g in groups:
        d = df[df[group_col] == g]
        kmf = KaplanMeierFitter()
        kmf.fit(d[time_col], d[event_col])
        rmst = restricted_mean_survival_time(kmf, t=t_star)
        print(f"  {g}: RMST = {rmst:.2f}")


def main():
    parser = argparse.ArgumentParser(description="Survival analysis")
    parser.add_argument("--input", required=True)
    parser.add_argument("--time", required=True, help="Time column name")
    parser.add_argument("--event", required=True, help="Event column name (1=event, 0=censored)")
    parser.add_argument("--group", default=None)
    parser.add_argument("--covariates", nargs="+", default=None)
    parser.add_argument("--cluster", default=None,
                        help="ID column for nested units → robust cluster SE in Cox")
    parser.add_argument("--time-unit", default="Months")
    parser.add_argument("--rmst-t", type=float, default=None,
                        help="t* for RMST (in same units as time)")
    parser.add_argument("--output", default="survival_analysis")

    args = parser.parse_args()

    df = load_data(args.input, args.time, args.event)
    print(f"\nN = {len(df)} | Events: {df[args.event].sum()} "
          f"({df[args.event].mean()*100:.1f}%)")

    print("\n── Kaplan-Meier Analysis ────────────────────")
    km_analysis(df, args.time, args.event, args.group,
                args.time_unit, args.output)

    if args.covariates:
        cox_analysis(df, args.time, args.event, args.covariates, args.output,
                     cluster_col=args.cluster)

    if args.rmst_t and args.group:
        rmst_analysis(df, args.time, args.event, args.group,
                      args.rmst_t, args.output)

    print(f"\n── Session Info ─────────────────────────────")
    print(f"Python: {sys.version.split()[0]}")
    for pkg in ["lifelines", "pandas", "numpy", "scipy"]:
        try:
            import importlib
            m = importlib.import_module(pkg)
            print(f"  {pkg}: {m.__version__}")
        except Exception:
            pass


if __name__ == "__main__":
    if len(sys.argv) == 1:
        print("Usage: python survival_analysis.py --input data.csv "
              "--time time_months --event event --group treatment")
    else:
        main()

"""
Template: Table 1 — Baseline Demographics
Generates a publication-ready demographics table from tabular data.

Usage:
    Modify the CONFIGURATION section below, then run:
        python table1_demographics.py

Input:  CSV or Excel file with one row per subject
Output: table1_demographics.csv, table1_demographics.md (console), summary text
"""

# === REPRODUCIBILITY HEADER ===
import sys
import datetime
import numpy as np
import pandas as pd
from scipy import stats

np.random.seed(42)
print(f"Date: {datetime.date.today()}")
print(f"Python: {sys.version}")
print(f"numpy: {np.__version__}, pandas: {pd.__version__}, scipy: {stats.scipy.__version__}")
print()

# === CONFIGURATION (modify for your study) ===
INPUT_FILE = "data.csv"           # Path to input data
OUTPUT_DIR = "."                  # Output directory
GROUP_COL = None                  # Column name for group comparison (None = no comparison)
CONTINUOUS_VARS = []              # List of continuous variable column names
CATEGORICAL_VARS = []            # List of categorical variable column names
VAR_LABELS = {}                   # Optional display labels: {"col_name": "Display Name (unit)"}
DECIMAL_PLACES = 1                # Decimal places for continuous variables
# ==============================================


def load_data(filepath: str) -> pd.DataFrame:
    """Load CSV or Excel file."""
    if filepath.endswith((".xlsx", ".xls")):
        return pd.read_excel(filepath)
    return pd.read_csv(filepath)


def test_normality(series: pd.Series, alpha: float = 0.05) -> tuple:
    """Test normality using Shapiro-Wilk (n<50) or Kolmogorov-Smirnov (n>=50)."""
    clean = series.dropna()
    if len(clean) < 3:
        return False, np.nan
    if len(clean) < 50:
        stat, p = stats.shapiro(clean)
    else:
        stat, p = stats.kstest(clean, "norm", args=(clean.mean(), clean.std()))
    return p >= alpha, p


def format_continuous(series: pd.Series, is_normal: bool, dp: int = 1) -> str:
    """Format continuous variable as mean +/- SD or median (IQR)."""
    clean = series.dropna()
    if is_normal:
        return f"{clean.mean():.{dp}f} +/- {clean.std():.{dp}f}"
    else:
        q1, median, q3 = clean.quantile([0.25, 0.5, 0.75])
        return f"{median:.{dp}f} ({q1:.{dp}f}-{q3:.{dp}f})"


def format_categorical(series: pd.Series) -> dict:
    """Format categorical variable as n (%)."""
    counts = series.value_counts(dropna=False)
    total = len(series)
    result = {}
    for cat, n in counts.items():
        pct = 100.0 * n / total if total > 0 else 0.0
        label = str(cat) if pd.notna(cat) else "Missing"
        result[label] = f"{n} ({pct:.1f}%)"
    return result


def compare_continuous(groups: list, is_normal: bool) -> tuple:
    """Compare continuous variable between groups."""
    clean_groups = [g.dropna() for g in groups]
    if len(clean_groups) == 2:
        if is_normal:
            stat, p = stats.ttest_ind(*clean_groups)
            return "t-test", p
        else:
            stat, p = stats.mannwhitneyu(*clean_groups, alternative="two-sided")
            return "Mann-Whitney U", p
    else:
        if is_normal:
            stat, p = stats.f_oneway(*clean_groups)
            return "ANOVA", p
        else:
            stat, p = stats.kruskal(*clean_groups)
            return "Kruskal-Wallis", p


def compare_categorical(contingency: pd.DataFrame) -> tuple:
    """Compare categorical variable between groups using chi-square or Fisher's exact."""
    table = contingency.values
    # Use Fisher's exact test for 2x2 tables with any expected count < 5
    if table.shape == (2, 2):
        expected = stats.chi2_contingency(table, correction=False)[3]
        if (expected < 5).any():
            odds_ratio, p = stats.fisher_exact(table)
            return "Fisher's exact", p
    # Otherwise chi-square (with Yates' correction for 2x2)
    chi2, p, dof, expected = stats.chi2_contingency(table)
    if (expected < 5).any():
        print(f"  Warning: expected counts < 5 in some cells; consider collapsing categories.")
    return "Chi-square", p


def format_p(p: float) -> str:
    """Format p-value per reporting rules."""
    if pd.isna(p):
        return ""
    if p < 0.001:
        return "<0.001"
    return f"{p:.3f}"


def build_table1(df: pd.DataFrame) -> pd.DataFrame:
    """Build the full Table 1 dataframe."""
    rows = []
    test_footnotes = set()

    if GROUP_COL and GROUP_COL in df.columns:
        groups = df[GROUP_COL].dropna().unique()
        groups = sorted(groups, key=str)
        group_dfs = {str(g): df[df[GROUP_COL] == g] for g in groups}
    else:
        groups = []
        group_dfs = {}

    n_row = {"Variable": "n"}
    if group_dfs:
        for g_name, g_df in group_dfs.items():
            n_row[g_name] = str(len(g_df))
    n_row["Overall"] = str(len(df))
    if group_dfs:
        n_row["p-value"] = ""
        n_row["Test"] = ""
    rows.append(n_row)

    # --- Continuous variables ---
    for var in CONTINUOUS_VARS:
        if var not in df.columns:
            print(f"Warning: '{var}' not found in data, skipping.")
            continue
        label = VAR_LABELS.get(var, var)
        is_normal, norm_p = test_normality(df[var])
        stat_type = "mean +/- SD" if is_normal else "median (IQR)"

        row = {"Variable": f"{label}, {stat_type}"}
        row["Overall"] = format_continuous(df[var], is_normal, DECIMAL_PLACES)

        if group_dfs:
            for g_name, g_df in group_dfs.items():
                row[g_name] = format_continuous(g_df[var], is_normal, DECIMAL_PLACES)
            grp_series = [g_df[var] for g_df in group_dfs.values()]
            test_name, p_val = compare_continuous(grp_series, is_normal)
            row["p-value"] = format_p(p_val)
            row["Test"] = test_name
            test_footnotes.add(test_name)

        rows.append(row)

    # --- Categorical variables ---
    for var in CATEGORICAL_VARS:
        if var not in df.columns:
            print(f"Warning: '{var}' not found in data, skipping.")
            continue
        label = VAR_LABELS.get(var, var)

        # Header row for this variable
        header_row = {"Variable": f"{label}, n (%)"}

        if group_dfs:
            contingency = pd.crosstab(df[var], df[GROUP_COL])
            test_name, p_val = compare_categorical(contingency)
            header_row["p-value"] = format_p(p_val)
            header_row["Test"] = test_name
            test_footnotes.add(test_name)
            for g_name in group_dfs:
                header_row[g_name] = ""
        header_row["Overall"] = ""
        rows.append(header_row)

        # Category rows
        all_cats = format_categorical(df[var])
        for cat_label, cat_val in all_cats.items():
            cat_row = {"Variable": f"  {cat_label}"}
            cat_row["Overall"] = cat_val
            if group_dfs:
                for g_name, g_df in group_dfs.items():
                    cat_counts = format_categorical(g_df[var])
                    cat_row[g_name] = cat_counts.get(cat_label, "0 (0.0%)")
                cat_row["p-value"] = ""
                cat_row["Test"] = ""
            rows.append(cat_row)

    # Build dataframe
    table_df = pd.DataFrame(rows)

    # Reorder columns
    col_order = ["Variable"]
    if group_dfs:
        col_order += [str(g) for g in sorted(groups, key=str)]
    col_order += ["Overall"]
    if group_dfs:
        col_order += ["p-value", "Test"]
    table_df = table_df[col_order]

    print("\n--- Test Methods Used ---")
    for t in sorted(test_footnotes):
        print(f"  - {t}")
    print()

    return table_df


def save_outputs(table_df: pd.DataFrame) -> None:
    """Save table as CSV and print markdown version."""
    import os
    csv_path = os.path.join(OUTPUT_DIR, "table1_demographics.csv")

    # Save CSV (without Test column for manuscript)
    export_df = table_df.drop(columns=["Test"], errors="ignore")
    export_df.to_csv(csv_path, index=False)
    print(f"Saved: {csv_path}")

    # Print markdown table
    print("\n--- Table 1. Baseline Characteristics ---\n")
    print(export_df.to_markdown(index=False))


def print_results_text(df: pd.DataFrame, table_df: pd.DataFrame) -> None:
    """Print copy-paste-ready text for Results section."""
    print("\n--- Results Text (copy-paste ready) ---\n")
    n_total = len(df)

    if GROUP_COL and GROUP_COL in df.columns:
        groups = df[GROUP_COL].dropna().unique()
        group_counts = df[GROUP_COL].value_counts()
        parts = [f"{g} (n = {group_counts[g]})" for g in sorted(groups, key=str)]
        print(f"A total of {n_total} subjects were included and divided into "
              f"{len(groups)} groups: {', '.join(parts)}. "
              f"Baseline characteristics are summarized in Table 1.")
    else:
        print(f"A total of {n_total} subjects were included. "
              f"Baseline characteristics are summarized in Table 1.")

    # Report missing data
    missing = df[CONTINUOUS_VARS + CATEGORICAL_VARS].isnull().sum()
    missing = missing[missing > 0]
    if len(missing) > 0:
        print(f"\nMissing data: ", end="")
        parts = [f"{VAR_LABELS.get(col, col)}: {n} ({100*n/len(df):.1f}%)"
                 for col, n in missing.items()]
        print("; ".join(parts) + ".")


# === MAIN ===
if __name__ == "__main__":
    print("=" * 60)
    print("Table 1: Baseline Demographics")
    print("=" * 60)

    df = load_data(INPUT_FILE)
    print(f"\nLoaded: {INPUT_FILE} ({df.shape[0]} rows, {df.shape[1]} columns)")

    # Auto-detect if CONTINUOUS_VARS and CATEGORICAL_VARS are empty
    if not CONTINUOUS_VARS and not CATEGORICAL_VARS:
        print("\nNo variables specified. Auto-detecting from data types...")
        for col in df.columns:
            if col == GROUP_COL:
                continue
            if pd.api.types.is_numeric_dtype(df[col]) and df[col].nunique() > 10:
                CONTINUOUS_VARS.append(col)
                print(f"  Continuous: {col}")
            elif df[col].nunique() <= 20:
                CATEGORICAL_VARS.append(col)
                print(f"  Categorical: {col}")

    table_df = build_table1(df)
    save_outputs(table_df)
    print_results_text(df, table_df)

# Analysis execution record challenge

Run `bash verify.sh` from this directory. Requires Python 3.10+, numpy, pandas,
SciPy, scikit-learn and matplotlib; no network or private data is used.

The bundled demo creates original synthetic observations in a temporary project.
The normal control checks exact confusion counts, two matching calculations and
a current byte audit. The defect control adds one synthetic input row: two audits
must report drift without rewriting the old record. The temporary project is removed.
The complete regression suite under `../../tests/test_analysis_run.py` additionally
covers units, denominators, undefined metrics, rounded-input refusal and output/code
changes. These controls test software behavior, not clinical validity.

#!/usr/bin/env bash
# Standalone original synthetic normal and stale-result controls.
set -euo pipefail
HERE="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
python3 - "$HERE/.." <<'PY'
import json
from pathlib import Path
import subprocess
import sys
import tempfile

scripts = Path(sys.argv[1]).resolve()
with tempfile.TemporaryDirectory() as temp:
    project = Path(temp) / "project"
    demo = subprocess.run([sys.executable, str(scripts / "demo_analysis_run.py"),
                           "--out", str(project)], capture_output=True, text=True, check=True)
    result = json.loads(demo.stdout)
    assert result["synthetic_counts"] == {"TP": 8, "FP": 3, "TN": 7, "FN": 2}
    assert result["comparison"]["status"] == "compared"
    assert result["comparison"]["recorded_numeric_results_equal"]
    runner = [sys.executable, str(scripts / "run_analysis.py")]
    audit = runner + ["audit", "--project-root", str(project), "--out", "runs/first"]
    assert subprocess.run(audit, capture_output=True).returncode == 0
    manifest = project / "runs/first/_analysis_outputs.md"
    original = manifest.read_bytes()
    with (project / "data.csv").open("a") as stream:
        stream.write("synthetic_changed,0,1\n")
    for _ in range(2):
        check = subprocess.run(audit, capture_output=True, text=True)
        assert check.returncode == 1
        assert "input:data" in json.loads(check.stdout)["changes"]
        assert manifest.read_bytes() == original
print("PASS: synthetic calculation, repeat, stale-input detection and read-only audit")
PY

"""
Analysis: synthetic BAD fixture for the generated-code quality gate.
Date: 2020-01-01
Random seed: (intentionally missing)
"""
import numpy as np
import pandas as pd
import json  # unused import (dead dependency)

# absolute path literal + source read
df = pd.read_csv("/Users/researcher/data/cohort.csv")

# hand-typed tabular data instead of read_csv + subset
ref = pd.DataFrame({
    "auc": [0.81, 0.83, 0.79, 0.88, 0.84, 0.77, 0.82, 0.86, 0.80, 0.85, 0.78, 0.87, 0.83, 0.81],
})

# randomness with no seed set -> non-reproducible
boot = np.random.choice(df["auc"].values, size=1000, replace=True)

breakpoint()  # debugger left in

# writes back to the source path -> overwrites raw data
df.to_csv("/Users/researcher/data/cohort.csv")

"""
Analysis: synthetic CLEAN fixture for the generated-code quality gate.
Date: 2026-01-01
Random seed: 42
"""
import numpy as np
import pandas as pd

np.random.seed(42)

# portable relative path; no hand-typed data
df = pd.read_csv("cohort.csv")

# seeded randomness -> reproducible
boot = np.random.choice(df["auc"].values, size=1000, replace=True)
boot_mean = float(np.mean(boot))

# derived output written to a NEW path, not the source
summary = df["auc"].describe()
summary.to_csv("auc_summary.csv")
print(f"bootstrap mean AUC = {boot_mean:.3f}")

"""
Analysis: synthetic CLEAN fixture exercising the colorblind-safe palette.
A hex-color list must NOT be flagged HARDCODED_DATA_LITERAL even alongside a
data-file read (regression for the WONG-palette false positive).
Date: 2026-01-01
Random seed: 42
"""
import numpy as np
import pandas as pd

np.random.seed(42)

# The Wong (2011) colorblind-safe palette that make-figures recommends. Eight
# string literals; the digits live inside the hex codes, not in tabular data.
WONG = ["#000000", "#E69F00", "#56B4E9", "#009E73",
        "#F0E442", "#0072B2", "#D55E00", "#CC79A7"]

# portable relative path; no hand-typed numeric data
df = pd.read_csv("cohort.csv")

means = df.groupby("group")["auc"].mean()
print(WONG[0], float(means.iloc[0]))