1"""F6-dev-01: Preston curve, replication entry point.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 np, pandas as pd7import statsmodels.api as sm8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112OUT = "out"13os.makedirs(OUT, exist_ok=True)14rd = dict(keep_default_na=False, na_values=[""])15ind = pd.read_csv("data/wb_indicators.csv", **rd)16ent = pd.read_csv("data/wb_entities.csv", **rd)17n_input_rows = {"wb_indicators.csv": int(len(ind)), "wb_entities.csv": int(len(ent))}1819excl = []20def log(step, rule, before, after, why):21 excl.append(dict(step=step, rule=rule, n_before=before, n_after=after,22 n_removed=before - after, rationale=why))2324assert not ind.duplicated(["iso3", "year"]).any()25unmatched = sorted(set(ind.iso3) - set(ent.iso3))26assert not unmatched, unmatched27df = ind.merge(ent[["iso3", "region", "income_level"]], on="iso3", how="left")2829n0 = len(df)30df = df[df.region != "Aggregates"].copy()31log(1, "drop entities with region == 'Aggregates'", n0, len(df),32 "analysis population is countries/economies only (TASK rule)")3334df["complete"] = df.gdp_pc_ppp.notna() & df.life_expectancy.notna()35counts = df.groupby("year").complete.sum()36eligible = counts[counts >= 150]37if eligible.empty:38 raise SystemExit("no year with >=150 complete countries")39year = int(eligible.index.max())4041n1 = len(df)42d = df[df.year == year].copy()43log(2, f"keep analysis year {year} (most recent year with >=150 complete countries)", n1, len(d),44 "estimand is a single cross-section in the analysis year")45n2 = len(d)46d = d[d.complete].copy()47n_excl_missing = n2 - len(d)48log(3, "drop countries with missing gdp_pc_ppp or life_expectancy", n2, len(d),49 "no imputation (TASK rule)")50assert (d.gdp_pc_ppp > 0).all()51pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)5253d["log_gdp_pc"] = np.log(d.gdp_pc_ppp)54y, x = d.life_expectancy.to_numpy(), d.log_gdp_pc.to_numpy()55n = len(d)5657# Baseline: mean-only prediction58baseline_rmse = float(np.sqrt(np.mean((y - y.mean()) ** 2)))5960# OLS with HC3 robust SE61X = sm.add_constant(x)62fit = sm.OLS(y, X).fit(cov_type="HC3")63alpha, beta = map(float, fit.params)64se = float(fit.bse[1])65ci_lo, ci_hi = map(float, fit.conf_int(alpha=0.05)[1]) # normal quantile with HC366model_rmse = float(np.sqrt(np.mean(fit.resid ** 2)))67r = float(np.corrcoef(x, y)[0, 1])68r2 = float(fit.rsquared)6970if ci_lo > 0: concl = "positive_effect"71elif ci_hi < 0: concl = "negative_effect"72else: concl = "no_detectable_effect"7374missing_list = sorted(df[(df.year == year) & ~df.complete].iso3)75caf = df[(df.iso3 == "CAF")].set_index("year").life_expectancy76data_issues = [77 {"issue": f"wb_indicators.csv mixes countries with {int((ind.merge(ent, on='iso3').region == 'Aggregates').sum() // 9)} aggregate entities (regions/income groups)",78 "handling": "Excluded every entity whose region is 'Aggregates' in wb_entities.csv"},79 {"issue": f"{n_excl_missing} countries/economies lack gdp_pc_ppp in {year}: {', '.join(missing_list)} (small territories plus conflict/sanctioned economies; missingness likely not random)",80 "handling": "Excluded and counted in n_excluded_missing; no imputation"},81 {"issue": "Central African Republic life-expectancy series is implausibly volatile in earlier years (e.g. "82 f"{caf[2019]:.1f} in 2019, {caf[2022]:.1f} in 2022, {caf[2023]:.1f} in 2023)",83 "handling": f"Its {year} value is plausible and kept (no task rule allows dropping it); flagged as a data-quality concern"},84 {"issue": "Life expectancy is reported for every entity-year including tiny territories, indicating modelled (UN WPP) estimates; latest-year values may be provisional",85 "handling": "Used as published; noted as a limitation"},86]8788res = {89 "task_id": "F6-dev-01",90 "feasibility": "complete",91 "feasibility_reason": None,92 "data_issues": data_issues,93 "n_input_rows": n_input_rows,94 "n_analysed": int(n),95 "random_seed": None,96 "analysis_year": year,97 "n_excluded_missing": int(n_excl_missing),98 "baseline_metric": baseline_rmse,99 "model_metric": model_rmse,100 "pearson_r": r,101 "r_squared": r2,102 "effect_estimate": beta,103 "intercept": alpha,104 "effect_se": se,105 "effect_se_type": "HC3",106 "effect_ci_low": ci_lo,107 "effect_ci_high": ci_hi,108 "conclusion": concl,109 "n_countries_total": int(df.iso3.nunique()),110 "complete_counts_by_year": {str(k): int(v) for k, v in counts.items()},111 "effect_se_classical": float(sm.OLS(y, X).fit().bse[1]),112}113with open(f"{OUT}/results.json", "w") as f:114 json.dump(res, f, indent=2)115116# Table 1117order = ["Low income", "Lower middle income", "Upper middle income", "High income"]118t = (d.groupby("income_level")119 .agg(n_countries=("iso3", "size"), mean_log_gdp_pc=("log_gdp_pc", "mean"),120 mean_life_expectancy=("life_expectancy", "mean"))121 .reset_index())122t["_o"] = t.income_level.map({k: i for i, k in enumerate(order)}).fillna(99)123t = t.sort_values(["_o", "income_level"]).drop(columns="_o")124t.round(4).to_csv(f"{OUT}/table1.csv", index=False)125126# Figure 1127fig, ax = plt.subplots(figsize=(7, 5))128colors = dict(zip(order, ["tab:red", "tab:orange", "tab:green", "tab:blue"]))129for lvl, g in d.groupby("income_level"):130 ax.scatter(g.log_gdp_pc, g.life_expectancy, s=14, alpha=0.75, label=lvl, color=colors.get(lvl, "grey"))131xs = np.linspace(x.min(), x.max(), 100)132ax.plot(xs, alpha + beta * xs, "k-", lw=1.5, label=f"OLS: LE = {alpha:.1f} + {beta:.2f} ln(GDPpc)")133ax.set_xlabel("ln(GDP per capita, PPP, constant 2021 int. $)")134ax.set_ylabel("Life expectancy at birth (years)")135ax.set_title(f"Preston curve, {year} (n = {n} countries/economies)")136ax.legend(fontsize=8, loc="lower right")137fig.tight_layout()138fig.savefig(f"{OUT}/figure1.png", dpi=150)139140print(f"year={year} n={n} excluded_missing={n_excl_missing} beta={beta:.3f} se(HC3)={se:.3f} "141 f"CI=[{ci_lo:.3f},{ci_hi:.3f}] R2={r2:.3f} RMSE base={baseline_rmse:.3f} model={model_rmse:.3f} -> {concl}")142