1"""Replication entry point for F6-dev-01 (Preston curve).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/; writes out/results.json, out/figure1.png, out/table1.csv, out/exclusions.csv.4No randomness is used (random_seed recorded as null)."""5import json, os6import numpy as np7import pandas as pd8import statsmodels.api as sm9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213DATA, OUT = "data", "out"14os.makedirs(OUT, exist_ok=True)15read = lambda f: pd.read_csv(os.path.join(DATA, f), keep_default_na=False, na_values=[""])16ind, ent = read("wb_indicators.csv"), read("wb_entities.csv")17n_input_rows = {"wb_indicators.csv": int(len(ind)), "wb_entities.csv": int(len(ent))}18issues, excl = [], []1920def log_excl(step, rule, before, after, why):21 excl.append(dict(step=step, rule=rule, n_before=int(before), n_after=int(after),22 n_removed=int(before - after), rationale=why))2324# integrity checks25assert not ind.duplicated(["iso3", "year"]).any() and not ent.iso3.duplicated().any()26assert set(ind.iso3) <= set(ent.iso3)27assert (ind.gdp_pc_ppp.dropna() > 0).all()2829# Step 1: drop aggregates30m = ind.merge(ent[["iso3", "region", "income_level"]], on="iso3", how="left")31n0 = len(m)32n_agg_ent = int(m.loc[m.region == "Aggregates", "iso3"].nunique())33c = m[m.region != "Aggregates"].copy()34log_excl(1, "drop entity-years whose region == 'Aggregates' in wb_entities.csv", n0, len(c),35 "information rule: countries/economies only")36issues.append({"issue": f"wb_indicators.csv contains {n_agg_ent} regional/income/other aggregate entities "37 f"({n0 - len(c)} entity-year rows)", "handling": "excluded via region == 'Aggregates' (rule)"})3839# Step 2: analysis year selection40cov = c.assign(both=c.gdp_pc_ppp.notna() & c.life_expectancy.notna()).groupby("year").both.sum()41elig = cov[cov >= 150]42analysis_year = int(elig.index.max())43n1 = len(c)44y = c[c.year == analysis_year].copy()45log_excl(2, f"keep analysis year {analysis_year} (latest year with >=150 complete countries)", n1, len(y),46 "information rule: single cross-section in analysis year")4748# Step 3: drop missing49miss = y[y.gdp_pc_ppp.isna() | y.life_expectancy.isna()]50n_excluded_missing = int(len(miss))51n2 = len(y)52d = y.dropna(subset=["gdp_pc_ppp", "life_expectancy"]).copy()53log_excl(3, "drop countries with gdp_pc_ppp or life_expectancy missing in analysis year", n2, len(d),54 "information rule: no imputation")55issues.append({"issue": f"{n_excluded_missing} countries/economies lack an indicator in {analysis_year} "56 f"(gdp missing: {int(miss.gdp_pc_ppp.isna().sum())}, LE missing: {int(miss.life_expectancy.isna().sum())}): "57 + ", ".join(sorted(miss.iso3)),58 "handling": "excluded, counted in n_excluded_missing; possible selection bias discussed in report"})59issues.append({"issue": "Central African Republic life expectancy series is implausibly volatile in 2015-2022 "60 "(e.g. 18.8 years in 2022, 31.5 in 2019); 2023 value 57.4 is plausible",61 "handling": "affects non-analysis years only; 2023 value kept; noted as data-quality caveat"})62issues.append({"issue": "West Bank and Gaza life expectancy drops from 76.7 (2022) to 65.2 (2023), consistent with "63 "conflict mortality but a large outlier relative to its income",64 "handling": "kept (genuine reported value, no rule permits dropping it); a leave-5-most-influential-out (Cook distance) sensitivity is reported"})6566# Estimation67d["ln_gdp"] = np.log(d.gdp_pc_ppp)68X = sm.add_constant(d.ln_gdp.values)69yv = d.life_expectancy.values70fit = sm.OLS(yv, X).fit(cov_type="HC3")71alpha, beta = float(fit.params[0]), float(fit.params[1])72se = float(fit.bse[1])73z = 1.95996398454005474ci_lo, ci_hi = beta - z * se, beta + z * se75resid = yv - (alpha + beta * d.ln_gdp.values)76baseline_rmse = float(np.sqrt(np.mean((yv - yv.mean()) ** 2)))77model_rmse = float(np.sqrt(np.mean(resid ** 2)))78r = float(np.corrcoef(d.ln_gdp, yv)[0, 1])79conclusion = "positive_effect" if ci_lo > 0 else ("negative_effect" if ci_hi < 0 else "no_detectable_effect")8081# Supplementary (not estimand): other robust SEs, sensitivity dropping the 5 most influential obs (Cook's D)82fit_ols = sm.OLS(yv, X).fit()83se_alt = {k: float(sm.OLS(yv, X).fit(cov_type=k).bse[1]) for k in ["HC0", "HC1", "HC2", "HC3"]}84cooks = fit_ols.get_influence().cooks_distance[0]85top = np.argsort(cooks)[::-1][:5]86keep = np.setdiff1d(np.arange(len(d)), top)87fit_s = sm.OLS(yv[keep], X[keep]).fit(cov_type="HC3")88d["resid"] = resid8990results = {91 "task_id": "F6-dev-01",92 "feasibility": "complete",93 "feasibility_reason": None,94 "analysis_year": analysis_year,95 "n_excluded_missing": n_excluded_missing,96 "n_analysed": int(len(d)),97 "baseline_metric": baseline_rmse,98 "model_metric": model_rmse,99 "pearson_r": r,100 "r_squared": float(fit.rsquared),101 "effect_estimate": beta,102 "intercept": alpha,103 "effect_se": se,104 "effect_ci_low": float(ci_lo),105 "effect_ci_high": float(ci_hi),106 "conclusion": conclusion,107 "data_issues": issues,108 "n_input_rows": n_input_rows,109 "random_seed": None,110 "se_type": "HC3",111 "complete_countries_by_year": {str(k): int(v) for k, v in cov.items()},112 "n_countries_in_analysis_year": int(n2),113 "supp_effect_se_by_type": se_alt,114 "supp_effect_se_classical": float(fit_ols.bse[1]),115 "supp_mean_life_expectancy": float(yv.mean()),116 "supp_sens_drop5_cooks_iso3": [str(s) for s in d.iso3.values[top]],117 "supp_sens_drop5_effect_estimate": float(fit_s.params[1]),118 "supp_sens_drop5_effect_se": float(fit_s.bse[1]),119 "supp_sens_drop5_r_squared": float(fit_s.rsquared),120 "supp_largest_neg_residual_iso3": str(d.iso3.values[np.argmin(resid)]),121 "supp_largest_neg_residual": float(resid.min()),122}123with open(os.path.join(OUT, "results.json"), "w") as f:124 json.dump(results, f, indent=2)125126pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(127 os.path.join(OUT, "exclusions.csv"), index=False)128129order = ["Low income", "Lower middle income", "Upper middle income", "High income"]130t = (d.groupby("income_level").agg(n_countries=("iso3", "size"), mean_log_gdp_pc=("ln_gdp", "mean"),131 mean_life_expectancy=("life_expectancy", "mean")).reset_index())132t["o"] = t.income_level.map({k: i for i, k in enumerate(order)}).fillna(99)133t = t.sort_values(["o", "income_level"]).drop(columns="o")134t.to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.4f")135136fig, ax = plt.subplots(figsize=(7, 5))137cols = dict(zip(order, ["#d7191c", "#fdae61", "#abd9e9", "#2c7bb6"]))138for lvl in order:139 s = d[d.income_level == lvl]140 ax.scatter(s.ln_gdp, s.life_expectancy, s=18, c=cols[lvl], label=f"{lvl} (n={len(s)})", edgecolor="k", lw=0.3)141xs = np.linspace(d.ln_gdp.min(), d.ln_gdp.max(), 100)142ax.plot(xs, alpha + beta * xs, "k-", lw=1.5, label=f"OLS: LE = {alpha:.1f} + {beta:.2f} ln(GDPpc)")143ax.set_xlabel("ln(GDP per capita, PPP, constant 2021 intl $)")144ax.set_ylabel("Life expectancy at birth (years)")145ax.set_title(f"Preston curve, {analysis_year} (n={len(d)} countries/economies)")146ax.legend(fontsize=8, loc="lower right")147fig.tight_layout()148fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120, metadata={"Software": None})149plt.close(fig)150151print(f"year={analysis_year} n={len(d)} excluded_missing={n_excluded_missing} beta={beta:.4f} "152 f"HC3 se={se:.4f} CI=[{ci_lo:.4f},{ci_hi:.4f}] R2={fit.rsquared:.4f} conclusion={conclusion}")153