Appendix B — F, Wald, and likelihood ratio tests

Chapter 9 showed that the F test and the likelihood ratio test can give different p-values for the same comparison, and explained why: \(\chi^2\) is what the \(F\) distribution becomes as the denominator degrees of freedom go to infinity, so it has slightly less mass in its tail.

That was one example. This appendix looks at the pattern systematically, and adds a third statistic (the Wald test) that you will meet in the output of most generalized linear models.

B.1 Three ways to test the same thing

Consider the simplest possible comparison: an intercept-only model against one with a single predictor. Three statistics are available.

The F statistic compares reductions in sum of squared error, scaled by the error that remains:

\[F = \frac{(\text{SSE}_C - \text{SSE}_A)/1}{\text{SSE}_A/(n-2)}\]

referred to an \(F\) distribution with 1 and \(n - 2\) degrees of freedom.

The Wald statistic takes that same quantity and refers it to a \(\chi^2\) distribution instead, ignoring the fact that the denominator was estimated.

The likelihood ratio statistic compares log-likelihoods, which for the normal model reduces to a function of the ratio of the two error sums:

\[G^2 = n \log\left(\frac{\text{SSE}_C}{\text{SSE}_A}\right)\]

also referred to \(\chi^2\) with one degree of freedom.

All three test the same null hypothesis. None of them agrees exactly with the others.

B.2 Watching them converge

We can hold the association constant and vary only the sample size. MASS::mvrnorm() with empirical = TRUE generates data whose sample correlation is exactly what we ask for, which removes sampling noise and isolates the effect of \(n\).

# create dataset with exact known correlation
# then extract F, G^2, and the F, Wald, and LRT chi-square p-values

get_stats <- function(n, r) {

  Sigma <- matrix(c(1, r,
                    r, 1), 2, 2)

  mat <- MASS::mvrnorm(n, mu = c(0, 0), Sigma = Sigma, empirical = TRUE)

  x = mat[,1]
  y = mat[,2]

  m0 <- lm(y ~ 1)
  m1 <- lm(y ~ x)

  rss0 <- deviance(m0)
  rss1 <- deviance(m1)

  f_stat <- ((rss0 - rss1) / 1) / (rss1 / (n-2))
  g2 <- n * log(rss0 / rss1)

  fp <- 1 - pf(f_stat, 1, n-2)
  waldp <- 1 - pchisq(f_stat, 1)
  lrtp <- 1 - pchisq(g2, 1)

  log_fp <- log(fp)
  log_waldp <- log(waldp)
  log_lrtp <- log(lrtp)

  tmp <- tibble::tibble(
    f = f_stat,
    g2 = g2,
    fp = fp,
    waldp = waldp,
    lrtp = lrtp,
    log_fp = log_fp,
    log_waldp = log_waldp,
    log_lrtp = log_lrtp
  )

  return(tmp)

}

# for a given PRE (r2), get the data and compare stats for increasing N
mytests <- expand_grid(
  n = 30:80,
  r = sqrt(.1)) |>
  rowwise() |>
  mutate(stats = list(get_stats(n, r))) |>
  unnest_wider(stats)

mytests |>
  select(n, f, g2, fp, waldp, lrtp) |>
  head(4) |>
  mutate(across(everything(), \(z) round(z, 4)))

# A tibble: 4 x 6
      n     f    g2     fp  waldp   lrtp
  <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl>
1    30  3.11  3.16 0.0887 0.0778 0.0754
2    31  3.22  3.27 0.0831 0.0726 0.0707
3    32  3.33  3.37 0.0779 0.0679 0.0663
4    33  3.44  3.48 0.073  0.0635 0.0622

Note that f and g2 are not the same statistic. The F statistic scales the reduction in SSE by the error that remains; \(G^2\) is a function of the ratio of the two error sums. They happen to be close here, and they get closer as \(n\) grows, but they are computed differently and referred to different distributions.

Every one of these datasets has exactly the same association (\(R^2 = 0.1\)), so any difference between the columns is a property of the test, not of the data.

tinyplot
ggplot2

