"""
Meta-regression Level 2 for predicting ypN0 (axillary pathologic complete
response) after neoadjuvant systemic therapy in breast cancer.

Author: Nivel 2 predictor pipeline.
Inputs:
    data/ypN0_rates_by_subtype.csv
    data/predictor_effect_sizes.csv
    data/regimen_effect.csv

Outputs:
    data/nomogram_coefficients.csv
    nivel2_metaregresion/results_pooled_subtype.csv
    nivel2_metaregresion/results_validation.csv
    nivel2_metaregresion/results_meta_or_subtype.csv

Method:
    - DerSimonian-Laird random-effects pooling of logit-transformed proportions
      for per-subtype ypN0 rates.
    - DL inverse-variance random-effects pooling of log-odds-ratios for
      subtype effects, with reference HR+/HER2-.
    - Multivariable logistic prediction equation built from pooled coefficients,
      complemented with literature priors for predictors with sparse meta-data
      (Ki67, age, cN-stage, regimens, genomic risk).
    - Validation by reconstructing predicted ypN0 in pivotal studies (Z1071,
      RISAS, Montagna, van Amstel KEYNOTE-522).

Anchoring references:
    van Amstel FJG et al. BJS 2024. DOI: 10.1093/bjs/znae203
    Samiei S et al. JAMA Surg 2021. DOI: 10.1001/jamasurg.2021.0891
"""

from __future__ import annotations

import os
import math
from pathlib import Path
from typing import Iterable

import numpy as np
import pandas as pd
from scipy import stats

# -------------------- paths --------------------
HERE = Path(__file__).resolve().parent
ROOT = HERE.parent
DATA = ROOT / "data"
OUT_DATA = DATA
OUT_HERE = HERE


# ======================================================================
# Generic utilities
# ======================================================================
def logit(p: float) -> float:
    """Logit transform with clamp to avoid +/-inf when p in {0,1}."""
    p = max(min(p, 1 - 1e-6), 1e-6)
    return math.log(p / (1 - p))


def inv_logit(x: float) -> float:
    return 1.0 / (1.0 + math.exp(-x))


def logit_var(k: float, n: float) -> float:
    """Variance of logit(k/n) using delta method: 1/k + 1/(n-k)."""
    if k <= 0 or k >= n:
        # continuity correction
        k = max(min(k, n - 0.5), 0.5)
    return 1.0 / k + 1.0 / (n - k)


def wilson_ci(k: int, n: int, conf: float = 0.95) -> tuple[float, float]:
    if n == 0:
        return (float("nan"), float("nan"))
    z = stats.norm.ppf(1 - (1 - conf) / 2)
    phat = k / n
    denom = 1 + z * z / n
    centre = (phat + z * z / (2 * n)) / denom
    half = (z * math.sqrt(phat * (1 - phat) / n + z * z / (4 * n * n))) / denom
    return (max(0.0, centre - half), min(1.0, centre + half))


# ======================================================================
# DerSimonian-Laird random-effects pooling
# ======================================================================
def dersimonian_laird(y: np.ndarray, v: np.ndarray) -> dict:
    """Pool log-effects y with within-study variances v using DL random effects.

    Returns dict with keys: pooled, se, ci_low, ci_high, tau2, q, df, i2, k.
    All on the same scale as the inputs (logit for proportions, log-OR for ORs).
    """
    y = np.asarray(y, dtype=float)
    v = np.asarray(v, dtype=float)
    k = len(y)
    if k == 0:
        return dict(pooled=np.nan, se=np.nan, ci_low=np.nan, ci_high=np.nan,
                    tau2=np.nan, q=np.nan, df=0, i2=np.nan, k=0)
    if k == 1:
        se = math.sqrt(v[0])
        return dict(pooled=float(y[0]), se=se,
                    ci_low=float(y[0] - 1.96 * se),
                    ci_high=float(y[0] + 1.96 * se),
                    tau2=0.0, q=0.0, df=0, i2=0.0, k=1)
    # Fixed-effect weights
    w_fe = 1.0 / v
    mu_fe = float(np.sum(w_fe * y) / np.sum(w_fe))
    q = float(np.sum(w_fe * (y - mu_fe) ** 2))
    df = k - 1
    c = float(np.sum(w_fe) - np.sum(w_fe ** 2) / np.sum(w_fe))
    tau2 = max(0.0, (q - df) / c) if c > 0 else 0.0
    # Random-effects weights
    w_re = 1.0 / (v + tau2)
    mu_re = float(np.sum(w_re * y) / np.sum(w_re))
    se_re = float(math.sqrt(1.0 / np.sum(w_re)))
    i2 = max(0.0, (q - df) / q) * 100.0 if q > 0 else 0.0
    return dict(pooled=mu_re, se=se_re,
                ci_low=mu_re - 1.96 * se_re,
                ci_high=mu_re + 1.96 * se_re,
                tau2=tau2, q=q, df=df, i2=i2, k=k)


