Chapter 1
Prediction, inference and data

Lectures 1–2 · why we collect data, how we collect it, and why that decides what it can tell us

Count the predictions that touched your life today. Your phone guessed the next word you would type. A streaming service guessed what you'd want to watch. A fast-food app offered you a deal it predicted would bring you back this week. Somewhere, an airline sold more tickets for your flight than the plane has seats, predicting that some passengers wouldn't show up.

None of these systems can see the future. Each one found a pattern in past data and bet that the pattern would continue. And each one was designed knowing it would sometimes be wrong: the deal you ignore, the word you delete, the passenger who has to be bumped. The goal is never to be right every time. It is to keep the cost of being wrong at a level everyone can live with.

That, in one sentence, is what statistics is for. This chapter sets up the whole course: what prediction and inference are, where data come from, and why the way data are collected decides what they can and cannot tell you.

There is very little mathematics in this chapter, and that is deliberate. The ideas here are the ones that decide whether a calculation means anything at all. A perfectly computed average of badly collected data is still a bad answer, and no formula later in the course can rescue it. So read this chapter the way a detective reads a crime scene: before asking what the numbers say, ask where the numbers came from.

Here is the route through the chapter.

• Sections 1.1 and 1.2 separate the two big goals of statistics: prediction (guessing what will happen) and inference (learning the mechanism that makes it happen).
• Sections 1.3 and 1.4 are about the raw material: who or what we measure (populations and samples), and what we measure about them (variables).
• Section 1.5 asks the hardest question in the chapter: when does an association in data reveal a cause?
• Section 1.6 introduces the assumption behind almost every method in the course, independence.
• Section 1.7 puts it all on one map.

Lecture 1 — Prediction, inference and sampling

1.1 Predicting under uncertainty

The seat that isn't empty

Airlines know that some ticket holders don't show up. Fly with empty seats and you lose money you could have earned. So they sell a few more tickets than there are seats. Sell too many, and more people show up than can fit: someone is bumped, and the airline pays compensation (and gets bad press). The airline must decide today how many tickets to sell, without knowing who will turn up.

Two things make this decision possible:

1. A pattern. From years of past flights, the airline knows roughly what fraction of ticket holders show up.
2. A cost for each kind of mistake. An empty seat costs one lost fare. A bumped passenger costs compensation, a hotel, perhaps a lost customer.

With those two ingredients, the airline can choose the number of tickets that keeps the average cost of its mistakes as small as possible. Try it.

Notice three things.

• Even the best choice still makes mistakes on many flights. Uncertainty doesn't go away; it gets managed.
• When bumping becomes more expensive, the best number of tickets goes down. The right prediction depends on the cost of being wrong, not just on the data.
• If passengers become less predictable, the airline needs a bigger safety margin. Uncertainty itself has a price.

The same structure is everywhere. A pharmacy predicts how many doses of a vaccine to stock. A hospital predicts how many nurses it needs on a Friday night. A forensic lab predicts how long its backlog of DNA samples will take to clear. In every case there is a pattern learned from data, an unknown future, and a cost for each way of being wrong.

The definition has three ingredients, and it helps to name them in the airline story.

• The unknown $y$ is the number of passengers who will show up.
• The information $x$ is everything known when tickets go on sale: the route, the day of the week, past show-up rates.
• The loss measures how bad each outcome is. Here the airline doesn't just guess $y$; it makes a decision (how many tickets to sell) whose cost depends on $y$. With $c_E$ the cost of an empty seat and $c_B$ the cost of bumping one passenger, the loss on a 100-seat flight is $$L(y)=\{\begin{array}{cc}{c}_E(100-y),\hfill & y<100,\hfill \\ 0,\hfill & y=100,\hfill \\ c_B(y-100),\hfill & y>100.\hfill \end{array}$$ The top line counts the cost of the $100-y$ empty seats, the bottom line the cost of the $y-100$ bumped passengers, and the middle line is a perfect fit.

For instance, with $c_E=\$250$ and $c_B=\$1000$, a flight on which 97 people show up costs $3$$\times \$250$$=\$750$, and a flight on which 102 show up costs $2$$\times \$1000$$=\$2000$. The loss is lopsided: being two over is much worse than being three under. That lopsidedness is exactly why the best number of tickets is not simply "whatever makes the average turnout equal 100". We will make "small on average" precise in Chapter 4, where it becomes the expected value of the loss.

