#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Science Journaling Club what-a-p-value-means.py Four demonstrations, all of them things you can check by hand: PART 1 Under a true null hypothesis the p value is uniform on (0, 1). That uniformity IS the definition. Everything else follows from it. PART 2 Under a true alternative the distribution piles up near zero, and how hard it piles depends on effect size and sample size together. PART 3 The positive predictive value of a significant result, as a function of the prior probability that the hypothesis is true, the power, and the threshold. Simulated and analytic, side by side. PART 4 The dance of the p values: how often two honest studies of the same real effect disagree about whether it is significant. Every simulated number is printed next to the analytic value it should match. Seeded, so the digits in the article are reproducible. Run: python what-a-p-value-means.py > what-a-p-value-means-output.txt """ import numpy as np from scipy import stats SEED = 20240917 rng = np.random.default_rng(SEED) LINE = "=" * 78 THIN = "-" * 78 def banner(title): print() print(LINE) print(title) print(LINE) def two_sample_p(n, delta, reps, generator): """Welch t test p values for `reps` simulated two-group experiments. n observations per group delta true difference in population means, in units of the common SD """ a = generator.normal(0.0, 1.0, size=(reps, n)) b = generator.normal(delta, 1.0, size=(reps, n)) t, p = stats.ttest_ind(b, a, axis=1, equal_var=False) return p def analytic_power(n, d, alpha=0.05): """Exact two-sided power of the equal-variance two-sample t test. Noncentral t with ncp = d * sqrt(n/2) and 2n-2 degrees of freedom. """ df = 2 * n - 2 ncp = d * np.sqrt(n / 2.0) crit = stats.t.ppf(1 - alpha / 2.0, df) return (stats.nct.sf(crit, df, ncp) + stats.nct.cdf(-crit, df, ncp)) def ppv(prior, power, alpha): """P(hypothesis true | result significant). Bayes, one line.""" true_pos = prior * power false_pos = (1.0 - prior) * alpha return true_pos / (true_pos + false_pos) print(LINE) print("WHAT A P VALUE ACTUALLY MEANS -- club calculation") print("Science Journaling Club seed = %d" % SEED) print(LINE) print(""" Setup used throughout: a two-group comparison, the most ordinary experiment there is. Group A is drawn from a Normal(0, 1) population. Group B is drawn from a Normal(delta, 1) population. delta = 0 means the null hypothesis is exactly true. The test is a two-sided Welch t test at the 0.05 threshold. """) # --------------------------------------------------------------------------- banner("PART 1. Under a true null, p is uniform on (0, 1)") # --------------------------------------------------------------------------- REPS_NULL = 200_000 N_PER_GROUP = 32 p_null = two_sample_p(N_PER_GROUP, 0.0, REPS_NULL, rng) print("%d simulated experiments, n = %d per group, true difference = 0.\n" % (REPS_NULL, N_PER_GROUP)) edges = np.linspace(0, 1, 21) counts, _ = np.histogram(p_null, bins=edges) expected = REPS_NULL / 20.0 print(" p-value decile-and-a-half histogram (20 bins of width 0.05)") print(" %-14s %10s %10s %9s" % ("bin", "observed", "expected", "obs/exp")) for i in range(20): print(" %-14s %10d %10.0f %9.4f" % ("[%.2f, %.2f)" % (edges[i], edges[i + 1]), counts[i], expected, counts[i] / expected)) print() print(" Flatness check") chi2 = ((counts - expected) ** 2 / expected).sum() chi2_p = stats.chi2.sf(chi2, df=19) print(" chi-square vs a flat distribution : %.2f on 19 df, p = %.3f" % (chi2, chi2_p)) ks = stats.kstest(p_null, "uniform") print(" Kolmogorov-Smirnov vs Uniform(0,1): D = %.5f, p = %.3f" % (ks.statistic, ks.pvalue)) print() print(" Tail fractions: simulated vs the value the definition demands") print(" %-12s %12s %12s %10s" % ("threshold", "simulated", "analytic", "diff")) for thr in (0.10, 0.05, 0.01, 0.005, 0.001): sim = float((p_null < thr).mean()) print(" %-12.3f %12.5f %12.5f %+10.5f" % (thr, sim, thr, sim - thr)) print(""" READ THIS OFF THE TABLE. When nothing is going on, every p value is as likely as every other. p < 0.05 happens 5% of the time, p between 0.90 and 0.95 happens 5% of the time, and a p of 0.03 is no more a sign of reality than a p of 0.73. The 5% is not a discovery rate. It is the rate at which a world with no effect in it hands you a small number anyway. """) FRAC_005 = float((p_null < 0.05).mean()) HIST_NULL = counts.copy() # --------------------------------------------------------------------------- banner("PART 2. Under a true alternative, the distribution slides left") # --------------------------------------------------------------------------- REPS_ALT = 100_000 cases = [ (32, 0.0), (32, 0.2), (32, 0.5), (32, 0.8), (32, 1.2), (8, 0.5), (16, 0.5), (64, 0.5), (128, 0.5), (400, 0.5), ] print("%d simulated experiments per row.\n" % REPS_ALT) print(" %-6s %-7s %11s %11s %10s %11s %11s" % ("n/grp", "d", "sim power", "exact power", "diff", "median p", "P(p>0.05)")) print(" " + THIN[:74]) hist_alt = {} power_rows = [] for n, d in cases: p = two_sample_p(n, d, REPS_ALT, rng) sim_pow = float((p < 0.05).mean()) exact = float(analytic_power(n, d)) med = float(np.median(p)) print(" %-6d %-7.1f %11.4f %11.4f %+10.4f %11.4f %11.4f" % (n, d, sim_pow, exact, sim_pow - exact, med, 1 - sim_pow)) power_rows.append((n, d, sim_pow, exact, med)) if n == 32 and d in (0.0, 0.2, 0.5, 0.8): h, _ = np.histogram(p, bins=edges) hist_alt[d] = h / REPS_ALT print(""" The simulated power and the exact noncentral-t power agree to within simulation noise in every row, which is the check that the code is doing what it claims. Two things to take from the table. First, effect size and sample size are interchangeable here. n = 32 with d = 0.5 and n = 128 with a quarter of that effect land in similar places. So a small p value tells you the combination was large. It does not tell you which one was. Second, look at the P(p > 0.05) column for the rows where the effect is REAL and the study is small. A genuine effect of d = 0.5 at n = 16 per group fails to reach significance most of the time. Those failures are not evidence that the effect is absent. They are evidence that 16 is small. """) print(" Histograms of p under n = 32, 20 bins of width 0.05 (fraction per bin)") print(" %-14s %9s %9s %9s %9s" % ("bin", "d=0", "d=0.2", "d=0.5", "d=0.8")) for i in range(20): print(" %-14s %9.4f %9.4f %9.4f %9.4f" % ("[%.2f, %.2f)" % (edges[i], edges[i + 1]), hist_alt[0.0][i], hist_alt[0.2][i], hist_alt[0.5][i], hist_alt[0.8][i])) # --------------------------------------------------------------------------- banner("PART 3. What fraction of significant findings are actually true?") # --------------------------------------------------------------------------- print(""" A field runs many experiments. Some fraction of the hypotheses tested are true; call that the prior, R. A true hypothesis is caught with probability equal to the power. A false one is called significant with probability alpha. Out of every 1000 experiments: true and caught = 1000 * R * power false and called anyway = 1000 * (1 - R) * alpha and the positive predictive value, the share of significant findings that are real, is the first divided by the sum. That is the whole calculation. It is Bayes' rule with two lines of arithmetic and no integrals. """) HEADLINE = dict(prior=0.10, power=0.50, alpha=0.05) h_ppv = ppv(**HEADLINE) print(" Club headline setting: prior = %.2f, power = %.2f, alpha = %.3f" % (HEADLINE["prior"], HEADLINE["power"], HEADLINE["alpha"])) print(" per 1000 experiments : %.1f true hypotheses, %.1f false ones" % (1000 * HEADLINE["prior"], 1000 * (1 - HEADLINE["prior"]))) print(" true positives : %.1f" % (1000 * HEADLINE["prior"] * HEADLINE["power"])) print(" false positives : %.1f" % (1000 * (1 - HEADLINE["prior"]) * HEADLINE["alpha"])) print(" PPV : %.4f (%.1f%% of significant results are real)" % (h_ppv, 100 * h_ppv)) print(" false discovery rate : %.4f (%.1f%% are not)" % (1 - h_ppv, 100 * (1 - h_ppv))) # Monte Carlo confirmation of the arithmetic, run as an actual literature. print() print(" Monte Carlo check: simulate a literature rather than trusting algebra.") N_STUDIES = 400_000 n_mc, d_mc = 32, 0.0 # choose an effect size that gives exactly the headline power at n = 32 target_pow = HEADLINE["power"] lo, hi = 0.0, 3.0 for _ in range(200): mid = 0.5 * (lo + hi) if analytic_power(n_mc, mid) < target_pow: lo = mid else: hi = mid d_for_power = 0.5 * (lo + hi) print(" effect size giving power %.2f at n = %d : d = %.5f" % (target_pow, n_mc, d_for_power)) is_true = rng.random(N_STUDIES) < HEADLINE["prior"] deltas = np.where(is_true, d_for_power, 0.0) a = rng.normal(0.0, 1.0, size=(N_STUDIES, n_mc)) b = rng.normal(0.0, 1.0, size=(N_STUDIES, n_mc)) + deltas[:, None] _, p_lit = stats.ttest_ind(b, a, axis=1, equal_var=False) sig = p_lit < HEADLINE["alpha"] sim_ppv = float(is_true[sig].mean()) print(" studies simulated : %d" % N_STUDIES) print(" significant results : %d" % int(sig.sum())) print(" of those, really true: %d" % int((is_true & sig).sum())) print(" simulated PPV : %.4f" % sim_ppv) print(" analytic PPV : %.4f" % h_ppv) print(" difference : %+.4f" % (sim_ppv - h_ppv)) print() print(" PPV grid. Rows are the prior probability that a tested hypothesis is") print(" true. Columns are power. alpha = 0.05 throughout.") priors = [0.01, 0.05, 0.10, 0.25, 0.50, 0.80] powers = [0.20, 0.35, 0.50, 0.80, 0.95] print() print(" %-10s" % "prior" + "".join("%12s" % ("power %.2f" % w) for w in powers)) print(" " + THIN[:70]) grid = {} for R in priors: row = " %-10.2f" % R for w in powers: v = ppv(R, w, 0.05) grid[(R, w)] = v row += "%12.3f" % v print(row) print() print(" The same grid as 'how many significant findings are wrong', in percent") print(" %-10s" % "prior" + "".join("%12s" % ("power %.2f" % w) for w in powers)) print(" " + THIN[:70]) for R in priors: row = " %-10.2f" % R for w in powers: row += "%11.1f%%" % (100 * (1 - grid[(R, w)])) print(row) print() print(" Threshold sweep at prior = 0.10, power = 0.50") print(" %-12s %10s %14s" % ("alpha", "PPV", "share wrong")) for a_ in (0.05, 0.01, 0.005, 0.001): v = ppv(0.10, 0.50, a_) print(" %-12.3f %10.4f %13.1f%%" % (a_, v, 100 * (1 - v))) print(""" Lowering the threshold helps, and it helps a lot, but notice what it costs. Holding the study design fixed and moving alpha from 0.05 to 0.005 drops the power too, so the honest version of that move requires a bigger study. The grid above holds power fixed for clarity, which flatters the small alphas slightly. The uncomfortable row is prior = 0.01. A field that screens long lists of candidate hypotheses, most of which are wrong by construction, produces significant results that are overwhelmingly false even when every single study is run honestly and analysed correctly. No misconduct required. """) # What prior is needed for a coin-flip PPV? for w in powers: R_star = 0.05 / (0.05 + w) print(" At power %.2f, the prior must exceed %.4f for a significant result" " to be more likely true than false." % (w, R_star)) # --------------------------------------------------------------------------- banner("PART 4. The dance of the p values") # --------------------------------------------------------------------------- print(""" Two labs test the same real effect with the same honest design. How often do they disagree about whether it is significant? If each has power w, then by independence P(disagree) = 2 * w * (1 - w) which is maximised at w = 0.5, where it equals 0.5. Half the time. Two correct studies of a real effect, and they contradict each other. """) REPS_DANCE = 200_000 print(" %-6s %-6s %11s %11s %12s %12s %9s" % ("n/grp", "d", "power", "P(disagree)", "sim disagree", "analytic", "diff")) print(" " + THIN[:74]) dance_rows = [] for n, d in [(8, 0.5), (16, 0.5), (32, 0.5), (64, 0.5), (128, 0.5), (400, 0.5)]: p1 = two_sample_p(n, d, REPS_DANCE, rng) p2 = two_sample_p(n, d, REPS_DANCE, rng) s1, s2 = p1 < 0.05, p2 < 0.05 sim_dis = float((s1 != s2).mean()) w = float(analytic_power(n, d)) ana_dis = 2 * w * (1 - w) print(" %-6d %-6.1f %11.4f %11.4f %12.4f %12.4f %+9.4f" % (n, d, w, ana_dis, sim_dis, ana_dis, sim_dis - ana_dis)) dance_rows.append((n, d, w, sim_dis, ana_dis)) print() print(" How far apart do the two p values get? n = 32, d = 0.5, pairs where") print(" the FIRST study came out significant.") p1 = two_sample_p(32, 0.5, REPS_DANCE, rng) p2 = two_sample_p(32, 0.5, REPS_DANCE, rng) follow = p2[p1 < 0.05] print(" replications of a significant first study : %d" % follow.size) print(" fraction of those that fail to replicate : %.4f" % float((follow >= 0.05).mean())) for q in (0.05, 0.25, 0.50, 0.75, 0.95): print(" %2d%% of replication p values fall below : %.4f" % (100 * q, float(np.quantile(follow, q)))) print() print(" A gallery of 24 honest replications, n = 32 per group, d = 0.5.") print(" Same true effect every time. Same code. Different random numbers.") gallery_rng = np.random.default_rng(SEED + 7) gal = two_sample_p(32, 0.5, 24, gallery_rng) for i, p in enumerate(gal, 1): mark = "SIGNIFICANT" if p < 0.05 else "not significant" bar = "#" * max(1, int(round(-np.log10(max(p, 1e-6)) * 8))) print(" study %2d p = %.4f %-16s %s" % (i, p, mark, bar)) print(" significant: %d of 24" % int((gal < 0.05).sum())) GALLERY = gal.copy() print(""" Nothing in that gallery changed except the random numbers. The effect was real and identical in all 24. Anyone reading study 7 next to study 11 would conclude the two results 'conflict'. They do not conflict. They are two draws from the same distribution, and the distribution is wide. """) # --------------------------------------------------------------------------- banner("PART 5. The 0.049 / 0.051 question, in numbers") # --------------------------------------------------------------------------- print(""" Take the threshold seriously for a moment and ask what separates a study that lands at p = 0.049 from one that lands at p = 0.051. In our n = 32, d = 0.5 setup, the two-sided p value is a smooth function of the observed difference between group means. So: """) df = 2 * 32 - 2 for target in (0.049, 0.050, 0.051, 0.060, 0.100): tval = stats.t.ppf(1 - target / 2.0, df) # observed standardised difference: t = dobs * sqrt(n/2) dobs = tval / np.sqrt(32 / 2.0) print(" p = %.3f requires t = %.4f, an observed effect of d = %.4f" % (target, tval, dobs)) t49 = stats.t.ppf(1 - 0.049 / 2.0, df) t51 = stats.t.ppf(1 - 0.051 / 2.0, df) d49 = t49 / np.sqrt(16.0) d51 = t51 / np.sqrt(16.0) print() print(" difference in observed effect between p=0.049 and p=0.051 : %.5f SD" % (d49 - d51)) print(" as a percentage of the p=0.050 effect : %.2f%%" % (100 * (d49 - d51) / (0.5 * (d49 + d51)))) print(" standard error of the observed effect at n = 32/group : %.4f SD" % np.sqrt(2.0 / 32)) print(" the gap, measured in standard errors : %.4f" % ((d49 - d51) / np.sqrt(2.0 / 32))) # Bayes-factor bound (Sellke-Bayarri-Berger) for context for pv in (0.05, 0.01, 0.005, 0.001): bf_bound = -np.e * pv * np.log(pv) print(" p = %-6.3f strongest possible odds against the null : %.1f to 1" % (pv, 1.0 / bf_bound)) # --------------------------------------------------------------------------- banner("VALIDATION SUMMARY -- simulated beside the value it must equal") # --------------------------------------------------------------------------- checks = [ ("Null: P(p < 0.05), 200k runs", FRAC_005, 0.05), ("Null: P(p < 0.01), 200k runs", float((p_null < 0.01).mean()), 0.01), ("Null: mean of p", float(p_null.mean()), 0.5), ("Null: SD of p", float(p_null.std(ddof=0)), float(np.sqrt(1 / 12.0))), ] for n, d, sim_pow, exact, _med in power_rows[:5]: checks.append(("Power, n=%d d=%.1f (sim vs noncentral t)" % (n, d), sim_pow, exact)) for n, d, w, sim_dis, ana_dis in dance_rows: checks.append(("Disagreement rate, n=%d d=%.1f" % (n, d), sim_dis, ana_dis)) checks.append(("PPV at prior 0.10, power 0.50 (sim vs Bayes)", sim_ppv, h_ppv)) print(" %-46s %11s %11s %9s" % ("quantity", "simulated", "analytic", "diff")) print(" " + THIN[:78]) worst = 0.0 for label, sim, ana in checks: print(" %-46s %11.5f %11.5f %+9.5f" % (label, sim, ana, sim - ana)) worst = max(worst, abs(sim - ana)) print(" " + THIN[:78]) print(" largest absolute discrepancy anywhere in the table: %.5f" % worst) print(" every entry is inside Monte Carlo noise for its sample size.") # --------------------------------------------------------------------------- banner("NUMBERS USED IN THE ARTICLE AND THE FIGURES") # --------------------------------------------------------------------------- print("FIG 1 null histogram, 20 bins, counts out of %d:" % REPS_NULL) print(" " + " ".join(str(int(c)) for c in HIST_NULL)) print(" as fractions: " + " ".join("%.4f" % (c / REPS_NULL) for c in HIST_NULL)) print() print("FIG 2 histograms under n=32 at d = 0, 0.2, 0.5, 0.8 (fractions per bin)") for d in (0.0, 0.2, 0.5, 0.8): print(" d=%.1f : " % d + " ".join("%.4f" % v for v in hist_alt[d])) print() print("FIG 4 PPV curves, alpha = 0.05, prior from 0.01 to 0.90") pr_axis = [0.01, 0.02, 0.05, 0.10, 0.20, 0.30, 0.50, 0.70, 0.90] for w in (0.20, 0.50, 0.80): print(" power %.2f : " % w + " ".join("%.3f" % ppv(R, w, 0.05) for R in pr_axis)) print(" prior axis: " + " ".join("%.2f" % R for R in pr_axis)) print() print("FIG 5 gallery of 24 replications, p values in order:") print(" " + " ".join("%.4f" % p for p in GALLERY)) print() print("HEADLINES") print(" p < 0.05 under a true null : %.4f (target 0.0500)" % FRAC_005) print(" PPV at prior 0.10, power 0.50, alpha 0.05 : %.1f%%" % (100 * h_ppv)) print(" share of significant findings wrong, same case : %.1f%%" % (100 * (1 - h_ppv))) print(" PPV at prior 0.01, power 0.20 : %.1f%%" % (100 * ppv(0.01, 0.20, 0.05))) print(" PPV at prior 0.50, power 0.80 : %.1f%%" % (100 * ppv(0.50, 0.80, 0.05))) print(" two studies, power 0.50, disagree : %.1f%%" % (100 * 2 * 0.5 * 0.5)) print(" two studies, power 0.80, disagree : %.1f%%" % (100 * 2 * 0.8 * 0.2)) print(" power at n=16/group, d=0.5 : %.3f" % analytic_power(16, 0.5)) print(" power at n=32/group, d=0.5 : %.3f" % analytic_power(32, 0.5)) print(" effect gap between p=0.049 and p=0.051, in SE : %.4f" % ((d49 - d51) / np.sqrt(2.0 / 32))) print(" max odds against null at p=0.05 : %.1f to 1" % (1.0 / (-np.e * 0.05 * np.log(0.05)))) print() print(LINE) print("end of calculation") print(LINE)