Statistics

The last chapter ended with a promise: a sample can tell you about the whole. The law of large numbers guarantees that the mean of many independent observations ends up close to the expectation, and the central limit theorem tells you exactly how it wobbles around it. What’s left is to learn to use these theorems without fooling ourselves.

That’s harder than it sounds. “The first principle is that you must not fool yourself, and you are the easiest person to fool,” Richard Feynman told Caltech graduates in 1974. Mark Twain attributed to the British prime minister Benjamin Disraeli a line about three kinds of lies: “lies, damned lies, and statistics”. It has never been found in Disraeli’s speeches or books, but the line lives on, because everyone has seen true numbers propping up false conclusions.

The chapter is set out like a myth-busting show. We have five newspaper headlines. They’re made up, but you’ve met each of them in dozens of versions, and the last two have real stories behind them. Behind each headline we’ll find the number it rests on, work out what that number really says, and borrow a tool from mathematical statistics to test the myth.

The hospital’s average temperature

At Acme Ltd the average salary is 145 thousand a month, and it rose by 10 % over the year

Let’s look at the payroll. Acme has ten people. Nine employees earn (in thousands) $40$, $42$, $45$, $48$, $50$, $52$, $55$, $58$ and $60$, and the director earns $1000$. The total is $1450$; divide by ten and you get exactly $145$. The headline isn’t lying. But nobody in the company earns anything like $145$ thousand: nine people earn at most $60$, and the director almost seven times the average.

Now the second half of the headline. Over the year the director raised his own salary to $1145$ thousand, and the others didn’t get a penny more. The total became $1595$ and the mean $159.5$, exactly $10$ % more than before. Again the headline is honest, though nine employees out of ten got nothing from this rise. In the same way, the average temperature of a hospital can be perfectly normal when half the patients have a fever and the other half have gone cold for good.

The arithmetic mean $\bar x = \frac{x_1 + \dots + x_n}{n}$ answers the question of how much everyone would get if everything were shared equally. The question “how much does an ordinary employee earn” is answered by another number.

It has little in common with the median of a triangle from chapter 17: that one halves the area of a triangle, this one a list of numbers. At Acme the median is the average of the fifth and sixth numbers, $(50 + 52)/2 = 51$. The director can raise his salary to a billion and the median won’t budge: the director was last on the list and stays last. Statisticians call this property robustness. One wild value can’t drag the median far, while it can drag the mean anywhere. There’s an old joke: when Bill Gates walks into a bar, the average customer becomes a billionaire.

That’s why it’s more sensible to ask for the median when it comes to incomes. National statistics offices, for example, publish both mean and median earnings, and the median is noticeably lower: the distribution of incomes, as at Acme, is stretched out to the right, and a few large salaries pull the mean up.

The mean and the median balance the data in different ways. Put identical weights on a light plank at the points $x_1, \dots, x_n$ and place a fulcrum at the point $c$. A weight at $x_i$ turns the plank with a force proportional to its lever arm $x_i - c$, and the plank tips in the direction given by the sign of the sum $\sum (x_i - c)$. This sum is zero exactly when $c = \bar x$, because $\sum (x_i - \bar x) = \sum x_i - n\bar x = 0$. The mean is the balance point, the centre of mass of the weights. The median, by contrast, is found by voting: the fulcrum is at the median when there are as many weights on the left as on the right, however far away they lie.

There’s another way to look at these two numbers. Suppose the whole set must be represented by a single number $c$, and for each $x_i$ we pay a penalty for missing. If the penalty is the distance $|x_i - c|$, the best choice is the median. If it grows like the square of the distance, the best choice is the mean. Let’s prove both claims.

The trick is to follow how the sum changes when $c$ moves a tiny bit.

Mark the numbers on a line, and a trial point $c$. The sum $S(c)$ is the total length of the segments from $c$ to all the numbers. Move $c$ to the right by a small $h$, without jumping over any number. Let $L$ numbers lie to the left of $c$ and $R$ to the right. Each segment to a number on the left gets longer by $h$, each segment to a number on the right gets shorter by $h$. So the sum changes by $(L - R)\,h$. While there are more numbers to the right of $c$ ($R > L$), moving right reduces the sum; when there are more on the left, moving left reduces it. It pays to move towards where there are more numbers. So the graph of $S(c)$ is a polygonal line, and on each of its segments the slope is $L - R$. To the left of all the numbers the slope is $-n$, and each time it passes a number it goes up by $2$: one number has moved from the right-hand ones to the left-hand ones. The slope is negative while $L < R$, and positive when $L > R$. The polygonal line takes its smallest value where the slope changes sign. For odd $n = 2k + 1$ that’s the point $x_{k+1}$: $k$ numbers to its left and $k$ to its right; it’s the median. For even $n = 2k$ the slope between $x_k$ and $x_{k+1}$ is zero, the line runs horizontally, and every point of that segment gives the minimum, including the midpoint, the median.

