1"""H2 (and H4/H5) on Seshat: Y1 crisis severity (primary), Y1b severe binary, Y2 polity duration (Cox), Y3 succession.2This is the first module that reads outcomes together with P. Usage: python -m wp35.econometrics.h2_seshat3Everything follows research/preregistration.md §5.2–5.3, §5.5–5.7, §6.4"""5from __future__ import annotations6import json7from pathlib import Path8import numpy as np9import pandas as pd10import statsmodels.api as sm11import statsmodels.formula.api as smf12from lifelines import CoxPHFitter13from lifelines.statistics import proportional_hazard_test14from wp35.econometrics.common import fit_formula, coef_row, cv_compare, fold_win_share, hinge_fit, binary_metrics, oos_loglik_binomial1516ROOT = Path(__file__).resolve().parents[3]17PROC = ROOT / "data" / "processed"; RES = ROOT / "results" / "empirical"18RES.mkdir(parents=True, exist_ok=True)19BIN = sm.families.Binomial()20FE = "C(macro_region) + C(era)"21XS = "X_warfare + X_prod_infra" # epidemic_v dropped (A-006): it is itself a CrisisDB consequence item222324def load_crisis():25 cr = pd.read_parquet(PROC / "seshat_crisis_indices.parquet")26 y = pd.read_parquet(PROC / "seshat_crisis.parquet")[["crisis_id", "y1_severity", "y1b_severe", "y_reform", "epidemic_v", "y1_n_coded"]]27 d = cr.merge(y, on="crisis_id")28 meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet")[["polity_id", "centralization", "homeNGA"]]29 d = d.merge(meta, on="polity_id", how="left", suffixes=("", "_m"))30 d["homeNGA"] = d["homeNGA"].fillna(d.get("homeNGA_m"))31 # persistence baseline: an earlier crisis record in the same NGA32 d = d.sort_values("crisis_start")33 d["prior_crisis_nga"] = 0.034 for i, r in d.iterrows():35 d.loc[i, "prior_crisis_nga"] = float(((d.homeNGA == r.homeNGA) & (d.crisis_start < r.crisis_start)).any())36 d["quasi"] = d.centralization.eq("quasi-polity")37 d["P_placebo"] = (d.D_pc1 - d.D_pc1.mean()) / d.D_pc1.std() - (d.X_prod_infra - d.X_prod_infra.mean()) / d.X_prod_infra.std()38 for c in ["P_main", "D_pc1", "K_mean", "K_irt", "complexity_pc1", "log_pop"]:39 d[c + "_z"] = (d[c] - d[c].mean()) / d[c].std()40 return d414243def wald_sum_zero(res, a: str, b: str) -> dict:44 """Wald test of beta_a + beta_b = 0 (the restriction imposed by P = z(D) - z(K))."""45 import scipy.stats as st46 V = res.cov_params(); s_ = res.params[a] + res.params[b]; var = V.loc[a, a] + V.loc[b, b] + 2 * V.loc[a, b]47 w = s_ ** 2 / var; return {"sum_coef": float(s_), "sum_se": float(np.sqrt(var)), "wald_p": float(1 - st.chi2.cdf(w, 1))}484950def primary_models(d: pd.DataFrame):51 rows = []52 # common estimation sample for every specification that uses covariates/baselines (A-006)53 common = d.dropna(subset=["y1_severity", "P_main_z", "D_pc1_z", "K_mean_z", "complexity_pc1_z", "prior_crisis_nga", "X_warfare", "X_prod_infra", "P_placebo", "macro_region", "era"])54 specs = {55 "A1 primary: Y1 ~ P + FE": (f"y1_severity ~ P_main_z + {FE}", "P_main_z"),56 "A2: Y1 ~ P + X + FE": (f"y1_severity ~ P_main_z + {XS} + {FE}", "P_main_z"),57 "A3: Y1 ~ P + complexity PC1 + X + FE": (f"y1_severity ~ P_main_z + complexity_pc1_z + {XS} + {FE}", "P_main_z"),58 "B1 baseline D only": (f"y1_severity ~ D_pc1_z + {XS} + {FE}", "D_pc1_z"),59 "B2 baseline K only": (f"y1_severity ~ K_mean_z + {XS} + {FE}", "K_mean_z"),60 "B3 baseline complexity PC1": (f"y1_severity ~ complexity_pc1_z + {XS} + {FE}", "complexity_pc1_z"),61 "B4 baseline persistence": (f"y1_severity ~ prior_crisis_nga + {XS} + {FE}", "prior_crisis_nga"),62 "B5 D and K separately": (f"y1_severity ~ D_pc1_z + K_mean_z + {XS} + {FE}", "D_pc1_z"),63 "Placebo: P built with production infrastructure instead of K (X = warfare only)": (f"y1_severity ~ P_placebo + X_warfare + {FE}", "P_placebo"),64 "Placebo comparison: P with the same controls (X = warfare only)": (f"y1_severity ~ P_main_z + X_warfare + {FE}", "P_main_z"),65 }66 for label, (f, term) in specs.items():67 sample = d if label.startswith("A1") else common68 r = fit_formula(f, sample, BIN, groups="seshat_region"); rows.append(coef_row(r, term, label))69 if "B5" in label:70 rows.append(coef_row(r, "K_mean_z", label)); wt = wald_sum_zero(r, "D_pc1_z", "K_mean_z")71 rows.append({"model": "B5 Wald test of beta_D + beta_K = 0 (restriction implied by P)", "term": "D+K", "coef": wt["sum_coef"], "se": wt["sum_se"], "p": wt["wald_p"], "n": r._n_used})72 # binary secondary73 # binary secondary outcome: several region/era cells are all-severe (quasi-separation) -> region FE only, HC1 SE74 for label, f in {"A1b: Y1b ~ P + region FE (logit, HC1; era FE dropped: separation)": "y1b_severe ~ P_main_z + C(macro_region)",75 "A2b: Y1b ~ P + X + region FE (logit, HC1)": f"y1b_severe ~ P_main_z + {XS} + C(macro_region)",76 "A3b: Y1b ~ P (logit, HC1, no FE)": "y1b_severe ~ P_main_z"}.items():77 r = fit_formula(f, d, BIN, groups=None); rows.append(coef_row(r, "P_main_z", label))78 # exploratory: reforms79 r = fit_formula(f"y_reform ~ P_main_z + {XS} + {FE}", d, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "Exploratory: reform consequences ~ P + X + FE"))80 return pd.DataFrame(rows)818283def robustness(d: pd.DataFrame):84 rows = []85 variants = {"P_D_pc1_K_irt": "K = IRT latent", "P_log_pop_K_mean": "D = log population only", "P_D_pc1_K_pc1": "K = PC1 of blocks",86 "P_main_lo": "D at lower bounds", "P_main_hi": "D at upper bounds", "P_main_noinf": "inferred codes dropped",87 "P_within_era": "standardized within era", "P_within_region": "standardized within macro-region", "P_rank": "rank-based P",88 "P_nosh": "D without settlement hierarchy", "P_noadmlev": "K without administrative levels", "P_nosh_noadmlev": "both hierarchy indicators removed"}89 for v, lab in variants.items():90 if v not in d.columns:91 # crisis-level file lacks some polity-level variants → merge from polity file92 continue93 dd = d.copy(); dd["_P"] = (dd[v] - dd[v].mean()) / dd[v].std()94 r = fit_formula(f"y1_severity ~ _P + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "_P", f"R: {lab}"))95 dd = d[~d.quasi]; r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: quasi-polities dropped"))96 dd = d[d.years_since_start >= 25]; r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: crises starting ≥25 years after polity start"))97 dd = d[d.y1_n_coded >= 12]; r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: ≥12 of 14 items coded"))98 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + C(macro_region)", d, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: no era FE"))99 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + C(era)", d, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: no region FE"))100 r = fit_formula(f"y1_severity ~ P_main_z + {XS}", d, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: no FE"))101 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", d, BIN, groups="homeNGA"); rows.append(coef_row(r, "P_main_z", "R: SE clustered by NGA"))102 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", d, BIN, groups="polity_id"); rows.append(coef_row(r, "P_main_z", "R: SE clustered by polity"))103 return pd.DataFrame(rows)104105106def out_of_sample(d: pd.DataFrame):107 formulas = {"A2 P": f"y1_severity ~ P_main_z + {XS} + {FE}", "B0 FE only": f"y1_severity ~ {XS} + {FE}", "B1 D only": f"y1_severity ~ D_pc1_z + {XS} + {FE}",108 "B2 K only": f"y1_severity ~ K_mean_z + {XS} + {FE}", "B3 complexity": f"y1_severity ~ complexity_pc1_z + {XS} + {FE}", "B4 persistence": f"y1_severity ~ prior_crisis_nga + {XS} + {FE}",109 "B5 D+K": f"y1_severity ~ D_pc1_z + K_mean_z + {XS} + {FE}"}110 # era FE with leave-one-region-out can leave unseen levels; macro-region FE cannot be used when the region is held out111 f_loro = {k: v.replace("C(macro_region) + ", "") for k, v in formulas.items()}112 loro = cv_compare(d, f_loro, "y1_severity", BIN, fold_col="macro_region", fractional=True)113 kf = cv_compare(d, formulas, "y1_severity", BIN, n_repeats=10, k=5, fractional=True, block_col="polity_id")114 loro["scheme"] = "leave-one-macro-region-out"; kf["scheme"] = "10x5-fold"115 cv = pd.concat([loro, kf]); cv.to_csv(RES / "h2_seshat_cv_folds.csv", index=False)116 wins = []117 for scheme, c in [("loro", loro), ("kfold", kf)]:118 for b in ["B0 FE only", "B1 D only", "B2 K only", "B3 complexity", "B4 persistence", "B5 D+K"]:119 w = fold_win_share(c, "A2 P", b); w["scheme"] = scheme; wins.append(w)120 wins = pd.DataFrame(wins); wins.to_csv(RES / "h2_seshat_cv_wins.csv", index=False)121 # binary Y1b out-of-sample metrics122 fb = {k: v.replace("y1_severity", "y1b_severe") for k, v in formulas.items()}123 kfb = cv_compare(d, fb, "y1b_severe", BIN, n_repeats=10, k=5, block_col="polity_id")124 kfb.to_csv(RES / "h2_seshat_cv_folds_binary.csv", index=False)125 summ = kfb.groupby("model")[["oos_loglik", "auc", "brier"]].mean().reset_index(); summ.to_csv(RES / "h2_seshat_cv_binary_summary.csv", index=False)126 return wins, summ127128129def h4_nonlinear(d: pd.DataFrame):130 """Linear vs hinge vs spline, out of sample (10x5-fold); hinge knot chosen on the training fold only."""131 rng = np.random.default_rng(20261101)132 base = f"y1_severity ~ P_main_z + {XS} + {FE}"133 dd = d.dropna(subset=["y1_severity", "P_main_z", "X_warfare", "X_prod_infra", "macro_region", "era"]).reset_index(drop=True)134 rows = []; pols = dd.polity_id.unique()135 from wp35.econometrics.common import map_unseen_levels136 for rep in range(10):137 perm = rng.permutation(len(pols)); parts = np.array_split(perm, 5)138 for i, te in enumerate(parts):139 m = dd.polity_id.isin(pols[te]).values; tr, ts = dd[~m], map_unseen_levels(dd[~m], dd[m], base)140 try:141 lin = smf.glm(base, tr, family=BIN).fit(); ll_lin = oos_loglik_binomial(ts.y1_severity, lin.predict(ts))142 h = hinge_fit(tr, "y1_severity", base, "P_main_z", BIN)143 ts["_hinge"] = np.clip(ts.P_main_z - h[1], 0, None); ll_h = oos_loglik_binomial(ts.y1_severity, h[2].predict(ts))144 sp = smf.glm(f"y1_severity ~ cr(P_main_z, df=3) + {XS} + {FE}", tr, family=BIN).fit(); ll_s = oos_loglik_binomial(ts.y1_severity, sp.predict(ts))145 rows.append({"fold": f"{rep}_{i}", "ll_linear": ll_lin, "ll_hinge": ll_h, "ll_spline": ll_s, "knot": h[1]})146 except Exception as e:147 rows.append({"fold": f"{rep}_{i}", "error": str(e)[:80]})148 r = pd.DataFrame(rows); r.to_csv(RES / "h4_seshat_folds.csv", index=False)149 ok = r.dropna(subset=["ll_linear"])150 summ = {"folds": int(len(ok)), "hinge_beats_linear_share": float((ok.ll_hinge > ok.ll_linear).mean()), "spline_beats_linear_share": float((ok.ll_spline > ok.ll_linear).mean()),151 "mean_ll_linear": float(ok.ll_linear.mean()), "mean_ll_hinge": float(ok.ll_hinge.mean()), "mean_ll_spline": float(ok.ll_spline.mean()), "median_knot": float(ok.knot.median())}152 # in-sample hinge for description153 h = hinge_fit(dd, "y1_severity", base, "P_main_z", BIN); summ["insample_knot"] = float(h[1]); summ["insample_hinge_coef"] = float(h[2].params["_hinge"]); summ["insample_P_coef"] = float(h[2].params["P_main_z"])154 return summ155156157def item_level(d: pd.DataFrame):158 """Per-item logits (A-002): does P (and D) predict each CrisisDB consequence? Plus dimensionality diagnostics."""159 from wp35.harmonization.seshat import INSTABILITY_ITEMS160 cols = [i + "_v" for i in INSTABILITY_ITEMS]161 full = pd.read_parquet(PROC / "seshat_crisis.parquet").set_index("crisis_id")[cols]162 dd = d.set_index("crisis_id").join(full, rsuffix="_full")163 rows = []164 for it in INSTABILITY_ITEMS:165 y = dd[it + "_v_full"] if it + "_v_full" in dd else dd[it + "_v"]166 sub = dd.assign(_y=y).dropna(subset=["_y", "P_main_z", "D_pc1_z"]); sub = sub[sub._y.isin([0.0, 1.0])]167 row = {"item": it, "n": int(len(sub)), "share_present": float(sub._y.mean())}168 for term in ["P_main_z", "D_pc1_z", "K_mean_z"]:169 try:170 m = smf.glm(f"_y ~ {term} + C(era)", sub, family=BIN).fit(cov_type="HC1"); row[f"coef_{term}"] = m.params[term]; row[f"se_{term}"] = m.bse[term]171 except Exception:172 row[f"coef_{term}"] = np.nan; row[f"se_{term}"] = np.nan173 rows.append(row)174 out = pd.DataFrame(rows); out.to_csv(RES / "h2_seshat_items.csv", index=False)175 M = full.replace({0.5: np.nan}); C = M.corr(min_periods=30); off = C.values[np.triu_indices(len(cols), 1)]176 ev = np.linalg.eigvalsh(C.fillna(0).values)177 diag = {"mean_interitem_corr": float(np.nanmean(off)), "median_interitem_corr": float(np.nanmedian(off)), "first_eigen_share": float(ev[-1] / ev.sum()), "n_items": len(cols)}178 (RES / "h2_seshat_items_dimensionality.json").write_text(json.dumps(diag, indent=2))179 return out, diag180181182def h5_blocks(d: pd.DataFrame):183 rows = []184 for col, name in [("macro_region", "region"), ("era", "era"), ("cultural_sphere", "sphere")]:185 for g in sorted(d[col].dropna().unique()):186 held = d[d[col] != g]; inb = d[d[col] == g]187 other = "C(era)" if col in ("macro_region", "cultural_sphere") else "C(macro_region)"188 try:189 try:190 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {other}", held, BIN, groups="seshat_region"); cov = "cluster"191 except Exception:192 try:193 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {other}", held, BIN, groups=None); cov = "HC1 (cluster cov singular)"194 except Exception:195 r = fit_formula(f"y1_severity ~ P_main_z + {XS}", held, BIN, groups="seshat_region"); cov = "cluster, FE dropped (design singular)"196 row = coef_row(r, "P_main_z", f"H5 leave-{name}-out: {g}"); row["held_out_n"] = int(inb.y1_severity.notna().sum()); row["cov"] = cov197 if inb.P_main_z.notna().sum() >= 12:198 try:199 rb = fit_formula(f"y1_severity ~ P_main_z + {XS}", inb, BIN, groups="seshat_region"); row["within_block_coef"] = rb.params["P_main_z"]; row["within_block_se"] = rb.bse["P_main_z"]200 except Exception:201 try:202 rb = fit_formula("y1_severity ~ P_main_z", inb, BIN, groups=None); row["within_block_coef"] = rb.params["P_main_z"]; row["within_block_se"] = rb.bse["P_main_z"]; row["within_block_note"] = "P only, HC1"203 except Exception as e2:204 row["within_block_note"] = str(e2)[:60]205 rows.append(row)206 except Exception as e:207 rows.append({"model": f"H5 {name} {g}", "error": str(e)[:80]})208 out = pd.DataFrame(rows); out.to_csv(RES / "h5_seshat_blocks.csv", index=False)209 return out210211212def y2_cox_and_y3(dp: pd.DataFrame):213 rows = []214 d = dp.dropna(subset=["P_main", "duration", "macro_region"]).copy()215 d["event"] = (~d.censored).astype(int); d["P_z"] = (d.P_main - d.P_main.mean()) / d.P_main.std()216 for c in ["X_warfare", "X_prod_infra", "X_military", "D_pc1", "K_mean", "complexity_pc1", "n_polities_in_nga"]:217 d[c + "_z"] = (d[c] - d[c].mean()) / d[c].std()218 d["era_num"] = d.midpoint219 specs = {"Y2-1 Cox: P + strata(macro-region)": ["P_z"], "Y2-2 Cox: P + X + midpoint + strata": ["P_z", "X_warfare_z", "X_prod_infra_z", "midpoint"],220 "Y2-3 Cox: D + K + X + midpoint + strata": ["D_pc1_z", "K_mean_z", "X_warfare_z", "X_prod_infra_z", "midpoint"], "Y2-4 Cox: complexity PC1 + X + strata": ["complexity_pc1_z", "X_warfare_z", "X_prod_infra_z", "midpoint"],221 "Y2-5 Cox: D + K + X + midpoint + polities per NGA (periodization) + strata": ["D_pc1_z", "K_mean_z", "X_warfare_z", "X_prod_infra_z", "midpoint", "n_polities_in_nga_z"]}222 for label, cols in specs.items():223 dd = d.dropna(subset=cols + ["seshat_region"])[cols + ["duration", "event", "macro_region", "seshat_region"]].copy()224 try:225 cph = CoxPHFitter(); cph.fit(dd, duration_col="duration", event_col="event", strata=["macro_region"], cluster_col="seshat_region", robust=True)226 for t in (cols[:2] + (["n_polities_in_nga_z"] if "Y2-5" in label else [])) if ("Y2-3" in label or "Y2-5" in label) else cols[:1]:227 s = cph.summary.loc[t]228 rows.append({"model": label, "term": t, "coef": s["coef"], "se": s["se(coef)"], "ci_lo": s["coef lower 95%"], "ci_hi": s["coef upper 95%"], "p": s["p"], "n": len(dd), "events": int(dd.event.sum()), "hr": np.exp(s["coef"])})229 if "Y2-2" in label:230 ph = proportional_hazard_test(cph, dd, time_transform="rank"); rows[-1]["schoenfeld_p_P"] = float(ph.summary.loc["P_z", "p"])231 except Exception as e:232 rows.append({"model": label, "error": str(e)[:100]})233 # Y3 discontinuous succession234 meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet")[["polity_id", "y3_discontinuous"]]235 d3 = dp.merge(meta, on="polity_id").dropna(subset=["y3_discontinuous", "P_main"]).copy(); d3["P_z"] = (d3.P_main - d3.P_main.mean()) / d3.P_main.std()236 for label, f in {"Y3-1 logit: discontinuous succession ~ P + FE": "y3_discontinuous ~ P_z + C(macro_region) + C(era)", "Y3-2 logit: + X": "y3_discontinuous ~ P_z + X_warfare + X_prod_infra + C(macro_region) + C(era)"}.items():237 try:238 r = fit_formula(f, d3, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_z", label, {"events": int(r._data_used.y3_discontinuous.sum())}))239 except Exception as e:240 rows.append({"model": label, "error": str(e)[:100]})241 out = pd.DataFrame(rows); out.to_csv(RES / "h2_seshat_y2_y3.csv", index=False)242 return out243244245def main():246 d = load_crisis(); dp = pd.read_parquet(PROC / "seshat_indices.parquet")247 d.to_parquet(PROC / "seshat_h2_sample.parquet", index=False)248 prim = primary_models(d); prim.to_csv(RES / "h2_seshat_primary.csv", index=False)249 rob = robustness(d); rob.to_csv(RES / "h2_seshat_robustness.csv", index=False)250 wins, summ = out_of_sample(d)251 h4 = h4_nonlinear(d); (RES / "h4_seshat_summary.json").write_text(json.dumps(h4, indent=2))252 h5 = h5_blocks(d)253 items, idiag = item_level(d)254 y23 = y2_cox_and_y3(dp)255 print("\n=== ITEM LEVEL ===\n", items.round(3).to_string(index=False), "\n", idiag)256 print("corr(K_mean, polities per NGA) =", round(float(dp.K_mean.corr(dp.n_polities_in_nga)), 3))257 pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 70)258 cols = ["model", "term", "coef", "se", "p", "n"]259 print("=== PRIMARY (Y1 fractional logit, cluster region) ===\n", prim[cols].round(3).to_string(index=False))260 print("\n=== ROBUSTNESS ===\n", rob[cols].round(3).to_string(index=False))261 print("\n=== OOS wins of A2 over baselines ===\n", wins.round(3).to_string(index=False))262 print("\n=== Y1b OOS summary ===\n", summ.round(3).to_string(index=False))263 print("\n=== H4 ===", json.dumps(h4, indent=1))264 print("\n=== H5 ===\n", h5[[c for c in h5.columns if c in ("model", "coef", "se", "p", "n", "held_out_n", "within_block_coef", "within_block_se")]].round(3).to_string(index=False))265 print("\n=== Y2 Cox / Y3 ===\n", y23[[c for c in y23.columns if c in ("model", "term", "coef", "se", "p", "n", "events", "hr", "schoenfeld_p_P", "error")]].round(3).to_string(index=False))266267268if __name__ == "__main__":269 main()270