""" soil-carbon-review.py -- the Science Journaling Club's pooling code for "How Much Carbon Can Farmland Actually Hold? A Review of the Field Trial Evidence" (Review, Meta-analysis, Soil science). WHAT THIS IS. A random-effects meta-analysis, written from the formulas, of published estimates of soil organic carbon (SOC) change under practices that are sold as carbon sequestration: no-till, reduced tillage, cover cropping, crop residue / straw retention, crop rotation. Three carbon-ADDING practices (manure, manure plus fertiliser, organic substitution for mineral fertiliser) are carried alongside as a contrast, because they are the control group for the whole argument: if a practice that genuinely imports carbon behaves differently with depth from a practice that only rearranges it, that difference is the evidence. WHAT THIS IS NOT. A systematic review. We are a school club. We searched, we read abstracts and open-access results sections, and one person extracted each number while another checked it against the source text. There was no duplicate independent screening, no protocol registration, no grey literature, no contact with authors for unpublished data. Every number below is traceable to a sentence in a published paper, and PART 1 of the printed output lists where each one came from so a reader can go and check it. That is the standard we can actually meet, and it is a lower standard than the word "meta-analysis" usually implies. ONE FURTHER WARNING, STATED UP FRONT. Most of our effect sizes come from published syntheses rather than from individual field trials, because syntheses report an effect size with an interval and individual trials usually do not. Syntheses of the same literature share primary studies. Our entries are therefore NOT independent, which is the assumption every pooling formula below rests on. We report the pooled number because the brief asked for one. We do not think it should be used to price anything. ============================================================================= THE ESTIMATOR ============================================================================= Effects are log response ratios. For a reported percentage change p in SOC under the practice relative to its control, y = ln(1 + p/100) and for a reported 95% interval [lo, hi] on the same percentage scale, se(y) = [ ln(1 + hi/100) - ln(1 + lo/100) ] / (2 * 1.959964) The log scale is used because response ratios are multiplicative and their sampling distribution is much closer to normal in logs. Back-transform with 100*(exp(y) - 1). Fixed-effect weights are inverse variance, w_i = 1/v_i. The DerSimonian-Laird estimate of the between-study variance is tau^2 = max(0, (Q - (k-1)) / C), C = sum(w) - sum(w^2)/sum(w) with Cochran's Q = sum w_i (y_i - theta_FE)^2, and I^2 = max(0, (Q - (k-1)) / Q) Random-effects weights are w*_i = 1/(v_i + tau^2). The prediction interval uses t on k-2 degrees of freedom and the total variance tau^2 + var(theta). Small-study bias is assessed with Egger's regression of the standard normal deviate y_i/se_i on precision 1/se_i; a non-zero intercept indicates funnel asymmetry. ============================================================================= PART 1 VALIDATES THE MACHINERY BEFORE IT TOUCHES REAL DATA ============================================================================= Two checks, both printed: (a) Simulate k studies from a known true effect with known tau^2, pool them, confirm the truth lands inside the stated interval, and repeat the whole thing 2000 times to confirm coverage is near 95%. (b) Build a dataset with no between-study heterogeneity at all, confirm DerSimonian-Laird returns tau^2 = 0, and confirm that the fixed-effect and random-effects estimates then coincide to machine precision. This is an algebraic identity, not a coincidence: at tau^2 = 0 the two weight vectors are the same vector. """ import math import warnings import numpy as np try: with warnings.catch_warnings(): warnings.simplefilter("ignore") from scipy import stats as _st def t_crit(df, p=0.975): return float(_st.t.ppf(p, df)) def t_sf(x, df): return float(_st.t.sf(abs(x), df)) * 2.0 def z_sf(x): return float(_st.norm.sf(abs(x))) * 2.0 def chi2_sf(x, df): return float(_st.chi2.sf(x, df)) except Exception: # pragma: no cover def t_crit(df, p=0.975): return 1.959964 + 2.5 / max(df, 1) def t_sf(x, df): return math.erfc(abs(x) / math.sqrt(2.0)) def z_sf(x): return math.erfc(abs(x) / math.sqrt(2.0)) def chi2_sf(x, df): return float("nan") Z95 = 1.959964 SEP = "=" * 74 def head(t): print("\n" + SEP) print(t) print(SEP) def sub(t): print("\n" + t) print("-" * len(t)) # ========================================================================= # THE DATA # ========================================================================= # Fields: # key, short citation, practice, percent change, 95% low, 95% high, # reported sample size (as the source states it), maximum sampling depth # in cm (None where the source does not state one), family: # "move" practice that mainly redistributes existing carbon # "add" practice that imports carbon from outside the field # provenance: where in the source the number was read. STUDIES = [ dict(key="Du17-ESM", cite="Du et al. 2017", ref=1, practice="no-till (equivalent soil mass)", pct=3.8, lo=1.4, hi=6.3, n="95 comparisons, 57 sites", depth=40, family="move", prov="Abstract, verbatim: '3.8% [95% CI: 1.4%-6.3%; P=0.005]' (ESM)"), dict(key="Du17-FD", cite="Du et al. 2017", ref=1, practice="no-till (fixed depth)", pct=5.1, lo=2.5, hi=7.7, n="95 comparisons, 57 sites", depth=30, family="move", prov="Abstract, verbatim: '5.1% (95%CI: 2.5%-7.7%; P<0.001)' (FD)"), dict(key="Bei23-NT", cite="Beillouin et al. 2023", ref=2, practice="no-till", pct=9.3, lo=5.6, hi=13.0, n="second-order: 230 meta-analyses", depth=None, family="move", prov="Results text: 'no-till farming (+9,3%, CI [+5,6, +13])'"), dict(key="Bei23-RT", cite="Beillouin et al. 2023", ref=2, practice="reduced tillage", pct=12.0, lo=0.1, hi=24.0, n="second-order", depth=None, family="move", prov="Results text: 'reduced tillage intensity (+12%, CI [+0.1, +24])'"), dict(key="Bei23-RES", cite="Beillouin et al. 2023", ref=2, practice="crop residue retention", pct=13.0, lo=9.8, hi=16.0, n="second-order", depth=None, family="move", prov="Results text: 'crop residue retention (+13%, CI [+9.8, +16])'"), dict(key="Bei23-ROT", cite="Beillouin et al. 2023", ref=2, practice="crop rotation", pct=6.5, lo=-0.9, hi=14.0, n="second-order", depth=None, family="move", prov="Results text: 'crop rotation resulted ... 6.5% (CI [-0.9, +14])'"), dict(key="Jian20-CC", cite="Jian et al. 2020", ref=3, practice="cover cropping", pct=15.5, lo=13.8, hi=17.3, n="not stated in abstract", depth=30, family="move", prov="Abstract, verbatim: 'mean change of 15.5% (95% CI 13.8%-17.3%)'"), dict(key="Liu14-STR", cite="Liu et al. 2014", ref=4, practice="straw return", pct=12.8, lo=12.0, hi=13.6, n="176 field studies", depth=20, family="move", prov="Abstract: '12.8 +/- 0.4%'; we read +/- as one standard error " "and widened it to 95% (12.8 +/- 1.96*0.4). FLAGGED: derived."), dict(key="Han16-CFS", cite="Han et al. 2016", ref=5, practice="straw + chemical fertiliser", pct=19.5, lo=18.5, hi=21.5, n="global compilation", depth=20, family="move", prov="Abstract: '2.0 (1.9-2.2) g/kg (19.5%)'. We scaled the reported " "absolute interval by 19.5/2.0. FLAGGED: derived."), dict(key="Han16-CFM", cite="Han et al. 2016", ref=5, practice="manure + chemical fertiliser", pct=36.2, lo=33.1, hi=39.3, n="global compilation", depth=20, family="add", prov="Abstract: '3.5 (3.2-3.8) g/kg (36.2%)'. Scaled by 36.2/3.5. " "FLAGGED: derived."), dict(key="GG21-MAN", cite="Gross & Glaser 2021", ref=6, practice="manure application", pct=35.0, lo=32.0, hi=39.0, n="592 pairwise comparisons, 101 studies", depth=30, family="add", prov="Results text, verbatim: 'increase of SOC stocks of 35% " "(95% CI 32-39%)'"), dict(key="Bei23-ORG", cite="Beillouin et al. 2023", ref=2, practice="organic for mineral substitution", pct=34.0, lo=20.0, hi=49.0, n="second-order", depth=None, family="add", prov="Results text: 'increase in SOC of +34%, CI [+20, +49]'"), ] # Effect sizes we found, read, and could NOT pool, listed so the reader can # see the size of the discard pile rather than only the survivors. NOT_POOLED = [ ("Luo et al. 2010", 7, "no-till, 0-10 cm", "+3.15 +/- 2.42 Mg C/ha", "absolute stock, no control stock given -> no response ratio"), ("Luo et al. 2010", 7, "no-till, 20-40 cm", "-3.30 +/- 1.61 Mg C/ha", "absolute stock, same reason"), ("Luo et al. 2010", 7, "no-till, 0-40 cm", "no significant increase", "reported as a null, no point estimate given in the abstract"), ("Haddaway et al. 2017", 8, "no-till, 0-30 cm, >=10 yr", "+4.6 Mg/ha (0.78-8.43)", "absolute stock, no control stock given"), ("Haddaway et al. 2017", 8, "no-till, full profile", "no effect detected", "reported as a null"), ("Angers & Eriksen-Hamel 2008", 9, "no-till vs full-inversion, profile", "+4.9 Mg/ha", "no interval retrievable from the sources we could read"), ("Virto et al. 2012", 10, "no-till, 0-30 cm ESM", "+6.7% (3.4 Mg C/ha)", "no interval in the abstract"), ("McClelland et al. 2021", 11, "cover crops, 0-30 cm", "+12% (1.11 Mg C/ha)", "interval is in a figure we could not read numerically"), ("Prairie et al. 2023", 12, "no-till, 0-20 cm", "+11.3%", "interval is in Fig. 2 only"), ("Prairie et al. 2023", 12, "no-till, >20 cm", "not significant", "reported as a null"), ("Bai et al. 2019", 13, "conservation tillage", "+5%", "no interval in the abstract"), ("Bai et al. 2019", 13, "cover crops", "+6%", "no interval in the abstract"), ("Crystal-Ornelas et al. 2021", 14, "conservation tillage, organic systems", "+14%", "no interval in the abstract"), ("Emde et al. 2021", 15, "irrigation", "+5.9%", "intervals in Fig. 3 only"), ("Sun et al. 2020", 16, "no-till, 260 paired studies", "no global summary figure", "regional results only"), ] # Rate-scale entries (Mg C per hectare per year). Kept separate because a # rate and a percentage are not the same quantity and must not share a pool. RATES = [ dict(key="WP02", cite="West & Post 2002", ref=17, practice="no-till", est=0.57, se=0.14, n="276 paired treatments, 67 experiments", depth="mostly 0-30 cm", prov="'57 +/- 14 g C m-2 yr-1', converted: 1 g/m2 = 0.01 Mg/ha"), dict(key="PD15", cite="Poeplau & Don 2015", ref=18, practice="cover crops", est=0.32, se=0.08, n="139 plots, 37 sites", depth="22 cm mean", prov="Abstract: '0.32 +/- 0.08 Mg ha-1 yr-1 in a mean soil depth " "of 22 cm'"), dict(key="Du17-FDr", cite="Du et al. 2017 (fixed depth)", ref=1, practice="no-till", est=0.300, se=(0.547 - 0.054) / (2 * Z95), n="95 comparisons", depth="0-30 cm, fixed depth", prov="Abstract: '0.300 Mg ha-1 yr-1, 95%CI: 0.054-0.547'"), dict(key="Du17-ESMr", cite="Du et al. 2017 (equivalent soil mass)", ref=1, practice="no-till", est=0.141, se=(0.384 - (-0.102)) / (2 * Z95), n="95 comparisons", depth="0-30 cm, equivalent soil mass", prov="Abstract: '0.141 Mg ha-1 yr-1, 95%CI: -0.102-0.384'"), ] # Within-study depth contrasts. Not pooled. These are paired comparisons # inside a single synthesis, so everything except depth is held roughly # constant, which makes them stronger evidence on depth than any regression # across studies could be. DEPTH_PAIRS = [ ("Luo et al. 2010", 7, "no-till", "0-10 cm", "+3.15 Mg C/ha", "20-40 cm", "-3.30 Mg C/ha", "0-40 cm: no increase"), ("Du et al. 2017", 1, "no-till", "0-20 cm", "accumulation", "30-40 cm", "depletion", "FD 5.1% vs ESM 3.8%"), ("Prairie et al. 2023", 12, "no-till", "0-20 cm", "+11.3%", ">20 cm", "not significant", "118 studies, 157 experiments"), ("Jian et al. 2020", 3, "cover crops", "<=30 cm", "significant increase", ">30 cm", "not significant", "-"), ("Haddaway et al. 2017", 8, "no-till", "0-30 cm", "+4.6 Mg/ha", "full profile", "no effect", ">=10 yr comparisons"), ("Gross & Glaser 2021", 6, "manure", "<=15 cm", "+40% (conv. tillage)", ">30 cm", "+23% (conv. tillage)", "CARBON-ADDING CONTROL"), ] # ========================================================================= # CORE MACHINERY # ========================================================================= def pct_to_lnrr(pct, lo, hi): """Percentage change and its 95% interval -> log response ratio and se.""" y = math.log1p(pct / 100.0) se = (math.log1p(hi / 100.0) - math.log1p(lo / 100.0)) / (2 * Z95) return y, se def back(y): """Log response ratio -> percentage change.""" return 100.0 * (math.exp(y) - 1.0) def dersimonian_laird(y, se): """Random-effects pooling. Returns a dictionary of everything printed.""" y = np.asarray(y, dtype=float) se = np.asarray(se, dtype=float) v = se ** 2 k = len(y) w = 1.0 / v theta_fe = float(np.sum(w * y) / np.sum(w)) var_fe = float(1.0 / np.sum(w)) Q = float(np.sum(w * (y - theta_fe) ** 2)) df = k - 1 C = float(np.sum(w) - np.sum(w ** 2) / np.sum(w)) tau2 = max(0.0, (Q - df) / C) if C > 0 else 0.0 I2 = max(0.0, (Q - df) / Q) if Q > 0 else 0.0 H2 = Q / df if df > 0 else float("nan") ws = 1.0 / (v + tau2) theta_re = float(np.sum(ws * y) / np.sum(ws)) var_re = float(1.0 / np.sum(ws)) se_re = math.sqrt(var_re) ci = (theta_re - Z95 * se_re, theta_re + Z95 * se_re) if k > 2: tc = t_crit(k - 2) pi_half = tc * math.sqrt(tau2 + var_re) pi = (theta_re - pi_half, theta_re + pi_half) else: pi = (float("nan"), float("nan")) z = theta_re / se_re return dict(k=k, theta_fe=theta_fe, se_fe=math.sqrt(var_fe), theta_re=theta_re, se_re=se_re, ci=ci, pi=pi, Q=Q, df=df, p_Q=chi2_sf(Q, df), tau2=tau2, tau=math.sqrt(tau2), I2=I2, H2=H2, z=z, p=z_sf(z), weights=ws / np.sum(ws)) def egger(y, se): """Egger's regression test for funnel asymmetry. Regress the standard normal deviate y/se on precision 1/se by ordinary least squares. Under symmetry the intercept is zero; a large intercept means small (imprecise) studies sit systematically off to one side. """ y = np.asarray(y, float) se = np.asarray(se, float) snd = y / se prec = 1.0 / se k = len(y) X = np.column_stack([np.ones(k), prec]) beta, *_ = np.linalg.lstsq(X, snd, rcond=None) resid = snd - X @ beta dof = k - 2 s2 = float(resid @ resid) / dof cov = s2 * np.linalg.inv(X.T @ X) se_int = math.sqrt(cov[0, 0]) t = beta[0] / se_int return dict(intercept=float(beta[0]), se=se_int, slope=float(beta[1]), t=float(t), df=dof, p=t_sf(t, dof)) def metareg(y, se, x): """Weighted meta-regression with a method-of-moments residual tau^2. Fit by weighted least squares with w = 1/(v + tau2_res), where tau2_res is found by matching the weighted residual sum of squares to its expectation under the model. Three iterations is plenty here. """ y = np.asarray(y, float) v = np.asarray(se, float) ** 2 x = np.asarray(x, float) k = len(y) X = np.column_stack([np.ones(k), x]) p = X.shape[1] tau2 = 0.0 for _ in range(40): w = 1.0 / (v + tau2) W = np.diag(w) XtWX = X.T @ W @ X beta = np.linalg.solve(XtWX, X.T @ W @ y) resid = y - X @ beta Qres = float(resid @ W @ resid) # trace term for the method-of-moments update P = W - W @ X @ np.linalg.solve(XtWX, X.T @ W) tr = float(np.trace(P @ np.diag(v))) denom = float(np.trace(P)) new = max(0.0, tau2 + (Qres - (k - p)) / denom) if denom > 0 else 0.0 # the direct moment solution, guarded new2 = max(0.0, (Qres - tr) / denom) if denom > 0 else 0.0 cand = new2 if abs(new2 - tau2) < abs(new - tau2) else new if abs(cand - tau2) < 1e-12: tau2 = cand break tau2 = cand w = 1.0 / (v + tau2) W = np.diag(w) XtWX = X.T @ W @ X cov = np.linalg.inv(XtWX) beta = cov @ (X.T @ W @ y) resid = y - X @ beta Qres = float(resid @ W @ resid) se_b = np.sqrt(np.diag(cov)) z = beta / se_b return dict(beta=beta, se=se_b, z=z, p=[z_sf(t) for t in z], tau2=tau2, Qres=Qres, df=k - p, k=k, p_Qres=chi2_sf(Qres, k - p)) # ========================================================================= # PART 0 -- HEADER # ========================================================================= print(SEP) print("HOW MUCH CARBON CAN FARMLAND ACTUALLY HOLD?") print("Random-effects pooling of published soil organic carbon effect sizes") print("Science Journaling Club -- our own code, written from the formulas") print(SEP) print(""" Read PART 1 before PART 3. The machinery is checked against a known answer first, because a pooling routine that has never been asked to recover a number you already know is a routine you have no reason to trust. """.strip()) # ========================================================================= # PART 1 -- VALIDATION # ========================================================================= head("PART 1 -- DOES THE MACHINERY WORK? (two checks, before any real data)") sub("Check A: recover a known true effect from synthetic data") rng = np.random.default_rng(20250913) TRUE_MU = math.log(1.12) # a true 12% increase TRUE_TAU2 = 0.0025 # real between-study spread K_SIM = 12 theta_i = rng.normal(TRUE_MU, math.sqrt(TRUE_TAU2), K_SIM) se_i = rng.uniform(0.012, 0.075, K_SIM) y_i = rng.normal(theta_i, se_i) r = dersimonian_laird(y_i, se_i) print(f" true mu = {TRUE_MU:.6f} lnRR (= {back(TRUE_MU):+.2f}% change)") print(f" true tau^2 = {TRUE_TAU2:.6f}") print(f" k = {K_SIM} synthetic studies") print(f" pooled (random) = {r['theta_re']:.6f} " f"[{r['ci'][0]:.6f}, {r['ci'][1]:.6f}]") print(f" ... as a percentage = {back(r['theta_re']):+.2f}% " f"[{back(r['ci'][0]):+.2f}%, {back(r['ci'][1]):+.2f}%]") print(f" estimated tau^2 = {r['tau2']:.6f} (true {TRUE_TAU2:.6f})") print(f" I^2 = {100*r['I2']:.1f}%") inside = r['ci'][0] <= TRUE_MU <= r['ci'][1] print(f" TRUTH INSIDE THE 95% INTERVAL: {'YES' if inside else 'NO'}") # coverage over many replicates -- one run landing inside proves nothing cover = 0 tau2s = [] NREP = 2000 for _ in range(NREP): th = rng.normal(TRUE_MU, math.sqrt(TRUE_TAU2), K_SIM) s = rng.uniform(0.012, 0.075, K_SIM) yy = rng.normal(th, s) rr = dersimonian_laird(yy, s) tau2s.append(rr["tau2"]) if rr["ci"][0] <= TRUE_MU <= rr["ci"][1]: cover += 1 print(f" coverage over {NREP} replicates = {100*cover/NREP:.1f}% " f"(nominal 95%)") print(f" mean estimated tau^2 over replicates = {np.mean(tau2s):.6f} " f"(true {TRUE_TAU2:.6f})") print(" Two things to read off that pair of lines. The estimator recovers") print(" tau^2 almost exactly on average, so it is not biased. Coverage") print(" still lands below the nominal 95%, because the interval treats the") print(" estimated tau^2 as if it were known exactly, and it is not. Our") print(" intervals on the real data are therefore a little too narrow, and") print(" we would rather print that than hide it.") sub("Check B: at tau^2 = 0 the fixed and random estimates must coincide") # Build data with literally no heterogeneity: all studies share one true # effect, and we hand the routine the exact fitted values so Q <= k-1. mu0 = math.log(1.07) se0 = np.array([0.02, 0.03, 0.05, 0.04, 0.025, 0.06]) y0 = np.full(len(se0), mu0) # identical effects -> Q = 0 exactly r0 = dersimonian_laird(y0, se0) print(f" Q = {r0['Q']:.12f} (df = {r0['df']})") print(f" tau^2 (DL) = {r0['tau2']:.12f}") print(f" I^2 = {100*r0['I2']:.4f}%") print(f" fixed-effect theta = {r0['theta_fe']:.15f}") print(f" random-effects theta = {r0['theta_re']:.15f}") gap = abs(r0["theta_fe"] - r0["theta_re"]) print(f" |difference| = {gap:.3e}") print(f" IDENTITY HOLDS: {'YES' if gap < 1e-12 else 'NO'}") print(" Why it must: at tau^2 = 0 the random-effects weight 1/(v+tau^2)") print(" IS the fixed-effect weight 1/v, so the two weighted means are the") print(" same weighted mean. Any code that fails this is wired wrong.") # a second, less trivial tau^2 = 0 case: real spread, but less than sampling # error alone would predict, so DL truncates at zero se1 = np.array([0.03, 0.05, 0.04, 0.06, 0.035]) y1 = np.array([0.0700, 0.0690, 0.0705, 0.0695, 0.0702]) r1 = dersimonian_laird(y1, se1) print(f"\n second case (real but tiny spread): Q = {r1['Q']:.6f}, " f"df = {r1['df']}, tau^2 = {r1['tau2']:.12f}") print(f" fixed {r1['theta_fe']:.15f} vs random {r1['theta_re']:.15f} " f"-> gap {abs(r1['theta_fe']-r1['theta_re']):.3e}") # ========================================================================= # PART 2 -- THE DATASET # ========================================================================= head("PART 2 -- THE DATASET, WITH PROVENANCE FOR EVERY NUMBER") for s in STUDIES: s["y"], s["se"] = pct_to_lnrr(s["pct"], s["lo"], s["hi"]) print(f"{'key':<11s}{'source':<25s}{'practice':<34s}" f"{'pct':>7s}{'lo':>7s}{'hi':>7s}{'lnRR':>8s}{'se':>8s}{'depth':>7s}") print("-" * 114) for s in STUDIES: d = f"{s['depth']}" if s["depth"] else "n/s" print(f"{s['key']:<11s}{s['cite']:<25s}{s['practice']:<34s}" f"{s['pct']:>7.1f}{s['lo']:>7.1f}{s['hi']:>7.1f}" f"{s['y']:>8.4f}{s['se']:>8.4f}{d:>7s}") sub("Where each number came from") for s in STUDIES: print(f" [{s['ref']}] {s['key']:<11s} {s['prov']}") sub("Effect sizes we found but could NOT pool") print(f" {len(NOT_POOLED)} entries discarded. The discard pile matters: it is") print(" larger than the pool, and it is not a random sample of the evidence.") print(" Nulls and absolute-stock results are over-represented in it, which") print(" biases whatever survives UPWARD.\n") for c, ref, what, val, why in NOT_POOLED: print(f" [{ref}] {c:<30s} {what:<34s} {val:<26s} {why}") # ========================================================================= # PART 3 -- POOLING # ========================================================================= head("PART 3 -- RANDOM-EFFECTS POOL (all 12 effect sizes)") Y = [s["y"] for s in STUDIES] SE = [s["se"] for s in STUDIES] res = dersimonian_laird(Y, SE) print(f" k = {res['k']}") print(f" fixed-effect pooled = {res['theta_fe']:.5f} lnRR " f"= {back(res['theta_fe']):+.2f}%") print(f" RANDOM-EFFECTS POOL = {res['theta_re']:.5f} lnRR " f"= {back(res['theta_re']):+.2f}%") print(f" 95% CI = [{res['ci'][0]:.5f}, {res['ci'][1]:.5f}] " f"= [{back(res['ci'][0]):+.2f}%, {back(res['ci'][1]):+.2f}%]") print(f" 95% PREDICTION INT. = [{back(res['pi'][0]):+.2f}%, " f"{back(res['pi'][1]):+.2f}%]") print(f" z = {res['z']:.2f}, p = {res['p']:.3g}") print() print(f" Cochran's Q = {res['Q']:.2f} on {res['df']} df, " f"p = {res['p_Q']:.3g}") print(f" tau^2 = {res['tau2']:.5f} (tau = {res['tau']:.4f} " f"on the log scale)") print(f" I^2 = {100*res['I2']:.1f}%") print(f" H^2 = {res['H2']:.1f}") print() print(" READ THE PREDICTION INTERVAL, NOT THE CONFIDENCE INTERVAL. The") print(" confidence interval describes how well we know the AVERAGE of these") print(" studies. The prediction interval describes the range a NEW field") print(" would plausibly fall in, and it is the one a farmer or a credit") print(" buyer actually faces.") sub("Study weights in the random-effects pool") for s, w in zip(STUDIES, res["weights"]): print(f" {s['key']:<11s} {100*w:>6.2f}% {s['cite']} / {s['practice']}") print("\n Note how flat these are. With tau^2 this large, tau^2 swamps the") print(" sampling variances and every study ends up weighted almost equally,") print(" so a huge synthesis of 592 comparisons counts about the same as one") print(" line lifted from a second-order review.") # ========================================================================= # PART 4 -- SUBGROUPS # ========================================================================= head("PART 4 -- SPLIT BY WHAT THE PRACTICE ACTUALLY DOES") for fam, label in [("move", "REARRANGES carbon already in the field " "(tillage, residues, cover crops, rotation)"), ("add", "IMPORTS carbon from outside the field " "(manure, organic substitution)")]: sel = [s for s in STUDIES if s["family"] == fam] rr = dersimonian_laird([s["y"] for s in sel], [s["se"] for s in sel]) print(f"\n {label}") print(f" k = {rr['k']}, pooled = {back(rr['theta_re']):+.2f}% " f"[{back(rr['ci'][0]):+.2f}%, {back(rr['ci'][1]):+.2f}%], " f"I^2 = {100*rr['I2']:.1f}%, tau^2 = {rr['tau2']:.4f}") for s in sel: print(f" {s['key']:<11s} {s['pct']:>6.1f}% {s['practice']}") move = [s for s in STUDIES if s["family"] == "move"] add = [s for s in STUDIES if s["family"] == "add"] rm = dersimonian_laird([s["y"] for s in move], [s["se"] for s in move]) ra = dersimonian_laird([s["y"] for s in add], [s["se"] for s in add]) diff = ra["theta_re"] - rm["theta_re"] sed = math.sqrt(ra["se_re"] ** 2 + rm["se_re"] ** 2) print(f"\n Difference (add - move) = {diff:.4f} lnRR, se {sed:.4f}, " f"z = {diff/sed:.2f}, p = {z_sf(diff/sed):.3g}") print(f" On the percentage scale the adding practices run about " f"{back(ra['theta_re'])/back(rm['theta_re']):.1f}x the rearranging ones.") # ========================================================================= # PART 5 -- THE DEPTH QUESTION # ========================================================================= head("PART 5 -- SAMPLING DEPTH") withd = [s for s in STUDIES if s["depth"] is not None] print(f" Of {len(STUDIES)} pooled effect sizes, {len(withd)} state a maximum") print(" sampling depth. That is the first result of this section and it is") print(" not a good one.\n") print(f" {'key':<11s}{'depth cm':>9s}{'pct':>8s}{'se(lnRR)':>10s} source") for s in sorted(withd, key=lambda a: a["depth"]): print(f" {s['key']:<11s}{s['depth']:>9d}{s['pct']:>8.1f}" f"{s['se']:>10.4f} {s['cite']}") mr = metareg([s["y"] for s in withd], [s["se"] for s in withd], [s["depth"] for s in withd]) b0, b1 = mr["beta"] print(f"\n Meta-regression of lnRR on maximum sampling depth (k = {mr['k']}):") print(f" intercept = {b0:+.5f} (se {mr['se'][0]:.5f})") print(f" slope = {b1:+.6f} per cm (se {mr['se'][1]:.6f}), " f"z = {mr['z'][1]:.2f}, p = {mr['p'][1]:.3g}") print(f" residual tau^2 = {mr['tau2']:.5f}, Q_res = {mr['Qres']:.2f} " f"on {mr['df']} df") per10 = 100 * (math.exp(b1 * 10) - 1) print(f" A slope of {b1:+.6f} per cm means each extra 10 cm of required") print(f" sampling depth moves the estimated effect by {per10:+.2f}%.") print(" With k = 7 and one point beyond 30 cm, this regression is") print(" underpowered and we are not going to pretend otherwise. The") print(" number is printed so a reader can see how weak it is.") sub("Pooled effect under a minimum-depth requirement") print(" This is the calculation that matters for a credit protocol: if you") print(" refuse to count any study that did not sample at least D cm, what") print(" is left, and what does it say?\n") print(f" {'require >= D cm':<18s}{'k':>4s}{'pooled %':>11s}" f"{'95% CI':>24s}{'I^2':>8s}") DEPTH_GATES = [0, 20, 30, 40] depth_curve = [] for D in DEPTH_GATES: sel = [s for s in STUDIES if (s["depth"] or 0) >= D] if len(sel) == 0: print(f" {'>= ' + str(D):<18s}{0:>4d} nothing left") continue if len(sel) == 1: s = sel[0] lo, hi = s["y"] - Z95 * s["se"], s["y"] + Z95 * s["se"] print(f" {'>= ' + str(D):<18s}{1:>4d}{back(s['y']):>11.2f}" f"{'[' + f'{back(lo):+.2f}, {back(hi):+.2f}' + ']':>24s}" f"{'n/a':>8s}") depth_curve.append((D, 1, back(s["y"]), back(lo), back(hi))) continue rr = dersimonian_laird([s["y"] for s in sel], [s["se"] for s in sel]) ci_txt = f"[{back(rr['ci'][0]):+.2f}, {back(rr['ci'][1]):+.2f}]" print(f" {'>= ' + str(D):<18s}{rr['k']:>4d}{back(rr['theta_re']):>11.2f}" f"{ci_txt:>24s}{100*rr['I2']:>7.1f}%") depth_curve.append((D, rr["k"], back(rr["theta_re"]), back(rr["ci"][0]), back(rr["ci"][1]))) sub("Within-study depth contrasts (NOT pooled, and stronger for that)") print(" Each row is one synthesis comparing its own shallow layer with its") print(" own deep layer. Everything except depth is held constant, which no") print(" between-study regression can manage.\n") for cite, ref, prac, d1, v1, d2, v2, note in DEPTH_PAIRS: print(f" [{ref}] {cite:<26s} {prac:<12s} {d1:>10s}: {v1:<24s}" f" {d2:>12s}: {v2:<22s} {note}") print("\n Five of the six rows are practices that rearrange carbon, and all") print(" five lose the effect below the plough layer. The sixth is manure,") print(" which is carbon carted in from somewhere else, and it holds a") print(" significant effect below 30 cm. That contrast is the argument.") # ========================================================================= # PART 6 -- SMALL-STUDY BIAS # ========================================================================= head("PART 6 -- FUNNEL ASYMMETRY AND EGGER'S TEST") eg = egger(Y, SE) print(f" Egger intercept = {eg['intercept']:+.3f} (se {eg['se']:.3f})") print(f" t = {eg['t']:.2f} on {eg['df']} df, p = {eg['p']:.3g}") print(f" slope (the bias-adjusted effect) = {eg['slope']:+.4f} lnRR " f"= {back(eg['slope']):+.2f}%") print() print(" Funnel coordinates (for Figure 4):") print(f" {'key':<11s}{'lnRR':>9s}{'se':>9s}{'1/se':>9s}{'y/se':>9s}") for s in STUDIES: print(f" {s['key']:<11s}{s['y']:>9.4f}{s['se']:>9.4f}" f"{1/s['se']:>9.2f}{s['y']/s['se']:>9.2f}") print() print(" CAUTION. Egger's test is close to meaningless at k = 12 with I^2") print(" this high: heterogeneity alone generates funnel asymmetry, and with") print(" a dozen points the test has almost no power. Treat the number as a") print(" description of the plot, not as a hypothesis test.") # ========================================================================= # PART 7 -- LEAVE ONE OUT # ========================================================================= head("PART 7 -- LEAVE-ONE-OUT SENSITIVITY") print(f" {'omitted':<11s}{'k':>3s}{'pooled %':>11s}{'95% CI':>24s}" f"{'I^2':>8s}{'tau^2':>9s}") loo = [] for i, s in enumerate(STUDIES): sel = [t for j, t in enumerate(STUDIES) if j != i] rr = dersimonian_laird([t["y"] for t in sel], [t["se"] for t in sel]) ci_s = f"[{back(rr['ci'][0]):+.2f}, {back(rr['ci'][1]):+.2f}]" print(f" {s['key']:<11s}{rr['k']:>3d}{back(rr['theta_re']):>11.2f}" f"{ci_s:>24s}{100*rr['I2']:>7.1f}%{rr['tau2']:>9.4f}") loo.append(back(rr["theta_re"])) print(f"\n Range across leave-one-out fits: {min(loo):.2f}% to {max(loo):.2f}%") print(f" Full-data pool: {back(res['theta_re']):.2f}%") print(f" Widest single-study swing: {max(abs(x - back(res['theta_re'])) for x in loo):.2f} " f"percentage points") sub("Pre-planned sensitivity: one effect size per publication") print(" Beillouin et al. contribute four rows and Han et al. two, which") print(" double-counts those papers. Keeping only the first row from each") print(" publication gives:") seen = set() onepp = [] for s in STUDIES: if s["cite"] in seen: continue seen.add(s["cite"]) onepp.append(s) r_one = dersimonian_laird([s["y"] for s in onepp], [s["se"] for s in onepp]) print(f" k = {r_one['k']}, pooled = {back(r_one['theta_re']):+.2f}% " f"[{back(r_one['ci'][0]):+.2f}%, {back(r_one['ci'][1]):+.2f}%], " f"I^2 = {100*r_one['I2']:.1f}%") print(" kept:", ", ".join(s["key"] for s in onepp)) sub("Pre-planned sensitivity: drop the two derived intervals") noderived = [s for s in STUDIES if s["key"] not in ("Liu14-STR", "Han16-CFS", "Han16-CFM")] r_nd = dersimonian_laird([s["y"] for s in noderived], [s["se"] for s in noderived]) print(f" k = {r_nd['k']}, pooled = {back(r_nd['theta_re']):+.2f}% " f"[{back(r_nd['ci'][0]):+.2f}%, {back(r_nd['ci'][1]):+.2f}%], " f"I^2 = {100*r_nd['I2']:.1f}%") # ========================================================================= # PART 8 -- THE RATE POOL # ========================================================================= head("PART 8 -- THE RATE POOL (Mg C per hectare per year)") print(" A percentage is not a credit. Credits are issued in tonnes per") print(" hectare per year, so here is the small set of estimates reported in") print(" that unit, kept in a separate pool because you cannot average a") print(" percentage with a rate.\n") print(f" {'key':<12s}{'practice':<14s}{'rate':>7s}{'se':>7s}" f"{'95% CI':>20s} depth") for r_ in RATES: lo, hi = r_["est"] - Z95 * r_["se"], r_["est"] + Z95 * r_["se"] print(f" {r_['key']:<12s}{r_['practice']:<14s}{r_['est']:>7.3f}" f"{r_['se']:>7.3f}{'[' + f'{lo:+.3f}, {hi:+.3f}' + ']':>20s} " f"{r_['depth']}") rr = dersimonian_laird([r_["est"] for r_ in RATES], [r_["se"] for r_ in RATES]) print(f"\n Pooled rate = {rr['theta_re']:.3f} Mg C/ha/yr " f"[{rr['ci'][0]:.3f}, {rr['ci'][1]:.3f}]") print(f" Q = {rr['Q']:.2f} on {rr['df']} df, tau^2 = {rr['tau2']:.5f}, " f"I^2 = {100*rr['I2']:.1f}%") print(" (No log transform here: these are already differences, not ratios.)") sub("The accounting pair -- the cleanest number in this whole review") fd = [r_ for r_ in RATES if r_["key"] == "Du17-FDr"][0] esm = [r_ for r_ in RATES if r_["key"] == "Du17-ESMr"][0] print(f" Same 95 comparisons, same 57 sites, same 0-30 cm, same authors.") print(f" fixed depth : {fd['est']:.3f} Mg C/ha/yr " f"[0.054, 0.547] significant") print(f" equivalent soil mass : {esm['est']:.3f} Mg C/ha/yr " f"[-0.102, 0.384] NOT significant") print(f" ratio : {fd['est']/esm['est']:.2f}x") print(f" difference : {fd['est']-esm['est']:.3f} Mg C/ha/yr") print(" Nothing changed in the field. Only the arithmetic changed. No-till") print(" loosens soil, bulk density falls, and a fixed 30 cm core therefore") print(" scoops up less soil under no-till than under the plough. Correct") print(" for that by comparing equal masses of soil rather than equal depths") print(" and more than half the sequestration rate goes away.") # ========================================================================= # PART 9 -- SATURATION # ========================================================================= head("PART 9 -- SATURATION: HOW LONG DOES A PRACTICE ACTUALLY BUY?") print(""" A soil is not a bucket you can keep filling. Carbon comes in from plant residues and roots and leaves again through microbial respiration, and the loss term grows with the stock itself. The simplest honest model is one line: dC/dt = I - k*C with C the stock in Mg C/ha, I the annual carbon input reaching the stable pool, and k the fractional loss rate per year. It has a closed solution: C(t) = C_eq + (C0 - C_eq) * exp(-k*t), C_eq = I / k Raise the input and the equilibrium rises proportionally, but the approach to it decays exponentially. That exponential is the whole argument about permanence: the practice does not stop working because the farmer stops trying, it stops working because the stock has climbed near its new ceiling and respiration has caught up with the extra input. """.rstrip()) C0 = 45.0 # Mg C/ha in 0-30 cm, a mid-range arable stock K = 0.03 # per year I0 = C0 * K # input that holds the stock where it already is DI = 0.50 # extra input from the practice, Mg C/ha/yr def stock(t, c0=C0, k=K, inp=I0 + DI): ceq = inp / k return ceq + (c0 - ceq) * math.exp(-k * t) CEQ = (I0 + DI) / K print(f"\n C0 = {C0:.1f} Mg C/ha (0-30 cm), k = {K:.3f}/yr") print(f" baseline input I0 = C0*k = {I0:.3f} Mg C/ha/yr (holds C0 steady)") print(f" practice adds dI = {DI:.2f} Mg C/ha/yr") print(f" new equilibrium C_eq = (I0+dI)/k = {CEQ:.2f} Mg C/ha") print(f" total gain available = {CEQ - C0:.2f} Mg C/ha, and not one tonne more") print(f"\n {'year':>6s}{'stock':>10s}{'cumulative gain':>18s}" f"{'that year':>12s}{'% of total':>12s}") for t in [0, 5, 10, 20, 30, 40, 50, 75, 100, 150]: c = stock(t) rate = DI * math.exp(-K * t) print(f" {t:>6d}{c:>10.2f}{c - C0:>18.2f}{rate:>12.3f}" f"{100*(c-C0)/(CEQ-C0):>11.1f}%") for floor in [0.30, 0.20, 0.10, 0.05]: yrs = math.log(DI / floor) / K print(f"\n Years until the annual gain falls below {floor:.2f} " f"Mg C/ha/yr: {yrs:.1f}") half = math.log(2) / K print(f" Half the total gain is banked by year {half:.1f}.") print(f" Ninety per cent by year {math.log(10)/K:.1f}.") print("\n Published durations, for comparison with the model:") print(" Liu et al. 2014 straw return saturates after about 12 years") print(" Han et al. 2016 28-73 yr (straw+NPK), 26-117 yr (manure+NPK)") print(" Poeplau & Don 2015 new steady state after ~155 yr, total gain") print(" 16.7 +/- 1.5 Mg C/ha at 22 cm depth") print(" Poulton et al. 2018 Rothamsted farmyard manure gave 18 and 43 per") print(" mille per year in the first 20 years, still") print(" above 7 per mille at 40-60 years, then fell") print(" Georgiou et al. 2022 accrual runs ~3x faster in soils at a tenth") print(" of their mineral capacity than at a half") sub("The Georgiou capacity check, done arithmetically") print(" Georgiou et al. report global mineral-associated carbon at 42% of") print(" mineralogical capacity in surface soil and 21% at depth, across 1144") print(" profiles, and accrual roughly 3x faster at one tenth of capacity") print(" than at one half. In our one-line model the accrual rate at stock C") print(" is I - k*C, so the ratio of rates at two stocks C_a and C_b is:") for ca, cb in [(0.1, 0.5)]: Cmax = 100.0 Ia = (I0 + DI) ra_ = Ia - K * (ca * Cmax) rb_ = Ia - K * (cb * Cmax) print(f" at {ca:.0%} of a {Cmax:.0f} Mg/ha capacity: " f"{ra_:.3f} Mg C/ha/yr") print(f" at {cb:.0%} of the same capacity: " f"{rb_:.3f} Mg C/ha/yr") print(f" ratio = {ra_/rb_:.2f}x, against their reported ~3x.") print(" The one-line model gets the sign and the rough size of the") print(" capacity effect without being told about it. That is a weak check,") print(" not a validation, and we are calling it weak.") # ========================================================================= # PART 10 -- WHAT A CREDIT WOULD BE WORTH # ========================================================================= head("PART 10 -- THE CREDIT ARITHMETIC, IN PLAIN NUMBERS") CO2_PER_C = 44.0 / 12.0 for label, rate in [("pooled rate, all four entries", rr["theta_re"]), ("Du et al. fixed depth", fd["est"]), ("Du et al. equivalent soil mass", esm["est"]), ("West & Post no-till", 0.57), ("Poeplau & Don cover crops", 0.32)]: co2 = rate * CO2_PER_C print(f" {label:<34s}{rate:>7.3f} Mg C/ha/yr = " f"{co2:>6.3f} t CO2e/ha/yr") for price in [20, 50, 100]: print(f" at ${price:>3d}/t CO2e: ${co2*price:>7.2f} per hectare " f"per year") print("\n A 200-hectare arable farm, at the equivalent-soil-mass rate and") esm_co2 = esm["est"] * CO2_PER_C print(f" $50/t CO2e: {esm_co2:.3f} t CO2e/ha/yr x 200 ha x $50 = " f"${esm_co2*200*50:,.0f} per year, gross.") print(" Soil sampling to 30 cm on a grid, laboratory analysis, verification") print(" and registry fees are not free. We are not going to invent a cost") print(" figure, but the reader can see the size of the pot being divided.") # ========================================================================= # PART 11 -- HEADLINE NUMBERS FOR THE ARTICLE # ========================================================================= head("PART 11 -- HEADLINE NUMBERS QUOTED IN THE ARTICLE") print(f" studies screened and read : 22") print(f" effect sizes extracted : {len(STUDIES) + len(NOT_POOLED)}") print(f" effect sizes with usable uncertainty: {len(STUDIES)}") print(f" distinct publications in the pool : " f"{len(set(s['cite'] for s in STUDIES))}") print(f" pooled SOC change (random effects) : {back(res['theta_re']):+.1f}% " f"[{back(res['ci'][0]):+.1f}%, {back(res['ci'][1]):+.1f}%]") print(f" 95% prediction interval : " f"[{back(res['pi'][0]):+.1f}%, {back(res['pi'][1]):+.1f}%]") print(f" I^2 : {100*res['I2']:.0f}%") print(f" tau^2 : {res['tau2']:.4f}") print(f" Q : {res['Q']:.1f} on " f"{res['df']} df") print(f" rearranging practices pooled : {back(rm['theta_re']):+.1f}% " f"[{back(rm['ci'][0]):+.1f}%, {back(rm['ci'][1]):+.1f}%] (k={rm['k']})") print(f" importing practices pooled : {back(ra['theta_re']):+.1f}% " f"[{back(ra['ci'][0]):+.1f}%, {back(ra['ci'][1]):+.1f}%] (k={ra['k']})") print(f" depth meta-regression slope : {b1:+.6f} lnRR per cm " f"({per10:+.2f}% per 10 cm), p = {mr['p'][1]:.2f}") print(f" effect sizes stating a max depth : {len(withd)} of {len(STUDIES)}") print(f" Egger intercept : {eg['intercept']:+.2f}, " f"p = {eg['p']:.2f}") print(f" leave-one-out range : {min(loo):.1f}% to {max(loo):.1f}%") print(f" one-row-per-paper pool : {back(r_one['theta_re']):+.1f}% " f"(k={r_one['k']})") print(f" fixed depth vs equivalent soil mass : {fd['est']:.3f} vs " f"{esm['est']:.3f} Mg C/ha/yr ({fd['est']/esm['est']:.1f}x)") print(f" saturation ceiling, default settings: {CEQ:.1f} Mg C/ha " f"(gain {CEQ-C0:.1f})") print(f" gain banked by year 20 : {stock(20)-C0:.2f} Mg C/ha") print(f" gain banked by year 50 : {stock(50)-C0:.2f} Mg C/ha") print(f" years until gain < 0.10 Mg C/ha/yr : " f"{math.log(DI/0.10)/K:.1f}") print("\n" + SEP) print("WHAT THIS NUMBER IS WORTH") print(SEP) print(""" Not much, as a number. I^2 near the ceiling means the studies are not estimating one common effect, so the pooled percentage is an average over a set of things that are not the same thing. The prediction interval is the honest summary and it is wide enough to include no change at all. What the review does establish, and we think it establishes it firmly: 1. Practices that rearrange existing carbon lose their effect below the plough layer. Practices that import carbon do not. 2. The measured effect depends on the accounting convention by a factor of two, in the same dataset, with nothing happening in the field. 3. Most of the literature does not sample deep enough to tell the difference, and a protocol that accepts 30 cm cores cannot either. 4. Whatever the rate is, it is temporary, and the model that says so is one line long. A club extracting effect sizes from abstracts is not a systematic review with two independent screeners. Weight this accordingly, and weight it less than the individual long-term experiments it is built from. Those trials ran for decades on budgets nobody wanted to renew, and they are the reason any of this is arguable at all. """.strip()) print(SEP)