# ======================================================================
# Data loading + filtering
# ======================================================================
def load_subtype_rates() -> pd.DataFrame:
    df = pd.read_csv(DATA / "ypN0_rates_by_subtype.csv")
    # Map subtype strings into canonical buckets
    canon = {
        "HR-/HER2+": "HR-/HER2+",
        "HER2+": "HER2+_combined",
        "HR+/HER2+": "HR+/HER2+",
        "TNBC": "TNBC",
        "HR+/HER2-": "HR+/HER2-",
        "Luminal A": "HR+/HER2-",
        "Luminal B": "HR+/HER2-",
    }
    df["subtype_canon"] = df["subtype"].map(canon)
    return df


def load_predictor_effects() -> pd.DataFrame:
    df = pd.read_csv(DATA / "predictor_effect_sizes.csv")
    # Keep only logistic models for the nomogram (drop Cox / KM rows)
    df["model_type_l"] = df["model_type"].fillna("").str.lower()
    return df


def load_regimen() -> pd.DataFrame:
    return pd.read_csv(DATA / "regimen_effect.csv")


# ======================================================================
# Part A - pooled rates by subtype
# ======================================================================
# Approximate per-subtype N from Samiei 2021 pooled (57,531 across 33 studies)
# Distribution typical of NAC trials: ~50% HR+/HER2-, ~25% HER2+, ~20% TNBC.
SAMIEI_TOTAL = 57531
SAMIEI_SUBTYPE_N = {
    "HR+/HER2-": int(SAMIEI_TOTAL * 0.50),
    "HER2+_combined": int(SAMIEI_TOTAL * 0.25),
    "HR-/HER2+": int(SAMIEI_TOTAL * 0.12),
    "HR+/HER2+": int(SAMIEI_TOTAL * 0.13),
    "TNBC": int(SAMIEI_TOTAL * 0.20),
}


def _study_n(row: pd.Series, default: int = 100) -> int:
    """Pick a reasonable N for this row to enable inverse-variance pooling."""
    for col in ("n_subgroup", "n_total"):
        v = row.get(col)
        if pd.notna(v) and v and float(v) > 0:
            return int(v)
    # If row comes from Samiei pooled, use proportional split
    if row.get("study_id", "").startswith("Samiei2021"):
        return SAMIEI_SUBTYPE_N.get(row["subtype_canon"], 5000)
    return default


