""" insect-decline-meta.py -- the Science Journaling Club's own random-effects meta-analysis for "The Insect Decline Literature Disagrees With Itself, and We Wanted to Know How Badly" (Review, Meta-analysis, Ecology). WHAT THIS IS ------------ A club-built DerSimonian-Laird random-effects pooling of published long-term insect abundance and biomass trends, plus the diagnostics that decide whether pooling was a reasonable thing to do at all: Cochran's Q, I-squared, a funnel plot, Egger's regression test for small-study asymmetry, and leave-one-out. WHAT THIS IS NOT ---------------- A systematic review. We are a school club. We did not run two independent screeners over a pre-registered search string, we did not contact authors for unpublished data, and our search was "read the obvious papers and chase their reference lists." Every number below is traceable to a printed sentence or a printed table in the cited paper; the file records where. But the SET of papers is a convenience sample, and that is the single largest weakness of the pooled estimate. The article says so at length. THE COMMON SCALE ---------------- Every effect is converted to rho, the natural-log annual rate of change of total abundance or biomass: N(t) = N(0) * exp(rho * t) so rho = ln(N(t)/N(0)) / t and the percent-per-year figure people quote is 100 * (exp(rho) - 1). Per decade it is 100 * (exp(10*rho) - 1). A log-linear rate is the only scale on which a 27-year German biomass series and a 21-year Ohio butterfly count are even in principle comparable. Conversions used, with the source of each raw number recorded in STUDIES: log10 slope b (per year), SE s -> rho = b*ln(10), SE = s*ln(10) percent per year p with 95% CI [lo, hi] (also in percent per year): rho = ln(1 + p/100) SE = [ln(1+hi/100) - ln(1+lo/100)] / (2 * 1.959964) RANDOM EFFECTS, DERSIMONIAN-LAIRD --------------------------------- w_i = 1 / v_i (fixed-effect weights, v_i = SE_i^2) Q = sum w_i (y_i - y_FE)^2 C = sum w_i - (sum w_i^2)/(sum w_i) tau2 = max(0, (Q - (k-1)) / C) w*_i = 1 / (v_i + tau2) y_RE = sum w*_i y_i / sum w*_i SE(y_RE)= 1 / sqrt(sum w*_i) I2 = max(0, (Q - (k-1)) / Q) Egger's test is the weighted regression of y_i/SE_i on 1/SE_i; the intercept is the asymmetry statistic. Egger detects asymmetry, not bias. Standard library only. Python 3. No numpy, no scipy, no external packages, so that anyone with a bare Python install can run it and get the same digits. """ import math LN10 = math.log(10.0) Z95 = 1.959963984540054 # --------------------------------------------------------------------------- # small numerical helpers (no scipy in the room) # --------------------------------------------------------------------------- def norm_cdf(x): """Standard normal CDF via the error function.""" return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0))) def norm_sf2(x): """Two-sided normal tail probability.""" return 2.0 * (1.0 - norm_cdf(abs(x))) def gammap_series(a, x, itmax=500, eps=1e-14): ap = a total = 1.0 / a delta = total for _ in range(itmax): ap += 1.0 delta *= x / ap total += delta if abs(delta) < abs(total) * eps: break return total * math.exp(-x + a * math.log(x) - math.lgamma(a)) def gammaq_cf(a, x, itmax=500, eps=1e-14): tiny = 1e-300 b = x + 1.0 - a c = 1.0 / tiny d = 1.0 / b h = d for i in range(1, itmax + 1): an = -i * (i - a) b += 2.0 d = an * d + b if abs(d) < tiny: d = tiny c = b + an / c if abs(c) < tiny: c = tiny d = 1.0 / d delta = d * c h *= delta if abs(delta - 1.0) < eps: break return math.exp(-x + a * math.log(x) - math.lgamma(a)) * h def gammaq(a, x): if x < a + 1.0: return 1.0 - gammap_series(a, x) return gammaq_cf(a, x) def chi2_sf(x, k): """ Upper tail of the chi-square distribution with k degrees of freedom, via the regularised incomplete gamma Q(k/2, x/2). """ if x <= 0: return 1.0 return gammaq(k / 2.0, x / 2.0) def betacf(a, b, x, itmax=300, eps=1e-14): tiny = 1e-300 qab, qap, qam = a + b, a + 1.0, a - 1.0 c = 1.0 d = 1.0 - qab * x / qap if abs(d) < tiny: d = tiny 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) < tiny: d = tiny c = 1.0 + aa / c if abs(c) < tiny: c = tiny 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) < tiny: d = tiny c = 1.0 + aa / c if abs(c) < tiny: c = tiny d = 1.0 / d delta = d * c h *= delta if abs(delta - 1.0) < eps: break return h def betainc(a, b, x): 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) if x < (a + 1.0) / (a + b + 2.0): front = math.exp(lbeta + a * math.log(x) + b * math.log(1.0 - x)) return front * betacf(a, b, x) / a back = math.exp(lbeta + b * math.log(1.0 - x) + a * math.log(x)) return 1.0 - back * betacf(b, a, 1.0 - x) / b def t_sf2(t, df): """ Two-sided p-value for Student's t, as the regularised incomplete beta I_{df/(df+t^2)}(df/2, 1/2). """ if df <= 0: return float('nan') x = df / (df + t * t) return betainc(df / 2.0, 0.5, x) # --------------------------------------------------------------------------- # THE EXTRACTED STUDIES # --------------------------------------------------------------------------- # Each row is: key, short label, locality, span, taxa/metric, sample size, # rho (ln per year), SE(rho), and PROVENANCE -- the sentence or table the # number came out of. If you cannot check the provenance, do not trust the row. # # Rule we set before extracting: one effect per independent monitoring # dataset, chosen a priori as the paper's own headline series for that # dataset, never the largest effect on offer. STUDIES = [ dict( key="hallmann17", label="Hallmann 2017", locality="63 reserves, Germany", span="1989-2016 (27 yr)", taxa="flying insects, Malaise biomass", n="1,503 trap samples / 63 sites", rho=-0.0630, se=0.0020, source=("PLoS ONE 12(10):e0185809, Results: 'annual trend " "coefficient = -0.063, sd = 0.002, i.e. 6.1% annual decline'"), ), dict( key="shortall_here", label="Shortall 2009 Hereford", locality="Hereford, England", span="1973-2002 (30 yr)", taxa="aerial biomass, 12.2 m suction trap", n="30 annual biomass indices", rho=-0.01885 * LN10, se=0.00410 * LN10, source=("Insect Conserv. Divers. 2:251-260, Table 1 'Biomass " "Hereford': slope -0.01885, SE 0.00410 (log10 g per sample)"), ), dict( key="shortall_roth", label="Shortall 2009 Rothamsted", locality="Harpenden, England", span="1973-2002 (30 yr)", taxa="aerial biomass, 12.2 m suction trap", n="30 annual biomass indices", rho=0.00557 * LN10, se=0.00409 * LN10, source=("Insect Conserv. Divers. 2:251-260, Table 1 'Biomass " "Rothamsted': slope +0.00557, SE 0.00409 (log10)"), ), dict( key="shortall_star", label="Shortall 2009 Starcross", locality="Starcross, Devon, England", span="1973-2002 (30 yr)", taxa="aerial biomass, 12.2 m suction trap", n="30 annual biomass indices", rho=-0.00217 * LN10, se=0.00265 * LN10, source=("Insect Conserv. Divers. 2:251-260, Table 1 'Biomass " "Starcross': slope -0.00217, SE 0.00265 (log10)"), ), dict( key="shortall_wye", label="Shortall 2009 Wye", locality="Wye, Kent, England", span="1973-2002 (30 yr)", taxa="aerial biomass, 12.2 m suction trap", n="30 annual biomass indices", rho=-0.00129 * LN10, se=0.00265 * LN10, source=("Insect Conserv. Divers. 2:251-260, Table 1 'Biomass " "Wye': slope -0.00129, SE 0.00265 (log10)"), ), dict( key="hallmann20_lep", label="Hallmann 2020 Kaaistoep moths", locality="De Kaaistoep, Netherlands", span="1997-2017 (21 yr)", taxa="macro-moths at light", n="497 trapping nights / 54,492 individuals", rho=-0.040, se=0.006, source=("Insect Conserv. Divers. 13:127-139, Table 2 Lepidoptera: " "rho -0.040 (0.006)"), ), dict( key="hallmann20_col", label="Hallmann 2020 Kaaistoep beetles", locality="De Kaaistoep, Netherlands", span="1997-2017 (21 yr)", taxa="beetles at light", n="572 trapping nights / 257,793 individuals", rho=-0.048, se=0.010, source=("Insect Conserv. Divers. 13:127-139, Table 2 Coleoptera: " "rho -0.048 (0.010)"), ), dict( key="hallmann20_car", label="Hallmann 2020 Wijster carabids", locality="Wijster, Drenthe, Netherlands", span="1986-2016 (26 sampled yr)", taxa="ground beetles, 48 pitfall traps", n="264,986 individuals / 156 species", rho=-0.044, se=0.006, source=("Insect Conserv. Divers. 13:127-139, Results: 'rho = -0.044, " "se = 0.006, P < 0.001, 4% decline per year'"), ), dict( key="wepprich19", label="Wepprich 2019 Ohio", locality="Ohio, USA (104 sites)", span="1996-2016 (21 yr)", taxa="butterflies, 81 species, transects", n="24,405 surveys", rho=-0.020, se=0.005, source=("PLoS ONE 14(7):e0216270, Results: 'declined at an annual " "rate of 2.0% (b1 = -0.020, std. err. 0.005, p < 0.001)'"), ), dict( key="edwards25", label="Edwards 2025 contiguous US", locality="contiguous USA (2,478 sites)", span="2000-2020 (21 yr)", taxa="butterflies, 554 species, 35 programmes", n="76,957 surveys / 12.6 M individuals", pct=-1.3, pct_lo=-2.3, pct_hi=-0.2, source=("Science 387:1090-1094: 'a rate of 1.3% annually [95% " "confidence interval: -2.3%, -0.2%]'"), ), dict( key="weiss24", label="Weiss 2024 NE Germany", locality="beech forest, NE Germany", span="1999-2022 (24 yr)", taxa="carabid beetles, abundance", n="24-year pitfall series", pct=-3.1, pct_lo=-5.3, pct_hi=-1.0, source=("Ecography 2024:e07020, Abstract: 'significant linear " "declines in abundance and biomass with annual rates of " "-3.1% (95% CI [-5.3, -1]) and -4.9% (95% CI [-9.4, -1.6])'"), ), ] def finish_rows(rows): """Fill in rho/SE for rows given as percent-per-year with a 95% CI.""" out = [] for r in rows: r = dict(r) if "rho" not in r: r["rho"] = math.log(1.0 + r["pct"] / 100.0) hi = math.log(1.0 + r["pct_hi"] / 100.0) lo = math.log(1.0 + r["pct_lo"] / 100.0) r["se"] = (hi - lo) / (2.0 * Z95) out.append(r) return out STUDIES = finish_rows(STUDIES) # --------------------------------------------------------------------------- # the machinery # --------------------------------------------------------------------------- def fixed_effect(ys, ses): ws = [1.0 / (s * s) for s in ses] sw = sum(ws) est = sum(w * y for w, y in zip(ws, ys)) / sw se = math.sqrt(1.0 / sw) return est, se, ws def dersimonian_laird(ys, ses): """ Returns the fixed-effect estimate, Q, tau^2, the random-effects estimate and its interval, I^2 and H^2. """ k = len(ys) fe, fe_se, ws = fixed_effect(ys, ses) Q = sum(w * (y - fe) ** 2 for w, y in zip(ws, ys)) df = k - 1 sw = sum(ws) sw2 = sum(w * w for w in ws) C = sw - sw2 / sw tau2 = max(0.0, (Q - df) / C) if C > 0 else 0.0 wstar = [1.0 / (s * s + tau2) for s in ses] swstar = sum(wstar) re = sum(w * y for w, y in zip(wstar, ys)) / swstar re_se = math.sqrt(1.0 / swstar) I2 = max(0.0, (Q - df) / Q) * 100.0 if Q > 0 else 0.0 H2 = Q / df if df > 0 else float('nan') return dict( k=k, fe=fe, fe_se=fe_se, Q=Q, df=df, Q_p=chi2_sf(Q, df), tau2=tau2, tau=math.sqrt(tau2), C=C, re=re, re_se=re_se, re_lo=re - Z95 * re_se, re_hi=re + Z95 * re_se, I2=I2, H2=H2, pi_lo=re - Z95 * math.sqrt(tau2 + re_se ** 2), pi_hi=re + Z95 * math.sqrt(tau2 + re_se ** 2), z=re / re_se, p=norm_sf2(re / re_se), weights=wstar, swstar=swstar, ) def eggers_test(ys, ses): """ Egger's regression: the standard normal deviate y/se regressed on the precision 1/se. The intercept, and its t-test, are the asymmetry statistic. Ordinary least squares on (x = 1/se, z = y/se). """ k = len(ys) xs = [1.0 / s for s in ses] zs = [y / s for y, s in zip(ys, ses)] mx = sum(xs) / k mz = sum(zs) / k sxx = sum((x - mx) ** 2 for x in xs) sxz = sum((x - mx) * (z - mz) for x, z in zip(xs, zs)) slope = sxz / sxx intercept = mz - slope * mx resid = [z - (intercept + slope * x) for x, z in zip(xs, zs)] df = k - 2 s2 = sum(r * r for r in resid) / df se_int = math.sqrt(s2 * (1.0 / k + mx * mx / sxx)) t = intercept / se_int return dict(intercept=intercept, se=se_int, t=t, df=df, p=t_sf2(t, df), slope=slope) def leave_one_out(rows): out = [] for i in range(len(rows)): sub = rows[:i] + rows[i + 1:] r = dersimonian_laird([s["rho"] for s in sub], [s["se"] for s in sub]) out.append((rows[i]["label"], r)) return out def pct_yr(rho): return 100.0 * (math.exp(rho) - 1.0) def pct_dec(rho): return 100.0 * (math.exp(10.0 * rho) - 1.0) # --------------------------------------------------------------------------- # VALIDATION -- run the machinery on data whose answer we already know # --------------------------------------------------------------------------- class LCG: """ Tiny seeded linear congruential generator plus Box-Muller, so the synthetic check is reproducible on any Python build without depending on the implementation details of the random module. """ def __init__(self, seed): self.s = seed & 0xFFFFFFFFFFFF def u(self): self.s = (25214903917 * self.s + 11) & 0xFFFFFFFFFFFF return (self.s >> 16) / 4294967296.0 def normal(self): u1 = max(self.u(), 1e-12) u2 = self.u() return math.sqrt(-2.0 * math.log(u1)) * math.cos(2.0 * math.pi * u2) def validate_recovery(true_mu=-0.025, true_tau=0.015, k=14, seed=20260913, trials=2000): """ Generate k studies from a known random-effects model, pool them, and check that the stated interval covers the truth. Then repeat many times and report coverage, which should land near 95%, and the mean tau estimate. """ rng = LCG(seed) ses = [0.002 + 0.010 * rng.u() for _ in range(k)] ys = [true_mu + true_tau * rng.normal() + s * rng.normal() for s in ses] res = dersimonian_laird(ys, ses) cover = 0 tau_hats = [] rng2 = LCG(seed + 1) for _ in range(trials): s2 = [0.002 + 0.010 * rng2.u() for _ in range(k)] y2 = [true_mu + true_tau * rng2.normal() + s * rng2.normal() for s in s2] r2 = dersimonian_laird(y2, s2) if r2["re_lo"] <= true_mu <= r2["re_hi"]: cover += 1 tau_hats.append(r2["tau"]) return res, ys, ses, cover / trials, sum(tau_hats) / trials def validate_tau_zero(): """ When observed dispersion is no larger than sampling error, DL sets tau^2 = 0 and the random-effects weights collapse onto the fixed-effect weights. The two estimates and their standard errors must then agree to machine precision. Build such a dataset by construction: put every study exactly on the fixed-effect mean, so Q = 0 exactly. """ ses = [0.004, 0.009, 0.006, 0.012, 0.003, 0.007] ys = [-0.02] * len(ses) r = dersimonian_laird(ys, ses) return r, ys, ses # --------------------------------------------------------------------------- # SITE SELECTION BIAS -- regression to the mean with no real trend # --------------------------------------------------------------------------- def selection_bias_sim(sel_strength, n_sites=400, n_years=27, seed=7, site_sd=1.0, year_sd=0.6, true_rho=0.0): """ Monitoring programmes are not sited at random. They are sited where there is something worth watching. Model site i's log abundance in year t as x[i,t] = mu_i + e[i,t] + true_rho * t with mu_i ~ N(0, site_sd^2) a permanent site quality and e[i,t] ~ N(0, year_sd^2) independent year-to-year noise. Nothing in this model declines unless true_rho < 0. Now select sites for monitoring on their FIRST year's observed value: rank sites by x[i,0] and keep the top fraction (1 - sel_strength). At sel_strength = 0 every site is kept and there is no selection; at 0.9 only the top tenth of first-year observations gets monitored. Then fit an ordinary least squares log-linear trend to each kept site's series and report the mean slope. With true_rho = 0 the honest answer is zero. What comes out is not zero, because e[i,0] was large by construction and has nowhere to go but down. """ rng = LCG(seed) mus = [site_sd * rng.normal() for _ in range(n_sites)] series = [] for i in range(n_sites): row = [mus[i] + year_sd * rng.normal() + true_rho * t for t in range(n_years)] series.append(row) keep_n = max(8, int(round(n_sites * (1.0 - sel_strength)))) order = sorted(range(n_sites), key=lambda i: series[i][0], reverse=True) kept = order[:keep_n] ts = list(range(n_years)) mt = sum(ts) / n_years stt = sum((t - mt) ** 2 for t in ts) slopes = [] for i in kept: y = series[i] my = sum(y) / n_years sty = sum((t - mt) * (v - my) for t, v in zip(ts, y)) slopes.append(sty / stt) mean_slope = sum(slopes) / len(slopes) sd = math.sqrt(sum((s - mean_slope) ** 2 for s in slopes) / max(1, len(slopes) - 1)) return dict(sel=sel_strength, kept=keep_n, slope=mean_slope, se=sd / math.sqrt(len(slopes)), pct_yr=pct_yr(mean_slope), pct_dec=pct_dec(mean_slope), total=100.0 * (math.exp(mean_slope * (n_years - 1)) - 1.0)) # --------------------------------------------------------------------------- # report # --------------------------------------------------------------------------- def rule(ch="-", n=78): return ch * n def main(): print(rule("=")) print("THE INSECT DECLINE LITERATURE DISAGREES WITH ITSELF") print("Science Journaling Club -- random-effects pooling and diagnostics") print("Effect scale: rho = natural-log annual rate of change in total") print("abundance or biomass. percent per year = 100*(exp(rho)-1).") print(rule("=")) print() # ---- 0. validation ---------------------------------------------------- print("PART 0 -- VALIDATING THE MACHINERY BEFORE POINTING IT AT ANYTHING") print(rule()) print("Check 0a. Recovery of a known effect from synthetic data.") true_mu, true_tau = -0.025, 0.015 res0, ys0, ses0, coverage, mean_tau = validate_recovery(true_mu, true_tau) print(" true mu = %+.6f (%.3f%% per year)" % (true_mu, pct_yr(true_mu))) print(" true tau = %.6f" % true_tau) print(" k = %d synthetic studies" % res0["k"]) print(" pooled RE estimate = %+.6f 95%% CI [%+.6f, %+.6f]" % (res0["re"], res0["re_lo"], res0["re_hi"])) print(" tau-hat = %.6f (I2 = %.1f%%)" % (res0["tau"], res0["I2"])) inside = res0["re_lo"] <= true_mu <= res0["re_hi"] print(" truth inside the stated 95%% interval? %s" % ("YES" if inside else "NO")) print(" 2000-replicate coverage of the 95%% interval = %.1f%%" % (coverage * 100.0)) print(" 2000-replicate mean tau-hat = %.6f (true %.6f)" % (mean_tau, true_tau)) print(" [DL is known to under-cover a little when k is small and tau is") print(" large. High 80s to mid 90s is the documented behaviour of this") print(" estimator, not a coding error.]") print() print("Check 0b. Fixed and random effects must coincide when tau^2 = 0.") r0, y0, s0 = validate_tau_zero() print(" k = %d, all effects identical, so Q must be exactly 0." % r0["k"]) print(" Q = %.3e" % r0["Q"]) print(" tau^2 = %.3e" % r0["tau2"]) print(" FE estimate = %+.12f SE %.12f" % (r0["fe"], r0["fe_se"])) print(" RE estimate = %+.12f SE %.12f" % (r0["re"], r0["re_se"])) d_est = abs(r0["fe"] - r0["re"]) d_se = abs(r0["fe_se"] - r0["re_se"]) print(" |FE - RE| = %.3e |SE_FE - SE_RE| = %.3e" % (d_est, d_se)) ok = d_est < 1e-12 and d_se < 1e-12 print(" identical to machine precision? %s" % ("YES" if ok else "NO")) print() print(" Same check with dispersion small but not zero, so that tau^2 is") print(" clamped to 0 by the max(0, .) rather than landing there by") print(" construction:") s1 = [0.004, 0.009, 0.006, 0.012, 0.003, 0.007] y1 = [-0.0205, -0.0198, -0.0202, -0.0199, -0.0201, -0.0200] r1 = dersimonian_laird(y1, s1) print(" Q = %.4f on %d df -> tau^2 = %.3e" % (r1["Q"], r1["df"], r1["tau2"])) print(" |FE - RE| = %.3e" % abs(r1["fe"] - r1["re"])) print(" identical to machine precision? %s" % ("YES" if abs(r1["fe"] - r1["re"]) < 1e-12 else "NO")) print() # ---- 1. the extracted table ------------------------------------------ print("PART 1 -- THE EXTRACTED STUDIES") print(rule()) print("%-31s %-11s %10s %9s %8s" % ("study", "span", "rho", "SE", "%/yr")) for s in STUDIES: print("%-31s %-11s %+10.5f %9.5f %+8.2f" % (s["label"], s["span"].split(" ")[0], s["rho"], s["se"], pct_yr(s["rho"]))) print() print("Provenance of every number above:") for s in STUDIES: print(" [%s] %s" % (s["key"], s["source"])) print() ys = [s["rho"] for s in STUDIES] ses = [s["se"] for s in STUDIES] # ---- 2. pooling ------------------------------------------------------- print("PART 2 -- POOLING") print(rule()) res = dersimonian_laird(ys, ses) print(" k = %d effects" % res["k"]) print(" Fixed-effect estimate = %+.5f (%.2f%%/yr) SE %.5f" % (res["fe"], pct_yr(res["fe"]), res["fe_se"])) print(" Cochran's Q = %.2f on %d df p = %.3e" % (res["Q"], res["df"], res["Q_p"])) print(" DerSimonian-Laird tau^2 = %.6f tau = %.5f" % (res["tau2"], res["tau"])) print(" I-squared = %.1f%%" % res["I2"]) print(" H-squared (Q/df) = %.2f" % res["H2"]) print() print(" RANDOM-EFFECTS POOLED ESTIMATE") print(" rho = %+.5f 95%% CI [%+.5f, %+.5f]" % (res["re"], res["re_lo"], res["re_hi"])) print(" per year = %+.2f%% 95%% CI [%+.2f%%, %+.2f%%]" % (pct_yr(res["re"]), pct_yr(res["re_lo"]), pct_yr(res["re_hi"]))) print(" per decade = %+.1f%% 95%% CI [%+.1f%%, %+.1f%%]" % (pct_dec(res["re"]), pct_dec(res["re_lo"]), pct_dec(res["re_hi"]))) print(" z = %.2f, p = %.2e" % (res["z"], res["p"])) print() print(" 95% PREDICTION INTERVAL for the next comparable study") print(" rho in [%+.5f, %+.5f] = [%+.2f%%, %+.2f%%] per year" % (res["pi_lo"], res["pi_hi"], pct_yr(res["pi_lo"]), pct_yr(res["pi_hi"]))) print(" [The confidence interval describes the mean. The prediction") print(" interval describes a study. They are different objects and") print(" the second is the honest one to quote at a reader who wants") print(" to know what their local woodland is doing.]") print() print(" Random-effects weights, percent of total:") for s, w in zip(STUDIES, res["weights"]): print(" %-31s %6.2f%%" % (s["label"], 100.0 * w / res["swstar"])) print() print(" For contrast, the FIXED-effect weights, which is why nobody") print(" should be using a fixed-effect model on this literature:") _, _, wfix = fixed_effect(ys, ses) swf = sum(wfix) for s, w in zip(STUDIES, wfix): print(" %-31s %6.2f%%" % (s["label"], 100.0 * w / swf)) print() # ---- 3. funnel and Egger --------------------------------------------- print("PART 3 -- SMALL-STUDY ASYMMETRY") print(rule()) egg = eggers_test(ys, ses) print(" Egger's regression intercept = %+.4f (SE %.4f)" % (egg["intercept"], egg["se"])) print(" t = %+.3f on %d df, p = %.4f" % (egg["t"], egg["df"], egg["p"])) print(" regression slope = %+.5f (%+.2f%%/yr): the effect this test" % (egg["slope"], pct_yr(egg["slope"]))) print(" extrapolates to a study of infinite precision.") print(" [A significant intercept says the funnel is lopsided. It does") print(" NOT say why. With k=11 and one study an order of magnitude") print(" more precise than the rest, this test has almost no power and") print(" almost no interpretability. We report it because we said we") print(" would, and then we do not lean on it.]") print() print(" Funnel coordinates, most precise first:") print(" %-31s %10s %10s %8s" % ("study", "rho", "SE", "1/SE")) for s in sorted(STUDIES, key=lambda r: r["se"]): print(" %-31s %+10.5f %10.5f %8.1f" % (s["label"], s["rho"], s["se"], 1.0 / s["se"])) print() # ---- 4. leave one out ------------------------------------------------- print("PART 4 -- LEAVE-ONE-OUT") print(rule()) print(" %-31s %10s %22s %8s" % ("omitted", "rho", "95% CI", "I2")) loo = leave_one_out(STUDIES) for label, r in loo: print(" %-31s %+10.5f [%+.5f,%+.5f] %7.1f%%" % (label, r["re"], r["re_lo"], r["re_hi"], r["I2"])) lows = min(r["re"] for _, r in loo) highs = max(r["re"] for _, r in loo) print() print(" range of the pooled estimate across omissions: %+.5f to %+.5f" % (lows, highs)) print(" i.e. %.2f%%/yr to %.2f%%/yr" % (pct_yr(lows), pct_yr(highs))) print(" No single study moves the pooled estimate outside the interval") print(" reported in Part 2." if (res["re_lo"] <= lows and highs <= res["re_hi"]) else " At least one omission moves the estimate outside the " "interval reported in Part 2.") print() # ---- 5. subgroup sanity checks --------------------------------------- print("PART 5 -- SENSITIVITY TO THE DEPENDENCE WE KNOW WE HAVE") print(rule()) print(" Two papers contribute more than one effect. Shortall 2009 gives") print(" four traps; Hallmann 2020 gives three datasets from two") print(" localities. Re-pool keeping one effect per PAPER, then one per") print(" LOCALITY, then split by metric.") print() one_per_paper = ["hallmann17", "shortall_here", "hallmann20_lep", "wepprich19", "edwards25", "weiss24"] sub = [s for s in STUDIES if s["key"] in one_per_paper] r = dersimonian_laird([s["rho"] for s in sub], [s["se"] for s in sub]) print(" one per paper (k=%2d): rho %+.5f [%+.5f,%+.5f] = %+.2f%%/yr" " I2 %.1f%%" % (r["k"], r["re"], r["re_lo"], r["re_hi"], pct_yr(r["re"]), r["I2"])) one_per_loc = ["hallmann17", "shortall_here", "shortall_roth", "shortall_star", "shortall_wye", "hallmann20_lep", "hallmann20_car", "wepprich19", "edwards25", "weiss24"] sub2 = [s for s in STUDIES if s["key"] in one_per_loc] r2 = dersimonian_laird([s["rho"] for s in sub2], [s["se"] for s in sub2]) print(" one per locality (k=%2d): rho %+.5f [%+.5f,%+.5f] = %+.2f%%/yr" " I2 %.1f%%" % (r2["k"], r2["re"], r2["re_lo"], r2["re_hi"], pct_yr(r2["re"]), r2["I2"])) biomass = ["hallmann17", "shortall_here", "shortall_roth", "shortall_star", "shortall_wye"] sub3 = [s for s in STUDIES if s["key"] in biomass] r3 = dersimonian_laird([s["rho"] for s in sub3], [s["se"] for s in sub3]) sub4 = [s for s in STUDIES if s["key"] not in biomass] r4 = dersimonian_laird([s["rho"] for s in sub4], [s["se"] for s in sub4]) print(" biomass only (k=%2d): rho %+.5f [%+.5f,%+.5f] = %+.2f%%/yr" % (r3["k"], r3["re"], r3["re_lo"], r3["re_hi"], pct_yr(r3["re"]))) print(" abundance only (k=%2d): rho %+.5f [%+.5f,%+.5f] = %+.2f%%/yr" % (r4["k"], r4["re"], r4["re_lo"], r4["re_hi"], pct_yr(r4["re"]))) europe = ["hallmann17", "shortall_here", "shortall_roth", "shortall_star", "shortall_wye", "hallmann20_lep", "hallmann20_col", "hallmann20_car", "weiss24"] sub6 = [s for s in STUDIES if s["key"] in europe] sub7 = [s for s in STUDIES if s["key"] not in europe] r6 = dersimonian_laird([s["rho"] for s in sub6], [s["se"] for s in sub6]) r7 = dersimonian_laird([s["rho"] for s in sub7], [s["se"] for s in sub7]) print(" Europe only (k=%2d): rho %+.5f [%+.5f,%+.5f] = %+.2f%%/yr" % (r6["k"], r6["re"], r6["re_lo"], r6["re_hi"], pct_yr(r6["re"]))) print(" North America (k=%2d): rho %+.5f [%+.5f,%+.5f] = %+.2f%%/yr" % (r7["k"], r7["re"], r7["re_lo"], r7["re_hi"], pct_yr(r7["re"]))) print() fw = 100.0 * (1 / 0.0020 ** 2) / sum(1 / s ** 2 for s in ses) print(" Excluding the single most precise study (Hallmann 2017, SE") print(" 0.0020, which carries %.1f%% of the FIXED-effect weight and" % fw) print(" %.1f%% of the random-effects weight):" % (100.0 * res["weights"][0] / res["swstar"])) sub5 = [s for s in STUDIES if s["key"] != "hallmann17"] r5 = dersimonian_laird([s["rho"] for s in sub5], [s["se"] for s in sub5]) print(" rho %+.5f [%+.5f,%+.5f] = %+.2f%%/yr, I2 %.1f%%" % (r5["re"], r5["re_lo"], r5["re_hi"], pct_yr(r5["re"]), r5["I2"])) print() # ---- 6. external benchmarks ------------------------------------------ print("PART 6 -- HOW THE POOL COMPARES WITH THE PUBLISHED SYNTHESES") print(rule()) print(" These are NOT in the pool. They are syntheses of overlapping") print(" primary data; including them would double-count. They are the") print(" external check on whether our number is sane.") print() bench = [ ("van Klink 2020 terrestrial, as published", -0.92), ("van Klink 2020 terrestrial, erratum", -1.11), ("van Klink 2020 terrestrial, outliers removed", -0.66), ("van Klink 2020 terrestrial, North America out", -0.49), ("van Klink 2020 freshwater", 1.08), ("Haase 2023 European freshwater abundance", 1.17), ("Sockman 2025 subalpine Colorado, one meadow", -6.6), ] print(" %-47s %8s %9s" % ("synthesis or benchmark", "%/yr", "%/decade")) for name, p in bench: rho = math.log(1 + p / 100.0) print(" %-47s %+8.2f %+9.1f" % (name, p, pct_dec(rho))) print(" %-47s %+8.2f %+9.1f" % ("THIS POOL, k=%d, random effects" % res["k"], pct_yr(res["re"]), pct_dec(res["re"]))) print() vk = math.log(1 - 0.0111) print(" Difference between our pool and van Klink's corrected figure:") print(" %+.5f in log units, a factor of %.2f on the annual rate" % (res["re"] - vk, res["re"] / vk)) print(" Over 30 years: %.0f%% loss on our number versus %.0f%% loss" % (-100 * (math.exp(30 * res["re"]) - 1), -100 * (math.exp(30 * vk) - 1))) print(" on theirs. Same literature. Different rooms.") print() print(" Is van Klink's corrected estimate inside our 95%% CI? %s" % ("YES" if res["re_lo"] <= vk <= res["re_hi"] else "NO")) print(" Is it inside our 95%% PREDICTION interval? %s" % ("YES" if res["pi_lo"] <= vk <= res["pi_hi"] else "NO")) print() # ---- 7. selection bias ------------------------------------------------ print("PART 7 -- SITE SELECTION BIAS WITH NO REAL TREND UNDERNEATH") print(rule()) print(" 400 sites, 27 years, permanent site quality SD 1.0, annual noise") print(" SD 0.6, TRUE trend exactly zero. Sites are chosen on their first") print(" year's observed abundance. Mean fitted log-linear slope:") print() print(" %-14s %6s %12s %9s %13s" % ("kept", "n", "slope", "%/yr", "27-yr total")) for sel in [0.0, 0.25, 0.50, 0.75, 0.90, 0.95, 0.975]: s = selection_bias_sim(sel) print(" top %-10.1f%% %6d %+12.5f %+9.2f %+12.1f%%" % (100 * (1 - sel), s["kept"], s["slope"], s["pct_yr"], s["total"])) print() print(" Repeat with a REAL decline of 1%/yr underneath, to see how much") print(" selection adds on top of a true signal:") print(" %-14s %6s %12s %9s %13s" % ("kept", "n", "slope", "%/yr", "27-yr total")) for sel in [0.0, 0.50, 0.90, 0.975]: s = selection_bias_sim(sel, true_rho=math.log(0.99)) print(" top %-10.1f%% %6d %+12.5f %+9.2f %+12.1f%%" % (100 * (1 - sel), s["kept"], s["slope"], s["pct_yr"], s["total"])) print() print(" And again with year-to-year noise turned down to SD 0.2, which") print(" is what you would have if traps were quiet and weather did not") print(" matter. Regression to the mean needs noise to work with:") print(" %-14s %6s %12s %9s" % ("kept", "n", "slope", "%/yr")) for sel in [0.0, 0.90, 0.975]: s = selection_bias_sim(sel, year_sd=0.2) print(" top %-10.1f%% %6d %+12.5f %+9.2f" % (100 * (1 - sel), s["kept"], s["slope"], s["pct_yr"])) print() print(" And once more with the noise turned UP to SD 1.0, which is what") print(" a light trap in a bad summer looks like:") print(" %-14s %6s %12s %9s %13s" % ("kept", "n", "slope", "%/yr", "27-yr total")) for sel in [0.0, 0.50, 0.90, 0.975]: s = selection_bias_sim(sel, year_sd=1.0) print(" top %-10.1f%% %6d %+12.5f %+9.2f %+12.1f%%" % (100 * (1 - sel), s["kept"], s["slope"], s["pct_yr"], s["total"])) print() print(" Full grid, for the figure in the article. Rows are the fraction") print(" of sites kept; columns are the year-to-year noise SD. Entries") print(" are the mean fitted trend in percent per year, true trend zero.") print() print(" Each entry is averaged over 24 independent seeds, because a") print(" single run of 10 surviving sites is itself a noisy estimate and") print(" we are not going to draw a wobble and call it a mechanism.") print() keeps = [1.0, 0.75, 0.50, 0.35, 0.25, 0.15, 0.10, 0.05, 0.025] print(" %-10s %10s %10s %10s" % ("kept", "sd=0.2", "sd=0.6", "sd=1.0")) for kp in keeps: row = [] for nsd in (0.2, 0.6, 1.0): vals = [selection_bias_sim(1.0 - kp, year_sd=nsd, seed=1 + 37 * j) ["pct_yr"] for j in range(24)] row.append(sum(vals) / len(vals)) print(" %-10s %+10.3f %+10.3f %+10.3f" % ("%.1f%%" % (100 * kp), row[0], row[1], row[2])) print() print(" How much of OUR pooled estimate could this account for?") base = selection_bias_sim(0.90) print(" pooled estimate %+.5f (%.2f%% per year)" % (res["re"], pct_yr(res["re"]))) print(" top-10%% selection %+.5f (%.2f%% per year)" % (base["slope"], base["pct_yr"])) print(" share explained %.1f%%" % (100.0 * base["slope"] / res["re"])) print(" [That share is an upper bound under a deliberately unkind") print(" model. It is not nothing, and it is not the whole thing.]") print() print(rule("=")) print("END. Every printed number is reproducible by running this file.") print(rule("=")) if __name__ == "__main__": main()