set.seed(301)        # makes the simulation repeatable
# number who show up on each of 10,000 simulated flights
shows <- rbinom(10000, size = 110, prob = 0.9)
mean(shows)                  # average number who show up
mean(shows > 100)            # fraction with someone bumped
mean(pmax(100 - shows, 0))   # average empty seats
mean(pmax(shows - 100, 0))   # average bumped passengers

[1] 99.0012
[1] 0.3273
[1] 1.7871
[1] 0.7883

flights <- data.frame(shows = shows, bumped = shows > 100)
ggplot(flights, aes(x = shows, fill = bumped)) +
  geom_bar() +
  geom_vline(xintercept = 100.5, linetype = "dashed") +
  scale_fill_manual(values = c("grey60", "orange"),
                    labels = c("everyone seated", "someone bumped")) +
  labs(x = "passengers who show up", y = "number of flights", fill = NULL)

Prediction as classification

Sometimes the quantity to predict isn't a number but a category: spam or not spam, guilty or not guilty, headlamp glass or window glass. A prediction rule for a category is called a classifier. There are now two ways to be wrong, and they rarely cost the same.

library(MASS)          # comes with R; contains the fgl data
dim(fgl)               # 214 fragments, 10 variables
table(fgl$type)        # how many fragments of each type
# average barium, by type
aggregate(Ba ~ type, data = fgl, FUN = mean)

[1] 214  10

 WinF WinNF   Veh   Con  Tabl  Head 
   70    76    17    13     9    29 
   type          Ba
1  WinF 0.012714286
2 WinNF 0.050263158
3   Veh 0.008823529
4   Con 0.187692308
5  Tabl 0.000000000
6  Head 1.040000000

ggplot(fgl, aes(x = Ba, y = type)) +
  geom_jitter(height = 0.15, width = 0, alpha = 0.6) +
  geom_vline(xintercept = 0.2, linetype = "dashed") +
  labs(x = "barium oxide (% by weight)", y = NULL)

table(predicted_headlamp = fgl$Ba > 0.2,
      actually_headlamp  = fgl$type == "Head")

                  actually_headlamp
predicted_headlamp FALSE TRUE
             FALSE   180    3
             TRUE      5   26

Three lessons from the widget are worth stating plainly. First, every cut-off makes mistakes; moving it only changes which mistakes. Second, the data limit what any rule can do: three headlamp fragments contain no barium at all, so no cut-off catches them, and once the line is left of the lowest barium-bearing headlamp (0.4%), moving it further left only adds false alarms. Third, the best cut-off depends on the relative cost as well as on the data, although here it is stable: 0.3% is cheapest unless a miss costs less than a tenth of a false alarm, and then the cheapest rule jumps to 0.7%, accepting ten more misses to avoid one false alarm. This is the same trade-off a medical screening test or a fraud filter faces, and we will meet it again as the trade-off between Type I and Type II errors in Chapter 9.

1.2 Inference: learning the mechanism

The comet that came back

Prediction asks what will happen? There is a deeper question: why does it happen? What is the mechanism that produced the data, and how does it work?

The universe seems to run on rules. Planets keep to their orbits, drugs act on the body in consistent ways, coins land heads about half the time no matter who flips them. The data we collect are produced by these rules, blurred by randomness. If we can learn the rule from the data, we gain something prediction alone never gives us: we can explain the past, predict the future, and predict what would happen if we changed something.

Most mechanisms we study are less tidy than gravity. A drug lowers blood pressure, but by different amounts in different people. A spam filter's errors depend on what spammers try next. So we think of data as the output of a machine with two parts: a systematic part, the rule, and a random part, the noise.

To make the picture concrete, here is one way to write a data-generating process for the birth weights we meet in Example 1.4: $$\underset{\mathrm{\text{what we observe}}}{\underbrace{\mathrm{\text{birth weight}}}}=\underset{\mathrm{\text{rule}}}{\underbrace{{\mu }_{\mathrm{\text{group}}}}}+\underset{\mathrm{\text{noise}}}{\underbrace{\epsilon }}.$$ Here ${\mu }_{\mathrm{\text{group}}}$ is the average birth weight for babies in a given group (say, mothers who smoked), a fixed but unknown parameter, and $\epsilon$ ("epsilon") is everything else that makes one baby heavier or lighter than its group's average: genes, the length of the pregnancy, luck. We never see $\mu$ or $\epsilon$ separately, only their sum. Inference is the art of pulling them apart using many observations: the rule is what stays the same from baby to baby, and the noise is what doesn't.