def pool_subtype_rates(df: pd.DataFrame) -> pd.DataFrame:
    """Pool baseline (chemo +/- HER2 targeting) ypN0 rates per subtype.

    Exclusions for the BASELINE pool:
      - Pembrolizumab arms (Kuemmel/Woodfin/Perry) -> these are modelled as a
        REGIMEN EFFECT, not the baseline rate. Including them would conflate
        rate with treatment.
      - Barron 2018 NCDB rows -> ypN0 *conditional on breast pCR*, not marginal.
      - Sharma NeoSTOP -> combined breast+axillary pCR, not pure ypN0.
      - Montagna 2024 -> ypN0 by inclusion (rate = 1.0 by design).
      - Boughey 2022 Ki67 sub-cells of <50 ages -> they double-count the parent
        315-pt cohort (we retain only the parent rate 0.248).
      - Boughey 2022 genomic sub-cells -> stratified beyond the marginal rate.
      - Lobular subset rates (Mukhtar) -> retained because they don't overlap.
    """
    EXCLUDE = {
        "Kuemmel2026_KEYNOTE522",     # pembro/placebo + chemo
        "Woodfin2023",                # pembro TNBC
        "Perry2025_MSK",              # pembro TNBC
        "Barron2018_NCDB",            # conditional on breast pCR
        "Sharma2021_NeoSTOP",         # combined breast+axillary pCR
        "Montagna2024_JAMAOnc",       # ypN0 by design
    }
    EXCLUDE_SUBTYPE_PREFIX = (
        "HR+/HER2- (age",         # Boughey 2022 Ki67 sub-cells (parent retained)
        "HR+/HER2- (high",        # Boughey 2022 genomic high (parent retained)
        "HR+/HER2- (low",         # Boughey 2022 genomic low (parent retained)
    )
    df = df[~df["study_id"].isin(EXCLUDE)].copy()
    df = df[~df["subtype"].astype(str).str.startswith(EXCLUDE_SUBTYPE_PREFIX)]
    rows = []
    subtype_order = ["HR-/HER2+", "HR+/HER2+", "HER2+_combined",
                     "TNBC", "HR+/HER2-"]
    for s in subtype_order:
        sub = df[df["subtype_canon"] == s].copy()
        if sub.empty:
            continue
        # Drop rows with no rate
        sub = sub[pd.notna(sub["ypN0_rate"])]
        if sub.empty:
            continue
        y = []
        v = []
        ns = []
        for _, r in sub.iterrows():
            p = float(r["ypN0_rate"])
            n_eff = _study_n(r)
            k_eff = max(0.5, p * n_eff)
            y.append(logit(p))
            v.append(logit_var(k_eff, n_eff))
            ns.append(n_eff)
        res = dersimonian_laird(np.array(y), np.array(v))
        rows.append({
            "subtype": s,
            "k_studies": res["k"],
            "n_total": int(sum(ns)),
            "pooled_rate": inv_logit(res["pooled"]),
            "ci_low": inv_logit(res["ci_low"]),
            "ci_high": inv_logit(res["ci_high"]),
            "tau2_logit": res["tau2"],
            "i2_percent": res["i2"],
            "q_stat": res["q"],
        })
    return pd.DataFrame(rows)