The trick is to see a set of $n$ numbers as a single point of $n$-dimensional space. For two numbers it’s the ordinary plane.

Two numbers $x_1, x_2$ are a point $X = (x_1, x_2)$, and the “one-number answer” $c$ is a point $C = (c, c)$ on the diagonal $y = x$. By the distance formula (chapter 18) $|XC|^2 = (x_1 - c)^2 + (x_2 - c)^2$, exactly our sum of squares. Drop a perpendicular from $X$ to the diagonal. Its foot is the point $M = (\bar x, \bar x)$. Let’s check: the vector $\overrightarrow{MX} = (x_1 - \bar x,\ x_2 - \bar x)$ has coordinate sum $x_1 + x_2 - 2\bar x = 0$, and that’s its dot product with the direction of the diagonal $(1, 1)$ (chapter 22). The product is zero, so the vectors are perpendicular. The triangle $XMC$ is right-angled, with the right angle at $M$. By Pythagoras’ theorem $|XC|^2 = |XM|^2 + |MC|^2$. The leg $|XM|^2 = (x_1 - \bar x)^2 + (x_2 - \bar x)^2$ doesn’t depend on $c$, and the other leg $|MC|^2 = 2(\bar x - c)^2$ vanishes only when $C = M$. The point of the diagonal nearest to $X$ is $M$, and the sum of squares is smallest at $c = \bar x$. For $n$ numbers the argument is the same, only in $n$-dimensional space, where Pythagoras is one line of algebra. Write $x_i - c = (x_i - \bar x) + (\bar x - c)$ and square: $(x_i - c)^2 = (x_i - \bar x)^2 + 2(x_i - \bar x)(\bar x - c) + (\bar x - c)^2$. Add over all $i$. The middle terms give $2(\bar x - c)\sum(x_i - \bar x) = 0$, which is the perpendicularity, and the last ones give $n(\bar x - c)^2$. We have the identity.

Which penalty to choose is decided by the problem, not by mathematics. But squares are easier to compute with: sums of squares have formulas like the one just proved, while sums of absolute values don’t. That’s why classical statistics rests on the mean and squared deviations, and a measure of spread is built from the same squares. It’s needed: an average salary of $50$ thousand can occur in a company where everyone earns between $48$ and $52$, and in one where half earn $10$ and half $90$.

The sum of squared deviations from the mean. Squaring removes the sign: deviations up and down don’t cancel each other, as they would in the sum of the deviations themselves, which is always zero. The square root at the end takes us back from “squared money” to money. We divide not by $n$ but by $n - 1$. Why will become clear in the second investigation; for large $n$ the difference is negligible. Example: for the nine Acme employees without the director the mean is $50$, the deviations are $-10, -8, -5, -2, 0, 2, 5, 8, 10$, the sum of their squares is $386$, and $s = \sqrt{386/8} \approx 6.9$ thousand. With the director $s \approx 300$ thousand, twice the mean itself. Such a spread is a sure sign that the mean alone says little here.

Find the sample standard deviation of the numbers $3$, $7$, $7$, $9$, $14$.

The mean is $\bar x = 40/5 = 8$. The deviations are $-5$, $-1$, $-1$, $1$, $6$ (their sum, as it should be, is zero). The sum of squares is $25 + 1 + 1 + 1 + 36 = 64$; divide by $n - 1 = 4$: $s^2 = 16$, $s = 4$.

Over a year the mean income of a city’s residents rose by $8$ %, while the median didn’t change. What can be said for certain?

Total income rose, and these numbers say nothing about most residents Most residents’ incomes rose The incomes of those richer than the median rose Everyone’s income rose by 8 %

The mean rose, so the total rose, and someone got a raise. The median stayed put, so the person in the middle of the list earns the same as whoever stood there a year ago. Who got the raise can’t be told from these two numbers.

Both numbers in the headline are correct. But the mean answers the question “how much would everyone get if it were shared equally”, while the reader hears “how much does an ordinary person earn”. That question is answered by the median, and with any mean you should also ask about the spread.

A thousand people for a whole city

Poll: 52 % of residents for building a new park, 48 % against. 1000 people surveyed

