"""WP21 — Python meta-analytic toolkit (independent cross-check of metafor + Bayesian + ML meta-regression). Functions: reml_re (DerSimonian-Laird start, REML by Newton), pet_peese, bayes_hierarchical (Gibbs with half-Cauchy tau via slice sampling; R-hat, ESS over 4 chains), ml_meta_regression (RF / GBM with GroupKFold by study), leave_one_out. Entry point: python scripts/econometrics/python_meta.py -> results/econometrics/python_*.csv """ from __future__ import annotations import json, math, sys from pathlib import Path import numpy as np, pandas as pd from scipy import optimize, stats ROOT = Path(__file__).resolve().parents[2]; OUT = ROOT / "results" / "econometrics"; OUT.mkdir(parents=True, exist_ok=True) def reml_re(y, v): y = np.asarray(y, float); v = np.asarray(v, float); k = len(y) if k < 2: 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) def nll(lt): t2 = math.exp(lt); w = 1 / (v + t2); mu = (w * y).sum() / w.sum() return 0.5 * (np.log(v + t2).sum() + math.log(w.sum()) + (w * (y - mu) ** 2).sum()) r = optimize.minimize_scalar(nll, bounds=(-20, 10), method="bounded"); t2 = math.exp(r.x) if r.x > -19 else 0.0 w = 1 / (v + t2); mu = (w * y).sum() / w.sum(); se = math.sqrt(1 / w.sum()) 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.0 pi = stats.t.ppf(0.975, max(k - 2, 1)) * math.sqrt(se ** 2 + t2) 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) def pet_peese(y, se): y = np.asarray(y, float); se = np.asarray(se, float); v = se ** 2; w = 1 / v def wls(X): XtW = X.T * w; beta = np.linalg.solve(XtW @ X, XtW @ y); resid = y - X @ beta cov = np.linalg.inv(XtW @ X) * max((w * resid ** 2).sum() / (len(y) - X.shape[1]), 1e-12) return beta, np.sqrt(np.diag(cov)) X1 = np.column_stack([np.ones_like(y), se]); b1, s1 = wls(X1) X2 = np.column_stack([np.ones_like(y), v]); b2, s2 = wls(X2) pet_t = b1[0] / s1[0] est = b2[0] if abs(pet_t) > stats.t.ppf(0.975, len(y) - 2) else b1[0] 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) def bayes_hierarchical(y, v, n_iter=6000, burn=2000, chains=4, prior_sd=10.0, tau_scale=5.0, seed=0): """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; slice sampling for tau. Returns posterior summaries and R-hat/ESS.""" y = np.asarray(y, float); v = np.asarray(v, float); k = len(y); rng = np.random.default_rng(seed) draws_mu = np.zeros((chains, n_iter - burn)); draws_tau = np.zeros((chains, n_iter - burn)) for c in range(chains): 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.0 theta = y.copy() for it in range(n_iter): prec = 1 / v + 1 / tau ** 2; theta = rng.normal((y / v + mu / tau ** 2) / prec, np.sqrt(1 / prec)) 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)) # slice sampling on log tau def logp(lt): t = math.exp(lt); return -k * lt - ((theta - mu) ** 2).sum() / (2 * t * t) + math.log(1 / (1 + (t / tau_scale) ** 2)) + lt lt = math.log(tau); lz = logp(lt) - rng.exponential(); L, R = lt - 1, lt + 1 while logp(L) > lz: L -= 1 while logp(R) > lz: R += 1 for _ in range(100): prop = rng.uniform(L, R) if logp(prop) > lz: lt = prop; break if prop < lt: L = prop else: R = prop tau = math.exp(lt) if it >= burn: draws_mu[c, it - burn] = mu; draws_tau[c, it - burn] = tau def rhat_ess(x): 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 / n rhat = math.sqrt(var / W) if W > 0 else np.nan # ESS via autocorrelation (Geyer initial positive sequence, per chain, summed) ess = 0.0 for c in range(m): xc = x[c] - x[c].mean(); ac = np.correlate(xc, xc, "full")[n - 1:] / (xc @ xc + 1e-300); s = 1.0 for lag in range(1, n // 2): if ac[lag] + ac[lag + 1] < 0: break s += 2 * ac[lag] ess += n / max(s, 1e-9) return rhat, ess rm, em = rhat_ess(draws_mu); rt, et = rhat_ess(np.log(draws_tau)) mu_all = draws_mu.ravel(); tau_all = draws_tau.ravel() 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), 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) def leave_one_out(df, group): rows = [] for g in df[group].unique(): s = df[df[group] != g] if len(s) >= 3: r = reml_re(s.es, s.es_se ** 2); rows.append({group + "_left_out": g, "mu": r["mu"], "se": r["se"], "k": r["k"]}) return pd.DataFrame(rows) def ml_meta_regression(df, features, target="es", group="study_id", seed=20261020): from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor from sklearn.model_selection import GroupKFold from sklearn.inspection import permutation_importance X = pd.get_dummies(df[features], dummy_na=True).fillna(0).astype(float); y = df[target].values; groups = df[group].values res = {} 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))): gkf = GroupKFold(n_splits=min(5, len(np.unique(groups)))); pred = np.zeros_like(y, dtype=float) for tr, te in gkf.split(X, y, groups): mdl.fit(X.iloc[tr], y[tr]); pred[te] = mdl.predict(X.iloc[te]) ss_res = ((y - pred) ** 2).sum(); ss_tot = ((y - y.mean()) ** 2).sum() mdl.fit(X, y); imp = permutation_importance(mdl, X, y, n_repeats=10, random_state=seed) 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))} return res def main(): d = pd.read_csv(ROOT / "data" / "final" / "effects.csv"); d = d[np.isfinite(d.es) & np.isfinite(d.es_se) & (d.es_se > 0)] rows = []; bayes = []; loo = [] for (fam, tx), s in d.groupby(["family", "taxonomy_id"]): if s.study_id.nunique() < 5: continue 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)) if len(s) >= 10: r.update(pet_peese(s.es, s.es_se)) rows.append(r) if s.study_id.nunique() >= 10: # study-level aggregation (mean within study, SE pooled) for the Bayesian model to respect independence 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) b = bayes_hierarchical(agg.es.values, agg.v.values, seed=20261020); b.update(family=fam, taxonomy_id=tx); bayes.append(b) l = leave_one_out(s, "study_id"); l["family"] = fam; l["taxonomy_id"] = tx; loo.append(l) if "country" in s: l2 = leave_one_out(s, "country"); l2["family"] = fam; l2["taxonomy_id"] = tx; loo.append(l2) pd.DataFrame(rows).to_csv(OUT / "python_pooled.csv", index=False) if bayes: pd.DataFrame(bayes).to_csv(OUT / "python_bayes.csv", index=False) if loo: pd.concat(loo).to_csv(OUT / "python_leave_one_out.csv", index=False) 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] big = d[d.groupby(["family", "taxonomy_id"]).study_id.transform("nunique") >= 10] if len(big) > 50: ml = ml_meta_regression(big, feats); (OUT / "python_ml_meta_regression.json").write_text(json.dumps(ml, indent=2, default=float)) print("pooled", len(rows), "bayes", len(bayes)) if __name__ == "__main__": main()