1"""Replication entry point for F6-dev-01 (Preston curve, WDI).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"14MIN_COUNTRIES = 15015HC = "HC3"16Z = 1.959963984540054 # 97.5% standard-normal quantile1718ind = pd.read_csv(os.path.join(DATA, "wb_indicators.csv"), keep_default_na=False, na_values=[""])19ent = pd.read_csv(os.path.join(DATA, "wb_entities.csv"), keep_default_na=False, na_values=[""])20n_input_rows = {"wb_indicators.csv": int(len(ind)), "wb_entities.csv": int(len(ent))}21data_issues = []22excl = []2324def log_ex(step, rule, before, after, why):25 excl.append(dict(step=step, rule=rule, n_before=int(before), n_after=int(after),26 n_removed=int(before - after), rationale=why))2728# --- integrity checks ---29assert not ind.duplicated(["iso3", "year"]).any(), "duplicate iso3-year rows"30assert not ent.duplicated("iso3").any(), "duplicate entity rows"31unmatched = sorted(set(ind.iso3) - set(ent.iso3))32if unmatched:33 data_issues.append({"issue": f"{len(unmatched)} indicator iso3 codes absent from entity metadata: {unmatched}",34 "handling": "excluded (country status cannot be verified)"})35assert (ind.gdp_pc_ppp.dropna() > 0).all(), "non-positive GDP per capita"3637m = ind.merge(ent[["iso3", "region", "income_level"]], on="iso3", how="left")38n0 = len(m)39m = m[m.region.notna()]40log_ex(1, "drop indicator rows whose iso3 is not in wb_entities.csv", n0, len(m),41 "country/aggregate status cannot be established")42n1 = len(m)43n_agg = m.loc[m.region == "Aggregates", "iso3"].nunique()44c = m[m.region != "Aggregates"].copy()45log_ex(2, "drop entities with region == 'Aggregates' (all years)", n1, len(c),46 f"information rule: countries/economies only ({n_agg} aggregate entities)")47data_issues.append({"issue": f"wb_indicators.csv contains {n_agg} regional/income/other aggregate entities "48 "mixed with countries/economies",49 "handling": "excluded via region == 'Aggregates' in wb_entities.csv"})5051# --- analysis year selection ---52c["complete"] = c.gdp_pc_ppp.notna() & c.life_expectancy.notna()53counts = c.groupby("year")["complete"].sum().astype(int)54eligible = counts[counts >= MIN_COUNTRIES]55if eligible.empty:56 raise SystemExit("No year with >= 150 complete countries: premise invalid")57year = int(eligible.index.max())5859n2 = len(c)60y = c[c.year == year].copy()61log_ex(3, f"keep analysis year {year} only", n2, len(y),62 f"most recent year with >= {MIN_COUNTRIES} countries having both indicators")63n3 = len(y)64miss = y[~y.complete]65d = y[y.complete].copy()66n_excluded_missing = int(len(miss))67log_ex(4, "drop countries with GDPpc or LE missing in analysis year (no imputation)", n3, len(d),68 "information rule; missing: " + "; ".join(miss.iso3))69miss_g = int(miss.gdp_pc_ppp.isna().sum()); miss_l = int(miss.life_expectancy.isna().sum())70data_issues.append({"issue": f"{n_excluded_missing} countries/economies in {year} lack an indicator "71 f"(GDPpc missing: {miss_g}, LE missing: {miss_l}): {', '.join(miss.iso3)}; "72 "mostly small territories and conflict/sanctioned economies, so missingness is not random",73 "handling": "excluded from the estimate, counted in n_excluded_missing; no imputation"})7475# anomalies outside the analysis year (documented, not used)76le_p = c.pivot(index="iso3", columns="year", values="life_expectancy")77jumps = le_p.diff(axis=1).abs().max(axis=1)78big = jumps[jumps > 10].sort_values(ascending=False)79if len(big):80 data_issues.append({"issue": "implausibly large year-to-year life-expectancy swings (>10 years) in: "81 + ", ".join(f"{k} (max |dLE|={v:.1f})" for k, v in big.items())82 + "; e.g. CAF LE 18.8 in 2022 vs 57.4 in 2023; PSE 2023 LE 65.2 vs 76.7 in 2022 (plausibly war-related; not verified)",83 "handling": f"only {year} values enter the estimate; values kept as published (no editing); "84 "influence checked with leave-one-out slopes"})85aggname = ent[ent.name != ent.name.str.strip()]86if len(aggname):87 data_issues.append({"issue": f"{len(aggname)} entity names in wb_entities.csv carry trailing whitespace",88 "handling": "cosmetic; joins use iso3 only"})8990# --- estimation ---91d["log_gdp_pc"] = np.log(d.gdp_pc_ppp)92x = d.log_gdp_pc.to_numpy(); le = d.life_expectancy.to_numpy()93n = len(d)94baseline_rmse = float(np.sqrt(np.mean((le - le.mean()) ** 2)))95X = sm.add_constant(x)96fit = sm.OLS(le, X).fit(cov_type=HC)97alpha, beta = (float(v) for v in fit.params)98se = float(fit.bse[1])99resid = le - (alpha + beta * x)100model_rmse = float(np.sqrt(np.mean(resid ** 2)))101r = float(np.corrcoef(x, le)[0, 1])102r2 = float(fit.rsquared)103lo, hi = beta - Z * se, beta + Z * se104se_hc = {h: float(sm.OLS(le, X).fit(cov_type=h).bse[1]) for h in ["HC0", "HC1", "HC2", "HC3"]}105se_classical = float(sm.OLS(le, X).fit().bse[1])106107# leave-one-out influence108loo = []109for i in range(n):110 k = np.arange(n) != i111 b = np.polyfit(x[k], le[k], 1)[0]112 loo.append(b)113loo = np.array(loo)114iso = d.iso3.to_numpy()115116if lo > 0:117 concl = "positive_effect"118elif hi < 0:119 concl = "negative_effect"120else:121 concl = "no_detectable_effect"122123results = {124 "task_id": "F6-dev-01",125 "analysis_year": year,126 "n_excluded_missing": n_excluded_missing,127 "baseline_metric": baseline_rmse,128 "model_metric": model_rmse,129 "pearson_r": r,130 "r_squared": r2,131 "effect_estimate": beta,132 "intercept": alpha,133 "effect_se": se,134 "effect_ci_low": float(lo),135 "effect_ci_high": float(hi),136 "conclusion": concl,137 "feasibility": "complete",138 "feasibility_reason": None,139 "data_issues": data_issues,140 "n_input_rows": n_input_rows,141 "n_analysed": int(n),142 "random_seed": None,143 # supplementary diagnostics144 "se_type": HC,145 "ci_method": "beta +/- 1.959964 * HC3 SE (normal approximation)",146 "complete_countries_by_year": {str(k): int(v) for k, v in counts.items()},147 "n_countries_in_year": int(n3),148 "n_aggregate_entities_excluded": int(n_agg),149 "effect_se_hc0": se_hc["HC0"], "effect_se_hc1": se_hc["HC1"],150 "effect_se_hc2": se_hc["HC2"], "effect_se_hc3": se_hc["HC3"],151 "effect_se_classical": se_classical,152 "loo_slope_min": float(loo.min()), "loo_slope_min_dropped": str(iso[loo.argmin()]),153 "loo_slope_max": float(loo.max()), "loo_slope_max_dropped": str(iso[loo.argmax()]),154 "rmse_reduction_share": float(1 - model_rmse / baseline_rmse),155 "le_gain_per_doubling_gdp": float(beta * np.log(2)),156 "mean_life_expectancy": float(le.mean()),157 "mean_log_gdp_pc": float(x.mean()),158}159160with open(os.path.join(OUT, "results.json"), "w") as f:161 json.dump(results, f, indent=2)162163# --- table ---164order = ["Low income", "Lower middle income", "Upper middle income", "High income"]165t = (d.groupby("income_level")166 .agg(n_countries=("iso3", "size"), mean_log_gdp_pc=("log_gdp_pc", "mean"),167 mean_life_expectancy=("life_expectancy", "mean"))168 .reset_index())169t["k"] = t.income_level.map({v: i for i, v in enumerate(order)}).fillna(99)170t = t.sort_values(["k", "income_level"]).drop(columns="k")171assert t.n_countries.sum() == n172t.to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.6f")173174# --- exclusions ledger ---175pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(176 os.path.join(OUT, "exclusions.csv"), index=False)177178# --- figure ---179fig, ax = plt.subplots(figsize=(7.5, 5.2))180cols = dict(zip(order, ["#d7191c", "#fdae61", "#abd9e9", "#2c7bb6"]))181for lvl in order:182 s = d[d.income_level == lvl]183 ax.scatter(s.log_gdp_pc, s.life_expectancy, s=18, alpha=0.8, color=cols[lvl], label=f"{lvl} (n={len(s)})")184xx = np.linspace(x.min(), x.max(), 100)185ax.plot(xx, alpha + beta * xx, color="k", lw=1.5,186 label=f"OLS: LE = {alpha:.1f} + {beta:.2f} ln(GDPpc), R² = {r2:.2f}")187ax.set_xlabel("ln(GDP per capita, PPP, constant 2021 int'l $)")188ax.set_ylabel("Life expectancy at birth (years)")189ax.set_title(f"Preston curve, {year}: {n} countries/economies (WDI)")190ax.legend(fontsize=8, loc="lower right")191fig.tight_layout()192fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)193194print(f"year={year} n={n} excluded_missing={n_excluded_missing} beta={beta:.4f} se({HC})={se:.4f} "195 f"CI=[{lo:.4f},{hi:.4f}] r={r:.4f} R2={r2:.4f} RMSE base={baseline_rmse:.4f} model={model_rmse:.4f} -> {concl}")196