The first objection you hear to any poll: how can a thousand people speak for a city of millions? The second, the opposite: if a poll got it wrong, too few people were asked. Both are wrong, and the best illustration is a story with which many statistics courses begin.

In 1936 the USA was electing a president: the incumbent Franklin Roosevelt against the Republican Alf Landon. The Literary Digest magazine, which had correctly predicted the winners of several previous elections, sent out about ten million questionnaires and got more than two million replies. The magazine’s forecast: Landon would win about $57$ % of the vote. George Gallup surveyed dozens of times fewer people and predicted a Roosevelt win. Roosevelt got about $61$ % and won $46$ states out of $48$. Two years later the magazine ceased publication.

The Literary Digest’s trouble wasn’t the number of replies but whose replies they were. The addresses came from telephone directories, lists of car owners and the magazine’s subscribers: at the height of the Great Depression these were people noticeably better off than average. On top of that, only about one questionnaire in four came back, and Roosevelt’s opponents were keener to reply. Two million replies collected with a bias turned out worse than a few thousand collected fairly.

A cook tasting soup doesn’t eat the whole pot: he stirs it and tastes a spoonful. Stirring is the randomness of the sample, and it matters more than the size of the spoon. And the size of the pot hardly matters at all: well-stirred soup can be tasted with a spoon just as well from a saucepan as from a cauldron. The Literary Digest was scooping from one edge.

Let’s compute how far off an honest poll is. Let the share of park supporters in the city be $p$. The answer of respondent number $i$ is a random variable $X_i$: $1$ if they’re for it and $0$ if against. Its expectation is $p$ and its variance $p(1 - p)$ (chapter 48). The share in favour in the poll is the mean $\bar X = \frac{X_1 + \dots + X_n}{n}$. How widely does it scatter around $p$?

The idea: the square of a sum is the area of a square made of pieces, and for independent variables only the pieces on the diagonal survive on average.

The expectation of the mean is $\mu$ by linearity: $\mathbb E\bar X = \frac1n(\mathbb E X_1 + \dots + \mathbb E X_n) = \mu$. For the variance introduce the deviations $Y_i = X_i - \mu$, with $\mathbb E Y_i = 0$ and $\mathbb E Y_i^2 = \sigma^2$. Then $\bar X - \mu = \frac{Y_1 + \dots + Y_n}{n}$, and $\operatorname{Var}\bar X = \mathbb E(\bar X - \mu)^2 = \frac{1}{n^2}\,\mathbb E\,(Y_1 + \dots + Y_n)^2$. The square of the sum is the area of a square with side $Y_1 + \dots + Y_n$. Cut it into $n^2$ rectangles $Y_iY_j$, as in chapter 7 (in the picture $n = 3$): $(Y_1 + \dots + Y_n)^2 = \sum_{i,j} Y_iY_j$. The expectation of a sum is the sum of the expectations, so we can work piece by piece. On the diagonal are the squares $Y_i^2$. The expectation of each is $\sigma^2$: that’s the definition of variance. There are $n$ of them, $n\sigma^2$ together. Off the diagonal are the rectangles $Y_iY_j$ with $i \ne j$. The variables $Y_i$ and $Y_j$ are independent, and the expectation of a product of independent variables is the product of the expectations: for discrete variables $\sum_{a,b} ab\,P(Y_i = a)\,P(Y_j = b) = \sum_a a\,P(Y_i = a) \cdot \sum_b b\,P(Y_j = b)$. So $\mathbb E(Y_iY_j) = 0 \cdot 0 = 0$. In a single experiment such a rectangle can have a negative “area”, when $Y_i$ and $Y_j$ have different signs; on average the pluses and minuses cancel out. What’s left is $\mathbb E(Y_1 + \dots + Y_n)^2 = n\sigma^2$, so $\operatorname{Var}\bar X = \frac{n\sigma^2}{n^2} = \frac{\sigma^2}{n}$, and the standard deviation is $\sigma/\sqrt n$. The sum of $n$ independent deviations grows not like $n$ but like $\sqrt n$: they partly cancel each other.

For a poll $\sigma^2 = p(1 - p)$, and with $p$ around a half the standard error of a thousand replies is $\sqrt{0.25/1000} \approx 0.016$: one and a half percentage points. The size of the city doesn’t appear in the formula: that’s where the cook’s calm comes from. Strictly speaking, sampling without replacement brings in a factor $\sqrt{(N - n)/(N - 1)}$, where $N$ is the size of the population, but for a thousand people out of a million it differs from one by five ten-thousandths. The formula does have a square root, though: to halve the error you have to ask four times as many people. That’s why a typical national poll is one or two thousand people; beyond that, precision gets expensive too fast.

