22 Inference

Chapter 21 produced a number. This chapter asks the harder question: what does that number license you to say about the world?

The gap is real. A sample mean of 10.4 is a fact about your data; “the population mean is about 10.4, give or take 0.6” is a claim about everyone you did not measure. Inference is the machinery for making that leap honestly, and most of the chapter is about the ways it gets overstated.

22.1 From sample to population

Two quantities, easy to conflate:

Symbol	Name	Status
\(\mu\), \(\sigma\), \(p\)	parameters	fixed, unknown
\(\bar{x}\), \(s\), \(\hat{p}\)	statistics	computed, random

The population parameter does not move; your estimate of it does. Every inferential statement is really a statement about how much your estimate would move if you did the study again.

Which is why the sampling scheme matters more than any subsequent calculation. A biased sample gives a precise answer to the wrong question, and no amount of statistics repairs it — the bootstrap warning in Section 21.9 applies to inference generally.

22.2 Sampling distributions

The sampling distribution of a statistic is its distribution across repeated samples from the population. You almost never see it — you have one sample — but everything in this chapter is a claim about its shape.

set.seed(3)
draws <- replicate(20000, mean(rnorm(25, mean = 5, sd = 2)))
c(
  centre_of_sampling_dist = mean(draws),
  sd_of_sampling_dist = sd(draws),
  theory = 2 / sqrt(25)
)

centre_of_sampling_dist     sd_of_sampling_dist                  theory 
              4.9999493               0.4010631               0.4000000

That standard deviation has its own name — the standard error — because it measures error in an estimate, not spread in the data:

\[ \text{SE}(\bar{X}) = \frac{\sigma}{\sqrt{n}} \tag{22.1}\]

The \(\sqrt{n}\) is the central economic fact of statistics: to halve your uncertainty you need four times the data.

22.3 The law of large numbers

As the sample grows, the sample mean converges to the population mean:

\[ \bar{X}_n \longrightarrow \mu \qquad \text{as } n \to \infty \]

This is the guarantee that collecting more data helps. It is also weaker than people assume: it says the average settles down eventually, not that it does so quickly, and it requires the mean to exist at all (Section 19.7).

set.seed(5)
x <- rnorm(10000, mean = 3, sd = 1)
running <- cumsum(x) / seq_along(x)
round(running[c(10, 100, 1000, 10000)], 4)

[1] 2.9211 3.0316 3.0174 3.0018

idx <- seq(10, 10000, by = 10)
draw_line(
  x = idx,
  y = running[idx],
  points = FALSE,
  xlab = "sample size",
  ylab = "running mean"
)

Figure 22.1: The running mean of a growing sample. Early values swing widely; later ones settle toward 3. The convergence is real but slow — the envelope narrows as \(1/\sqrt{n}\), not as \(1/n\).

Watch out

The gambler’s fallacy is a misreading of this theorem. The average converges because early deviations are diluted by later data, not because they are corrected. After a run of heads, the coin does not owe you tails — the imbalance simply becomes a smaller fraction of a larger total.

22.4 The central limit theorem

The law of large numbers says where \(\bar{X}\) goes. The central limit theorem says what its distribution looks like on the way:

\[ \bar{X}_n \;\approx\; \mathcal{N}\!\left(\mu, \frac{\sigma^2}{n}\right) \qquad\text{for large } n \tag{22.2}\]

Whatever the population looks like. Skewed, discrete, bounded, bimodal — the distribution of the mean tends to normal regardless. That is why the normal appears everywhere, and why so much inference works without knowing the population’s shape.

set.seed(7)
clt_sd <- function(n) {
  m <- replicate(8000, mean(rexp(n, rate = 1)))
  c(
    sd = sd(m), theory = 1 / sqrt(n),
    skewness = mean(scale(m)^3)
  )
}
sapply(c(1, 5, 30, 100), clt_sd)

             [,1]      [,2]      [,3]      [,4]
sd       1.013712 0.4459902 0.1832286 0.1013721
theory   1.000000 0.4472136 0.1825742 0.1000000
skewness 2.029239 0.8055666 0.4315786 0.2202400

