""" edna-detection-review.py Science Journaling Club - systematic review / random-effects meta-analysis QUESTION Across published studies that surveyed the same water bodies with both environmental DNA and a conventional method, how much better (or worse) is one eDNA survey unit than one conventional survey unit at detecting a species that is genuinely there? EFFECT SIZE log odds ratio of per-unit detection probability, eDNA relative to the conventional comparator, conditional on the site being occupied. This is exactly the p of a MacKenzie-style occupancy model, on the logit scale, differenced between two methods. y_i = logit(p_eDNA) - logit(p_conv) v_i = SE[logit(p_eDNA)]^2 + SE[logit(p_conv)]^2 y > 0 means eDNA wins per unit of effort. Every number in DATA below is traceable to a published source, recorded in `source` on each row. MODEL DerSimonian-Laird random effects, with Cochran's Q, I-squared, a prediction interval, Egger's regression test for small-study asymmetry, leave-one-out sensitivity, a subgroup split by survey unit, and a meta-regression on log10 water volume. VALIDATION Before any of that runs on real data, two checks: (1) recovery of a known true effect from synthetic studies, (2) the algebraic identity that fixed-effect and random-effects estimates coincide when tau^2 comes out at zero. Run: python edna-detection-review.py > edna-detection-review-output.txt """ import math import sys import random # --------------------------------------------------------------------------- # small numeric helpers (no scipy: keep the file readable and portable) # --------------------------------------------------------------------------- def logit(p): return math.log(p / (1.0 - p)) def inv_logit(x): return 1.0 / (1.0 + math.exp(-x)) def norm_cdf(z): return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0))) def two_sided_p(z): return 2.0 * (1.0 - norm_cdf(abs(z))) def chi2_sf(x, k): """Upper tail of a chi-square with k degrees of freedom. Exact series for even k, and the erfc-based series for odd k. Good to better than 1e-12 for the sizes used here (k <= 40). """ if x <= 0: return 1.0 if k % 2 == 0: term = math.exp(-x / 2.0) total = term for i in range(1, k // 2): term *= (x / 2.0) / i total += term return min(1.0, total) else: s = math.sqrt(x) total = math.erfc(s / math.sqrt(2.0)) term = math.exp(-x / 2.0) * math.sqrt(2.0 / math.pi) * s for i in range(1, (k + 1) // 2): total += term term *= x / (2.0 * i + 1.0) return min(1.0, total) def t_sf(t, df): """Upper tail of Student's t. Continued-fraction free: uses the incomplete beta via a simple continued fraction (Lentz).""" x = df / (df + t * t) ib = _betainc(df / 2.0, 0.5, x) p = 0.5 * ib return p if t > 0 else 1.0 - p def _betacf(a, b, x): 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, 300): 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 de = d * c h *= de if abs(de - 1.0) < 3e-14: 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) + a * math.log(x) + b * math.log(1.0 - x)) front = math.exp(lbeta) if x < (a + 1.0) / (a + b + 2.0): return front * _betacf(a, b, x) / a return 1.0 - math.exp( math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) + b * math.log(1.0 - x) + a * math.log(x)) * _betacf(b, a, 1.0 - x) / b Z95 = 1.959963985 # --------------------------------------------------------------------------- # converting whatever a paper reported into logit p and its standard error # --------------------------------------------------------------------------- def from_ci(p, lo, hi): """Point estimate with a 95% interval, all on the probability scale.""" return logit(p), (logit(hi) - logit(lo)) / (2.0 * Z95) def from_se(p, se_p): """Point estimate with an SE on the probability scale (delta method).""" return logit(p), se_p / (p * (1.0 - p)) def from_n(p, n): """Point estimate with no reported uncertainty. We reconstruct a binomial SE from the number of independent survey units the paper reports, which is conservative because it ignores replicate sub-samples within a unit.""" return logit(p), math.sqrt(1.0 / (n * p * (1.0 - p))) def from_counts(a, n1, c, n2): """Detections out of survey units, both methods. Haldane-Anscombe 0.5 correction applied throughout so that zero cells stay finite.""" b, d = n1 - a, n2 - c a_, b_, c_, d_ = a + 0.5, b + 0.5, c + 0.5, d + 0.5 y = math.log((a_ * d_) / (b_ * c_)) v = 1.0 / a_ + 1.0 / b_ + 1.0 / c_ + 1.0 / d_ return y, v # --------------------------------------------------------------------------- # THE DATA # Each row: one study, one target species, one paired contrast. # Where a study reported several species, we take the MEDIAN contrast by log # odds ratio, decided before looking at the pooled result, so that no study # gets represented by whichever of its species flatters eDNA most. # --------------------------------------------------------------------------- DATA = [ dict( key="Biggs 2015", ref=8, species="Great crested newt (Triturus cristatus)", group="amphibian", unit="visit", volume_L=0.09, # 600 mL pooled, 90 mL actually precipitated pore="none (ethanol precipitation)", n=140, kind="counts", edna=(139, 140), conv=(105, 140), conv_name="torch count (median of 3 field methods)", source="Text: eDNA detected newts on 139 of 140 visits (99.3%); " "bottle traps 76%, torch counts 75%, egg searches 44%, over " "35 ponds visited 4 times.", ), dict( key="Valentini 2016", ref=7, species="Amphibian assemblage, 39 French water bodies", group="amphibian", unit="visit", volume_L=0.09, # Biggs protocol, ponds pore="none (ethanol precipitation)", n=39, kind="ci", edna=(0.97, 0.90, 0.99), conv=(0.58, 0.50, 0.63), conv_name="single traditional visit (netting, visual, calling)", source="Results: 'the detection probability with eDNA metabarcoding " "was 0.97 (CI = 0.90-0.99) vs. 0.58 (CI = 0.50-0.63) for " "traditional surveys.'", ), dict( key="Schmelzle 2016", ref=9, species="Tidewater goby (Eucyclogobius newberryi)", group="fish", unit="sample", volume_L=None, pore="not recorded by us", n=29, kind="derived", edna=(0.74, 29), conv=(0.39, 29), conv_name="seine haul", source="Abstract/summary: detection probability 0.74 for eDNA vs 0.39 " "for seining across 29 paired locations. No uncertainty " "reported in the text available to us; SE reconstructed from " "n = 29 sites.", ), dict( key="Hinlo 2017", ref=10, species="Redfin perch (Perca fluviatilis)", group="fish", unit="site-season", volume_L=12.0, # six 2 L samples per site pore="1.2 um glass fibre", n=14, kind="counts", edna=(3, 14), conv=(0, 14), conv_name="fyke net set", source="Table 1 detection matrix, 8 sites x 2 seasons, restricted to " "the 14 site-seasons where both methods ran. Median of three " "species contrasts (carp -1.17, redfin +2.18, weatherloach " "+4.36 on the log odds scale).", ), dict( key="Eiler 2018", ref=11, species="Pool frog (Pelophylax lessonae)", group="amphibian", unit="visit", volume_L=0.85, # 200-1500 mL, midpoint pore="0.45 um", n=49, kind="derived", edna=(0.38, 49), conv=(0.40, 49), conv_name="calling-male count plus visual", source="Results: 'detection probability for P. lessonae was 0.38 per " "sample with the eDNA method compared to 0.40 when detecting " "calling males and silent observed individuals'; 49 matched " "samples. SE reconstructed from n = 49.", ), dict( key="Rose 2019", ref=12, species="Northern watersnake (Nerodia sipedon)", group="reptile", unit="sample", volume_L=0.5, pore="0.45 um", n=61, kind="ci", edna=(0.44, 0.15, 0.75), conv=(0.74, 0.48, 0.95), conv_name="one trap-night (30 traps)", source="Reported per-unit detection with 95% credible intervals: eDNA " "0.44 (0.15-0.75), trapping 0.74 (0.48-0.95). Median of two " "species contrasts; N. sipedon has the larger site count (61).", ), dict( key="Akre 2019", ref=13, species="Wood turtle (Glyptemys insculpta)", group="reptile", unit="sample", volume_L=2.0, pore="0.45 um cellulose nitrate", n=37, kind="ci", edna=(0.55, 0.38, 0.71), conv=(0.88, 0.58, 0.98), conv_name="one-hour visual encounter survey of 1 km", source="Reported detection probabilities with 95% CI: eDNA per sample " "0.55 (0.38-0.71); visual encounter survey 0.88 (0.58-0.98); " "37 stream reaches.", ), dict( key="Moss 2022", ref=14, species="California newt (Taricha torosa)", group="amphibian", unit="visit", volume_L=0.25, # 500 mL split across two filters pore="0.45 or 5 um", n=20, kind="ci", edna=(0.65, 0.44, 0.81), conv=(0.69, 0.48, 0.84), conv_name="seine haul", source="Results: metabarcoding p = 0.65 (95% CI 0.44-0.81), seining " "p = 0.69 (0.48-0.84), 20 ponds. Median of the three species " "with full intervals reported (red-legged frog +3.24, newt " "-0.18, bullfrog -0.18).", ), dict( key="Quilumbaquin 2023", ref=15, species="Amphibian assemblage, Ecuadorian Amazon", group="amphibian", unit="visit", volume_L=1.0, pore="0.45 um nitrocellulose", n=36, kind="ci", edna=(0.42, 0.40, 0.45), conv=(0.17, 0.14, 0.20), conv_name="visual encounter survey", source="Abstract: 'eDNA detected 28 species and had a detection " "probability (DP) of 0.42 CI [0.40-0.45], while VES recorded " "20 species with a DP of 0.17 CI [0.14-0.20].'", ), dict( key="Li 2024", ref=16, species="Amphibian assemblage, Zhoushan Archipelago", group="amphibian", unit="site", volume_L=2.0, pore="1.5 um glass fibre", n=21, kind="derived", edna=(0.54, 21), conv=(0.24, 21), conv_name="traditional line transect method", source="Results: 'the mean detection probability of eDNA is 0.54, " "while the mean detection probability of TLTM is 0.24'; 21 " "islands. No uncertainty reported; SE reconstructed from " "n = 21 islands.", ), dict( key="Di Girolamo 2024", ref=17, species="American mink (Neogale vison)", group="mammal", unit="visit", volume_L=0.5, pore="0.45 um nitrocellulose", n=21, kind="se", edna=(0.25, 0.08), conv=(0.36, 0.16), conv_name="one camera-trap week", source="Results: camera traps rho = 0.36 (SE 0.16), eDNA rho = 0.25 " "(SE 0.08), 7 sites over 21 survey weeks.", ), dict( key="Dougherty 2025", ref=18, species="Alewife (Alosa pseudoharengus)", group="fish", unit="site-season", volume_L=3.55, # 10 replicates of 355 mL pore="1.2 um glass microfibre", n=7, kind="counts", edna=(4, 7), conv=(7, 7), conv_name="purse seine survey (6 sets)", source="Results: 'we observed a 51.7% detection rate with 3 false " "negatives out of a sample size of 7 lakes'; seining detected " "alewife in every lake where the species was present.", ), ] # Studies found, read, and deliberately NOT pooled, with the reason. EXCLUDED = [ ("Smart 2015 (ref 19)", "smooth newt, Melbourne", "reports per-site RANGES only (eDNA 0.29-1.0 per sample, traps 0.01-0.26 " "per trap) with no pooled estimate or interval. Nothing to pool."), ("Hunter 2015 (ref 20)", "Burmese python, Florida", "conventional effort returned zero detections in 5,935 trap-nights and " "zero standardised visual sightings. The odds ratio is not finite and a " "continuity correction would invent the comparison."), ("Lopes 2017 (ref 21)", "tropical stream anurans, Brazil", "eDNA detection probabilities only (no paired conventional p), but the " "20 L vs 60 L contrast is used below as direct evidence on volume."), ("Dougherty 2016 crayfish", "rusty crayfish, Wisconsin", "eDNA detection modelled against trap CPUE as a covariate; no detection " "probability estimated for trapping itself."), ("de Souza 2016", "waterdog and musk turtle, Alabama", "eDNA seasonal detection only; no conventional comparator fitted."), ("Lyet 2021", "terrestrial mammals via stream eDNA", "per-taxon eDNA detection given for 30-80 L samples, camera comparison " "expressed as taxa counts rather than a paired per-unit probability."), ("NCDOT RP2023-18", "Roanoke logperch, North Carolina", "grey literature: eDNA 6/8 vs electrofishing 1/8 at historically " "occupied sites. Excluded because it is not peer reviewed."), ] # Direct evidence on water volume, from a single study that varied it. # Lopes et al. 2017, per-sample detection p11 at 20 L and at 60 L. LOPES_VOLUME = [ ("Hylodes phyllodes", 0.614, 0.761), ("Hylodes asper", 0.570, 0.649), ("Cycloramphus boraceiensis", 0.154, 0.596), ] # --------------------------------------------------------------------------- # effect sizes # --------------------------------------------------------------------------- def effect(row): k = row["kind"] if k == "counts": a, n1 = row["edna"] c, n2 = row["conv"] return from_counts(a, n1, c, n2) if k == "ci": p, lo, hi = row["edna"] le, se_e = from_ci(p, lo, hi) p, lo, hi = row["conv"] lc, se_c = from_ci(p, lo, hi) elif k == "se": le, se_e = from_se(*row["edna"]) lc, se_c = from_se(*row["conv"]) elif k == "derived": le, se_e = from_n(*row["edna"]) lc, se_c = from_n(*row["conv"]) else: raise ValueError(k) return le - lc, se_e ** 2 + se_c ** 2 # --------------------------------------------------------------------------- # DerSimonian-Laird random effects # --------------------------------------------------------------------------- def fixed_effect(y, v): w = [1.0 / vi for vi in v] sw = sum(w) est = sum(wi * yi for wi, yi in zip(w, y)) / sw se = math.sqrt(1.0 / sw) return est, se, w def dersimonian_laird(y, v): k = len(y) fe, fe_se, w = fixed_effect(y, v) Q = sum(wi * (yi - fe) ** 2 for wi, yi in zip(w, y)) df = k - 1 sw = sum(w) sw2 = sum(wi ** 2 for wi in w) C = sw - sw2 / sw tau2 = max(0.0, (Q - df) / C) if C > 0 else 0.0 I2 = max(0.0, (Q - df) / Q) * 100.0 if Q > 0 else 0.0 H2 = Q / df if df > 0 else float("nan") wr = [1.0 / (vi + tau2) for vi in v] swr = sum(wr) est = sum(wi * yi for wi, yi in zip(wr, y)) / swr se = math.sqrt(1.0 / swr) z = est / se return dict(k=k, fe=fe, fe_se=fe_se, est=est, se=se, z=z, p=two_sided_p(z), tau2=tau2, tau=math.sqrt(tau2), Q=Q, df=df, Qp=chi2_sf(Q, df) if df > 0 else float("nan"), I2=I2, H2=H2, lo=est - Z95 * se, hi=est + Z95 * se) def prediction_interval(res): """Where a 13th study's true effect would be expected to land.""" if res["df"] < 1: return (float("nan"), float("nan")) t = 2.200985 # t(0.975, 11 df); recomputed below if k changes df = res["k"] - 2 t = _t_crit(df) half = t * math.sqrt(res["tau2"] + res["se"] ** 2) return res["est"] - half, res["est"] + half def _t_crit(df, alpha=0.025): lo, hi = 0.0, 100.0 for _ in range(200): mid = 0.5 * (lo + hi) if t_sf(mid, df) > alpha: lo = mid else: hi = mid return 0.5 * (lo + hi) def egger(y, v): """Classic Egger regression: standard normal deviate on precision. SND_i = y_i / se_i regressed on 1/se_i; the intercept is the test.""" snd = [yi / math.sqrt(vi) for yi, vi in zip(y, v)] prec = [1.0 / math.sqrt(vi) for vi in v] n = len(y) mx = sum(prec) / n my = sum(snd) / n sxx = sum((x - mx) ** 2 for x in prec) sxy = sum((x - mx) * (yy - my) for x, yy in zip(prec, snd)) slope = sxy / sxx intercept = my - slope * mx resid = [yy - (intercept + slope * x) for x, yy in zip(prec, snd)] s2 = sum(r * r for r in resid) / (n - 2) se_int = math.sqrt(s2 * (1.0 / n + mx * mx / sxx)) t = intercept / se_int return dict(intercept=intercept, se=se_int, t=t, df=n - 2, p=2.0 * t_sf(abs(t), n - 2), slope=slope) def meta_regress(y, v, x, tau2): """Weighted least squares of y on a single moderator, random effects weights 1/(v + tau2).""" w = [1.0 / (vi + tau2) for vi in v] sw = sum(w) mx = sum(wi * xi for wi, xi in zip(w, x)) / sw my = sum(wi * yi for wi, yi in zip(w, y)) / sw sxx = sum(wi * (xi - mx) ** 2 for wi, xi in zip(w, x)) sxy = sum(wi * (xi - mx) * (yi - my) for wi, xi, yi in zip(w, x, y)) slope = sxy / sxx intercept = my - slope * mx se_slope = math.sqrt(1.0 / sxx) z = slope / se_slope return dict(slope=slope, se=se_slope, z=z, p=two_sided_p(z), intercept=intercept) # --------------------------------------------------------------------------- # reporting helpers # --------------------------------------------------------------------------- def rule(ch="-", n=78): print(ch * n) def head(title): print() rule("=") print(title) rule("=") def show(res, label, as_or=True): print(" %-22s %7.3f [%7.3f, %7.3f]" % (label, res["est"], res["lo"], res["hi"])) if as_or: print(" %-22s %7.3f [%7.3f, %7.3f]" % (" as an odds ratio", math.exp(res["est"]), math.exp(res["lo"]), math.exp(res["hi"]))) # --------------------------------------------------------------------------- # PART 0. VALIDATE THE MACHINERY BEFORE TRUSTING IT # --------------------------------------------------------------------------- def validate(): head("PART 0. VALIDATION OF THE POOLING CODE") print(""" Check 1. Recovery of a known true effect. We simulate 14 studies drawn from a random-effects world with a true mean log odds ratio of MU = 0.800 and a true between-study SD of TAU = 0.450. Each study gets its own within-study variance. If the code is right, the pooled estimate should sit near 0.800 and its 95% interval should contain it in roughly 95% of repeats. """) random.seed(20260913) MU, TAU = 0.800, 0.450 ses = [0.18, 0.22, 0.25, 0.30, 0.33, 0.36, 0.40, 0.45, 0.50, 0.55, 0.62, 0.70, 0.80, 0.95] y, v = [], [] for s in ses: theta = random.gauss(MU, TAU) y.append(random.gauss(theta, s)) v.append(s * s) res = dersimonian_laird(y, v) print(" single synthetic data set, k = %d" % res["k"]) show(res, "pooled log OR", as_or=False) print(" true value MU %7.3f" % MU) inside = res["lo"] <= MU <= res["hi"] print(" interval contains MU: %s" % ("YES" if inside else "NO")) print(" tau-hat %.3f against true TAU %.3f" % (res["tau"], TAU)) print(" I-squared %.1f%%, Q = %.2f on %d df (p = %.4f)" % (res["I2"], res["Q"], res["df"], res["Qp"])) hits, reps, taus, ests = 0, 4000, [], [] for _ in range(reps): yy, vv = [], [] for s in ses: theta = random.gauss(MU, TAU) yy.append(random.gauss(theta, s)) vv.append(s * s) r = dersimonian_laird(yy, vv) ests.append(r["est"]) taus.append(r["tau2"]) if r["lo"] <= MU <= r["hi"]: hits += 1 cov = 100.0 * hits / reps mean_est = sum(ests) / reps print() print(" %d repeats:" % reps) print(" mean pooled estimate %7.4f (true %.4f, bias %+.4f)" % (mean_est, MU, mean_est - MU)) print(" mean tau-squared %7.4f (true %.4f)" % (sum(taus) / reps, TAU * TAU)) print(" 95%% interval coverage %6.1f%% (nominal 95.0%%)" % cov) ok1 = abs(mean_est - MU) < 0.02 and 91.0 < cov < 97.5 print(" VERDICT: %s" % ("PASS" if ok1 else "FAIL")) print(""" Check 2. The tau-squared = 0 identity. When Q <= df the DerSimonian-Laird estimator truncates tau-squared to zero, and the random-effects weights 1/(v + 0) become the fixed-effect weights 1/v exactly. The two estimates and their standard errors must then be identical, not merely close. We build a data set with no real heterogeneity and check it to machine precision. """) y2 = [0.30, 0.32, 0.28, 0.31, 0.29, 0.30, 0.31] v2 = [0.09, 0.16, 0.12, 0.10, 0.14, 0.11, 0.13] r2 = dersimonian_laird(y2, v2) print(" Q = %.4f on %d df, tau-squared = %.6f" % (r2["Q"], r2["df"], r2["tau2"])) print(" fixed effect %.12f (SE %.12f)" % (r2["fe"], r2["fe_se"])) print(" random effects %.12f (SE %.12f)" % (r2["est"], r2["se"])) d1 = abs(r2["fe"] - r2["est"]) d2 = abs(r2["fe_se"] - r2["se"]) print(" absolute differences: estimate %.3e, SE %.3e" % (d1, d2)) ok2 = r2["tau2"] == 0.0 and d1 < 1e-12 and d2 < 1e-12 print(" VERDICT: %s" % ("PASS" if ok2 else "FAIL")) print(""" Check 3. Egger's test on symmetric data. Feed the regression a funnel that is symmetric by construction. The intercept should be statistically indistinguishable from zero. """) random.seed(7) ys, vs = [], [] for s in [0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.75, 0.9, 1.1]: ys.append(random.gauss(0.4, s)) vs.append(s * s) eg = egger(ys, vs) print(" intercept %+.3f (SE %.3f), t = %+.2f on %d df, p = %.3f" % (eg["intercept"], eg["se"], eg["t"], eg["df"], eg["p"])) ok3 = eg["p"] > 0.05 print(" VERDICT: %s" % ("PASS" if ok3 else "FAIL")) print() print(" ALL VALIDATION CHECKS: %s" % ("PASS" if (ok1 and ok2 and ok3) else "FAIL")) return ok1 and ok2 and ok3 # --------------------------------------------------------------------------- # main # --------------------------------------------------------------------------- def main(): print("=" * 78) print("eDNA versus conventional survey detection: a club meta-analysis") print("Science Journaling Club | analysis/edna-detection-review.py") print("=" * 78) validate() # ---- the extracted table ------------------------------------------- head("PART 1. THE EXTRACTED STUDIES") print() print(" %-18s %-34s %8s %8s %9s" % ("study", "target", "p(eDNA)", "p(conv)", "log OR")) rule() rows = [] for r in DATA: y, v = effect(r) r["y"], r["v"] = y, v r["se"] = math.sqrt(v) if r["kind"] == "counts": pe = r["edna"][0] / r["edna"][1] pc = r["conv"][0] / r["conv"][1] else: pe = r["edna"][0] pc = r["conv"][0] r["pe"], r["pc"] = pe, pc rows.append(r) print(" %-18s %-34s %8.2f %8.2f %+9.3f" % (r["key"], r["species"][:34], pe, pc, y)) rule() print() print(" Per-study detail, with the uncertainty we used and where it came from") for r in rows: print() print(" %s [ref %d]" % (r["key"], r["ref"])) print(" target %s" % r["species"]) print(" comparator %s" % r["conv_name"]) print(" survey unit %s; n = %d; volume %s L; filter %s" % (r["unit"], r["n"], ("%.2f" % r["volume_L"]) if r["volume_L"] else "not recorded", r["pore"])) print(" y = %+.3f, var = %.4f, SE = %.3f (input type: %s)" % (r["y"], r["v"], r["se"], r["kind"])) print(" source %s" % r["source"]) print() print(" Studies read and deliberately excluded from pooling:") for name, what, why in EXCLUDED: print(" - %s, %s" % (name, what)) print(" %s" % why) y = [r["y"] for r in rows] v = [r["v"] for r in rows] # ---- the pooled result --------------------------------------------- head("PART 2. POOLED RANDOM-EFFECTS ESTIMATE (DerSimonian-Laird)") res = dersimonian_laird(y, v) print() print(" k = %d contrasts from %d studies" % (res["k"], res["k"])) show(res, "random effects") print(" z = %+.3f, p = %.4f" % (res["z"], res["p"])) print() print(" fixed effect %7.3f [%7.3f, %7.3f]" % (res["fe"], res["fe"] - Z95 * res["fe_se"], res["fe"] + Z95 * res["fe_se"])) print() print(" Cochran's Q %7.2f on %d df, p = %.5f" % (res["Q"], res["df"], res["Qp"])) print(" tau-squared %7.4f (tau = %.3f on the log odds scale)" % (res["tau2"], res["tau"])) print(" I-squared %6.1f%%" % res["I2"]) print(" H-squared %7.2f" % res["H2"]) pi_lo, pi_hi = prediction_interval(res) print(" 95%% prediction interval [%7.3f, %7.3f] (odds ratio %.2f to %.2f)" % (pi_lo, pi_hi, math.exp(pi_lo), math.exp(pi_hi))) print() print(" Read on the probability scale: if one conventional survey has a") print(" detection probability of 0.50, the pooled odds ratio puts one eDNA") print(" survey unit at %.2f, and the prediction interval runs %.2f to %.2f." % (inv_logit(logit(0.5) + res["est"]), inv_logit(logit(0.5) + pi_lo), inv_logit(logit(0.5) + pi_hi))) # ---- small-study bias ---------------------------------------------- head("PART 3. SMALL-STUDY BIAS: FUNNEL AND EGGER") eg = egger(y, v) print() print(" Egger regression of the standard normal deviate on precision") print(" intercept %+.3f (SE %.3f)" % (eg["intercept"], eg["se"])) print(" t %+.3f on %d df" % (eg["t"], eg["df"])) print(" p %.4f" % eg["p"]) print(" slope %+.3f" % eg["slope"]) print() print(" Funnel coordinates (x = log OR, y = standard error):") print(" %-18s %9s %9s" % ("study", "log OR", "SE")) for r in sorted(rows, key=lambda q: q["se"]): print(" %-18s %+9.3f %9.3f" % (r["key"], r["y"], r["se"])) semax = max(r["se"] for r in rows) print() print(" Pseudo 95%% funnel at SE = %.2f spans %+.3f to %+.3f around the pooled mean" % (semax, res["est"] - Z95 * semax, res["est"] + Z95 * semax)) # ---- leave one out -------------------------------------------------- head("PART 4. LEAVE-ONE-OUT SENSITIVITY") print() print(" %-18s %9s %19s %8s %8s" % ("omitted", "pooled", "95% CI", "I2 (%)", "tau2")) rule() loo = [] for i in range(len(rows)): yy = [y[j] for j in range(len(y)) if j != i] vv = [v[j] for j in range(len(v)) if j != i] r = dersimonian_laird(yy, vv) loo.append((rows[i]["key"], r)) print(" %-18s %+9.3f [%+6.3f, %+6.3f] %8.1f %8.4f" % (rows[i]["key"], r["est"], r["lo"], r["hi"], r["I2"], r["tau2"])) rule() lo_est = min(r["est"] for _, r in loo) hi_est = max(r["est"] for _, r in loo) print(" range of the pooled estimate under leave-one-out: %+.3f to %+.3f" % (lo_est, hi_est)) print(" as odds ratios: %.2f to %.2f" % (math.exp(lo_est), math.exp(hi_est))) sign_flips = sum(1 for _, r in loo if (r["est"] > 0) != (res["est"] > 0)) crosses = sum(1 for _, r in loo if r["lo"] <= 0 <= r["hi"]) print(" sign changes: %d of %d; intervals containing zero: %d of %d" % (sign_flips, len(loo), crosses, len(loo))) # ---- subgroups ------------------------------------------------------ head("PART 5. SUBGROUPS AND MODERATORS") print() print(" (a) By survey unit. A 'sample' contrast compares one bottle of") print(" water with one net haul. A 'visit' or 'site-season' contrast") print(" compares whole survey rounds, which usually bundle several") print(" water samples against a single conventional effort.") print() groups = {} for r in rows: groups.setdefault("sample" if r["unit"] == "sample" else "round", []).append(r) for name, g in sorted(groups.items()): rg = dersimonian_laird([q["y"] for q in g], [q["v"] for q in g]) print(" %-8s k = %2d pooled %+.3f [%+.3f, %+.3f] OR %.2f I2 %.1f%%" % (name, rg["k"], rg["est"], rg["lo"], rg["hi"], math.exp(rg["est"]), rg["I2"])) if len(groups) == 2: (n1, g1), (n2, g2) = sorted(groups.items()) r1 = dersimonian_laird([q["y"] for q in g1], [q["v"] for q in g1]) r2 = dersimonian_laird([q["y"] for q in g2], [q["v"] for q in g2]) d = r2["est"] - r1["est"] sd = math.sqrt(r1["se"] ** 2 + r2["se"] ** 2) print(" difference (%s minus %s) %+.3f (SE %.3f), z = %+.2f, p = %.3f" % (n2, n1, d, sd, d / sd, two_sided_p(d / sd))) print() print(" (b) By broad taxon.") tg = {} for r in rows: tg.setdefault(r["group"], []).append(r) for name, g in sorted(tg.items()): if len(g) == 1: print(" %-10s k = 1 single contrast %+.3f (SE %.3f)" % (name, g[0]["y"], g[0]["se"])) else: rg = dersimonian_laird([q["y"] for q in g], [q["v"] for q in g]) print(" %-10s k = %2d pooled %+.3f [%+.3f, %+.3f] I2 %.1f%%" % (name, rg["k"], rg["est"], rg["lo"], rg["hi"], rg["I2"])) print() print(" (c) Meta-regression on log10 of the water volume processed,") print(" restricted to the %d contrasts that report it." % sum(1 for r in rows if r["volume_L"])) sub = [r for r in rows if r["volume_L"]] mr = meta_regress([r["y"] for r in sub], [r["v"] for r in sub], [math.log10(r["volume_L"]) for r in sub], res["tau2"]) print(" slope %+.3f per decade of litres (SE %.3f), z = %+.2f, p = %.3f" % (mr["slope"], mr["se"], mr["z"], mr["p"])) print(" intercept %+.3f" % mr["intercept"]) print(" Volumes in the set span %.2f L to %.2f L, a factor of %.0f." % (min(r["volume_L"] for r in sub), max(r["volume_L"] for r in sub), max(r["volume_L"] for r in sub) / min(r["volume_L"] for r in sub))) print() print(" (d) Within-study evidence on volume, from the one study that") print(" varied it experimentally (Lopes 2017, 20 L vs 60 L samples):") print(" %-32s %8s %8s %10s" % ("species", "p(20 L)", "p(60 L)", "log OR")) tot = [] for sp, p20, p60 in LOPES_VOLUME: lor = logit(p60) - logit(p20) tot.append(lor) print(" %-32s %8.3f %8.3f %+10.3f" % (sp, p20, p60, lor)) print(" mean log OR for tripling volume: %+.3f (odds ratio %.2f)" % (sum(tot) / len(tot), math.exp(sum(tot) / len(tot)))) # ---- false positives ------------------------------------------------ head("PART 6. THE FALSE POSITIVE PROBLEM") print(""" Everything above treats a positive as truth. It is not. Write p11 for the probability that an occupied site gives a positive (true positive rate) and p10 for the probability that an unoccupied site gives one anyway, through contamination, primer cross-reactivity, DNA transported from upstream, or a tag-jump in the sequencing run. With prior occupancy psi, one positive sample carries a posterior probability of occupancy PPV = psi*p11 / ( psi*p11 + (1-psi)*p10 ) The method's whole selling point is rare species. Rare means small psi. And PPV falls off a cliff as psi falls, at any fixed p10. """) p11 = 0.60 print(" With p11 = %.2f, one positive sample:" % p11) print() print(" %-10s" % "psi", end="") p10s = [0.001, 0.005, 0.01, 0.02, 0.05] for q in p10s: print(" p10=%-6.3f" % q, end="") print() for psi in [0.50, 0.20, 0.10, 0.05, 0.02, 0.01, 0.005]: print(" %-10.3f" % psi, end="") for p10 in p10s: ppv = psi * p11 / (psi * p11 + (1 - psi) * p10) print(" %-10.3f" % ppv, end="") print() print() print(" Requiring r positive replicates out of r before calling a") print(" detection, with independent replicates:") print() print(" %-4s %-12s %-12s %-12s" % ("r", "p11^r", "p10^r", "PPV at psi=0.02")) for r in range(1, 6): a = p11 ** r b = 0.02 ** r ppv = 0.02 * a / (0.02 * a + 0.98 * b) print(" %-4d %-12.5f %-12.8f %-12.4f" % (r, a, b, ppv)) print() print(" The cost of that rule is sensitivity: p11^r falls from %.3f to %.3f" % (p11, p11 ** 5)) print(" as r goes from 1 to 5, so the false negative rate climbs from") print(" %.2f to %.2f per sample." % (1 - p11, 1 - p11 ** 5)) print() print(" What a false positive rate does to the pooled comparison. Suppose") print(" eDNA carries p10 = 0.02 and the conventional method carries none.") print(" At psi = 0.05, the apparent eDNA detection rate at a randomly") print(" chosen site is psi*p11 + (1-psi)*p10, and the share of those") print(" positives that are real is only:") for psi in [0.30, 0.10, 0.05, 0.02]: app = psi * p11 + (1 - psi) * 0.02 real = psi * p11 / app print(" psi = %.2f -> apparent positive rate %.4f, of which %.1f%% real" % (psi, app, 100 * real)) # ---- what the article quotes --------------------------------------- head("PART 7. NUMBERS QUOTED IN THE ARTICLE") print() print(" pooled log OR %+.3f" % res["est"]) print(" pooled odds ratio %.2f [%.2f, %.2f]" % (math.exp(res["est"]), math.exp(res["lo"]), math.exp(res["hi"]))) print(" p value %.4f" % res["p"]) print(" Q %.2f on %d df, p = %.5f" % (res["Q"], res["df"], res["Qp"])) print(" I-squared %.1f%%" % res["I2"]) print(" tau %.3f log odds" % res["tau"]) print(" 95%% prediction interval %.2f to %.2f as odds ratios" % (math.exp(pi_lo), math.exp(pi_hi))) print(" Egger intercept %+.3f, p = %.3f" % (eg["intercept"], eg["p"])) print(" leave-one-out OR range %.2f to %.2f" % (math.exp(lo_est), math.exp(hi_est))) print(" studies pooled %d" % len(rows)) print(" studies found and excluded %d" % len(EXCLUDED)) print(" PPV, psi=0.02, p10=0.01, r=1 %.3f" % (0.02 * 0.6 / (0.02 * 0.6 + 0.98 * 0.01))) print(" PPV, psi=0.02, p10=0.02, r=2 %.3f" % (0.02 * 0.36 / (0.02 * 0.36 + 0.98 * 0.0004))) print() rule("=") print("end of output") rule("=") if __name__ == "__main__": main()