"""H2 (and H4/H5) on Seshat: Y1 crisis severity (primary), Y1b severe binary, Y2 polity duration (Cox), Y3 succession. This is the first module that reads outcomes together with P. Usage: python -m wp35.econometrics.h2_seshat Everything follows research/preregistration.md §5.2–5.3, §5.5–5.7, §6. """ from __future__ import annotations import json from pathlib import Path import numpy as np import pandas as pd import statsmodels.api as sm import statsmodels.formula.api as smf from lifelines import CoxPHFitter from lifelines.statistics import proportional_hazard_test from wp35.econometrics.common import fit_formula, coef_row, cv_compare, fold_win_share, hinge_fit, binary_metrics, oos_loglik_binomial ROOT = Path(__file__).resolve().parents[3] PROC = ROOT / "data" / "processed"; RES = ROOT / "results" / "empirical" RES.mkdir(parents=True, exist_ok=True) BIN = sm.families.Binomial() FE = "C(macro_region) + C(era)" XS = "X_warfare + X_prod_infra" # epidemic_v dropped (A-006): it is itself a CrisisDB consequence item def load_crisis(): cr = pd.read_parquet(PROC / "seshat_crisis_indices.parquet") y = pd.read_parquet(PROC / "seshat_crisis.parquet")[["crisis_id", "y1_severity", "y1b_severe", "y_reform", "epidemic_v", "y1_n_coded"]] d = cr.merge(y, on="crisis_id") meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet")[["polity_id", "centralization", "homeNGA"]] d = d.merge(meta, on="polity_id", how="left", suffixes=("", "_m")) d["homeNGA"] = d["homeNGA"].fillna(d.get("homeNGA_m")) # persistence baseline: an earlier crisis record in the same NGA d = d.sort_values("crisis_start") d["prior_crisis_nga"] = 0.0 for i, r in d.iterrows(): d.loc[i, "prior_crisis_nga"] = float(((d.homeNGA == r.homeNGA) & (d.crisis_start < r.crisis_start)).any()) d["quasi"] = d.centralization.eq("quasi-polity") 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() for c in ["P_main", "D_pc1", "K_mean", "K_irt", "complexity_pc1", "log_pop"]: d[c + "_z"] = (d[c] - d[c].mean()) / d[c].std() return d def wald_sum_zero(res, a: str, b: str) -> dict: """Wald test of beta_a + beta_b = 0 (the restriction imposed by P = z(D) - z(K)).""" import scipy.stats as st 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] w = s_ ** 2 / var; return {"sum_coef": float(s_), "sum_se": float(np.sqrt(var)), "wald_p": float(1 - st.chi2.cdf(w, 1))} def primary_models(d: pd.DataFrame): rows = [] # common estimation sample for every specification that uses covariates/baselines (A-006) 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"]) specs = { "A1 primary: Y1 ~ P + FE": (f"y1_severity ~ P_main_z + {FE}", "P_main_z"), "A2: Y1 ~ P + X + FE": (f"y1_severity ~ P_main_z + {XS} + {FE}", "P_main_z"), "A3: Y1 ~ P + complexity PC1 + X + FE": (f"y1_severity ~ P_main_z + complexity_pc1_z + {XS} + {FE}", "P_main_z"), "B1 baseline D only": (f"y1_severity ~ D_pc1_z + {XS} + {FE}", "D_pc1_z"), "B2 baseline K only": (f"y1_severity ~ K_mean_z + {XS} + {FE}", "K_mean_z"), "B3 baseline complexity PC1": (f"y1_severity ~ complexity_pc1_z + {XS} + {FE}", "complexity_pc1_z"), "B4 baseline persistence": (f"y1_severity ~ prior_crisis_nga + {XS} + {FE}", "prior_crisis_nga"), "B5 D and K separately": (f"y1_severity ~ D_pc1_z + K_mean_z + {XS} + {FE}", "D_pc1_z"), "Placebo: P built with production infrastructure instead of K (X = warfare only)": (f"y1_severity ~ P_placebo + X_warfare + {FE}", "P_placebo"), "Placebo comparison: P with the same controls (X = warfare only)": (f"y1_severity ~ P_main_z + X_warfare + {FE}", "P_main_z"), } for label, (f, term) in specs.items(): sample = d if label.startswith("A1") else common r = fit_formula(f, sample, BIN, groups="seshat_region"); rows.append(coef_row(r, term, label)) if "B5" in label: rows.append(coef_row(r, "K_mean_z", label)); wt = wald_sum_zero(r, "D_pc1_z", "K_mean_z") 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}) # binary secondary # binary secondary outcome: several region/era cells are all-severe (quasi-separation) -> region FE only, HC1 SE for label, f in {"A1b: Y1b ~ P + region FE (logit, HC1; era FE dropped: separation)": "y1b_severe ~ P_main_z + C(macro_region)", "A2b: Y1b ~ P + X + region FE (logit, HC1)": f"y1b_severe ~ P_main_z + {XS} + C(macro_region)", "A3b: Y1b ~ P (logit, HC1, no FE)": "y1b_severe ~ P_main_z"}.items(): r = fit_formula(f, d, BIN, groups=None); rows.append(coef_row(r, "P_main_z", label)) # exploratory: reforms 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")) return pd.DataFrame(rows) def robustness(d: pd.DataFrame): rows = [] 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", "P_main_lo": "D at lower bounds", "P_main_hi": "D at upper bounds", "P_main_noinf": "inferred codes dropped", "P_within_era": "standardized within era", "P_within_region": "standardized within macro-region", "P_rank": "rank-based P", "P_nosh": "D without settlement hierarchy", "P_noadmlev": "K without administrative levels", "P_nosh_noadmlev": "both hierarchy indicators removed"} for v, lab in variants.items(): if v not in d.columns: # crisis-level file lacks some polity-level variants → merge from polity file continue dd = d.copy(); dd["_P"] = (dd[v] - dd[v].mean()) / dd[v].std() r = fit_formula(f"y1_severity ~ _P + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "_P", f"R: {lab}")) 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")) 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")) 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")) 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")) 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")) 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")) 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")) 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")) return pd.DataFrame(rows) def out_of_sample(d: pd.DataFrame): 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}", "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}", "B5 D+K": f"y1_severity ~ D_pc1_z + K_mean_z + {XS} + {FE}"} # era FE with leave-one-region-out can leave unseen levels; macro-region FE cannot be used when the region is held out f_loro = {k: v.replace("C(macro_region) + ", "") for k, v in formulas.items()} loro = cv_compare(d, f_loro, "y1_severity", BIN, fold_col="macro_region", fractional=True) kf = cv_compare(d, formulas, "y1_severity", BIN, n_repeats=10, k=5, fractional=True, block_col="polity_id") loro["scheme"] = "leave-one-macro-region-out"; kf["scheme"] = "10x5-fold" cv = pd.concat([loro, kf]); cv.to_csv(RES / "h2_seshat_cv_folds.csv", index=False) wins = [] for scheme, c in [("loro", loro), ("kfold", kf)]: for b in ["B0 FE only", "B1 D only", "B2 K only", "B3 complexity", "B4 persistence", "B5 D+K"]: w = fold_win_share(c, "A2 P", b); w["scheme"] = scheme; wins.append(w) wins = pd.DataFrame(wins); wins.to_csv(RES / "h2_seshat_cv_wins.csv", index=False) # binary Y1b out-of-sample metrics fb = {k: v.replace("y1_severity", "y1b_severe") for k, v in formulas.items()} kfb = cv_compare(d, fb, "y1b_severe", BIN, n_repeats=10, k=5, block_col="polity_id") kfb.to_csv(RES / "h2_seshat_cv_folds_binary.csv", index=False) summ = kfb.groupby("model")[["oos_loglik", "auc", "brier"]].mean().reset_index(); summ.to_csv(RES / "h2_seshat_cv_binary_summary.csv", index=False) return wins, summ def h4_nonlinear(d: pd.DataFrame): """Linear vs hinge vs spline, out of sample (10x5-fold); hinge knot chosen on the training fold only.""" rng = np.random.default_rng(20261101) base = f"y1_severity ~ P_main_z + {XS} + {FE}" dd = d.dropna(subset=["y1_severity", "P_main_z", "X_warfare", "X_prod_infra", "macro_region", "era"]).reset_index(drop=True) rows = []; pols = dd.polity_id.unique() from wp35.econometrics.common import map_unseen_levels for rep in range(10): perm = rng.permutation(len(pols)); parts = np.array_split(perm, 5) for i, te in enumerate(parts): m = dd.polity_id.isin(pols[te]).values; tr, ts = dd[~m], map_unseen_levels(dd[~m], dd[m], base) try: lin = smf.glm(base, tr, family=BIN).fit(); ll_lin = oos_loglik_binomial(ts.y1_severity, lin.predict(ts)) h = hinge_fit(tr, "y1_severity", base, "P_main_z", BIN) ts["_hinge"] = np.clip(ts.P_main_z - h[1], 0, None); ll_h = oos_loglik_binomial(ts.y1_severity, h[2].predict(ts)) 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)) rows.append({"fold": f"{rep}_{i}", "ll_linear": ll_lin, "ll_hinge": ll_h, "ll_spline": ll_s, "knot": h[1]}) except Exception as e: rows.append({"fold": f"{rep}_{i}", "error": str(e)[:80]}) r = pd.DataFrame(rows); r.to_csv(RES / "h4_seshat_folds.csv", index=False) ok = r.dropna(subset=["ll_linear"]) 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()), "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())} # in-sample hinge for description 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"]) return summ def item_level(d: pd.DataFrame): """Per-item logits (A-002): does P (and D) predict each CrisisDB consequence? Plus dimensionality diagnostics.""" from wp35.harmonization.seshat import INSTABILITY_ITEMS cols = [i + "_v" for i in INSTABILITY_ITEMS] full = pd.read_parquet(PROC / "seshat_crisis.parquet").set_index("crisis_id")[cols] dd = d.set_index("crisis_id").join(full, rsuffix="_full") rows = [] for it in INSTABILITY_ITEMS: y = dd[it + "_v_full"] if it + "_v_full" in dd else dd[it + "_v"] sub = dd.assign(_y=y).dropna(subset=["_y", "P_main_z", "D_pc1_z"]); sub = sub[sub._y.isin([0.0, 1.0])] row = {"item": it, "n": int(len(sub)), "share_present": float(sub._y.mean())} for term in ["P_main_z", "D_pc1_z", "K_mean_z"]: try: 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] except Exception: row[f"coef_{term}"] = np.nan; row[f"se_{term}"] = np.nan rows.append(row) out = pd.DataFrame(rows); out.to_csv(RES / "h2_seshat_items.csv", index=False) M = full.replace({0.5: np.nan}); C = M.corr(min_periods=30); off = C.values[np.triu_indices(len(cols), 1)] ev = np.linalg.eigvalsh(C.fillna(0).values) 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)} (RES / "h2_seshat_items_dimensionality.json").write_text(json.dumps(diag, indent=2)) return out, diag def h5_blocks(d: pd.DataFrame): rows = [] for col, name in [("macro_region", "region"), ("era", "era"), ("cultural_sphere", "sphere")]: for g in sorted(d[col].dropna().unique()): held = d[d[col] != g]; inb = d[d[col] == g] other = "C(era)" if col in ("macro_region", "cultural_sphere") else "C(macro_region)" try: try: r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {other}", held, BIN, groups="seshat_region"); cov = "cluster" except Exception: try: r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {other}", held, BIN, groups=None); cov = "HC1 (cluster cov singular)" except Exception: r = fit_formula(f"y1_severity ~ P_main_z + {XS}", held, BIN, groups="seshat_region"); cov = "cluster, FE dropped (design singular)" 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"] = cov if inb.P_main_z.notna().sum() >= 12: try: 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"] except Exception: try: 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" except Exception as e2: row["within_block_note"] = str(e2)[:60] rows.append(row) except Exception as e: rows.append({"model": f"H5 {name} {g}", "error": str(e)[:80]}) out = pd.DataFrame(rows); out.to_csv(RES / "h5_seshat_blocks.csv", index=False) return out def y2_cox_and_y3(dp: pd.DataFrame): rows = [] d = dp.dropna(subset=["P_main", "duration", "macro_region"]).copy() d["event"] = (~d.censored).astype(int); d["P_z"] = (d.P_main - d.P_main.mean()) / d.P_main.std() for c in ["X_warfare", "X_prod_infra", "X_military", "D_pc1", "K_mean", "complexity_pc1", "n_polities_in_nga"]: d[c + "_z"] = (d[c] - d[c].mean()) / d[c].std() d["era_num"] = d.midpoint 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"], "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"], "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"]} for label, cols in specs.items(): dd = d.dropna(subset=cols + ["seshat_region"])[cols + ["duration", "event", "macro_region", "seshat_region"]].copy() try: cph = CoxPHFitter(); cph.fit(dd, duration_col="duration", event_col="event", strata=["macro_region"], cluster_col="seshat_region", robust=True) 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]: s = cph.summary.loc[t] 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"])}) if "Y2-2" in label: ph = proportional_hazard_test(cph, dd, time_transform="rank"); rows[-1]["schoenfeld_p_P"] = float(ph.summary.loc["P_z", "p"]) except Exception as e: rows.append({"model": label, "error": str(e)[:100]}) # Y3 discontinuous succession meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet")[["polity_id", "y3_discontinuous"]] 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() 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(): try: 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())})) except Exception as e: rows.append({"model": label, "error": str(e)[:100]}) out = pd.DataFrame(rows); out.to_csv(RES / "h2_seshat_y2_y3.csv", index=False) return out def main(): d = load_crisis(); dp = pd.read_parquet(PROC / "seshat_indices.parquet") d.to_parquet(PROC / "seshat_h2_sample.parquet", index=False) prim = primary_models(d); prim.to_csv(RES / "h2_seshat_primary.csv", index=False) rob = robustness(d); rob.to_csv(RES / "h2_seshat_robustness.csv", index=False) wins, summ = out_of_sample(d) h4 = h4_nonlinear(d); (RES / "h4_seshat_summary.json").write_text(json.dumps(h4, indent=2)) h5 = h5_blocks(d) items, idiag = item_level(d) y23 = y2_cox_and_y3(dp) print("\n=== ITEM LEVEL ===\n", items.round(3).to_string(index=False), "\n", idiag) print("corr(K_mean, polities per NGA) =", round(float(dp.K_mean.corr(dp.n_polities_in_nga)), 3)) pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 70) cols = ["model", "term", "coef", "se", "p", "n"] print("=== PRIMARY (Y1 fractional logit, cluster region) ===\n", prim[cols].round(3).to_string(index=False)) print("\n=== ROBUSTNESS ===\n", rob[cols].round(3).to_string(index=False)) print("\n=== OOS wins of A2 over baselines ===\n", wins.round(3).to_string(index=False)) print("\n=== Y1b OOS summary ===\n", summ.round(3).to_string(index=False)) print("\n=== H4 ===", json.dumps(h4, indent=1)) 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)) 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)) if __name__ == "__main__": main()