The parent is exponential — strongly right-skewed, with skewness 2. The skewness of the sample mean falls steadily toward zero as \(n\) grows, which is normality arriving.

set.seed(8)
draw_density(
  list(
    `n = 1` = rexp(8000, 1),
    `n = 5` = replicate(8000, mean(rexp(5, 1))),
    `n = 30` = replicate(8000, mean(rexp(30, 1)))
  ),
  xlab = "sample mean"
)

Figure 22.2: Sampling distributions of the mean of \(n\) exponential observations. At \(n=1\) it is the parent’s sharp right skew; by \(n=30\) it is visibly symmetric. Normality is a property of the average, not of the data.

Watch out

“Large \(n\)” is not a fixed number. The more skewed the parent, the slower the convergence — and for a distribution with infinite variance the theorem does not apply at all, however much data you have.

The common advice of \(n > 30\) is a rule of thumb for mild skew, not a law.

22.5 Confidence intervals

A point estimate without a range is an overclaim. A confidence interval attaches one:

\[ \bar{x} \pm t^*\,\frac{s}{\sqrt{n}} \tag{22.3}\]

where \(t^*\) is the critical value from Student’s \(t\) with \(n-1\) degrees of freedom (Section 20.8).

set.seed(9)
samp <- rnorm(20, mean = 5, sd = 2)
se <- sd(samp) / sqrt(20)
tstar <- qt(0.975, df = 19)
c(
  estimate = mean(samp),
  se = se,
  lower = mean(samp) - tstar * se,
  upper = mean(samp) + tstar * se
)

 estimate        se     lower     upper 
4.6914550 0.4385854 3.7734853 5.6094248

The \(t\) rather than the normal matters at this sample size, and the cost of getting it wrong is measurable:

set.seed(10)
coverage <- function(crit_fn) {
  hits <- replicate(20000, {
    x <- rnorm(20, mean = 5, sd = 2)
    e <- sd(x) / sqrt(20)
    abs(mean(x) - 5) < crit_fn() * e
  })
  mean(hits)
}
c(
  with_t = coverage(function() qt(0.975, 19)),
  with_z = coverage(function() 1.96)
)

 with_t  with_z 
0.95095 0.93350

The \(t\) interval covers about 95% of the time, as advertised. Using \(1.96\) instead covers only about 93% — an interval that quietly claims more precision than it has.

Watch out

“95% confidence” describes the procedure, not this interval. It means: if you repeated the whole study many times, about 95% of the intervals produced would contain the true value.

For the one interval in front of you, the parameter is either in it or not. The statement you probably want — “there is a 95% chance \(\mu\) is in here” — is a credible interval and requires the Bayesian machinery of Section 21.8.

22.6 Hypothesis testing

The formal apparatus for “could this have been chance?”

State a null hypothesis \(H_0\) — usually “no effect” — and an alternative \(H_1\).
Choose a test statistic whose distribution under \(H_0\) you know.
Compute how extreme the observed value is under that distribution.
Reject \(H_0\) if it is extreme enough.

The logic is proof by contradiction, softened: assume no effect, and see whether the data would be surprising if that were true.

set.seed(12)
group_a <- rnorm(30, mean = 0, sd = 1)
group_b <- rnorm(30, mean = 0.8, sd = 1)
tt <- t.test(group_b, group_a)
c(
  difference = unname(diff(rev(tt$estimate))),
  t_statistic = unname(tt$statistic),
  p_value = tt$p.value
)

  difference  t_statistic      p_value 
0.9073787904 3.9814109700 0.0001932539

Watch out

Failing to reject \(H_0\) is not evidence that \(H_0\) is true. It means the data were consistent with it — which a tiny, noisy study guarantees regardless of the truth.

“No significant difference” and “no difference” are entirely different claims, and conflating them is how underpowered studies get read as evidence of absence.

22.7 p-values and what they are not

\[ p = P(\text{a statistic at least as extreme as observed} \mid H_0) \tag{22.4}\]

Note what is conditioned on what. The \(p\)-value is a probability about the data, given the null — and by Section 18.5, reversing a conditional requires a prior, which no \(p\)-value has.

