"""Independent re-implementation and robustness checks for Zenodo 10.5281/zenodo.22677480.
Usage: python3 -I replicate.py path/to/CORT_IT_open_data.csv
"""
import sys
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf

d = pd.read_csv(sys.argv[1], na_values=[""])
for c in ["is_verified", "stringent_verif", "no_or_wrong_memory"]:
    d[c] = d[c].map({True: True, False: False, "TRUE": True, "FALSE": False})
v = d[d.is_verified == True].copy()
v["lit"] = np.log10(v.intermp + 4)

def r(df, x, y="lit"):
    s = df[[x, y]].dropna()
    rr, p = stats.pearsonr(s[x], s[y])
    rs, ps = stats.spearmanr(s[x], s[y])
    return f"r={rr:.3f} p={p:.2g} rho={rs:.3f} N={len(s)}"

print("== Reproduction (Python, independent of authors' R) ==")
print("all cases", len(d), "verified", len(v))
print("unexp x logIT  ", r(v, "death_unexp"))
print("natural x logIT", r(v, "NATURAL_DEATH"))
print("age x logIT    ", r(v, "agedhp"))
m = smf.ols("lit ~ death_unexp + NATURAL_DEATH", v).fit()
print("joint OLS: b_unexp=%.4f p=%.3f | b_nat=%.4f p=%.3f | N=%d" % (
    m.params.death_unexp, m.pvalues.death_unexp, m.params.NATURAL_DEATH, m.pvalues.NATURAL_DEATH, m.nobs))
med = v.dropna(subset=["death_unexp", "intermp"]).groupby("death_unexp").intermp.agg(["median", "size"])
print(med)
print("drop 0->4: %.0f%%" % (100 * (1 - med.loc[4, "median"] / med.loc[0, "median"])))

print("\n== Robustness checks (not in authors' script) ==")
s = v.dropna(subset=["death_unexp", "intermp"])
print("IT range months:", s.intermp.min(), s.intermp.max(), "| negative IT n =", (s.intermp < 0).sum(),
      "| IT>20y n =", (s.intermp > 240).sum())
print("excl. negative IT        ", r(v[v.intermp >= 0], "death_unexp"))
print("excl. IT > 20 years      ", r(v[v.intermp <= 240], "death_unexp"))
print("excl. same family (4)    ", r(v[v.family_assoc != 4], "death_unexp"))
print("same family only         ", r(v[v.family_assoc == 4], "death_unexp"))
print("no family link (0) only  ", r(v[v.family_assoc == 0], "death_unexp"))
for k in [1, 10]:
    v["alt"] = np.log10(v.intermp + 3 + k)
    print(f"log constant +{3+k:<3}          ", r(v, "death_unexp", "alt"))

a = s[s.death_unexp == 4].intermp; b = s[s.death_unexp < 4].intermp
u = stats.mannwhitneyu(a, b)
print("Mann-Whitney IT, rating 4 vs 0-3: U=%.0f p=%.2g, medians %g vs %g" % (u.statistic, u.pvalue, a.median(), b.median()))
rng = np.random.default_rng(2026)
boots = [np.median(rng.choice(b, len(b))) - np.median(rng.choice(a, len(a))) for _ in range(5000)]
print("bootstrap 95%% CI for median difference (0-3 minus 4): [%.1f, %.1f] months" % tuple(np.percentile(boots, [2.5, 97.5])))
q = smf.quantreg("intermp ~ death_unexp", s).fit(q=0.5)
print("median regression slope: %.2f months per rating step, p=%.3g" % (q.params.death_unexp, q.pvalues.death_unexp))

print("\n== Coding/missingness checks ==")
v["unexp_missing"] = v.death_unexp.isna()
g = v.dropna(subset=["intermp"]).groupby("unexp_missing").intermp.agg(["median", "size"])
print("IT by whether unexpectedness was coded:\n", g)
print("unexpectedness rating distribution (%):\n", (100 * v.death_unexp.value_counts(normalize=True).sort_index()).round(0))
print("crosstab naturalness x unexpectedness:\n", pd.crosstab(v.NATURAL_DEATH, v.death_unexp))