Now we can pay the debt from the first investigation and explain the division by $n - 1$.

The number $n - 1$ is called the number of degrees of freedom: the deviations from the mean add up to zero, so only $n - 1$ of them can vary freely, and the last is determined by the others. Note that $s$ itself is on average slightly smaller than $\sigma$ (the square root of a mean isn’t the mean of the square roots), but the correction is small and usually neglected.

The standard error tells you how far the sample mean usually strays from the truth. The central limit theorem adds the shape: for large $n$ the distribution of $\bar X$ is close to normal with mean $\mu$ and standard deviation $\sigma/\sqrt n$. And a normal variable deviates from its mean by more than $1.96$ standard deviations in only $5$ % of cases. Hence the recipe every poll lives by.

The trick is to look at the same event from the side of $\mu$.

The sample gave a mean $\bar x$. Mark off $1.96\,\mathrm{SE}$ from it in both directions: that gives an interval centred on $\bar x$. The interval covers $\mu$ exactly when the distance from $\bar x$ to $\mu$ is at most $1.96\,\mathrm{SE}$. Distance is symmetric, so this is the same as $\bar x$ falling in a band of the same width around $\mu$. Move $\bar x$: while the point is in the band, its interval reaches $\mu$; once it leaves the band, it doesn’t. The probability that $\bar X$ falls in the band is the area under the bell above the band. The number $1.96$ is chosen so that this area is $0.95$: $\Phi(1.96) - \Phi(-1.96) \approx 0.95$, where $\Phi$ is the standard normal distribution function. So the interval covers $\mu$ with probability $95$ %. What’s random here isn’t $\mu$ but $\bar X$, and with it the interval itself. For a $99$ % level one takes $2.58$ instead of $1.96$, for $90$ % one takes $1.64$.

The sample mean: the best guess about the unknown $\mu$. How many standard errors to mark off. For the $95$ % level it’s the number $z$ with $\Phi(z) = 0.975$: $2.5$ % is left for each tail. The standard error, with the unknown $\sigma$ replaced by the sample $s$. With at least a few dozen observations this hardly matters; for small samples $1.96$ is replaced by a slightly larger number from Student’s distribution. Example: at a dairy $100$ cartons labelled “1 litre” were weighed. The mean is $996$ ml, $s = 12$ ml. The standard error is $12/\sqrt{100} = 1.2$ ml, the interval $996 \pm 2.35$, that is, from $993.6$ to $998.4$ ml. A litre isn’t in it: the filling machine needs adjusting.

For a proportion the formula is the same, only $s$ is computed from the proportion itself: $s = \sqrt{\hat p\,(1 - \hat p)}$, where $\hat p$ is the share in favour in the poll:

$$\hat p \pm 1.96\,\sqrt{\frac{\hat p\,(1 - \hat p)}{n}}.$$

In our poll $\hat p = 0.52$, $n = 1000$, and $1.96\sqrt{0.52 \cdot 0.48/1000} \approx 0.031$. The interval runs from $48.9$ to $55.1$ %, hence the famous “plus or minus three per cent”. Hence the verdict too: a half, $50$ %, lies inside the interval, and the poll hasn’t proved that park supporters are a majority. The “52 against 48” gap is worse still. If the share for is a point higher, the share against is a point lower, so the difference between the shares wobbles twice as much as each of them, and its margin of error is about $6$ points, more than the gap itself.

One last thing. Plus or minus three per cent is only the random error: who happened to be in the sample and who didn’t. Systematic errors like the ones that sank the Literary Digest (who didn’t pick up the phone, who fibbed, how the question was worded) aren’t included in those three per cent and don’t shrink as the sample grows. Statisticians call such an error sampling bias.

How many people must be surveyed so that, with a proportion of about $50$ %, the poll’s margin of error (half the width of the 95 per cent interval) is at most one percentage point?

We need $1.96\sqrt{0.25/n} \le 0.01$, that is, $\sqrt n \ge 1.96 \cdot 0.5/0.01 = 98$, and $n \ge 98^2 = 9604$. Three times more precise than a “standard” poll takes almost ten times as many people.

A thousand people is enough to learn a proportion to within about three points, and not enough to tell $52$ from $50$. And if the sample is biased, two million won’t help either.

The lady tasting tea

Scientists prove new drug lowers blood pressure. The result is statistically significant, p < 0.05

