1 Motivation

Students (and, sadly, many practitioners) are taught to chase an \(R^2\) as close to one as possible. This note shows, with a single and very familiar data set, that this is the wrong target. As we make a linear regression more and more flexible:

  • every in-sample measure of fit improves monotonically: \(R^2 \uparrow\), \(\text{MSE} \downarrow\), \(\text{RMSE} \downarrow\), \(\text{MAD} \downarrow\);
  • their out-of-sample counterparts eventually go the other way: \(\text{MSE}_{\text{out}} \uparrow\), \(\text{RMSE}_{\text{out}} \uparrow\), \(\text{MAD}_{\text{out}} \uparrow\).

The model with \(R^2 \approx 1\) will turn out to be, by far, the worst predictive model of the lot.

2 Data and the train/test split

We use the Boston housing data (MASS::Boston): \(n = 506\) census tracts, response medv (median house value, in $1000s) and \(d = 13\) predictors, all standardized. Half of the tracts are used for estimation (training) and the other half are kept aside to evaluate predictions (test).

library(tidyverse)
library(rsample)

set.seed(42)

db <- MASS::Boston |>
  select(-medv) |>
  scale() |>
  as_tibble()

db$medv <- MASS::Boston$medv

split <- initial_split(db, prop = 0.5)
trn   <- training(split)
tst   <- testing(split)

d <- ncol(db) - 1
c(n_train = nrow(trn), n_test = nrow(tst), d = d)
## n_train  n_test       d 
##     253     253      13

3 A nested family of increasingly flexible models

For \(k = 1, 2, \ldots, 5\), model \(M_k\) contains an intercept, all main effects and all interactions of up to \(k\) distinct predictors. In R this is the formula medv ~ .^k (for \(k=1\) it is simply medv ~ .). For instance, \(M_2\) is

\[ y_i = \beta_0 + \sum_{j=1}^{d} \beta_j x_{ij} + \sum_{j<l} \beta_{jl}\, x_{ij} x_{il} + \varepsilon_i , \]

and, in general, the number of regression coefficients in \(M_k\) is

\[ p_k \;=\; 1 + \sum_{j=1}^{k} \binom{d}{j}. \]

The models are nested, \(M_1 \subset M_2 \subset \cdots \subset M_5\).

orders <- 1:5
p_k    <- 1 + cumsum(choose(d, orders))
tibble(k = orders, p_k = p_k, n_train = nrow(trn), p_over_n = round(p_k / nrow(trn), 2))
## # A tibble: 5 × 4
##       k   p_k n_train p_over_n
##   <int> <dbl>   <int>    <dbl>
## 1     1    14     253     0.06
## 2     2    92     253     0.36
## 3     3   378     253     1.49
## 4     4  1093     253     4.32
## 5     5  2380     253     9.41

With \(n = 253\) training observations, \(M_1\) and \(M_2\) have \(p < n\), while from \(k = 3\) onwards we have \(p > n\): there are more coefficients than data points.

4 Estimation: (almost) ordinary least squares

When \(p > n\), \(X^\top X\) is singular and OLS is not unique. We therefore use a ridge estimator with a tiny penalty \(\lambda = 0.001\) (the intercept is not penalized):

\[ \widehat\beta_\lambda \;=\; \left(X^\top X + \lambda D\right)^{-1} X^\top y, \qquad D = \operatorname{diag}(0, 1, \ldots, 1). \]

For \(p < n\) this is numerically indistinguishable from OLS. For \(p > n\) and \(\lambda \to 0\) it converges to the minimum-norm least-squares solution, which interpolates the training data, i.e. \(X\widehat\beta = y\). The penalty is there only to make the linear system solvable, not to regularize in any meaningful way.

5 Measures of fit

Let \(\widehat y_i\) be the fitted value of observation \(i\). On the training set \(T\) (size \(n\)) and on the test set \(\mathcal{V}\) (size \(m\)) we compute

\[ R^2_{\text{in}} = \operatorname{cor}^2(y, \widehat y), \qquad \text{RMSE}_{\text{in}} = \sqrt{\frac{1}{n}\sum_{i\in T} (y_i-\widehat y_i)^2}, \qquad \text{MAD}_{\text{in}} = \frac{1}{n}\sum_{i\in T} |y_i-\widehat y_i|, \]

