#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ forking-paths.py -- Science Journaling Club, Volume 2 Issue 2, Winter 2026. QUESTION -------- Analysing a data set involves dozens of small decisions, each of which is defensible on its own. Suppose every one of them is made in good faith, but made *after* looking at the data. How often does a study of pure noise then produce a result significant at the conventional five percent level? WHAT THIS FILE IS ----------------- This is a computational study. There are no participants, no laboratory, no measurements of anything physical. Every "study" below is two lists of numbers drawn from the same normal distribution by a pseudo-random number generator, and every "finding" is therefore false by construction. The computation is the experiment. The only thing being measured is the behaviour of a decision procedure applied to noise. THE MODEL --------- One simulated study: * Two groups, "treatment" and "control", 40 participants each at full enrolment. Group membership has no effect on anything whatsoever. * Three outcome measures per participant, jointly normal, mean 0, variance 1, pairwise correlation r (default 0.5). Think of three related questionnaire scales. * One covariate per participant, standard normal, independent of everything including the outcomes. Think of a baseline measure someone thinks might matter. * One irrelevant binary attribute per participant, a fair coin. Think of a demographic split with no bearing on the question. Five researcher degrees of freedom, each applied alone and in every combination (2^5 = 32 combinations): S optional stopping test at n = 20, 25, 30, 35, 40 per group instead of only at n = 40, and stop at the first significant look M outcome choice test all three outcome measures, not only the first G subgroup analysis test the whole sample, then within each level of the irrelevant binary attribute C covariate inclusion test unadjusted, and adjusted for the covariate by ordinary least squares O outlier exclusion test with all data, and again after dropping any observation more than TRIM_K = 2.0 sample standard deviations from the mean of the sample being analysed A combination "finds" an effect if ANY analysis path it permits yields a two-sided p below 0.05. Every path is a legitimate analysis. The dishonesty lies entirely in the selection, not in any single test. The full combination permits 5 x 3 x 3 x 2 x 2 = 180 analysis paths. Because those paths share data they are far from independent, so we also report the effective number of independent tests each combination represents: k_eff = ln(1 - FPR) / ln(1 - 0.05) which is the number of genuinely independent coin flips at alpha = 0.05 that would give the same overall false positive rate. VALIDATION ---------- 1. The club's own two-sample t statistic and p value are printed beside scipy.stats.ttest_ind computed on identical data. 2. The club's covariate-adjusted (ANCOVA) t statistic is printed beside an independent LAPACK least-squares solve on identical data, with the tail probability taken from scipy rather than from this file. 3. With no researcher degrees of freedom at all the false positive rate must equal 0.05 within Monte Carlo error, printed with its confidence interval. That number is the check on everything else in the file. 4. Optional stopping alone, with K equally spaced looks at accumulating data, is compared look by look with the published table of Armitage, McPherson and Rowe (1969), J. R. Statist. Soc. A 132, 235-244. 5. Three rows of Table 1 of Simmons, Nelson & Simonsohn (2011) are reproduced exactly as that table's note specifies them, and the club's rates are printed beside the published ones at all three significance levels. 6. A convergence trace shows the running estimate settling as trials accumulate. ASSUMPTIONS ----------- * Normal errors, equal variances, independent observations. Real data are none of these things reliably. * The researcher always takes the smallest p value available. Real researchers are less systematic, so the rates here are an upper bound on the damage from these five freedoms specifically, and a lower bound on the damage from the full set of freedoms a real analysis contains. * The null is exactly true. Nothing here says anything about power, about bias in the estimate of a real effect, or about what happens when a small genuine effect is present. * Sample sizes are fixed by the peek schedule. The optional-stopping numbers depend strongly on that schedule; see the sensitivity section. LIMITATIONS, NAMED ------------------ * The five freedoms modelled here are a small sample of the published lists. Wicherts et al. (2016) enumerate 34. Transformation choice, exclusion of whole conditions, choice of test, HARKing and the file drawer are all absent from this model. * The 32 combinations are evaluated on one shared pool of simulated studies. Each rate is an unbiased estimate with a valid binomial standard error, but the rates are correlated across rows, so differences between rows are more precise than independent runs would give and the rows are not independent replications of each other. * Subgroup analyses with fewer than MIN_N = 3 observations in either arm are treated as unavailable (p = 1). At n = 20 per group this is rare but it is a modelling choice, and it makes the subgroup rate very slightly conservative. * Real analyses are not exhaustive searches. The rate at which a real researcher blunders into a false positive depends on how many paths they actually walk, which nothing here measures. REPRODUCING ----------- python forking-paths.py > forking-paths-output.txt Python 3.12, numpy required. scipy is used only for the validation print-outs and the script says so plainly if it is missing. Master seed 20260214, hard coded below; every stream is spawned from it through numpy SeedSequence, so the output is deterministic. """ import sys import time import platform import warnings import numpy as np warnings.filterwarnings("ignore") # -------------------------------------------------------------------------- # Constants. Every one of these is a modelling choice. # -------------------------------------------------------------------------- SEED = 20260214 # master seed, Winter 2026 issue ALPHA = 0.05 # nominal two-sided significance level N_MAX = 40 # participants per group at full enrolment LOOKS = (20, 25, 30, 35, 40) # the peek schedule, per group N_OUTCOMES = 3 # correlated outcome measures R_OUTCOME = 0.5 # pairwise correlation between outcomes TRIM_K = 2.0 # outlier rule: drop |z| > TRIM_K MIN_N = 3 # smallest arm a subgroup test will accept N_STUDIES = 50000 # studies in the main shared pool N_SENS = 30000 # studies per sensitivity cell N_STOP = 200000 # studies for the Armitage comparison STOP_STEP = 25 # per-group increment between equally spaced looks STOP_MAX_K = 20 # most looks in the Armitage comparison CHUNK = 10000 # studies held in memory at once N_SIMMONS = 120000 # studies for the 2011 Table 1 replication N_SIMMONS_PUB = 15000 # simulations behind the published table FREEDOMS = ("S", "M", "G", "C", "O") FREEDOM_NAME = { "S": "optional stopping (5 looks)", "M": "three outcome measures", "G": "subgroup on irrelevant split", "C": "covariate inclusion", "O": "outlier exclusion at 2 SD", } # -------------------------------------------------------------------------- # The club's own statistics. No scipy inside these: they are what is being # validated. gammaln is Lanczos; the incomplete beta is Lentz's continued # fraction, both standard textbook constructions written out here so that the # t-test has no library behind it. # -------------------------------------------------------------------------- _LANCZOS_G = 7.0 _LANCZOS_C = ( 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7, ) def gammaln(z): """log Gamma(z) for z > 0, Lanczos approximation, vectorised.""" z = np.asarray(z, dtype=float) zz = z - 1.0 x = np.full(zz.shape, _LANCZOS_C[0]) if zz.shape else np.array(_LANCZOS_C[0]) for i in range(1, 9): x = x + _LANCZOS_C[i] / (zz + i) t = zz + _LANCZOS_G + 0.5 return 0.5 * np.log(2.0 * np.pi) + (zz + 0.5) * np.log(t) - t + np.log(x) def _betacf(a, b, x, iters=160): """Continued fraction for the incomplete beta, modified Lentz.""" tiny = 1e-300 qab = a + b qap = a + 1.0 qam = a - 1.0 c = np.ones_like(x) d = 1.0 - qab * x / qap d = np.where(np.abs(d) < tiny, tiny, d) d = 1.0 / d h = d.copy() for m in range(1, iters + 1): m2 = 2.0 * m aa = m * (b - m) * x / ((qam + m2) * (a + m2)) d = 1.0 + aa * d d = np.where(np.abs(d) < tiny, tiny, d) d = 1.0 / d c = 1.0 + aa / c c = np.where(np.abs(c) < tiny, tiny, c) h = h * d * c aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)) d = 1.0 + aa * d d = np.where(np.abs(d) < tiny, tiny, d) d = 1.0 / d c = 1.0 + aa / c c = np.where(np.abs(c) < tiny, tiny, c) h = h * d * c return h def betainc_reg(a, b, x): """Regularised incomplete beta I_x(a, b), vectorised, no library call.""" a, b, x = np.broadcast_arrays( np.asarray(a, float), np.asarray(b, float), np.asarray(x, float)) a = np.array(a, dtype=float) b = np.array(b, dtype=float) x = np.clip(np.array(x, dtype=float), 0.0, 1.0) out = np.zeros_like(x) interior = (x > 0.0) & (x < 1.0) if not np.any(interior): return np.where(x >= 1.0, 1.0, 0.0) ai, bi, xi = a[interior], b[interior], x[interior] lnbt = (gammaln(ai + bi) - gammaln(ai) - gammaln(bi) + ai * np.log(xi) + bi * np.log1p(-xi)) bt = np.exp(lnbt) lower = xi < (ai + 1.0) / (ai + bi + 2.0) res = np.empty_like(xi) if np.any(lower): res[lower] = bt[lower] * _betacf(ai[lower], bi[lower], xi[lower]) / ai[lower] upper = ~lower if np.any(upper): res[upper] = 1.0 - bt[upper] * _betacf( bi[upper], ai[upper], 1.0 - xi[upper]) / bi[upper] out[interior] = np.clip(res, 0.0, 1.0) out[x >= 1.0] = 1.0 return out def t_sf2(t, df): """Two-sided tail probability P(|T_df| >= |t|).""" t = np.asarray(t, float) df = np.asarray(df, float) bad = ~np.isfinite(t) | (df <= 0) tt = np.where(bad, 0.0, t) dfs = np.where(df <= 0, 1.0, df) x = dfs / (dfs + tt * tt) p = betainc_reg(0.5 * dfs, 0.5, x) return np.where(bad, 1.0, p) def t_crit(df, alpha=ALPHA): """Two-sided critical value(s), bisection on the club's own tail function. Vectorised over df so a whole lookup table costs one bisection.""" df = np.atleast_1d(np.asarray(df, float)) lo = np.zeros_like(df) hi = np.full_like(df, 200.0) for _ in range(80): mid = 0.5 * (lo + hi) gt = t_sf2(mid, df) > alpha lo = np.where(gt, mid, lo) hi = np.where(gt, hi, mid) out = 0.5 * (lo + hi) return out if out.size > 1 else float(out[0]) MAX_DF = 240 _dfs = np.arange(MAX_DF + 1, dtype=float) _dfs[0] = 1.0 TCRIT = np.asarray(t_crit(_dfs, ALPHA)) TCRIT[0] = np.inf # df = 0 can never be significant def _sig_from(t, df_int): """Significant at ALPHA, by table lookup rather than a tail evaluation. Degrees of freedom in this study are always whole numbers, so the table is exact and the continued fraction stays out of the inner loop.""" idx = np.clip(df_int, 0, MAX_DF) return np.abs(t) >= TCRIT[idx] # -------------------------------------------------------------------------- # Masked tests. Every analysis path in this study is "take the sample selected # by this boolean mask and run one test on it", so a single masked t test and # a single masked regression cover all 180 paths. Both are written as weighted # sums of powers, which is the arithmetic any textbook gives, and is what lets # 180 paths per study run at a sensible speed. # -------------------------------------------------------------------------- def masked_ttest(y, grp, mask, want_p=False): """Pooled-variance two-sample t test, one per row, over a boolean mask.""" g0 = mask & (grp == 0) g1 = mask & (grp == 1) w0 = g0.astype(np.float64) w1 = g1.astype(np.float64) n0 = w0.sum(axis=1) n1 = w1.sum(axis=1) ok = (n0 >= MIN_N) & (n1 >= MIN_N) n0s = np.where(ok, n0, 2.0) n1s = np.where(ok, n1, 2.0) s0 = (w0 * y).sum(axis=1) s1 = (w1 * y).sum(axis=1) q0 = (w0 * y * y).sum(axis=1) q1 = (w1 * y * y).sum(axis=1) ss0 = q0 - s0 * s0 / n0s ss1 = q1 - s1 * s1 / n1s df = n0s + n1s - 2.0 sp2 = np.maximum((ss0 + ss1) / df, 1e-300) with np.errstate(divide="ignore", invalid="ignore"): t = (s0 / n0s - s1 / n1s) / np.sqrt(sp2 * (1.0 / n0s + 1.0 / n1s)) t = np.where(ok & np.isfinite(t), t, 0.0) dfi = np.where(ok, df, 0.0).astype(np.int64) if want_p: return t, dfi, np.where(ok, t_sf2(t, np.maximum(df, 1.0)), 1.0) return t, dfi def masked_ancova(y, grp, cov, mask, want_p=False): """OLS of y on [1, group, covariate]; t on the group coefficient. Solved from weighted cross-products after sweeping out the intercept, which is algebraically identical to inverting the 3x3 normal-equation matrix and a good deal faster. Validation 2 checks it against a LAPACK least-squares solve on the same data. """ w = mask.astype(np.float64) g = grp.astype(np.float64)[None, :] wg = w * g wc = w * cov wy = w * y Sw = w.sum(axis=1) Sg = wg.sum(axis=1) Sc = wc.sum(axis=1) Sy = wy.sum(axis=1) Sgg = (wg * g).sum(axis=1) Scc = (wc * cov).sum(axis=1) Syy = (wy * y).sum(axis=1) Sgc = (wg * cov).sum(axis=1) Sgy = (wg * y).sum(axis=1) Scy = (wc * y).sum(axis=1) n0 = (mask & (grp == 0)).sum(axis=1) n1 = (mask & (grp == 1)).sum(axis=1) ok = (n0 >= MIN_N) & (n1 >= MIN_N) & (Sw >= 5) Sws = np.where(ok, Sw, 5.0) a = Sgg - Sg * Sg / Sws # centred sum of squares, group b = Sgc - Sg * Sc / Sws # centred cross product d = Scc - Sc * Sc / Sws # centred sum of squares, covariate u = Sgy - Sg * Sy / Sws v = Scy - Sc * Sy / Sws yy = Syy - Sy * Sy / Sws det = a * d - b * b det = np.where(np.abs(det) < 1e-300, 1e-300, det) bg = (d * u - b * v) / det bc = (a * v - b * u) / det rss = np.maximum(yy - bg * u - bc * v, 1e-300) df = Sws - 3.0 df = np.where(df <= 0, 1.0, df) s2 = rss / df with np.errstate(divide="ignore", invalid="ignore"): t = bg / np.sqrt(np.maximum(s2 * d / det, 1e-300)) t = np.where(ok & np.isfinite(t), t, 0.0) dfi = np.where(ok, df, 0.0).astype(np.int64) if want_p: return t, dfi, np.where(ok, t_sf2(t, np.maximum(df, 1.0)), 1.0) return t, dfi # -------------------------------------------------------------------------- # One pool of simulated studies, every analysis path evaluated. # -------------------------------------------------------------------------- def chol_equicorr(k, r): R = np.full((k, k), float(r)) np.fill_diagonal(R, 1.0) return np.linalg.cholesky(R) def simulate_pool(seedseq, n_studies, r=R_OUTCOME, trim_k=TRIM_K, looks=LOOKS, n_max=N_MAX, chunk=CHUNK, progress=None): """Run a pool of studies and return sig bool array (n_studies, L, M, G, A, O), True where that single analysis path reached p < ALPHA base_p the p values of the baseline path, kept so the uniformity of the null distribution can be checked Axis order: look, outcome measure, subgroup selection (0 = whole sample, 1 and 2 = the two levels of the irrelevant split), covariate adjustment (0 = no, 1 = yes), outlier exclusion (0 = keep all, 1 = trim). """ L = len(looks) N = 2 * n_max grp = np.concatenate([np.zeros(n_max, int), np.ones(n_max, int)]) within = np.concatenate([np.arange(n_max), np.arange(n_max)]) chol = chol_equicorr(N_OUTCOMES, r) rng = np.random.default_rng(seedseq) sig = np.zeros((n_studies, L, N_OUTCOMES, 3, 2, 2), dtype=bool) base_p = np.ones(n_studies, dtype=np.float64) done = 0 while done < n_studies: c = min(chunk, n_studies - done) z = rng.standard_normal((c, N, N_OUTCOMES)) y = z @ chol.T # (c, N, 3) correlated outcomes cov = rng.standard_normal((c, N)) # irrelevant covariate sub = rng.integers(0, 2, size=(c, N)) # irrelevant binary split for li in range(L): nl = looks[li] enrolled = np.broadcast_to(within < nl, (c, N)) for gi in range(3): sel = enrolled if gi == 0 else (enrolled & (sub == (gi - 1))) wsel = sel.astype(np.float64) nsel = wsel.sum(axis=1) nsafe = np.where(nsel >= 2, nsel, 2.0) for mi in range(N_OUTCOMES): ym = np.ascontiguousarray(y[:, :, mi]) for oi in range(2): if oi == 0: mask = sel else: sm = (wsel * ym).sum(axis=1) sq = (wsel * ym * ym).sum(axis=1) mu = sm / nsafe var = (sq - sm * sm / nsafe) / (nsafe - 1.0) sd = np.sqrt(np.maximum(var, 1e-300)) mask = sel & (np.abs(ym - mu[:, None]) <= trim_k * sd[:, None]) is_base = (li == L - 1 and gi == 0 and mi == 0 and oi == 0) if is_base: t0, d0, p0 = masked_ttest(ym, grp, mask, want_p=True) base_p[done:done + c] = p0 else: t0, d0 = masked_ttest(ym, grp, mask) sig[done:done + c, li, mi, gi, 0, oi] = _sig_from(t0, d0) t1, d1 = masked_ancova(ym, grp, cov, mask) sig[done:done + c, li, mi, gi, 1, oi] = _sig_from(t1, d1) done += c if progress is not None: print(progress.format(done=done, total=n_studies)) sys.stdout.flush() return sig, base_p def combo_axes(combo, n_looks=len(LOOKS)): """Which indices along each axis a combination of freedoms permits.""" looks = list(range(n_looks)) if "S" in combo else [n_looks - 1] outs = list(range(N_OUTCOMES)) if "M" in combo else [0] grps = [0, 1, 2] if "G" in combo else [0] adjs = [0, 1] if "C" in combo else [0] trims = [0, 1] if "O" in combo else [0] return looks, outs, grps, adjs, trims def combo_hits(sig, combo): """True where a combination's permitted paths contain a significant one.""" looks, outs, grps, adjs, trims = combo_axes(combo, sig.shape[1]) sub = sig[:, looks][:, :, outs][:, :, :, grps] sub = sub[:, :, :, :, adjs][:, :, :, :, :, trims] return sub.reshape(sub.shape[0], -1).any(axis=1) def n_paths(combo, n_looks=len(LOOKS)): looks, outs, grps, adjs, trims = combo_axes(combo, n_looks) return len(looks) * len(outs) * len(grps) * len(adjs) * len(trims) def wilson(k, n, z=1.959964): """Wilson score interval, which behaves near the edges where the normal approximation does not.""" if n == 0: return (0.0, 0.0) ph = k / n d = 1.0 + z * z / n centre = (ph + z * z / (2 * n)) / d half = z * np.sqrt(ph * (1 - ph) / n + z * z / (4 * n * n)) / d return (max(0.0, centre - half), min(1.0, centre + half)) def k_effective(fpr, alpha=ALPHA): """Number of independent alpha-level tests giving the same overall rate.""" if fpr >= 1.0: return float("inf") return np.log1p(-fpr) / np.log1p(-alpha) def k_eff_se(fpr, se, alpha=ALPHA): """Delta method standard error on k_eff.""" return se / ((1.0 - fpr) * abs(np.log1p(-alpha))) # -------------------------------------------------------------------------- # The Armitage comparison: K equally spaced looks at accumulating data. # -------------------------------------------------------------------------- ARMITAGE_PUBLISHED = { 1: 0.05, 2: 0.083, 3: 0.107, 4: 0.126, 5: 0.142, 10: 0.193, 20: 0.246, } def optional_stopping_curve(seedseq, n_studies, step=STOP_STEP, kmax=STOP_MAX_K, chunk=5000, alpha=ALPHA): """False positive rate after 1..kmax equally spaced looks.""" rng = np.random.default_rng(seedseq) crit = np.asarray(t_crit( np.array([2.0 * (k + 1) * step - 2.0 for k in range(kmax)]), alpha)) hits = np.zeros(kmax, dtype=np.int64) done = 0 while done < n_studies: c = min(chunk, n_studies - done) s0 = np.zeros(c); q0 = np.zeros(c) s1 = np.zeros(c); q1 = np.zeros(c) alive = np.ones(c, dtype=bool) # not yet stopped cum = np.zeros(kmax, dtype=np.int64) stopped_at = np.full(c, -1, dtype=np.int64) for k in range(kmax): a = rng.standard_normal((c, step)) b = rng.standard_normal((c, step)) s0 += a.sum(axis=1); q0 += (a * a).sum(axis=1) s1 += b.sum(axis=1); q1 += (b * b).sum(axis=1) n = float((k + 1) * step) m0 = s0 / n; m1 = s1 / n ss0 = q0 - n * m0 * m0 ss1 = q1 - n * m1 * m1 sp2 = (ss0 + ss1) / (2 * n - 2.0) t = (m0 - m1) / np.sqrt(sp2 * (2.0 / n)) sig = np.abs(t) >= crit[k] newly = alive & sig stopped_at[newly] = k alive = alive & ~sig cum[k] = int((stopped_at >= 0).sum()) hits += cum done += c return hits / float(n_studies) # -------------------------------------------------------------------------- # An exact replication of Table 1 of Simmons, Nelson & Simonsohn (2011). # The operationalisations below are copied from that table's note, not # invented here: Situation A is three t tests, one on each of two dependent # variables correlated at r = .50 and a third on their average; Situation B is # a t test after 20 observations per cell and another after 10 more; Situation # C is a t test, an analysis of covariance with a gender main effect, and an # analysis of covariance with a gender interaction, counted significant if the # effect of condition is significant in any of them or if the gender by # condition interaction is significant. Situation D needs three conditions, # which this study's model does not have, so it is not attempted. # -------------------------------------------------------------------------- SIMMONS_PUBLISHED = { # p<.10 p<.05 p<.01 "A": (0.178, 0.095, 0.022), "B": (0.145, 0.077, 0.016), "C": (0.216, 0.117, 0.027), "AB": (0.260, 0.144, 0.033), "ABC": (0.509, 0.309, 0.084), } SIMMONS_LABEL = { "A": "two dependent variables (r = .50)", "B": "addition of 10 more observations per cell", "C": "controlling for gender or its interaction", "AB": "combine A and B", "ABC": "combine A, B and C", } SIMMONS_ALPHAS = (0.10, 0.05, 0.01) def _ols_t(X, y): """t statistics for every coefficient of a stacked OLS. X is (S, n, k).""" XtX = np.einsum("sni,snj->sij", X, X) Xty = np.einsum("sni,sn->si", X, y) beta = np.linalg.solve(XtX, Xty[:, :, None])[:, :, 0] resid = y - np.einsum("sni,si->sn", X, beta) n, kk = X.shape[1], X.shape[2] df = float(n - kk) s2 = (resid * resid).sum(axis=1) / df inv = np.linalg.inv(XtX) se = np.sqrt(np.maximum(s2[:, None] * np.diagonal(inv, axis1=1, axis2=2), 1e-300)) return beta / se, df def simmons_replication(seedseq, n_sims, chunk=10000, centred_gender=False): """Return {row: (rate at .10, rate at .05, rate at .01)}.""" n1, n2 = 20, 30 # observations per cell, before and after rho = 0.5 N = 2 * n2 idx20 = np.concatenate([np.arange(n1), np.arange(n2, n2 + n1)]) idx30 = np.arange(N) cond_full = np.concatenate([np.zeros(n2), np.ones(n2)]) chol = chol_equicorr(2, rho) rng = np.random.default_rng(seedseq) crit = {} for a in SIMMONS_ALPHAS: for df in (2 * n1 - 2, 2 * n1 - 3, 2 * n1 - 4, 2 * n2 - 2, 2 * n2 - 3, 2 * n2 - 4): crit[(a, df)] = t_crit(float(df), a) rows = ("A", "B", "C", "AB", "ABC") hits = {r: np.zeros(len(SIMMONS_ALPHAS), dtype=np.int64) for r in rows} done = 0 while done < n_sims: c = min(chunk, n_sims - done) z = rng.standard_normal((c, N, 2)) yy = z @ chol.T dv = np.stack([yy[:, :, 0], yy[:, :, 1], 0.5 * (yy[:, :, 0] + yy[:, :, 1])], axis=2) gender = rng.integers(0, 2, size=(c, N)).astype(float) # t statistics, indexed [dv variant][sample size][statistic] store = {} for si, idx in enumerate((idx20, idx30)): nn = len(idx) cond = cond_full[idx] gsub = gender[:, idx] gcen = gsub - gsub.mean(axis=1, keepdims=True) if centred_gender else gsub ones = np.ones((c, nn)) condb = np.broadcast_to(cond, (c, nn)) X2 = np.stack([ones, condb, gcen], axis=2) X3 = np.stack([ones, condb, gcen, condb * gcen], axis=2) for vi in range(3): ysub = dv[:, idx, vi] g0 = ysub[:, :nn // 2] g1 = ysub[:, nn // 2:] m = nn // 2 sp2 = (g0.var(axis=1, ddof=1) + g1.var(axis=1, ddof=1)) / 2.0 t_plain = (g0.mean(axis=1) - g1.mean(axis=1)) / np.sqrt( sp2 * (2.0 / m)) t2, df2 = _ols_t(X2, ysub) t3, df3 = _ols_t(X3, ysub) store[(vi, si)] = [ (t_plain, float(nn - 2)), (t2[:, 1], df2), (t3[:, 1], df3), (t3[:, 3], df3), ] for ai, a in enumerate(SIMMONS_ALPHAS): def sig(vi, si, k): t, df = store[(vi, si)][k] return np.abs(t) >= crit[(a, int(df))] hits["A"][ai] += int(np.sum( sig(0, 0, 0) | sig(1, 0, 0) | sig(2, 0, 0))) hits["B"][ai] += int(np.sum(sig(0, 0, 0) | sig(0, 1, 0))) hits["C"][ai] += int(np.sum( sig(0, 0, 0) | sig(0, 0, 1) | sig(0, 0, 2) | sig(0, 0, 3))) ab = np.zeros(c, dtype=bool) abc = np.zeros(c, dtype=bool) for vi in range(3): for si in range(2): ab |= sig(vi, si, 0) for kk in range(4): abc |= sig(vi, si, kk) hits["AB"][ai] += int(np.sum(ab)) hits["ABC"][ai] += int(np.sum(abc)) done += c return {r: hits[r] / float(n_sims) for r in rows} # -------------------------------------------------------------------------- # Reporting helpers # -------------------------------------------------------------------------- def rule(ch="-", n=78): print(ch * n) def head(title): print() rule("=") print(title) rule("=") def main(): t_start = time.time() print("=" * 78) print("HOW MANY REASONABLE CHOICES DOES IT TAKE TO FIND A RESULT THAT IS NOT THERE") print("Science Journaling Club | Volume 2, Issue 2, Winter 2026") print("Original research: a Monte Carlo study of researcher degrees of freedom") print("=" * 78) print("Python : %s" % sys.version.split()[0]) print("numpy : %s" % np.__version__) print("platform : %s" % platform.platform()) print("MASTER SEED : %d (every stream spawned from it via SeedSequence)" % SEED) print("nominal alpha : %.3f two-sided" % ALPHA) print("group size : %d per group at full enrolment" % N_MAX) print("peek schedule : %s per group" % (", ".join(str(x) for x in LOOKS))) print("outcomes : %d, pairwise correlation r = %.2f" % (N_OUTCOMES, R_OUTCOME)) print("outlier rule : drop |z| > %.1f SD of the sample being analysed" % TRIM_K) print("studies (pool): %d, shared by all 32 combinations" % N_STUDIES) print() print("NOTE ON WHAT IS BEING SIMULATED: there are no participants. Both groups") print("are drawn from the same distribution, so every significant result counted") print("below is false by construction. Nothing here was observed in anything.") root = np.random.SeedSequence(SEED) (ss_valid, ss_pool, ss_stop, ss_sens_r0, ss_sens_r3, ss_sens_r5, ss_sens_r8, ss_sens_t20, ss_sens_t25, ss_sens_t30, ss_simmons, ss_simmons_c) = root.spawn(12) # ------------------------------------------------------------------ head("VALIDATION 1. The club's t test against scipy.stats.ttest_ind") # ------------------------------------------------------------------ rng_v = np.random.default_rng(ss_valid) va = rng_v.standard_normal((20000, 40)) vb = rng_v.standard_normal((20000, 40)) ymix = np.concatenate([va, vb], axis=1) grp_v = np.concatenate([np.zeros(40, int), np.ones(40, int)]) mask_v = np.ones_like(ymix, dtype=bool) tc, dfc, pc = masked_ttest(ymix, grp_v, mask_v, want_p=True) try: from scipy import stats as spstats have_scipy = True except Exception as exc: # pragma: no cover have_scipy = False print("scipy not importable (%s); the library comparison is skipped and" % exc) print("this file's t test therefore stands unvalidated. Say so if you quote it.") if have_scipy: ts, ps = spstats.ttest_ind(va, vb, axis=1, equal_var=True) print("First five studies, identical data, both routines:") print() print(" study club t scipy t diff " "club p scipy p diff") for i in range(5): print(" %-5d %12.8f %12.8f %+10.2e %12.9f %12.9f %+10.2e" % (i, tc[i], ts[i], tc[i] - ts[i], pc[i], ps[i], pc[i] - ps[i])) print() print(" Over all %d studies:" % len(tc)) print(" max |t_club - t_scipy | = %.3e" % np.max(np.abs(tc - ts))) print(" max |p_club - p_scipy | = %.3e" % np.max(np.abs(pc - ps))) print(" max relative p error = %.3e" % np.max(np.abs(pc - ps) / np.maximum(ps, 1e-12))) print(" disagreements at alpha = %d" % int(np.sum((pc < ALPHA) != (ps < ALPHA)))) agree = np.max(np.abs(pc - ps)) < 1e-9 print(" VERDICT: %s" % ("agreement to better than 1e-9." if agree else "DISAGREEMENT larger than 1e-9. Do not trust this file.")) # ------------------------------------------------------------------ head("VALIDATION 2. The club's ANCOVA against an independent LAPACK solve") # ------------------------------------------------------------------ covv = rng_v.standard_normal((20000, 80)) ta, dfa, pa = masked_ancova(ymix, grp_v, covv, mask_v, want_p=True) print("The club's covariate-adjusted test solves the normal equations directly.") print("Here it is checked against numpy.linalg.lstsq, which is LAPACK's") print("SVD-based least-squares driver and shares no code with it, with the tail") print("probability taken from scipy.stats.t rather than from this file.") print("(statsmodels is installed in this environment but will not import against") print("numpy %s, so the comparison is made this way instead.)" % np.__version__) print() if have_scipy: print(" study club t_group LAPACK t diff " "club p scipy p diff") worst_t = 0.0 worst_p = 0.0 for i in range(500): X = np.column_stack([np.ones(80), grp_v.astype(float), covv[i]]) beta, _, _, _ = np.linalg.lstsq(X, ymix[i], rcond=None) resid = ymix[i] - X @ beta dfr = 80 - 3 s2 = float(resid @ resid) / dfr XtXi = np.linalg.pinv(X.T @ X) se_b = np.sqrt(s2 * XtXi[1, 1]) t_ref = beta[1] / se_b p_ref = 2.0 * spstats.t.sf(abs(t_ref), dfr) worst_t = max(worst_t, abs(ta[i] - t_ref)) worst_p = max(worst_p, abs(pa[i] - p_ref)) if i < 5: print(" %-5d %12.8f %12.8f %+10.2e %12.9f %12.9f %+10.2e" % (i, ta[i], t_ref, ta[i] - t_ref, pa[i], p_ref, pa[i] - p_ref)) print() print(" Over 500 studies: max |t diff| = %.3e, max |p diff| = %.3e" % (worst_t, worst_p)) print(" degrees of freedom used by both: %d" % int(dfa[0])) print(" VERDICT: %s" % ("agreement to better than 1e-9." if max(worst_t, worst_p) < 1e-9 else "DISAGREEMENT larger than 1e-9.")) else: print(" scipy unavailable; comparison skipped.") # ------------------------------------------------------------------ head("VALIDATION 2b. A t test small enough to check by hand") # ------------------------------------------------------------------ ha = np.array([[4.0, 7.0, 5.0, 6.0, 8.0]]) hb = np.array([[6.0, 9.0, 7.0, 10.0, 8.0]]) hy = np.concatenate([ha, hb], axis=1) hg = np.concatenate([np.zeros(5, int), np.ones(5, int)]) hm = np.ones((1, 10), dtype=bool) ht, hdf, hp = masked_ttest(hy, hg, hm, want_p=True) print(" group A: %s mean %.4f variance %.4f" % (ha[0].tolist(), ha.mean(), ha.var(ddof=1))) print(" group B: %s mean %.4f variance %.4f" % (hb[0].tolist(), hb.mean(), hb.var(ddof=1))) sp2 = (4 * ha.var(ddof=1) + 4 * hb.var(ddof=1)) / 8.0 print(" pooled variance = (4*%.4f + 4*%.4f)/8 = %.4f" % (ha.var(ddof=1), hb.var(ddof=1), sp2)) print(" standard error = sqrt(%.4f * (1/5 + 1/5)) = %.6f" % (sp2, np.sqrt(sp2 * 0.4))) print(" t = (%.1f - %.1f) / %.6f = %.6f" % (ha.mean(), hb.mean(), np.sqrt(sp2 * 0.4), ht[0])) print(" df = %d" % int(hdf[0])) print(" club p = %.9f" % hp[0]) if have_scipy: rt, rp = spstats.ttest_ind(ha[0], hb[0], equal_var=True) print(" scipy t, p = %.6f, %.9f" % (rt, rp)) print(" difference = %+.2e, %+.2e" % (ht[0] - rt, hp[0] - rp)) # ------------------------------------------------------------------ head("THE MAIN POOL") # ------------------------------------------------------------------ print("Simulating %d studies and evaluating all %d analysis paths in each." % (N_STUDIES, 5 * 3 * 3 * 2 * 2)) print("Every path is a defensible analysis of the same null data.") print() t0 = time.time() SIG, base_p = simulate_pool(ss_pool, N_STUDIES, progress=" ... {done} / {total} studies") print(" pool complete in %.1f s" % (time.time() - t0)) # ------------------------------------------------------------------ head("VALIDATION 3. The baseline. No researcher degrees of freedom.") # ------------------------------------------------------------------ k = int(np.sum(base_p < ALPHA)) n = len(base_p) fpr = k / n se = np.sqrt(fpr * (1 - fpr) / n) lo, hi = wilson(k, n) z = (fpr - ALPHA) / se print("One outcome measure, one test, at the planned sample size of %d per group." % N_MAX) print("Nothing is chosen after seeing the data. This is the check on everything") print("else in this file: if it does not come out at %.3f, nothing below means" % ALPHA) print("anything.") print() print(" significant studies : %d of %d" % (k, n)) print(" club false positive rate: %.5f" % fpr) print(" standard error : %.5f" % se) print(" 95%% Wilson interval : [%.5f, %.5f]" % (lo, hi)) print(" accepted (nominal) value: %.5f" % ALPHA) print(" difference : %+.5f" % (fpr - ALPHA)) print(" distance in SE : %+.2f sigma" % z) print(" interval covers 0.05 : %s" % ("yes" if lo <= ALPHA <= hi else "NO")) print() if abs(z) < 3: print(" VERDICT: baseline recovers the nominal level. Proceed.") else: print(" VERDICT: baseline is %.2f sigma from nominal. Something is wrong;" % z) print(" everything downstream is suspect until it is found.") # a second, distribution-level check that costs nothing print() print(" Further check on the same numbers: under the null the p values must") print(" be uniform on (0,1). Deciles of the baseline p values, expected %d each:" % (n // 10)) counts, _ = np.histogram(base_p, bins=10, range=(0.0, 1.0)) print(" ", " ".join("%5d" % c for c in counts)) chi2 = float(np.sum((counts - n / 10.0) ** 2) / (n / 10.0)) print(" chi-square on 9 df = %.2f (5%% critical value 16.92, 1%% is 21.67)" % chi2) # ------------------------------------------------------------------ head("VALIDATION 4. Optional stopping against Armitage, McPherson & Rowe (1969)") # ------------------------------------------------------------------ print("Equally spaced looks at accumulating data, %d new participants per group" % STOP_STEP) print("between looks, stop at the first significant result. This is the setting") print("the 1969 table describes, and it is NOT the peek schedule used in the main") print("grid, which adds only %d participants per group in total." % (LOOKS[-1] - LOOKS[0])) print("Studies: %d. Look-by-look critical values from the club's own t quantile." % N_STOP) print() t0 = time.time() curve = optional_stopping_curve(ss_stop, N_STOP) print(" %d studies in %.1f s" % (N_STOP, time.time() - t0)) print() print(" looks club rate SE 95% interval published diff sigma") for kk in sorted(ARMITAGE_PUBLISHED): r = curve[kk - 1] s = np.sqrt(r * (1 - r) / N_STOP) lo2, hi2 = wilson(int(round(r * N_STOP)), N_STOP) pub = ARMITAGE_PUBLISHED[kk] print(" %5d %.5f %.5f [%.5f,%.5f] %.3f %+.4f %+6.2f" % (kk, r, s, lo2, hi2, pub, r - pub, (r - pub) / s)) print() print(" Full curve, looks 1 to %d:" % STOP_MAX_K) for kk in range(1, STOP_MAX_K + 1): print(" K=%-3d %.5f" % (kk, curve[kk - 1])) print() worst = max(abs(curve[kk - 1] - ARMITAGE_PUBLISHED[kk]) for kk in ARMITAGE_PUBLISHED) print(" largest absolute deviation from the published table: %.4f" % worst) print(" NOTE: the published values are rounded to three decimal places and are") print(" the large-sample limit. At %d observations per group per look the" % STOP_STEP) print(" discrete t test is close to that limit but not identical to it, so an") print(" agreement of a few thousandths is the most this comparison can show.") # ------------------------------------------------------------------ head("VALIDATION 5. Reproducing Table 1 of Simmons, Nelson & Simonsohn (2011)") # ------------------------------------------------------------------ print("That paper is the reason this study exists. Its Table 1 gives false") print("positive rates for four researcher degrees of freedom, from 15,000") print("simulations each, at three significance levels. Three of the four can be") print("built out of the pieces already in this file, so we built them and ran") print("%d simulations of each. Situation D needs a three-condition design," % N_SIMMONS) print("which our model does not have, so it is not attempted.") print() print("Every operationalisation is taken from the note under their table, not") print("from our own reading of the text. Gender is a fair coin per observation.") print() t0 = time.time() sim_rates = simmons_replication(ss_simmons, N_SIMMONS) print(" %d simulations per row in %.1f s" % (N_SIMMONS, time.time() - t0)) print() print("Their table rests on 15,000 simulations per cell, so it carries Monte") print("Carlo error of its own. The sigma column below is the difference over") print("the standard error of the DIFFERENCE, which combines both runs.") print() print(" row situation alpha club " "published diff their SE sigma") worst_sigma = 0.0 worst_row = "" for rname in ("A", "B", "C", "AB", "ABC"): for ai, a in enumerate(SIMMONS_ALPHAS): club = sim_rates[rname][ai] pub = SIMMONS_PUBLISHED[rname][ai] sd_us = np.sqrt(club * (1 - club) / N_SIMMONS) sd_them = np.sqrt(pub * (1 - pub) / N_SIMMONS_PUB) sd = np.sqrt(sd_us * sd_us + sd_them * sd_them) zz = (club - pub) / sd if abs(zz) > worst_sigma: worst_sigma = abs(zz) worst_row = "%s at alpha %.2f" % (rname, a) print(" SIM %-3s %-38s %.2f %.5f %.3f %+.4f %.5f %+6.2f" % (rname if ai == 0 else "", SIMMONS_LABEL[rname] if ai == 0 else "", a, club, pub, club - pub, sd_them, zz)) print() print(" largest deviation: %.2f sigma, at %s." % (worst_sigma, worst_row)) print(" The published figures are also rounded to a tenth of a percentage") print(" point, worth up to %.4f on its own, so an agreement closer than that" % 0.0005) print(" cannot be demonstrated by this comparison however long we run.") if worst_sigma < 3.0: print(" VERDICT: our implementation reproduces the published table.") else: print(" VERDICT: one or more rows sit more than 3 sigma from the published") print(" value. We could not find an error in our code; the candidates are") print(" their Monte Carlo error, the rounding, and the coding ambiguity") print(" described below. Read the row above and judge for yourself.") print() print(" One ambiguity we could not resolve from the published description.") print(" In the model carrying a gender interaction, the condition coefficient") print(" depends on how gender is coded. With gender as a 0/1 dummy it is the") print(" condition effect among the gender coded 0; with gender centred it is") print(" the average condition effect. The rows above use the 0/1 dummy, which") print(" is what regression software does by default. Repeating Situation C") print(" with gender centred instead:") sim_c = simmons_replication(ss_simmons_c, N_SIMMONS // 2, centred_gender=True) for ai, a in enumerate(SIMMONS_ALPHAS): print(" SIMC alpha %.2f centred %.5f dummy %.5f published %.3f" % (a, sim_c["C"][ai], sim_rates["C"][ai], SIMMONS_PUBLISHED["C"][ai])) print() print(" Replicating a paper about analytic flexibility required us to make an") print(" analytic choice the paper did not pin down. We report both.") # ------------------------------------------------------------------ head("THE GRID. All 32 combinations of five researcher degrees of freedom") # ------------------------------------------------------------------ print("S = optional stopping, M = three outcome measures, G = subgroup on an") print("irrelevant split, C = covariate inclusion, O = outlier exclusion at 2 SD.") print("Every combination takes the smallest p value among the analyses it permits.") print("k_eff is the number of independent 5% tests that would give the same rate.") print() combos = [] for bits in range(32): combo = "".join(FREEDOMS[i] for i in range(5) if bits & (1 << i)) combos.append(combo) combos.sort(key=lambda c: (len(c), c)) results = {} print(" combo freedoms paths k_sig FPR " "SE 95% interval k_eff SE") for combo in combos: mp = combo_hits(SIG, combo) ks = int(np.sum(mp)) f = ks / N_STUDIES s = np.sqrt(f * (1 - f) / N_STUDIES) lo3, hi3 = wilson(ks, N_STUDIES) ke = k_effective(f) kse = k_eff_se(f, s) np_ = n_paths(combo) results[combo] = dict(k=ks, fpr=f, se=s, lo=lo3, hi=hi3, keff=ke, keff_se=kse, paths=np_) label = combo if combo else "(none)" names = "+".join(combo) if combo else "baseline" print(" %-8s %-32s %5d %6d %.5f %.5f [%.5f,%.5f] %6.2f %.2f" % (label, names, np_, ks, f, s, lo3, hi3, ke, kse)) print() print(" ONE AT A TIME, in order of damage:") singles = sorted([c for c in combos if len(c) == 1], key=lambda c: -results[c]["fpr"]) for c in singles: r = results[c] print(" %-2s %-34s %.5f +/- %.5f (%.1fx nominal, k_eff %.2f)" % (c, FREEDOM_NAME[c], r["fpr"], r["se"], r["fpr"] / ALPHA, r["keff"])) print() print(" CUMULATIVE, adopting freedoms in the order S, M, G, C, O:") order = "SMGCO" prev = results[""]["fpr"] for i in range(6): c = "".join(sorted(order[:i], key=lambda ch: FREEDOMS.index(ch))) r = results[c] print(" after %d freedom(s) %-6s paths %3d FPR %.5f +/- %.5f " "step %+.5f k_eff %.2f" % (i, c if c else "none", r["paths"], r["fpr"], r["se"], r["fpr"] - prev, r["keff"])) prev = r["fpr"] print() print(" The order matters for the intermediate rows and not for the last one.") print(" Adopting in the reverse order O, C, G, M, S:") order2 = "OCGMS" prev = results[""]["fpr"] for i in range(6): c = "".join(sorted(order2[:i], key=lambda ch: FREEDOMS.index(ch))) r = results[c] print(" after %d freedom(s) %-6s paths %3d FPR %.5f step %+.5f" % (i, c if c else "none", r["paths"], r["fpr"], r["fpr"] - prev)) prev = r["fpr"] full = results["SMGCO"] print() print(" HEADLINE: with all five freedoms in play, %d of %d studies of pure" % (full["k"], N_STUDIES)) print(" noise produced a significant result. That is %.4f, or %.1f%%, against a" % (full["fpr"], 100 * full["fpr"])) print(" nominal %.0f%%. The inflation factor is %.2f. The 180 analysis paths" % (100 * ALPHA, full["fpr"] / ALPHA)) print(" behave like %.2f independent tests, which is %.1f%% of their nominal" % (full["keff"], 100 * full["keff"] / full["paths"])) print(" count: the paths overlap heavily because they share data.") # ------------------------------------------------------------------ head("EFFECTIVE INDEPENDENT TESTS: nominal paths against k_eff") # ------------------------------------------------------------------ print(" combo paths k_eff k_eff/paths redundancy") for combo in combos: r = results[combo] ratio = r["keff"] / r["paths"] if r["paths"] else 0.0 print(" %-8s %5d %6.2f %6.3f %.2f paths per effective test" % (combo if combo else "(none)", r["paths"], r["keff"], ratio, (r["paths"] / r["keff"]) if r["keff"] > 0 else float("nan"))) # ------------------------------------------------------------------ head("CONVERGENCE") # ------------------------------------------------------------------ print("Running estimate as studies accumulate, for the baseline and for the") print("full combination. Monte Carlo error falls as 1/sqrt(n) and nothing else") print("changes, so a trace that is still wandering at the right-hand edge would") print("mean the run was too short.") print() base_hits = combo_hits(SIG, "").astype(np.int64) full_hits = combo_hits(SIG, "SMGCO").astype(np.int64) three_hits = combo_hits(SIG, "SMG").astype(np.int64) cb = np.cumsum(base_hits) cf = np.cumsum(full_hits) ct = np.cumsum(three_hits) marks = sorted(set(np.unique(np.round(np.logspace( np.log10(100), np.log10(N_STUDIES), 60)).astype(int)).tolist())) print(" CONVERGENCE TABLE (trials, baseline, +/-2SE, S+M+G, all five)") for m in marks: b = cb[m - 1] / m f = cf[m - 1] / m t3 = ct[m - 1] / m sb = np.sqrt(max(b * (1 - b), 1e-12) / m) sf = np.sqrt(max(f * (1 - f), 1e-12) / m) print(" CONV %7d %.5f %.5f %.5f %.5f %.5f" % (m, b, b - 2 * sb, b + 2 * sb, t3, f)) print() print(" last 10000 studies only: baseline %.5f, all five %.5f" % (base_hits[-10000:].mean(), full_hits[-10000:].mean())) print(" first 10000 studies only: baseline %.5f, all five %.5f" % (base_hits[:10000].mean(), full_hits[:10000].mean())) # how deep did a false positive need the researcher to go print() print("DEPTH. For each study, the first point in the order S, M, G, C, O at") print("which it became significant. 'never' means no permitted path reached") print("p < 0.05 even with all five freedoms.") seq = ["", "S", "MS", "GMS", "CGMS", "CGMOS"] seq_norm = [] for c in seq: seq_norm.append("".join(sorted(c, key=lambda ch: FREEDOMS.index(ch)))) hit_mat = np.stack([combo_hits(SIG, c) for c in seq_norm], axis=1) first = np.full(N_STUDIES, -1) for i in range(hit_mat.shape[1] - 1, -1, -1): first = np.where(hit_mat[:, i], i, first) labels = ["0 (honest test)", "1 (peeking)", "2 (+ outcomes)", "3 (+ subgroups)", "4 (+ covariate)", "5 (+ outliers)"] print() for i, lab in enumerate(labels): c = int(np.sum(first == i)) print(" DEPTH %d %-18s %6d %6.2f%%" % (i, lab, c, 100.0 * c / N_STUDIES)) nev = int(np.sum(first == -1)) print(" DEPTH - never %6d %6.2f%%" % (nev, 100.0 * nev / N_STUDIES)) # ------------------------------------------------------------------ head("SENSITIVITY. Where a different modelling choice changes the answer") # ------------------------------------------------------------------ print("Each cell below is an independent run of %d studies with its own stream." % N_SENS) print() print("A. Correlation between the three outcome measures. Choosing among") print(" outcomes costs least when the outcomes agree with each other.") print() print(" r FPR(M alone) SE k_eff FPR(all five) SE") r_rows = [] for r_val, ss in ((0.0, ss_sens_r0), (0.3, ss_sens_r3), (0.5, ss_sens_r5), (0.8, ss_sens_r8)): Ps, _bp = simulate_pool(ss, N_SENS, r=r_val) fM = float(np.mean(combo_hits(Ps, "M"))) sM = np.sqrt(fM * (1 - fM) / N_SENS) fA = float(np.mean(combo_hits(Ps, "SMGCO"))) sA = np.sqrt(fA * (1 - fA) / N_SENS) r_rows.append((r_val, fM, sM, fA, sA)) print(" SENSR %.2f %.5f %.5f %6.2f %.5f %.5f" % (r_val, fM, sM, k_effective(fM), fA, sA)) del Ps print() print(" The main grid uses r = 0.5. At r = 0 the three measures are three") print(" independent tests; at r = 0.8 they are nearly one. The whole result") print(" moves with a parameter nobody measured.") print() print("B. The outlier threshold. 2 SD, 2.5 SD and 3 SD are all rules people") print(" defend in print.") print() print(" threshold mean dropped FPR(O alone) SE k_eff FPR(all five)") t_rows = [] for tk, ss in ((2.0, ss_sens_t20), (2.5, ss_sens_t25), (3.0, ss_sens_t30)): Ps, _bp = simulate_pool(ss, N_SENS, trim_k=tk) fO = float(np.mean(combo_hits(Ps, "O"))) sO = np.sqrt(fO * (1 - fO) / N_SENS) fA = float(np.mean(combo_hits(Ps, "SMGCO"))) # expected fraction of a normal sample beyond tk SD drop = 2.0 * (1.0 - 0.5 * (1.0 + _erf(tk / np.sqrt(2.0)))) t_rows.append((tk, drop, fO, sO, fA)) print(" SENST %.1f %5.2f%% %.5f %.5f %5.2f %.5f" % (tk, 100 * drop, fO, sO, k_effective(fO), fA)) del Ps print() print(" A stricter rule drops fewer points, so it moves the p value less, so") print(" it buys less. The 2 SD rule in the main grid is the most generous of") print(" the three and our headline number depends on that choice.") print() print("C. The peek schedule. The main grid peeks four extra times while adding") print(" only %d participants per group. Armitage's setting doubles, triples" % (LOOKS[-1] - LOOKS[0])) print(" and quadruples the sample instead.") print() print(" main grid, S alone, 5 looks at n = 20..40 : %.5f" % results["S"]["fpr"]) print(" equally spaced, 5 looks at n = 25..125 : %.5f" % curve[4]) print(" ratio : %.2f" % (curve[4] / results["S"]["fpr"])) print(" The same word, 'we peeked five times', covers both. They are not the") print(" same act and they do not cost the same.") # ------------------------------------------------------------------ head("SUMMARY TABLE FOR THE ARTICLE") # ------------------------------------------------------------------ print(" TAB combo paths k_sig FPR SE lo hi " "k_eff keff_se inflation") for combo in combos: r = results[combo] print(" TAB %-5s %5d %6d %.5f %.5f %.5f %.5f %6.2f %6.2f %6.2f" % (combo if combo else "-", r["paths"], r["k"], r["fpr"], r["se"], r["lo"], r["hi"], r["keff"], r["keff_se"], r["fpr"] / ALPHA)) print() print(" ARM K club published") for kk in range(1, STOP_MAX_K + 1): pub = ARMITAGE_PUBLISHED.get(kk) print(" ARM %-3d %.5f %s" % (kk, curve[kk - 1], ("%.3f" % pub) if pub else "-")) print() elapsed = time.time() - t_start print("=" * 78) print("Total wall clock: %.1f s (%.1f min)" % (elapsed, elapsed / 60.0)) print("Master seed %d. Rerunning this file reproduces every number above." % SEED) print("=" * 78) def _erf(x): """Abramowitz & Stegun 7.1.26, good to 1.5e-7; used only for a printed expected-fraction column, never inside a test.""" s = np.sign(x) x = np.abs(x) t = 1.0 / (1.0 + 0.3275911 * x) y = 1.0 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t * np.exp(-x * x) return s * y if __name__ == "__main__": main()