A p-value is not	Why
\(P(H_0 \text{ is true})\)	that reverses the conditioning
\(P(\text{the result was chance})\)	same error, phrased casually
a measure of effect size	a tiny effect gives a small \(p\) at large \(n\)
a measure of importance	statistical and practical significance differ
replicable as such	\(p\) varies wildly across repeat studies

The effect-size point deserves a demonstration, because it is the one that misleads in practice:

set.seed(14)
tiny_effect <- function(n) {
  x <- rnorm(n, mean = 0.05, sd = 1)
  t.test(x)$p.value
}
sapply(c(50, 500, 5000, 50000), function(n) {
  mean(replicate(200, tiny_effect(n)) < 0.05)
})

[1] 0.040 0.185 0.970 1.000

The true effect is 0.05 standard deviations — negligible for any purpose. At \(n = 50\) it is almost never detected; at \(n = 50{,}000\) it is detected nearly always. A small \(p\)-value at large \(n\) says the effect is not exactly zero, not that it matters.

22.8 Type I and Type II errors

Two ways to be wrong, and they trade off:

	\(H_0\) true	\(H_0\) false
Reject \(H_0\)	Type I error (\(\alpha\))	correct
Do not reject	correct	Type II error (\(\beta\))

\(\alpha\) is the false positive rate you choose — conventionally 0.05. \(\beta\) is the false negative rate, which you do not choose directly; it follows from \(\alpha\), the sample size, and the true effect size.

Lowering \(\alpha\) makes false positives rarer and false negatives commoner. There is no setting that reduces both; only more data does that.

22.9 Power

Power is \(1 - \beta\): the probability of detecting an effect that is really there.

ns <- c(20, 50, 100, 200)
pw <- sapply(ns, function(n) {
  power.t.test(n = n, delta = 0.5, sd = 1)$power
})
rbind(n_per_group = ns, power = round(pw, 4))

               [,1]    [,2]     [,3]     [,4]
n_per_group 20.0000 50.0000 100.0000 200.0000
power        0.3377  0.6969   0.9404   0.9988

nn <- seq(10, 200, by = 5)
draw_line(
  x = nn,
  y = sapply(nn, function(n) {
    power.t.test(n = n, delta = 0.5, sd = 1)$power
  }),
  points = FALSE,
  xlab = "n per group",
  ylab = "power"
)

Figure 22.3: Power against sample size for a medium effect (\(d = 0.5\)) at \(\alpha = 0.05\). Reaching the conventional 80% target needs about 64 per group; 20 per group has only a one-in-three chance of finding a real effect.

power.t.test(power = 0.8, delta = 0.5, sd = 1)$n

[1] 63.76576

In machine learning

Underpowered studies are worse than useless, and the reason is subtle. In a low-power study, an effect must be large by chance to clear the significance threshold — so the effects that do get published are systematically overstated.

This is the winner’s curse, and it is a major contributor to the replication crisis: the literature fills with exaggerated estimates, each individually “significant”.

22.10 Multiple comparisons

Run one test at \(\alpha = 0.05\) and you have a 5% false positive rate. Run twenty independent tests on pure noise and:

\[ P(\text{at least one false positive}) = 1 - 0.95^{20} = 0.64 \]

m <- c(1, 5, 20, 100)
rbind(
  tests = m,
  P_any_false_positive = round(1 - 0.95^m, 4)
)

                     [,1]   [,2]    [,3]     [,4]
tests                1.00 5.0000 20.0000 100.0000
P_any_false_positive 0.05 0.2262  0.6415   0.9941

At 100 tests it is essentially certain. This is why exploratory analyses that try many things need correction:

Method	Controls	Idea
Bonferroni	family-wise error rate	test each at \(\alpha/m\)
Holm	family-wise error rate	Bonferroni, stepwise; strictly better
Benjamini–Hochberg	false discovery rate	allow a fixed share of discoveries to be false

set.seed(16)
pvals <- replicate(20, t.test(rnorm(30))$p.value) # all null
c(
  raw_below_0.05 = sum(pvals < 0.05),
  bonferroni_below_0.05 =
    sum(p.adjust(pvals, "bonferroni") < 0.05),
  bh_below_0.05 = sum(p.adjust(pvals, "BH") < 0.05)
)

       raw_below_0.05 bonferroni_below_0.05         bh_below_0.05 
                    1                     0                     0

