1"""H1 — association between coordination capacity and social scale (no direction claimed).2Seshat (polity), Ethnographic Atlas (society), modern panel (country-year). Usage: python -m wp35.econometrics.h13"""4from __future__ import annotations5import json6from pathlib import Path7import numpy as np8import pandas as pd9import statsmodels.api as sm10import statsmodels.formula.api as smf11from statsmodels.miscmodels.ordinal_model import OrderedModel12from scipy.stats import norm13from wp35.econometrics.common import fit_formula, coef_row, conley_ols1415ROOT = Path(__file__).resolve().parents[3]16PROC = ROOT / "data" / "processed"; RES = ROOT / "results" / "empirical"17RES.mkdir(parents=True, exist_ok=True)181920def seshat():21 d = pd.read_parquet(PROC / "seshat_indices.parquet")22 d["quasi"] = d.centralization.eq("quasi-polity")23 rows = []24 base = "C(macro_region) + C(era)"25 specs = {26 "H1-S1 primary: D_pc1 ~ K_mean + FE": ("D_pc1 ~ K_mean + " + base, "K_mean"),27 "H1-S2: D_pc1 ~ K_mean + X + FE": ("D_pc1 ~ K_mean + X_warfare + X_prod_infra + X_military + " + base, "K_mean"),28 "H1-S3: log_pop ~ K_mean + X + FE": ("log_pop ~ K_mean + X_warfare + X_prod_infra + X_military + " + base, "K_mean"),29 "H1-S4: log_terr ~ K_mean + X + FE": ("log_terr ~ K_mean + X_warfare + X_prod_infra + X_military + " + base, "K_mean"),30 "H1-S5: D_pc1 ~ K_irt + X + FE": ("D_pc1 ~ K_irt + X_warfare + X_prod_infra + X_military + " + base, "K_irt"),31 "H1-S6: D_pc1 ~ K blocks + X + FE": ("D_pc1 ~ K_admin + K_rules_money + K_comm_infra + K_records + X_warfare + X_prod_infra + X_military + " + base, "K_admin"),32 "H1-S7: D_pc1 ~ K_mean + X + FE, no quasi-polities": ("D_pc1 ~ K_mean + X_warfare + X_prod_infra + X_military + " + base, "K_mean"),33 }34 for label, (f, term) in specs.items():35 dd = d[~d.quasi] if "no quasi" in label else d36 r = fit_formula(f, dd, None, groups="seshat_region")37 rows.append(coef_row(r, term, label, {"r2": r.rsquared}))38 if "K blocks" in label:39 for t in ["K_rules_money", "K_comm_infra", "K_records"]:40 rows.append(coef_row(r, t, label, {"r2": r.rsquared}))41 # region random effects (MixedLM) on the primary spec42 dd = d.dropna(subset=["D_pc1", "K_mean", "macro_region", "era", "seshat_region"])43 try:44 mm = smf.mixedlm("D_pc1 ~ K_mean + C(era)", dd, groups=dd["seshat_region"]).fit(reml=True)45 rows.append({"model": "H1-S8: D_pc1 ~ K_mean + era, region random intercept", "term": "K_mean", "coef": mm.params["K_mean"], "se": mm.bse["K_mean"],46 "ci_lo": mm.params["K_mean"] - 1.96 * mm.bse["K_mean"], "ci_hi": mm.params["K_mean"] + 1.96 * mm.bse["K_mean"], "p": mm.pvalues["K_mean"], "n": len(dd), "llf": mm.llf})47 except Exception as e:48 rows.append({"model": "H1-S8 mixedlm", "error": str(e)[:100]})49 # standardized effect: K_mean has range 0-1; report SD-scaled coefficient50 out = pd.DataFrame(rows); out["coef_per_sd_K"] = out.coef * d.K_mean.std()51 out.to_csv(RES / "h1_seshat.csv", index=False)52 return out535455def ea():56 d = pd.read_parquet(PROC / "dplace_ea_societies.parquet")57 d["hier"] = d.jurisdictional_hierarchy_beyond; d["cls"] = d.EA066; d["agri"] = d.EA005; d["intens"] = d.EA028; d["subs"] = d.EA04258 rows = []59 # (a) OLS of log community size on hierarchy levels with Conley SE and family-cluster SE60 f = "log_community_size ~ hier + C(subs) + intens + cls + C(region)"61 for cut in [500, 1000]:62 res, se, n = conley_ols(f, d, "Latitude", "Longitude", cut)63 rows.append({"model": f"H1-E1 OLS log community size ~ hierarchy + X + region, Conley {cut} km", "term": "hier", "coef": res.params["hier"], "se": se["hier"],64 "ci_lo": res.params["hier"] - 1.96 * se["hier"], "ci_hi": res.params["hier"] + 1.96 * se["hier"], "p": 2 * (1 - norm.cdf(abs(res.params["hier"] / se["hier"]))), "n": n, "r2": res.rsquared})65 r = fit_formula(f, d, None, groups="Family_ID"); rows.append(coef_row(r, "hier", "H1-E2 same, SE clustered by language family", {"r2": r.rsquared}))66 f2 = "log_population ~ hier + C(subs) + intens + cls + C(region)"67 r = fit_formula(f2, d, None, groups="Family_ID"); rows.append(coef_row(r, "hier", "H1-E3 log population (EA202) ~ hierarchy + X + region, family-clustered", {"r2": r.rsquared}))68 # (b) ordered logit of hierarchy on log community size69 dd = d.dropna(subset=["hier", "log_community_size", "intens", "cls", "region", "Family_ID"]).copy()70 X = pd.get_dummies(dd[["log_community_size", "intens", "cls", "region"]], columns=["region"], drop_first=True, dtype=float)71 try:72 om = OrderedModel(dd.hier.astype(int), X, distr="logit").fit(method="bfgs", disp=0, maxiter=500)73 # cluster-robust SE by family via sandwich on scores74 groups = pd.factorize(dd.Family_ID)[0]75 from statsmodels.stats.sandwich_covariance import cov_cluster76 try:77 V = cov_cluster(om, groups); se = float(np.sqrt(V[0, 0]))78 except Exception:79 se = float(om.bse.iloc[0])80 b = float(om.params.iloc[0])81 rows.append({"model": "H1-E4 ordered logit hierarchy ~ log community size + X + region (family-clustered)", "term": "log_community_size", "coef": b, "se": se, "ci_lo": b - 1.96 * se, "ci_hi": b + 1.96 * se,82 "p": 2 * (1 - norm.cdf(abs(b / se))), "n": len(dd), "llf": om.llf})83 except Exception as e:84 rows.append({"model": "H1-E4 ordered logit", "error": str(e)[:120]})85 out = pd.DataFrame(rows); out.to_csv(RES / "h1_ea.csv", index=False)86 return out878889def modern():90 p = pd.read_parquet(PROC / "modern_indices.parquet")91 d = p.dropna(subset=["log_pop", "K_modern", "log_gdppc", "v2x_polyarchy"]).copy()92 # within-country (demeaned) regression with year FE93 for c in ["log_pop", "K_modern", "log_gdppc", "v2x_polyarchy"]:94 d[c + "_dm"] = d[c] - d.groupby("country_text_id")[c].transform("mean")95 rows = []96 r = fit_formula("log_pop_dm ~ K_modern_dm + log_gdppc_dm + v2x_polyarchy_dm + C(year)", d, None, groups="country_text_id")97 rows.append(coef_row(r, "K_modern_dm", "H1-M1 within-country: log pop ~ K + log gdppc + polyarchy + year FE", {"r2": r.rsquared}))98 r = fit_formula("log_pop ~ K_modern + log_gdppc + v2x_polyarchy + C(year) + C(e_regionpol_6C)", d, None, groups="country_text_id")99 rows.append(coef_row(r, "K_modern", "H1-M2 pooled: log pop ~ K + log gdppc + polyarchy + year FE + region FE", {"r2": r.rsquared}))100 # long differences 50 years101 d50 = d.merge(d[["country_text_id", "year", "log_pop", "K_modern", "log_gdppc", "v2x_polyarchy"]].assign(year=lambda x: x.year + 50), on=["country_text_id", "year"], suffixes=("", "_lag50"))102 for c in ["log_pop", "K_modern", "log_gdppc", "v2x_polyarchy"]:103 d50["d" + c] = d50[c] - d50[c + "_lag50"]104 r = fit_formula("dlog_pop ~ dK_modern + dlog_gdppc + dv2x_polyarchy + C(year)", d50, None, groups="country_text_id")105 rows.append(coef_row(r, "dK_modern", "H1-M3 50-year differences: Δlog pop ~ ΔK + Δlog gdppc + Δpolyarchy + year FE", {"r2": r.rsquared}))106 out = pd.DataFrame(rows); out.to_csv(RES / "h1_modern.csv", index=False)107 return out108109110def main():111 a = seshat(); b = ea(); c = modern()112 pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 90)113 for name, t in [("SESHAT", a), ("EA", b), ("MODERN", c)]:114 print(f"\n=== H1 {name} ===\n", t[[c for c in t.columns if c in ("model", "term", "coef", "se", "p", "n", "r2", "coef_per_sd_K", "error")]].round(4).to_string(index=False))115116117if __name__ == "__main__":118 main()119