# ======================================================================
# Part B - meta-regression of subtype log-ORs
# ======================================================================
def pool_subtype_logors(df: pd.DataFrame) -> pd.DataFrame:
    """Pool per-subtype log-ORs vs HR+/HER2-. Combine direct OR rows from
    predictor_effect_sizes plus rates-derived ORs from ypN0_rates_by_subtype."""
    rates_df = load_subtype_rates()
    # Derive log-OR per study where both ref and target subtype exist
    derived = []
    grouped = rates_df[pd.notna(rates_df["ypN0_rate"])].groupby("study_id")
    for sid, g in grouped:
        sub_map = dict(zip(g["subtype_canon"], g["ypN0_rate"]))
        if "HR+/HER2-" not in sub_map:
            continue
        p_ref = sub_map["HR+/HER2-"]
        if p_ref <= 0 or p_ref >= 1:
            continue
        # Use a stand-in n for variance: study-specific n_subgroup if any, else
        # n_total / # subtypes
        n_default = max(100, int(g["n_total"].max() / max(1, len(g))))
        for target in ("TNBC", "HR-/HER2+", "HR+/HER2+", "HER2+_combined"):
            if target not in sub_map:
                continue
            p_t = sub_map[target]
            if p_t <= 0 or p_t >= 1:
                continue
            row_t = g[g["subtype_canon"] == target].iloc[0]
            row_r = g[g["subtype_canon"] == "HR+/HER2-"].iloc[0]
            n_t = _study_n(row_t, default=n_default)
            n_r = _study_n(row_r, default=n_default)
            k_t = max(0.5, p_t * n_t)
            k_r = max(0.5, p_ref * n_r)
            log_or = math.log((k_t / (n_t - k_t)) / (k_r / (n_r - k_r)))
            v = 1 / k_t + 1 / (n_t - k_t) + 1 / k_r + 1 / (n_r - k_r)
            derived.append({"study": sid, "target": target,
                            "log_or": log_or, "var": v})
    # Add direct ORs reported in predictor_effect_sizes.csv (van Amstel, An2025)
    eff = load_predictor_effects()
    direct = eff[(eff["predictor"] == "subtype") &
                 (eff["model_type_l"].str.contains("logistic"))].copy()
    for _, r in direct.iterrows():
        comp = str(r["comparison"])
        OR = r["OR"]
        if pd.isna(OR) or OR <= 0:
            continue
        # Map An2025 "ER-/HER2+ vs TNBC" by converting to TNBC vs HR+/HER2-
        # via assumption is too speculative; only retain vs HR+/HER2- contrasts.
        if "vs HR+/HER2-" not in comp:
            continue
        target = comp.split(" vs ")[0].strip()
        # Map labels
        rename = {"HR-/HER2+": "HR-/HER2+",
                  "HR+/HER2+": "HR+/HER2+",
                  "TNBC": "TNBC",
                  "HER2+": "HER2+_combined"}
        target_c = rename.get(target)
        if target_c is None:
            continue
        # If CI is missing, approximate SE conservatively from p-value or assume
        # log(OR) +/- 0.5
        if pd.notna(r["ci_low"]) and pd.notna(r["ci_high"]) and r["ci_low"] > 0:
            se = (math.log(r["ci_high"]) - math.log(r["ci_low"])) / (2 * 1.96)
        else:
            # Use a conservative wide SE: 0.5 on the log-OR scale
            se = 0.5
        derived.append({"study": r["study_id"], "target": target_c,
                        "log_or": math.log(OR), "var": se ** 2})
    d = pd.DataFrame(derived)
    out = []
    for tgt in ("TNBC", "HR-/HER2+", "HR+/HER2+", "HER2+_combined"):
        sub = d[d["target"] == tgt]
        if sub.empty:
            continue
        res = dersimonian_laird(sub["log_or"].values, sub["var"].values)
        out.append({
            "contrast": f"{tgt} vs HR+/HER2-",
            "k_studies": res["k"],
            "pooled_log_or": res["pooled"],
            "pooled_or": math.exp(res["pooled"]),
            "ci_low_or": math.exp(res["ci_low"]),
            "ci_high_or": math.exp(res["ci_high"]),
            "tau2": res["tau2"],
            "i2_percent": res["i2"],
        })
    return pd.DataFrame(out), d