The word “prove” in such headlines almost always means one thing: the result passed a statistical hypothesis test. The logic of such a test is best shown by a story that Ronald Fisher made famous in his book The Design of Experiments (1935).

In the 1920s Fisher worked at the Rothamsted agricultural research station in England. As the story goes, one day he poured a cup of tea for a colleague, the algae specialist Muriel Bristol, and she remarked that she preferred the milk to be poured into the cup first, and could taste the difference. How can such a claim be tested?

Fisher described this experiment. Eight cups are prepared: four with the milk poured first, four with the tea first. The cups are served in random order. The lady knows there are four of each and has to pick out the four “milk-first” ones. If she can’t tell the difference at all, her choice is just one of $\binom84 = 70$ equally likely sets of four (chapter 45), and she picks all four cups correctly with probability $1/70 \approx 0.014$. If that happens anyway, there are two versions: either the lady really can tell, or an event with a chance of one and a half per cent has occurred. It’s more reasonable to believe the lady. According to one witness’s recollection, she didn’t make a single mistake.

And if she got three cups out of four right? That seems pretty good too. But with a random choice exactly three are right in $\binom43\binom41 = 16$ ways out of $70$, and the probability of getting at least three right is $\frac{16 + 1}{70} \approx 0.24$. One person in four who tastes nothing would do as well. That’s not proof.

For the lady who got all the cups right, $p = 1/70 \approx 0.014$; for the one who got three, $p \approx 0.24$. The threshold $\alpha = 0.05$ caught on largely thanks to Fisher: back in 1925 he called it a convenient boundary, not a law of nature. Hence his remark that six cups aren’t enough: with three “milk-first” cups and three ordinary ones, guessing everything by chance has probability $1/\binom63 = 1/20$, exactly $0.05$.

A drug trial works the same way. The null hypothesis is that the drug doesn’t work, and the difference between the groups given the drug and the placebo is explained by who ended up in which group. The p-value says how rare the observed difference would be if that were so.

A drug trial gave $p = 0.03$. Which statement is correct?

If the drug doesn’t work, a difference between the groups as large or larger would arise by chance in about 3 % of trials The drug doesn’t work with probability 3 % A repeat trial will show the effect again with probability 97 % The drug’s effect is large

A p-value is the probability of the data given the hypothesis, not of the hypothesis given the data. It’s the same trap as the prosecutor’s fallacy from the chapter on probability.

A decision based on a threshold can go wrong in two ways.

The threshold $\alpha = 0.05$ keeps only false alarms under control, and only in a single test. That’s where the main trap lies.

Twenty jelly beans

In 2011 the webcomic xkcd published an episode called “Significant”. Scientists are asked to check whether jelly beans cause acne. They find no link. Then they test twenty colours of jelly bean separately, and for the green ones $p < 0.05$. Next morning in the paper: “Green jelly beans linked to acne! Only 5 % chance of coincidence”.

The joke is accurate down to the arithmetic. If jelly beans have nothing to do with it, each of the twenty tests gives a false alarm with probability $0.05$. But the probability that at least one of the twenty gives one is much greater.

The trick is to see probability as area. Let’s draw the case of two tests.

Let the outcome of each test be determined by a random number from the interval $[0, 1]$: an alarm if it’s less than $\alpha$. For a continuous test the p-value itself serves as this number. Two independent such numbers are a random point of the unit square, and the probability of any event is the area of the corresponding part of the square. An alarm in the first test is a vertical strip of width $\alpha$, in the second a horizontal strip of height $\alpha$. Each takes up a share $\alpha$ of the square. No alarms at all is the remaining square with side $1 - \alpha$ and area $(1 - \alpha)^2$. At least one alarm is the corner made of the two strips, with area $1 - (1 - \alpha)^2$. At $\alpha = 0.05$ that’s $0.0975$, almost twice $\alpha$. For $m$ independent tests it’s the same in an $m$-dimensional cube: the probabilities of “no alarm” multiply, and there are no alarms at all with probability $(1 - \alpha)^m$. At $m = 20$ at least one alarm happens with probability $1 - 0.95^{20} \approx 0.64$. If the tests are dependent, the alarm regions can lie any way at all, but the area of their union is at most the sum of their areas: shared pieces are counted twice in the sum. So the probability of at least one alarm is at most $m\alpha$. Hence the Bonferroni correction: to keep the overall probability of a false alarm in $m$ tests at most $\alpha$, each test is run with threshold $\alpha/m$.

