#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ acidification-calcifiers.py Science Journaling Club -- random-effects meta-analysis of calcification and growth responses of marine calcifiers to experimentally lowered pH / raised pCO2. WHAT THIS DOES 1. Validates the pooling machinery on synthetic data with a known true effect, and checks the mathematical identity that random-effects collapses onto fixed-effect when the between-study variance comes out at zero. 2. Pools 23 effect sizes drawn from 15 published experiments using the DerSimonian-Laird random-effects estimator. 3. Reports Cochran's Q, I-squared, tau-squared and a prediction interval. 4. Splits by taxon and tests whether between-taxon differences exceed sampling variation (Q_between). 5. Assesses small-study bias with a funnel plot (printed as coordinates) and Egger's regression test. 6. Runs leave-one-out sensitivity. EFFECT SIZE The log response ratio, y = ln(mean_treatment / mean_control), with the delta-method variance v = sd_t^2 / (n_t * mean_t^2) + sd_c^2 / (n_c * mean_c^2). y = 0 is no effect. y = -0.223 is a 20 percent reduction. Back-transform with percent = 100 * (exp(y) - 1). PROVENANCE Every row carries the archive or table it came from. Rows marked IMPUTED had no dispersion reported anywhere we could find it; those variances are set to the 75th percentile of the variances we could compute, which deliberately gives them little weight. Rows are single-study x single-taxon units; where a paper reported several strains, sites, colony sizes or species inside one taxon we combined them by inverse-variance weighting FIRST, so that no paper contributes more than one row per taxon. Pure standard library. No numpy, no scipy. Runs on any Python 3.8+. """ import math import random # -------------------------------------------------------------------------- # 0. Small statistical helpers, written out so nothing is hidden in a library # -------------------------------------------------------------------------- def norm_cdf(z): """Standard normal cumulative distribution function.""" return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0))) def _betacf(a, b, x, itmax=200, eps=3.0e-12): """Continued fraction for the incomplete beta function (Lentz's method).""" qab, qap, qam = a + b, a + 1.0, a - 1.0 c = 1.0 d = 1.0 - qab * x / qap if abs(d) < 1.0e-30: d = 1.0e-30 d = 1.0 / d h = d for m in range(1, itmax + 1): m2 = 2 * m aa = m * (b - m) * x / ((qam + m2) * (a + m2)) d = 1.0 + aa * d if abs(d) < 1.0e-30: d = 1.0e-30 c = 1.0 + aa / c if abs(c) < 1.0e-30: c = 1.0e-30 d = 1.0 / d h *= d * c aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)) d = 1.0 + aa * d if abs(d) < 1.0e-30: d = 1.0e-30 c = 1.0 + aa / c if abs(c) < 1.0e-30: c = 1.0e-30 d = 1.0 / d delta = d * c h *= delta if abs(delta - 1.0) < eps: break return h def betai(a, b, x): """Regularised incomplete beta function I_x(a, b).""" if x <= 0.0: return 0.0 if x >= 1.0: return 1.0 lbeta = (math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) + a * math.log(x) + b * math.log(1.0 - x)) if x < (a + 1.0) / (a + b + 2.0): return math.exp(lbeta) * _betacf(a, b, x) / a return 1.0 - math.exp(lbeta) * _betacf(b, a, 1.0 - x) / b def t_two_sided_p(t, df): """Two-sided p-value for Student's t with df degrees of freedom.""" if df <= 0: return float("nan") return betai(0.5 * df, 0.5, df / (df + t * t)) def chi2_upper_p(q, df): """Upper tail probability of chi-squared, via a series/continued fraction.""" if df <= 0: return float("nan") if q <= 0: return 1.0 a, x = 0.5 * df, 0.5 * q if x < a + 1.0: # series for P(a, x) ap, s, delta = a, 1.0 / a, 1.0 / a for _ in range(500): ap += 1.0 delta *= x / ap s += delta if abs(delta) < abs(s) * 1e-14: break return 1.0 - s * math.exp(-x + a * math.log(x) - math.lgamma(a)) b, c, d = x + 1.0 - a, 1.0e30, 1.0 / (x + 1.0 - a) # continued fraction h = d for i in range(1, 500): an = -i * (i - a) b += 2.0 d = an * d + b if abs(d) < 1.0e-30: d = 1.0e-30 c = b + an / c if abs(c) < 1.0e-30: c = 1.0e-30 d = 1.0 / d delta = d * c h *= delta if abs(delta - 1.0) < 1e-14: break return math.exp(-x + a * math.log(x) - math.lgamma(a)) * h def t_crit_95(df): """Two-sided 95 percent critical value of Student's t, found by bisection.""" if df <= 0: return float("nan") lo, hi = 1.0, 200.0 for _ in range(200): mid = 0.5 * (lo + hi) if t_two_sided_p(mid, df) > 0.05: lo = mid else: hi = mid return 0.5 * (lo + hi) def pct(y): """Log response ratio to percent change.""" return 100.0 * (math.exp(y) - 1.0) # -------------------------------------------------------------------------- # 1. The pooling engine # -------------------------------------------------------------------------- def fixed_effect(ys, vs): """Inverse-variance fixed-effect pool. Returns (estimate, variance).""" w = [1.0 / v for v in vs] sw = sum(w) return sum(wi * y for wi, y in zip(w, ys)) / sw, 1.0 / sw def cochran_q(ys, vs): """Cochran's Q about the fixed-effect mean, with its degrees of freedom.""" theta_fe, _ = fixed_effect(ys, vs) q = sum((y - theta_fe) ** 2 / v for y, v in zip(ys, vs)) return q, len(ys) - 1 def dersimonian_laird(ys, vs): """ Random-effects meta-analysis, DerSimonian-Laird moment estimator of tau^2. tau^2 = max(0, (Q - (k-1)) / C), C = sum(w) - sum(w^2)/sum(w) Returns a dict with everything downstream code needs. """ k = len(ys) w = [1.0 / v for v in vs] sw, sw2 = sum(w), sum(wi * wi for wi in w) q, df = cochran_q(ys, vs) c = sw - sw2 / sw tau2 = max(0.0, (q - df) / c) if c > 0 else 0.0 ws = [1.0 / (v + tau2) for v in vs] sws = sum(ws) theta = sum(wi * y for wi, y in zip(ws, ys)) / sws se = math.sqrt(1.0 / sws) i2 = max(0.0, 100.0 * (q - df) / q) if q > 0 else 0.0 # H^2: the ratio of observed to expected variation, 1 means none extra h2 = q / df if df > 0 else float("nan") theta_fe, var_fe = fixed_effect(ys, vs) # 95 percent prediction interval: where the NEXT experiment's true effect # should fall. Needs k >= 3 to be meaningful. if k >= 3: tcrit = t_crit_95(k - 2) spread = tcrit * math.sqrt(tau2 + se * se) pi = (theta - spread, theta + spread) else: pi = (float("nan"), float("nan")) return { "k": k, "theta": theta, "se": se, "lo": theta - 1.96 * se, "hi": theta + 1.96 * se, "z": theta / se, "p": 2.0 * (1.0 - norm_cdf(abs(theta / se))), "Q": q, "df": df, "Qp": chi2_upper_p(q, df), "tau2": tau2, "I2": i2, "H2": h2, "theta_fe": theta_fe, "se_fe": math.sqrt(var_fe), "pi_lo": pi[0], "pi_hi": pi[1], } def eggers_test(ys, vs): """ Egger's regression test for funnel asymmetry. Regress the standard normal deviate y/se on precision 1/se; a non-zero intercept is the asymmetry. Weighted-least-squares form, equivalent to regressing y on se with weights 1/v (Sterne & Egger 2001). """ k = len(ys) if k < 3: return None xs = [math.sqrt(v) for v in vs] # standard errors ws = [1.0 / v for v in vs] sw = sum(ws) mx = sum(w * x for w, x in zip(ws, xs)) / sw my = sum(w * y for w, y in zip(ws, ys)) / sw sxx = sum(w * (x - mx) ** 2 for w, x in zip(ws, xs)) sxy = sum(w * (x - mx) * (y - my) for w, x, y in zip(ws, xs, ys)) slope = sxy / sxx intercept = my - slope * mx resid = [y - (intercept + slope * x) for x, y in zip(xs, ys)] dof = k - 2 s2 = sum(w * r * r for w, r in zip(ws, resid)) / dof se_int = math.sqrt(s2 * (1.0 / sw + mx * mx / sxx)) t = intercept / se_int return {"intercept": intercept, "se": se_int, "t": t, "df": dof, "p": t_two_sided_p(t, dof), "slope": slope} def leave_one_out(rows): """Drop each row in turn, repool, and report the movement.""" out = [] for i in range(len(rows)): keep = rows[:i] + rows[i + 1:] r = dersimonian_laird([x["y"] for x in keep], [x["v"] for x in keep]) out.append((rows[i]["label"], r)) return out def residual_tau2(group_lists): """ DerSimonian-Laird estimate of the residual between-study variance, meaning the variance left over AFTER the group means are allowed to differ. tau2_res = max(0, (Q_within - df_within) / C), C = sum(w) - sum_g( sum_g(w^2) / sum_g(w) ) """ all_w, q_within, df_within, c_sub = [], 0.0, 0, 0.0 for g in group_lists: gw = [1.0 / r["v"] for r in g] all_w.extend(gw) if len(g) > 1: q, df = cochran_q([r["y"] for r in g], [r["v"] for r in g]) q_within += q df_within += df c_sub += sum(w * w for w in gw) / sum(gw) c = sum(all_w) - c_sub tau2 = max(0.0, (q_within - df_within) / c) if c > 0 else 0.0 return tau2, q_within, df_within def qm_test(group_lists): """ Mixed-effects moderator test. Group means are compared using weights 1/(v + tau2_residual), so that within-group heterogeneity is charged to the error term instead of being mistaken for a group difference. Q_M = sum_g W_g (theta_g - theta_grand)^2, df = G - 1 This is the test that answers "does taxon explain anything?". The fixed-effect Q_between below does NOT answer that question whenever within-group heterogeneity is present; see the calibration check. """ tau2, q_within, df_within = residual_tau2(group_lists) means = [] for g in group_lists: w = sum(1.0 / (r["v"] + tau2) for r in g) th = sum((1.0 / (r["v"] + tau2)) * r["y"] for r in g) / w means.append((th, w)) wtot = sum(w for _, w in means) grand = sum(w * th for th, w in means) / wtot qm = sum(w * (th - grand) ** 2 for th, w in means) df = len(group_lists) - 1 return {"QM": qm, "df": df, "p": chi2_upper_p(qm, df), "tau2_res": tau2, "Q_within": q_within, "df_within": df_within} def subgroup_test(rows, key="taxon"): """ Split by taxon, pool each group, and run BOTH moderator tests: the naive fixed-effect decomposition and the mixed-effects Q_M. They disagree violently here, and which of them you believe decides the article. """ groups = {} for r in rows: groups.setdefault(r[key], []).append(r) q_total, _ = cochran_q([r["y"] for r in rows], [r["v"] for r in rows]) q_within, df_within = 0.0, 0 per_group = {} for g, rs in groups.items(): ys = [r["y"] for r in rs] vs = [r["v"] for r in rs] if len(rs) > 1: qg, dfg = cochran_q(ys, vs) q_within += qg df_within += dfg per_group[g] = dersimonian_laird(ys, vs) q_between = q_total - q_within df_between = len(groups) - 1 qm = qm_test(list(groups.values())) return { "groups": per_group, "Q_total": q_total, "Q_within": q_within, "Q_between": q_between, "df_between": df_between, "p_between": chi2_upper_p(q_between, df_between), "qm": qm, } def validate_subgroup_calibration(seed=8675309, trials=1500): """ CHECK 3. Simulate study programmes in which every taxon has exactly the same true mean effect, so any "between-taxon difference" found is a false positive. A well-calibrated test rejects 5 percent of the time. Run both moderator tests at three levels of within-group heterogeneity. """ out = [] for tau_w in (0.00, 0.05, 0.28): rng = random.Random(seed + int(tau_w * 1000)) fe_hits, qm_hits = 0, 0 G, n, v = 8, 3, 0.03 for _ in range(trials): gl = [] for _g in range(G): gl.append([{"y": rng.gauss(0.0, math.sqrt(tau_w)) + rng.gauss(0.0, math.sqrt(v)), "v": v} for _ in range(n)]) allr = [r for g in gl for r in g] qt, _ = cochran_q([r["y"] for r in allr], [r["v"] for r in allr]) qw = sum(cochran_q([r["y"] for r in g], [r["v"] for r in g])[0] for g in gl if len(g) > 1) if chi2_upper_p(qt - qw, G - 1) < 0.05: fe_hits += 1 if qm_test(gl)["p"] < 0.05: qm_hits += 1 out.append((tau_w, fe_hits / trials, qm_hits / trials)) return out, trials # -------------------------------------------------------------------------- # 2. Machinery validation, run BEFORE the real data is touched # -------------------------------------------------------------------------- def validate_recovery(seed=20240917): """ Simulate k studies from a known random-effects model and check that the pooled estimate's 95 percent interval covers the true effect, and that the DL tau^2 lands near the true between-study variance. """ rng = random.Random(seed) true_theta, true_tau2, k = -0.300, 0.040, 40 ys, vs = [], [] for _ in range(k): v = rng.uniform(0.004, 0.090) # within-study variance mu_i = rng.gauss(true_theta, math.sqrt(true_tau2)) ys.append(rng.gauss(mu_i, math.sqrt(v))) vs.append(v) r = dersimonian_laird(ys, vs) covered = r["lo"] <= true_theta <= r["hi"] # repeat the whole simulation 500 times and count coverage hits = 0 trials = 500 for t in range(trials): rg = random.Random(seed + 1 + t) yy, vv = [], [] for _ in range(k): v = rg.uniform(0.004, 0.090) mu_i = rg.gauss(true_theta, math.sqrt(true_tau2)) yy.append(rg.gauss(mu_i, math.sqrt(v))) vv.append(v) rr = dersimonian_laird(yy, vv) if rr["lo"] <= true_theta <= rr["hi"]: hits += 1 return r, true_theta, true_tau2, covered, hits, trials def validate_zero_tau(): """ When Q <= k-1 the DL estimator returns tau^2 = 0 exactly, and then the random-effects weights 1/(v + 0) ARE the fixed-effect weights, so the two estimates must be identical to machine precision. This is an algebraic identity, not an empirical result, so a failure here means a coding bug. """ ys = [0.10, 0.11, 0.09, 0.10, 0.105, 0.095] vs = [0.02, 0.03, 0.025, 0.02, 0.028, 0.022] r = dersimonian_laird(ys, vs) diff_theta = abs(r["theta"] - r["theta_fe"]) diff_se = abs(r["se"] - r["se_fe"]) return r, diff_theta, diff_se # -------------------------------------------------------------------------- # 3. The extracted literature # -------------------------------------------------------------------------- # mean_c, sd_c, n_c, mean_t, sd_t, n_t are shown in the provenance string where # a single contrast produced the row. Where several contrasts inside one paper # were combined, the provenance says so and gives the combining rule. # # Inclusion window: treatment pCO2 between roughly 700 and 1500 uatm, or a pH # drop of 0.2 to 0.5 units, against that paper's own ambient control. That is # the end-of-century range, and holding to it is what makes the rows comparable. ROWS = [ dict(label="Riebesell 2000 - coccolithophore", taxon="Coccolithophore", y=-0.45213, v=0.053663, ref=1, imputed=False, prov="PANGAEA 728092 (raw PIC production, pmol/cell/day); our bins " "250-450 uatm control (n=22, mean 0.5649, sd 0.1732) vs " "650-850 uatm (n=8, mean 0.3594, sd 0.2259). Pools E. huxleyi " "and G. oceanica."), dict(label="Gazeau 2007 - mussel", taxon="Mollusc", y=-0.75083, v=0.034583, ref=2, imputed=False, prov="Gazeau et al. 2007 Table 1, incubation-level G in mmol CaCO3 " "g_FW^-1 h^-1. Control = 10 lowest-pCO2 incubations (<=756 " "uatm), mean 0.3390 sd 0.0946; treatment = 9 incubations at " "981-1461 uatm, mean 0.1600 sd 0.0786."), dict(label="Gazeau 2007 - oyster", taxon="Mollusc", y=+0.03709, v=0.002842, ref=2, imputed=False, prov="Gazeau et al. 2007 Table 2. Control = 7 incubations at 698-861 " "uatm, mean 0.2457 sd 0.0315; treatment = 6 at 1063-1258 uatm, " "mean 0.2550 sd 0.0138."), dict(label="Iglesias-Rodriguez 2008 - coccolithophore", taxon="Coccolithophore", y=+0.67940, v=0.000214, ref=3, imputed=False, prov="PANGAEA 718841, PIC production pmol/cell/day. Six independent " "experiments, each 300 uatm vs 750 uatm (Exp6: 280 vs 600), " "n=3 per cell; inverse-variance combined within the paper."), dict(label="Comeau 2009 - pteropod", taxon="Pteropod", y=-0.32542, v=0.001042, ref=4, imputed=False, prov="Comeau et al. 2009 Results: 45Ca calcification 0.36 +/- 0.027 " "(n=10) at pH_T 8.09 vs 0.26 +/- 0.018 (n=10) at pH_T 7.78; " "the paper states '28% lower'."), dict(label="Maier 2009 - cold-water coral", taxon="Coral", y=-0.78945, v=0.036437, ref=5, imputed=False, prov="Maier et al. 2009 Table 2, Skagerrak. 15-Mar: 0.046 +/- 0.010 " "SE (n=8) ambient vs 0.020 +/- 0.003 (n=8) at -0.30 pH. 17-Mar: " "0.021 +/- 0.004 vs 0.010 +/- 0.002. Two runs combined."), dict(label="Langer 2009 - coccolithophore", taxon="Coccolithophore", y=-0.16665, v=0.000493, ref=6, imputed=False, prov="Langer et al. 2009 Table 3, PIC production pg/cell/day, n=3. " "Four E. huxleyi strains, control ~400 uatm vs the level " "nearest 1000 uatm; strains combined. Strain lnRRs: -0.007, " "-0.329, -0.317, -0.439."), dict(label="Ries 2009 - crustaceans", taxon="Crustacean", y=+0.24508, v=0.009435, ref=7, imputed=False, prov="PANGAEA 733947, calcification %/day, 903 vs 409 uatm. " "Callinectes sapidus +0.326, Homarus americanus +0.062, " "Penaeus plebejus +0.584; combined."), dict(label="Ries 2009 - molluscs", taxon="Mollusc", y=-0.39998, v=0.013347, ref=7, imputed=False, prov="PANGAEA 733947, 903 vs 409 uatm, nine species (Argopecten, " "Crassostrea, Crepidula, Littorina, Mercenaria, Mya, Mytilus, " "Strombus, Urosalpinx); combined. Species range -4.07 to +0.96."), dict(label="Ries 2009 - echinoderms", taxon="Echinoderm", y=+0.87247, v=0.078156, ref=7, imputed=False, prov="PANGAEA 733947, 903 vs 409 uatm. Arbacia punctulata +1.652, " "Eucidaris tribuloides -0.442; combined. See the note on the " "archived control-block light metadata."), dict(label="Ries 2009 - temperate coral", taxon="Coral", y=-0.06230, v=0.002550, ref=7, imputed=False, prov="PANGAEA 733947, Oculina arbuscula, 0.1962 (n=11) at 409 uatm " "vs 0.1843 (n=12) at 903 uatm."), dict(label="Ries 2009 - coralline alga", taxon="Coralline alga", y=+0.61984, v=0.044330, ref=7, imputed=False, prov="PANGAEA 733947, Neogoniolithon sp., 0.0955 (n=10) at 409 uatm " "vs 0.1775 (n=10) at 903 uatm."), dict(label="Findlay 2010 - barnacle", taxon="Crustacean", y=+0.03326, v=None, ref=8, imputed=True, prov="PANGAEA 737438, Semibalanus balanoides calcification " "mg g^-1 d^-1, tank means only: pH 8.1 -> 7.7 at both " "temperatures (4.06->4.35 and 4.22->4.21). No dispersion " "archived; variance IMPUTED."), dict(label="Lischka 2011 - pteropod", taxon="Pteropod", y=-0.14409, v=0.003552, ref=9, imputed=False, prov="PANGAEA 761910, shell increment / shell diameter. 350 uatm " "(n=18, mean 0.7361, sd 0.0941) vs 1100 uatm (n=15, mean " "0.6373, sd 0.1269), the three temperatures pooled because " "temperature is crossed and balanced."), dict(label="Chauvin 2011 - Acropora muricata", taxon="Coral", y=-0.04359, v=0.008438, ref=10, imputed=False, prov="PANGAEA 771294, net calcification umol/nubbin/h, 413 vs 1134 " "uatm, at two sites (reef flat -0.179, back reef +0.051), " "n=5 each; combined."), dict(label="Uthicke 2012 - foraminifer", taxon="Foraminifer", y=-0.22966, v=0.003864, ref=11, imputed=False, prov="PANGAEA 831207, Marginopora vertebralis, %/day. TA1 pH 8.1 " "(n=11) vs 7.7 (n=12); TA2 pH 8.1 (n=12) vs 7.8 (n=12); " "combined."), dict(label="Courtney 2013 - tropical urchin", taxon="Echinoderm", y=-1.40765, v=1.562482, ref=12, imputed=False, prov="PANGAEA 824707, Echinometra viridis, percent buoyant-mass " "gain computed by us per individual: 448 uatm 15.86% (n=7, " "sd 10.88) vs 783 uatm 3.88% (n=8, sd 13.42)."), dict(label="Long 2013 - Tanner crab", taxon="Crustacean", y=-0.09531, v=None, ref=13, imputed=True, prov="Long et al. 2013 Results: 'Percent calcium was higher in " "Control crabs than in pH 7.8 or pH 7.5 crabs by 10% and 11%'. " "We use the pH 7.8 figure. No SE given; variance IMPUTED."), dict(label="Schoepf 2013 - Acropora millepora", taxon="Coral", y=-0.75502, v=None, ref=14, imputed=True, prov="Schoepf et al. 2013: calcification 'decreased by 53%' at 741 " "uatm in the second experimental half. Three congeneric species " "in the same paper showed no decrease. Variance IMPUTED."), dict(label="Comeau 2018 - Acropora yongei", taxon="Coral", y=-0.47765, v=0.040656, ref=15, imputed=False, prov="PANGAEA 892655, mg CaCO3 cm^-2 d^-1. Ambient 1.6955 (n=7, " "sd 0.4326) vs high-DIC low-pH 1.0516 (n=7, sd 0.4927)."), dict(label="Comeau 2018 - Pocillopora damicornis", taxon="Coral", y=+0.22028, v=0.073089, ref=15, imputed=False, prov="PANGAEA 892655. Ambient 0.6996 (n=5, sd 0.3930) vs high-DIC " "low-pH 0.8720 (n=6, sd 0.2136)."), dict(label="Comeau 2018 - coralline algae", taxon="Coralline alga", y=-1.33038, v=0.269418, ref=15, imputed=False, prov="PANGAEA 892655. Sporolithon durum 0.0376 (n=5) -> 0.0088 " "(n=5); Neogoniolithon sp. 0.0128 (n=4) -> 0.0040 (n=4); " "combined."), dict(label="Edmunds 2020 - Acropora hyacinthus", taxon="Coral", y=-0.03024, v=0.027932, ref=16, imputed=False, prov="PANGAEA 926042, whole colonies, mg/day, low-temperature arm. " "Large 1046.5 (n=4) -> 989.0 (n=4); small 188.2 (n=4) -> 196.0 " "(n=4); combined."), ] REFS = { 1: "Riebesell et al. 2000, Nature 407:364", 2: "Gazeau et al. 2007, Geophys. Res. Lett. 34:L07603", 3: "Iglesias-Rodriguez et al. 2008, Science 320:336", 4: "Comeau et al. 2009, Biogeosciences 6:1877", 5: "Maier et al. 2009, Biogeosciences 6:1671", 6: "Langer et al. 2009, Biogeosciences 6:2637", 7: "Ries, Cohen & McCorkle 2009, Geology 37:1131", 8: "Findlay et al. 2010, Estuar. Coast. Shelf Sci. 86:675", 9: "Lischka et al. 2011, Biogeosciences 8:919", 10: "Chauvin, Denis & Cuet 2011, Coral Reefs 30:911", 11: "Uthicke & Fabricius 2012, Glob. Change Biol. 18:2781", 12: "Courtney, Westfield & Ries 2013, J. Exp. Mar. Biol. Ecol. 440:169", 13: "Long et al. 2013, PLoS ONE 8:e60959", 14: "Schoepf et al. 2013, PLoS ONE 8:e75049", 15: "Comeau et al. 2018, Glob. Change Biol. 24:4857", 16: "Edmunds & Burgess 2020, J. Exp. Biol. 223:jeb217000", } def impute_variances(rows): """ Fill missing variances with the 75th percentile of the ones we could compute. Conservative on purpose: a fabricated variance should never buy a row much weight. """ known = sorted(r["v"] for r in rows if r["v"] is not None) idx = int(math.ceil(0.75 * len(known))) - 1 fill = known[max(0, idx)] for r in rows: if r["v"] is None: r["v"] = fill return fill # -------------------------------------------------------------------------- # 4. Report # -------------------------------------------------------------------------- W = 78 def rule(ch="-"): print(ch * W) def head(title): print() rule("=") print(title) rule("=") def fmt_row(r): return ("%-42s %+7.3f %+7.1f%% [%+6.1f, %+6.1f]" % (r["label"][:42], r["y"], pct(r["y"]), pct(r["y"] - 1.96 * math.sqrt(r["v"])), pct(r["y"] + 1.96 * math.sqrt(r["v"])))) def report(res, name): print("%-26s k = %d" % (name, res["k"])) print(" pooled lnRR %+.4f (SE %.4f)" % (res["theta"], res["se"])) print(" 95%% CI %+.4f to %+.4f" % (res["lo"], res["hi"])) print(" as percent %+.1f%% [%+.1f%%, %+.1f%%]" % (pct(res["theta"]), pct(res["lo"]), pct(res["hi"]))) print(" z = %.3f, p = %.4g" % (res["z"], res["p"])) print(" Q = %.2f on %d df, p = %.4g" % (res["Q"], res["df"], res["Qp"])) print(" tau^2 = %.4f, I^2 = %.1f%%, H^2 = %.2f" % (res["tau2"], res["I2"], res["H2"])) if not math.isnan(res["pi_lo"]): print(" 95%% prediction %+.1f%% to %+.1f%%" % (pct(res["pi_lo"]), pct(res["pi_hi"]))) def main(): print("SCIENCE JOURNALING CLUB") print("Ocean acidification and calcifiers: sorting the response literature") print("Random-effects meta-analysis of the log response ratio") rule("=") # ---------------- validation ---------------- head("CHECK 1 -- can the machinery recover a known effect?") r, true_theta, true_tau2, covered, hits, trials = validate_recovery() print("Simulated k = %d studies from a known random-effects model." % r["k"]) print(" true effect %+.4f (%.1f%%)" % (true_theta, pct(true_theta))) print(" true tau^2 %.4f" % true_tau2) print(" recovered %+.4f 95%% CI %+.4f to %+.4f" % (r["theta"], r["lo"], r["hi"])) print(" recovered tau^2 %.4f" % r["tau2"]) print(" true value inside its own interval: %s" % ("YES" if covered else "NO")) print(" coverage over %d repeat simulations: %.1f%% (nominal 95%%)" % (trials, 100.0 * hits / trials)) print(" VERDICT: %s" % ("PASS" if covered and 0.90 <= hits / trials <= 0.99 else "FAIL")) head("CHECK 2 -- does random-effects collapse onto fixed-effect at tau^2 = 0?") r0, d_theta, d_se = validate_zero_tau() print("Homogeneous synthetic set, k = %d." % r0["k"]) print(" Q = %.4f on %d df (Q <= df, so DL must return tau^2 = 0)" % (r0["Q"], r0["df"])) print(" tau^2 %.10f" % r0["tau2"]) print(" random-effects %+.10f (SE %.10f)" % (r0["theta"], r0["se"])) print(" fixed-effect %+.10f (SE %.10f)" % (r0["theta_fe"], r0["se_fe"])) print(" |difference| %.3e (estimate), %.3e (SE)" % (d_theta, d_se)) print(" VERDICT: %s" % ("PASS" if d_theta < 1e-12 and d_se < 1e-12 else "FAIL")) head("CHECK 3 -- is the subgroup test calibrated?") print("Every taxon given the SAME true effect, so every rejection is a false") print("positive. A calibrated test rejects 5 percent of the time. G = 8 groups,") print("n = 3 experiments each, within-study variance 0.03.") print() print("%-16s %22s %22s" % ("within-group tau^2", "fixed-effect Q_between", "mixed-effects Q_M")) rule() cal, trials = validate_subgroup_calibration() for tau_w, fe_rate, qm_rate in cal: print("%-16.2f %21.1f%% %21.1f%%" % (tau_w, 100 * fe_rate, 100 * qm_rate)) rule() print("%d simulated programmes per row." % trials) print() print("The fixed-effect Q_between decomposition is fine when within-group") print("heterogeneity is zero and catastrophically anticonservative when it is not.") print("At the level of within-taxon disagreement this literature actually shows,") print("it calls a difference significant in almost every simulated programme where") print("no difference exists. Q_M stays usable, though it is still somewhat") print("anticonservative at n = 3 because the residual variance is itself estimated") print("from very few studies. VERDICT: report Q_M, and treat Q_between as a warning.") # ---------------- data ---------------- fill = impute_variances(ROWS) n_imp = sum(1 for r in ROWS if r["imputed"]) studies = sorted({r["ref"] for r in ROWS}) head("THE EXTRACTED DATA") print("%d effect sizes from %d published experiments." % (len(ROWS), len(studies))) print("%d of %d rows carry a variance computed from reported dispersion; " "%d were imputed at %.5f." % (len(ROWS) - n_imp, len(ROWS), n_imp, fill)) print() print("%-42s %7s %8s %s" % ("study / unit", "lnRR", "percent", "95% CI, percent")) rule() for r in sorted(ROWS, key=lambda x: x["y"]): print(fmt_row(r)) rule() print() print("Provenance") rule() for r in ROWS: print("%s" % r["label"]) print(" source: %s" % REFS[r["ref"]]) words, line = r["prov"].split(), " " for wd in words: if len(line) + len(wd) + 1 > W - 2: print(line) line = " " line += " " + wd print(line) ys = [r["y"] for r in ROWS] vs = [r["v"] for r in ROWS] # ---------------- overall pool ---------------- head("OVERALL POOL") overall = dersimonian_laird(ys, vs) report(overall, "all taxa") print() print("For contrast, the fixed-effect pool is %+.4f (SE %.4f). The random-" % (overall["theta_fe"], overall["se_fe"])) print("effects estimate is wider because tau^2 is not zero. Nothing here is") print("a substitute for looking at the spread.") # ---------------- taxa ---------------- head("BY TAXON") sg = subgroup_test(ROWS) order = sorted(sg["groups"], key=lambda g: sg["groups"][g]["theta"]) print("%-18s %2s %9s %9s %22s" % ("taxon", "k", "lnRR", "percent", "95% CI (percent)")) rule() for g in order: res = sg["groups"][g] print("%-18s %2d %+9.3f %+8.1f%% [%+7.1f%%, %+7.1f%%]" % (g, res["k"], res["theta"], pct(res["theta"]), pct(res["lo"]), pct(res["hi"]))) rule() print() print("Does taxon explain anything? Two tests, two answers.") print(" Q_total = %.2f" % sg["Q_total"]) print(" Q_within = %.2f" % sg["Q_within"]) print() print(" (a) NAIVE fixed-effect decomposition") print(" Q_between = %.2f on %d df, p = %.4g" % (sg["Q_between"], sg["df_between"], sg["p_between"])) print(" Taken at face value this says taxon matters enormously.") print(" Check 3 says this test rejects a true null %.0f%% of the time at the" % (100 * cal[2][1])) print(" level of within-taxon heterogeneity we actually have. Discard it.") print() qm = sg["qm"] print(" (b) MIXED-EFFECTS moderator test, the one that answers the question") print(" residual tau^2 = %.4f (within-taxon disagreement, charged to error)" % qm["tau2_res"]) print(" Q_M = %.2f on %d df, p = %.4g" % (qm["QM"], qm["df"], qm["p"])) if qm["p"] < 0.05: print(" Taxon is a real moderator once within-taxon spread is accounted for.") else: print(" Once within-taxon disagreement is charged to the error term, the") print(" differences BETWEEN taxa are indistinguishable from the differences") print(" WITHIN them. Taxon explains essentially none of the heterogeneity.") print() print("Within-taxon heterogeneity (is the taxon label enough on its own?)") for g in order: res = sg["groups"][g] if res["k"] > 1: print(" %-18s k=%d Q=%6.2f df=%d p=%.4g I^2=%5.1f%% tau^2=%.4f" % (g, res["k"], res["Q"], res["df"], res["Qp"], res["I2"], res["tau2"])) else: print(" %-18s k=1 (single experiment, no within-taxon test)" % g) # ---------------- small-study bias ---------------- head("SMALL-STUDY BIAS") eg = eggers_test(ys, vs) print("Egger's regression test") print(" intercept %+.4f (SE %.4f), t = %.3f on %d df, p = %.4g" % (eg["intercept"], eg["se"], eg["t"], eg["df"], eg["p"])) if eg["p"] < 0.10: print(" The funnel is asymmetric. Small, imprecise experiments report") print(" systematically different effects from large precise ones.") else: print(" No detectable funnel asymmetry. That is weak evidence at this k;") print(" Egger's test has very little power below about 20 studies.") print() print("Egger's test is fragile here, so drop each row in turn and refit:") eg_int, eg_p = [], [] for i in range(len(ROWS)): keep = ROWS[:i] + ROWS[i + 1:] e = eggers_test([x["y"] for x in keep], [x["v"] for x in keep]) eg_int.append(e["intercept"]) eg_p.append(e["p"]) worst = max(range(len(ROWS)), key=lambda i: eg_p[i]) print(" intercept ranges %+.3f to %+.3f; p ranges %.4g to %.4g" % (min(eg_int), max(eg_int), min(eg_p), max(eg_p))) print(" the single most influential row is %s" % ROWS[worst]["label"]) print(" (dropping it: intercept %+.3f, p = %.4g)" % (eg_int[worst], eg_p[worst])) print() print("Funnel plot coordinates (x = lnRR, y = standard error, small at top)") print("%-42s %8s %8s" % ("study / unit", "lnRR", "SE")) rule() for r in sorted(ROWS, key=lambda x: math.sqrt(x["v"])): print("%-42s %+8.3f %8.3f" % (r["label"][:42], r["y"], math.sqrt(r["v"]))) rule() se_max = max(math.sqrt(r["v"]) for r in ROWS) print("Pseudo 95%% limits at the pooled estimate %+.3f:" % overall["theta"]) for frac in (0.25, 0.50, 0.75, 1.00): s = se_max * frac print(" SE %.3f -> %+.3f to %+.3f" % (s, overall["theta"] - 1.96 * s, overall["theta"] + 1.96 * s)) # ---------------- leave one out ---------------- head("LEAVE-ONE-OUT SENSITIVITY") print("%-42s %9s %9s %8s" % ("dropped", "lnRR", "percent", "I^2")) rule() loo = leave_one_out(ROWS) for label, res in loo: print("%-42s %+9.4f %+8.1f%% %7.1f%%" % (label[:42], res["theta"], pct(res["theta"]), res["I2"])) rule() print("%-42s %+9.4f %+8.1f%% %7.1f%%" % ("(none dropped)", overall["theta"], pct(overall["theta"]), overall["I2"])) swing_lo = min(res["theta"] for _, res in loo) swing_hi = max(res["theta"] for _, res in loo) print() print("Pooled estimate ranges from %+.4f to %+.4f (%.1f%% to %.1f%%) across" % (swing_lo, swing_hi, pct(swing_lo), pct(swing_hi))) print("the leave-one-out set. The sign of the pooled effect is %s." % ("stable" if swing_lo * swing_hi > 0 else "NOT stable")) # ---------------- independence check ---------------- head("NON-INDEPENDENCE CHECK") print("Three papers contribute more than one row. Dropping the multi-row") print("papers down to a single inverse-variance-combined row each:") by_ref = {} for r in ROWS: by_ref.setdefault(r["ref"], []).append(r) collapsed = [] for ref, rs in sorted(by_ref.items()): if len(rs) == 1: collapsed.append((rs[0]["y"], rs[0]["v"], REFS[ref])) else: th, vv = fixed_effect([x["y"] for x in rs], [x["v"] for x in rs]) collapsed.append((th, vv, REFS[ref] + " (%d rows collapsed)" % len(rs))) cres = dersimonian_laird([c[0] for c in collapsed], [c[1] for c in collapsed]) for th, vv, nm in sorted(collapsed): print(" %-52s %+7.3f" % (nm[:52], th)) print() report(cres, "one row per paper") # ---------------- verdict ---------------- head("WHAT THE NUMBERS SAY") print("Pooled effect on calcification and growth across all calcifying taxa:") print(" %+.1f%% (95%% CI %+.1f%% to %+.1f%%), k = %d" % (pct(overall["theta"]), pct(overall["lo"]), pct(overall["hi"]), overall["k"])) print("Heterogeneity: I^2 = %.1f%%, tau^2 = %.4f, Q = %.1f on %d df (p = %.3g)" % (overall["I2"], overall["tau2"], overall["Q"], overall["df"], overall["Qp"])) print("Prediction interval for the next experiment: %+.1f%% to %+.1f%%" % (pct(overall["pi_lo"]), pct(overall["pi_hi"]))) print() print("Does taxonomic group explain that spread?") print(" naive fixed-effect Q_between = %.1f on %d df, p = %.3g (miscalibrated)" % (sg["Q_between"], sg["df_between"], sg["p_between"])) print(" mixed-effects Q_M = %.2f on %d df, p = %.3g (the answer: no)" % (sg["qm"]["QM"], sg["qm"]["df"], sg["qm"]["p"])) print() print("Read the prediction interval before the confidence interval. The") print("confidence interval says where the AVERAGE sits. The prediction") print("interval says what the next experiment is likely to find, and it") print("spans both signs.") print() print("Then read Q_M. The spread between taxa looks large, but it is no larger") print("than the spread inside a single taxon, and the honest summary of this") print("literature is that experiments on the same organisms disagree with each") print("other about as much as experiments on different phyla do. That is a") print("statement about the experiments, and it is the finding.") print() rule("=") if __name__ == "__main__": main()