\[ R^2_{\text{out}} = \operatorname{cor}^2(y^{\ast}, \widehat y^{\ast}), \qquad \text{RMSE}_{\text{out}} = \sqrt{\frac{1}{m}\sum_{i\in V} (y^\ast_i-\widehat y^\ast_i)^2}, \qquad \text{MAD}_{\text{out}} = \frac{1}{m}\sum_{i\in V} |y^\ast_i-\widehat y^\ast_i|, \]

where \(\widehat y^\ast_i = x_i^{\ast\top}\widehat\beta\) uses coefficients estimated only on the training set. (MSE is simply RMSE squared.)

6 Fitting the five models

fit_order <- function(k, trn, tst, lambda = 0.001) {

  model_formula <- if (k == 1) "medv ~ ." else sprintf("medv ~ .^%d", k)

  # training
  X <- model.matrix(formula(model_formula), data = trn)
  y <- trn$medv

  D <- diag(ncol(X))
  D[1, 1] <- 0

  beta_hat  <- solve(crossprod(X) + lambda * D, crossprod(X, y)) # (near-)OLS via ridge
  y_hat_trn <- drop(X %*% beta_hat)

  # testing
  X_tst <- model.matrix(formula(model_formula), data = tst)
  y_tst <- tst$medv
  y_hat <- drop(X_tst %*% beta_hat)

  tibble(
    k        = k,
    formula  = model_formula,
    p        = ncol(X),
    R2_in    = cor(y, y_hat_trn)^2,
    RMSE_in  = sqrt(mean((y - y_hat_trn)^2)),
    MAD_in   = mean(abs(y - y_hat_trn)),
    R2_out   = cor(y_tst, y_hat)^2,
    RMSE_out = sqrt(mean((y_tst - y_hat)^2)),
    MAD_out  = mean(abs(y_tst - y_hat))
  )
}

results <- map_dfr(orders, fit_order, trn = trn, tst = tst)
results |>
  mutate(across(starts_with("R2"), ~ round(.x, 4)),
         across(c(starts_with("RMSE"), starts_with("MAD")), ~ round(.x, 2))) |>
  knitr::kable(caption = "In-sample versus out-of-sample performance, k = 1, ..., 5")
In-sample versus out-of-sample performance, k = 1, …, 5
k formula p R2_in RMSE_in MAD_in R2_out RMSE_out MAD_out
1 medv ~ . 14 0.7303 4.97 3.51 0.7469 4.48 3.41
2 medv ~ .^2 92 0.9307 2.52 1.82 0.5854 6.83 4.06
3 medv ~ .^3 378 0.9946 0.70 0.36 0.1407 25.41 13.66
4 medv ~ .^4 1093 0.9998 0.13 0.05 0.1381 36.82 16.02
5 medv ~ .^5 2380 1.0000 0.03 0.01 0.0519 47.13 16.71 
results |>
  select(k, p, RMSE_in, RMSE_out, MAD_in, MAD_out) |>
  pivot_longer(-c(k, p),
               names_to  = c("measure", "sample"),
               names_sep = "_") |>
  mutate(sample = factor(sample, levels = c("in", "out"),
                         labels = c("In-sample (training)", "Out-of-sample (test)"))) |>
  ggplot(aes(k, value, colour = sample)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  facet_wrap(~ measure, scales = "free_y") +
  scale_y_log10() +
  scale_x_continuous(breaks = orders,
                     labels = paste0("k=", orders, "\np=", results$p)) +
  labs(x = "Highest interaction order (number of coefficients)",
       y = "Error (log scale)", colour = NULL,
       title = "In-sample error keeps falling; out-of-sample error keeps rising") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

results |>
  select(k, R2_in, R2_out) |>
  pivot_longer(-k, names_to = "sample", values_to = "R2") |>
  mutate(sample = recode(sample, R2_in = "In-sample R^2", R2_out = "Out-of-sample R^2")) |>
  ggplot(aes(k, R2, colour = sample)) +
  geom_hline(yintercept = 1, linetype = 2, colour = "grey50") +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "Highest interaction order k", y = expression(R^2), colour = NULL,
       title = "The 'holy grail' R^2 = 1 is reached exactly where prediction collapses") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

7 Why does this happen?

7.1 In-sample fit can only improve

Because the models are nested, the least-squares problem for \(M_{k+1}\) is a minimization over a larger set than that for \(M_k\) (just set the extra coefficients to zero to recover \(M_k\)). Hence

