1"""Phase 3 measurement: build D, K, P (and rival covariates) exactly as pre-registered, with alternatives,2measurement error, cross-region comparability diagnostics and an a-priori power computation.3Usage: python -m wp35.measurement.build_indices4No outcome (Y1, Y2, Y3, Y4-6) is read or used anywhere in this module.5"""6from __future__ import annotations7import json8from pathlib import Path9import numpy as np10import pandas as pd11import matplotlib12matplotlib.use("Agg")13import matplotlib.pyplot as plt14import statsmodels.api as sm15from wp35.measurement.irt import fit_2pl16from wp35.harmonization.build_panels import K_BLOCKS, STANDARDS, D_VARS, PROD_INFRA, MILITARY1718ROOT = Path(__file__).resolve().parents[3]19PROC = ROOT / "data" / "processed"; RES = ROOT / "results" / "empirical"; FIG = ROOT / "paper" / "figures"20RES.mkdir(parents=True, exist_ok=True); FIG.mkdir(parents=True, exist_ok=True)21SEED = 2026110122K_ITEMS_MODERN = ["v2clrspct", "v2stcritrecadm", "v2cltrnslw", "v2stfisccap", "v2svstterr", "v2elembcap"]23BIN_K = [v for b in K_BLOCKS.values() for v in b if v != "administrative_level"]242526def z(x: pd.Series, ref: pd.Series | None = None) -> pd.Series:27 ref = x if ref is None else ref28 return (x - ref.mean()) / ref.std(ddof=0)293031def pc1_pairwise(df: pd.DataFrame, min_coded: int = 2):32 """First principal component of the pairwise-complete correlation matrix; scores = loading-weighted mean of the33 available standardized variables (requires >= min_coded coded). Loadings signed so that the sum is positive."""34 Zs = (df - df.mean()) / df.std(ddof=0)35 C = Zs.corr(min_periods=10).fillna(0).values36 w, V = np.linalg.eigh(C); v = V[:, -1]; v = v if v.sum() > 0 else -v37 load = pd.Series(v / np.abs(v).sum(), index=df.columns)38 avail = Zs.notna()39 num = (Zs.fillna(0) * load).sum(axis=1); den = (avail * load.abs()).sum(axis=1)40 score = num / den.where(den > 0)41 score[avail.sum(axis=1) < min_coded] = np.nan42 return score, load, float(w[-1] / w.sum())434445def block_score(wide: pd.DataFrame, items: list[str], min_share: float = 0.5) -> pd.Series:46 sub = wide.reindex(columns=items).copy()47 if "administrative_level" in sub.columns:48 r = sub["administrative_level"].rank(method="average")49 sub["administrative_level"] = (r - 1) / (r.max() - 1)50 coded = sub.notna().mean(axis=1)51 return sub.mean(axis=1).where(coded >= min_share)525354def build_seshat():55 wide = pd.read_parquet(PROC / "seshat_polity_wide.parquet")56 agg = pd.read_parquet(PROC / "seshat_polity_agg.parquet")57 meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet").set_index("polity_id")58 logw = agg.pivot(index="polity_id", columns="variable", values="logv_mean")59 lo = agg.pivot(index="polity_id", columns="variable", values="lo_min"); hi = agg.pivot(index="polity_id", columns="variable", values="hi_max")60 out = pd.DataFrame(index=wide.index)61 # ---- D ----62 Dmat = pd.DataFrame({"log_pop": logw["polity_population"], "log_terr": logw["polity_territory"],63 "log_cap": logw["population_of_the_largest_settlement"], "settl_hier": wide["settlement_hierarchy"]})64 out["D_pc1"], D_load, D_var = pc1_pairwise(Dmat)65 out["log_pop"] = Dmat.log_pop; out["log_terr"] = Dmat.log_terr; out["log_cap"] = Dmat.log_cap; out["settl_hier"] = Dmat.settl_hier66 # bounds for measurement-error sensitivity (geometric variables only)67 with np.errstate(divide="ignore"):68 out["log_pop_lo"] = np.log(lo["polity_population"].where(lo["polity_population"] > 0)); out["log_pop_hi"] = np.log(hi["polity_population"].where(hi["polity_population"] > 0))69 out["log_terr_lo"] = np.log(lo["polity_territory"].where(lo["polity_territory"] > 0)); out["log_terr_hi"] = np.log(hi["polity_territory"].where(hi["polity_territory"] > 0))70 # ---- K blocks ----71 for b, items in K_BLOCKS.items():72 out[f"K_{b}"] = block_score(wide, items)73 blocks = out[[f"K_{b}" for b in K_BLOCKS]]74 out["K_mean"] = blocks.mean(axis=1).where(blocks.notna().sum(axis=1) >= 3)75 out["K_n_blocks"] = blocks.notna().sum(axis=1)76 out["K_pc1"], K_load, K_var = pc1_pairwise(blocks, min_coded=3)77 # IRT on binary items (disputed 0.5 -> missing)78 Xb = wide.reindex(columns=BIN_K).copy(); Xb = Xb.where(Xb != 0.5)79 irt = fit_2pl(Xb.values)80 out["K_irt"] = irt["eap"]; out["K_irt_psd"] = irt["psd"]; out["K_irt_nitems"] = irt["n_items"]81 item_params = pd.DataFrame({"item": BIN_K, "a": irt["a"], "b": irt["b"], "block": [next(k for k, v in K_BLOCKS.items() if i in v) for i in BIN_K]})82 item_params.to_csv(RES / "measurement_irt_item_parameters.csv", index=False)83 # standards sub-index (sparse; latent only)84 out["K_standards"] = block_score(wide, STANDARDS)85 # ---- P ----86 for dname in ["D_pc1", "log_pop"]:87 for kname in ["K_mean", "K_irt", "K_pc1"]:88 out[f"P_{dname}_{kname}"] = z(out[dname]) - z(out[kname])89 out["P_main"] = out["P_D_pc1_K_mean"]90 # standardization windows (sensitivity)91 df = out.join(meta[["macro_region", "era"]])92 out["P_within_era"] = df.groupby("era").D_pc1.transform(z) - df.groupby("era").K_mean.transform(z)93 out["P_within_region"] = df.groupby("macro_region").D_pc1.transform(z) - df.groupby("macro_region").K_mean.transform(z)94 out["P_rank"] = out.D_pc1.rank(pct=True) - out.K_mean.rank(pct=True)95 # bounds96 Dlo = pd.DataFrame({"log_pop": out.log_pop_lo, "log_terr": out.log_terr_lo, "log_cap": Dmat.log_cap, "settl_hier": Dmat.settl_hier})97 Dhi = pd.DataFrame({"log_pop": out.log_pop_hi, "log_terr": out.log_terr_hi, "log_cap": Dmat.log_cap, "settl_hier": Dmat.settl_hier})98 out["D_pc1_lo"] = ((Dlo - Dmat.mean()) / Dmat.std(ddof=0)).fillna(0).mul(D_load).sum(axis=1) / (Dlo.notna() * D_load.abs()).sum(axis=1).replace(0, np.nan)99 out["D_pc1_hi"] = ((Dhi - Dmat.mean()) / Dmat.std(ddof=0)).fillna(0).mul(D_load).sum(axis=1) / (Dhi.notna() * D_load.abs()).sum(axis=1).replace(0, np.nan)100 out.loc[out.D_pc1.isna(), ["D_pc1_lo", "D_pc1_hi"]] = np.nan101 out["P_main_lo"] = z(out.D_pc1_lo, out.D_pc1) - z(out.K_mean); out["P_main_hi"] = z(out.D_pc1_hi, out.D_pc1) - z(out.K_mean)102 # drop inferred codes sensitivity: recompute K_mean using only non-inferred facts103 facts = pd.read_parquet(PROC / "seshat_facts.parquet")104 ni = facts[(~facts.inferred) & facts.value.notna()].groupby(["polity_id", "variable"]).value.mean().unstack()105 out["K_mean_noinf"] = pd.concat([block_score(ni, it) for it in K_BLOCKS.values()], axis=1).pipe(lambda b: b.mean(axis=1).where(b.notna().sum(axis=1) >= 3)).reindex(out.index)106 out["P_main_noinf"] = z(out.D_pc1) - z(out.K_mean_noinf)107 # ---- rival covariates (excluded from K by design) ----108 wf = facts[facts.sheet == "Warfare"].groupby(["polity_id", "variable"]).value.mean().unstack()109 out["X_warfare"] = wf.mean(axis=1).where(wf.notna().mean(axis=1) >= 0.5).reindex(out.index)110 out["X_prod_infra"] = block_score(wide, PROD_INFRA)111 mil = wide.reindex(columns=MILITARY).copy(); r = mil["military_level"].rank(); mil["military_level"] = (r - 1) / (r.max() - 1)112 out["X_military"] = mil.mean(axis=1).where(mil.notna().mean(axis=1) >= 0.5)113 out["X_moralizing"] = wide.get("moralizing_supernatural_punishment_and_reward")114 rt = facts[facts.sheet == "Religion"].groupby(["polity_id", "variable"]).value.mean().unstack()115 out["X_moralizing"] = rt.get("moralizing_supernatural_punishment_and_reward").reindex(out.index) if "moralizing_supernatural_punishment_and_reward" in rt else np.nan116 # Turchin-style single complexity axis over ALL social-complexity indicators (baseline B3)117 sc_vars = [v for v in facts[facts.sheet == "Social complexity"].variable.unique() if v in wide.columns]118 scmat = wide[sc_vars].copy()119 for v in ["polity_population", "polity_territory", "population_of_the_largest_settlement"]:120 scmat[v] = logw[v]121 scmat = scmat.loc[:, scmat.notna().sum() >= 100]122 out["complexity_pc1"], C_load, C_var = pc1_pairwise(scmat, min_coded=5)123 out = out.join(meta[["seshat_region", "homeNGA", "macro_region", "era", "start_year", "end_year", "midpoint", "duration", "censored", "centralization", "utm_zone"]])124 out.index.name = "polity_id"125 out.reset_index().to_parquet(PROC / "seshat_indices.parquet", index=False)126 # ---- pre-crisis versions (same loadings/standardization as the full sample) ----127 pre = pd.read_parquet(PROC / "seshat_precrisis_agg.parquet")128 pw = pre.pivot(index="crisis_id", columns="variable", values="value_mean"); plog = pre.pivot(index="crisis_id", columns="variable", values="logv_mean")129 cr = pd.DataFrame(index=pw.index)130 Dp = pd.DataFrame({"log_pop": plog.get("polity_population"), "log_terr": plog.get("polity_territory"), "log_cap": plog.get("population_of_the_largest_settlement"), "settl_hier": pw.get("settlement_hierarchy")}, index=pw.index)131 Zp = (Dp - Dmat.mean()) / Dmat.std(ddof=0)132 cr["D_pc1"] = (Zp.fillna(0) * D_load).sum(axis=1) / (Zp.notna() * D_load.abs()).sum(axis=1).replace(0, np.nan)133 cr.loc[Zp.notna().sum(axis=1) < 2, "D_pc1"] = np.nan134 cr["log_pop"] = Dp.log_pop135 for b, items in K_BLOCKS.items():136 cr[f"K_{b}"] = block_score(pw, items)137 bl = cr[[f"K_{b}" for b in K_BLOCKS]]; cr["K_mean"] = bl.mean(axis=1).where(bl.notna().sum(axis=1) >= 3)138 Xp = pw.reindex(columns=BIN_K).where(lambda d: d != 0.5)139 # EAP with the full-sample item parameters (no refit): posterior under fixed a, b140 nodes, weights = np.polynomial.hermite_e.hermegauss(41); weights /= weights.sum()141 from scipy.special import expit142 Pq = expit(irt["a"][None, :] * (nodes[:, None] - irt["b"][None, :]))143 Xv = Xp.values; obs = ~np.isnan(Xv); Xz_ = np.where(obs, Xv, 0)144 L = Xz_ @ np.log(Pq + 1e-300).T + ((1 - Xz_) * obs) @ np.log(1 - Pq + 1e-300).T + np.log(weights)[None, :]145 post = np.exp(L - L.max(axis=1, keepdims=True)); post /= post.sum(axis=1, keepdims=True)146 cr["K_irt"] = np.where(obs.sum(axis=1) >= 3, post @ nodes, np.nan)147 cr["K_pc1"] = (((bl - blocks.mean()) / blocks.std(ddof=0)).fillna(0) * K_load).sum(axis=1) / (bl.notna() * K_load.abs()).sum(axis=1).replace(0, np.nan)148 cr.loc[bl.notna().sum(axis=1) < 3, "K_pc1"] = np.nan149 for dname in ["D_pc1", "log_pop"]:150 for kname in ["K_mean", "K_irt", "K_pc1"]:151 cr[f"P_{dname}_{kname}"] = z(cr[dname], out[dname]) - z(cr[kname], out[kname])152 cr["P_main"] = cr["P_D_pc1_K_mean"]153 wfp = pre[pre.variable.isin(wf.columns)].pivot(index="crisis_id", columns="variable", values="value_mean")154 cr["X_warfare"] = wfp.mean(axis=1).where(wfp.notna().mean(axis=1) >= 0.5).reindex(cr.index)155 cr["X_prod_infra"] = block_score(pw, PROD_INFRA)156 milp = pw.reindex(columns=MILITARY).copy(); milp["military_level"] = milp["military_level"].rank().pipe(lambda r: (r - 1) / (r.max() - 1)); cr["X_military"] = milp.mean(axis=1).where(milp.notna().mean(axis=1) >= 0.5)157 scp = pw.reindex(columns=scmat.columns).copy()158 for v in ["polity_population", "polity_territory", "population_of_the_largest_settlement"]:159 scp[v] = plog.get(v)160 Zc = (scp - scmat.mean()) / scmat.std(ddof=0)161 cr["complexity_pc1"] = (Zc.fillna(0) * C_load).sum(axis=1) / (Zc.notna() * C_load.abs()).sum(axis=1).replace(0, np.nan)162 cr.loc[Zc.notna().sum(axis=1) < 5, "complexity_pc1"] = np.nan163 crisis = pd.read_parquet(PROC / "seshat_crisis.parquet").set_index("crisis_id")164 cr["polity_id"] = crisis.polity_id; cr["crisis_start"] = crisis.crisis_start165 cr = cr.join(meta[["macro_region", "era", "seshat_region", "homeNGA", "start_year"]], on="polity_id")166 cr["years_since_start"] = cr.crisis_start - cr.start_year167 cr.reset_index().to_parquet(PROC / "seshat_crisis_indices.parquet", index=False)168 # ---- diagnostics: collinearity, comparability, measurement error ----169 diag = {"D_pc1_loadings": D_load.round(4).to_dict(), "D_pc1_var_share": round(D_var, 4),170 "K_pc1_loadings": K_load.round(4).to_dict(), "K_pc1_var_share": round(K_var, 4),171 "complexity_pc1_var_share": round(C_var, 4), "complexity_pc1_n_vars": int(scmat.shape[1]),172 "irt_loglik": irt["loglik"], "irt_iterations": irt["iterations"],173 "n_D_pc1": int(out.D_pc1.notna().sum()), "n_K_mean": int(out.K_mean.notna().sum()), "n_K_irt": int(out.K_irt.notna().sum()),174 "n_P_main": int(out.P_main.notna().sum()), "n_P_main_crisis": int(cr.P_main.notna().sum()),175 "corr_zD_zK": {"D_pc1~K_mean": round(out.D_pc1.corr(out.K_mean), 4), "D_pc1~K_irt": round(out.D_pc1.corr(out.K_irt), 4),176 "log_pop~K_mean": round(out.log_pop.corr(out.K_mean), 4), "K_mean~K_irt": round(out.K_mean.corr(out.K_irt), 4),177 "K_mean~K_pc1": round(out.K_mean.corr(out.K_pc1), 4), "P_main~complexity_pc1": round(out.P_main.corr(out.complexity_pc1), 4),178 "P_main~D_pc1": round(out.P_main.corr(out.D_pc1), 4), "P_main~K_mean": round(out.P_main.corr(out.K_mean), 4)},179 "P_variance_share_of_D_plus_K": round(float(out.P_main.var() / (z(out.D_pc1).var() + z(out.K_mean).var())), 4),180 "block_correlations": blocks.corr().round(3).to_dict(),181 "block_alpha_by_macro_region": {}, "block_corr_mean_by_era": {},182 "measurement_error": {"D_pc1_halfwidth_median": float(((out.D_pc1_hi - out.D_pc1_lo) / 2).median()),183 "D_pc1_halfwidth_q90": float(((out.D_pc1_hi - out.D_pc1_lo) / 2).quantile(0.9)),184 "K_irt_psd_median": float(out.K_irt_psd.median()), "K_irt_psd_q90": float(out.K_irt_psd.quantile(0.9)),185 "sd_K_irt": float(out.K_irt.std()), "P_main_range_median": float((out.P_main_hi - out.P_main_lo).abs().median())}}186 for r, g in df.groupby("macro_region"):187 b = g[[f"K_{k}" for k in K_BLOCKS]].dropna()188 if len(b) >= 8:189 k = b.shape[1]; alpha = k / (k - 1) * (1 - b.var(ddof=1).sum() / b.sum(axis=1).var(ddof=1))190 diag["block_alpha_by_macro_region"][r] = {"n": int(len(b)), "alpha": round(float(alpha), 3)}191 for e, g in df.groupby("era"):192 b = g[[f"K_{k}" for k in K_BLOCKS]]193 c = b.corr().values; diag["block_corr_mean_by_era"][e] = {"n": int(b.dropna().shape[0]), "mean_offdiag_corr": round(float(c[np.triu_indices(4, 1)].mean()), 3)}194 # IRT invariance: refit on Old World vs Americas/Oceania and compare difficulties195 new_world = df.macro_region.isin(["North_America", "South_America", "Southeast_Asia_Oceania"])196 inv = {}197 for name, mask in [("old_world", ~new_world), ("new_world_oceania", new_world)]:198 f = fit_2pl(Xb[mask.values].values); inv[name] = {"n": int(mask.sum()), "b": f["b"].round(3).tolist(), "a": f["a"].round(3).tolist()}199 bo, bn = np.array(inv["old_world"]["b"]), np.array(inv["new_world_oceania"]["b"])200 ok = np.isfinite(bo) & np.isfinite(bn) & (np.abs(bo) < 5.9) & (np.abs(bn) < 5.9)201 diag["irt_invariance"] = {"items_compared": int(ok.sum()), "corr_difficulty_old_vs_new": round(float(np.corrcoef(bo[ok], bn[ok])[0, 1]), 3),202 "mean_abs_diff_difficulty": round(float(np.abs(bo[ok] - bn[ok]).mean()), 3), "items": BIN_K, "fits": inv}203 # leave-one-block-out sensitivity of K and P204 loo = {}205 for b in K_BLOCKS:206 others = [f"K_{k}" for k in K_BLOCKS if k != b]; kk = out[others].mean(axis=1).where(out[others].notna().sum(axis=1) >= 2)207 loo[b] = {"corr_K_with_full": round(float(kk.corr(out.K_mean)), 4), "corr_P_with_main": round(float((z(out.D_pc1) - z(kk)).corr(out.P_main)), 4)}208 diag["leave_one_block_out"] = loo209 (RES / "measurement_diagnostics.json").write_text(json.dumps(diag, indent=2))210 # ---- figures ----211 fig, axes = plt.subplots(1, 3, figsize=(12, 3.8))212 axes[0].scatter(out.K_mean, out.D_pc1, s=8, alpha=0.5, c="#444"); axes[0].set_xlabel("K (mean of 4 blocks)"); axes[0].set_ylabel("D (PC1 of log scale variables)"); axes[0].set_title(f"r = {out.D_pc1.corr(out.K_mean):.2f}, n = {int((out.D_pc1.notna() & out.K_mean.notna()).sum())}", fontsize=9)213 axes[1].scatter(out.K_mean, out.K_irt, s=8, alpha=0.5, c="#444"); axes[1].set_xlabel("K mean"); axes[1].set_ylabel("K IRT (2PL EAP)"); axes[1].set_title(f"r = {out.K_mean.corr(out.K_irt):.2f}", fontsize=9)214 axes[2].errorbar(item_params.b, range(len(item_params)), xerr=None, fmt="o", ms=3, color="#444"); axes[2].set_yticks(range(len(item_params))); axes[2].set_yticklabels(item_params.item, fontsize=5); axes[2].set_xlabel("IRT difficulty b"); axes[2].set_title("Item difficulties", fontsize=9)215 fig.tight_layout(); fig.savefig(FIG / "fig_measurement_DK.pdf"); plt.close(fig)216 fig, ax = plt.subplots(figsize=(4.5, 4.2)); ax.scatter(bo[ok], bn[ok], s=12, c="#444"); lim = [min(bo[ok].min(), bn[ok].min()) - .3, max(bo[ok].max(), bn[ok].max()) + .3]; ax.plot(lim, lim, "k--", lw=.8)217 ax.set_xlabel("difficulty, Old World fit"); ax.set_ylabel("difficulty, Americas/Oceania fit"); ax.set_title(f"IRT item invariance (r = {diag['irt_invariance']['corr_difficulty_old_vs_new']:.2f})", fontsize=9)218 fig.tight_layout(); fig.savefig(FIG / "fig_irt_invariance.pdf"); plt.close(fig)219 return out, cr, diag220221222def build_modern():223 p = pd.read_parquet(PROC / "modern_panel.parquet")224 items = p[K_ITEMS_MODERN]225 Zi = (items - items.mean()) / items.std(ddof=0)226 p["K_modern"] = Zi.mean(axis=1).where(Zi.notna().sum(axis=1) >= 4)227 p["K_modern_pc1"], Km_load, Km_var = pc1_pairwise(items, min_coded=4)228 p["D_modern"] = p.log_pop229 p["P_modern"] = z(p.D_modern) - z(p.K_modern)230 p["P_modern_pc1"] = z(p.D_modern) - z(p.K_modern_pc1)231 p["P_modern_within_year"] = p.groupby("year").D_modern.transform(z) - p.groupby("year").K_modern.transform(z)232 p["P_modern_area"] = z(((z(p.log_pop) + z(p.log_area)) / 2)) - z(p.K_modern)233 p.to_parquet(PROC / "modern_indices.parquet", index=False)234 diag = {"K_modern_pc1_loadings": Km_load.round(4).to_dict(), "K_modern_pc1_var_share": round(Km_var, 4),235 "n_P_modern": int(p.P_modern.notna().sum()), "corr_D_K_modern": round(float(p.D_modern.corr(p.K_modern)), 4),236 "corr_K_modern_polyarchy": round(float(p.K_modern.corr(p.v2x_polyarchy)), 4),237 "item_corr": items.corr().round(3).to_dict(),238 "corr_K_by_region": {int(r): round(float(g.D_modern.corr(g.K_modern)), 3) for r, g in p.groupby("e_regionpol_6C") if g.P_modern.notna().sum() > 50}}239 (RES / "measurement_diagnostics_modern.json").write_text(json.dumps(diag, indent=2))240 return p, diag241242243def power_h2a(cr: pd.DataFrame, n_rep: int = 400):244 """A-priori MDE for H2-A: random outcomes assigned to the actual P distribution of crisis records (no real Y used)."""245 rng = np.random.default_rng(SEED)246 d = cr.dropna(subset=["P_main"]).copy()247 P = z(d.P_main).values; reg = pd.factorize(d.seshat_region)[0]; n = len(P)248 base = 0.77 # Y1b positive share read from coverage_summary (122/158); used only to set the intercept249 res = []250 for beta in [0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0]:251 rej = 0252 for _ in range(n_rep):253 a0 = np.log(base / (1 - base))254 y = rng.random(n) < 1 / (1 + np.exp(-(a0 + beta * P)))255 try:256 m = sm.Logit(y.astype(float), sm.add_constant(P)).fit(disp=0, cov_type="cluster", cov_kwds={"groups": reg})257 rej += m.pvalues[1] < 0.05258 except Exception:259 pass260 res.append({"beta": beta, "power": rej / n_rep, "n": n})261 res = pd.DataFrame(res); res.to_csv(RES / "power_h2a_binary.csv", index=False)262 # fractional outcome: Binomial(14, logistic(a0 + beta P)) with base severity 0.29263 res2 = []264 for beta in [0.1, 0.15, 0.2, 0.3, 0.4, 0.5]:265 rej = 0266 for _ in range(n_rep):267 a0 = np.log(0.29 / 0.71); pr = 1 / (1 + np.exp(-(a0 + beta * P))); y = rng.binomial(14, pr) / 14268 try:269 m = sm.GLM(y, sm.add_constant(P), family=sm.families.Binomial()).fit(cov_type="cluster", cov_kwds={"groups": reg})270 rej += m.pvalues[1] < 0.05271 except Exception:272 pass273 res2.append({"beta": beta, "power": rej / n_rep, "n": n})274 res2 = pd.DataFrame(res2); res2.to_csv(RES / "power_h2a_fractional.csv", index=False)275 return res, res2276277278def main():279 out, cr, diag = build_seshat()280 p, dm = build_modern()281 pw1, pw2 = power_h2a(cr)282 print(json.dumps({k: v for k, v in diag.items() if k in ("D_pc1_loadings", "D_pc1_var_share", "K_pc1_loadings", "n_D_pc1", "n_K_mean", "n_K_irt", "n_P_main", "n_P_main_crisis", "corr_zD_zK", "P_variance_share_of_D_plus_K", "measurement_error", "irt_invariance", "leave_one_block_out", "block_alpha_by_macro_region")}, indent=1, default=str)[:6000])283 print(json.dumps({k: v for k, v in dm.items() if k != "item_corr"}, indent=1))284 print("power binary:\n", pw1.to_string(index=False)); print("power fractional:\n", pw2.to_string(index=False))285286287if __name__ == "__main__":288 main()289