Every one of these twenty tests is null by construction, so every rejection is a false positive. Correction removes them.

Watch out

The count that matters is every test you could have run, not every test you report. Trying several outcomes, subgroups or model specifications and reporting the one that reached significance is multiplicity whether or not the others appear in the write-up.

Pre-registration exists to make that count auditable.

22.11 Permutation tests

If the null hypothesis is “the group labels are irrelevant”, you can test it without any distributional assumption at all: shuffle the labels and see how often chance produces a difference as large as the one you observed.

set.seed(1)
a <- rnorm(30, mean = 0)
b <- rnorm(30, mean = 0.8)
observed <- mean(b) - mean(a)
pooled <- c(a, b)
set.seed(99)
perm <- replicate(5000, {
  shuffled <- sample(pooled)
  mean(shuffled[31:60]) - mean(shuffled[1:30])
})
c(
  observed = observed,
  permutation_p = mean(abs(perm) >= abs(observed)),
  t_test_p = t.test(b, a)$p.value
)

     observed permutation_p      t_test_p 
 0.8503164130  0.0002000000  0.0003331962

draw_histogram(
  perm,
  xlab = "difference under label shuffling"
)

Figure 22.4: The permutation null distribution of the group difference, built by shuffling labels 5,000 times. The observed difference sits far outside it, which is what a small p-value means — stated without assuming any distribution.

The permutation \(p\)-value closely matches the \(t\)-test’s, but it assumed nothing about normality — the null distribution was constructed from the data. Where a parametric assumption is doubtful, or the statistic is unusual, this is often the more trustworthy route.

22.12 From inference to prediction

Classical inference and machine learning ask different questions of the same data, and much confusion comes from not noticing which one is being asked.

	Inference	Prediction
Question	is there an effect, and how big?	how accurate is the forecast?
Target	a parameter	a future outcome
Validated by	sampling theory, \(p\)-values	held-out data
Prefers	interpretable, unbiased	accurate, bias accepted
Fears	false positives	overfitting

Both are legitimate. The mistake is importing habits across the divide: reading a \(p\)-value from a model selected by cross-validation, or judging a causal claim by its test-set accuracy.

The bias-variance tradeoff (Section 21.10) is where they part company. Inference has traditionally prized unbiasedness; prediction happily accepts bias for lower variance, which is what regularization does. Neither is wrong — they are optimizing different things.

In machine learning

The two are converging. Conformal prediction attaches finite-sample coverage guarantees to any model’s predictions; causal inference borrows machine learning to estimate nuisance functions while keeping valid inference on the effect of interest.

The dividing line is not between fields but between questions, and knowing which one you are asking is the whole skill.

22.13 Summary

Concept	Statement
Standard error	\(\sigma/\sqrt{n}\); four times the data halves it
Law of large numbers	\(\bar{X}_n \to \mu\)
Central limit theorem	\(\bar{X}_n \approx \mathcal{N}(\mu, \sigma^2/n)\), any parent
Confidence interval	\(\bar{x} \pm t^*s/\sqrt{n}\); a claim about the procedure
p-value	\(P(\text{data} \mid H_0)\), never \(P(H_0 \mid \text{data})\)
Type I / II	false positive \(\alpha\) / false negative \(\beta\)
Power	\(1-\beta\); 80% needs \(\approx 64\) per group at \(d = 0.5\)
Multiplicity	20 tests give a 64% chance of a false positive
Permutation test	build the null by shuffling; assumes nothing

22.14 Exercises

1. A sample of 50 has mean 12 and standard deviation 4. Construct a 95% confidence interval.

Solution

n1 <- 50
se1 <- 4 / sqrt(n1)
t1 <- qt(0.975, n1 - 1)
c(
  se = se1,
  margin = t1 * se1,
  lower = 12 - t1 * se1,
  upper = 12 + t1 * se1
)

        se     margin      lower      upper 
 0.5656854  1.1367874 10.8632126 13.1367874

