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 = null)."""5import json, os6import numpy as np7import pandas as pd8import statsmodels.api as sm9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213OUT = "out"14os.makedirs(OUT, exist_ok=True)15TASK_ID = "F6-dev-01"16MIN_COUNTRIES = 1501718ind = pd.read_csv("data/wb_indicators.csv", keep_default_na=False, na_values=[""],19 dtype={"iso3": str, "name": str})20ent = pd.read_csv("data/wb_entities.csv", keep_default_na=False, na_values=[""], dtype=str)21n_input_rows = {"wb_indicators.csv": int(len(ind)), "wb_entities.csv": int(len(ent))}2223assert not ind.duplicated(["iso3", "year"]).any()24assert not ent.duplicated("iso3").any()25assert set(ind.iso3) <= set(ent.iso3)2627ex = []28def log_ex(step, rule, nb, na, why):29 ex.append(dict(step=step, rule=rule, n_before=nb, n_after=na, n_removed=nb - na, rationale=why))3031# Step 1: remove aggregates (entity-year rows)32df = ind.merge(ent[["iso3", "region", "income_level"]], on="iso3", how="left")33n0 = len(df)34df = df[df.region != "Aggregates"].copy()35log_ex(1, "drop region == 'Aggregates'", n0, len(df), "analysis population is countries/economies only (TASK info rule)")3637# Analysis year: latest with both indicators for >= 150 countries38complete = df.gdp_pc_ppp.notna() & df.life_expectancy.notna()39counts = df[complete].groupby("year").size()40eligible = counts[counts >= MIN_COUNTRIES]41if eligible.empty:42 raise SystemExit("No year with >=150 complete countries: premise invalid")43year = int(eligible.index.max())4445# Step 2: keep analysis year46n1 = len(df)47d = df[df.year == year].copy()48log_ex(2, f"keep year == {year}", n1, len(d), "latest year 2015-2023 with both indicators for >=150 countries/economies")49n_year = len(d)5051# Step 3: drop missing indicator52d = d[d.gdp_pc_ppp.notna() & d.life_expectancy.notna()].copy()53n_excl_missing = n_year - len(d)54log_ex(3, "drop rows with missing gdp_pc_ppp or life_expectancy", n_year, len(d), "no imputation (TASK info rule)")55missing_iso = sorted(set(df[df.year == year].iso3) - set(d.iso3))56pd.DataFrame(ex).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)5758assert (d.gdp_pc_ppp > 0).all()59d["log_gdp"] = np.log(d.gdp_pc_ppp)60y = d.life_expectancy.to_numpy(); x = d.log_gdp.to_numpy(); n = len(d)6162# Baseline: mean-only prediction63baseline_rmse = float(np.sqrt(np.mean((y - y.mean()) ** 2)))6465# OLS with HC3 robust SE (primary); t(n-2) critical value for 95% CI66X = sm.add_constant(x)67fit = sm.OLS(y, X).fit(cov_type="HC3")68beta, alpha = float(fit.params[1]), float(fit.params[0])69se = float(fit.bse[1])70from scipy import stats71tcrit = float(stats.t.ppf(0.975, n - 2))72ci_low, ci_high = beta - tcrit * se, beta + tcrit * se73resid = fit.resid74model_rmse = float(np.sqrt(np.mean(resid ** 2)))75r = float(np.corrcoef(x, y)[0, 1])76r2 = float(fit.rsquared)7778# Sensitivity: other HC variants and classical SE79se_alt = {f"sensitivity_se_{k.lower()}": float(sm.OLS(y, X).fit(cov_type=k).bse[1]) for k in ("HC0", "HC1", "HC2")}80se_alt["sensitivity_se_classical"] = float(sm.OLS(y, X).fit().bse[1])81# Sensitivity: drop the 5 largest |residual| observations (influence check, descriptive only)82order = np.argsort(-np.abs(resid))83keep = np.ones(n, bool); keep[order[:5]] = False84fit_trim = sm.OLS(y[keep], X[keep]).fit(cov_type="HC3")85top_resid = d.iloc[order[:5]].iso3.tolist()8687if ci_low > 0:88 conclusion = "positive_effect"89elif ci_high < 0:90 conclusion = "negative_effect"91else:92 conclusion = "no_detectable_effect"9394# Table 195lvl_order = ["Low income", "Lower middle income", "Upper middle income", "High income"]96t1 = (d.groupby("income_level")97 .agg(n_countries=("iso3", "size"), mean_log_gdp_pc=("log_gdp", "mean"), mean_life_expectancy=("life_expectancy", "mean"))98 .reset_index())99t1["_o"] = t1.income_level.map({k: i for i, k in enumerate(lvl_order)}).fillna(99)100t1 = t1.sort_values(["_o", "income_level"]).drop(columns="_o")101t1.to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.4f")102103# Figure 1104fig, ax = plt.subplots(figsize=(7, 5))105colors = {"Low income": "tab:red", "Lower middle income": "tab:orange", "Upper middle income": "tab:green", "High income": "tab:blue"}106for lvl in lvl_order + sorted(set(d.income_level) - set(lvl_order)):107 s = d[d.income_level == lvl]108 if len(s):109 ax.scatter(s.log_gdp, s.life_expectancy, s=14, alpha=0.75, color=colors.get(lvl, "gray"), label=f"{lvl} (n={len(s)})")110xx = np.linspace(x.min(), x.max(), 100)111ax.plot(xx, alpha + beta * xx, "k-", lw=1.5, label=f"OLS: LE = {alpha:.1f} + {beta:.2f} ln(GDPpc), R² = {r2:.2f}")112ax.set_xlabel("ln(GDP per capita, PPP, constant 2021 int'l $)")113ax.set_ylabel("Life expectancy at birth (years)")114ax.set_title(f"Preston curve, {year} ({n} countries/economies, WDI)")115ax.legend(fontsize=7, loc="lower right")116fig.tight_layout()117fig.savefig(os.path.join(OUT, "figure1.png"), dpi=150)118plt.close(fig)119120caf = d[d.iso3 == "CAF"]121data_issues = [122 {"issue": "wb_indicators.csv contains regional/income-group aggregates (79 entities with region 'Aggregates')",123 "handling": "Excluded via wb_entities.csv region == 'Aggregates' before choosing the analysis year."},124 {"issue": f"{n_excl_missing} countries/economies lack GDP per capita PPP in {year} (" + ", ".join(missing_iso) + "); missingness is non-random (small territories, conflict/sanctioned states)",125 "handling": "Excluded without imputation as required; counted in n_excluded_missing; possible selection bias discussed in report."},126 {"issue": "Central African Republic life expectancy is erratic across 2015-2022 (e.g. 18.8 in 2022, 31.5 in 2019)",127 "handling": f"Only the analysis-year value is used ({float(caf.life_expectancy.iloc[0]) if len(caf) else 'n/a'} in {year}), which is plausible; kept. Influence assessed by the drop-5-largest-residuals sensitivity."},128 {"issue": "West Bank and Gaza life expectancy falls ~11.5 years in 2023 (conflict)",129 "handling": "Kept as reported; real data feature, covered by the residual-trimming sensitivity if influential."},130 {"issue": "Recent-year WDI values may be modelled estimates/projections rather than measurements",131 "handling": "Not verifiable from the provided files; noted as a limitation."},132]133134res = {135 "task_id": TASK_ID,136 "feasibility": "complete",137 "feasibility_reason": None,138 "data_issues": data_issues,139 "n_input_rows": n_input_rows,140 "n_analysed": int(n),141 "random_seed": None,142 "analysis_year": year,143 "n_countries_in_year": int(n_year),144 "n_excluded_missing": int(n_excl_missing),145 "complete_countries_by_year": {str(k): int(v) for k, v in counts.items()},146 "baseline_metric": baseline_rmse,147 "model_metric": model_rmse,148 "pearson_r": r,149 "r_squared": r2,150 "effect_estimate": beta,151 "intercept": alpha,152 "effect_se": se,153 "effect_se_type": "HC3",154 "effect_ci_low": float(ci_low),155 "effect_ci_high": float(ci_high),156 "ci_critical_value_t": tcrit,157 "effect_t_stat": beta / se,158 "conclusion": conclusion,159 **se_alt,160 "sensitivity_trim5_effect": float(fit_trim.params[1]),161 "sensitivity_trim5_se_hc3": float(fit_trim.bse[1]),162 "sensitivity_trim5_dropped_iso3": top_resid,163 "mean_life_expectancy": float(y.mean()),164}165with open(os.path.join(OUT, "results.json"), "w") as fh:166 json.dump(res, fh, indent=2)167print(f"year={year} n={n} excluded_missing={n_excl_missing} beta={beta:.4f} se_HC3={se:.4f} "168 f"CI=[{ci_low:.4f},{ci_high:.4f}] R2={r2:.4f} rmse_base={baseline_rmse:.4f} rmse_ols={model_rmse:.4f} -> {conclusion}")169