# ======================================================================
# Part C - build nomogram coefficients
# ======================================================================
def build_nomogram(pooled_subtype: pd.DataFrame,
                   pooled_ors: pd.DataFrame) -> pd.DataFrame:
    """Build a multivariable logistic equation:
       logit(P(ypN0)) = a0
           + a_subtype[ subtype ]
           + a_ki67  * I(Ki67 high category)
           + a_age   * (age - 50) / 10
           + a_cN23  * I(cN2/3)
           + a_pembro* I(pembrolizumab)
           + a_pertuz* I(pertuzumab)
           + a_carbo * I(carboplatin)
           + a_genolow*I(low genomic risk, only HR+/HER2-)

    Intercept = logit of pooled rate for HR+/HER2- reference.
    Subtype offsets = pooled log-OR vs HR+/HER2-.
    Other coefficients = pooled or literature priors with documented sources.
    """
    coef = []

    # Intercept: pooled logit of HR+/HER2- rate
    ref_row = pooled_subtype[pooled_subtype["subtype"] == "HR+/HER2-"].iloc[0]
    a0 = logit(float(ref_row["pooled_rate"]))
    coef.append(dict(term="intercept",
                     description="log-odds of ypN0 for the reference: HR+/HER2-, Ki67<20%, age>=50, cN1, anthracycline-taxane without pembro/pertuzumab/carbo, average genomic risk",
                     beta=a0,
                     or_=math.exp(a0),
                     se=float(ref_row.get("tau2_logit", 0.0)) ** 0.5,
                     source="DL pooled subtype rate (Part A)"))

    # Subtype offsets
    pool_map = {row["contrast"].split(" vs ")[0]: row
                for _, row in pooled_ors.iterrows()}
    for label, key in (("subtype_HR-/HER2+", "HR-/HER2+"),
                       ("subtype_HR+/HER2+", "HR+/HER2+"),
                       ("subtype_TNBC", "TNBC")):
        r = pool_map.get(key)
        if r is None:
            continue
        coef.append(dict(term=label,
                         description=f"log-OR for {key} relative to HR+/HER2-",
                         beta=float(r["pooled_log_or"]),
                         or_=math.exp(float(r["pooled_log_or"])),
                         se=float(r.get("tau2", 0.0)) ** 0.5,
                         source="DL pooled log-OR (Part B, ratesderived + direct ORs)"))

    # Ki67 (categorical, HR+/HER2- specific)
    # Boughey 2022 Ki67: rates 28.4% (Ki67>=20%) vs 8.1% (<20%) -> OR
    or_ki67 = (0.284 / (1 - 0.284)) / (0.081 / (1 - 0.081))
    coef.append(dict(term="ki67_20_30",
                     description="Ki67 20-30% vs <20% (partial effect; half of full Ki67>=30% effect)",
                     beta=math.log(or_ki67) / 2,
                     or_=math.exp(math.log(or_ki67) / 2),
                     se=0.20,
                     source="Boughey 2022 Ann Surg Oncol (PMID 35569077), modelled as graded"))
    coef.append(dict(term="ki67_gt30",
                     description="Ki67 >=30% vs <20% (HR+/HER2- only)",
                     beta=math.log(or_ki67),
                     or_=or_ki67,
                     se=0.20,
                     source="Boughey 2022 (rates 28.4 vs 8.1%); approximate"))

    # Age per decade (continuous, centred at 50)
    # An2025 NCDB: OR=0.99/year -> log-OR = log(0.99) per year, *10 per decade.
    coef.append(dict(term="age_per_decade",
                     description="OR per decade above 50 (younger -> higher ypN0)",
                     beta=10 * math.log(0.99),
                     or_=math.exp(10 * math.log(0.99)),
                     se=0.02,
                     source="An 2025 NCDB (PMID 40684017); OR 0.99/yr"))

    # cN stage cN2/3 vs cN1
    # Z1071 (Boughey 2013): cN1 41.0% (215/525) vs cN2 46.1% (12/26)
    p_cN1 = 215 / 525
    p_cN2 = 12 / 26
    or_cN23 = (p_cN2 / (1 - p_cN2)) / (p_cN1 / (1 - p_cN1))
    coef.append(dict(term="cN_stage_cN23",
                     description="cN2/3 vs cN1 (slight positive effect)",
                     beta=math.log(or_cN23),
                     or_=or_cN23,
                     se=0.30,
                     source="Boughey 2013 ACOSOG Z1071; small subgroup"))

    # Regimens
    # Pembrolizumab in TNBC: KEYNOTE-522 76.7% (pembro) vs 69.9% (control)
    or_pembro = (0.767 / (1 - 0.767)) / (0.699 / (1 - 0.699))
    coef.append(dict(term="regimen_pembro",
                     description="Pembrolizumab vs no IO (TNBC only, additive effect)",
                     beta=math.log(or_pembro),
                     or_=or_pembro,
                     se=0.20,
                     source="Kuemmel 2026 KEYNOTE-522 (PMID 41739398)"))
    # Pertuzumab dual HER2 effect: literature consensus +10-15 ypN0% over single
    # HER2 blockade. We anchor with NeoSphere-type effect: OR ~ 1.7 (assumed).
    coef.append(dict(term="regimen_pertuzumab",
                     description="Dual HER2 (T+P) vs trastuzumab alone",
                     beta=math.log(1.7),
                     or_=1.7,
                     se=0.20,
                     source="Literature prior; NeoSphere/TRYPHAENA pattern"))
    # Carboplatin in TNBC: BrighTNess/CALGB40603 +13% pCR -> OR ~ 1.5 axillary
    coef.append(dict(term="regimen_carbo",
                     description="Carboplatin addition in TNBC vs no carbo",
                     beta=math.log(1.5),
                     or_=1.5,
                     se=0.20,
                     source="NeoSTOP / BrighTNess pattern; conservative"))

    # Genomic risk (HR+/HER2- only): high vs low
    # Boughey 2022 genomic: high genomic risk 25% vs low 10.4%
    or_geno = (0.25 / 0.75) / (0.104 / 0.896)
    coef.append(dict(term="genomic_high",
                     description="High genomic risk (RS>25 or MammaPrint high) vs low (HR+/HER2- only)",
                     beta=math.log(or_geno),
                     or_=or_geno,
                     se=0.25,
                     source="Boughey 2022 genomic (PMID 35876925)"))

    return pd.DataFrame(coef)