Roughly \(12 \pm 1.14\). Note \(t^* = 2.01\) against the normal’s \(1.96\) — at \(n = 50\) the difference is small, but at \(n = 10\) it would be \(2.26\), a 15% wider interval.

2. Demonstrate the central limit theorem starting from a uniform population.

Solution

set.seed(18)
unif_means <- function(n) {
  m <- replicate(6000, mean(runif(n)))
  c(sd = sd(m), theory = sqrt(1 / 12) / sqrt(n))
}
sapply(c(1, 2, 5, 30), unif_means)

            [,1]      [,2]      [,3]       [,4]
sd     0.2904030 0.2019125 0.1308414 0.05268300
theory 0.2886751 0.2041241 0.1290994 0.05270463

The uniform is flat and bounded — as unlike a normal as a well-behaved distribution gets — yet the standard deviation of its sample mean follows \(\sigma/\sqrt{n}\) exactly, and the shape becomes normal remarkably fast. Even \(n = 5\) is close.

3. Compute the probability of at least one false positive across 10, 50 and 200 independent tests at \(\alpha = 0.05\).

Solution

mm <- c(10, 50, 200)
rbind(
  tests = mm,
  P_any = round(1 - 0.95^mm, 6),
  bonferroni_alpha = round(0.05 / mm, 6)
)

                      [,1]      [,2]       [,3]
tests            10.000000 50.000000 200.000000
P_any             0.401263  0.923055   0.999965
bonferroni_alpha  0.005000  0.001000   0.000250

At 200 tests a false positive is a certainty to six decimal places. Bonferroni’s threshold becomes brutal — \(0.00025\) — which is why it is often replaced by false discovery rate control when many tests are genuinely expected to be non-null.

4. How many observations per group are needed for 80% power at \(d = 0.3\)? At \(d = 0.8\)?

Solution

c(
  d_0.3 = power.t.test(power = 0.8, delta = 0.3, sd = 1)$n,
  d_0.5 = power.t.test(power = 0.8, delta = 0.5, sd = 1)$n,
  d_0.8 = power.t.test(power = 0.8, delta = 0.8, sd = 1)$n
)

    d_0.3     d_0.5     d_0.8 
175.38510  63.76576  25.52463

About 176, 64 and 26. Required \(n\) scales as \(1/d^2\), so halving the effect you want to detect quadruples the sample — the same \(\sqrt{n}\) arithmetic as Equation 22.1, seen from the design side.

5. Run a permutation test where the null is true and confirm the \(p\)-value behaves correctly.

Solution

set.seed(20)
a0 <- rnorm(25)
b0 <- rnorm(25) # same distribution: null is true
obs0 <- mean(b0) - mean(a0)
pool0 <- c(a0, b0)
perm0 <- replicate(5000, {
  s <- sample(pool0)
  mean(s[26:50]) - mean(s[1:25])
})
c(
  observed = obs0,
  permutation_p = mean(abs(perm0) >= abs(obs0))
)

     observed permutation_p 
  -0.09443933    0.75140000

A large \(p\)-value, as it should be: the observed difference is unremarkable among differences produced by shuffling. Under a true null the \(p\)-value is approximately uniform on \([0,1]\) — which is exactly why 5% of null tests cross 0.05, and hence why Section 22.10 matters.

6. Show that a \(p\)-value below 0.05 says little about effect size by comparing a small effect at large \(n\) with a large effect at small \(n\).

Solution

set.seed(22)
small_big_n <- rnorm(10000, mean = 0.04, sd = 1)
big_small_n <- rnorm(15, mean = 1.2, sd = 1)
rbind(
  small_effect_large_n = c(
    effect = mean(small_big_n),
    p = t.test(small_big_n)$p.value
  ),
  large_effect_small_n = c(
    effect = mean(big_small_n),
    p = t.test(big_small_n)$p.value
  )
)

                         effect           p
small_effect_large_n 0.02640107 0.008716254
large_effect_small_n 1.10068051 0.001162055

Both are significant, and the effects differ by a factor of about forty. Report the estimate and its interval, not just the \(p\)-value — the \(p\) tells you the effect is detectable, and the interval tells you whether it is worth anything.