1"""WP21 — Python meta-analytic toolkit (independent cross-check of metafor + Bayesian + ML meta-regression).2Functions: reml_re (DerSimonian-Laird start, REML by Newton), pet_peese, bayes_hierarchical (Gibbs with half-Cauchy tau via3slice sampling; R-hat, ESS over 4 chains), ml_meta_regression (RF / GBM with GroupKFold by study), leave_one_out.4Entry point: python scripts/econometrics/python_meta.py -> results/econometrics/python_*.csv5"""6from __future__ import annotations7import json, math, sys8from pathlib import Path9import numpy as np, pandas as pd10from scipy import optimize, stats1112ROOT = Path(__file__).resolve().parents[2]; OUT = ROOT / "results" / "econometrics"; OUT.mkdir(parents=True, exist_ok=True)131415def reml_re(y, v):16 y = np.asarray(y, float); v = np.asarray(v, float); k = len(y)17 if k < 2:18 return dict(mu=y.mean() if k else np.nan, se=np.sqrt(v.mean()) if k else np.nan, tau2=0.0, I2=np.nan, k=k)19 def nll(lt):20 t2 = math.exp(lt); w = 1 / (v + t2); mu = (w * y).sum() / w.sum()21 return 0.5 * (np.log(v + t2).sum() + math.log(w.sum()) + (w * (y - mu) ** 2).sum())22 r = optimize.minimize_scalar(nll, bounds=(-20, 10), method="bounded"); t2 = math.exp(r.x) if r.x > -19 else 0.023 w = 1 / (v + t2); mu = (w * y).sum() / w.sum(); se = math.sqrt(1 / w.sum())24 wf = 1 / v; muf = (wf * y).sum() / wf.sum(); Q = (wf * (y - muf) ** 2).sum(); I2 = max(0.0, (Q - (k - 1)) / Q) * 100 if Q > 0 else 0.025 pi = stats.t.ppf(0.975, max(k - 2, 1)) * math.sqrt(se ** 2 + t2)26 return dict(mu=mu, se=se, ci_lo=mu - 1.96 * se, ci_hi=mu + 1.96 * se, pi_lo=mu - pi, pi_hi=mu + pi, tau2=t2, I2=I2, Q=Q, k=k)272829def pet_peese(y, se):30 y = np.asarray(y, float); se = np.asarray(se, float); v = se ** 2; w = 1 / v31 def wls(X):32 XtW = X.T * w; beta = np.linalg.solve(XtW @ X, XtW @ y); resid = y - X @ beta33 cov = np.linalg.inv(XtW @ X) * max((w * resid ** 2).sum() / (len(y) - X.shape[1]), 1e-12)34 return beta, np.sqrt(np.diag(cov))35 X1 = np.column_stack([np.ones_like(y), se]); b1, s1 = wls(X1)36 X2 = np.column_stack([np.ones_like(y), v]); b2, s2 = wls(X2)37 pet_t = b1[0] / s1[0]38 est = b2[0] if abs(pet_t) > stats.t.ppf(0.975, len(y) - 2) else b1[0]39 return dict(pet_b0=b1[0], pet_se0=s1[0], pet_slope=b1[1], pet_slope_se=s1[1], peese_b0=b2[0], peese_se0=s2[0], pet_peese=est)404142def bayes_hierarchical(y, v, n_iter=6000, burn=2000, chains=4, prior_sd=10.0, tau_scale=5.0, seed=0):43 """y_i ~ N(theta_i, v_i); theta_i ~ N(mu, tau^2); mu ~ N(0, prior_sd^2); tau ~ half-Cauchy(0, tau_scale). Gibbs for theta, mu;44 slice sampling for tau. Returns posterior summaries and R-hat/ESS."""45 y = np.asarray(y, float); v = np.asarray(v, float); k = len(y); rng = np.random.default_rng(seed)46 draws_mu = np.zeros((chains, n_iter - burn)); draws_tau = np.zeros((chains, n_iter - burn))47 for c in range(chains):48 mu = y.mean() + rng.normal() * y.std(ddof=1) if k > 1 else y.mean(); tau = max(y.std(ddof=1), 1e-3) * rng.uniform(0.5, 1.5) if k > 1 else 1.049 theta = y.copy()50 for it in range(n_iter):51 prec = 1 / v + 1 / tau ** 2; theta = rng.normal((y / v + mu / tau ** 2) / prec, np.sqrt(1 / prec))52 pm = (theta.sum() / tau ** 2) / (k / tau ** 2 + 1 / prior_sd ** 2); pv = 1 / (k / tau ** 2 + 1 / prior_sd ** 2); mu = rng.normal(pm, math.sqrt(pv))53 # slice sampling on log tau54 def logp(lt):55 t = math.exp(lt); return -k * lt - ((theta - mu) ** 2).sum() / (2 * t * t) + math.log(1 / (1 + (t / tau_scale) ** 2)) + lt56 lt = math.log(tau); lz = logp(lt) - rng.exponential(); L, R = lt - 1, lt + 157 while logp(L) > lz: L -= 158 while logp(R) > lz: R += 159 for _ in range(100):60 prop = rng.uniform(L, R)61 if logp(prop) > lz: lt = prop; break62 if prop < lt: L = prop63 else: R = prop64 tau = math.exp(lt)65 if it >= burn:66 draws_mu[c, it - burn] = mu; draws_tau[c, it - burn] = tau67 def rhat_ess(x):68 m, n = x.shape; means = x.mean(1); W = x.var(1, ddof=1).mean(); B = n * means.var(ddof=1); var = (n - 1) / n * W + B / n69 rhat = math.sqrt(var / W) if W > 0 else np.nan70 # ESS via autocorrelation (Geyer initial positive sequence, per chain, summed)71 ess = 0.072 for c in range(m):73 xc = x[c] - x[c].mean(); ac = np.correlate(xc, xc, "full")[n - 1:] / (xc @ xc + 1e-300); s = 1.074 for lag in range(1, n // 2):75 if ac[lag] + ac[lag + 1] < 0: break76 s += 2 * ac[lag]77 ess += n / max(s, 1e-9)78 return rhat, ess79 rm, em = rhat_ess(draws_mu); rt, et = rhat_ess(np.log(draws_tau))80 mu_all = draws_mu.ravel(); tau_all = draws_tau.ravel()81 return dict(mu_mean=mu_all.mean(), mu_sd=mu_all.std(), mu_q025=np.quantile(mu_all, .025), mu_q975=np.quantile(mu_all, .975), tau_median=np.median(tau_all),82 tau_q025=np.quantile(tau_all, .025), tau_q975=np.quantile(tau_all, .975), rhat_mu=rm, ess_mu=em, rhat_tau=rt, ess_tau=et, k=k)838485def leave_one_out(df, group):86 rows = []87 for g in df[group].unique():88 s = df[df[group] != g]89 if len(s) >= 3:90 r = reml_re(s.es, s.es_se ** 2); rows.append({group + "_left_out": g, "mu": r["mu"], "se": r["se"], "k": r["k"]})91 return pd.DataFrame(rows)929394def ml_meta_regression(df, features, target="es", group="study_id", seed=20261020):95 from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor96 from sklearn.model_selection import GroupKFold97 from sklearn.inspection import permutation_importance98 X = pd.get_dummies(df[features], dummy_na=True).fillna(0).astype(float); y = df[target].values; groups = df[group].values99 res = {}100 for name, mdl in (("random_forest", RandomForestRegressor(n_estimators=500, min_samples_leaf=5, random_state=seed)), ("gradient_boosting", GradientBoostingRegressor(n_estimators=300, max_depth=3, learning_rate=0.03, random_state=seed))):101 gkf = GroupKFold(n_splits=min(5, len(np.unique(groups)))); pred = np.zeros_like(y, dtype=float)102 for tr, te in gkf.split(X, y, groups):103 mdl.fit(X.iloc[tr], y[tr]); pred[te] = mdl.predict(X.iloc[te])104 ss_res = ((y - pred) ** 2).sum(); ss_tot = ((y - y.mean()) ** 2).sum()105 mdl.fit(X, y); imp = permutation_importance(mdl, X, y, n_repeats=10, random_state=seed)106 res[name] = {"cv_r2_grouped": 1 - ss_res / ss_tot, "cv_rmse": math.sqrt(ss_res / len(y)), "importance": dict(zip(X.columns, imp.importances_mean))}107 return res108109110def main():111 d = pd.read_csv(ROOT / "data" / "final" / "effects.csv"); d = d[np.isfinite(d.es) & np.isfinite(d.es_se) & (d.es_se > 0)]112 rows = []; bayes = []; loo = []113 for (fam, tx), s in d.groupby(["family", "taxonomy_id"]):114 if s.study_id.nunique() < 5:115 continue116 r = reml_re(s.es, s.es_se ** 2); r.update(family=fam, taxonomy_id=tx, k_studies=s.study_id.nunique(), n_estimates=len(s))117 if len(s) >= 10:118 r.update(pet_peese(s.es, s.es_se))119 rows.append(r)120 if s.study_id.nunique() >= 10:121 # study-level aggregation (mean within study, SE pooled) for the Bayesian model to respect independence122 agg = s.groupby("study_id").apply(lambda g: pd.Series({"es": np.average(g.es, weights=1 / g.es_se ** 2), "v": 1 / (1 / g.es_se ** 2).sum()}), include_groups=False)123 b = bayes_hierarchical(agg.es.values, agg.v.values, seed=20261020); b.update(family=fam, taxonomy_id=tx); bayes.append(b)124 l = leave_one_out(s, "study_id"); l["family"] = fam; l["taxonomy_id"] = tx; loo.append(l)125 if "country" in s:126 l2 = leave_one_out(s, "country"); l2["family"] = fam; l2["taxonomy_id"] = tx; loo.append(l2)127 pd.DataFrame(rows).to_csv(OUT / "python_pooled.csv", index=False)128 if bayes: pd.DataFrame(bayes).to_csv(OUT / "python_bayes.csv", index=False)129 if loo: pd.concat(loo).to_csv(OUT / "python_leave_one_out.csv", index=False)130 feats = [c for c in ["taxonomy_id", "family", "continent", "income_group", "data_mid_year", "property_type", "price_type", "estimator_class", "time_fe", "location_fe", "pub_type", "n_obs"] if c in d]131 big = d[d.groupby(["family", "taxonomy_id"]).study_id.transform("nunique") >= 10]132 if len(big) > 50:133 ml = ml_meta_regression(big, feats); (OUT / "python_ml_meta_regression.json").write_text(json.dumps(ml, indent=2, default=float))134 print("pooled", len(rows), "bayes", len(bayes))135136137if __name__ == "__main__":138 main()139