"""Part 2: do Phanerozoic innovations with a physiological oxygen requirement originate after rises in oxygen (or other environmental steps) more often than expected if origin times were independent of the curve? For each environmental reconstruction c(age) on a 1-Myr grid we tabulate, for every age, rise_W : c(age) - c(age + W) change over the preceding W Myr (W = 20 by default) level_q : the quantile of c(age) among all ages of the domain record : 1 if c(age) is within 5 per cent of its running maximum over all older ages of the record (a 'gate-opening epoch': a threshold is being crossed for the first time) minlag : age_first(c >= c(age)) - age, the time since the curve FIRST reached the level it has at `age`. If an innovation needs oxygen above some threshold and appears at `age`, its threshold cannot exceed c(age), so the gate opened no later than age + minlag: minlag is the smallest gate-to- innovation lag compatible with ANY threshold, for that reconstruction. Test statistics are the means of these over the innovations (date uncertainty averaged over draws of first appearance plus range extension). Null distributions: shift : all origin times shifted together by s Myr, wrapping around the domain (every integer s). Keeps the clustering of the origins and the autocorrelation of the curve. jitter : each origin time moved independently by up to +/- 50 Myr (tests association beyond the broad trend) pbdb : origin times redrawn from the distribution of first appearances of all animal and plant genera in the Paleobiology Database (innovations compared with ordinary originations, which share the pull of the recent and the unevenness of the rock record). Effect sizes are reported in natural units (per cent O2 per 20 Myr, rate ratios) with their null intervals. Added in revision (referee): a run with literature dates for all innovations; a count of all 140 oxygen tests in the direction gating predicts and in the opposite direction (part2_multiplicity.json); the look-back-window sensitivity on both the common and the largest innovation set, with n; and power for 14 and 9 independent origins, because the real origins are clustered and the 27-origin power is an upper bound. """ import json import os import numpy as np import pandas as pd import pbdb from config import DATA, RAW, SEED from curves import build as build_curves W = 20 REC_TOL = 0.05 R_DATES = 300 def curve_stats(c, w=W): """c: array indexed by age (Ma), length L+1. Returns dict of arrays on ages 0..L-w.""" L = len(c) - 1 n = L - w + 1 rise = c[:n] - c[w:w + n] order = np.argsort(np.argsort(c[:n])) level_q = (order + 0.5) / n # running maximum over all older ages (from the oldest end of the record towards the present) runmax = np.maximum.accumulate(c[::-1])[::-1] record = (c >= (1 - REC_TOL) * runmax)[:n].astype(float) # first (oldest) age at which the curve was >= its value at each age minlag = np.empty(n) for a in range(n): older = np.nonzero(c[a:] >= c[a])[0] minlag[a] = older[-1] # offset of the oldest age with c >= c[a]; 0 if never before return dict(rise=rise, level_q=level_q, record=record, minlag=minlag, rising=(rise > 0).astype(float), n=n) def lookup(arr, ages): idx = np.clip(np.rint(ages).astype(int), 0, len(arr) - 1) return arr[idx] def test_curve(c, origins, rng, w=W, null_times=None, label=""): """origins: (R, n) origin ages (Ma) for the innovations inside the domain. Returns list of result rows.""" st = curve_stats(c, w) n_dom = st["n"] O = np.clip(origins, 0, n_dom - 1) rows = [] shifts = np.arange(1, n_dom) def stat_fns(): yield "rise_%dMyr" % w, st["rise"], "mean change of the curve over the %d Myr before origin" % w yield "frac_rising", st["rising"], "fraction of origins at times when the curve had risen over %d Myr" % w yield "level_quantile", st["level_q"], "mean quantile of the curve at origin (0.5 = typical)" yield "frac_record", st["record"], "fraction of origins within 5%% of the running maximum" yield "minlag_le50", (st["minlag"] <= 50).astype(float), "fraction with minimum gate lag <= 50 Myr" yield "minlag_median", st["minlag"], "median over innovations of the minimum gate lag (Myr)" for name, arr, desc in stat_fns(): agg = np.median if name == "minlag_median" else np.mean def T(ages): v = lookup(arr, ages) # (R, n) return float(np.mean(agg(v, axis=1))) obs = T(O) null_shift = np.array([T((O + s) % n_dom) for s in shifts]) jit = np.array([T(np.abs(((O + rng.uniform(-50, 50, size=(1, O.shape[1]))) + n_dom) % (2 * n_dom) - n_dom) .clip(0, n_dom - 1)) for _ in range(500)]) row = dict(curve=label, statistic=name, description=desc, n_innovations=O.shape[1], domain_Ma=n_dom - 1, observed=obs, time_average=float(agg(arr)) if name != "minlag_median" else float(np.median(arr)), null_shift_mean=float(null_shift.mean()), null_shift_sd=float(null_shift.std()), null_shift_p05=float(np.quantile(null_shift, 0.05)), null_shift_p95=float(np.quantile(null_shift, 0.95)), p_shift_ge=float((np.sum(null_shift >= obs) + 1) / (len(null_shift) + 1)), p_shift_le=float((np.sum(null_shift <= obs) + 1) / (len(null_shift) + 1)), null_jitter_mean=float(jit.mean()), null_jitter_sd=float(jit.std()), p_jitter_ge=float((np.sum(jit >= obs) + 1) / (len(jit) + 1))) if null_times is not None: nb = np.array([T(rng.choice(null_times, size=(1, O.shape[1]))) for _ in range(1000)]) row.update(null_pbdb_mean=float(nb.mean()), null_pbdb_sd=float(nb.std()), p_pbdb_ge=float((np.sum(nb >= obs) + 1) / (len(nb) + 1))) row["z_shift"] = (obs - row["null_shift_mean"]) / row["null_shift_sd"] if row["null_shift_sd"] > 0 else np.nan rows.append(row) # rate ratio: origins per Myr when the curve is rising versus falling fr = st["rising"].mean() def RR(ages): k = lookup(st["rising"], ages).mean(axis=1) # fraction of innovations in rising times k = np.clip(k, 1e-3, 1 - 1e-3) return float(np.exp(np.mean(np.log((k / fr) / ((1 - k) / (1 - fr)))))) obs = RR(O) null = np.array([RR((O + s) % n_dom) for s in shifts]) rows.append(dict(curve=label, statistic="rate_ratio_rising_vs_falling", description="origins per Myr while the curve is rising divided by origins per Myr while falling", n_innovations=O.shape[1], domain_Ma=n_dom - 1, observed=obs, time_average=float(fr), null_shift_mean=float(np.exp(np.mean(np.log(null)))), null_shift_sd=float(np.std(np.log(null))), null_shift_p05=float(np.quantile(null, 0.05)), null_shift_p95=float(np.quantile(null, 0.95)), p_shift_ge=float((np.sum(null >= obs) + 1) / (len(null) + 1)), p_shift_le=float((np.sum(null <= obs) + 1) / (len(null) + 1)))) return rows, st def pbdb_origination_times(rng, n=200000): """Ages drawn from the first appearances of all animal and plant genera (stage resolution, uniform within stage), from the PBDB diversity tables (X_Ft + X_FL = genera first appearing in the bin).""" out = [] for f in ("div_animalia_stage.csv", "div_plantae_stage.csv"): d = pbdb.read_pbdb_csv(os.path.join(RAW, "pbdb", f)) wgt = (d["X_Ft"] + d["X_FL"]).to_numpy(float) p = wgt / wgt.sum() k = rng.choice(len(d), size=n // 2, p=p) out.append(rng.uniform(d["min_ma"].to_numpy()[k], d["max_ma"].to_numpy()[k])) return np.concatenate(out) def _shift_p(arr, O, n_dom): """One-sided circular-shift p-value (large values = association) for one set of origin ages O (n,).""" s = np.arange(1, n_dom) idx = (np.rint(O).astype(int)[None, :] + s[:, None]) % n_dom null = arr[idx].mean(axis=1) obs = arr[np.clip(np.rint(O).astype(int), 0, n_dom - 1)].mean() return (np.sum(null >= obs) + 1) / (len(null) + 1) def simulate_gated(rng, c, n, frac, tau, n_dom): """n origin ages on the domain 0..n_dom-1. A fraction `frac` are gated: an oxygen threshold is drawn uniformly between the level at the old end of the record and the maximum of the curve, the gate opens the first time the curve reaches it, and the origin follows after an exponential lag of mean tau (Myr). The rest are uniform over the domain.""" out = np.empty(n) c = c[:n_dom] lo, hi = c[-1], c.max() for i in range(n): if rng.uniform() < frac: while True: th = rng.uniform(lo, hi) gate = np.nonzero(c >= th)[0][-1] # oldest age at which the curve is >= threshold a = gate - rng.exponential(tau) if a >= 0: out[i] = a break else: out[i] = rng.uniform(0, n_dom - 1) return out def power_analysis(rng, curves, n=27, nsim=300): """Power of the circular-shift tests (one-sided, alpha = 0.05) to detect threshold gating, when the curve used in the test is the true one and when it is a different reconstruction.""" g = curves["O2 GEOCARBSULFOR median"] k = curves["O2 COPSE standard"][:len(g)] stg, stk = curve_stats(g), curve_stats(k) n_dom = stg["n"] rows = [] for true_name, c_true in (("GEOCARBSULFOR median", g), ("COPSE standard", k)): for test_name, st in (("GEOCARBSULFOR median", stg), ("COPSE standard", stk)): for frac in (0.0, 0.25, 0.5, 0.75, 1.0): for tau in (10.0, 30.0): if frac == 0.0 and tau != 10.0: continue hits = {"rise_20Myr": 0, "frac_record": 0, "minlag_le50": 0} for _ in range(nsim): O = simulate_gated(rng, c_true, n, frac, tau, n_dom) for name, arr in (("rise_20Myr", st["rise"]), ("frac_record", st["record"]), ("minlag_le50", (st["minlag"] <= 50).astype(float))): hits[name] += _shift_p(arr, O, n_dom) < 0.05 rows.append(dict(true_curve=true_name, test_curve=test_name, gated_fraction=frac, mean_lag_Myr=tau, n_innovations=n, n_sim=nsim, **{"power_" + a: b / nsim for a, b in hits.items()})) df = pd.DataFrame(rows) df.to_csv(os.path.join(DATA, "part2_power.csv"), index=False) return df def effective_n(med, domain, gap=10.0): """How many separate episodes the origins inside a curve's domain form. Two counts: clusters of origins separated by more than `gap` Myr, and the number of distinct W-Myr bins that contain an origin.""" a = np.sort(med[med <= domain - W - 2]) clusters = 1 + int(np.sum(np.diff(a) > gap)) if len(a) else 0 bins = int(len(np.unique(np.floor(a / W)))) if len(a) else 0 return {"n": int(len(a)), "n_clusters_gap_gt_%dMyr" % gap: clusters, "n_occupied_%dMyr_bins" % W: bins} def power_vs_n(rng, curves, ns=(27, 14, 9), nsim=300): """The quoted power assumes n independent origins. The real origins are clustered, so the power for fewer, independent origins bounds it from below. Half the origins gated, test curve = true curve.""" rows = [] g = curves["O2 GEOCARBSULFOR median"] k = curves["O2 COPSE standard"][:len(g)] for name, c in (("GEOCARBSULFOR median", g), ("COPSE standard", k)): st = curve_stats(c) for tau in (10.0, 30.0): for n in ns: hits = 0 for _ in range(nsim): O = simulate_gated(rng, c, n, 0.5, tau, st["n"]) hits += _shift_p(st["rise"], O, st["n"]) < 0.05 rows.append(dict(curve=name, gated_fraction=0.5, mean_lag_Myr=tau, n_independent_origins=n, n_sim=nsim, power_rise_20Myr=hits / nsim)) df = pd.DataFrame(rows) df.to_csv(os.path.join(DATA, "part2_power_vs_n.csv"), index=False) return df def curve_agreement(curves): """Do the oxygen reconstructions agree on WHEN oxygen rose? Correlations of level and of 20-Myr change.""" names = [n for n in curves if n.startswith("O2")] rows = [] for i, a in enumerate(names): for b in names[i + 1:]: L = min(len(curves[a]), len(curves[b])) - 1 ca, cb = curves[a][:L + 1], curves[b][:L + 1] ra, rb = ca[:L + 1 - W] - ca[W:], cb[:L + 1 - W] - cb[W:] rows.append(dict(curve_a=a, curve_b=b, overlap_Ma=L, corr_level=float(np.corrcoef(ca, cb)[0, 1]), corr_rise_20Myr=float(np.corrcoef(ra, rb)[0, 1]), frac_same_sign_rise=float(np.mean(np.sign(ra) == np.sign(rb))))) df = pd.DataFrame(rows) df.to_csv(os.path.join(DATA, "part2_curve_agreement.csv"), index=False) return df def run_all(): rng = np.random.default_rng(SEED + 10) cv = build_curves() tab, draws = pbdb.build() ids = list(tab["id"]) fad = np.column_stack([draws[i][0][:R_DATES] for i in ids]) ext = np.column_stack([draws[i][1][:R_DATES] for i in ids]) origin = fad + ext # (R, n_innov), Ma alt = tab["alt_old"].to_numpy(float) origin_alt = np.where(np.isnan(alt)[None, :], origin, np.maximum(origin, alt[None, :])) # every innovation dated from the literature bracket (no PBDB dates), same range extensions origin_lit = np.column_stack([draws[i][2][:R_DATES] for i in ids]) + ext med = np.median(origin, axis=0) cls = tab["o2_class"].to_numpy() null_times = pbdb_origination_times(rng) t = cv["table"] gc = cv["gcso"] curves = { "O2 GEOCARBSULFOR median": t["gcso_median"].to_numpy()[:551], "O2 COPSE standard": t["copse_standard"].to_numpy()[:801], "O2 COPSE no feedbacks": t["copse_no_feedbacks"].to_numpy()[:801], "O2 COPSE both feedbacks": t["copse_both_feedbacks"].to_numpy()[:801], "O2 charcoal midpoint": t["charcoal_mid"].to_numpy()[:401], "T PhanDA GMST": t["phanda_gmst"].to_numpy()[:486], "T Scotese GAT": t["scotese_gat"].to_numpy()[:541], "CO2 PhanDA (log)": np.log10(t["phanda_co2"].to_numpy()[:486]), } subsets = { "all": np.ones(len(ids), bool), "high-O2-demand": cls == "high", "younger than 500 Ma (no Ediacaran-Cambrian cluster)": med < 500, "PBDB-dated only": tab["fad_basis"].to_numpy() == "PBDB", } rows, per_innov = [], {} for cname, c in curves.items(): L = len(c) - 1 for sname, mask in subsets.items(): inside = mask & (med <= L - W - 2) if inside.sum() < 5: continue for dates_label, Omat in (("fossil first appearance + range extension", origin), ("older contested claims adopted", origin_alt), ("literature dates for all innovations", origin_lit)): if not dates_label.startswith("fossil") and sname != "all": continue ins = inside if dates_label.startswith("fossil") else mask & (np.median(Omat, axis=0) <= L - W - 2) r, st = test_curve(c, Omat[:, ins], rng, null_times=null_times[null_times <= L - W], label=cname) for x in r: x.update(subset=sname, dates=dates_label) rows += r st = curve_stats(c) per_innov[cname] = {i: (float(np.median(lookup(st["minlag"], np.clip(origin[:, k], 0, st["n"] - 1)))), float(np.mean(lookup(st["rise"], np.clip(origin[:, k], 0, st["n"] - 1)))), float(np.mean(lookup(st["level_q"], np.clip(origin[:, k], 0, st["n"] - 1))))) if med[k] <= L - W - 2 else (np.nan, np.nan, np.nan) for k, i in enumerate(ids)} res = pd.DataFrame(rows) # one-sided p-value in the direction gating predicts (origins after rises, at high levels, at record # highs, with short minimum lags) and in the opposite direction low_is_gating = res["statistic"] == "minlag_median" res["p_gating_direction"] = np.where(low_is_gating, res["p_shift_le"], res["p_shift_ge"]) res["p_opposite_direction"] = np.where(low_is_gating, res["p_shift_ge"], res["p_shift_le"]) res.to_csv(os.path.join(DATA, "part2_tests.csv"), index=False) ox = res[res["curve"].str.startswith("O2") & res["dates"].str.startswith("fossil")] full = ox[ox["subset"] == "all"] mult = dict( statistics=sorted(ox["statistic"].unique()), curves=sorted(ox["curve"].unique()), subsets=sorted(ox["subset"].unique()), n_tests=int(len(ox)), n_p_lt_05_gating_direction=int((ox["p_gating_direction"] < 0.05).sum()), n_p_lt_05_opposite_direction=int((ox["p_opposite_direction"] < 0.05).sum()), expected_by_chance_each_direction=round(0.05 * len(ox), 1), by_statistic={st: dict(gating=int((g["p_gating_direction"] < 0.05).sum()), opposite=int((g["p_opposite_direction"] < 0.05).sum()), n=int(len(g))) for st, g in ox.groupby("statistic")}, full_set=dict(n_tests=int(len(full)), n_p_lt_05_gating_direction=int((full["p_gating_direction"] < 0.05).sum()), n_p_lt_05_opposite_direction=int((full["p_opposite_direction"] < 0.05).sum()), min_p_gating_direction=float(full["p_gating_direction"].min()), min_p_opposite_direction=float(full["p_opposite_direction"].min())), significant_gating_direction=ox[ox["p_gating_direction"] < 0.05][ ["curve", "subset", "statistic", "p_gating_direction"]].round(4).to_dict(orient="records"), significant_opposite_direction=ox[ox["p_opposite_direction"] < 0.05][ ["curve", "subset", "statistic", "p_opposite_direction"]].round(4).to_dict(orient="records")) with open(os.path.join(DATA, "part2_multiplicity.json"), "w") as f: json.dump(mult, f, indent=1) # sensitivity to the oxygen reconstruction: every COPSE variant and 200 GEOCARBSULFOR realisations sens = [] inside_all = med <= 550 - W - 2 for k, v in cv["copse"].items(): r, _ = test_curve(v[:801], origin[:100][:, med <= 800 - W - 2], rng, label="COPSE " + k) for x in r: if x["statistic"] in ("rise_%dMyr" % W, "rate_ratio_rising_vs_falling", "minlag_le50", "frac_record"): sens.append(dict(family="COPSE", member=k, statistic=x["statistic"], observed=x["observed"], null_mean=x["null_shift_mean"], z=x.get("z_shift", np.nan), p_ge=x["p_shift_ge"])) pick = rng.choice(gc.shape[1], size=200, replace=False) for j in pick: st = curve_stats(gc[:, j]) O = np.clip(origin[:60][:, inside_all], 0, st["n"] - 1) for name, arr in (("rise_%dMyr" % W, st["rise"]), ("minlag_le50", (st["minlag"] <= 50).astype(float)), ("frac_record", st["record"])): obs = float(lookup(arr, O).mean()) null = np.array([lookup(arr, (O + s) % st["n"]).mean() for s in range(1, st["n"], 3)]) sens.append(dict(family="GEOCARBSULFOR MC", member=int(j), statistic=name, observed=obs, null_mean=float(null.mean()), z=float((obs - null.mean()) / null.std()), p_ge=float((np.sum(null >= obs) + 1) / (len(null) + 1)))) sens = pd.DataFrame(sens) sens.to_csv(os.path.join(DATA, "part2_reconstruction_sensitivity.csv"), index=False) # per-innovation table pi = tab[["id", "name", "o2_class", "fad_old", "fad_young", "fad_basis", "ext50", "ext95", "origin50", "alt_old", "verified", "lit_age", "pbdb_fad_taxon", "pbdb_younger_than_literature", "pbdb_minus_lit_Myr"]].copy() for cname in ("O2 GEOCARBSULFOR median", "O2 COPSE standard", "O2 charcoal midpoint"): short = cname.split()[1].lower()[:6] pi["minlag_" + short] = [per_innov[cname][i][0] for i in ids] pi["rise20_" + short] = [per_innov[cname][i][1] for i in ids] pi["levelq_" + short] = [per_innov[cname][i][2] for i in ids] pi.to_csv(os.path.join(DATA, "part2_per_innovation.csv"), index=False) # how many innovations could have been promptly oxygen-gated under every reconstruction lagcols = [c for c in pi.columns if c.startswith("minlag_")] have = pi[lagcols].notna().any(axis=1) summ = dict( n_innovations=int(len(pi)), n_with_any_curve=int(have.sum()), n_minlag_gt50_all_curves=int(((pi[lagcols] > 50) | pi[lagcols].isna()).all(axis=1)[have].sum()), n_minlag_gt100_all_curves=int(((pi[lagcols] > 100) | pi[lagcols].isna()).all(axis=1)[have].sum()), n_minlag_le50_any_curve=int((pi[lagcols] <= 50).any(axis=1).sum()), n_minlag_le50_all_curves=int(((pi[lagcols] <= 50) | pi[lagcols].isna()).all(axis=1)[have].sum()), ids_minlag_gt100_all=list(pi["id"][have & ((pi[lagcols] > 100) | pi[lagcols].isna()).all(axis=1)]), ids_minlag_le50_all=list(pi["id"][have & ((pi[lagcols] <= 50) | pi[lagcols].isna()).all(axis=1)]), ) with open(os.path.join(DATA, "part2_summary.json"), "w") as f: json.dump(summ, f, indent=1) # sensitivity of the primary statistic to the length of the look-back window wrows = [] for cname in ("O2 GEOCARBSULFOR median", "O2 COPSE standard", "O2 charcoal midpoint"): c = curves[cname] L = len(c) - 1 for w in (10, 20, 40): for set_label, ins in (("common set (inside the 40-Myr domain)", med <= L - 40 - 2), ("largest set for this window", med <= L - w - 2)): r, _ = test_curve(c, origin[:, ins], rng, w=w, label=cname) x = r[0] wrows.append(dict(curve=cname, window_Myr=w, innovation_set=set_label, n_innovations=x["n_innovations"], observed=x["observed"], null_shift_mean=x["null_shift_mean"], null_shift_p05=x["null_shift_p05"], null_shift_p95=x["null_shift_p95"], p_shift_ge=x["p_shift_ge"], z_shift=x["z_shift"])) wsens = pd.DataFrame(wrows) wsens.to_csv(os.path.join(DATA, "part2_window_sensitivity.csv"), index=False) # null distribution of the primary statistic, saved for the figure nrows = [] for cname in ("O2 GEOCARBSULFOR median", "O2 COPSE standard", "O2 charcoal midpoint"): c = curves[cname] st = curve_stats(c) ins = med <= len(c) - 1 - W - 2 O = np.clip(origin[:, ins], 0, st["n"] - 1) nrows.append(dict(curve=cname, shift_Myr=0, mean_rise=float(lookup(st["rise"], O).mean()))) for sft in range(1, st["n"]): nrows.append(dict(curve=cname, shift_Myr=sft, mean_rise=float(lookup(st["rise"], (O + sft) % st["n"]).mean()))) pd.DataFrame(nrows).to_csv(os.path.join(DATA, "part2_null_primary.csv"), index=False) power = power_analysis(rng, curves) agree = curve_agreement(curves) eff = {"GEOCARBSULFOR": effective_n(med, 550), "COPSE": effective_n(med, 800), "charcoal": effective_n(med, 400)} power_n = power_vs_n(rng, curves) summ["effective_n"] = eff with open(os.path.join(DATA, "part2_summary.json"), "w") as f: json.dump(summ, f, indent=1) return dict(tests=res, sens=sens, per_innovation=pi, summary=summ, origin=origin, ids=ids, tab=tab, curves=curves, power=power, power_n=power_n, agreement=agree, null_times=null_times, window=wsens, med=med, multiplicity=mult) if __name__ == "__main__": import time t0 = time.time() pd.set_option("display.width", 260); pd.set_option("display.max_colwidth", 60) r = run_all() q = r["tests"] print(q[(q["dates"].str.startswith("fossil"))][["curve", "subset", "statistic", "n_innovations", "observed", "time_average", "null_shift_mean", "null_shift_p05", "null_shift_p95", "p_shift_ge", "p_shift_le", "p_jitter_ge", "p_pbdb_ge"]].round(3).to_string()) print(r["per_innovation"].round(2).to_string()) print(json.dumps(r["summary"], indent=1)) s = r["sens"] print(s.groupby(["family", "statistic"]).agg(obs=("observed", "median"), null=("null_mean", "median"), z_med=("z", "median"), z_lo=("z", "min"), z_hi=("z", "max"), frac_p_lt_05=("p_ge", lambda x: float(np.mean(x < 0.05)))).round(3)) print(r["window"].round(3).to_string()) print(r["power"].round(2).to_string()) print(r["agreement"].round(2).to_string()) print("seconds", time.time() - t0)