The last four words of the definition matter. A conclusion without a measure of its uncertainty is not inference, it's a guess. Here is a mechanism you can't see inside. Learn about it the only way scientists ever can: by watching what it produces.

The widget shows the central fact of inference. With a handful of spins, a wide range of mechanisms could have produced what you saw. As data accumulate, the plausible range narrows, but it never shrinks to a single point. Chapter 8 shows how to compute that range.

Prediction or inference?

The same data can serve both goals, but they ask different questions.

library(MASS)
table(birthwt$smoke)     # 0 = did not smoke, 1 = smoked
aggregate(bwt ~ smoke, data = birthwt, FUN = mean)

  0   1 
115  74 
  smoke      bwt
1     0 3055.696
2     1 2771.919

bw <- birthwt
bw$smoker <- factor(bw$smoke, labels = c("did not smoke", "smoked"))
ggplot(bw, aes(x = smoker, y = bwt)) +
  geom_boxplot(outlier.shape = NA, width = 0.5) +
  geom_jitter(width = 0.12, alpha = 0.5) +
  labs(x = "mother during pregnancy", y = "birth weight (g)")

There is a third kind of question, the most ambitious of all.

1.3 Populations, samples and how we choose them

We can't measure everything

A police department seizes 2,000 small bags of white powder in a raid. To charge the suspects, the lab must show the bags contain an illegal drug. Testing every bag would take weeks. Instead, the lab tests a few dozen and draws a conclusion about all 2,000.

This is the basic move of statistics: look at a part to learn about the whole. Whether it works depends almost entirely on which part we look at.

A memory aid: parameter goes with population, statistic goes with sample.

mean(USArrests$Murder)     # all 50 states

[1] 7.788

set.seed(1973)
chosen <- sample(rownames(USArrests), 10)   # 10 random states
chosen
mean(USArrests[chosen, "Murder"])            # the sample mean

 [1] "New Hampshire"  "Delaware"       "Mississippi"   
 [4] "Maine"          "Kentucky"       "North Carolina"
 [7] "West Virginia"  "New York"       "Iowa"          
[10] "Rhode Island"  
[1] 7.13

means <- replicate(1000, mean(sample(USArrests$Murder, 10)))
mean(means)    # the average of the 1,000 sample means
sd(means)      # how much they typically vary (Chapter 2)

[1] 7.80708
[1] 1.223097

ggplot(data.frame(means), aes(x = means)) +
  geom_histogram(binwidth = 0.4, colour = "white") +
  geom_vline(xintercept = mean(USArrests$Murder), linetype = "dashed") +
  labs(x = "sample mean of 10 states (murder arrests per 100,000)",
       y = "number of samples")

Choosing the sample

The simplest fair method treats every unit identically.

Drawing names from a well-shaken hat produces a simple random sample. So does R's sample() function. One consequence is the fairness property we want:

A tiny case makes the theorem concrete. Take a population of $N=4$ units, $\{A,B,C,D\}$, and a simple random sample of $n=2$. The possible samples are $\{A,$$B\},$$\{A,$$C\},$$\{A,$$D\},$$\{B,$$C\},$$\{B,$$D\},$$\{C,$$D\},$ so $M=6$, each with probability $1/6$. Unit $A$ appears in $k=3$ of them, as does every other unit, so each unit is included with probability $3/6$$=1/2$$=n/N$. Check the double count: $nM$$=2$$\times 6$$=12$ and $Nk$$=4$$\times 3$$=12$.

The theorem says every unit has the same chance. That's necessary for a fair sample but not enough on its own: for a list of 20 units, a method that tosses a coin to pick either the first 10 or the last 10 also gives every unit a 50% chance, but it never picks a mix. A simple random sample makes every combination equally likely, which is a much stronger guarantee.

Other designs are often more practical:

• Stratified sample. Split the population into groups (strata) that differ in important ways, then take a simple random sample from each. A survey of NYC residents might sample separately from each of the five boroughs, so that every borough is represented.
• Cluster sample. Randomly choose whole groups, then measure everyone in them. A health survey might randomly choose 30 apartment buildings and interview every household in each. Cheaper, but people in the same building tend to be similar, so you learn less per person.
• Systematic sample. Take every $k$th unit from a list, starting at a random point. Easy, and usually fine, unless the list has a pattern that repeats every $k$ units.
• Convenience sample. Take whoever is easy to reach: your friends, the first 50 people through the door, an online poll. Cheap, fast, and often badly misleading.

