import sys import os #!/usr/bin/env python3 """K16a meta_check2.py - independent recomputation of the 2026 update (update2026.json) and of the wcrf2007 Fig 4.3.3 runs in all.json. Written from scratch by the K16a controller. Does NOT read or import scripts/update2026.py, pool.py, mvmeta.py, multiverse.py or catsets.py. Inputs: kontroll/reanalyser/update2026.json, all.json (published results + per-study rows) database/estimates.csv (raw rows, to rebuild rows that the json lacks) Methods (own implementation) y = ln RR, se = (ln hi - ln lo)/(2*1.959964) FE: inverse variance. Q, df=k-1, I2 = max(0,(Q-df)/Q). DL: tau2 = max(0,(Q-df)/(sum w - sum w^2/sum w)); RE weights 1/(v+tau2). HKSJ: q = sum w*(y-mu)^2/(k-1); se = sqrt(q/sum w*); t(k-1) quantile (no max(1,q) truncation; the truncated variant is also reported). REML: Fisher-scoring/fixed-point iteration on tau2. Egger: OLS of y/se on 1/se; intercept t-test, df = k-2. Greenland-Longnecker (GL): pseudo-counts fitted to the reported RRs with the observed total cases and category denominators; cov(L_i,L_j) = 1/A0 (+1/B0 for cumulative-incidence / case-control); GLS slope through the origin on dose - ref dose. Rules used when denominators/cases are missing are stated in RAPPORT.md (they are this checker's assumptions, not the pipeline's). Usage: python3 meta_check2.py (writes nothing; prints all comparisons; RAPPORT.md was built from the output) """ import csv import itertools import json import math import statistics import numpy as np from scipy import optimize, stats BASE = sys.argv[1] if len(sys.argv) > 1 else os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) U = json.load(open(BASE + "/kontroll/reanalyser/update2026.json")) ALL = json.load(open(BASE + "/kontroll/reanalyser/all.json")) EST = {r["estimate_id"]: r for r in csv.DictReader(open(BASE + "/database/estimates.csv"))} Z = 1.959964 TOL4 = 5e-5 # ----------------------------------------------------------------------------- core pooling def yv(rr, lo, hi): se = (math.log(hi) - math.log(lo)) / (2 * Z) return math.log(rr), se * se def dl(rows): """rows: list of (rr, lo, hi) or (y, v, 'yv').""" ys, vs = [], [] for r in rows: if len(r) == 3 and r[2] == "yv": raise ValueError y, v = yv(*r[:3]) ys.append(y) vs.append(v) return dl_yv(ys, vs) def dl_yv(ys, vs): ys, vs = np.array(ys, float), np.array(vs, float) k = len(ys) w = 1 / vs sw = w.sum() mfe = (w * ys).sum() / sw Q = float((w * (ys - mfe) ** 2).sum()) df = k - 1 C = sw - (w ** 2).sum() / sw tau2 = max(0.0, (Q - df) / C) if k > 1 else 0.0 ws = 1 / (vs + tau2) mre = (ws * ys).sum() / ws.sum() sere = math.sqrt(1 / ws.sum()) I2 = max(0.0, (Q - df) / Q) * 100 if Q > 0 else 0.0 out = dict(k=k, mu=mre, se=sere, re=[math.exp(mre), math.exp(mre - Z * sere), math.exp(mre + Z * sere)], fe=[math.exp(mfe), math.exp(mfe - Z / math.sqrt(sw)), math.exp(mfe + Z / math.sqrt(sw))], Q=Q, tau2=tau2, I2=I2, ys=ys, vs=vs, wpct=100 * ws / ws.sum()) if k > 1: q = float((ws * (ys - mre) ** 2).sum()) / df t = stats.t.ppf(0.975, df) se_h = math.sqrt(q / ws.sum()) out["hksj"] = [math.exp(mre), math.exp(mre - t * se_h), math.exp(mre + t * se_h)] se_h1 = math.sqrt(max(1.0, q) / ws.sum()) out["hksj_trunc"] = [math.exp(mre), math.exp(mre - t * se_h1), math.exp(mre + t * se_h1)] return out def reml(ys, vs): ys, vs = np.array(ys), np.array(vs) tau2 = dl_yv(ys, vs)["tau2"] for _ in range(1000): w = 1 / (vs + tau2) mu = (w * ys).sum() / w.sum() # REML fixed-point: tau2 = sum w^2((y-mu)^2 - v) / sum w^2 + 1/sum w new = max(0.0, (w ** 2 * ((ys - mu) ** 2 - vs)).sum() / (w ** 2).sum() + 1 / w.sum()) if abs(new - tau2) < 1e-12: tau2 = new break tau2 = new w = 1 / (vs + tau2) mu = (w * ys).sum() / w.sum() se = math.sqrt(1 / w.sum()) return [math.exp(mu), math.exp(mu - Z * se), math.exp(mu + Z * se), tau2] def egger(ys, vs): ys, se = np.array(ys), np.sqrt(np.array(vs)) Y, X = ys / se, 1 / se A = np.column_stack([np.ones_like(X), X]) beta, res, *_ = np.linalg.lstsq(A, Y, rcond=None) resid = Y - A @ beta df = len(Y) - 2 s2 = (resid ** 2).sum() / df cov = s2 * np.linalg.inv(A.T @ A) t = beta[0] / math.sqrt(cov[0, 0]) return float(beta[0]), float(2 * stats.t.sf(abs(t), df)) # ----------------------------------------------------------------------------- reporting LINES = [] def out(s=""): print(s) LINES.append(s) def f3(t): return "%.4f (%.4f-%.4f)" % tuple(t[:3]) def cmp(name, pub, rec, note=""): d = max(abs(a - b) for a, b in zip(pub[:3], rec[:3])) m2 = all(round(a + 1e-12, 2) == round(b + 1e-12, 2) for a, b in zip(pub[:3], rec[:3])) m4 = d < TOL4 out("| %s | %s | %s | %s | %s | %s |" % (name, f3(pub), f3(rec), "yes" if m2 else "NO", "yes" if m4 else "NO (max diff %.5f)" % d, note)) return m4 def cmp_scalar(name, pub, rec, note="", tol=TOL4): m2 = round(pub, 2) == round(rec, 2) m4 = abs(pub - rec) < tol out("| %s | %.4f | %.4f | %s | %s | %s |" % (name, pub, rec, "yes" if m2 else "NO", "yes" if m4 else "NO (diff %.5f)" % (rec - pub), note)) return m4 HDR = "| Item | Published | Recomputed | Match 2 dp | Match 4 dp | Note |\n|---|---|---|---|---|---|" # ----------------------------------------------------------------------------- data SET = {r["key"]: (r["rr"], r["lo"], r["hi"]) for r in U["set"]} KEYS = [r["key"] for r in U["set"]] SLUG = dict(zip(U["main"]["loo"].keys(), KEYS)) # loo keys are in set order MAIN = [SET[k] for k in KEYS] def e(eid): r = EST[eid] return float(r["point"]), float(r["ci_low"]), float(r["ci_high"]) def powr(t, p): return tuple(x ** p for x in t) def comb_fixed(rows): p = dl(rows) return tuple(p["fe"]) # ----------------------------------------------------------------------------- 0. main pool out("## 0. Main pool (k=20) - baseline") out(HDR) P = dl(MAIN) m = U["main"] cmp("DL random effects", m["re"], P["re"]) cmp("Fixed effect", m["fe"], P["fe"]) cmp_scalar("tau2", m["tau2"], P["tau2"]) cmp_scalar("Q", m["Q"], P["Q"]) cmp_scalar("I2 (%)", m["I2"], P["I2"]) for r, wp in zip(U["set"], P["wpct"]): if abs(r["w"] - wp) > 1e-3: out("| weight %s | %.4f | %.4f | | NO | |" % (r["key"], r["w"], wp)) out("| per-study RE weights (20) | json `w` | recomputed | | %s | max diff %.2e |" % ( "yes" if max(abs(r["w"] - wp) for r, wp in zip(U["set"], P["wpct"])) < 1e-6 else "NO", max(abs(r["w"] - wp) for r, wp in zip(U["set"], P["wpct"])))) out() # ----------------------------------------------------------------------------- 1. subgroups out("## 1. Subgroups (membership from json `studies`)") out(HDR) for name, g in U["subgroups"].items(): rows = [SET[s] for s in g["studies"]] p = dl(rows) ok = cmp(name + " (k=%d)" % g["k"], g["re"], p["re"], "I2 pub %.2f / rec %.2f" % (g["I2"], p["I2"])) assert p["k"] == g["k"] # membership audit against row attributes attr = {r["key"]: r for r in U["set"]} pub_year = {k: int(''.join(c for c in k.split("(")[-1] if c.isdigit())[:4]) for k in KEYS} checks = { "region: Europe": [k for k in KEYS if attr[k]["region"] == "Europe"], "region: N America": [k for k in KEYS if attr[k]["region"] == "N America"], "region: Asia": [k for k in KEYS if attr[k]["region"] == "Asia"], "women-only cohorts": [k for k in KEYS if attr[k]["sex"] == "women"], "mixed/men cohorts": [k for k in KEYS if attr[k]["sex"] != "women"], "published ≥ 2015": [k for k in KEYS if pub_year[k] >= 2015], "published < 2015": [k for k in KEYS if pub_year[k] < 2015], "≥ 1000 cases": [k for k in KEYS if attr[k]["cases"] >= 1000], "< 1000 cases": [k for k in KEYS if attr[k]["cases"] < 1000], } out() out("Membership audit (subgroup lists rebuilt from the row attributes region/sex/year/cases):") for n, lst in checks.items(): same = set(lst) == set(U["subgroups"][n]["studies"]) out("- %s: %s" % (n, "identical" if same else "DIFFERS: rebuilt %s vs json %s" % ( sorted(set(lst) - set(U["subgroups"][n]["studies"])), sorted(set(U["subgroups"][n]["studies"]) - set(lst))))) part = [k for k in KEYS if "partial" in attr[k]["expo"] or "undefined" in attr[k]["expo"] or "EPIC definition" in attr[k]["expo"] or "incl. hamburger" in attr[k]["expo"] or k.startswith("NIH-AARP")] out("- partial/undefined (rebuilt from exposure text + NIH-AARP 'partial' flag in estimates.csv): %s" % ( "identical" if set(part) == set(U["subgroups"]["partial / undefined exposure rows only"]["studies"]) else "DIFFERS %s" % sorted(set(part) ^ set(U["subgroups"]["partial / undefined exposure rows only"]["studies"])))) out("- unprocessed + partial lists are complementary and cover all 20: %s" % ( set(U["subgroups"]["explicitly unprocessed (database flag = no; excludes partial/undefined rows)"]["studies"]) | set(U["subgroups"]["partial / undefined exposure rows only"]["studies"]) == set(KEYS))) out() # ----------------------------------------------------------------------------- 2. sensitivity out("## 2. Sensitivity analyses") # rebuild Wang 2020 and Vulcan 2017 rows from estimates.csv ex = 100 / (3 * 120 / 7) wang = comb_fixed([powr(e("wang2020_075"), ex), powr(e("wang2020_086"), ex)]) def gl_slope(doses, rr, lo, hi, cases, denom, kind): """Greenland-Longnecker. doses incl. reference (first). Returns (b, var_b).""" L = np.log(np.array(rr[1:], float)) v = np.array([yv(r, l, h)[1] for r, l, h in zip(rr[1:], lo[1:], hi[1:])]) A = np.array(cases, float) N = np.array(denom, float) Ap = A.sum() x = np.array(doses[1:], float) - doses[0] if kind == "ir": # person-time / large cohort: A_j = A0 * RR_j * N_j/N_0, sum A = Ap RR = np.concatenate([[1.0], np.exp(L)]) A0 = Ap / (RR * N / N[0]).sum() c = 1 / A0 C = np.full((len(L), len(L)), c) else: # 'ci' (cumulative incidence, N_j fixed) or 'cc' (case-control, B_j = controls fixed) Bp = (N - A).sum() if kind == "ci" else N.sum() def fit(a0): if kind == "ci": # A_j/(N_j-A_j) = OR... use RR on risk scale: A_j/N_j = RR_j * A0/N0 Aj = np.concatenate([[a0], np.exp(L) * a0 * N[1:] / N[0]]) return Aj, N - Aj # case-control: B_j fixed to controls; A_j = OR_j * A0 * B_j/B_0 B = N Aj = np.concatenate([[a0], np.exp(L) * a0 * B[1:] / B[0]]) return Aj, B a0 = optimize.brentq(lambda a: fit(a)[0].sum() - Ap, 1e-9, Ap) Aj, Bj = fit(a0) if kind == "ci": c = 1 / Aj[0] - 1 / N[0] # covariance of log RRs on the risk scale else: c = 1 / Aj[0] + 1 / Bj[0] C = np.full((len(L), len(L)), c) np.fill_diagonal(C, v) Ci = np.linalg.inv(C) var = 1 / (x @ Ci @ x) b = var * (x @ Ci @ L) return float(b), float(var) def rr100(b, var, per=100.0): s = math.sqrt(var) return (math.exp(per * b), math.exp(per * (b - Z * s)), math.exp(per * (b + Z * s))) def cat_rows(set_id): rs = sorted([r for r in EST.values() if r["category_set_id"] == set_id], key=lambda r: int(r["category_rank"])) return rs def vulcan_sex(sex): # ids exactly as listed in the json provenance string for the Vulcan row ids = {"men": [96, 102, 108, 114, 120], "women": [94, 100, 106, 112, 118]}[sex] rs = [EST["vulcan2017_%03d" % i] for i in ids] assert all(r["sex"] == sex and r["exposure"] == "red_unprocessed" for r in rs) d = [float(r["dose_point"]) for r in rs] rr = [float(r["point"]) for r in rs] lo = [float(r["ci_low"]) if r["ci_low"] else 1 for r in rs] hi = [float(r["ci_high"]) if r["ci_high"] else 1 for r in rs] return d, rr, lo, hi, [float(r["n_cases"]) for r in rs], [float(r["person_years"]) for r in rs], rs vrows = [] for sx in ("men", "women"): d, rr, lo, hi, A, N, rs = vulcan_sex(sx) b, var = gl_slope(d, rr, lo, hi, A, N, "ir") vrows.append(rr100(b, var)) vulcan = comb_fixed(vrows) out(HDR) cmp("Wang 2020 row (rebuilt from wang2020_075/086, ^%.3f, fixed-combined)" % ex, U["wang"]["row"], wang) cmp("Vulcan 2017 row (rebuilt: GL per sex from 'Full model' categories, fixed-combined)", U["vulcan"]["row"], vulcan, "ids used: " + ";".join(r["estimate_id"] for r in vulcan_sex("men")[6] + vulcan_sex("women")[6])) specs = { "without Million Women (calibrated)": [k for k in KEYS if not k.startswith("Million")], "without NIH-AARP (energy-density assumption)": [k for k in KEYS if not k.startswith("NIH")], "without nested case-control (Tiemersma)": [k for k in KEYS if "Tiemersma" not in k], "without the two largest (MWS, NIH-AARP)": [k for k in KEYS if not (k.startswith("NIH") or k.startswith("Million"))], } for name, g in U["sensitivity"].items(): if name in specs: rows = [SET[k] for k in specs[name]] elif name.startswith("with Malmö"): rows = MAIN + [tuple(U["vulcan"]["row"])] elif name.startswith("NHS/HPFS"): rows = [SET[k] if not k.startswith("NHS + HPFS") else tuple(U["wang"]["row"]) for k in KEYS] p = dl(rows) assert p["k"] == g["k"], name cmp(name[:60] + " (k=%d)" % g["k"], g["re"], p["re"], "I2 pub %.2f / rec %.2f" % (g["I2"], p["I2"])) out() # ----------------------------------------------------------------------------- 3. negative controls out("## 3. Negative controls (rows rebuilt from estimates.csv; json lists study ids only)") NEG = { "poultry": [("dankner2024", "dankner2024_040", 1), ("english2004", "english2004_035", 700 / 120), ("gilsing2015", "gilsing2015_256", 2), ("mejborn2021", "mejborn2021_018", 1), ("norat2005", "norat2005_048", 1), ("spencer2010", "spencer2010_050", 2), ("bradbury2020", "bradbury2020_020", 4)], "fish": [("english2004", "english2004_050", 700 / 120), ("gilsing2015", "gilsing2015_316", 4), ("norat2005", "norat2005_049", 1), ("spencer2010", "spencer2010_065", 2), ("bradbury2020", "bradbury2020_025", 4)], } out("Row choice: the main-model continuous (per-increment) row per cohort, converted to per 100 g/day by " "exponentiation (per 50 g ^2, per 25 g ^4, per serving/week ^(700/120)=5.833 at 120 g/serving, per 100 g ^1).") out() out("| Control | Study | estimate_id | exposure | as printed | per 100 g |") out("|---|---|---|---|---|---|") for ctl, lst in NEG.items(): for st, eid, p_ in lst: r = EST[eid] assert r["study_id"] == st and r["exposure"] == ctl and r["outcome"] == "colorectal" out("| %s | %s | %s | %s | %s %s (%s-%s) | %.4f (%.4f-%.4f) |" % ( ctl, st, eid, r["exposure"], r["dose_unit"], r["point"], r["ci_low"], r["ci_high"], *powr(e(eid), p_))) out() out(HDR) for ctl, lst in NEG.items(): g = U["negative_controls"][ctl] assert [s for s, _, _ in lst] == g["studies"] p = dl([powr(e(eid), p_) for _, eid, p_ in lst]) cmp("%s (k=%d)" % (ctl, g["k"]), g["re"], p["re"], "I2 pub %.2f / rec %.2f" % (g["I2"], p["I2"])) out() # ----------------------------------------------------------------------------- 4. LOO, Egger, HKSJ, FE, REML out("## 4. Leave-one-out, Egger, HKSJ, fixed effect, REML") out(HDR) ok_loo = True loo_re = [] for slug, key in SLUG.items(): p = dl([SET[k] for k in KEYS if k != key]) loo_re.append(p["re"][0]) d = max(abs(a - b) for a, b in zip(m["loo"][slug], p["re"])) if d >= TOL4: ok_loo = False cmp("LOO without " + slug, m["loo"][slug], p["re"]) out("| LOO, all 20 runs (point and CI) | json `loo` | recomputed | %s | %s | |" % ( "yes", "yes, max diff < 5e-5" if ok_loo else "NO (see rows above)")) out("| LOO range of point estimates | 1.09-1.11 | %.4f-%.4f | %s | | |" % ( min(loo_re), max(loo_re), "yes" if (round(min(loo_re), 2), round(max(loo_re), 2)) == (1.09, 1.11) else "NO")) eg = egger(P["ys"], P["vs"]) cmp_scalar("Egger intercept", m["egger"][0], eg[0]) cmp_scalar("Egger p", m["egger"][1], eg[1]) cmp("HKSJ (Hartung-Knapp, t df=19)", m["hksj"], P["hksj"], "truncated max(1,q) variant: " + f3(P["hksj_trunc"])) cmp("Fixed effect", m["fe"], P["fe"]) rm = reml(P["ys"], P["vs"]) cmp("REML", m["reml"], rm, "tau2 pub %.6f / rec %.6f" % (m["reml"][3], rm[3])) out() # ----------------------------------------------------------------------------- 6. GL (needed for multiverse Wada option) SP = U["spline"] def gl_set(s, mode): """mode = dict of assumptions for missing data. Returns (b, var, rule).""" rs = cat_rows(s["set"]) rr = [float(r["point"]) for r in rs] lo = [float(r["ci_low"]) if r["ci_low"] else 1.0 for r in rs] hi = [float(r["ci_high"]) if r["ci_high"] else 1.0 for r in rs] doses = s["doses"] cases = [float(r["n_cases"]) if r["n_cases"] else None for r in rs] den = None rule = [] if s["denom"] == "person_years" and all(r["person_years"] for r in rs): den, kind = [float(r["person_years"]) for r in rs], "ir" elif s["denom"] == "n_participants": den, kind = [float(r["n_in_category"]) for r in rs], "ir" elif s["denom"] == "n_controls": den, kind = [float(r["n_in_category"]) for r in rs], "cc" else: kind = "ir" rule.append("denominators missing -> equal per category") if any(c is None for c in cases): tot = mode["total_cases"].get(s["set"]) cases = [tot / len(rs)] * len(rs) rule.append("cases missing -> total %s split equally" % tot) if den is None: den = [mode["N"]] * len(rs) b, var = gl_slope(doses, rr, lo, hi, cases, den, kind) return b, var, kind, "; ".join(rule), [r["estimate_id"] for r in rs] TOT = {"larsson2005_mv_colorectal_women_beefpork": 733, "alrajabi2021_m2_colorectal_men_red": 232 * 118 / 265, "alrajabi2021_m2_colorectal_women_red": 232 * 147 / 265} MODE = dict(total_cases=TOT, N=1e6) GL = {} for s in SP["sets"]: GL[s["set"]] = gl_set(s, MODE) # per-study GL rows vs the main-set rows (same category data, 100 g) out("## 6a. GL slopes per study vs the per-study rows in the main pool (same category sets)") out(HDR) pairs = [("pietinen1999_mv_colorectal_men_red_unprocessed", "ATBC (Pietinen 1999)"), ("tiemersma2002_mv_colorectal_all_red", "Netherlands MP (Tiemersma 2002)"), ("larsson2005_mv_colorectal_women_beefpork", "SMC (Larsson 2005)"), ("ollberding2012_redunp_mv_colorectal_all", "MEC (Ollberding 2012)"), ("knuppel2020_m2_colorectal_all_red_unprocessed", "UK Biobank (Knuppel 2020)"), ("yiannakou2022_mv_colorectal_women_red_unprocessed", "BWHS (Yiannakou 2022)")] for sid, key in pairs: b, var, kind, rule, ids = GL[sid] cmp("GL %s" % key, SET[key], rr100(b, var), "%s %s" % (kind, rule)) al = comb_fixed([rr100(*GL["alrajabi2021_m2_colorectal_men_red"][:2]), rr100(*GL["alrajabi2021_m2_colorectal_women_red"][:2])]) cmp("GL Alberta's Tomorrow (Al Rajabi 2021), sexes fixed-combined", SET["Alberta's Tomorrow (Al Rajabi 2021)"], al, "cases assumed: 232 split by model-1 sex ratio, then equally across quartiles") out() out("## 6b. Linear two-stage slope on the 12 spline sets (GL per set, DL pooling)") out(HDR) ys, vs = [], [] for s in SP["sets"]: b, var = GL[s["set"]][:2] ys.append(100 * b) vs.append(1e4 * var) pl = dl_yv(ys, vs) cmp("linear slope per 100 g, 12 sets", SP["linear_same_sets"], pl["re"], "Q %.2f, I2 %.1f; assumptions for missing data listed in RAPPORT" % (pl["Q"], pl["I2"])) out("| per-set slopes (per 100 g) | | " + "; ".join("%s %.3f" % (s["set"].split("_")[0] + "/" + s["sex"], math.exp(y)) for s, y in zip(SP["sets"], ys)) + " | | | |") out("## 6c. Impact: main pool with this checker's GL rows in place of the json GL rows") out(HDR) alt = dict(SET) for sid, key in pairs: alt[key] = rr100(*GL[sid][:2]) alt["Alberta's Tomorrow (Al Rajabi 2021)"] = al pa = dl([alt[k] for k in KEYS]) cmp("main pool, 7 GL rows replaced", U["main"]["re"], pa["re"], "I2 %.1f" % pa["I2"]) ycases = sum(float(r["n_cases"]) for r in cat_rows("yiannakou2022_mv_colorectal_women_red_unprocessed")) out("| BWHS cases: json row %d vs sum of category cases in estimates.csv %d | | | | | data note |" % ( attr["BWHS (Yiannakou 2022)"]["cases"], ycases)) out() # ----------------------------------------------------------------------------- 5. multiverse out("## 5. Multiverse (alternative rows rebuilt from estimates.csv / cup2016.json; json gives option names only)") cup16 = json.load(open(BASE + "/kontroll/reanalyser/cup2016.json")) def find_cup_row(label): for run in cup16["runs"] if "runs" in cup16 else []: for r in run.get("rows", []): if r["label"] == label and r.get("note", "").startswith("CUP"): return (r["rr"], r["lo"], r["hi"]) s = json.dumps(cup16) raise KeyError(label) def first_cup(label): def walk(o): if isinstance(o, dict): if o.get("label") == label and str(o.get("note", "")).startswith("CUP"): yield (o["rr"], o["lo"], o["hi"]) for v in o.values(): yield from walk(v) elif isinstance(o, list): for v in o: yield from walk(v) vals = set(walk(cup16)) return sorted(vals) wada = comb_fixed([rr100(*GL["wada2017_mv_colorectal_men_red_unprocessed"][:2]), rr100(*GL["wada2017_mv_colorectal_women_red_unprocessed"][:2])]) nih = e("zouiouich2025_134") OPT = { "Netherlands MP (Tiemersma 2002)": {"main": SET["Netherlands MP (Tiemersma 2002)"], "CUP value": first_cup("tiemersma2002")[0]}, "MCCS (English 2004)": {"main": SET["MCCS (English 2004)"], "adjusted for processed meat": powr(e("english2004_076"), 700 / 120)}, "EPIC (Papadimitriou 2022)": {"main": SET["EPIC (Papadimitriou 2022)"], "Norat 2005 observed": e("norat2005_031"), "Norat 2005 calibrated": e("norat2005_032")}, "SMC (Larsson 2005)": {"main": SET["SMC (Larsson 2005)"], "total red meat per 100 g (paper)": e("larsson2005_085")}, "MEC (Ollberding 2012)": {"main": SET["MEC (Ollberding 2012)"], "total red meat (CUP's row)": first_cup("ollberding2012")[0]}, "JPHC + Takayama (Islam 2019)": {"main": SET["JPHC + Takayama (Islam 2019)"], "Wada 2017 (Takayama only)": wada}, "UK Biobank (Knuppel 2020)": {"main": SET["UK Biobank (Knuppel 2020)"], "Bradbury 2020 (same cohort)": powr(e("bradbury2020_010"), 2)}, "NIH-AARP (Zouiouich 2025)": {"main": SET["NIH-AARP (Zouiouich 2025)"], "at 1800 kcal": powr(nih, 100 / 90), "at 2200 kcal": powr(nih, 100 / 110)}, } # NIH-AARP: per 50 g/1000 kcal. At E kcal/day, 50 g/1000 kcal = 0.05*E g/day, so per 100 g/day the exponent # is 100/(0.05*E): 1800 kcal -> 1.111, 2200 kcal -> 0.909 (used above). A run with the exponents the other # way round (0.9 / 1.1) is reported as a diagnostic, because that is what the published medians match. OPT_SWAP = json.loads(json.dumps(OPT)) OPT_SWAP["NIH-AARP (Zouiouich 2025)"]["at 1800 kcal"] = powr(nih, 0.9) OPT_SWAP["NIH-AARP (Zouiouich 2025)"]["at 2200 kcal"] = powr(nih, 1.1) out("CUP rows found in cup2016.json: tiemersma2002 %s; ollberding2012 %s" % (first_cup("tiemersma2002"), first_cup("ollberding2012"))) out() out("| Dimension | Option | Row used (per 100 g) | Source |") out("|---|---|---|---|") SRC = {"CUP value": "cup2016.json CUP row", "adjusted for processed meat": "english2004_076 ^5.833", "Norat 2005 observed": "norat2005_031", "Norat 2005 calibrated": "norat2005_032", "total red meat per 100 g (paper)": "larsson2005_085", "total red meat (CUP's row)": "cup2016.json CUP row", "Wada 2017 (Takayama only)": "GL wada2017_013-016 / 085-088, sexes fixed-combined", "Bradbury 2020 (same cohort)": "bradbury2020_010 ^2", "at 1800 kcal": "zouiouich2025_134 ^(100/90)", "at 2200 kcal": "zouiouich2025_134 ^(100/110)"} for dim, opts in OPT.items(): assert list(opts) == U["multiverse"]["dims"][dim], dim for o, row in opts.items(): out("| %s | %s | %.4f (%.4f-%.4f) | %s |" % (dim, o, *row, SRC.get(o, "json main row"))) out() dims = list(OPT) def run_mv(OPT): res = [] for combo in itertools.product(*[list(OPT[d]) for d in dims]): rows = dict((k, SET[k]) for k in KEYS) for d, o in zip(dims, combo): rows[d] = tuple(OPT[d][o]) p = dl(list(rows.values())) for model in ("random", "fixed"): est = p["re"] if model == "random" else p["fe"] res.append((dict(zip(dims, combo), model=model), est)) return res def report_mv(res): n = len(res) sig = sum(1 for _, r in res if r[1] > 1) pts = [r[0] for _, r in res] mv = U["multiverse"] out(HDR) cmp_scalar("n combinations", mv["n"], n) cmp_scalar("n with lower CI > 1", mv["significant"], sig) cmp_scalar("median RR", mv["median"], statistics.median(pts)) cmp_scalar("min RR", mv["min"], min(pts)) cmp_scalar("max RR", mv["max"], max(pts)) for lab, bc in mv["by_choice"].items(): dim, opt = lab.split(" = ") sub = [r[0] for c, r in res if (c["model"] == opt if dim == "model" else c[dim] == opt)] sg = sum(1 for c, r in res if (c["model"] == opt if dim == "model" else c[dim] == opt) and r[1] > 1) cmp_scalar("by_choice " + lab + " median (n %d/%d, sig %d/%d)" % (bc["n"], len(sub), bc["significant"], sg), bc["median"], statistics.median(sub)) out() for tag, OO in (("A. NIH exponents as derived (1800 kcal ^1.111, 2200 kcal ^0.909)", OPT), ("B. NIH exponents reversed (1800 kcal ^0.9, 2200 kcal ^1.1) - diagnostic", OPT_SWAP)): out("### " + tag) report_mv(run_mv(OO)) # ----------------------------------------------------------------------------- 7. wcrf2007 Fig 4.3.3 out("## 7. wcrf2007 'Per 100 g/day - Fig 4.3.3' (all.json)") out(HDR) runs = ALL["wcrf2007"] runs = runs["runs"] if isinstance(runs, dict) and "runs" in runs else runs cand = [] def walk_runs(o): if isinstance(o, dict): if str(o.get("name", "")).startswith("Per 100 g/day") and "rows" in o: cand.append(o) for v in o.values(): walk_runs(v) elif isinstance(o, list): for v in o: walk_runs(v) walk_runs(ALL["wcrf2007"]) for run in cand: rows = [(r["rr"], r["lo"], r["hi"]) for r in run["rows"]] p = dl(rows) out("| run: %s | | | | | |" % run["name"]) cmp(" json re", run["re"], p["re"]) cmp(" json fe", run["fe"], p["fe"]) cmp_scalar(" Q", run["Q"], p["Q"]) cmp_scalar(" tau2", run["tau2"], p["tau2"]) cmp_scalar(" I2", run["I2"], p["I2"]) cmp(" published figure 1.29 (1.05-1.59) vs FE", [1.29, 1.05, 1.59], p["fe"], "2-dp target only") cmp(" published text 1.29 (1.04-1.60) vs RE", [1.29, 1.04, 1.60], p["re"], "2-dp target only") for slug, v in run.get("loo", {}).items(): q = dl([(r["rr"], r["lo"], r["hi"]) for r in run["rows"] if r["label"] != slug]) cmp(" LOO without " + slug, v, q["re"])