In real life twenty tests rarely line up so honestly. A researcher tries different ways of analysing the data: exclude outliers or keep them, adjust for age or not, measure after a week or after a month, stop recruiting participants as soon as p dips below $0.05$. Each choice on its own looks innocent, but together they’re twenty jelly beans, only invisible. This is called p-hacking.

The consequences turned out to be serious. In 2005 the epidemiologist John Ioannidis published a paper with a provocative title, “Why Most Published Research Findings Are False”. In 2015 a large group of psychologists, the Open Science Collaboration, repeated a hundred experiments from papers in leading journals. In the original papers $97$ % of the results were significant; in the repetitions, $36$ %. After this replication crisis many journals began accepting papers with pre-registration: before collecting data the researcher announces which hypothesis will be tested and how.

A researcher tested $10$ independent hypotheses at level $\alpha = 0.05$, and all of them are in fact false (there are no effects). With what probability will they announce at least one “discovery”? Round the answer to three decimal places.

No false alarms at all: probability $0.95^{10} \approx 0.599$. At least one: $1 - 0.599 = 0.401$, almost forty per cent. With the Bonferroni correction each hypothesis would be tested at level $0.005$, and the probability of a false discovery would be $1 - 0.995^{10} \approx 0.049$.

“p < 0.05” means only that such a result rarely arises by chance if there’s no effect. Ask how many hypotheses were tested before this one was found, and whether anyone has replicated the result.

Chocolate and Nobel prizes

More chocolate, more Nobel laureates: correlation 0.79

This headline is almost real. In 2012 the cardiologist Franz Messerli published a note in the New England Journal of Medicine: countries that eat more chocolate have more Nobel laureates per head of population. Over $23$ countries the correlation coefficient came out at $0.79$. The note was written with a good deal of irony, but the news reported it in earnest.

What does the number $0.79$ mean? The correlation coefficient measures how closely the points of a scatter plot hug a straight line.

The sum of products of deviations. A point to the right of and above the centre $(\bar x, \bar y)$, or to the left and below, gives a positive term; a point to the right and below, or to the left and above, a negative one. If the points stretch along a rising line, the pluses win. Normalisation: we divide by the product of the “lengths” of the deviations. After this $r$ doesn’t depend on the units (whether chocolate is counted in grams or kilograms) and, as we’ll now prove, lies between $-1$ and $1$. Example: for the points $(1, 1)$, $(2, 3)$, $(3, 2)$ the means are $\bar x = \bar y = 2$, the deviations in $x$ are $-1, 0, 1$ and in $y$ are $-1, 1, 0$. The numerator is $1 + 0 + 0 = 1$, the denominator $\sqrt2 \cdot \sqrt2 = 2$, and $r = 0.5$.

The same trick as with the sum of squares: the whole column of deviations is one arrow in a space of $n$ dimensions.

Collect the deviations into two vectors: $\mathbf u = (x_1 - \bar x,\ \dots,\ x_n - \bar x)$ and $\mathbf v = (y_1 - \bar y,\ \dots,\ y_n - \bar y)$. The numerator of $r$ is their dot product, the denominator the product of their lengths: $r = \frac{\langle\mathbf u, \mathbf v\rangle}{|\mathbf u|\,|\mathbf v|}$. Two arrows starting from one point always lie in one plane, and that’s the plane we draw. First, the line. For any slope $b$, choosing the shift of the line $y = a + bx$ means representing the numbers $y_i - bx_i$ by a single number $a$, and by the theorem on sums of squares the best $a$ is their mean, $\bar y - b\bar x$. So the best line passes through $(\bar x, \bar y)$, and its vertical deviations are $(y_i - \bar y) - b(x_i - \bar x)$: the coordinates of the vector $\mathbf v - b\mathbf u$. Expand the squared length: $|\mathbf v - b\mathbf u|^2 = |\mathbf v|^2 - 2b\,\langle\mathbf u, \mathbf v\rangle + b^2|\mathbf u|^2$. This is a quadratic in $b$; it takes its smallest value at the vertex of the parabola (chapter 10), at $b = \frac{\langle\mathbf u, \mathbf v\rangle}{|\mathbf u|^2}$. In the picture $b\mathbf u$ is the projection of $\mathbf v$ onto the line of the vector $\mathbf u$, and the remainder $\mathbf v - b\mathbf u$ is perpendicular to it. Substitute this $b$: the smallest value is $|\mathbf v|^2 - \frac{\langle\mathbf u, \mathbf v\rangle^2}{|\mathbf u|^2} = |\mathbf v|^2\,(1 - r^2)$. A sum of squares is never negative, so $1 - r^2 \ge 0$, that is, $|r| \le 1$. Zero occurs only when the remainder is zero, $\mathbf v = b\mathbf u$: all the points lie on the line. By Pythagoras’ theorem the length of the projection is $|\mathbf v|\,|r|$, so $r$ is the cosine of the angle between $\mathbf u$ and $\mathbf v$. The slope of the best line is $b = \frac{\langle\mathbf u, \mathbf v\rangle}{|\mathbf u|^2} = r\,\frac{|\mathbf v|}{|\mathbf u|} = r\,\frac{s_y}{s_x}$: the factor $\sqrt{n - 1}$ in $s_x$ and $s_y$ cancels.