def predict_logit(features: dict, coef_df: pd.DataFrame) -> float:
    """Compute logit(P(ypN0)) from a feature dict and coefficient table."""
    coef = dict(zip(coef_df["term"], coef_df["beta"]))
    z = coef["intercept"]
    # Subtype offset
    st = features.get("subtype", "HR+/HER2-")
    if st == "HR-/HER2+" and "subtype_HR-/HER2+" in coef:
        z += coef["subtype_HR-/HER2+"]
    elif st == "HR+/HER2+" and "subtype_HR+/HER2+" in coef:
        z += coef["subtype_HR+/HER2+"]
    elif st == "TNBC" and "subtype_TNBC" in coef:
        z += coef["subtype_TNBC"]
    # Ki67 (HR+/HER2- only)
    if st == "HR+/HER2-":
        ki67_cat = features.get("ki67_cat", "lt20")
        if ki67_cat == "20_30":
            z += coef.get("ki67_20_30", 0)
        elif ki67_cat == "gt30":
            z += coef.get("ki67_gt30", 0)
        # Genomic risk (only HR+/HER2-)
        if features.get("genomic_high", False):
            z += coef.get("genomic_high", 0)
    # Age (per decade above 50; positive coef means older -> lower ypN0,
    # but we anchored beta as log(0.99)*10 which is negative)
    age = features.get("age", 50)
    z += coef.get("age_per_decade", 0) * (age - 50) / 10
    # cN stage
    if features.get("cN23", False):
        z += coef.get("cN_stage_cN23", 0)
    # Regimens
    if features.get("pembro", False) and st == "TNBC":
        z += coef.get("regimen_pembro", 0)
    if features.get("pertuzumab", False) and ("HER2+" in st or "HER2+_combined" in st):
        z += coef.get("regimen_pertuzumab", 0)
    if features.get("carbo", False) and st == "TNBC":
        z += coef.get("regimen_carbo", 0)
    return z


def predict_prob(features: dict, coef_df: pd.DataFrame) -> float:
    return inv_logit(predict_logit(features, coef_df))


