1"""Preston curve, F6-dev-01. Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""2import json, os3import numpy as np, pandas as pd4import statsmodels.api as sm5import matplotlib6matplotlib.use("Agg")7import matplotlib.pyplot as plt89SEED = 2026102910rng = np.random.default_rng(SEED)11OUT = "out"12os.makedirs(OUT, exist_ok=True)13kw = dict(keep_default_na=False, na_values=[""]) # 'NA' is not a missing marker14ind = pd.read_csv("data/wb_indicators.csv", **kw)15ent = pd.read_csv("data/wb_entities.csv", **kw)16n_rows = {"wb_indicators.csv": len(ind), "wb_entities.csv": len(ent)}17issues = []1819# ---- audit ----20assert not ind.duplicated(["iso3", "year"]).any()21assert set(ind.iso3) <= set(ent.iso3)22d = ind.merge(ent[["iso3", "name", "region", "income_level"]], on="iso3", suffixes=("", "_meta"))23nm = d.loc[d.name != d.name_meta, "iso3"].unique().tolist()24if nm:25 issues.append({"issue": f"Entity names differ between files for {nm} (trailing whitespace in metadata names; all are aggregates)",26 "handling": "Joined on iso3 only; names ignored."})27n_agg = int(ind.iso3.isin(ent.iso3[ent.region == "Aggregates"]).sum())28issues.append({"issue": f"wb_indicators.csv mixes countries with regional/income aggregates ({n_agg} aggregate rows, {ind.iso3[ind.iso3.isin(ent.iso3[ent.region=='Aggregates'])].nunique()} aggregate entities)",29 "handling": "Excluded all entities with region == 'Aggregates' per the task rule."})30c = d[d.region != "Aggregates"].copy()31c["complete"] = c.gdp_pc_ppp.notna() & c.life_expectancy.notna()32per_year = c.groupby("year").complete.sum()33year = int(per_year[per_year >= 150].index.max())3435# implausible LE volatility (audit, not used for estimation unless it hits the analysis year)36s = c.sort_values(["iso3", "year"]).copy()37s["dle"] = s.groupby("iso3").life_expectancy.diff()38vol = s[s.dle.abs() > 10]39for iso, g in vol.groupby("iso3"):40 ser = s[s.iso3 == iso].set_index("year").life_expectancy.round(1).to_dict()41 issues.append({"issue": f"{iso}: life expectancy jumps by >10 years between consecutive years (ending in {g.year.tolist()}); series 2015-2023 = {ser}",42 "handling": f"Only the {year} cross-section is used; values kept as published (no basis to overwrite official data). Robustness: slope re-estimated excluding all such entities (see out/sensitivity.json)."})4344y = c[c.year == year]45miss = y[~y.complete]46n_excl = int(len(miss))47issues.append({"issue": f"{n_excl} countries/economies lack an indicator in {year} (all missing GDP pc PPP): {', '.join(miss.iso3)}; missingness is not random (small territories/microstates and conflict/sanctioned economies)",48 "handling": "Listwise exclusion, no imputation (task rule); selection bias discussed in report."})49a = y[y.complete].copy()50a["lny"] = np.log(a.gdp_pc_ppp)51assert (a.gdp_pc_ppp > 0).all()5253# ---- estimate ----54X = sm.add_constant(a.lny.values)55le = a.life_expectancy.values56fit = sm.OLS(le, X).fit(cov_type="HC3")57alpha, beta = fit.params58se = fit.bse[1]59ci = fit.conf_int(0.05)[1] # normal-based (z) interval for robust cov60resid = fit.resid61baseline = float(np.sqrt(np.mean((le - le.mean()) ** 2)))62model = float(np.sqrt(np.mean(resid ** 2)))63r = float(np.corrcoef(a.lny, le)[0, 1])64concl = "positive_effect" if ci[0] > 0 else ("negative_effect" if ci[1] < 0 else "no_detectable_effect")6566# ---- sensitivity (reported in report, saved separately) ----67sens = {}68for hc in ["HC0", "HC1", "HC2", "HC3"]:69 sens[f"se_{hc}"] = float(sm.OLS(le, X).fit(cov_type=hc).bse[1])70B = 200071bs = np.empty(B)72n = len(a)73for b in range(B):74 i = rng.integers(0, n, n)75 bs[b] = np.polyfit(a.lny.values[i], le[i], 1)[0]76sens["pairs_bootstrap_ci"] = [float(np.percentile(bs, 2.5)), float(np.percentile(bs, 97.5))]77sub = ~a.iso3.isin(vol.iso3.unique())78sens["beta_excl_volatile"] = float(np.polyfit(a.lny[sub], le[sub], 1)[0])79infl = sm.OLS(le, X).fit().get_influence()80cook = infl.cooks_distance[0]81keep = cook < 4 / n82sens["beta_excl_cooks_gt_4_over_n"] = float(np.polyfit(a.lny[keep], le[keep], 1)[0])83sens["n_cooks_gt_4_over_n"] = int((~keep).sum())84sens["top_influential"] = a.iso3.values[np.argsort(-cook)[:5]].tolist()85q = sm.OLS(le, sm.add_constant(np.column_stack([a.lny, a.lny ** 2]))).fit(cov_type="HC3")86sens["quadratic_coef_lny2"] = float(q.params[2]); sens["quadratic_coef_lny2_p"] = float(q.pvalues[2])87sens["largest_neg_residuals"] = a.assign(e=resid).nsmallest(5, "e")[["iso3", "e"]].round(2).values.tolist()88for name, f in [("beta_year_2019_prepandemic", 2019)]:89 z = c[(c.year == f) & c.complete]90 sens[name] = float(np.polyfit(np.log(z.gdp_pc_ppp), z.life_expectancy, 1)[0])91with open(f"{OUT}/sensitivity.json", "w") as fh:92 json.dump(sens, fh, indent=2)9394# ---- table & figure ----95order = ["Low income", "Lower middle income", "Upper middle income", "High income"]96t = (a.groupby("income_level").agg(n_countries=("iso3", "size"), mean_log_gdp_pc=("lny", "mean"),97 mean_life_expectancy=("life_expectancy", "mean"))98 .reindex([o for o in order if o in set(a.income_level)]).reset_index())99t.round(4).to_csv(f"{OUT}/table1.csv", index=False)100101fig, ax = plt.subplots(figsize=(7, 5))102cols = dict(zip(order, ["#d62728", "#ff7f0e", "#2ca02c", "#1f77b4"]))103for lev in order:104 m = a.income_level == lev105 ax.scatter(a.lny[m], le[m], s=14, alpha=.75, color=cols[lev], label=f"{lev} (n={m.sum()})")106xs = np.linspace(a.lny.min(), a.lny.max(), 100)107ax.plot(xs, alpha + beta * xs, "k-", lw=1.5, label=f"OLS: LE = {alpha:.1f} + {beta:.2f}·ln(GDPpc)")108ax.set_xlabel("ln(GDP per capita, PPP, constant 2021 int'l $)")109ax.set_ylabel("Life expectancy at birth (years)")110ax.set_title(f"Preston curve, {year} (n={n} countries/economies; R²={fit.rsquared:.2f})")111ax.legend(fontsize=8, loc="lower right")112fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150)113114res = {115 "task_id": "F6-dev-01", "analysis_year": year, "n_excluded_missing": n_excl,116 "baseline_metric": baseline, "model_metric": model, "pearson_r": r,117 "r_squared": float(fit.rsquared), "effect_estimate": float(beta), "intercept": float(alpha),118 "effect_se": float(se), "effect_ci_low": float(ci[0]), "effect_ci_high": float(ci[1]),119 "conclusion": concl, "feasibility": "complete", "feasibility_reason": None,120 "data_issues": issues, "n_input_rows": n_rows, "n_analysed": int(n), "random_seed": SEED,121}122with open(f"{OUT}/results.json", "w") as fh:123 json.dump(res, fh, indent=2)124print(json.dumps(res, indent=2)); print(json.dumps(sens, indent=2)); print(t)125