Praise hurts, scolding helps?

The last line of the theorem implies a phenomenon that confuses newspapers and scientists alike. Convert $x$ and $y$ into standard units: how many sample standard deviations a value is from the mean. Then $\bar x = \bar y = 0$, $s_x = s_y = 1$, and the best line is simply $y = r\,x$. If the correlation isn’t perfect, $|r| < 1$, and the forecast of $y$ is always closer to the mean than $x$ itself.

Galton was the first to notice this. In 1886 he compared the heights of adult children with those of their parents and saw that the children of very tall parents are also tall on average, but shorter than their parents, while the children of very short parents are taller on average. Galton called this “regression towards mediocrity”; hence the word “regression”, which is now used for any fitting of a line to points.

Daniel Kahneman recalled telling Israeli Air Force flight instructors in the 1960s that reward teaches better than punishment. One instructor objected: when he praised a cadet for an excellent landing, the next one usually went worse, and when he scolded one for a bad landing, the next was better. He saw this as the effect of his words. In fact he was seeing regression to the mean: a landing is made up of skill and luck, and after an exceptionally good attempt the next one will most likely be ordinary, praise or no praise.

Regression to the mean explains much of what looks like an effect. A patient goes to the doctor when they feel at their worst, and will get better under almost any treatment. Cameras are put up at the junctions that had the most accidents last year, and the next year there are fewer accidents there, partly for no reason at all. That’s why a serious trial has a control group: it was recruited the same way, and regression to the mean is the same in it as in the main group.

The correlation between the marks on two tests is $r = 0.6$. A pupil scored two standard deviations above the mean on the first. By how many standard deviations above the mean is the best forecast of their second mark?

In standard units the best forecast is $r\,x = 0.6 \cdot 2 = 1.2$. The pupil will most likely do well again, but not quite so brilliantly: part of their first success was luck, and luck doesn’t repeat to order.

Whose fault is it

Back to the chocolate. The correlation of $0.79$ is real, but it doesn’t say that chocolate makes people cleverer. There are at least three explanations. The first is chance: among a multitude of pairs of indicators some are bound to turn out linked, like the green jelly beans. The second is reverse causation: maybe Nobel laureates celebrate their prizes with chocolate. The third, and the most plausible, is a common cause: rich countries both eat more chocolate and spend more on science. A hidden common cause is called a confounder.

One more warning: $r$ captures only straight-line relationships. For points on the parabola $y = x^2$, symmetric about the axis, $r = 0$, even though $y$ is completely determined by $x$. And the four data sets that Francis Anscombe devised in 1973 have the same means, variances, $r \approx 0.82$ and the same best line $y = 3 + 0.5x$, yet on a chart they look nothing like one another. Before computing a correlation, plot the points.

So how do you find out whether chocolate works? Only by an experiment in which chocolate is handed out by lot: participants are randomly divided into two groups, one gets chocolate and the other doesn’t. The lot knows nothing about wealth, age or habits, so all confounders are shared equally between the groups on average. This idea, randomisation, was also brought into science by Fisher. The tea experiment starts with it too: the order of the cups there was chosen by lot.

Correlation measures how closely the points hug a line and says nothing about what causes what. Ask what else could have driven both quantities, and whether there was an experiment with random allocation.

Berkeley, 1973

University weeds out women: 44 % of men admitted to graduate school, only 35 % of women

This case is real. In autumn 1973, $8442$ men and $4321$ women applied to graduate school at the University of California, Berkeley. $44$ % of the men and $35$ % of the women were admitted. With numbers like these a nine-point difference can’t be explained by chance: the margin of error of the difference is under two points. The university grew worried, and the statistician Peter Bickel and colleagues set out to investigate.

Admission decisions at Berkeley were made not for the university as a whole but by departments. Bickel, Hammel and O’Connell looked at each department separately. Here are the six largest; in their 1975 paper the departments are labelled with letters.

