1"""Replication entry point for F6-dev-01 (Preston curve). Run from workspace root."""2import json, os3import numpy as np, pandas as pd4import statsmodels.api as sm5import matplotlib6matplotlib.use("Agg")7import matplotlib.pyplot as plt89SEED = None # no randomness used10os.makedirs("out", exist_ok=True)11# keep_default_na=False so that e.g. Namibia's ISO2-like strings / "NA" are never read as missing12rd = dict(keep_default_na=False, na_values=[""])13ind = pd.read_csv("data/wb_indicators.csv", **rd)14ent = pd.read_csv("data/wb_entities.csv", **rd)15issues = []1617dup = int(ind.duplicated(["iso3", "year"]).sum())18issues.append({"issue": f"wb_indicators.csv mixes countries with regional/income aggregates ({ent.region.eq('Aggregates').sum()} aggregate codes in metadata)",19 "handling": "Dropped every iso3 whose region is 'Aggregates' in wb_entities.csv"})20unmatched = sorted(set(ind.iso3) - set(ent.iso3))21issues.append({"issue": f"{len(unmatched)} iso3 codes in indicators not found in entity metadata; {dup} duplicate iso3-year rows",22 "handling": "Checked; nothing to handle" if not unmatched and dup == 0 else "Unmatched codes excluded"})23issues.append({"issue": "Literal country codes/names such as 'NA' could be parsed as missing by default CSV readers",24 "handling": "Read CSVs with only empty strings treated as missing"})2526countries = set(ent.loc[ent.region != "Aggregates", "iso3"])27c = ind[ind.iso3.isin(countries)].copy()28both = c.gdp_pc_ppp.notna() & c.life_expectancy.notna()29counts = both.groupby(c.year).sum()30year = int(max(y for y, n in counts.items() if n >= 150))3132caf = c[(c.iso3 == "CAF") & (c.life_expectancy < 35)]33if len(caf):34 issues.append({"issue": "Implausible life-expectancy values for Central African Republic in non-analysis years ("35 + ", ".join(f"{int(r.year)}: {r.life_expectancy:.1f}" for r in caf.itertuples()) + ")",36 "handling": f"Not in analysis year {year}; left unchanged, flagged only"})3738y = c[c.year == year].merge(ent[["iso3", "income_level", "region"]], on="iso3")39miss = y[y.gdp_pc_ppp.isna() | y.life_expectancy.isna()]40issues.append({"issue": f"{len(miss)} countries/economies in {year} lack an indicator (all GDP pc PPP missing: "41 + ", ".join(miss.iso3) + "); missingness concentrated in small territories and conflict/sanctioned economies",42 "handling": "Excluded listwise, no imputation; counted in n_excluded_missing"})43a = y.drop(miss.index).copy()44nonpos = int((a.gdp_pc_ppp <= 0).sum())45issues.append({"issue": f"{nonpos} non-positive GDP values (log undefined)", "handling": "Checked; none found" if nonpos == 0 else "Excluded"})46a = a[a.gdp_pc_ppp > 0]47a["lngdp"] = np.log(a.gdp_pc_ppp)4849X = sm.add_constant(a.lngdp.values)50fit = sm.OLS(a.life_expectancy.values, X).fit(cov_type="HC3")51fit0 = sm.OLS(a.life_expectancy.values, X).fit()52le = a.life_expectancy.values53base = float(np.sqrt(np.mean((le - le.mean()) ** 2)))54rmse = float(np.sqrt(np.mean(fit.resid ** 2)))55r = float(np.corrcoef(a.lngdp, le)[0, 1])56beta, se = float(fit.params[1]), float(fit.bse[1])57lo, hi = [float(v) for v in fit.conf_int(0.05)[1]]58concl = "positive_effect" if lo > 0 else "negative_effect" if hi < 0 else "no_detectable_effect"5960# sensitivity (reported in report, not in required keys)61sens = {}62for hc in ["HC0", "HC1", "HC2", "HC3"]:63 f = sm.OLS(le, X).fit(cov_type=hc); sens[hc] = float(f.bse[1])64sens["classical_se"] = float(fit0.bse[1])65yrs = {}66for yy in sorted(c.year.unique()):67 s = c[(c.year == yy) & c.gdp_pc_ppp.notna() & c.life_expectancy.notna()]68 f = sm.OLS(s.life_expectancy.values, sm.add_constant(np.log(s.gdp_pc_ppp.values))).fit(cov_type="HC3")69 yrs[int(yy)] = {"n": len(s), "beta": float(f.params[1]), "se": float(f.bse[1])}70infl = fit0.get_influence(); cd = infl.cooks_distance[0]71top = a.assign(cook=cd, resid=fit0.resid).sort_values("cook", ascending=False).head(5)72trim = a[cd <= 4 / len(a)]73ft = sm.OLS(trim.life_expectancy.values, sm.add_constant(trim.lngdp.values)).fit(cov_type="HC3")74quad = sm.OLS(le, np.column_stack([np.ones(len(a)), a.lngdp, a.lngdp ** 2])).fit(cov_type="HC3")75sensitivity = {"se_by_estimator": sens, "slope_by_year": yrs,76 "top5_cooks": top[["iso3", "name", "gdp_pc_ppp", "life_expectancy", "resid", "cook"]].to_dict("records"),77 "trimmed_cooks_4_over_n": {"n": len(trim), "beta": float(ft.params[1]), "se": float(ft.bse[1])},78 "quadratic_in_log": {"coef": [float(v) for v in quad.params], "se": [float(v) for v in quad.bse]}}79with open("out/sensitivity.json", "w") as fh:80 json.dump(sensitivity, fh, indent=2)8182tab = (a.groupby("income_level").agg(n_countries=("iso3", "size"), mean_log_gdp_pc=("lngdp", "mean"),83 mean_life_expectancy=("life_expectancy", "mean")).reset_index())84order = ["Low income", "Lower middle income", "Upper middle income", "High income"]85tab["k"] = tab.income_level.map({k: i for i, k in enumerate(order)}).fillna(99)86tab.sort_values("k").drop(columns="k").round(4).to_csv("out/table1.csv", index=False)8788fig, ax = plt.subplots(figsize=(7, 5))89for lvl in order:90 s = a[a.income_level == lvl]91 ax.scatter(s.lngdp, s.life_expectancy, s=14, alpha=.75, label=f"{lvl} (n={len(s)})")92xs = np.linspace(a.lngdp.min(), a.lngdp.max(), 100)93ax.plot(xs, fit.params[0] + beta * xs, "k-", lw=2, label=f"OLS: LE = {fit.params[0]:.1f} + {beta:.2f}·ln(GDPpc)")94ax.set_xlabel("ln(GDP per capita, PPP, constant 2021 int'l $)")95ax.set_ylabel("Life expectancy at birth (years)")96ax.set_title(f"Preston curve, {year} (n={len(a)} countries/economies)")97ax.legend(fontsize=8); fig.tight_layout(); fig.savefig("out/figure1.png", dpi=150)9899res = {"task_id": "F6-dev-01", "analysis_year": year, "n_excluded_missing": int(len(miss)),100 "baseline_metric": base, "model_metric": rmse, "pearson_r": r, "r_squared": float(fit.rsquared),101 "effect_estimate": beta, "intercept": float(fit.params[0]), "effect_se": se,102 "effect_ci_low": lo, "effect_ci_high": hi, "conclusion": concl,103 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,104 "n_input_rows": {"wb_indicators.csv": int(len(ind)), "wb_entities.csv": int(len(ent))},105 "n_analysed": int(len(a)), "random_seed": SEED}106with open("out/results.json", "w") as fh:107 json.dump(res, fh, indent=2)108print(json.dumps(res, indent=2)); print(json.dumps(sensitivity, indent=1)); print(tab)109