\[ \text{RSS}_{k+1} \;=\; \min_{\beta \in \mathbb{R}^{p_{k+1}}} \|y - X_{k+1}\beta\|^2 \;\le\; \min_{\beta \in \mathbb{R}^{p_{k}}} \|y - X_{k}\beta\|^2 \;=\; \text{RSS}_{k}, \]

and since \(R^2 = 1 - \text{RSS}/\text{TSS}\) (under OLS with an intercept, this equals \(\operatorname{cor}^2(y,\widehat y)\)), \(R^2\) is non-decreasing in model size. Adding any regressor, even pure noise, cannot lower \(R^2\). In the extreme case \(p \ge n\) with \(\operatorname{rank}(X) = n\), the column space of \(X\) is all of \(\mathbb{R}^n\), so \(\widehat y = y\), \(\text{RSS} = 0\) and

\[ R^2_{\text{in}} = 1, \qquad \text{RMSE}_{\text{in}} = \text{MAD}_{\text{in}} = 0 . \]

The tiny ridge penalty is the only reason the table shows \(R^2_{\text{in}}\) slightly below one and in-sample errors slightly above zero. A “perfect fit” here says nothing about the data and everything about the number of parameters.

7.2 Out-of-sample error is a bias–variance trade-off

Write \(y = f(x) + \varepsilon\) with \(\operatorname{E}(\varepsilon)=0\) and \(\operatorname{V}(\varepsilon)=\sigma^2\). For a new point \(x^\ast\), the expected squared prediction error of \(\widehat f\) decomposes as

\[ \operatorname{E}\big[(y^\ast - \widehat f(x^\ast))^2\big] \;=\; \underbrace{\sigma^2}_{\text{irreducible}} \;+\; \underbrace{\big(\operatorname{E}\widehat f(x^\ast) - f(x^\ast)\big)^2}_{\text{bias}^2} \;+\; \underbrace{\operatorname{V}\big(\widehat f(x^\ast)\big)}_{\text{variance}} . \]

Richer models reduce the bias but increase the variance. For OLS with a correct linear model, \(p < n\) and Gaussian errors, the classical results make the contrast explicit:

\[ \operatorname{E}\big[\text{MSE}_{\text{in}}\big] \;=\; \sigma^2\left(1 - \frac{p}{n}\right) \qquad\text{whereas}\qquad \operatorname{E}\big[\text{MSE}_{\text{out}}\big] \;\approx\; \sigma^2\left(1 + \frac{p}{n}\right). \]

The in-sample error is optimistic by roughly \(2\sigma^2 p/n\), and the gap widens linearly with \(p\). As \(p \to n\) the in-sample MSE goes to zero while the out-of-sample MSE grows (for random designs it actually blows up like \(\sigma^2 (n-1)/(n-p-1)\)). Every extra coefficient is fitting noise that will not be there in the next sample.

7.3 Reading the table

  • \(M_1\) (14 coefficients): honest model; in-sample and out-of-sample RMSE are of the same order (4.97 vs 4.48).
  • \(M_2\) (92 coefficients): \(R^2_{\text{in}}\) jumps to 0.931, but the test RMSE already worsens to 6.83.
  • \(M_3\) to \(M_5\) (\(p > n\)): the training data are essentially interpolated (\(R^2_{\text{in}} \to 1\), in-sample RMSE 0.034 for \(k=5\)) while the test RMSE explodes to 47.1, roughly 11 times worse than the simplest model, and the out-of-sample \(R^2\) drops to 0.052.

The same story is told by RMSE and MAD alike: it is not an artifact of a particular loss function.

8 Take-away messages

  1. In-sample \(R^2\), MSE, RMSE, MAD, etc. are monotone in model complexity for nested least-squares models. They cannot be used to choose among them.
  2. \(R^2 \approx 1\) is not evidence of a good model; with \(p \ge n\) it is guaranteed and therefore carries no information.
  3. Models must be judged on data they have not seen: hold-out samples, cross-validation, or criteria that penalize complexity (adjusted \(R^2\), AIC, BIC, marginal likelihoods / Bayes factors, LOO-CV, WAIC).
  4. Flexibility is not the enemy per se: what is missing above is genuine regularization or shrinkage (a properly tuned ridge/lasso penalty, or priors that shrink interactions towards zero), which trades a little bias for a large reduction in variance.

9 Appendix: counting coefficients

# number of non-intercept columns of the k = 5 design matrix
sum(choose(d, 1:5))
## [1] 2379
results$p[5] - 1
## [1] 2379