1"""Replication entry point for F6-dev-01 (Preston curve). Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np6import pandas as pd7import statsmodels.api as sm8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112SEED = 0 # no randomness used; fixed for completeness13np.random.seed(SEED)14DATA, OUT = "data", "out"15os.makedirs(OUT, exist_ok=True)1617# keep_default_na=False so that ISO code "NAM" (Namibia) etc. is never parsed as NaN18rd = dict(keep_default_na=False, na_values=[""])19ind = pd.read_csv(os.path.join(DATA, "wb_indicators.csv"), **rd)20ent = pd.read_csv(os.path.join(DATA, "wb_entities.csv"), **rd)21n_input_rows = {"wb_indicators.csv": int(len(ind)), "wb_entities.csv": int(len(ent))}22issues = []2324# ---- audit ----25assert not ind.duplicated(["iso3", "year"]).any()26assert not ent.duplicated("iso3").any()27unmatched = sorted(set(ind.iso3) - set(ent.iso3))28assert not unmatched, unmatched29ent["region"] = ent["region"].str.strip()30ent["income_level"] = ent["income_level"].str.strip()31ent["name"] = ent["name"].str.strip()3233m = ind.merge(ent[["iso3", "region", "income_level"]], on="iso3", how="left")34n_agg_ent = int(m.loc[m.region == "Aggregates", "iso3"].nunique())35issues.append({"issue": f"wb_indicators.csv mixes countries/economies with {n_agg_ent} regional/income/other aggregates (region == 'Aggregates').",36 "handling": "All aggregate entities excluded before any computation, per TASK rules."})37issues.append({"issue": "Some entity names carry trailing whitespace in wb_entities.csv (e.g. 'Latin America & Caribbean ', 'Sub-Saharan Africa ').",38 "handling": "Strings stripped; join is on iso3 so no effect on the analysis."})39issues.append({"issue": "ISO3 code 'NAM' (Namibia) would be read as missing by pandas' default NA parsing.",40 "handling": "CSV read with keep_default_na=False and only empty strings treated as missing."})41c = m[m.region != "Aggregates"].copy()4243# ---- analysis year ----44c["complete"] = c.gdp_pc_ppp.notna() & c.life_expectancy.notna()45per_year = c.groupby("year")["complete"].sum()46eligible = per_year[per_year >= 150]47analysis_year = int(eligible.index.max())4849y = c[c.year == analysis_year]50missing = y[~y.complete]51n_excluded_missing = int(len(missing))52issues.append({"issue": f"In {analysis_year}, {n_excluded_missing} countries/economies lack an indicator (GDP pc PPP missing for: "53 + ", ".join(missing.iso3) + "; LE missing for "54 + str(int(missing.life_expectancy.isna().sum())) + "). Missingness is non-random (small territories, conflict/sanctioned economies such as SSD, SYR, YEM, VEN, PRK, CUB, ERI).",55 "handling": "Listwise exclusion, no imputation, counted in n_excluded_missing; discussed as a selection threat."})5657# implausible LE swings in the panel (outside analysis year)58swing = c.groupby("iso3").life_expectancy.agg(lambda s: s.max() - s.min())59big = swing[swing > 15].index.tolist()60issues.append({"issue": "Very large year-to-year swings in life expectancy for " + ", ".join(big)61 + " (e.g. CAF 18.8 years in 2022 vs 57.4 in 2023; SSD 35-40 in 2015-2017), reflecting crisis-mortality estimates/revisions.",62 "handling": f"Only {analysis_year} values enter the estimate; CAF {analysis_year} value kept (plausible); a sensitivity excluding CAF is reported in the report."})6364d = y[y.complete].copy()65d["log_gdp"] = np.log(d.gdp_pc_ppp)66assert (d.gdp_pc_ppp > 0).all()67n = len(d)6869# ---- baseline and model ----70le = d.life_expectancy.values71baseline_rmse = float(np.sqrt(np.mean((le - le.mean()) ** 2)))72X = sm.add_constant(d.log_gdp.values)73fit = sm.OLS(le, X).fit(cov_type="HC3")74alpha, beta = map(float, fit.params)75resid = fit.resid76model_rmse = float(np.sqrt(np.mean(resid ** 2)))77se = float(fit.bse[1])78ci = fit.conf_int(alpha=0.05)[1] # normal-based 95% CI with HC3 SE79r = float(np.corrcoef(d.log_gdp, le)[0, 1])80r2 = float(fit.rsquared)8182lo, hi = float(ci[0]), float(ci[1])83conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")8485# ---- sensitivity (printed, saved to out/sensitivity.json) ----86sens = {}87for k in ["nonrobust", "HC0", "HC1", "HC2", "HC3"]:88 f = sm.OLS(le, X).fit(cov_type=k)89 sens[f"se_{k}"] = float(f.bse[1])90noc = d[d.iso3 != "CAF"]91f = sm.OLS(noc.life_expectancy.values, sm.add_constant(noc.log_gdp.values)).fit(cov_type="HC3")92sens["beta_excl_CAF"] = float(f.params[1]); sens["se_excl_CAF_HC3"] = float(f.bse[1])93infl = sm.OLS(le, X).fit().get_influence()94cooks = infl.cooks_distance[0]95top = np.argsort(cooks)[::-1][:5]96sens["top5_cooks"] = {d.iso3.iloc[i]: float(cooks[i]) for i in top}97keep = cooks < 4 / n98f = sm.OLS(le[keep], X[keep]).fit(cov_type="HC3")99sens["beta_excl_cooks_gt_4_over_n"] = float(f.params[1]); sens["n_excl_cooks"] = int((~keep).sum())100Xq = np.column_stack([X, (d.log_gdp.values - d.log_gdp.mean()) ** 2])101fq = sm.OLS(le, Xq).fit(cov_type="HC3")102sens["quadratic_term"] = float(fq.params[2]); sens["quadratic_term_p_HC3"] = float(fq.pvalues[2])103sens["rho_spearman"] = float(d[["log_gdp", "life_expectancy"]].corr(method="spearman").iloc[0, 1])104# slope by year (all eligible years) for stability105sens["beta_by_year"] = {}106for yr in eligible.index:107 dy = c[(c.year == yr) & c.complete]108 fy = sm.OLS(dy.life_expectancy.values, sm.add_constant(np.log(dy.gdp_pc_ppp.values))).fit()109 sens["beta_by_year"][int(yr)] = float(fy.params[1])110sens["complete_per_year"] = {int(k): int(v) for k, v in per_year.items()}111sens["largest_neg_resid"] = {d.iso3.iloc[i]: float(resid[i]) for i in np.argsort(resid)[:5]}112sens["largest_pos_resid"] = {d.iso3.iloc[i]: float(resid[i]) for i in np.argsort(resid)[::-1][:5]}113with open(os.path.join(OUT, "sensitivity.json"), "w") as fh:114 json.dump(sens, fh, indent=2)115116# ---- table ----117order = ["Low income", "Lower middle income", "Upper middle income", "High income"]118tab = (d.groupby("income_level")119 .agg(n_countries=("iso3", "size"), mean_log_gdp_pc=("log_gdp", "mean"),120 mean_life_expectancy=("life_expectancy", "mean"))121 .reindex([o for o in order if o in set(d.income_level)] +122 sorted(set(d.income_level) - set(order)))123 .reset_index())124tab.round(4).to_csv(os.path.join(OUT, "table1.csv"), index=False)125126# ---- figure ----127fig, ax = plt.subplots(figsize=(7, 5))128cols = {"Low income": "tab:red", "Lower middle income": "tab:orange",129 "Upper middle income": "tab:green", "High income": "tab:blue"}130for lvl, g in d.groupby("income_level"):131 ax.scatter(g.log_gdp, g.life_expectancy, s=16, alpha=0.75, label=lvl, color=cols.get(lvl, "grey"))132xs = np.linspace(d.log_gdp.min(), d.log_gdp.max(), 100)133ax.plot(xs, alpha + beta * xs, "k-", lw=1.5,134 label=f"OLS: LE = {alpha:.1f} + {beta:.2f}·ln(GDPpc), R²={r2:.2f}")135ax.set_xlabel("ln(GDP per capita, PPP, constant 2021 int'l $)")136ax.set_ylabel("Life expectancy at birth (years)")137ax.set_title(f"Preston curve, {analysis_year} (n={n} countries/economies)")138ax.legend(fontsize=8, loc="lower right")139fig.tight_layout()140fig.savefig(os.path.join(OUT, "figure1.png"), dpi=150)141142results = {143 "task_id": "F6-dev-01",144 "analysis_year": analysis_year,145 "n_excluded_missing": n_excluded_missing,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_ci_low": lo,154 "effect_ci_high": hi,155 "conclusion": conclusion,156 "feasibility": "complete",157 "feasibility_reason": None,158 "data_issues": issues,159 "n_input_rows": n_input_rows,160 "n_analysed": int(n),161 "random_seed": SEED,162}163with open(os.path.join(OUT, "results.json"), "w") as fh:164 json.dump(results, fh, indent=2)165print(json.dumps(results, indent=2)); print(json.dumps(sens, indent=2)); print(tab)166