1"""Replication entry point for F6-dev-01 (Preston curve). Run from workspace root:2 ~/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 0 for completeness; nothing depends on it)."""5import json, os6import numpy as np7import pandas as pd8import statsmodels.api as sm9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213SEED = 014np.random.seed(SEED)15TASK_ID = "F6-dev-01"16OUT = "out"17os.makedirs(OUT, exist_ok=True)1819# keep_default_na=False so codes such as "NA" are never parsed as missing; only empty cells are NaN20ind = pd.read_csv("data/wb_indicators.csv", keep_default_na=False, na_values=[""],21 dtype={"iso3": str, "name": str})22ent = pd.read_csv("data/wb_entities.csv", keep_default_na=False, na_values=[""], dtype=str)23n_input_rows = {"wb_indicators.csv": int(len(ind)), "wb_entities.csv": int(len(ent))}24data_issues = []25excl = []2627assert not ind.duplicated(["iso3", "year"]).any()28assert set(ind.iso3) <= set(ent.iso3), "indicator codes missing from metadata"2930# Step 1: drop aggregates31agg = set(ent.loc[ent.region == "Aggregates", "iso3"])32n0 = len(ind)33c = ind[~ind.iso3.isin(agg)].copy()34excl.append(("1", "drop entities with region == 'Aggregates'", n0, len(c),35 "Analysis population is countries/economies only (TASK rule)"))36data_issues.append({"issue": f"wb_indicators.csv mixes {ind.iso3.isin(agg).groupby(ind.iso3).first().sum()} aggregate entities ({n0-len(c)} rows) with countries/economies",37 "handling": "Excluded all entities with region == 'Aggregates' in wb_entities.csv"})3839# Step 2: analysis year40both = c.gdp_pc_ppp.notna() & c.life_expectancy.notna()41cover = both.groupby(c.year).sum()42eligible = cover[(cover >= 150) & (cover.index >= 2015) & (cover.index <= 2023)]43if eligible.empty:44 raise SystemExit("premise_invalid: no year with >=150 complete countries")45year = int(eligible.index.max())46n1 = len(c)47y = c[c.year == year].copy()48excl.append(("2", f"keep analysis year {year} only", n1, len(y),49 "Latest 2015-2023 year with both indicators for >=150 countries/economies"))5051# Step 3: missing indicators52miss = y.gdp_pc_ppp.isna() | y.life_expectancy.isna()53n_excluded_missing = int(miss.sum())54missing_names = sorted(y.loc[miss, "name"])55n2 = len(y)56d = y[~miss].copy()57excl.append(("3", "drop countries with missing GDPpc PPP or life expectancy in analysis year", n2, len(d),58 "No imputation (TASK rule); dropped: " + "; ".join(missing_names)))59data_issues.append({"issue": f"{n_excluded_missing} countries/economies lack GDP per capita PPP (or LE) in {year}: " + ", ".join(missing_names),60 "handling": "Excluded from the estimate, counted in n_excluded_missing; no imputation (possible selection: many are small territories or conflict/sanctioned economies)"})6162# anomaly scan (outside analysis year: documented only)63cs = c.sort_values(["iso3", "year"]).copy()64cs["dle"] = cs.groupby("iso3").life_expectancy.diff()65jumps = cs.loc[cs.dle.abs() > 10, ["iso3", "year", "life_expectancy"]]66if len(jumps):67 data_issues.append({"issue": "Implausible year-to-year life-expectancy swings (>10 years) in: " +68 ", ".join(f"{r.iso3} {r.year} ({r.life_expectancy:.1f})" for r in jumps.itertuples()),69 "handling": f"Only the {year} cross-section is used; values in {year} for these entities were inspected and kept as reported (no editing of source data)"})70assert (d.gdp_pc_ppp > 0).all()7172pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "rationale"]).assign(73 n_removed=lambda t: t.n_before - t.n_after)[74 ["step", "rule", "n_before", "n_after", "n_removed", "rationale"]].to_csv(f"{OUT}/exclusions.csv", index=False)7576# Baseline and model77d["lngdp"] = np.log(d.gdp_pc_ppp)78le = d.life_expectancy.to_numpy(); x = d.lngdp.to_numpy(); n = len(d)79baseline_rmse = float(np.sqrt(np.mean((le - le.mean()) ** 2)))80X = sm.add_constant(x)81fit = sm.OLS(le, X).fit(cov_type="HC3")82alpha, beta = map(float, fit.params)83resid = le - fit.fittedvalues84model_rmse = float(np.sqrt(np.mean(resid ** 2)))85se = float(fit.bse[1])86ci_low, ci_high = map(float, fit.conf_int(alpha=0.05)[1]) # t(n-2) critical value87r = float(np.corrcoef(x, le)[0, 1])88r2 = float(fit.rsquared)89# sensitivity: other HC variants and classical SE (reported, not used for conclusion)90sens = {f"sensitivity_se_{k}": float(sm.OLS(le, X).fit(cov_type=k).bse[1]) for k in ["HC0", "HC1", "HC2"]}91sens["sensitivity_se_classical"] = float(sm.OLS(le, X).fit().bse[1])92# sensitivity: drop entities with erratic LE series (flagged above) present in the analysis sample93flag = set(jumps.iso3) & set(d.iso3)94keep = ~d.iso3.isin(flag).to_numpy()95fs = sm.OLS(le[keep], X[keep]).fit(cov_type="HC3")96sens["sensitivity_beta_drop_erratic_LE"] = float(fs.params[1])97sens["sensitivity_n_drop_erratic_LE"] = int(keep.sum())98# sensitivity: max absolute change in beta from leaving out any single country99loo = [float(sm.OLS(np.delete(le, i), np.delete(X, i, axis=0)).fit().params[1]) for i in range(len(le))]100sens["sensitivity_beta_loo_min"] = float(min(loo)); sens["sensitivity_beta_loo_max"] = float(max(loo))101102if ci_low > 0: concl = "positive_effect"103elif ci_high < 0: concl = "negative_effect"104else: concl = "no_detectable_effect"105106# Table 1107order = ["Low income", "Lower middle income", "Upper middle income", "High income"]108d = d.merge(ent[["iso3", "income_level"]], on="iso3", how="left")109t1 = d.groupby("income_level").agg(n_countries=("iso3", "size"), mean_log_gdp_pc=("lngdp", "mean"),110 mean_life_expectancy=("life_expectancy", "mean")).reset_index()111t1["_o"] = t1.income_level.map({k: i for i, k in enumerate(order)}).fillna(99)112t1 = t1.sort_values(["_o", "income_level"]).drop(columns="_o")113t1.round(4).to_csv(f"{OUT}/table1.csv", index=False)114115# Figure 1116fig, ax = plt.subplots(figsize=(7, 5))117cols = {"Low income": "tab:red", "Lower middle income": "tab:orange", "Upper middle income": "tab:green", "High income": "tab:blue"}118for lvl in order:119 s = d[d.income_level == lvl]120 ax.scatter(s.lngdp, s.life_expectancy, s=14, alpha=0.75, color=cols[lvl], label=f"{lvl} (n={len(s)})")121xx = np.linspace(x.min(), x.max(), 100)122ax.plot(xx, alpha + beta * xx, "k-", lw=1.5, label=f"OLS: LE = {alpha:.1f} + {beta:.2f}·ln(GDPpc), R²={r2:.2f}")123ax.set_xlabel("ln(GDP per capita, PPP, constant 2021 intl $)")124ax.set_ylabel("Life expectancy at birth (years)")125ax.set_title(f"Preston curve, {year} (n={n} countries/economies)")126ax.legend(fontsize=7, loc="lower right")127fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150); plt.close(fig)128129res = {130 "task_id": TASK_ID, "feasibility": "complete", "feasibility_reason": None,131 "analysis_year": year, "n_excluded_missing": n_excluded_missing,132 "baseline_metric": baseline_rmse, "model_metric": model_rmse,133 "pearson_r": r, "r_squared": r2, "effect_estimate": beta, "intercept": alpha,134 "effect_se": se, "effect_se_type": "HC3", "effect_ci_low": ci_low, "effect_ci_high": ci_high,135 "conclusion": concl, "data_issues": data_issues, "n_input_rows": n_input_rows,136 "n_analysed": int(n), "random_seed": SEED,137 "n_countries_in_population": int(c.iso3.nunique()),138 "complete_countries_by_year": {str(k): int(v) for k, v in cover.items()},139 **sens,140}141with open(f"{OUT}/results.json", "w") as f:142 json.dump(res, f, indent=2)143print(f"year={year} n={n} excluded={n_excluded_missing} beta={beta:.4f} se(HC3)={se:.4f} "144 f"CI=[{ci_low:.4f},{ci_high:.4f}] r={r:.4f} R2={r2:.4f} rmse0={baseline_rmse:.3f} rmse={model_rmse:.3f} -> {concl}")145