See this for yourself. A researcher in New York who wants a quick estimate might only collect data from nearby states.

Watch two things as you draw more samples: where each method's estimates are centred, and how spread out they are. The convenience estimates cluster tightly, around 4.7, far from the truth of 7.788. In 1973 the Northeast had lower murder arrest rates than the South, so a sample drawn only from the Northeast is biased, however carefully it's drawn. Random and stratified samples land around the truth; the stratified samples vary a little less, because they never get unlucky and miss a whole region.

weighted.mean(USArrests$Murder, w = state.x77[, "Population"])

[1] 9.225131

Lecture 2 — Variables, studies and independence

1.4 Variables, data and measurement

Once we have a sample, we measure things about each unit. Data almost always arrive as a table: one row per unit, one column per measurement. Statisticians call this shape tidy data, and R stores it as a data frame (Section 0.5).

A quick test sorts most variables. Ask two questions in turn.

1. Does arithmetic make sense? If "the average value" means something, the variable is numerical. If the values are labels (even labels written as numbers, like a jersey number or a ZIP code), it is categorical.
2. If categorical: is there a natural order? Then ordinal, otherwise nominal. If numerical: can you list the possible values, with gaps between them? Then discrete, otherwise continuous.

The type of a variable decides which summaries and plots make sense (Chapter 2) and which probability models fit it (Chapters 4 and 5). An average ZIP code, or a histogram of blood types, is a sign that someone skipped this step. Getting it right is the first step of every analysis.

str(birthwt)

'data.frame':    189 obs. of  10 variables:
 $ low  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
 $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
 $ race : int  2 3 1 1 1 3 1 3 1 1 ...
 $ smoke: int  0 0 1 1 1 0 0 0 1 1 ...
 $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ ht   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ ui   : int  1 0 0 1 1 0 0 0 0 0 ...
 $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
 $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...

mean(birthwt$race)    # R computes this, but it is nonsense
birthwt$race <- factor(birthwt$race,
                       labels = c("white", "black", "other"))
table(birthwt$race)

[1] 1.846561

white black other 
   96    26    67 

The distribution is the central object of the whole course. In Chapter 2 we describe the distribution of a sample with summaries and pictures; from Chapter 3 on we describe the distribution of a population with probability.

both <- data.frame(arrests_1973 = USArrests$Murder,
                   offences_1976 = state.x77[, "Murder"],
                   row.names = state.abb)
head(both)
# how closely they agree (correlation, Chapter 6)
cor(both$arrests_1973, both$offences_1976)

   arrests_1973 offences_1976
AL         13.2          15.1
AK         10.0          11.3
AZ          8.1           7.8
AR          8.8          10.1
CA          9.0          10.3
CO          7.9           6.8
[1] 0.9336909

