import os #!/usr/bin/env python3 """K11 meta_check.py - independent DerSimonian-Laird re-pooling of database/synthesis_extraction.csv. Written from scratch by the K11 controller (no reuse of scripts/pool.py or kontroll/scripts/pool_synth.py). Method y_i = ln(RR_i) se_i = (ln(ci_high) - ln(ci_low)) / (2*1.96) fixed: w_i = 1/se_i^2 ; Q = sum w_i (y_i - y_FE)^2 ; df = k-1 DL: tau2 = max(0, (Q-df) / (sum w - sum w^2/sum w)) RE: w*_i = 1/(se_i^2 + tau2) ; y_RE = sum w* y / sum w* ; se_RE = sqrt(1/sum w*) I2 = max(0, (Q-df)/Q) Per-study RE weight % = 100 * w*_i / sum w*. Part A: every synth_id x outcome x analysis x increment group with usable rr/ci (naive grouping, as the table's own columns define it). Part B: the published analyses, with rows selected explicitly by CSV line number (header = line 1) and guarded by study_id, compared to the published pooled value (and printed weights where they exist). Usage: python3 meta_check.py [path/to/synthesis_extraction.csv] """ import csv import math import sys from collections import OrderedDict CSV = sys.argv[1] if len(sys.argv) > 1 else os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "database", "synthesis_extraction.csv") Z = 1.959964 def load(path): rows = {} with open(path, newline="") as f: for i, r in enumerate(csv.DictReader(f), start=2): rows[i] = r return rows def fnum(x): try: return float(x) except (TypeError, ValueError): return None def pool(items): """items: list of (label, rr, lo, hi). Returns dict or None.""" ys, vs, labs = [], [], [] for lab, rr, lo, hi in items: if rr is None or lo is None or hi is None or lo <= 0 or hi <= lo: continue ys.append(math.log(rr)) se = (math.log(hi) - math.log(lo)) / (2 * Z) vs.append(se * se) labs.append(lab) k = len(ys) if k == 0: return None w = [1 / v for v in vs] sw = sum(w) yfe = sum(wi * yi for wi, yi in zip(w, ys)) / sw Q = sum(wi * (yi - yfe) ** 2 for wi, yi in zip(w, ys)) df = k - 1 C = sw - sum(wi * wi for wi in w) / sw tau2 = max(0.0, (Q - df) / C) if (k > 1 and C > 0) else 0.0 ws = [1 / (v + tau2) for v in vs] sws = sum(ws) yre = sum(wi * yi for wi, yi in zip(ws, ys)) / sws sere = math.sqrt(1 / sws) I2 = max(0.0, (Q - df) / Q) * 100 if Q > 0 else 0.0 sefe = math.sqrt(1 / sw) return dict(k=k, rr=math.exp(yre), lo=math.exp(yre - Z * sere), hi=math.exp(yre + Z * sere), fe=math.exp(yfe), fe_lo=math.exp(yfe - Z * sefe), fe_hi=math.exp(yfe + Z * sefe), Q=Q, df=df, I2=I2, tau2=tau2, weights=[(lab, 100 * wi / sws) for lab, wi in zip(labs, ws)]) def fmt(p): if p is None: return "no usable rows" return ("k=%d RE %.2f (%.2f-%.2f) FE %.2f (%.2f-%.2f) Q=%.2f df=%d I2=%.1f%% tau2=%.4f" % (p["k"], p["rr"], p["lo"], p["hi"], p["fe"], p["fe_lo"], p["fe_hi"], p["Q"], p["df"], p["I2"], p["tau2"])) def item(rows, ln, study=None, override=None): r = rows[ln] if study is not None and r["study_id"] != study: raise SystemExit("line %d is %s, expected %s - table changed, update selection" % (ln, r["study_id"], study)) rr, lo, hi = fnum(r["rr"]), fnum(r["ci_low"]), fnum(r["ci_high"]) if override: rr, lo, hi = override return ("L%d %s %s" % (ln, r["study_id"], r["outcome"]), rr, lo, hi, fnum(r["weight_pct"])) def run_set(rows, name, published, spec, extra=None): """spec: list of (line, study_id[, override tuple or None[, printed weight for this analysis]]).""" its = [] for s in spec: ln, st = s[0], s[1] ov = s[2] if len(s) > 2 else None it = item(rows, ln, st, ov) if len(s) > 3: # printed weight for THIS analysis (row's weight_pct belongs to another figure) it = it[:4] + (s[3],) its.append(it) if extra: its += extra # (label, rr, lo, hi, printed_weight) p = pool([(a, b, c, d) for a, b, c, d, _ in its]) print("-" * 100) print(name) print(" published : %s" % published) print(" recomputed: %s" % fmt(p)) if p: pw = {a: e for a, b, c, d, e in its} printed = any(e is not None for e in pw.values()) for lab, w in p["weights"]: if printed and pw.get(lab) is not None: print(" %-40s RE wt %6.2f%% printed %6.2f%% diff %+5.2f" % (lab, w, pw[lab], w - pw[lab])) else: print(" %-40s RE wt %6.2f%%" % (lab, w)) dropped = [a for a, b, c, d, e in its if a not in dict(p["weights"])] if dropped: print(" (no usable CI, not pooled: %s)" % ", ".join(dropped)) return p def part_a(rows): print("=" * 100) print("PART A - naive groups: synth_id x outcome x analysis x increment (all rows with rr+ci)") print("NB: these groups mix analyses the table does not distinguish (see K11 report).") g = OrderedDict() for ln, r in rows.items(): key = (r["synth_id"], r["outcome"], r["analysis"], r["increment"]) g.setdefault(key, []).append(item(rows, ln)) for key, its in g.items(): usable = [x for x in its if x[1] is not None and x[2] is not None and x[3] is not None] if not usable: print("%-60s n_rows=%d (no usable rr/ci)" % (" | ".join(key), len(its))) continue p = pool([(a, b, c, d) for a, b, c, d, _ in usable]) print("%-60s n_rows=%d %s" % (" | ".join(key), len(its), fmt(p))) def part_b(rows): print() print("=" * 100) print("PART B - published analyses re-pooled from explicitly selected rows") # ---------------- Norat 2002 ----------------- norat = [(6, "willett1990"), (7, "giovannucci1994"), (10, "bostick1994"), (13, "kato1997"), (14, "chen1998"), (15, "hsing1998"), (16, "singh1998"), (18, "pietinen1999")] run_set(rows, "Norat 2002 Table I / Fig 2 - red meat HvL, cohort (9 studies; Thun has no CI -> 8 pooled)", "1.27 (1.11-1.45), n=9, p-het 0.45", norat) run_set(rows, "Norat 2002 sensitivity: + Thun 1992 as read off Fig 2 (approx 1.14 (0.93-1.39), digitised by K11 from norat_p-05.png pixels)", "1.27 (1.11-1.45)", norat, extra=[("Thun 1992 digitised", 1.14, 0.93, 1.39, None)]) norat_alt = [s for s in norat if s[1] != "singh1998"] + [(16, "singh1998", (1.50, 0.95, 2.36))] run_set(rows, "Norat 2002 sensitivity: + Thun digitised, Singh as plotted in Fig 2 (approx 1.50 (0.95-2.36), digitised)", "1.27 (1.11-1.45)", norat_alt, extra=[("Thun 1992 digitised", 1.14, 0.93, 1.39, None)]) # ---------------- Larsson & Wolk 2006 ----------------- lw_main = [(19, "bostick1994"), (20, "kato1997"), (21, "chen1998"), (22, "hsing1998"), (24, "singh1998"), (25, "pietinen1999"), (26, "jarvinen2001"), (29, "tiemersma2002"), (30, "flood2003"), (31, "wei2004"), (34, "wei2004"), (37, "english2004"), (40, "larsson2005"), (45, "chao2005"), (50, "norat2005")] run_set(rows, "Larsson 2006 Fig 1 / text p2660 - red meat HvL CRC (15 studies)", "1.28 (1.15-1.42), Q=4.86, p=0.99, I2=0%", lw_main) lw_colon10 = [(19, "bostick1994"), (23, "hsing1998"), (24, "singh1998"), (27, "jarvinen2001"), (32, "wei2004"), (35, "wei2004"), (38, "english2004"), (41, "larsson2005"), (46, "chao2005"), (51, "norat2005")] run_set(rows, "Larsson 2006 Table II colon - all 10 colon rows in Table I", "1.21 (1.05-1.40), n=9, Q=2.35", lw_colon10) for drop in lw_colon10: sub = [s for s in lw_colon10 if s != drop] p = pool([item(rows, s[0], s[1])[:4] for s in sub]) print(" leave out L%d %-15s -> %s" % (drop[0], drop[1], fmt(p))) run_set(rows, "Larsson 2006 Table II rectum (7 studies)", "1.56 (1.25-1.95), Q=4.82, I2=0%", [(28, "jarvinen2001"), (33, "wei2004"), (36, "wei2004"), (39, "english2004"), (44, "larsson2005"), (49, "chao2005"), (54, "norat2005")]) run_set(rows, "Larsson 2006 Table II proximal colon (3)", "1.15 (0.91-1.47), Q=0.72", [(42, "larsson2005"), (47, "chao2005"), (52, "norat2005")]) run_set(rows, "Larsson 2006 Table II distal colon (3)", "1.31 (0.62-2.79), Q=12.36, I2=83.8%", [(43, "larsson2005"), (48, "chao2005"), (53, "norat2005")]) # ---------------- WCRF 2007 ----------------- run_set(rows, "WCRF 2007 Fig 4.3.2 - per time/week (9 lines, 7 studies)", "1.43 (1.05-1.94); I2=58% p=0.01 (CUP 2010 Table 36)", [(82, "willett1990"), (83, "bostick1994"), (84, "giovannucci1994"), (85, "singh1998"), (86, "chen1998"), (87, "chen1998"), (88, "tiemersma2002"), (89, "tiemersma2002"), (90, "english2004")]) run_set(rows, "WCRF 2007 Fig 4.3.3 - per 100 g/day (Larsson CI printed 1.05-0.95 -> not usable)", "1.29 (1.05-1.59) fig / 1.29 (1.04-1.60) text / 1.29 (0.94-1.78) CUP T36; I2=32%", [(91, "pietinen1999"), (92, "norat2005"), (93, "larsson2005")]) run_set(rows, "WCRF 2007 Fig 4.3.3 sensitivity - Larsson upper CI imputed symmetric on log scale (1.43^2/1.05=1.95)", "1.29 (1.05-1.59)", [(91, "pietinen1999"), (92, "norat2005"), (93, "larsson2005", (1.43, 1.05, 1.43 ** 2 / 1.05))]) # ---------------- CUP 2010 SLR / Chan 2011 ----------------- run_set(rows, "CUP 2010 SLR Fig 70 (= Chan Fig 4 CRC) - DR per 100 g/day CRC", "1.17 (1.05-1.31), I2=0%, p=0.483", [(96, "lee2009"), (97, "nothlings2009"), (98, "larsson2005"), (99, "norat2005"), (100, "english2004"), (101, "tiemersma2002"), (102, "jarvinen2001"), (103, "pietinen1999")]) run_set(rows, "CUP 2010 SLR Fig 71 female subtotal (Lee, Larsson from Fig 70 rows; Tiemersma F)", "1.05 (0.78-1.42), I2=21.9%; printed wts Lee 35.70, Larsson 55.60, Tiemersma 8.71", [(96, "lee2009", None, 35.70), (98, "larsson2005", None, 55.60), (104, "tiemersma2002")]) run_set(rows, "CUP 2010 SLR Fig 71 male subtotal (Tiemersma M; Pietinen from Fig 70 row)", "1.28 (0.49-3.35), I2=64.3%; printed wts Tiemersma 43.79, Pietinen 56.21", [(105, "tiemersma2002"), (103, "pietinen1999", None, 56.21)]) run_set(rows, "CUP 2010 SLR Fig 74 - DR per 100 g/day colon", "1.12 (0.97-1.29), I2=0%, p=0.888", [(106, "lee2009"), (107, "oba2006"), (108, "norat2005"), (109, "english2004"), (110, "wei2004"), (111, "wei2004"), (112, "jarvinen2001"), (113, "singh1998"), (114, "bostick1994")]) run_set(rows, "CUP 2010 SLR Fig 75 female subtotal", "1.00 (0.72-1.38); printed wts Lee 26.64, Oba 4.66, Wei 50.94, Bostick 17.76", [(106, "lee2009", None, 26.64), (115, "oba2006"), (111, "wei2004", None, 50.94), (114, "bostick1994", None, 17.76)]) run_set(rows, "CUP 2010 SLR Fig 75 male subtotal", "1.06 (0.75-1.50); printed wts Oba 7.54, Wei 92.46", [(116, "oba2006"), (110, "wei2004", None, 92.46)]) run_set(rows, "Chan 2011 Fig 4 - DR per 100 g/day colon (10 rows incl. Larsson 2005)", "1.17 (1.02-1.33), I2=0%, p=0.646", [(155, "larsson2005"), (156, "bostick1994"), (157, "singh1998"), (158, "jarvinen2001"), (159, "english2004"), (160, "wei2004"), (161, "wei2004"), (162, "norat2005"), (163, "oba2006"), (164, "lee2009")]) run_set(rows, "CUP 2010 SLR Fig 78 (= Chan Fig 4 rectal) - DR per 100 g/day rectum", "1.18 (0.98-1.42), I2=0%, p=0.666", [(117, "lee2009"), (118, "larsson2005"), (119, "norat2005"), (120, "english2004"), (121, "wei2004"), (122, "wei2004"), (123, "jarvinen2001")]) run_set(rows, "Chan 2011 Table 1 - HvL CRC (= SLR Fig 73, 12 studies)", "1.10 (1.00-1.21), I2=22%, p=0.22", [(124, "cross2010"), (125, "lee2009"), (126, "nothlings2009"), (127, "butler2008"), (128, "chan2005"), (129, "larsson2005"), (130, "norat2005"), (131, "english2004"), (132, "tiemersma2002"), (133, "jarvinen2001"), (134, "pietinen1999"), (135, "kato1997")]) colon12 = [(136, "lee2009"), (137, "oba2006"), (138, "oba2006"), (139, "larsson2005"), (140, "larsson2005"), (141, "norat2005"), (142, "english2004"), (143, "wei2004"), (144, "wei2004"), (145, "jarvinen2001"), (146, "singh1998"), (147, "bostick1994")] run_set(rows, "Chan 2011 Table 1 - HvL colon (SLR Fig 77, all 12 lines entered separately)", "1.18 (1.04-1.35), n=10, I2=0%, p=0.70", colon12) # within-study fixed-effect combination of Oba M+F and Larsson prox+dist, then DL over 10 studies def fe_comb(lines): its = [item(rows, l)[:4] for l in lines] p = pool(its) se = (math.log(p["fe_hi"]) - math.log(p["fe_lo"])) / (2 * Z) return (p["fe"], math.exp(math.log(p["fe"]) - Z * se), math.exp(math.log(p["fe"]) + Z * se)) oba = fe_comb([137, 138]); lar = fe_comb([139, 140]) rest = [s for s in colon12 if s[0] not in (137, 138, 139, 140)] run_set(rows, "Chan 2011 HvL colon - Oba M+F and Larsson prox+dist combined (FE) within study -> 10 studies", "1.18 (1.04-1.35), n=10, I2=0%", rest, extra=[("Oba 2006 M+F FE %.2f (%.2f-%.2f)" % oba, oba[0], oba[1], oba[2], None), ("Larsson 2005 prox+dist FE %.2f (%.2f-%.2f)" % lar, lar[0], lar[1], lar[2], None)]) run_set(rows, "Chan 2011 Table 1 - HvL rectum (= SLR Fig 81, 7 lines)", "1.14 (0.83-1.56), I2=38%, p=0.14", [(148, "lee2009"), (149, "larsson2005"), (150, "norat2005"), (151, "english2004"), (152, "wei2004"), (153, "wei2004"), (154, "jarvinen2001")]) run_set(rows, "Chan 2011 Table 1 - HvL proximal colon (Larsson 2005 + Chao 2005)", "1.13 (0.83-1.54), n=2, I2=0%, p=0.83", [(168, "larsson2005"), (170, "chao2005")]) run_set(rows, "Chan 2011 Table 1 - HvL distal colon (Larsson 2005 + Chao 2005)", "1.57 (0.98-2.49), n=2, I2=53%, p=0.15", [(169, "larsson2005"), (171, "chao2005")]) print(" NB: no pair of Larsson/Chao/Norat proximal or distal values available in the sources reproduces Chan T1;") print(" Larsson+Norat distal = 1.83 (1.25-2.68), Larsson+Chao = 1.18 (0.43-3.25). Unresolved (see report).") if __name__ == "__main__": rows = load(CSV) print("meta_check.py on %s (%d data rows)" % (CSV, len(rows))) part_a(rows) part_b(rows)