#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ allometry-scaling.py Science Journaling Club, Volume 2, Issue 4, Summer 2026. "Does Metabolism Really Scale to the Three Quarter Power" ============================================================================= QUESTION ============================================================================= Kleiber's law states that whole-organism metabolic rate B scales with body mass M as B = a * M^b, with b = 3/4. Rubner's older surface-law argument gives b = 2/3. The two candidate exponents differ by 1/12 = 0.0833, and the argument over which is right has run since the 1930s. We ask three things of one published compilation: (1) When the data are fitted properly, is b = 3/4, b = 2/3, or neither? (2) How much does the answer move when you switch fitting method from ordinary least squares to reduced major axis, a switch that is itself one of the historical grounds of dispute? (3) How much does the answer move when you change which animals are in the sample: different taxonomic classes, and narrower mass windows? ============================================================================= MODEL ============================================================================= The whole study is a regression on logarithms. Write y = log10(B / watt), x = log10(M / gram), y = alpha + b*x + e. Two estimators of b are computed, by hand, from the sums of squares: OLS (ordinary least squares, y on x) b_ols = Sxy / Sxx This is the maximum-likelihood slope when all the error is in y and none in x. It is the standard choice in the metabolic-scaling literature and it is what Kleiber himself used. RMA (reduced major axis, also called standardised major axis, SMA) b_rma = sign(r) * sd(y) / sd(x) = b_ols / |r| This is the slope that minimises the sum of the areas of the triangles between points and line. It is symmetric in x and y, and it is the usual choice when mass is itself measured with error, which it is. The identity b_rma = b_ols / |r| is exact and is checked numerically below. Because |r| <= 1, RMA is ALWAYS at least as steep as OLS. The gap is entirely determined by the scatter: b_rma/b_ols - 1 = 1/|r| - 1. Inference: OLS slope: t = (b - beta0) / SE(b), df = n-2, with SE(b) = sqrt( SSE/(n-2) / Sxx ). RMA slope: confidence interval by the standard formula (Warton et al. 2006), b * (sqrt(B+1) +/- sqrt(B)) with B = F_{1,n-2} (1-r^2)/(n-2); point hypothesis b = beta0 tested by the exact correlation test: rescale y' = y/beta0, then H0 becomes "SMA slope = 1", which holds iff corr(y'-x, y'+x) = 0. Slope heterogeneity across taxa: an ANCOVA-style F test comparing the common-slope model to the separate-slopes model. Curvature: a quadratic term in x (Kolokotrones et al. 2010). The t and F distributions are implemented here from the regularised incomplete beta function (Lentz continued fraction) so that the script needs only numpy. Both are validated against published table values before use. ============================================================================= DATA ============================================================================= Source: AnAge, the Animal Ageing and Longevity Database, Build 15 (release date 3 July 2023), part of the Human Ageing Genomic Resources. URL: https://genomics.senescence.info/species/dataset.zip Retrieved: 14 September 2026, cached at analysis/data/anage_dataset.zip and extracted to analysis/data/anage_data.txt. The SHA-256 of the cached archive is printed at run time. Citation: Tacutu, R. et al. (2018) Human Ageing Genomic Resources: new and updated databases. Nucleic Acids Research 46(D1), D1083-D1090. doi:10.1093/nar/gkx1042 Fields used: "Body mass (g)", "Metabolic rate (W)", "Class", "Order", "Genus", "Species". AnAge's metabolic rate field is a compiled basal or standard metabolic rate in watts, drawn by AnAge's curators from the primary physiological literature; it is a secondary compilation and we did not go behind it to the original papers. Every species carrying BOTH a positive body mass and a positive metabolic rate is used. Nothing is excluded, no outlier is trimmed, and no value was entered by hand. NO ANIMAL WAS MEASURED BY US. The club has no laboratory, no respirometer and no animals. The computation is the experiment. Every number below is the output of arithmetic performed on somebody else's published measurements. ============================================================================= ASSUMPTIONS ============================================================================= A1. AnAge's compiled metabolic rates are comparable across species, i.e. that "basal" or "standard" means roughly the same thing in every source paper AnAge drew on. This is the weakest assumption in the study. A2. Error in log mass is small relative to error in log metabolic rate (required by OLS, not by RMA). Reporting both is the response to not knowing the error ratio. A3. Residuals are independent across species. THIS IS FALSE. Species share ancestry, so close relatives resemble each other and the effective sample size is smaller than n. All confidence intervals below are therefore too narrow, probably substantially. We do not have a phylogeny, so we cannot fix this; we quantify the direction of the problem by refitting at the level of taxonomic order, where the pseudoreplication is much reduced. A4. A single power law is the right functional form. Section 9 tests this by adding a quadratic term, and it fails. ============================================================================= LIMITATIONS, STATED PLAINLY ============================================================================= L1. No phylogenetic correction (A3). Published phylogenetic analyses of mammal BMR (Capellini et al. 2010; Sieg et al. 2009; White & Seymour 2003) shift the exponent and widen the intervals. L2. No body-temperature correction. AnAge carries a temperature field but it is sparse and inconsistently defined, so we left it alone. Temperature is known to matter (Clarke et al. 2010; Gillooly et al. 2001). L3. Ectotherms (16 reptiles, 18 amphibians here) are measured at ambient temperature and are not thermally comparable to endotherms. We report them, and we say so wherever a pooled fit includes them. L4. Compilation bias. AnAge reports one number per species with no measure of within-species uncertainty, so we cannot weight by precision and we cannot separate measurement error from biological variation. L5. Mass range. Even 6.7 decades is small next to the 20-odd decades that span bacteria to whales. Nothing here speaks to unicellular life or to invertebrates. L6. One number per species means intraspecific scaling is invisible. Heusner (1982) argued the interspecific 3/4 is an artefact of pooling intraspecific 2/3 relationships; our data cannot test that. Seed: 20260614 (used only for the bootstrap resampling; every other number in this file is deterministic and seed-independent). Run: python allometry-scaling.py > allometry-scaling-output.txt """ import csv import hashlib import math import os import sys import time import numpy as np # ---------------------------------------------------------------------------- # Constants # ---------------------------------------------------------------------------- SEED = 20260614 N_BOOT = 20000 ALPHA = 0.05 THREE_QUARTERS = 0.75 TWO_THIRDS = 2.0 / 3.0 HERE = os.path.dirname(os.path.abspath(__file__)) DATA_DIR = os.path.join(HERE, "data") DATA_TXT = os.path.join(DATA_DIR, "anage_data.txt") DATA_ZIP = os.path.join(DATA_DIR, "anage_dataset.zip") DATA_URL = "https://genomics.senescence.info/species/dataset.zip" RETRIEVED = "2026-09-14" if hasattr(sys.stdout, "reconfigure"): sys.stdout.reconfigure(encoding="utf-8") T0 = time.time() def rule(ch="=", n=78): print(ch * n) def head(title): print() rule() print(title) rule() # ============================================================================ # SECTION 0. Distribution functions, built from scratch, then validated. # ============================================================================ def betacf(a, b, x): """Continued fraction for the incomplete beta function (Lentz's method).""" MAXIT, EPS, FPMIN = 300, 3.0e-16, 1.0e-300 qab, qap, qam = a + b, a + 1.0, a - 1.0 c = 1.0 d = 1.0 - qab * x / qap if abs(d) < FPMIN: d = FPMIN d = 1.0 / d h = d for m in range(1, MAXIT + 1): m2 = 2 * m aa = m * (b - m) * x / ((qam + m2) * (a + m2)) d = 1.0 + aa * d if abs(d) < FPMIN: d = FPMIN c = 1.0 + aa / c if abs(c) < FPMIN: c = FPMIN 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) < FPMIN: d = FPMIN c = 1.0 + aa / c if abs(c) < FPMIN: c = FPMIN d = 1.0 / d dl = d * c h *= dl if abs(dl - 1.0) < EPS: break return h def betainc(a, b, x): """Regularised incomplete beta I_x(a,b).""" if x <= 0.0: return 0.0 if x >= 1.0: return 1.0 lbeta = (math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) + a * math.log(x) + b * math.log(1.0 - x)) 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 def t_sf2(t, df): """Two-sided tail probability of Student's t.""" t = abs(float(t)) return betainc(df / 2.0, 0.5, df / (df + t * t)) def t_ppf(p, df): """Inverse Student t CDF by bisection. p in (0,1).""" lo, hi = -400.0, 400.0 for _ in range(200): mid = 0.5 * (lo + hi) cdf = 1.0 - 0.5 * t_sf2(mid, df) if mid >= 0 else 0.5 * t_sf2(mid, df) if cdf < p: lo = mid else: hi = mid return 0.5 * (lo + hi) def f_sf(f, df1, df2): """Upper tail probability of the F distribution.""" if f <= 0: return 1.0 return betainc(df2 / 2.0, df1 / 2.0, df2 / (df2 + df1 * f)) def tcrit(df, alpha=ALPHA): return t_ppf(1.0 - alpha / 2.0, df) # ============================================================================ # SECTION 0b. The two estimators, written out longhand. # ============================================================================ def ols(x, y): """Ordinary least squares of y on x. Returns a dict of everything.""" x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) n = x.size mx, my = x.mean(), y.mean() dx, dy = x - mx, y - my Sxx = float(dx @ dx) Syy = float(dy @ dy) Sxy = float(dx @ dy) b = Sxy / Sxx a = my - b * mx resid = y - (a + b * x) SSE = float(resid @ resid) df = n - 2 s2 = SSE / df se_b = math.sqrt(s2 / Sxx) se_a = math.sqrt(s2 * (1.0 / n + mx * mx / Sxx)) r = Sxy / math.sqrt(Sxx * Syy) tc = tcrit(df) return dict(n=n, b=b, a=a, se_b=se_b, se_a=se_a, r=r, r2=r * r, SSE=SSE, Sxx=Sxx, Syy=Syy, Sxy=Sxy, df=df, s2=s2, lo=b - tc * se_b, hi=b + tc * se_b, resid=resid, mx=mx, my=my) def rma(x, y, alpha=ALPHA): """Reduced major axis (standardised major axis) slope and interval.""" x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) n = x.size sdx = x.std(ddof=1) sdy = y.std(ddof=1) r = float(np.corrcoef(x, y)[0, 1]) b = math.copysign(1.0, r) * sdy / sdx a = y.mean() - b * x.mean() df = n - 2 Fc = tcrit(df, alpha) ** 2 # F_{1,df} = t^2_{df} B = Fc * (1.0 - r * r) / df lo = b * (math.sqrt(B + 1.0) - math.sqrt(B)) hi = b * (math.sqrt(B + 1.0) + math.sqrt(B)) return dict(n=n, b=b, a=a, r=r, r2=r * r, lo=lo, hi=hi, df=df, B=B) def rma_test(x, y, beta0): """Exact test of H0: RMA slope = beta0, via the zero-correlation test on the rescaled residual and fitted axis scores.""" x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) / float(beta0) u = y - x v = y + x n = x.size r = float(np.corrcoef(u, v)[0, 1]) df = n - 2 t = r * math.sqrt(df) / math.sqrt(max(1.0 - r * r, 1e-300)) return t, t_sf2(t, df), df def verdict(lo, hi): """Which of the two candidate exponents does this interval exclude?""" e34 = not (lo <= THREE_QUARTERS <= hi) e23 = not (lo <= TWO_THIRDS <= hi) if e34 and e23: return "excludes BOTH 3/4 and 2/3" if e34: return "excludes 3/4, contains 2/3" if e23: return "excludes 2/3, contains 3/4" return "excludes NEITHER" # ============================================================================ # SECTION 1. VALIDATION # ============================================================================ head("SECTION 1. VALIDATION AGAINST CLOSED FORM") print(""" Nothing below this line is worth reading unless the two regression routines are right. So they are checked first, on a five-point dataset small enough to do by hand, against closed-form answers computed independently with exact rational arithmetic (Python's fractions module, no floating point at all). """) from fractions import Fraction as F # noqa: E402 (deliberately local to ยง1) vx = [0, 1, 2, 3, 4] vy = [1, 3, 2, 5, 4] # ---- exact closed form, rational arithmetic ------------------------------- n_v = len(vx) fx = [F(v) for v in vx] fy = [F(v) for v in vy] mxf = sum(fx) / n_v myf = sum(fy) / n_v Sxx_f = sum((v - mxf) ** 2 for v in fx) Syy_f = sum((v - myf) ** 2 for v in fy) Sxy_f = sum((a_ - mxf) * (b_ - myf) for a_, b_ in zip(fx, fy)) b_f = Sxy_f / Sxx_f a_f = myf - b_f * mxf res_f = [yy - (a_f + b_f * xx) for xx, yy in zip(fx, fy)] SSE_f = sum(rr ** 2 for rr in res_f) s2_f = SSE_f / (n_v - 2) seb_f = math.sqrt(float(s2_f / Sxx_f)) r2_f = Sxy_f ** 2 / (Sxx_f * Syy_f) r_f = math.sqrt(float(r2_f)) brma_f = math.sqrt(float(Syy_f / Sxx_f)) # sd(y)/sd(x), n-1 cancels print("Hand dataset: x = %s" % vx) print(" y = %s" % vy) print() print("Exact sums (rational): Sxx = %s Syy = %s Sxy = %s" % (Sxx_f, Syy_f, Sxy_f)) print("Exact OLS slope b = Sxy/Sxx = %s = %.12f" % (b_f, float(b_f))) print("Exact OLS intercept a = ybar - b*xbar = %s = %.12f" % (a_f, float(a_f))) print("Exact SSE = Syy - b*Sxy = %s = %.12f" % (SSE_f, float(SSE_f))) print("Exact r^2 = Sxy^2/(Sxx*Syy) = %s = %.12f" % (r2_f, float(r2_f))) print("Exact RMA slope = sqrt(Syy/Sxx) = sqrt(%s) = %.12f" % (Syy_f / Sxx_f, brma_f)) print("Exact RMA via identity b_ols/|r| = %.12f" % (float(b_f) / r_f)) V = ols(vx, vy) R = rma(vx, vy) print() print("%-34s %18s %18s %13s" % ("quantity", "club code", "closed form", "difference")) print("-" * 86) rows_val = [ ("OLS slope b", V["b"], float(b_f)), ("OLS intercept a", V["a"], float(a_f)), ("residual sum of squares", V["SSE"], float(SSE_f)), ("r^2", V["r2"], float(r2_f)), ("SE(b)", V["se_b"], seb_f), ("RMA slope", R["b"], brma_f), ("RMA slope via b_ols/|r|", R["b"], float(b_f) / r_f), ] maxdiff = 0.0 for name, got, want in rows_val: d = got - want maxdiff = max(maxdiff, abs(d)) print("%-34s %18.12f %18.12f %13.2e" % (name, got, want, d)) print("-" * 86) print("largest absolute discrepancy: %.3e" % maxdiff) assert maxdiff < 1e-12, "regression code disagrees with the closed form" print() print("Residuals, printed in full, club code against exact rational values:") print("%6s %10s %10s %16s %16s %13s" % ("i", "x", "y", "club resid", "exact resid", "difference")) for i in range(n_v): print("%6d %10.4f %10.4f %16.12f %16.12f %13.2e" % (i, vx[i], vy[i], V["resid"][i], float(res_f[i]), V["resid"][i] - float(res_f[i]))) print("sum of club residuals = %.3e (must be 0 for a fit with an intercept)" % float(V["resid"].sum())) print("sum of x*resid = %.3e (must be 0: the normal equations)" % float(np.asarray(vx, dtype=float) @ V["resid"])) # ---- independent numpy cross-check --------------------------------------- print() np_b, np_a = np.polyfit(vx, vy, 1) A = np.vstack([np.ones(n_v), np.asarray(vx, dtype=float)]).T ls = np.linalg.lstsq(A, np.asarray(vy, dtype=float), rcond=None)[0] print("numpy.polyfit slope %.12f intercept %.12f" % (np_b, np_a)) print("numpy.linalg.lstsq slope %.12f intercept %.12f" % (ls[1], ls[0])) print("club code slope %.12f intercept %.12f" % (V["b"], V["a"])) print("max difference vs numpy: %.3e" % max(abs(np_b - V["b"]), abs(ls[1] - V["b"]), abs(np_a - V["a"]))) # ---- degenerate case: perfect correlation must make OLS and RMA identical -- print() px = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]) py = 3.0 + 0.75 * px Vp, Rp = ols(px, py), rma(px, py) print("Perfectly collinear check, y = 3 + 0.75x exactly:") print(" OLS slope %.15f RMA slope %.15f difference %.3e" % (Vp["b"], Rp["b"], Vp["b"] - Rp["b"])) print(" r^2 = %.15f, so the OLS/RMA gap 1/|r| - 1 = %.3e" % (Vp["r2"], 1.0 / abs(Vp["r"]) - 1.0)) # ---- the t and F machinery ------------------------------------------------ print() print("Student t critical values, two-sided 95%, club code against published") print("tables (Fisher & Yates / any standard statistical table):") print("%8s %18s %18s %13s" % ("df", "club code", "table value", "difference")) tab = [(1, 12.706205), (2, 4.302653), (3, 3.182446), (5, 2.570582), (10, 2.228139), (20, 2.085963), (30, 2.042272), (60, 2.000298), (120, 1.979930)] tmax = 0.0 for df_, want in tab: got = tcrit(df_) tmax = max(tmax, abs(got - want)) print("%8d %18.6f %18.6f %13.2e" % (df_, got, want, got - want)) print("largest discrepancy: %.2e" % tmax) print("(table values are quoted to six decimal places, so a few units in the") print("last place is the closest agreement the comparison can show)") assert tmax < 1e-5 # large-df limit: t must approach the normal quantile 1.959964 from above, at # the rate given by the standard Cornish-Fisher term z(1 + (z^2+1)/(4v)). z975 = 1.9599639845400545 for df_ in [1000, 10000, 100000]: approx = z975 * (1.0 + (z975 * z975 + 1.0) / (4.0 * df_)) print(" df = %-7d club t = %.8f z(1+(z^2+1)/4v) = %.8f diff %.2e" % (df_, tcrit(df_), approx, tcrit(df_) - approx)) print(" the club value converges on the normal quantile %.6f from above," % z975) print(" which is the behaviour a correct t quantile must show.") print() print("F distribution check: for df1 = 1, F = t^2 exactly, so the upper tail") print("of F must equal the two-sided tail of t.") for tt, dd in [(2.0, 10), (3.182446, 3), (1.5, 420)]: print(" t = %-9.6f df = %-6d 2-sided t tail %.10f F=%.6f tail %.10f diff %.2e" % (tt, dd, t_sf2(tt, dd), tt * tt, f_sf(tt * tt, 1, dd), t_sf2(tt, dd) - f_sf(tt * tt, 1, dd))) print(" I_0.5(1,1) = %.15f (must be exactly 0.5)" % betainc(1.0, 1.0, 0.5)) print(" I_0.25(2,3) = %.12f (closed form 1 - (1-x)^3(1+3x) at x=0.25 " "= %.12f)" % (betainc(2.0, 3.0, 0.25), 1 - (0.75 ** 3) * (1 + 3 * 0.25))) # ---- RMA interval and RMA point test must agree --------------------------- print() print("Internal consistency: the RMA confidence interval and the RMA point") print("hypothesis test are computed by two different routes, so testing the") print("slope value sitting exactly on the interval endpoint must return") print("p = 0.05 to within rounding.") rng_chk = np.random.default_rng(11) cx = rng_chk.normal(0, 1, 60) cy = 0.7 * cx + rng_chk.normal(0, 0.35, 60) Rc = rma(cx, cy) for edge, nm in [(Rc["lo"], "lower"), (Rc["hi"], "upper")]: tt, pp, _ = rma_test(cx, cy, edge) print(" %s endpoint b = %.6f p = %.6f |p - 0.05| = %.2e" % (nm, edge, pp, abs(pp - 0.05))) print() print("VALIDATION PASSED. Both estimators reproduce closed-form answers to") print("better than 1e-12 (largest %.2e), and the t quantiles match published" % maxdiff) print("tables to better than 1e-5 (largest %.2e)." % tmax) # ============================================================================ # SECTION 2. DATA # ============================================================================ head("SECTION 2. THE DATA") if not os.path.exists(DATA_TXT): print("FATAL: cached data file missing at %s" % DATA_TXT) print("Download %s and extract anage_data.txt into that directory." % DATA_URL) sys.exit(1) if os.path.exists(DATA_ZIP): with open(DATA_ZIP, "rb") as fh: zsha = hashlib.sha256(fh.read()).hexdigest() else: zsha = "(archive not cached)" with open(DATA_TXT, "rb") as fh: tsha = hashlib.sha256(fh.read()).hexdigest() print("source : AnAge, Animal Ageing and Longevity Database, Build 15") print("release date : 3 July 2023") print("url : %s" % DATA_URL) print("retrieved : %s" % RETRIEVED) print("cached archive: analysis/data/anage_dataset.zip") print(" sha256 : %s" % zsha) print("cached table : analysis/data/anage_data.txt") print(" sha256 : %s" % tsha) print("citation : Tacutu et al. (2018) Nucleic Acids Research 46, D1083") print() recs = [] total_rows = 0 with open(DATA_TXT, encoding="utf-8", errors="replace") as fh: for row in csv.DictReader(fh, delimiter="\t"): total_rows += 1 mr = row["Metabolic rate (W)"].strip() bm = row["Body mass (g)"].strip() if not mr or not bm: continue try: m = float(bm) w = float(mr) except ValueError: continue if m <= 0 or w <= 0: continue recs.append(dict(cls=row["Class"].strip(), order=row["Order"].strip(), sp=(row["Genus"].strip() + " " + row["Species"].strip()), common=row["Common name"].strip(), mass=m, bmr=w)) print("rows in AnAge build 15 : %d" % total_rows) print("rows carrying BOTH mass and BMR > 0 : %d" % len(recs)) names = [r["sp"] for r in recs] print("distinct species names : %d" % len(set(names))) print("duplicate species names : %d" % (len(names) - len(set(names)))) mass = np.array([r["mass"] for r in recs]) bmr = np.array([r["bmr"] for r in recs]) X = np.log10(mass) Y = np.log10(bmr) CLS = np.array([r["cls"] for r in recs]) ORD = np.array([r["order"] for r in recs]) SP = np.array([r["sp"] for r in recs]) print() print("mass range : %.2f g to %.0f g" % (mass.min(), mass.max())) print(" : %.3f decades (a factor of %.3g)" % (X.max() - X.min(), mass.max() / mass.min())) print("BMR range : %.3g W to %.4g W (%.3f decades)" % (bmr.min(), bmr.max(), Y.max() - Y.min())) print() print("breakdown by class:") print("%-14s %6s %10s %14s %10s" % ("class", "n", "decades", "median mass g", "share")) CLASSES = ["Mammalia", "Aves", "Reptilia", "Amphibia"] for c in CLASSES: sel = CLS == c print("%-14s %6d %10.3f %14.1f %9.1f%%" % (c, sel.sum(), X[sel].max() - X[sel].min(), float(np.median(mass[sel])), 100.0 * sel.sum() / len(recs))) print("%-14s %6d %10.3f %14.1f %9.1f%%" % ("TOTAL", len(recs), X.max() - X.min(), float(np.median(mass)), 100.0)) print() print("lightest five and heaviest five, as a sanity check that the units are") print("what we think they are:") idx = np.argsort(mass) print("%-32s %-12s %14s %14s" % ("species", "class", "mass (g)", "BMR (W)")) for i in list(idx[:5]) + list(idx[-5:]): print("%-32s %-12s %14.4g %14.5g" % (SP[i], CLS[i], mass[i], bmr[i])) print() print("A human check: our mass-specific BMR for Homo sapiens should land near") print("the textbook 1.2 W/kg resting.") hi_ = np.where(SP == "Homo sapiens")[0] if hi_.size: i = hi_[0] print(" Homo sapiens: mass %.0f g, BMR %.2f W, that is %.2f W/kg and %d kcal/day" % (mass[i], bmr[i], 1000.0 * bmr[i] / mass[i], round(bmr[i] * 86400.0 / 4184.0))) # ============================================================================ # SECTION 3. THE HEADLINE FITS # ============================================================================ head("SECTION 3. EXPONENTS BY BOTH METHODS, WITH INTERVALS") GROUPS = [ ("Mammalia", CLS == "Mammalia"), ("Aves", CLS == "Aves"), ("Reptilia", CLS == "Reptilia"), ("Amphibia", CLS == "Amphibia"), ("Endotherms (M+A)", (CLS == "Mammalia") | (CLS == "Aves")), ("Ectotherms (R+A)", (CLS == "Reptilia") | (CLS == "Amphibia")), ("All four classes", np.ones(len(recs), dtype=bool)), ] fits = {} print("%-18s %5s %7s %7s %8s %8s %17s %7s" % ("group", "n", "dec", "r^2", "method", "b", "95% CI", "SE")) print("-" * 92) for name, sel in GROUPS: xs, ys = X[sel], Y[sel] o = ols(xs, ys) m = rma(xs, ys) fits[name] = dict(ols=o, rma=m, sel=sel, dec=float(xs.max() - xs.min())) print("%-18s %5d %7.2f %7.4f %8s %8.4f [%.4f, %.4f] %7.4f" % (name, o["n"], xs.max() - xs.min(), o["r2"], "OLS", o["b"], o["lo"], o["hi"], o["se_b"])) print("%-18s %5s %7s %7s %8s %8.4f [%.4f, %.4f] %7s" % ("", "", "", "", "RMA", m["b"], m["lo"], m["hi"], "-")) print("-" * 92) print() print("The raw sums behind the headline mammal fit, printed so the division") print("can be done on paper:") _mo = fits["Mammalia"]["ols"] print(" n = %d" % _mo["n"]) print(" xbar = %.6f (log10 grams)" % _mo["mx"]) print(" ybar = %.6f (log10 watts)" % _mo["my"]) print(" Sxx = %.4f" % _mo["Sxx"]) print(" Sxy = %.4f" % _mo["Sxy"]) print(" Syy = %.4f" % _mo["Syy"]) print(" b = Sxy/Sxx = %.4f / %.4f = %.6f" % (_mo["Sxy"], _mo["Sxx"], _mo["b"])) print(" r = Sxy/sqrt(Sxx*Syy) = %.6f" % _mo["r"]) print(" b_rma = sqrt(Syy/Sxx) = sqrt(%.4f/%.4f) = %.6f" % (_mo["Syy"], _mo["Sxx"], math.sqrt(_mo["Syy"] / _mo["Sxx"]))) print(" SSE = Syy - b*Sxy = %.4f" % (_mo["Syy"] - _mo["b"] * _mo["Sxy"])) print(" SE(b) = sqrt( SSE/(n-2) / Sxx ) = sqrt(%.6f/%.4f) = %.6f" % (_mo["SSE"] / _mo["df"], _mo["Sxx"], _mo["se_b"])) print(" t(3/4)= (%.6f - 0.75)/%.6f = %.4f on %d df" % (_mo["b"], _mo["se_b"], (_mo["b"] - 0.75) / _mo["se_b"], _mo["df"])) print() print("Intercepts (log10 a, with B in watts and M in grams), and the") print("predicted BMR of a 1 kg animal, which is what the intercept means:") print("%-18s %10s %10s %14s" % ("group", "OLS log a", "RMA log a", "BMR at 1 kg, W")) for name, _ in GROUPS: o, m = fits[name]["ols"], fits[name]["rma"] print("%-18s %10.4f %10.4f %14.4f" % (name, o["a"], m["a"], 10 ** (o["a"] + o["b"] * 3.0))) # ============================================================================ # SECTION 4. TESTING 3/4 AND 2/3 EXPLICITLY # ============================================================================ head("SECTION 4. DOES THE INTERVAL CONTAIN 3/4? DOES IT CONTAIN 2/3?") print("3/4 = %.6f 2/3 = %.6f difference = %.6f" % (THREE_QUARTERS, TWO_THIRDS, THREE_QUARTERS - TWO_THIRDS)) print() print("OLS, t = (b - beta0)/SE(b), df = n-2:") print("%-18s %8s %9s %8s %10s %9s %8s %10s" % ("group", "b", "SE", "t(3/4)", "p(3/4)", "t(2/3)", "p(2/3)", "verdict")) print("-" * 100) for name, _ in GROUPS: o = fits[name]["ols"] t34 = (o["b"] - THREE_QUARTERS) / o["se_b"] t23 = (o["b"] - TWO_THIRDS) / o["se_b"] print("%-18s %8.4f %9.4f %8.2f %10.3g %9.2f %8.3g %s" % (name, o["b"], o["se_b"], t34, t_sf2(t34, o["df"]), t23, t_sf2(t23, o["df"]), verdict(o["lo"], o["hi"]))) print() print("RMA, exact zero-correlation test on the rescaled axis scores:") print("%-18s %8s %8s %10s %9s %10s %8s" % ("group", "b", "t(3/4)", "p(3/4)", "t(2/3)", "p(2/3)", "verdict")) print("-" * 100) for name, sel in GROUPS: m = fits[name]["rma"] t34, p34, _ = rma_test(X[sel], Y[sel], THREE_QUARTERS) t23, p23, _ = rma_test(X[sel], Y[sel], TWO_THIRDS) print("%-18s %8.4f %8.2f %10.3g %9.2f %10.3g %s" % (name, m["b"], t34, p34, t23, p23, verdict(m["lo"], m["hi"]))) print() print("Distance of each fitted exponent from the two candidates, in standard") print("errors of that fit (OLS only, since OLS is the one with a closed-form SE):") print("%-18s %12s %12s" % ("group", "SEs from 3/4", "SEs from 2/3")) for name, _ in GROUPS: o = fits[name]["ols"] print("%-18s %12.2f %12.2f" % (name, (o["b"] - THREE_QUARTERS) / o["se_b"], (o["b"] - TWO_THIRDS) / o["se_b"])) # ============================================================================ # SECTION 5. THE SIZE OF THE METHOD EFFECT # ============================================================================ head("SECTION 5. HOW FAR APART THE TWO METHODS SIT, ON THE SAME DATA") print("Algebraically b_rma = b_ols / |r| exactly, so the gap depends only on") print("the scatter. Confirm the identity numerically, then read the size off.") print() print("%-18s %8s %9s %9s %9s %9s %10s" % ("group", "r^2", "b_ols", "b_rma", "b_ols/|r|", "gap", "gap %")) print("-" * 84) ident_max = 0.0 for name, _ in GROUPS: o, m = fits[name]["ols"], fits[name]["rma"] pred = o["b"] / abs(o["r"]) ident_max = max(ident_max, abs(pred - m["b"])) print("%-18s %8.4f %9.4f %9.4f %9.4f %9.4f %9.2f%%" % (name, o["r2"], o["b"], m["b"], pred, m["b"] - o["b"], 100.0 * (m["b"] / o["b"] - 1.0))) print("-" * 84) print("largest deviation from the identity b_rma = b_ols/|r| : %.3e" % ident_max) print() print("Read that column again. The OLS/RMA gap for mammals is %.4f, which is" % (fits["Mammalia"]["rma"]["b"] - fits["Mammalia"]["ols"]["b"])) print("%.1f%% of the whole distance between 3/4 and 2/3 (%.4f)." % (100.0 * (fits["Mammalia"]["rma"]["b"] - fits["Mammalia"]["ols"]["b"]) / (THREE_QUARTERS - TWO_THIRDS), THREE_QUARTERS - TWO_THIRDS)) gap_all = fits["All four classes"]["rma"]["b"] - fits["All four classes"]["ols"]["b"] print("For all four classes pooled it is %.4f, which is %.0f%% of that" % (gap_all, 100.0 * gap_all / (THREE_QUARTERS - TWO_THIRDS))) print("distance: bigger than the entire quantity under dispute.") print() print("What r^2 would a dataset need for the two methods to agree to within") print("a given tolerance, at b_ols = 0.70?") print("%12s %14s %14s" % ("tolerance", "required |r|", "required r^2")) for tol in [0.0833, 0.02, 0.01, 0.005, 0.001]: need_r = 0.70 / (0.70 + tol) print("%12.4f %14.5f %14.5f" % (tol, need_r, need_r ** 2)) # ============================================================================ # SECTION 6. BOOTSTRAP, AS A CHECK ON THE ANALYTIC INTERVALS # ============================================================================ head("SECTION 6. BOOTSTRAP INTERVALS (seed %d, %d resamples)" % (SEED, N_BOOT)) rng = np.random.default_rng(SEED) print("Non-parametric bootstrap over species, %d resamples per group." % N_BOOT) print("If the analytic intervals are right, these should agree closely. They") print("do NOT fix the phylogenetic non-independence: resampling species keeps") print("the same shared ancestry, so both intervals are too narrow together.") print() print("%-18s %9s %19s %19s %9s" % ("group", "method", "analytic 95% CI", "bootstrap 95% CI", "boot SE")) print("-" * 84) boot_store = {} for name, sel in GROUPS: xs, ys = X[sel], Y[sel] n = xs.size if n < 8: continue idx_b = rng.integers(0, n, size=(N_BOOT, n)) bx = xs[idx_b] by = ys[idx_b] mxb = bx.mean(axis=1, keepdims=True) myb = by.mean(axis=1, keepdims=True) dxb = bx - mxb dyb = by - myb Sxxb = (dxb * dxb).sum(axis=1) Syyb = (dyb * dyb).sum(axis=1) Sxyb = (dxb * dyb).sum(axis=1) good = Sxxb > 0 b_ols_b = Sxyb[good] / Sxxb[good] rb = Sxyb[good] / np.sqrt(Sxxb[good] * Syyb[good]) b_rma_b = np.sign(rb) * np.sqrt(Syyb[good] / Sxxb[good]) boot_store[name] = (b_ols_b, b_rma_b) o, m = fits[name]["ols"], fits[name]["rma"] qo = np.percentile(b_ols_b, [2.5, 97.5]) qm = np.percentile(b_rma_b, [2.5, 97.5]) print("%-18s %9s [%.4f, %.4f] [%.4f, %.4f] %9.4f" % (name, "OLS", o["lo"], o["hi"], qo[0], qo[1], b_ols_b.std(ddof=1))) print("%-18s %9s [%.4f, %.4f] [%.4f, %.4f] %9.4f" % ("", "RMA", m["lo"], m["hi"], qm[0], qm[1], b_rma_b.std(ddof=1))) print() print("Bootstrap probability that the true exponent exceeds each candidate,") print("read straight off the resample distribution:") print("%-18s %9s %14s %14s" % ("group", "method", "P(b > 3/4)", "P(b > 2/3)")) for name in ["Mammalia", "Aves", "Endotherms (M+A)", "All four classes"]: bo, bm = boot_store[name] print("%-18s %9s %14.4f %14.4f" % (name, "OLS", (bo > THREE_QUARTERS).mean(), (bo > TWO_THIRDS).mean())) print("%-18s %9s %14.4f %14.4f" % ("", "RMA", (bm > THREE_QUARTERS).mean(), (bm > TWO_THIRDS).mean())) # ============================================================================ # SECTION 7. DOES THE EXPONENT DIFFER BETWEEN TAXONOMIC GROUPS? # ============================================================================ head("SECTION 7. TAXONOMIC HETEROGENEITY") def ancova_slopes(xs, ys, groups): """Compare a common-slope model against a separate-slopes model.""" labs = sorted(set(groups)) k = len(labs) n = xs.size # separate slopes: sum of within-group SSE, df = n - 2k sse_sep = 0.0 for g in labs: s = groups == g sse_sep += ols(xs[s], ys[s])["SSE"] df_sep = n - 2 * k # common slope, separate intercepts (ANCOVA): df = n - (k+1) Sxx_w = 0.0 Sxy_w = 0.0 for g in labs: s = groups == g dx = xs[s] - xs[s].mean() dy = ys[s] - ys[s].mean() Sxx_w += float(dx @ dx) Sxy_w += float(dx @ dy) b_common = Sxy_w / Sxx_w sse_com = 0.0 for g in labs: s = groups == g pred = ys[s].mean() + b_common * (xs[s] - xs[s].mean()) rr = ys[s] - pred sse_com += float(rr @ rr) df_com = n - (k + 1) df1 = df_com - df_sep Fstat = ((sse_com - sse_sep) / df1) / (sse_sep / df_sep) return dict(labs=labs, k=k, b_common=b_common, F=Fstat, df1=df1, df2=df_sep, p=f_sf(Fstat, df1, df_sep), sse_sep=sse_sep, sse_com=sse_com, se_common=math.sqrt((sse_com / df_com) / Sxx_w)) sel4 = np.ones(len(recs), dtype=bool) A4 = ancova_slopes(X, Y, CLS) print("All four classes, separate slopes against one common slope:") print(" classes : %s" % ", ".join(A4["labs"])) print(" common slope : %.4f (SE %.4f)" % (A4["b_common"], A4["se_common"])) print(" SSE separate slopes: %.4f on %d df" % (A4["sse_sep"], A4["df2"])) print(" SSE common slope : %.4f on %d df" % (A4["sse_com"], A4["df2"] + A4["df1"])) print(" F(%d, %d) = %.3f, p = %.4g" % (A4["df1"], A4["df2"], A4["F"], A4["p"])) print(" verdict : %s" % ("slopes differ between classes" if A4["p"] < 0.05 else "no detectable difference in slope between classes")) selE = (CLS == "Mammalia") | (CLS == "Aves") AE = ancova_slopes(X[selE], Y[selE], CLS[selE]) print() print("Mammals against birds only (the two well-sampled classes):") print(" common slope : %.4f (SE %.4f)" % (AE["b_common"], AE["se_common"])) print(" F(%d, %d) = %.3f, p = %.4g" % (AE["df1"], AE["df2"], AE["F"], AE["p"])) print() print("Pairwise slope contrasts, OLS, Welch-style on the two standard errors:") print("%-22s %9s %9s %9s %8s %10s" % ("pair", "b1", "b2", "b1-b2", "t", "p")) print("-" * 74) pairs = [("Mammalia", "Aves"), ("Mammalia", "Reptilia"), ("Aves", "Reptilia"), ("Mammalia", "Amphibia"), ("Reptilia", "Amphibia")] for g1, g2 in pairs: o1, o2 = fits[g1]["ols"], fits[g2]["ols"] d = o1["b"] - o2["b"] se = math.sqrt(o1["se_b"] ** 2 + o2["se_b"] ** 2) dfw = ((o1["se_b"] ** 2 + o2["se_b"] ** 2) ** 2 / (o1["se_b"] ** 4 / o1["df"] + o2["se_b"] ** 4 / o2["df"])) print("%-22s %9.4f %9.4f %9.4f %8.2f %10.3g" % (g1 + " vs " + g2, o1["b"], o2["b"], d, d / se, t_sf2(d / se, dfw))) print() print("Mammalian orders with at least 8 species, fitted separately. This is") print("the same question asked one level down.") print("%-20s %5s %7s %9s %19s %8s" % ("order", "n", "decades", "b (OLS)", "95% CI", "r^2")) print("-" * 76) mam = CLS == "Mammalia" order_rows = [] for o_ in sorted(set(ORD[mam])): s = mam & (ORD == o_) if s.sum() < 8: continue xs, ys = X[s], Y[s] if xs.max() - xs.min() < 0.5: continue f_ = ols(xs, ys) order_rows.append((o_, f_, float(xs.max() - xs.min()))) print("%-20s %5d %7.2f %9.4f [%.4f, %.4f] %8.4f" % (o_, f_["n"], xs.max() - xs.min(), f_["b"], f_["lo"], f_["hi"], f_["r2"])) ob = np.array([r[1]["b"] for r in order_rows]) print("-" * 76) print("within-order slopes: n = %d orders, mean %.4f, sd %.4f, range %.4f to %.4f" % (ob.size, ob.mean(), ob.std(ddof=1), ob.min(), ob.max())) print("how many of those %d order-level intervals contain 3/4: %d" % (len(order_rows), sum(1 for r in order_rows if r[1]["lo"] <= 0.75 <= r[1]["hi"]))) print("how many contain 2/3 : %d" % sum(1 for r in order_rows if r[1]["lo"] <= TWO_THIRDS <= r[1]["hi"])) print("how many contain both : %d" % sum(1 for r in order_rows if r[1]["lo"] <= TWO_THIRDS and r[1]["hi"] >= 0.75)) print() print("Order-level means (one point per mammalian order, n >= 3 species).") print("This is the crudest possible answer to the pseudoreplication problem:") print("it throws away nearly all the data to buy near-independence.") oxs, oys, onames = [], [], [] for o_ in sorted(set(ORD[mam])): s = mam & (ORD == o_) if s.sum() < 3: continue oxs.append(float(X[s].mean())) oys.append(float(Y[s].mean())) onames.append((o_, int(s.sum()))) oxs = np.array(oxs) oys = np.array(oys) Fo = ols(oxs, oys) Ro = rma(oxs, oys) print(" orders used : %d (%s)" % (len(onames), ", ".join("%s n=%d" % t for t in onames))) print(" OLS slope : %.4f 95%% CI [%.4f, %.4f] %s" % (Fo["b"], Fo["lo"], Fo["hi"], verdict(Fo["lo"], Fo["hi"]))) print(" RMA slope : %.4f 95%% CI [%.4f, %.4f] %s" % (Ro["b"], Ro["lo"], Ro["hi"], verdict(Ro["lo"], Ro["hi"]))) print(" r^2 = %.4f, n = %d, so the interval is %.1f times wider than the" % (Fo["r2"], Fo["n"], (Fo["hi"] - Fo["lo"]) / (fits["Mammalia"]["ols"]["hi"] - fits["Mammalia"]["ols"]["lo"]))) print(" species-level one. Losing independence costs precision, and that is") print(" the honest trade.") # ============================================================================ # SECTION 8. RANGE RESTRICTION # ============================================================================ head("SECTION 8. WHAT HAPPENS WHEN YOU NARROW THE MASS WINDOW") print("The core statistical criticism of Kleiber's law is that the exponent") print("is estimated across an enormous mass range, and that a line fitted") print("across six decades will look straight whatever the data do locally.") print("So: fit the same mammal data inside sliding windows and watch.") print() xm, ym = X[mam], Y[mam] WIN = 2.0 step = 0.25 lo0 = math.floor(xm.min() * 4) / 4 hi0 = math.ceil(xm.max() * 4) / 4 print("Sliding window of width %.1f decades, step %.2f decades, mammals only." % (WIN, step)) print("%9s %9s %6s %9s %19s %8s %8s %s" % ("win lo", "win hi", "n", "b (OLS)", "95% CI", "b (RMA)", "r^2", "verdict")) print("-" * 104) slide = [] c = lo0 while c + WIN <= hi0 + 1e-9: s = (xm >= c) & (xm < c + WIN) if s.sum() >= 12: f_ = ols(xm[s], ym[s]) m_ = rma(xm[s], ym[s]) slide.append((c, c + WIN, f_, m_)) print("%9.2f %9.2f %6d %9.4f [%.4f, %.4f] %8.4f %8.4f %s" % (c, c + WIN, f_["n"], f_["b"], f_["lo"], f_["hi"], m_["b"], f_["r2"], verdict(f_["lo"], f_["hi"]))) c += step print("-" * 104) sb = np.array([s[2]["b"] for s in slide]) print("window slopes: n = %d windows, min %.4f, max %.4f, range %.4f, sd %.4f" % (sb.size, sb.min(), sb.max(), sb.max() - sb.min(), sb.std(ddof=1))) print("that range is %.2f times the whole 3/4-minus-2/3 gap of %.4f" % ((sb.max() - sb.min()) / (THREE_QUARTERS - TWO_THIRDS), THREE_QUARTERS - TWO_THIRDS)) print("windows whose interval excludes 3/4 : %d of %d" % (sum(1 for s in slide if not (s[2]["lo"] <= 0.75 <= s[2]["hi"])), len(slide))) print("windows whose interval excludes 2/3 : %d of %d" % (sum(1 for s in slide if not (s[2]["lo"] <= TWO_THIRDS <= s[2]["hi"])), len(slide))) print("windows whose interval excludes both: %d of %d" % (sum(1 for s in slide if not (s[2]["lo"] <= 0.75 <= s[2]["hi"]) and not (s[2]["lo"] <= TWO_THIRDS <= s[2]["hi"])), len(slide))) print("windows whose interval excludes neither: %d of %d" % (sum(1 for s in slide if (s[2]["lo"] <= 0.75 <= s[2]["hi"]) and (s[2]["lo"] <= TWO_THIRDS <= s[2]["hi"])), len(slide))) print() print("Nested truncation, mammals. Left column drops the heaviest species,") print("right column drops the lightest. Same data, same code, one cut.") print("%11s %6s %8s %18s %4s %11s %6s %8s %18s" % ("keep M <= g", "n", "b", "95% CI", "|", "keep M >= g", "n", "b", "95% CI")) print("-" * 108) lows = [2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0] highs = [-0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5] trunc_lo, trunc_hi = [], [] for cu, cd in zip(lows, highs): s1 = xm <= cu s2 = xm >= cd if s1.sum() >= 12: f1 = ols(xm[s1], ym[s1]) trunc_lo.append((cu, f1)) left = "%11.4g %6d %8.4f [%.4f, %.4f]" % (10.0 ** cu, f1["n"], f1["b"], f1["lo"], f1["hi"]) else: left = "%11s %6s %8s %18s" % ("-", "-", "-", "-") if s2.sum() >= 12: f2 = ols(xm[s2], ym[s2]) trunc_hi.append((cd, f2)) right = "%11.4g %6d %8.4f [%.4f, %.4f]" % (10.0 ** cd, f2["n"], f2["b"], f2["lo"], f2["hi"]) else: right = "%11s %6s %8s %18s" % ("-", "-", "-", "-") print("%s %4s %s" % (left, "|", right)) print("-" * 108) print("Dropping everything above 100 g gives b = %.4f. Dropping everything" % trunc_lo[0][1]["b"]) print("below 316 g gives b = %.4f. Same database, same estimator, same day." % trunc_hi[-1][1]["b"]) print("The two intervals [%.4f, %.4f] and [%.4f, %.4f] do not overlap." % (trunc_lo[0][1]["lo"], trunc_lo[0][1]["hi"], trunc_hi[-1][1]["lo"], trunc_hi[-1][1]["hi"])) # ============================================================================ # SECTION 9. HOW WIDE A MASS RANGE DO YOU ACTUALLY NEED? # ============================================================================ head("SECTION 9. THE SENSITIVITY CURVE: RANGE AGAINST RESOLVING POWER") print("To tell 3/4 from 2/3 you need a 95%% interval narrower than %.4f, so a" % (THREE_QUARTERS - TWO_THIRDS)) print("half-width below %.4f. How much mass range buys that?" % ((THREE_QUARTERS - TWO_THIRDS) / 2)) print() print("For each window width W we draw %d random windows of that width from" % 400) print("the mammal mass axis (seeded), fit OLS inside each, and report the") print("spread of the fitted exponent and the mean interval half-width.") print() rng2 = np.random.default_rng(SEED + 1) print("%7s %8s %9s %9s %9s %9s %10s %10s %10s" % ("W dec", "windows", "mean n", "mean b", "sd of b", "min b", "max b", "mean CIhw", "P(resolve)")) print("-" * 96) sens = [] for W in [0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0]: if W > (xm.max() - xm.min()): continue bs, hws, ns, res = [], [], [], [] tries = 0 while len(bs) < 400 and tries < 20000: tries += 1 c = rng2.uniform(xm.min(), xm.max() - W) s = (xm >= c) & (xm < c + W) if s.sum() < 12: continue if xm[s].max() - xm[s].min() < 0.5 * W: continue f_ = ols(xm[s], ym[s]) bs.append(f_["b"]) hws.append(0.5 * (f_["hi"] - f_["lo"])) ns.append(f_["n"]) res.append(1.0 if 0.5 * (f_["hi"] - f_["lo"]) < (THREE_QUARTERS - TWO_THIRDS) / 2 else 0.0) if not bs: continue bs = np.array(bs) hws = np.array(hws) sens.append((W, bs, hws, float(np.mean(ns)), float(np.mean(res)))) print("%7.2f %8d %9.1f %9.4f %9.4f %9.4f %10.4f %10.4f %10.3f" % (W, bs.size, np.mean(ns), bs.mean(), bs.std(ddof=1), bs.min(), bs.max(), hws.mean(), np.mean(res))) print("-" * 96) print() print("Read the 'sd of b' column as the answer to the whole ninety-year") print("argument's statistical half. Below about two decades of mass the") print("sampling scatter in the exponent is itself comparable to the quantity") print("being argued over, so a study on a narrow mass range cannot settle it") print("no matter how carefully the respirometry was done.") # ============================================================================ # SECTION 10. IS IT EVEN A STRAIGHT LINE? # ============================================================================ head("SECTION 10. CURVATURE") print("Kolokotrones et al. (2010) reported that mammalian metabolic scaling is") print("convex on log-log axes, so that no single exponent describes it. Test:") print("add a quadratic term and see whether it earns its degree of freedom.") print() def quad_fit(xs, ys): xc = xs - xs.mean() A = np.vstack([np.ones(xs.size), xc, xc * xc]).T coef, *_ = np.linalg.lstsq(A, ys, rcond=None) resid = ys - A @ coef dfq = xs.size - 3 s2q = float(resid @ resid) / dfq cov = s2q * np.linalg.inv(A.T @ A) se = np.sqrt(np.diag(cov)) return coef, se, dfq, float(resid @ resid), float(xs.mean()) for gname in ["Mammalia", "Aves", "Endotherms (M+A)"]: s = fits[gname]["sel"] xs, ys = X[s], Y[s] coef, se, dfq, sseq, xbar = quad_fit(xs, ys) lin = ols(xs, ys) Fq = ((lin["SSE"] - sseq) / 1.0) / (sseq / dfq) print("%s (n = %d, centred at log10 M = %.3f):" % (gname, xs.size, xbar)) print(" linear term %.4f +/- %.4f" % (coef[1], se[1])) print(" quadratic term %.5f +/- %.5f t = %.2f p = %.3g" % (coef[2], se[2], coef[2] / se[2], t_sf2(coef[2] / se[2], dfq))) print(" F(1, %d) for adding curvature = %.2f, p = %.3g" % (dfq, Fq, f_sf(Fq, 1, dfq))) print(" local slope d logB/d logM at 10 g : %.4f" % (coef[1] + 2 * coef[2] * (1.0 - xbar))) print(" local slope at 1 kg : %.4f" % (coef[1] + 2 * coef[2] * (3.0 - xbar))) print(" local slope at 100 kg : %.4f" % (coef[1] + 2 * coef[2] * (5.0 - xbar))) print() print("The local-slope columns are the point. If curvature is real then") print("asking whether 'the' exponent is 3/4 or 2/3 is asking a question the") print("data decline to answer, because the slope depends on where you stand.") # ============================================================================ # SECTION 11. COMPARISON TO PUBLISHED FITS # ============================================================================ head("SECTION 11. OUR NUMBERS BESIDE PUBLISHED ONES") mo = fits["Mammalia"]["ols"] print("Club mammal OLS exponent: %.4f, SE %.4f, 95%% CI [%.4f, %.4f], n = %d" % (mo["b"], mo["se_b"], mo["lo"], mo["hi"], mo["n"])) print() print("%-46s %8s %10s %12s" % ("published value", "b", "our b - theirs", "in our SEs")) print("-" * 80) published = [ ("Kleiber (1932/1947), the canonical 3/4", 0.75), ("Rubner surface law, 2/3", TWO_THIRDS), ("White & Seymour (2003), 619 mammals, OLS", 0.69), ("White & Seymour (2003), refined set n=469", 0.68), ("Savage et al. (2004), binned mammal BMR", 0.737), ("Hayssen & Lacy (1985), mammals", 0.70), ("Clarke, Rothery & Isaac (2010), mammals", 0.70), ("Capellini et al. (2010), phylogenetic", 0.75), ] for label, val in published: print("%-46s %8.4f %10.4f %12.2f" % (label, val, mo["b"] - val, (mo["b"] - val) / mo["se_b"])) print("-" * 80) print("Those 'in our SEs' figures assume our SE is the right one. It is not:") print("assumption A3 fails, so treat anything under about three as agreement.") # ============================================================================ # SECTION 12. THE ANSWER # ============================================================================ head("SECTION 12. SUMMARY") mm_o = fits["Mammalia"]["ols"] mm_r = fits["Mammalia"]["rma"] av_o = fits["Aves"]["ols"] al_o = fits["All four classes"]["ols"] al_r = fits["All four classes"]["rma"] print("1. Mammals, n = %d, %.2f decades of mass." % (mm_o["n"], fits["Mammalia"]["dec"])) print(" OLS b = %.4f, 95%% CI [%.4f, %.4f]. %s" % (mm_o["b"], mm_o["lo"], mm_o["hi"], verdict(mm_o["lo"], mm_o["hi"]))) print(" RMA b = %.4f, 95%% CI [%.4f, %.4f]. %s" % (mm_r["b"], mm_r["lo"], mm_r["hi"], verdict(mm_r["lo"], mm_r["hi"]))) _mid = 0.5 * (THREE_QUARTERS + TWO_THIRDS) print(" The exact midpoint of the two disputed values is %.6f." % _mid) print(" Our mammal OLS exponent sits %.6f from it, which is %.3f of one" % (mm_o["b"] - _mid, (mm_o["b"] - _mid) / mm_o["se_b"])) print(" standard error. We did not arrange that and we cannot explain it.") print() print("2. Birds, n = %d. OLS b = %.4f, 95%% CI [%.4f, %.4f]. %s" % (av_o["n"], av_o["b"], av_o["lo"], av_o["hi"], verdict(av_o["lo"], av_o["hi"]))) print() print("3. All four classes, n = %d, %.2f decades." % (al_o["n"], fits["All four classes"]["dec"])) print(" OLS b = %.4f [%.4f, %.4f] %s" % (al_o["b"], al_o["lo"], al_o["hi"], verdict(al_o["lo"], al_o["hi"]))) print(" RMA b = %.4f [%.4f, %.4f] %s" % (al_r["b"], al_r["lo"], al_r["hi"], verdict(al_r["lo"], al_r["hi"]))) print(" The two methods put the answer on OPPOSITE SIDES of 3/4 on the same") print(" %d species. The method is not a detail." % al_o["n"]) print() print("4. Slopes differ between classes: F(%d, %d) = %.2f, p = %.3g." % (A4["df1"], A4["df2"], A4["F"], A4["p"])) print(" Mammal minus bird slope = %.4f, which is %.0f%% of the 3/4-to-2/3 gap." % (mm_o["b"] - av_o["b"], 100.0 * (mm_o["b"] - av_o["b"]) / (THREE_QUARTERS - TWO_THIRDS))) print() print("5. Restricting the mass range moves the answer by %.4f across %d" % (sb.max() - sb.min(), len(slide))) print(" two-decade windows, which is %.1f times the quantity in dispute." % ((sb.max() - sb.min()) / (THREE_QUARTERS - TWO_THIRDS))) print() print("6. Neither 3/4 nor 2/3 survives as a universal constant in this") print(" compilation. Both survive as a description of some subset of it.") print() print("Runtime: %.2f s. numpy %s, Python %s" % (time.time() - T0, np.__version__, sys.version.split()[0])) print("Seed %d (bootstrap and window sampling only)." % SEED) # ============================================================================ # SECTION 13. FIGURE DATA # ============================================================================ head("SECTION 13. NUMBERS BEHIND THE FIGURES") print("[FIG2] forest plot: group, method, b, lo, hi") for name, _ in GROUPS: o, m = fits[name]["ols"], fits[name]["rma"] print("FIG2 %-18s OLS %.5f %.5f %.5f" % (name, o["b"], o["lo"], o["hi"])) print("FIG2 %-18s RMA %.5f %.5f %.5f" % (name, m["b"], m["lo"], m["hi"])) print() print("[FIG3] sliding windows: centre, n, b_ols, lo, hi, b_rma") for c0, c1, f_, m_ in slide: print("FIG3 %.3f %d %.5f %.5f %.5f %.5f" % (0.5 * (c0 + c1), f_["n"], f_["b"], f_["lo"], f_["hi"], m_["b"])) print() print("[FIG4] sensitivity: W, sd_b, mean_halfwidth, p_resolve, mean_n") for W, bs, hws, mn, pr in sens: print("FIG4 %.2f %.5f %.5f %.4f %.1f" % (W, bs.std(ddof=1), hws.mean(), pr, mn)) print() print("[FIG5] curvature, mammals: local slope at each decade of mass") s = fits["Mammalia"]["sel"] coef, se, dfq, sseq, xbar = quad_fit(X[s], Y[s]) xcq = np.vstack([np.ones(X[s].size), X[s] - xbar, (X[s] - xbar) ** 2]).T covq = (sseq / dfq) * np.linalg.inv(xcq.T @ xcq) tq = tcrit(dfq) for xv in np.arange(0.5, 6.51, 0.25): g = np.array([0.0, 1.0, 2.0 * (xv - xbar)]) sl = coef[1] + 2 * coef[2] * (xv - xbar) sg = math.sqrt(float(g @ covq @ g)) print("FIG5 %.2f %.5f %.5f %.5f" % (xv, sl, sl - tq * sg, sl + tq * sg)) print() print("[FIG1] scatter: the full point cloud, class, log10 mass, log10 BMR") print("FIG1N %d" % len(recs)) for i in range(len(recs)): print("FIG1 %s %.4f %.4f" % (CLS[i][:4], X[i], Y[i]))