# ======================================================================
# Part D - validation against pivotal studies
# ======================================================================
def validate(coef_df: pd.DataFrame) -> pd.DataFrame:
    """Reconstruct expected ypN0 from study-level case-mix and compare.

    Each entry encodes the published mix (subtype proportions, regimens)."""
    studies = [
        # Z1071 (Boughey 2013): cN1, anthracycline-taxane, pre-pembro era.
        # Published mix (JAMA 2013 Table 1): HR+/HER2- 56.0%, HER2+ 26.0%,
        # TNBC 18.0%. Pertuzumab not used (trial 2009-2011).
        {
            "study": "Z1071 (Boughey 2013)",
            "observed": 0.41,
            "mix": [
                ("TNBC", 0.18, {"pembro": False, "carbo": False}),
                ("HR-/HER2+", 0.08, {"pertuzumab": False}),
                ("HR+/HER2+", 0.18, {"pertuzumab": False}),
                ("HR+/HER2-", 0.56, {"ki67_cat": "lt20"}),
            ],
            "age": 50, "cN23": False,
        },
        # RISAS (Simons 2022 / van Amstel 2024): Dutch protocols 2017-19.
        # Mix per RISAS publication (HR+/HER2- ~57%, HER2+ ~25%, TNBC ~18%).
        # Pertuzumab routine in HER2+, no pembro in TNBC (pre-KEYNOTE-522 era).
        {
            "study": "RISAS (Simons 2022)",
            "observed": 0.304,
            "mix": [
                ("TNBC", 0.18, {"pembro": False, "carbo": False}),
                ("HR-/HER2+", 0.12, {"pertuzumab": True}),
                ("HR+/HER2+", 0.13, {"pertuzumab": True}),
                ("HR+/HER2-", 0.57, {"ki67_cat": "lt20"}),
            ],
            "age": 52, "cN23": False,
        },
        # Montagna 2024: cohort selected ypN0 (rate = 1.0 by design) but mix
        # informs HER2+ subset. Mix: 54% HER2+, modern era with pembro+pertuz.
        # Expected: predicted will be HIGH but <1 (selection-defined).
        {
            "study": "Montagna 2024 (modern mix)",
            "observed": 1.0,  # selection-defined
            "mix": [
                ("HR-/HER2+", 0.27, {"pertuzumab": True}),
                ("HR+/HER2+", 0.27, {"pertuzumab": True}),
                ("TNBC", 0.20, {"pembro": True, "carbo": True}),
                ("HR+/HER2-", 0.26, {"ki67_cat": "20_30"}),
            ],
            "age": 50, "cN23": False,
        },
        # van Amstel 2024 (RISAS overall): same trial as Simons, 29.7% observed
        {
            "study": "van Amstel 2024 (RISAS overall)",
            "observed": 0.297,
            "mix": [
                ("TNBC", 0.18, {"pembro": False, "carbo": False}),
                ("HR-/HER2+", 0.12, {"pertuzumab": True}),
                ("HR+/HER2+", 0.13, {"pertuzumab": True}),
                ("HR+/HER2-", 0.57, {"ki67_cat": "lt20"}),
            ],
            "age": 52, "cN23": False,
        },
        # KEYNOTE-522 (Kuemmel 2026) pembro arm, pure TNBC
        {
            "study": "KEYNOTE-522 pembro arm (TNBC)",
            "observed": 0.767,
            "mix": [
                ("TNBC", 1.0, {"pembro": True, "carbo": True}),
            ],
            "age": 49, "cN23": False,
        },
        # KEYNOTE-522 control arm
        {
            "study": "KEYNOTE-522 placebo arm (TNBC)",
            "observed": 0.699,
            "mix": [
                ("TNBC", 1.0, {"pembro": False, "carbo": True}),
            ],
            "age": 49, "cN23": False,
        },
    ]
    rows = []
    for s in studies:
        expected = 0.0
        for subtype, prop, regimen in s["mix"]:
            feats = {"subtype": subtype, "age": s["age"], "cN23": s["cN23"]}
            feats.update(regimen)
            p = predict_prob(feats, coef_df)
            expected += prop * p
        rows.append({
            "study": s["study"],
            "predicted_ypN0": expected,
            "observed_ypN0": s["observed"],
            "absolute_error": abs(expected - s["observed"]),
        })
    val = pd.DataFrame(rows)
    # Brier-like calibration: mean squared error on those with observed != 1
    mask = val["observed_ypN0"] < 1.0
    brier = float(np.mean((val.loc[mask, "predicted_ypN0"] -
                           val.loc[mask, "observed_ypN0"]) ** 2))
    val.attrs["brier"] = brier
    return val