ggplot(both, aes(x = arrests_1973, y = offences_1976)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_point() +
  labs(x = "murder arrests per 100,000 (1973)",
       y = "murders per 100,000 (1976)")

1.5 Where data come from: experiments, observation and causes

Association is not causation, and why

In Example 1.4 the babies of smokers were 284 g lighter on average. Does smoking cause lower birth weight? The difficulty is that smokers and non-smokers might differ in other ways. Suppose hardship (low income, stress, poor access to prenatal care) makes a mother more likely to smoke and also lowers birth weight by itself. Then smokers' babies would be lighter on average even if smoking itself did nothing.

The clearest way to think about this is to draw what influences what.

A confounder creates an association between $X$ and $Y$ that has nothing to do with $X$ causing $Y$. The diagram shows two routes connecting smoking and birth weight: the direct arrow (the effect we want), and the back route through hardship (confounding). A simple comparison of averages mixes the two together.

A made-up example shows how strong the back route can be. Suppose, in an imaginary town, smoking has no effect at all on birth weight. Hardship lowers the average birth weight by 300 g, and mothers facing hardship are much more likely to smoke:

(Each group has 100 mothers; "average weight" is the average birth weight of their babies.)

Within each row, smokers and non-smokers have the same average, because smoking does nothing. Now ignore the rows and compare all smokers with all non-smokers.

• Smokers: 60 at 2,900 g and 20 at 3,200 g, so their average is $\frac{60\times 2900+20\times 3200}{80}$$=\frac{\mathrm{238,000}}{80}$$=2975\mathrm{\text{g}}.$
• Non-smokers: 40 at 2,900 g and 80 at 3,200 g, so their average is $\frac{40\times 2900+80\times 3200}{120}$$=\frac{\mathrm{372,000}}{120}$$=3100\mathrm{\text{g}}.$

Smokers' babies appear 125 g lighter, even though smoking does nothing in this town. The whole difference travelled along the back route: smokers are mostly from the hardship group, and that is why their babies are lighter. Comparing within each row (like with like) gives the right answer, zero. The widget lets you build your own version of this and watch the apparent effect appear from nowhere.

Observational studies and experiments

Why does randomization matter so much? Because a coin can't be influenced by anything about the person.

salk <- data.frame(
  design = rep(c("randomized", "observed control"), each = 3),
  group  = rep(c("vaccinated", "not vaccinated", "no consent"), 2),
  rate   = c(28, 71, 46, 25, 54, 44))
salk$group <- factor(salk$group, levels = c("vaccinated", "not vaccinated", "no consent"))
ggplot(salk, aes(x = group, y = rate)) +
  geom_col(width = 0.6) +
  facet_wrap(~ design) +
  scale_x_discrete(labels = c("vaccine", "no\nvaccine", "no\nconsent")) +
  labs(x = NULL, y = "polio cases per 100,000")

1.6 Independence and IID data

Almost every method in this course assumes that each observation brings new information. If you measure the refractive index of the same glass fragment five times, you learn a lot about that fragment and the measuring instrument, but very little about glass in general. Five fragments from five different windows tell you far more.

Here is the idea in everyday language; Chapter 3 makes it precise with probability.

Independence is surprisingly hard to fake. People asked to write down "random" coin flips almost always switch between H and T too often, because a run of five heads feels non-random. Independent flips don't remember what came before, so they produce longer runs than our intuition expects. Try it.

This isn't just a party trick. Auditors and forensic accountants use the same idea to spot invented numbers in expense reports and fabricated data in scientific papers: made-up data are usually too regular.

Why are long runs so common? Each flip has a 1 in 2 chance of matching the one before it, so once a flip has landed, the chance that the next four all match it (making a run of five) is $$\frac{1}{2}\times \frac{1}{2}\times \frac{1}{2}\times \frac{1}{2}=\frac{1}{16}.$$ That is unlikely at any one place, but in 50 flips there are 46 places where a run of five could start. Simulation shows that about 4 out of every 5 sequences of 50 fair flips contain a run of five or more identical faces. Your intuition judges each run on its own ("five in a row, how unlikely!") instead of counting all the places one could have happened. We will be able to compute these chances properly in Chapter 3.

airquality$date <- as.Date(paste(1973, airquality$Month,
                                 airquality$Day, sep = "-"))
ggplot(airquality, aes(x = date, y = Ozone)) +
  geom_line() + geom_point(size = 1, na.rm = TRUE) +
  labs(x = NULL, y = "ozone (ppb)")

oz <- airquality$Ozone
today    <- oz[-length(oz)]   # days 1 to 152
tomorrow <- oz[-1]            # days 2 to 153
cor(today, tomorrow, use = "complete.obs")

[1] 0.5515109

ggplot(data.frame(today, tomorrow), aes(x = today, y = tomorrow)) +
  geom_point(alpha = 0.7, na.rm = TRUE) +
  labs(x = "ozone today (ppb)", y = "ozone tomorrow (ppb)")

1.7 The map of the course

We can now see the whole course in one picture. There is a population, or more generally a data-generating process, that we want to understand. We observe a sample of data from it.

• Describing the data (Chapter 2): summaries and pictures of what we actually observed.
• Probability (Chapters 3–7) reasons forwards: if we knew the mechanism, what kinds of data would it produce, and how much would they vary from sample to sample?
• Inference (Chapters 8–9) reasons backwards: having seen the data, what can we say about the mechanism, and how sure can we be?
• Prediction uses what we've learned to guess the next case.

Backwards reasoning is only possible because we first understand forwards reasoning. To judge whether 38 cocaine-positive bags out of 40 is consistent with "at least 90% of all 2,000 bags are cocaine", we need to know what samples from such a population typically look like. That is why the course spends five chapters on probability before it returns to data.

1.8 Exercises

Each exercise has a hint and an answer button. Try the exercise first; open the hint if you are stuck; and only look at the answer to check your own. The answers are brief: in a homework or exam you would be expected to explain your reasoning in full sentences.