3.2 Estimating the standard error
The standard error is the number we want. Where is it supposed to come from?
At the end of the previous section we identified the first obstacle to statistical inference. The standard error of the sample mean is
\[\mathrm{se}(\bar{X}) = \frac{\sigma}{\sqrt{n}}\]
but the population standard deviation \(\sigma\) is unknown.
This should not surprise us. Every quantity in that formula except \(n\) belongs to the population, while all we possess is a sample. If the population mean must be estimated, the population standard deviation must be estimated too.
The obvious strategy is the one that already worked. Replace the unknown population quantity by its sample analogue. For the mean this gave
\[\mu \quad\longrightarrow\quad \bar{X}\]
For the variance we might therefore try
\[\sigma^2 \quad\longrightarrow\quad \frac{1}{n}\sum_{i=1}^{n}(X_i - \bar{X})^2\]
It is an appealing estimator. It measures the average squared distance of the observations from their centre, exactly as the population variance does.
Unfortunately it is wrong. Not disastrously wrong, but systematically wrong: it underestimates the population variance.
Before trying to understand why, let us see the problem.
A statistical laboratory
As in the previous unit, we need a population where the truth is known. A fair die is ideal, because its variance is known exactly:
\[\mu = 3.5 \qquad \sigma^2 = \frac{35}{12} = 2.9167\]
Draw samples of ten rolls, compute the sum of squared deviations from the sample mean, and divide it by four different numbers. If any of these is a good estimator, its average across thousands of samples should land on 2.9167.
set.seed(11)
n <- 10
ss <- replicate(20000, {
d <- sample(1:6, n, replace = TRUE)
sum((d - mean(d))^2) # the sum of squared deviations
})
divisors <- c(n, n - 1, n - 2, n - 3)
tab <- data.frame(
divisor = divisors,
avg_variance = sapply(divisors, function(k) mean(ss / k)),
avg_sd = sapply(divisors, function(k) mean(sqrt(ss / k)))
)
round(tab, 4)#> divisor avg_variance avg_sd
#> 1 10 2.630 1.601
#> 2 9 2.922 1.688
#> 3 8 3.287 1.790
#> 4 7 3.757 1.914#> true_variance true_sd
#> 2.917 1.708Notice what the table shows. Dividing by \(n\) is always too small. Dividing by \(n - 2\) is too large, and \(n - 3\) is worse. Only one divisor consistently lands on the truth, and it is not \(n\).
It is \(n - 1\), which gives 2.9219 against a true 2.9167.
That observation raises two questions, and the rest of this section answers them in turn.
- Why does dividing by \(n\) underestimate the variance?
- Why is the correction exactly \(n - 1\), rather than \(n - 2\) or some other number?
The second question is the one that matters. Any downward adjustment would move the estimate in the right direction; only one gets the size right.
Why dividing by n fails
The reason is that \(\bar{x}\) is computed from the very observations whose spread we are measuring.
The sample mean sits, by construction, in the middle of the sample. It is the value that makes \(\sum(x_i - \bar{x})^2\) as small as it can possibly be — no other number, including the true \(\mu\), gives a smaller total. So when we measure the spread of the sample around its own mean, we measure it around the most flattering point available, and the answer comes out too small.
Once \(\bar{x}\) is known, the deviations are not free to take any values. They must sum to zero:
\[\sum_{i=1}^{n}(x_i - \bar{x}) = 0\]
Given the first \(n-1\) of them, the last is determined. Only \(n - 1\) of the \(n\) deviations carry independent information about the spread.
The quantity \(n - 1\) is called the degrees of freedom: the number of observations, less the number of quantities estimated from them along the way.
Why the correction is exactly n − 1
The argument above explains the direction of the error. It does not explain its size, and the size is what distinguishes \(n-1\) from \(n-2\). For that we need to take an expectation.
Write \(\widetilde{S}^2\) for the divide-by-\(n\) version, and simplify it first. Expanding the square,
\[ \begin{aligned} \widetilde{S}^2 &= \frac{1}{n}\sum_{i=1}^{n}(X_i - \bar{X})^2 \\[4pt] &= \frac{1}{n}\left(\sum_{i=1}^{n} X_i^2 - 2\bar{X}\sum_{i=1}^{n} X_i + n\bar{X}^2\right) \\[4pt] &= \frac{1}{n}\left(\sum_{i=1}^{n} X_i^2 - 2n\bar{X}^2 + n\bar{X}^2\right) && \text{since } \textstyle\sum X_i = n\bar{X} \\[4pt] &= \frac{1}{n}\sum_{i=1}^{n} X_i^2 \;-\; \bar{X}^2 \end{aligned} \]
Now take expectations:
\[E(\widetilde{S}^2) = \frac{1}{n}E\!\left(\sum_{i=1}^{n} X_i^2\right) - E(\bar{X}^2)\]
Both terms follow from rearranging the definition of variance, \(\mathrm{Var}(X) = E(X^2) - [E(X)]^2\), into \(E(X^2) = \mathrm{Var}(X) + [E(X)]^2\).
For a single observation, \(E(X_i) = \mu\) and \(\mathrm{Var}(X_i) = \sigma^2\), so
\[E(X_i^2) = \sigma^2 + \mu^2 \qquad\Longrightarrow\qquad E\!\left(\sum_{i=1}^{n} X_i^2\right) = n(\sigma^2 + \mu^2)\]
For the sample mean, the previous unit established \(E(\bar{X}) = \mu\) and \(\mathrm{Var}(\bar{X}) = \sigma^2/n\), so
\[E(\bar{X}^2) = \frac{\sigma^2}{n} + \mu^2\]
Substituting both,
\[ \begin{aligned} E(\widetilde{S}^2) &= \frac{1}{n}\,n(\sigma^2 + \mu^2) - \left(\frac{\sigma^2}{n} + \mu^2\right) \\[4pt] &= \sigma^2 + \mu^2 - \frac{\sigma^2}{n} - \mu^2 \\[4pt] &= \sigma^2\left(1 - \frac{1}{n}\right) \;=\; \sigma^2\,\frac{n-1}{n} \end{aligned} \]
So the divide-by-\(n\) estimator is too small by exactly the factor \((n-1)/n\) — not approximately, and not only for dice. Multiplying it by \(n/(n-1)\) removes the factor, and multiplying by \(n/(n-1)\) is the same as dividing by \(n-1\).
The sample variance and sample standard deviation are
\[s^2 = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2 \qquad\qquad s = \sqrt{s^2}\]
\(s^2\) is an unbiased estimator of \(\sigma^2\): \(E(s^2) = \sigma^2\). The adjustment from \(n\) to \(n-1\) is called Bessel’s correction.
In R these are var() and sd(), both of which already use \(n - 1\).
Check the theory against the simulation. With \(n = 10\) the predicted shortfall is \((n-1)/n = 0.9\), so dividing by \(n\) should return \(0.9 \times 2.9167 = 2.625\). The simulation returned 2.6297.
When the population mean is known
Occasionally \(\mu\) is known even though \(\sigma^2\) is not — a machine calibrated to a specification, or a process with a target value. Then no degrees of freedom are spent estimating the centre, nothing is flattered, and the divisor is \(n\):
\[\widehat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i - \mu)^2\]
Note that this uses \(\mu\), not \(\bar{x}\). Since \(\bar{x}\) minimises the sum of squared deviations, the total measured from \(\mu\) is larger than the total measured from \(\bar{x}\) — which is precisely why it can be divided by the larger \(n\) and still come out right.
The rule is easier to remember as a question than as two formulas. How many quantities did I have to estimate from this sample before I could measure its spread?
Estimated the centre from the data: divide by \(n - 1\). Knew the centre already: divide by \(n\).
Back to the standard error
We can now return to where the section began. The standard error of the sample mean is \(\sigma/\sqrt{n}\), and with \(\sigma\) unknown we substitute \(s\):
\[\widehat{\mathrm{se}}(\bar{X}) = \frac{s}{\sqrt{n}}\]
This is the quantity reported by every statistical package. Note the hat: it is itself an estimate, depending on the random sample through \(s\), and therefore carrying sampling variability of its own.
That substitution looks harmless. It is not, and Section 3.5 measures the damage.
\(s^2\) is unbiased for \(\sigma^2\), but \(s\) is not unbiased for \(\sigma\). Taking a square root does not preserve unbiasedness.
The evidence is already in the table above. Look at the avg_sd column: the
\(n-1\) row gives 1.6876 against a true \(\sigma\) of
1.7078 — close, but low, and it stays low however many
samples are run. On that column no divisor in the table is right; the best
would fall between \(n-1\) and \(n-2\), and would depend on the distribution.
The bias is small, shrinks quickly with \(n\), and is not worth correcting in practice. What is worth taking away is that unbiasedness attaches to a specific estimator of a specific parameter, and does not survive being passed through a nonlinear function.