9. Hypothesis tests¶
A new drug seems to add sleep, a coin seems to favour heads, a machine seems to overfill its bags. Each time the data show some effect, and each time the same doubt follows: could chance alone have produced a pattern this large? A hypothesis test is the standard way to answer that doubt with a number. In this chapter you will learn the logic of a test, what a p-value is and (just as important) what it is not, the two kinds of error a test can make and how to plan a study so that it has a fair chance of finding a real effect. You will then run the workhorse tests of science, the one-sample, paired and two-sample t tests, on a real dataset from 1905, see that every test is a confidence interval in disguise, and finish with the trap that catches whole fields: running many tests and reporting the lucky ones.
In this chapter
You will learn:
- the logic of a test: a null hypothesis, an alternative, a test statistic and a rule
- what a p-value is, why it is uniform when the null is true, and the six principles of the American Statistical Association's 2016 statement on p-values
- type I and type II errors, power, and how to choose a sample size before collecting data
- the one-sample, paired and two-sample t tests (pooled and Welch), with Student's own sleep data
- why a test and a confidence interval are two views of one object (with proof)
- multiple testing: the Bonferroni correction (proof), the false discovery rate and the Benjamini-Hochberg rule, and how p-hacking manufactures false results
You need: Chapter 6 (sampling distributions and the standard error), Chapter 8 (confidence intervals and Student's t), and the normal law from the Probability book (Chapter 10). From this chapter on, Part III uses some calculus (here only in 9.4, and lightly). A gentler version is the Topics lesson on hypothesis tests.
Levels: 9.1 and 9.2 are intermediate; 9.3 is intermediate with an advanced ending; 9.4 is intermediate, with an advanced proof; 9.5 is advanced; 9.6 is intermediate, with an advanced last part.
9.1 The logic: null and alternative¶
A question about a machine¶
A machine fills bags of flour labelled 500 grams. (This example is made up.) Long experience says that each bag's weight varies around the machine's setting with a standard deviation of \(\sigma = 4\) grams, and the factory wants the setting to be exactly 500. An inspector weighs 16 bags and finds an average of \(\bar x = 502.1\) grams.
Is the machine overfilling? The average is above 500, but averages wobble. If the setting were exactly 500, the average of 16 bags would still not be exactly 500; it would land somewhere near it. The real question is: how surprising is 502.1 if the setting is 500?
That reversal is the whole idea of a test. We cannot directly measure how likely it is that the machine is overfilling. What we can calculate is how likely data like ours would be if it were not. If such data would be very rare in that case, we have evidence against it.
The pieces of a test¶
A hypothesis is a statement about the population (or the process) that produced the data, usually about a parameter such as a mean \(\mu\).
- The null hypothesis \(H_0\) is the "nothing is going on" statement: the setting is 500, the drug does nothing, the coin is fair. It is the statement we try to find evidence against. The word "null" comes from the idea of zero effect.
- The alternative hypothesis \(H_1\) (also written \(H_A\)) is what we will believe instead if the evidence against \(H_0\) is strong. It can be two-sided (\(\mu \ne 500\): the machine is wrong in either direction) or one-sided (\(\mu > 500\): we only care about overfilling).
- A test statistic is a number computed from the data that measures how far the data are from what \(H_0\) predicts, in a unit where we know its distribution under \(H_0\).
- A rejection rule says which values of the statistic count as "too surprising to believe \(H_0\)".
For the flour, \(H_0: \mu = 500\) and, since both overfilling and underfilling are problems, \(H_1: \mu \ne 500\). The sample mean of \(n\) independent bags has standard error \(\sigma/\sqrt n\) (Chapter 6), here \(4/\sqrt{16} = 1\) gram. Bag weights are close to normal, so if \(H_0\) is true then \(\bar X \sim N(500, 1^2)\) and the z statistic
has the standard normal law \(N(0,1)\). Here \(\mu_0 = 500\) is the value named by the null hypothesis. With our data, \(z = (502.1 - 500)/1 = 2.1\). The average sits 2.1 standard errors above the null value.
Under \(N(0,1)\), values at least 2.1 away from 0 in either direction have probability \(2(1 - \Phi(2.1)) \approx 0.036\), where \(\Phi\) is the standard normal cumulative distribution function. So if the machine were set correctly, only about 3.6 percent of inspections would produce an average as far from 500 as this one. That is either a fairly rare accident or a sign that the setting has drifted.
The courtroom picture¶
A test works like a criminal trial. The defendant (\(H_0\)) is presumed innocent. The prosecution (the data) must produce evidence strong enough to convict "beyond reasonable doubt". Three features carry over exactly:
- The burden of proof is on the alternative. We never set out to prove \(H_0\).
- "Not guilty" is not "innocent". If the evidence is weak, we fail to reject \(H_0\). That does not show \(H_0\) is true; the effect might be real but too small for our sample to see. A study of 5 patients that finds no effect has shown very little.
- The standard of proof is fixed in advance. The significance level \(\alpha\) (often 0.05) is the risk of convicting an innocent defendant that we agree to accept. It is chosen before looking at the data.
The oldest test¶
The first published argument with this shape is usually credited to John Arbuthnot, physician to Queen Anne. In a 1710 paper in the Philosophical Transactions of the Royal Society he looked at the London christening records for 82 consecutive years and noted that in every one of them more boys than girls were christened. If each year were equally likely to favour boys or girls, like a fair coin, the chance that all 82 years favour boys would be
Arbuthnot concluded that chance was not the cause (he drew a theological conclusion from it, which is not a statistical one). In modern words: \(H_0\) was "each year is a fair coin", the statistic was the number of boy-majority years, and the data were absurdly unlikely under \(H_0\). This is also the first sign test, which returns in Chapter 15.
The recipe¶
Every test in this book follows the same five steps.
- State \(H_0\) and \(H_1\) in terms of a parameter.
- Choose the significance level \(\alpha\).
- Choose a test statistic whose distribution under \(H_0\) is known (exactly or approximately).
- Compute the statistic from the data, and the p-value (9.2).
- Reject \(H_0\) if the p-value is at most \(\alpha\); otherwise fail to reject. Then say what this means in the language of the problem, with an estimate of the effect and its uncertainty.
The last half of step 5 matters most and is skipped most often. "We rejected \(H_0\)" is never the end of a report.
9.2 p-values¶
The definition¶
Definition. The p-value of a test is the probability, computed assuming \(H_0\) is true, of getting a test statistic at least as extreme as the one actually observed. "Extreme" means "in the direction of \(H_1\)".
For a two-sided z test with observed \(z\), the p-value is \(P(|Z| \ge |z|) = 2(1 - \Phi(|z|))\). For a one-sided test of \(H_1: \mu > \mu_0\), it is \(P(Z \ge z) = 1 - \Phi(z)\).
A small p-value says: if \(H_0\) were true, data like these would be rare. So either a rare event happened, or \(H_0\) is false.
Example 1 (a coin). A coin lands heads 60 times in 100 tosses. Test \(H_0: p = \tfrac12\) against \(H_1: p \ne \tfrac12\), where \(p\) is the chance of heads.
Step 1: the exact p-value. Under \(H_0\) the number of heads is binomial with \(n = 100\) and \(p = \tfrac12\). By the symmetry of this law, the two-sided p-value is twice the upper tail:
Step 2: the normal approximation. Under \(H_0\) the count has mean 50 and standard deviation \(\sqrt{100 \cdot \tfrac12 \cdot \tfrac12} = 5\). The z statistic is \((60 - 50)/5 = 2\), with two-sided p-value \(2(1 - \Phi(2)) \approx 0.0455\).
Step 3: the continuity correction. The count is a whole number, so "at least 60" is better approximated by "above 59.5" (Probability, Chapter 10): \(z = (59.5 - 50)/5 = 1.9\), giving \(2(1 - \Phi(1.9)) \approx 0.0574\), very close to the exact 0.0569.
Check: the crude approximation says "significant at 5 percent" and the exact answer says "not quite". A result that flips between methods this easily is weak evidence either way, and that is the honest summary.
Try it in Python¶
from math import comb, erf, sqrt
def Phi(z): # standard normal cdf
return 0.5 * (1 + erf(z / sqrt(2)))
exact = 2 * sum(comb(100, k) for k in range(60, 101)) / 2**100
print(exact) # 0.0569
print(2 * (1 - Phi(2.0))) # 0.0455, no correction
print(2 * (1 - Phi(1.9))) # 0.0574, with the continuity correction
The p-value is uniform when the null is true¶
Here is the single most useful fact about p-values, and it explains why "reject when \(p \le \alpha\)" has error rate exactly \(\alpha\).
Theorem. Suppose the test statistic \(T\) has a continuous distribution under \(H_0\), and large values of \(T\) count as extreme. Let \(F\) be its cumulative distribution function under \(H_0\), and define the p-value as the random variable \(P = 1 - F(T)\). Then, if \(H_0\) is true, \(P\) is uniformly distributed on \([0, 1]\). In particular \(P(P \le \alpha) = \alpha\) for every \(\alpha\) between 0 and 1.
Proof. Assume for simplicity that \(F\) is strictly increasing where it matters (the general case uses the generalised inverse and goes the same way). For \(u \in (0, 1)\),
The fourth step uses that \(T\) has no atoms, so \(P(T \ge t) = 1 - F(t)\). A random variable with \(P(P \le u) = u\) for all \(u\) is uniform on \([0, 1]\). \(\square\)
Two consequences:
- The rule "reject when \(p \le \alpha\)" rejects a true null with probability exactly \(\alpha\). That is what the significance level promises.
- When nothing is going on, every p-value between 0 and 1 is equally likely. A p-value of 0.8 is not "evidence for \(H_0\)", and a p-value of 0.06 is not "almost significant evidence". Under \(H_0\), 0.06 is exactly as ordinary as 0.8. If 1000 researchers each test a true null, about 50 of them will get \(p \le 0.05\). This fact is the engine of the problems in 9.6.
For a discrete statistic, such as a binomial count, the p-value can only take certain values, and then \(P(P \le \alpha) \le \alpha\): the test is conservative, rejecting a true null at most (and usually less than) \(\alpha\) of the time.
What a p-value is not¶
The p-value is the most misread number in science. The common errors all confuse "the probability of the data, given the hypothesis" with "the probability of the hypothesis, given the data". These are different, as Bayes' theorem shows (Probability, Chapter 4).
- A p-value of 0.03 does not mean "there is a 3 percent chance that \(H_0\) is true". The p-value is computed assuming \(H_0\) is true, so it cannot also be the probability that \(H_0\) is true. To get that probability you need a prior, as in Chapter 14.
- It does not mean "there is a 3 percent chance the result is due to chance".
- A small p-value does not mean a large or important effect. With a million observations, a difference of 0.01 grams between two machines can have a tiny p-value and no practical meaning.
- A large p-value does not show that \(H_0\) is true. It may just mean the sample was small.
- \(p = 0.049\) and \(p = 0.051\) are practically the same evidence. Treating one as a discovery and the other as a failure is a habit, not a law.
The ASA statement of 2016¶
Worries about misuse grew so much that in 2016 the American Statistical Association (ASA) took the unusual step of publishing a formal statement on one specific matter of statistical practice: "The ASA's statement on p-values: context, process, and purpose", by Ronald Wasserstein and Nicole Lazar, in The American Statistician. It lists six principles, which in paraphrase say:
- P-values can indicate how incompatible the data are with a specified statistical model.
- P-values do not measure the probability that the studied hypothesis is true, or the probability that the data were produced by random chance alone.
- Scientific conclusions and business or policy decisions should not be based only on whether a p-value passes a specific threshold.
- Proper inference requires full reporting and transparency.
- A p-value, or statistical significance, does not measure the size of an effect or the importance of a result.
- By itself, a p-value does not provide a good measure of evidence regarding a model or hypothesis.
Note what principle 1 says carefully: "incompatible with a model". The p-value is computed under the whole model, including assumptions like independence and normality, not only under the claim \(\mu = \mu_0\). A small p-value may mean that \(\mu \ne \mu_0\), or that the observations were not independent, or that the data were chosen selectively. The test cannot tell which.
The statement did not ban p-values. A 2019 editorial in the same journal by Wasserstein, Allen Schirm and Lazar went further and urged dropping the phrase "statistically significant"; that editorial expressed its authors' views, not an official ASA position. In 2021 a task force appointed by the ASA's president published a separate statement saying that p-values and significance tests, properly applied and interpreted, remain valid tools. The working consensus that emerged: report the p-value as a number (not just "\(p < 0.05\)"), always alongside an estimate of the effect and a confidence interval, and never let a threshold do the thinking.
A short history¶
The p-value and the fixed-level test come from two rival schools that textbooks later glued together.
Ronald Fisher, in Statistical Methods for Research Workers (1925), treated the p-value as a continuous measure of evidence against a null hypothesis. He suggested 0.05 as a convenient point for judging whether a deviation is significant (it is about 2 standard deviations for a normal statistic), but he did not treat it as a rigid rule, and he had no notion of an alternative hypothesis or of power.
Jerzy Neyman and Egon Pearson, in papers culminating in "On the problem of the most efficient tests of statistical hypotheses" (1933), set things up as a decision problem: fix \(\alpha\) in advance, specify an alternative, and choose the test with the smallest chance of the other error (9.3). For them a single test result was not evidence about a hypothesis but one use of a procedure with known long-run error rates. Their main theorem, the Neyman-Pearson lemma, is proved in Chapter 13.
Fisher and Neyman disagreed bitterly for decades. Modern practice uses Fisher's p-value with Neyman and Pearson's \(\alpha\), power and alternatives, which works well enough if you remember which idea each piece came from.
9.3 Errors of type I and II, power¶
Two ways to be wrong¶
A test makes a binary decision about a binary truth, so there are four cases.
| \(H_0\) is true | \(H_0\) is false | |
|---|---|---|
| Reject \(H_0\) | type I error (false alarm), probability \(\alpha\) | correct, probability \(1 - \beta\) (the power) |
| Fail to reject | correct, probability \(1 - \alpha\) | type II error (a miss), probability \(\beta\) |
- A type I error is rejecting a true null: convicting the innocent, a smoke alarm going off with no fire, declaring a useless drug effective. Its probability is the significance level \(\alpha\).
- A type II error is failing to reject a false null: acquitting the guilty, a smoke alarm sleeping through a fire, missing a drug that works. Its probability is written \(\beta\).
- The power of a test is \(1 - \beta\): the probability of correctly rejecting \(H_0\) when it is false.
There is a tension. Making the alarm less sensitive (smaller \(\alpha\)) means more missed fires (larger \(\beta\)), for the same sample. The only way to reduce both at once is more information: a larger sample, or less noisy measurements.
Power is not one number: \(\beta\) depends on how false \(H_0\) is. A machine set at 500.1 grams is "not 500", but no reasonable inspection will notice. So we compute the power function: power as a function of the true parameter.
Power of the z test¶
Return to the flour. Test \(H_0: \mu = 500\) against \(H_1: \mu \ne 500\) at \(\alpha = 0.05\), with \(\sigma = 4\) and \(n = 16\), so the standard error is 1. We reject when \(|Z| \ge 1.96\). Suppose the true mean is 502. What is the chance we notice?
When \(\mu = 502\), the statistic \(Z = (\bar X - 500)/1\) is no longer standard normal: \(\bar X \sim N(502, 1)\), so \(Z \sim N(2, 1)\). The power is
So a drift of 2 grams, half a standard deviation of one bag, is caught only about half the time. The picture shows both laws, with \(\alpha\) and \(\beta\) as areas.
In general, for a two-sided z test at level \(\alpha\) with true mean \(\mu\), write \(\delta = \mu - \mu_0\) for the true effect and \(z_{\alpha/2}\) for the upper \(\alpha/2\) point of \(N(0,1)\) (1.96 when \(\alpha = 0.05\)). Then the same calculation gives
Read the formula: power depends on the effect only through \(\delta\sqrt n/\sigma\), the effect measured in standard errors. Doubling the effect has the same impact as quadrupling the sample. The ratio \(\delta/\sigma\) (effect divided by the spread of one observation) is called the standardised effect size; in the flour example it is \(2/4 = 0.5\).
| \(n\) | standard error | power at \(\mu = 502\) |
|---|---|---|
| 16 | 1.00 | 0.516 |
| 32 | 0.71 | 0.807 |
| 64 | 0.50 | 0.979 |
Choosing a sample size¶
Power analysis turns the formula around: decide the smallest effect worth detecting and the power you want, then solve for \(n\). When the effect is large compared with the noise, the second term in the power formula is tiny, so we need \(\Phi(\delta\sqrt n/\sigma - z_{\alpha/2}) = 1 - \beta\), that is \(\delta\sqrt n/\sigma - z_{\alpha/2} = z_\beta\). Solving,
Example 2. The factory wants an 80 percent chance of detecting a drift of 2 grams at \(\alpha = 0.05\). Then \(z_{0.025} = 1.960\), \(z_{0.2} = 0.842\), and
so weigh 32 bags. Check: the table above gives power 0.807 at \(n = 32\).
Power analysis should be done before collecting data. A study with low power has two problems, not one. It will usually miss real effects; and when it does find "significance", the estimated effect is usually exaggerated, because only the samples that happened to overshoot crossed the threshold. Andrew Gelman and John Carlin (2014) call this the "type M" (magnitude) error. Low power is one of the main causes of the replication problems discussed in Chapter 20.
Level: advanced. A post hoc power calculation, plugging the observed effect into the power formula after a non-significant result, adds nothing. For the z test the "observed power" is a fixed decreasing function of the p-value, so it repeats the p-value in other units. This point was made clearly by John Hoenig and Dennis Heisey (2001).
9.4 One- and two-sample t tests¶
When \(\sigma\) is unknown¶
The z test needs the population standard deviation \(\sigma\), which is almost never known. The natural fix is to replace it with the sample standard deviation \(s\), and this is exactly the problem William Sealy Gosset solved at the Guinness brewery in 1908, publishing as "Student" (Chapter 8, section 8.3).
Theorem (Student, 1908; proved rigorously by Fisher). If \(X_1, \dots, X_n\) are independent \(N(\mu, \sigma^2)\), then
has the t distribution with \(n - 1\) degrees of freedom, written \(t_{n-1}\). Here \(S^2 = \frac{1}{n-1}\sum (X_i - \bar X)^2\) is the sample variance.
Proof sketch. Level: advanced. Three facts combine.
- \(Z = \sqrt n(\bar X - \mu)/\sigma \sim N(0, 1)\), because a sum of independent normals is normal.
- \(V = (n-1)S^2/\sigma^2\) has the chi-squared distribution with \(n - 1\) degrees of freedom, written \(\chi^2_{n-1}\): the law of a sum of \(n - 1\) squares of independent standard normals.
- \(\bar X\) and \(S^2\) are independent.
Then \(T = Z / \sqrt{V/(n-1)}\), and the t distribution is defined as the law of \(Z/\sqrt{V/k}\) with \(Z \sim N(0,1)\) and \(V \sim \chi^2_k\) independent; its density follows by a change of variables and one integral (done in any of the books listed at the end).
Facts 2 and 3 come from one idea: rotation. Let \(Y_i = (X_i - \mu)/\sigma\), independent standard normals. Choose an orthogonal matrix \(Q\) (a rotation: its rows are perpendicular unit vectors) whose first row is \((1/\sqrt n, \dots, 1/\sqrt n)\). Set \(W = QY\). The joint density of independent standard normals, \((2\pi)^{-n/2} e^{-|y|^2/2}\), depends only on the length \(|y|\), and rotations preserve length and volume, so the \(W_i\) are again independent standard normals. Now \(W_1 = \sqrt n\,\bar Y\), and because rotation preserves length,
So \((n-1)S^2/\sigma^2\) is a sum of \(n - 1\) squared independent standard normals (fact 2), and it depends only on \(W_2, \dots, W_n\), which are independent of \(W_1\), which carries \(\bar X\) (fact 3). \(\square\)
The t density is bell-shaped and symmetric like the normal, but with heavier tails, because \(S\) sometimes underestimates \(\sigma\). As the degrees of freedom grow, \(S\) settles near \(\sigma\) and \(t_k\) approaches \(N(0,1)\). The two-sided 5 percent points show it:
| degrees of freedom | 1 | 2 | 5 | 9 | 20 | 30 | 100 | \(\infty\) (normal) |
|---|---|---|---|---|---|---|---|---|
| \(t_{0.025}\) | 12.71 | 4.30 | 2.571 | 2.262 | 2.086 | 2.042 | 1.984 | 1.960 |
The one-sample t test¶
To test \(H_0: \mu = \mu_0\) with normal data and unknown \(\sigma\), compute
and compare it with \(t_{n-1}\). The two-sided p-value is \(2\,P(T_{n-1} \ge |t|)\).
Robustness. The t test assumes normal data, but thanks to the central limit theorem (Chapter 6) it is close to correct for moderately large samples from most distributions, provided the observations are independent. Its real weaknesses are strong skewness with small \(n\), heavy outliers (which inflate \(s\) and pull \(\bar x\)), and above all dependence between observations, which no sample size cures. For small skewed samples, the methods of Chapter 15 are safer.
Student's sleep data and the paired test¶
In his 1908 paper, Gosset illustrated his method on data published in 1905 by the pharmacologists Arthur Cushny and
A. R. Peebles: the extra hours of sleep (compared with a night without a drug) of 10 patients, each of whom was
given each of two sleep-inducing drugs on different occasions. The data, as they appear in the dataset sleep that
ships with the statistics language R, are:
| patient | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| drug A | 0.7 | −1.6 | −0.2 | −1.2 | −0.1 | 3.4 | 3.7 | 0.8 | 0.0 | 2.0 |
| drug B | 1.9 | 0.8 | 1.1 | 0.1 | −0.1 | 4.4 | 5.5 | 1.6 | 4.6 | 3.4 |
| difference B − A | 1.2 | 2.4 | 1.3 | 1.3 | 0.0 | 1.0 | 1.8 | 0.8 | 4.6 | 1.4 |
(Historians have noted that the original labelling of the drugs in this table was muddled, so we call them A and B. The numbers are what matter here.)
Is drug B better than drug A? The same patients took both, so the natural unit of analysis is the difference for each patient. This removes the large differences between patients (patient 7 sleeps much more than patient 2 on either drug) and keeps only the comparison we care about. A test on the differences is called a paired t test; it is simply a one-sample t test on the differences, with \(H_0: \mu_d = 0\), where \(\mu_d\) is the mean difference in the population of patients.
Step 1: summaries. The 10 differences have mean \(\bar d = 15.8/10 = 1.58\) hours and standard deviation \(s_d \approx 1.230\).
Step 2: statistic. The standard error is \(1.230/\sqrt{10} \approx 0.389\), so
Step 3: p-value. \(2\,P(T_9 \ge 4.06) \approx 0.0028\).
Step 4: interval. With \(t_{0.025, 9} = 2.262\), a 95 percent confidence interval for \(\mu_d\) is \(1.58 \pm 2.262 \times 0.389\), that is, from 0.70 to 2.46 hours.
Step 5: in words. If the drugs had the same average effect, a mean difference this large would happen in fewer than 3 of every 1000 such experiments. Drug B added about 1.6 hours more sleep than drug A on average, and the data are consistent with any true advantage from about 0.7 to 2.5 hours. Check: 0 is outside the interval, matching \(p < 0.05\) (section 9.5 explains why these always agree).
Two independent samples¶
Often the two groups contain different subjects: a treatment group and a control group, chosen by randomisation (Chapter 5). Then there is no pairing. Suppose group 1 has \(n_1\) observations with mean \(\bar x_1\) and standard deviation \(s_1\), and group 2 likewise. We test \(H_0: \mu_1 = \mu_2\).
The difference \(\bar X_1 - \bar X_2\) has mean \(\mu_1 - \mu_2\) and, because the groups are independent,
There are two standard ways to estimate this.
Pooled t test. If we are willing to assume equal variances \(\sigma_1 = \sigma_2 = \sigma\), estimate \(\sigma^2\) by the pooled variance, a weighted average of the two sample variances:
Under \(H_0\) with normal data and equal variances, this has exactly a \(t_{n_1 + n_2 - 2}\) distribution (the same rotation argument, applied to each group, shows the numerator is normal and independent of a chi-squared with \(n_1 + n_2 - 2\) degrees of freedom).
Welch's t test. Without assuming equal variances, use
Its exact null distribution is not a t distribution (this is the "Behrens-Fisher problem", which has no exact solution of the simple kind). Bernard Welch (1947) showed it is very well approximated by \(t_\nu\), with the Welch-Satterthwaite degrees of freedom
which need not be a whole number. The idea behind it: match the mean and variance of the estimated variance to those of a scaled chi-squared variable. Many statisticians now use Welch's test by default: when the variances are equal it loses very little, and when they are not, the pooled test can be badly wrong, especially with unequal group sizes.
Example 3 (the wrong analysis). Suppose we forgot that the sleep data are paired and treated the two columns as independent groups of 10. Then \(\bar x_A = 0.75\), \(s_A \approx 1.789\), \(\bar x_B = 2.33\), \(s_B \approx 2.002\).
Step 1. \(s_A^2/10 + s_B^2/10 = 0.3201 + 0.4009 = 0.7210\), with square root \(0.849\).
Step 2. \(t = (2.33 - 0.75)/0.849 \approx 1.86\).
Step 3. The Welch degrees of freedom work out to \(\nu \approx 17.8\), and the two-sided p-value is about \(0.079\).
The pooled test gives the same \(t\) here (with equal group sizes the two formulas for the standard error agree, see exercise 11) on 18 degrees of freedom, and \(p \approx 0.079\).
Check the lesson: the same numbers give \(p = 0.003\) when analysed correctly and \(p = 0.08\) when the pairing is ignored. Ignoring the pairing throws away the fact that the patient-to-patient differences cancel, and the between-patient spread swamps the drug effect. The design of the study decides the test, not the other way round.
Try it in Python¶
import statistics as st
from math import sqrt
A = [0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0]
B = [1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4]
d = [b - a for a, b in zip(A, B)]
t_paired = st.mean(d) / (st.stdev(d) / sqrt(len(d)))
se = sqrt(st.variance(A) / 10 + st.variance(B) / 10)
t_welch = (st.mean(B) - st.mean(A)) / se
print(round(t_paired, 3), round(t_welch, 3)) # 4.062 and 1.861
To turn \(t\) into a p-value without a table you need the t distribution's cdf; the library scipy.stats has it, or
you can simulate it as in Chapter 6.
9.5 Tests and intervals as two views¶
Level: advanced.
In the sleep example, the 95 percent interval excluded 0 and the 5 percent test rejected \(\mu_d = 0\). This is not a coincidence. A confidence interval and a family of tests are the same object seen from two sides.
The duality theorem¶
Let \(\theta\) be a parameter. Suppose that for every possible value \(\theta_0\) we have a test of \(H_0: \theta = \theta_0\) at level \(\alpha\). Write \(A(\theta_0)\) for its acceptance region: the set of datasets for which this test does not reject.
Theorem (inverting tests).
- The set \(C(x) = \{\theta_0 : x \in A(\theta_0)\}\), all parameter values that the data do not reject, is a confidence set with confidence level at least \(1 - \alpha\).
- Conversely, if \(C(X)\) is a \(1 - \alpha\) confidence set, then the rule "reject \(H_0: \theta = \theta_0\) when \(\theta_0 \notin C(x)\)" is a test of level \(\alpha\).
Proof. The key is that the two events are the same event, by the definition of \(C\):
For part 1, suppose the true value is \(\theta_0\). The test of \(\theta_0\) has level \(\alpha\), so \(P_{\theta_0}(X \in A(\theta_0)) \ge 1 - \alpha\), and therefore \(P_{\theta_0}(\theta_0 \in C(X)) \ge 1 - \alpha\). Since this holds whatever the true \(\theta_0\) is, \(C\) covers the truth with probability at least \(1 - \alpha\).
For part 2, if \(H_0: \theta = \theta_0\) is true, the chance of rejecting is \(P_{\theta_0}(\theta_0 \notin C(X)) \le \alpha\), because \(C\) has coverage at least \(1 - \alpha\). \(\square\)
Example 4. The one-sample t test at level \(\alpha\) does not reject \(\mu_0\) when \(|\bar x - \mu_0| / (s/\sqrt n) < t_{\alpha/2, n-1}\). Rearranged, that is
the familiar t interval of Chapter 8. The set of values not rejected is the interval.
What duality buys you¶
- An interval says more than a test. The test says only whether one value, \(\mu_0\), is plausible. The interval lists all plausible values at once. "The drug adds between 0.7 and 2.5 hours" contains the test result (0 is not plausible) and also the size of the effect and its uncertainty. This is why the ASA statement and most journals now ask for intervals.
- Significance versus importance becomes visible. Picture four 95 percent intervals for the effect of a treatment where only effects above 1 unit matter in practice: \((0.1, 0.3)\) is significant but unimportant; \((2, 5)\) is significant and important; \((-0.2, 3)\) is not significant but cannot rule out an important effect (the study was too small); \((-0.1, 0.1)\) is not significant and does rule out any important effect. Two "non-significant" results can mean opposite things. Only the interval tells them apart.
- Equivalence testing. Sometimes the goal is to show that two treatments are practically the same (a cheaper generic drug against the brand). Failing to reject "no difference" does not show that. Instead, fix a margin \(\Delta\) of differences too small to matter, and conclude equivalence if a \(1 - 2\alpha\) confidence interval lies entirely inside \((-\Delta, \Delta)\). This is the "two one-sided tests" procedure, used by drug regulators.
- Caution with approximations. For proportions, the Wald interval and the score test are not exact duals of each other; the Wilson interval of Chapter 8 is exactly the inversion of the score test, which is one reason it behaves better.
9.6 Multiple testing and p-hacking¶
The problem of many tests¶
A test at level 0.05 has a 5 percent chance of a false alarm. Run 20 independent tests of true null hypotheses and the chance of at least one false alarm is
With 100 tests it is about 0.994. If you test 20 foods for a link with some illness, and none has any effect, you will probably still find "a significant link" with one of them. Report only that one, and you have published a false finding while every individual calculation was correct.
This is the multiple comparisons problem, and it appears everywhere: many outcomes measured in one trial, many subgroups (men, women, over 60, under 60, ...), many time points, many genes (Chapter 18 meets tests on 20 000 genes at once), many ways to analyse one dataset.
The family-wise error rate and Bonferroni¶
Suppose we test \(m\) null hypotheses \(H_1, \dots, H_m\). The family-wise error rate (FWER) is the probability of rejecting at least one true null hypothesis.
Theorem (Bonferroni correction). If each hypothesis is tested at level \(\alpha/m\) (reject \(H_i\) when \(p_i \le \alpha/m\)), then FWER \(\le \alpha\), whatever the dependence between the tests.
Proof. Let \(I_0\) be the set of true null hypotheses, with \(m_0 \le m\) members. A family-wise error is the union of the events "\(p_i \le \alpha/m\)" for \(i \in I_0\). The probability of a union is at most the sum of the probabilities (the union bound, also called Boole's inequality), and each true-null p-value satisfies \(P(p_i \le \alpha/m) \le \alpha/m\) (9.2). So
No independence was used, which is Bonferroni's great strength. Its weakness is power: with \(m = 20\,000\) genes and \(\alpha = 0.05\), each gene must reach \(p \le 2.5 \times 10^{-6}\), and real but modest effects are missed. The name honours Carlo Emilio Bonferroni, who proved inequalities of this kind in the 1930s; its use for multiple tests was popularised by Olive Jean Dunn around 1960.
A uniformly better method with the same guarantee is Holm's step-down procedure (Sture Holm, 1979): sort the p-values \(p_{(1)} \le \dots \le p_{(m)}\), compare \(p_{(1)}\) with \(\alpha/m\), \(p_{(2)}\) with \(\alpha/(m-1)\), and so on, rejecting until the first failure. It rejects everything Bonferroni rejects and sometimes more.
The false discovery rate¶
When we screen thousands of hypotheses, insisting on no false alarms at all is too strict. A biologist following up 100 candidate genes can live with 5 of them being duds, as long as she knows the proportion is small. This leads to a different target.
Let \(R\) be the number of rejections ("discoveries") and \(V\) the number of those that are false (true nulls rejected). The false discovery proportion is \(V/R\) (taken as 0 when \(R = 0\)), and the false discovery rate (FDR) is its expected value, \(E[V/R]\).
The Benjamini-Hochberg (BH) procedure (Yoav Benjamini and Yosef Hochberg, 1995) controls the FDR at level \(q\):
- Sort the p-values: \(p_{(1)} \le p_{(2)} \le \dots \le p_{(m)}\).
- Find the largest \(k\) with \(p_{(k)} \le \frac{k}{m} q\).
- Reject the hypotheses with the \(k\) smallest p-values (all of them, even any earlier ones that were above their own line).
Theorem (Benjamini and Hochberg, 1995). If the p-values of the true null hypotheses are independent of each other and of the false ones, the BH procedure has \(\text{FDR} = \frac{m_0}{m}q \le q\).
Idea of the proof. Level: advanced. Write \(V = \sum_{i \in I_0} \mathbf{1}[H_i \text{ rejected}]\), so \(E[V/R] = \sum_{i \in I_0} E\big[\mathbf{1}[H_i \text{ rejected}]/R\big]\). For one true null \(i\), split by the number of rejections \(R = k\): \(H_i\) is rejected with \(R = k\) exactly when \(p_i \le kq/m\) and the other \(m - 1\) p-values produce \(k - 1\) further rejections. Fixing the other p-values, that second event does not depend on \(p_i\) (once \(p_i\) is small enough to be rejected), so by independence and uniformity of \(p_i\) each term contributes \(\frac{1}{k}\cdot\frac{kq}{m}\) times the probability of the corresponding event for the others, and these probabilities add up to 1 over \(k\). Each true null therefore contributes exactly \(q/m\), and there are \(m_0\) of them. The full argument is a page; see Benjamini and Hochberg (1995). Benjamini and Daniel Yekutieli (2001) showed BH still controls the FDR under a common kind of positive dependence, and that dividing \(q\) by \(1 + \frac12 + \dots + \frac1m\) makes it safe under any dependence.
Example 5. Twelve hypotheses (made-up p-values) are tested with \(q = 0.05\). Sorted, with BH lines \(\frac{k}{12} \times 0.05\):
| \(k\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| \(p_{(k)}\) | 0.0004 | 0.003 | 0.009 | 0.011 | 0.019 | 0.04 | 0.07 | 0.15 | 0.33 | 0.62 | 0.71 | 0.94 |
| line \(0.05k/12\) | 0.0042 | 0.0083 | 0.0125 | 0.0167 | 0.0208 | 0.025 | 0.0292 | 0.0333 | 0.0375 | 0.0417 | 0.0458 | 0.05 |
| below line? | yes | yes | yes | yes | yes | no | no | no | no | no | no | no |
The largest \(k\) with \(p_{(k)}\) below its line is \(k = 5\), so BH rejects 5 hypotheses. Bonferroni (threshold \(0.05/12 \approx 0.0042\)) rejects only 2. Unadjusted testing at 0.05 would reject 6, including the one at 0.04.
Which correction to use depends on the cost of a false alarm. For a confirmatory drug trial with a few pre-named outcomes, control the FWER. For a screening study whose findings will be checked again, control the FDR.
p-hacking and the garden of forking paths¶
Multiple testing is not only about the tests a researcher reports. It is about the tests they could have run.
p-hacking is the practice, often innocent, of trying analyses until one gives \(p < 0.05\): measuring several outcomes and reporting the one that worked; adding subjects a few at a time and stopping as soon as \(p\) dips below 0.05; dropping "outliers" after seeing their effect; trying the analysis with and without a control variable. Each choice looks reasonable alone. Together they multiply the chances of a false positive.
In 2011, Joseph Simmons, Leif Nelson and Uri Simonsohn published "False-positive psychology" in Psychological Science. By simulation they showed that combining a few such common freedoms (two outcome measures, optional extra subjects, an optional covariate, dropping a condition) could push the real false-positive rate of a nominal 5 percent test far above 50 percent. They also ran a real (deliberately absurd) experiment that "showed" listening to a particular song made people younger, using only legal-looking flexibility, to make the point vivid.
The stopping-rule trap is easy to understand from 9.2. Under \(H_0\) the running p-value wanders randomly as data arrive; checking it after every new subject gives it many chances to dip below 0.05, and stopping at the first dip makes the dip certain to be reported. In fact, by the law of the iterated logarithm, if you check after every observation and are willing to wait forever, a true null will eventually be rejected with probability 1.
Andrew Gelman and Eric Loken coined the phrase garden of forking paths for the subtler version: a researcher who runs only one analysis, but would have chosen a different one had the data looked different, faces the same multiplicity, even without any conscious fishing.
The remedies are procedural more than mathematical, and they appear again in Chapter 20:
- Preregistration: write down the hypotheses, outcomes, sample size and analysis before seeing the data.
- Report everything: all outcomes, all analyses tried, all conditions (ASA principle 4).
- Correct for the tests you do run (Bonferroni, Holm, BH).
- Separate exploration from confirmation: an interesting pattern found while exploring is a hypothesis for the next study, not a conclusion of this one.
- Replicate.
Check yourself¶
1. A study reports \(p = 0.02\) for the effect of a fertiliser on yield. Which of these is correct? (a) There is a 2 percent chance the fertiliser has no effect. (b) If the fertiliser had no effect, results at least this extreme would occur about 2 percent of the time. (c) The fertiliser increases yield by a large amount.
Only (b). The p-value is a probability of data (at least this extreme) computed assuming the null is true; it is not the probability that the null is true, which (a) claims, and it says nothing about the size of the effect, which (c) claims. A large study can give \(p = 0.02\) for a tiny effect.
2. With \(\sigma = 10\) known, 25 measurements have mean 104. Test \(H_0: \mu = 100\) against \(\mu \ne 100\).
The standard error is \(10/\sqrt{25} = 2\), so \(z = (104 - 100)/2 = 2\). The two-sided p-value is \(2(1 - \Phi(2)) \approx 0.0455\). At \(\alpha = 0.05\) we reject (narrowly). The 95 percent interval \(104 \pm 1.96 \times 2 = (100.08, 107.92)\) just excludes 100, as duality requires.
3. A smoke alarm is a test of \(H_0\): there is no fire. Describe its type I and type II errors. A manufacturer makes it more sensitive. What happens to \(\alpha\) and \(\beta\)?
Type I: the alarm sounds with no fire (burnt toast). Type II: there is a fire and it stays silent. More sensitivity means it sounds more easily, so \(\alpha\) (false alarms) goes up and \(\beta\) (missed fires) goes down. For the same sensor you cannot lower both; a better sensor (more information) is the only way.
4. In exercise 2's setting (\(\sigma = 10\), \(n = 25\), two-sided, \(\alpha = 0.05\)), what is the power if the true mean is 105?
The effect is \(\delta = 5\), which is \(5/2 = 2.5\) standard errors. Power \(= \Phi(2.5 - 1.96) + \Phi(-2.5 - 1.96) = \Phi(0.54) + \Phi(-4.46) \approx 0.7054 + 0.0000 \approx 0.705\). So there is about a 30 percent chance of missing a true effect of 5.
5. How many measurements are needed for 90 percent power to detect an effect of 5 when \(\sigma = 10\), with a two-sided test at 5 percent?
Use \(n = ((z_{0.025} + z_{0.10})\sigma/\delta)^2\) with \(z_{0.025} = 1.960\) and \(z_{0.10} = 1.282\): \(n = (3.242 \times 10/5)^2 = 6.484^2 \approx 42.03\). Round up: 43 measurements.
6. Six made-up measurements are 12.1, 11.4, 12.8, 12.5, 13.2, 12.6. Test \(H_0: \mu = 12\) with a two-sided t test at 5 percent, and give the 95 percent interval.
The mean is \(74.6/6 \approx 12.433\); the deviations from it squared sum to about 1.933, so \(s = \sqrt{1.933/5} \approx 0.622\). The standard error is \(0.622/\sqrt 6 \approx 0.254\) and \(t = 0.433/0.254 \approx 1.71\) on 5 degrees of freedom. Since \(t_{0.025,5} = 2.571\) and \(1.71 < 2.571\), we do not reject; the p-value is about 0.15. The interval \(12.433 \pm 2.571 \times 0.254\) is about \((11.78, 13.09)\), which contains 12. With only six measurements the data are consistent with a mean anywhere from about 11.8 to 13.1.
7. A lab runs 1000 independent experiments, every one on a treatment that in fact does nothing, and tests each at \(\alpha = 0.05\). About how many 'significant' results appear? If the journal publishes only those, what fraction of its published findings in this area are false?
Under each true null the p-value is uniform, so each test is significant with probability 0.05, and about \(1000 \times 0.05 = 50\) are significant. All 50 published findings are false: 100 percent. Significance at 5 percent says nothing by itself about the fraction of published results that are true; that depends on how many of the tested hypotheses were true to begin with, and on the power.
8. A 95 percent confidence interval for a mean is \((1.2, 3.8)\). Without other information, what do two-sided tests at 5 percent conclude for \(H_0: \mu = 1\) and for \(H_0: \mu = 2\)?
By duality, the interval is the set of values not rejected at 5 percent. Since 1 is outside it, \(H_0: \mu = 1\) is rejected; since 2 is inside it, \(H_0: \mu = 2\) is not rejected.
9. A trial measures 8 outcomes. (a) If all nulls are true and the tests are independent, what is the chance of at least one \(p \le 0.05\)? (b) What is the Bonferroni threshold for a family-wise rate of 0.05?
(a) \(1 - 0.95^8 \approx 0.337\), about one chance in three. (b) \(0.05/8 = 0.00625\): reject an outcome's null only if its p-value is at most 0.00625.
10. Apply Bonferroni and Benjamini-Hochberg (both at 0.05) to the five p-values 0.001, 0.012, 0.014, 0.030, 0.20.
Bonferroni: threshold \(0.05/5 = 0.01\), so only 0.001 is rejected. BH: the lines are \(0.01, 0.02, 0.03, 0.04, 0.05\) for \(k = 1, \dots, 5\). Compare: \(0.001 \le 0.01\), \(0.012 \le 0.02\), \(0.014 \le 0.03\), \(0.030 \le 0.04\), but \(0.20 > 0.05\). The largest \(k\) below its line is 4, so BH rejects the four smallest. Holm, for comparison, checks \(0.001 \le 0.01\) (reject), then \(0.012 \le 0.05/4 = 0.0125\) (reject), then \(0.014 \le 0.05/3 \approx 0.0167\) (reject), then \(0.030 > 0.05/2 = 0.025\) (stop). Holm rejects three.
11. Show that when \(n_1 = n_2 = n\), the pooled and Welch t statistics are equal.
With equal sizes, \(s_p^2 = \frac{(n-1)(s_1^2 + s_2^2)}{2n - 2} = \frac{s_1^2 + s_2^2}{2}\), so the pooled standard error is \(s_p\sqrt{2/n} = \sqrt{\frac{s_1^2 + s_2^2}{2}\cdot\frac{2}{n}} = \sqrt{\frac{s_1^2}{n} + \frac{s_2^2}{n}}\), which is Welch's. The statistics are equal; only the degrees of freedom used for the p-value differ (\(2n - 2\) against Welch's \(\nu\), which is at most \(2n - 2\)).
12. (Harder.) For a one-sided z test of \(H_0: \mu = \mu_0\) against \(\mu > \mu_0\), the p-value is \(1 - \Phi(Z)\) with \(Z \sim N(\theta, 1)\), where \(\theta = \sqrt n(\mu - \mu_0)/\sigma\). Find \(P(p \le u)\) as a function of \(u\) and \(\theta\), check that it is \(u\) when \(\theta = 0\), and evaluate the power at \(u = 0.05\) when \(\theta = 2\).
\(p \le u\) means \(1 - \Phi(Z) \le u\), that is \(Z \ge z_u\), where \(z_u = \Phi^{-1}(1 - u)\). Since \(Z - \theta\) is standard normal, \(P(p \le u) = P(Z - \theta \ge z_u - \theta) = 1 - \Phi(z_u - \theta)\). When \(\theta = 0\) this is \(1 - \Phi(z_u) = u\): uniform, as proved in 9.2. With \(u = 0.05\), \(z_{0.05} = 1.645\) and \(\theta = 2\), the power is \(1 - \Phi(-0.355) = \Phi(0.355) \approx 0.639\). Under the alternative the p-values pile up near 0 instead of being spread evenly.
What you learned¶
- A test asks how surprising the data would be if the null hypothesis were true; evidence against \(H_0\) is not proof of \(H_1\), and failing to reject is not proof of \(H_0\).
- The p-value is the probability, under \(H_0\), of a result at least as extreme as the one seen; for a continuous statistic it is uniform under \(H_0\), so the rule "reject when \(p \le \alpha\)" has false-alarm rate \(\alpha\).
- The ASA's 2016 statement: p-values measure incompatibility with a model, not the probability of a hypothesis or the size of an effect; decisions should not rest on a threshold alone.
- Type I error (false alarm, \(\alpha\)), type II error (miss, \(\beta\)), power \(1 - \beta\); power depends on the effect in standard errors, \(\delta\sqrt n/\sigma\), and \(n = ((z_{\alpha/2} + z_\beta)\sigma/\delta)^2\) plans a sample size.
- t tests: one-sample, paired (a one-sample test on differences), pooled and Welch two-sample; the design decides the test, as the sleep data showed (\(p = 0.003\) paired against \(0.08\) unpaired).
- Tests and confidence intervals are duals: the interval is the set of null values the test does not reject, and it shows the size of the effect too.
- Many tests produce false alarms: Bonferroni controls the family-wise error under any dependence, Benjamini-Hochberg controls the false discovery rate, and p-hacking is multiple testing in disguise.
Next: Chapter 10, tests for counts in categories: chi-square tests, Mendel's peas, Fisher's exact test, and odds ratios.
Further reading¶
- Student (W. S. Gosset), "The probable error of a mean", Biometrika 6, 1908.
- R. A. Fisher, Statistical Methods for Research Workers, 1925 (many later editions).
- J. Neyman and E. S. Pearson, "On the problem of the most efficient tests of statistical hypotheses", Philosophical Transactions of the Royal Society A 231, 1933.
- B. L. Welch, "The generalization of 'Student's' problem when several different population variances are involved", Biometrika 34, 1947.
- Y. Benjamini and Y. Hochberg, "Controlling the false discovery rate: a practical and powerful approach to multiple testing", Journal of the Royal Statistical Society B 57, 1995.
- Y. Benjamini and D. Yekutieli, "The control of the false discovery rate in multiple testing under dependency", Annals of Statistics 29, 2001.
- S. Holm, "A simple sequentially rejective multiple test procedure", Scandinavian Journal of Statistics 6, 1979.
- R. L. Wasserstein and N. A. Lazar, "The ASA's statement on p-values: context, process, and purpose", The American Statistician 70, 2016.
- J. P. Simmons, L. D. Nelson and U. Simonsohn, "False-positive psychology", Psychological Science 22, 2011.
- J. M. Hoenig and D. M. Heisey, "The abuse of power: the pervasive fallacy of power calculations for data analysis", The American Statistician 55, 2001.
- A. Gelman and J. Carlin, "Beyond power calculations: assessing type S (sign) and type M (magnitude) errors", Perspectives on Psychological Science 9, 2014.
- E. L. Lehmann and J. P. Romano, Testing Statistical Hypotheses, 3rd ed., 2005 (the full theory).
- G. Casella and R. L. Berger, Statistical Inference, 2nd ed., 2002 (chapters 8 and 9: tests and their inversion into intervals).
- On this site: the Topics lessons on hypothesis tests and margin of error and confidence intervals; the Probability book on the normal distribution and Bayes' theorem; this book's index.