plt(fp ~ n, data = mytests, type = "l", lwd = 2,
    ylim = range(c(mytests$fp, mytests$lrtp)),
    xlab = "Sample size", ylab = "p-value")

lines(mytests$n, mytests$waldp, lty = 2, lwd = 2)
lines(mytests$n, mytests$lrtp,  lty = 3, lwd = 2)

abline(h = 0.05, col = "gray60")

legend("topright", legend = c("F", "Wald", "LRT"),
       lty = 1:3, lwd = 2, bty = "n")

Figure B.1: Three p-values for the same association, as sample size grows.

long <- mytests |>
  select(n, fp, waldp, lrtp) |>
  pivot_longer(-n, names_to = "test", values_to = "p") |>
  mutate(test = factor(test,
                       levels = c("fp", "waldp", "lrtp"),
                       labels = c("F", "Wald", "LRT")))

ggplot(long, aes(x = n, y = p, linetype = test)) +
  geom_line(linewidth = 0.8) +
  geom_hline(yintercept = 0.05, color = "gray60") +
  labs(x = "Sample size", y = "p-value", linetype = "")

Figure B.2: Three p-values for the same association, as sample size grows.

The same three curves are easier to compare on a log scale, which spreads out the small p-values instead of squashing them against zero.

tinyplot
ggplot2

Show code

# plot the results for log p-values

plt(log_fp ~ n, data = mytests, type = "l", lwd = 2,
    ylim = range(c(mytests$log_fp, mytests$log_lrtp)),
    xlab = "Sample size", ylab = "log p-value")

lines(mytests$n, mytests$log_waldp, lty = 2, lwd = 2)
lines(mytests$n, mytests$log_lrtp,  lty = 3, lwd = 2)

abline(h = log(0.05), col = "gray60")

legend("bottomleft", legend = c("F", "Wald", "LRT"),
       lty = 1:3, lwd = 2, bty = "n")

Figure B.3: The same three p-values, logged.

Show code

long_log <- mytests |>
  select(n, log_fp, log_waldp, log_lrtp) |>
  pivot_longer(-n, names_to = "test", values_to = "p") |>
  mutate(test = factor(test,
                       levels = c("log_fp", "log_waldp", "log_lrtp"),
                       labels = c("F", "Wald", "LRT")))

ggplot(long_log, aes(x = n, y = p, linetype = test)) +
  geom_line(linewidth = 0.8) +
  geom_hline(yintercept = log(0.05), color = "gray60") +
  labs(x = "Sample size", y = "log p-value", linetype = "")

Figure B.4: The same three p-values, logged.

Two things are visible.

The F test is the most conservative. Its p-value is always the largest of the three, because it alone accounts for the fact that the denominator of the F ratio was estimated from the data rather than known.

The gap closes as \(n\) grows. By the right edge of the plot the three curves are nearly indistinguishable. All three are asymptotically equivalent; they differ in small samples, where the difference can flip a conclusion.

mytests |>
  filter(n %in% c(30, 40, 60, 80)) |>
  select(n, f, g2, fp, waldp, lrtp) |>
  mutate(across(everything(), \(z) round(z, 4)))

# A tibble: 4 x 6
      n     f    g2     fp  waldp   lrtp
  <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl>
1    30  3.11  3.16 0.0887 0.0778 0.0754
2    40  4.22  4.21 0.0468 0.0399 0.0401
3    60  6.44  6.32 0.0138 0.0111 0.0119
4    80  8.67  8.43 0.0043 0.0032 0.0037

At \(n = 30\) none of the three clears 0.05. At \(n = 40\) all of them do. In between there is a stretch where the choice of test decides the verdict. That is a good reason not to treat 0.05 as the important thing, as Chapter 6 argued.

B.3 Which to use

For the normal linear model, use F. It is exact rather than asymptotic, and its conservatism in small samples is a feature.

For generalized linear models, the F test is unavailable and the choice is between Wald and likelihood ratio. Wald is what appears by default in summary() output because it is cheap. It needs only the fitted model. The likelihood ratio test requires fitting both models, and is generally better behaved, particularly when coefficients are large.

When they disagree noticeably, trust the likelihood ratio test.