# ======================================================================
# Main pipeline
# ======================================================================
def equation_text(coef_df: pd.DataFrame) -> tuple[str, str]:
    """Return a (latex, plain) representation of the equation."""
    by_term = dict(zip(coef_df["term"], coef_df["beta"]))
    latex_parts = [r"\mathrm{logit}\,P(ypN0) = ", f"{by_term['intercept']:+.3f}"]
    plain_parts = ["logit P(ypN0) =", f"{by_term['intercept']:+.3f}"]
    # subtype dummies
    for k, lab in (("subtype_HR-/HER2+", "HR-/HER2+"),
                   ("subtype_HR+/HER2+", "HR+/HER2+"),
                   ("subtype_TNBC", "TNBC")):
        if k in by_term:
            latex_parts.append(f"{by_term[k]:+.3f}\\cdot \\mathbb{{1}}[{lab}]")
            plain_parts.append(f"{by_term[k]:+.3f} * 1[{lab}]")
    for k, lab in (("ki67_20_30", "Ki67 20-30%"),
                   ("ki67_gt30", "Ki67 >=30%"),
                   ("genomic_high", "high genomic risk")):
        if k in by_term:
            latex_parts.append(f"{by_term[k]:+.3f}\\cdot \\mathbb{{1}}[{lab}]")
            plain_parts.append(f"{by_term[k]:+.3f} * 1[{lab}]")
    if "age_per_decade" in by_term:
        latex_parts.append(f"{by_term['age_per_decade']:+.3f}\\cdot (age-50)/10")
        plain_parts.append(f"{by_term['age_per_decade']:+.3f} * (age-50)/10")
    if "cN_stage_cN23" in by_term:
        latex_parts.append(f"{by_term['cN_stage_cN23']:+.3f}\\cdot \\mathbb{{1}}[cN2/3]")
        plain_parts.append(f"{by_term['cN_stage_cN23']:+.3f} * 1[cN2/3]")
    for k, lab in (("regimen_pembro", "pembro"),
                   ("regimen_pertuzumab", "pertuzumab"),
                   ("regimen_carbo", "carbo")):
        if k in by_term:
            latex_parts.append(f"{by_term[k]:+.3f}\\cdot \\mathbb{{1}}[{lab}]")
            plain_parts.append(f"{by_term[k]:+.3f} * 1[{lab}]")
    return " ".join(latex_parts), " ".join(plain_parts)


def main() -> dict:
    df_rates = load_subtype_rates()
    pooled = pool_subtype_rates(df_rates)
    pooled_path = OUT_HERE / "results_pooled_subtype.csv"
    pooled.to_csv(pooled_path, index=False)

    pooled_ors, derived = pool_subtype_logors(df_rates)
    or_path = OUT_HERE / "results_meta_or_subtype.csv"
    pooled_ors.to_csv(or_path, index=False)

    coef = build_nomogram(pooled, pooled_ors)
    coef_path = OUT_DATA / "nomogram_coefficients.csv"
    coef.to_csv(coef_path, index=False)

    val = validate(coef)
    val_path = OUT_HERE / "results_validation.csv"
    val.to_csv(val_path, index=False)

    latex_eq, plain_eq = equation_text(coef)

    # Print human-readable summary
    print("=" * 72)
    print("PART A - Pooled per-subtype ypN0 rates (DerSimonian-Laird RE)")
    print("=" * 72)
    print(pooled.to_string(index=False, float_format=lambda x: f"{x:.3f}"))
    print()
    print("=" * 72)
    print("PART B - Pooled log-OR by subtype (reference HR+/HER2-)")
    print("=" * 72)
    print(pooled_ors.to_string(index=False, float_format=lambda x: f"{x:.3f}"))
    print()
    print("=" * 72)
    print("PART C - Nomogram coefficients (logit scale)")
    print("=" * 72)
    print(coef.to_string(index=False, float_format=lambda x: f"{x:.3f}"))
    print()
    print("PLAIN EQUATION:")
    print(" ", plain_eq)
    print()
    print("LATEX EQUATION:")
    print(" ", latex_eq)
    print()
    print("=" * 72)
    print("PART D - Validation predicted vs observed")
    print("=" * 72)
    print(val.to_string(index=False, float_format=lambda x: f"{x:.3f}"))
    print(f"\nBrier score (excluding selection-defined Montagna): {val.attrs['brier']:.4f}")
    print()
    return {
        "pooled_subtype_csv": str(pooled_path),
        "pooled_or_csv": str(or_path),
        "nomogram_coef_csv": str(coef_path),
        "validation_csv": str(val_path),
        "equation_latex": latex_eq,
        "equation_plain": plain_eq,
        "brier": val.attrs["brier"],
    }


if __name__ == "__main__":
    main()