In four departments out of six women were admitted more often than men, in the other two slightly less often. Yet in total across the same six departments $45$ % of men and $30$ % of women were admitted. The answer lies in where people applied. More than half the men, $1385$ out of $2691$, went for departments A and B, which admitted almost two thirds. Women rarely applied there, $133$ applications out of $1835$, and mostly went for C to F, which are hard for everyone to get into: they admitted between six and thirty-seven per cent.

There’s no magic here, only the arithmetic of averages, the same as with the plank and weights.

The trick is that the overall share is again a centre of mass, as on the plank with weights. The picture is for Berkeley, with departments A and B merged into “easy” and C to F into “hard”.

Part $j$ had $n_jp_j$ successes, $n_1p_1 + \dots + n_kp_k$ in all. Divide by $N$: $p = \frac{n_1}{N}p_1 + \dots + \frac{n_k}{N}p_k$. The weights $w_j = n_j/N$ are non-negative and add up to $1$, so $p$ is no less than the smallest $p_j$ (replace all the $p_j$ by the smallest and the sum can only go down) and no greater than the largest. For two parts $p = w\,p_1 + (1 - w)\,p_2 = p_2 + w\,(p_1 - p_2)$: a point between $p_2$ and $p_1$, a fraction $w$ of the way from $p_2$. Put $w$, the share of applications to the easy departments, on the horizontal axis and the overall share admitted on the vertical. For men the overall share lies on the segment from $(0,\ 25.5\,\%)$, all applying to the hard ones, to $(1,\ 62.5\,\%)$, all applying to the easy ones. For women there’s a similar segment: from $26.5$ % to $79.7$ %. Women were admitted more often both in the easy and in the hard departments, so their segment lies entirely above the men’s. But the groups stand at different places on their segments. For men $w = 1385/2691 \approx 0.51$, and their point is at $44.5$ %. For women $w = 133/1835 \approx 0.07$, almost at the left end, and their point is at $30.4$ %. A point on the upper segment far to the left has ended up below a point on the lower segment further to the right. If the weights were the same, the points would stand on one vertical, and the upper segment would stay on top: $p_{\text{w}} - p_{\text{m}} = \sum_j w_j\,(p_{\text{w},j} - p_{\text{m},j}) > 0$. So Simpson’s paradox is possible only when the groups are distributed differently among the parts.

If women had spread their applications over the six departments the same way as men, then with the same chances in each department about $52$ % of them would have been admitted, more than the men. Bickel and his co-authors concluded that if the departments favoured anyone, it was rather women. The question shifted from the admissions committees to something else: why did women more often apply where places were scarce? But that one isn’t answered by arithmetic.

The paradox turns up in medicine too. In 1986 British urologists compared two operations for removing kidney stones: open surgery and the minimally invasive keyhole procedure. Overall the keyhole one succeeded more often: $83$ % against $78$ %. But for small stones open surgery worked in $93$ % of cases against $87$ %, and for large ones in $73$ % against $69$ %. Open surgery was better in both groups and worse in total: it was used more often in the difficult cases, with large stones.

Which comparison should you believe: by group or in total? The numbers themselves won’t tell you; you need to understand where they came from. The size of the stone affected both the choice of operation and its outcome, that is, it was a confounder, and the comparison has to be made within the groups. There are opposite cases too: if the groups were themselves formed under the influence of the treatment, splitting by them only confuses things, and the total is what to believe. Statistics tells you what to compute; which of the numbers answers your question is decided by understanding how the situation works.

Department by department, women were admitted no less often than men. Their overall percentage is lower because they more often applied where everyone was turned down. Before comparing totals, ask whether the groups are distributed the same way over the parts.

You can practise the calculations from all five investigations (median and spread, standard error, a poll’s margin of error, multiple testing, a forecast with regression to the mean, and overall shares) in the trainer.

Where next

In all five stories the observations were independent: the respondents didn’t confer, the cups knew nothing of one another, each jelly bean had its own group. Independence held up both the standard error and the formula $1 - 0.95^{20}$. But much in the world works differently. The letters of a text depend on the previous ones: in English, q is almost always followed by u. Tomorrow’s weather resembles today’s. A user clicks a link that’s on the page they have open. In statistics we always assumed observations were independent, and here each one depends on what came before. How do you describe such dependence, and can you measure how much information a message carries? In 1913 Andrey Markov counted the vowels and consonants in Eugene Onegin, and in 1948 Claude Shannon answered the second question with a single formula. That’s chapter 50.