17 Model Comparison Criteria

“When a measure becomes a target, it ceases to be a good measure.” - Marilyn Strathern

So far, we have discussed using the coefficient of determination, \(R^2\) and adjusted coefficient of determination, \(R^2_a\), as a means of determining how well the model fit the data and also for comparing different models. We will now discuss other ways to measure fit and compare models. These criteria can be used to help in the model building process.

17.1 Possible Models

From any set of \(p - 1\) predictors, \(2^{p-1}\) alternative models can be constructed.

This calculation is based on the fact that each predictor can be either included or excluded from the model.

Counting possible subset models

If there are 4 possible predictors, each predictor has two choices:

included in the model
excluded from the model

This gives

\[ 2^4 = 16 \]

possible subset models. If there are 10 possible predictors, there are

\[ 2^{10}=1024 \]

possible subset models. The number of possible models grows quickly, which is why model-selection criteria are useful.

For example, the \[ 2^4 = 16 \] different possible subset models that can be formed from the pool of four \(X\) variables are \[ \begin{align*} & \text{None}\\ & x_{1}\\ & x_{2}\\ & x_{3}\\ & x_{4}\\ & x_{1},x_{2}\\ & x_{1},x_{3}\\ & x_{1},x_{4}\\ & x_{2},x_{3}\\ & x_{2},x_{4}\\ & x_{3},x_{4}\\ & x_{1},x_{2},x_{3}\\ & x_{1},x_{2},x_{4}\\ & x_{1},x_{3},x_{4}\\ & x_{2},x_{3},x_{4}\\ & x_{1},x_{2},x_{3},x_{4} \end{align*} \]

Model selection procedures, also known as subset selection or variable selection procedures, have been developed to identify a small group of regression models that are “good” according to a specified criterion.

This limited number might consist of three to six “good” subsets according to the criteria specified, so the investigator can then carefully study these regression models for choosing the final model.

While many criteria for comparing the regression models have been developed, we will focus on six: \[ \begin{align*} & R_{p}^{2}\\ & R_{a,p}^{2}\\ & AIC_{p}\\ & SBC_{p}\\ & PRESS_{p}\\ & \text{pred }R_{p}^{2} \end{align*} \]

Which direction is better?

The criteria do not all move in the same direction.

Criterion	Preferred direction	Main idea
\(R_p^2\)	Larger	Measures in-sample fit, but never decreases when predictors are added.
\(R_{a,p}^2\)	Larger	Adjusts \(R^2\) for the number of parameters in the model.
\(AIC_p\)	Smaller	Rewards fit but penalizes model complexity.
\(BIC_p\) or \(SBC_p\)	Smaller	Similar to AIC, but usually penalizes complexity more strongly.
\(PRESS_p\)	Smaller	Measures leave-one-out prediction error.
pred \(R_p^2\)	Larger	Converts PRESS into an \(R^2\)-style prediction measure.

No single criterion should replace judgment. The best model is usually one that fits the scientific question, predicts well, and is not more complicated than it needs to be.

Exhaustive subset selection with four predictors

When the number of possible predictors is small, we can fit every possible subset model and compare the criteria. This is called exhaustive subset selection.

The following example uses four possible predictors from mtcars: disp, hp, wt, and qsec. Since there are 4 possible predictors, there are \(2^4=16\) subset models.

library(tidyverse)
library(broom)

press_stat <- function(model) {
  sum((residuals(model) / (1 - hatvalues(model)))^2)
}

candidate_predictors <- c("disp", "hp", "wt", "qsec")

all_subsets <- map(
  0:length(candidate_predictors),
  \(k) {
    if (k == 0) {
      list(character(0))
    } else {
      combn(candidate_predictors, k, simplify = FALSE)
    }
  }
) |>
  unlist(recursive = FALSE)

sst_mpg <- sum((mtcars$mpg - mean(mtcars$mpg))^2)

exhaustive_results <- tibble(predictors = all_subsets) |>
  mutate(
    model = map_chr(
      predictors,
      \(vars) {
        if (length(vars) == 0) {
          "Intercept only"
        } else {
          paste(vars, collapse = " + ")
        }
      }
    ),
    formula = map(
      predictors,
      \(vars) {
        right_side <- if (length(vars) == 0) {
          "1"
        } else {
          paste(vars, collapse = " + ")
        }

        as.formula(paste("mpg ~", right_side))
      }
    ),
    fit = map(formula, \(form) lm(form, data = mtcars)),
    p = map_int(predictors, length) + 1,
    PRESS = map_dbl(fit, press_stat),
    pred_r_squared = 1 - PRESS / sst_mpg,
    glance = map(fit, glance)
  ) |>
  unnest(glance) |>
  select(
    model,
    p,
    r.squared,
    adj.r.squared,
    AIC,
    BIC,
    PRESS,
    pred_r_squared
  )

exhaustive_results |>
  arrange(BIC) |>
  knitr::kable(digits = 4)

model	p	r.squared	adj.r.squared	AIC	BIC	PRESS	pred_r_squared
hp + wt	3	0.8268	0.8148	156.6523	162.5153	246.5063	0.7811
wt + qsec	3	0.8264	0.8144	156.7205	162.5834	249.3488	0.7786
hp + wt + qsec	4	0.8348	0.8171	157.1426	164.4713	245.4834	0.7820
disp + hp + wt	4	0.8268	0.8083	158.6430	165.9717	261.3609	0.7679
disp + wt + qsec	4	0.8264	0.8078	158.7203	166.0490	259.1200	0.7699
disp + hp + wt + qsec	5	0.8351	0.8107	159.0696	167.8640	257.8314	0.7710
disp + wt	3	0.7809	0.7658	164.1678	170.0307	309.3015	0.7253
wt	2	0.7528	0.7446	166.0294	170.4266	328.0228	0.7087
disp + hp	3	0.7482	0.7309	168.6186	174.4815	343.9581	0.6945
disp	2	0.7183	0.7090	170.2094	174.6066	365.8296	0.6751
disp + hp + qsec	4	0.7542	0.7279	169.8525	177.1812	356.9407	0.6830
disp + qsec	3	0.7216	0.7024	171.8418	177.7048	378.4050	0.6640
hp	2	0.6024	0.5892	181.2386	185.6358	552.1057	0.5097
hp + qsec	3	0.6369	0.6118	180.3391	186.2020	540.3406	0.5201
qsec	2	0.1753	0.1478	204.5881	208.9853	1041.1210	0.0754
Intercept only	1	0.0000	0.0000	208.7555	211.6870	1199.8671	-0.0656

Sorting by BIC places the simpler models with strong fit near the top. If we sorted by a different criterion, the order might change.

17.2 Notation

Before discussing the criteria, we will need to develop some notation. We will denote the number of potential \(X\) variables in the pool by \(P-1\).

We assume that all regression models contain an intercept term \(\beta_0\).

Hence, the regression function containing all potential \(X\) variables contains \(P\) parameters, and the function with no \(X\) variables contains one parameter (\(\beta_0\)).

The number of \(X\) variables in a subset will be denoted by \(p-1\), as always, so that there are \(p\) parameters in the regression function for this subset of \(X\) variables. Thus, we have: \[ 1\le p \le P < n \]

17.3 Coefficient of Determination and SSE

Clearly, we would like models that fit the data well. Thus we would like a coefficient of multiple determination, \(R^2\), to be high.

We will denote the number of parameters in the potential model as a subscript and write the coefficient of determination as \(R^2_p\).

The \(R^2_p\) criterion is equivalent to using the error sum of squares \(SSE_p\) as the criterion. The \(R^2_p\) criterion is not intended to identify the subsets that maximize this criterion.

We know that \(R^2_p\) can never decrease as additional \(X\) variables are included in the model. Hence, \(R^2_p\) will be a maximum when all \(P - 1\) potential X variables are included in the regression model.

The intent in using the \(R^2_p\) criterion is to find the point where adding more \(X\) variables is not worthwhile because it leads to a very small increase in \(R^2_p\).

Often, this point is reached when only a limited number of \(X\) variables are included in the regression model.

Do not choose a model using R-squared alone

Because \(R_p^2\) cannot decrease when a predictor is added, the model with the most predictors will always look at least as good by this criterion. This makes \(R_p^2\) useful for describing fit, but risky as the only model-selection tool.

When comparing candidate models, look for whether the added variables produce a meaningful improvement, not merely any improvement.

Overfitting by adding noise predictors

The following example adds randomly generated noise predictors to a useful model. These noise variables are not meaningful predictors of mpg, but ordinary least squares can still use them to reduce in-sample error slightly.

set.seed(3386)

press_stat <- function(model) {
  sum((residuals(model) / (1 - hatvalues(model)))^2)
}

noise_dat <- mtcars |>
  mutate(
    noise_1 = rnorm(n()),
    noise_2 = rnorm(n()),
    noise_3 = rnorm(n()),
    noise_4 = rnorm(n())
  )

sst_noise <- sum((noise_dat$mpg - mean(noise_dat$mpg))^2)

noise_models <- list(
  "Base model" = lm(mpg ~ wt + qsec, data = noise_dat),
  "Add 1 noise predictor" = lm(
    mpg ~ wt + qsec + noise_1,
    data = noise_dat
  ),
  "Add 2 noise predictors" = lm(
    mpg ~ wt + qsec + noise_1 + noise_2,
    data = noise_dat
  ),
  "Add 4 noise predictors" = lm(
    mpg ~ wt + qsec + noise_1 + noise_2 + noise_3 + noise_4,
    data = noise_dat
  )
)

noise_comparison <- purrr::imap_dfr(
  noise_models,
  \(fit_object, model_name) {
    glance(fit_object) |>
      mutate(
        model = model_name,
        predictors = length(coef(fit_object)) - 1,
        PRESS = press_stat(fit_object),
        pred_r_squared = 1 - PRESS / sst_noise
      )
  }
) |>
  select(
    model,
    predictors,
    r.squared,
    adj.r.squared,
    AIC,
    BIC,
    PRESS,
    pred_r_squared
  )

noise_comparison |>
  knitr::kable(digits = 4)

model	predictors	r.squared	adj.r.squared	AIC	BIC	PRESS	pred_r_squared
Base model	2	0.8264	0.8144	156.7205	162.5834	249.3488	0.7786
Add 1 noise predictor	3	0.8298	0.8116	158.0855	165.4142	257.2757	0.7715
Add 2 noise predictors	4	0.8298	0.8046	160.0854	168.8798	269.1829	0.7609
Add 4 noise predictors	6	0.8344	0.7946	163.2150	174.9408	327.8328	0.7089

Notice that \(R^2\) increases as the noise predictors are added. That does not mean the model is better. The other criteria warn us that the extra predictors are adding complexity without useful predictive value.

17.4 Adjusted Coefficient of Determination and MSE

Since \(R^2_{p}\) does not take account of the number of parameters in the regression model and since max(\(R^2_{p}\)) can never decrease as \(p\) increases, the adjusted coefficient of multiple determination \(R^2_{a,p}\) has been suggested as an alternative criterion.

This coefficient takes the number of parameters in the regression model into account through the degrees of freedom.

Users of the \(R^2_{a,p}\) criterion seek to find a few subsets for which \(R^2_{a,p}\) is at the maximum or so close to the maximum that adding more variables is not worthwhile.

The \(R^2_{a,p}\) criterion is equivalent to using the mean square error \(MSE_p\) as the criterion. This comes from the relationship \[ \begin{align*} R_{a}^{2} & =1-\frac{\frac{SSE}{n-p}}{\frac{SSTO}{n-1}}\\ & =1-\left(\frac{n-1}{n-p}\right)\frac{SSE}{SSTO}\\ & = 1-\frac{MSE}{\frac{SSTO}{n-1}} \end{align*} \]

17.5 AIC and BIC

We have seen that \(R_{a,p}^2\) is a criterion that penalizes models having large numbers of predictors.

Two popular alternatives that also provide penalties for adding predictors are Akaike’s information criterion (\(AIC_p\)) and Schwarz’s Bayesian criterion (\(SBC_p\)).

A more popular name for \(SBC_p\) is Bayesian information criterion (\(BIC_p\)).

We search for models that have small values of \(AIC_p\), or \(BIC_p\), where these criteria are given by: \[ \begin{align} AIC_{p} & =n\ln SSE_{p}-n\ln n+2p\\ BIC_{p} & =n\ln SSE_{p}-n\ln n+\left(\ln n\right)p \end{align} \tag{17.1}\]

Notice that for both of these measures, the first term is \[ n\ln SSE_{p} \] which decreases as \(p\) increases.

The second term is fixed (for a given sample size \(n\)), and the third term increases with the number of parameters, \(p\).

Models with small \(SSE_p\) will do well by these criteria as long as the penalties \[ \begin{align*} 2p & \text{ for }AIC_p \\ (\ln n)p & \text{ for }BIC_p \end{align*} \] are not too large.

If \(n \ge 8\) the penalty for \(BIC_p\) is larger than that for \(AIC_p\); hence the \(BIC_p\) criterion tends to favor more parsimonious models.

A parsimonious model is a model that accomplishes a desired level of explanation or prediction with as few predictor variables as possible.

We consider the “best” models to be the ones with the lowest \(AIC_p\) or \(BIC_p\).

Compare AIC and BIC only for models fit to the same response

AIC and BIC are useful for comparing candidate models, but the models should be fit to the same response variable using the same set of observations. They are not meant to compare, for example, a model for mpg against a model for log(mpg).

Software may report versions of AIC and BIC that differ by constants from textbook formulas. This is usually not a problem because model selection depends on differences among models fit to the same data.

17.6 Criteria for Prediction

The \(PRESS_p\) (prediction sum of squares) criterion is a measure of how well the use of the fitted values for a subset model can predict the observed responses \(y_i\).

The PRESS measure differs from SSE in that each fitted value \(\hat{y}_i\) for the PRESS criterion is obtained by

deleting the \(i\)th case from the data set
estimating the regression function for the subset model from the remaining \(n - 1\) cases, and
then using the fitted regression function to obtain the predicted value \(\hat{y}_{i(i)}\) for the \(i\)th case.

We use the notation \(\hat{y}_{i(i)}\) now for the fitted value to indicate, by the first subscript \(i\), that it is a predicted value for the \(i\)th case and, by the second subscript \((i)\), that the \(i\)th case was omitted when the regression function was fitted.

The PRESS prediction error for the \(i\)th case then is: \[ y_i - \hat{y}_{i(i)} \] and the \(PRESS_p\) criterion is the sum of the squared prediction errors over all \(n\) cases: \[ \begin{align} PRESS_p = \sum_{i=1}^n\left(y_i - \hat{y}_{i(i)} \right)^2 \end{align} \tag{17.2}\]

Models with small \(PRESS_p\) values are considered good candidate models. The reason is that when the prediction errors \(y_i - \hat{y}_{i(i)}\) are small, so are the squared prediction errors and the sum of the squared prediction errors.

Thus, models with small \(PRESS_p\) values fit well in the sense of having small prediction errors.

Another measure of prediction is the predicted \(R^{2}\). It is related to PRESS by \[ \begin{align} \text{pred }R_{p}^{2} & =1-\frac{PRESS_{p}}{SSTO} \end{align} \tag{17.3}\] Large values of pred \(R_{p}^{2}\) indicate a model that is “good” at prediction.

17.6.1 PRESS and Leave-One-Out Cross-Validation

PRESS is closely connected to leave-one-out cross-validation (LOOCV). In LOOCV, we repeatedly fit the model after leaving out one observation, then predict the observation that was left out. PRESS is the sum of the squared prediction errors from that process.

For ordinary least squares models, PRESS can be computed without refitting the model \(n\) separate times. If \(e_i\) is the ordinary residual from the full model and \(h_{ii}\) is the leverage value for observation \(i\), then

\[ \begin{align*} y_i-\hat{y}_{i(i)} = \frac{e_i}{1-h_{ii}}. \end{align*} \]

So for ordinary least squares,

\[ \begin{align*} PRESS = \sum_{i=1}^n \left(\frac{e_i}{1-h_{ii}}\right)^2. \end{align*} \]

This explains why PRESS is a prediction-focused criterion. Each observation is judged by how well the model predicts it when that observation was not used to fit the model.

Example 17.1

Comparing candidate models with several criteria

We can obtain \(R^2\), \(R^2_a\), AIC, and BIC using the glance() function after fitting a model. We can also calculate PRESS from the residuals and leverage values of an ordinary least squares model.

Let’s compare several candidate models for mtcars. These models are not meant to be the final answer; they are a small set of candidate models used to illustrate how the criteria work.

library(tidyverse)
library(broom)

press_stat <- function(model) {
  sum((residuals(model) / (1 - hatvalues(model)))^2)
}

sst_mpg <- sum((mtcars$mpg - mean(mtcars$mpg))^2)

candidate_models <- list(
  "Full untransformed" = lm(
    mpg ~ disp + hp + drat + wt + qsec,
    data = mtcars
  ),
  "Full transformed" = lm(
    mpg ~ log(disp) + log(hp) + drat + log(wt) + qsec,
    data = mtcars
  ),
  "Smaller transformed" = lm(
    mpg ~ log(hp) + log(wt) + qsec,
    data = mtcars
  ),
  "Weight and time" = lm(
    mpg ~ log(wt) + qsec,
    data = mtcars
  )
)

model_comparison <- purrr::imap_dfr(
  candidate_models,
  \(fit_object, model_name) {
    glance(fit_object) |>
      mutate(
        model = model_name,
        PRESS = press_stat(fit_object),
        pred_r_squared = 1 - PRESS / sst_mpg
      )
  }
) |>
  select(
    model,
    r.squared,
    adj.r.squared,
    AIC,
    BIC,
    PRESS,
    pred_r_squared
  )

model_comparison |>
  arrange(AIC) |>
  knitr::kable(digits = 4)

model	r.squared	adj.r.squared	AIC	BIC	PRESS	pred_r_squared
Smaller transformed	0.8910	0.8793	143.8376	151.1663	156.4152	0.8611
Weight and time	0.8780	0.8696	145.4393	151.3023	167.0926	0.8516
Full transformed	0.8922	0.8715	147.4688	157.7289	171.8721	0.8474
Full untransformed	0.8489	0.8199	158.2784	168.5385	262.7954	0.7666

The table shows that different criteria may prefer different models. A model with more predictors may have a larger \(R^2\), while AIC, BIC, PRESS, or predicted \(R^2\) may favor a simpler model if the extra predictors do not improve the model enough.

What if the criteria disagree?

Model-selection criteria often disagree because they value different things. When that happens, use the disagreement as a prompt for deeper comparison.

Start with the goal. If prediction is the main goal, give more weight to PRESS, predicted \(R^2\), validation error, or cross-validation. If explanation is the main goal, give more weight to interpretability and subject-matter justification.
Prefer simpler models when performance is similar. If two models have nearly identical performance, the more parsimonious model is usually easier to explain, teach, and defend.
Look for unstable variables. If a predictor only appears in the model favored by one criterion, ask whether it has a meaningful interpretation or whether it is mostly taking advantage of this sample.
Check diagnostics. A model with attractive criteria can still violate assumptions. Residual plots, influential observations, and multicollinearity still matter.
Avoid automatic selection. Criteria help narrow the search. They should not replace thinking about the research question, data quality, and consequences of model complexity.

17.7 Recap

In this chapter, we introduced several criteria for comparing candidate regression models.

Idea	Meaning
Subset model	A model formed by including only some of the available predictors.
Number of possible models	With \(p-1\) possible predictors, there are \(2^{p-1}\) possible subset models.
Exhaustive subset selection	Fitting every possible subset model from a chosen predictor pool. This is feasible only when the number of candidate predictors is small.
\(R_p^2\)	Measures in-sample fit for a model with \(p\) parameters; larger is better, but it never decreases when predictors are added.
\(R_{a,p}^2\)	Adjusted \(R^2\); penalizes unnecessary predictors through the degrees of freedom.
AIC	A model-selection criterion that balances fit and complexity; smaller is better.
BIC/SBC	Similar to AIC but usually applies a stronger penalty for complexity; smaller is better.
Parsimonious model	A model that achieves the goal with relatively few predictors.
PRESS	A leave-one-out prediction error criterion; smaller values indicate better predictive performance.
LOOCV	Leave-one-out cross-validation; fit the model leaving out one observation, then predict the omitted observation.
Predicted \(R^2\)	An \(R^2\)-style measure based on PRESS; larger values indicate better prediction.
Overfitting	A model fitting random or sample-specific noise rather than useful signal.
Criterion disagreement	Different criteria may select different models because they reward fit and penalize complexity differently.

17.8 Check your understanding

Problems

Why are there \(2^{p-1}\) possible subset models when there are \(p-1\) candidate predictors?
Why is \(R^2\) risky to use as the only model-selection criterion?
What problem is adjusted \(R^2\) trying to address?
For AIC and BIC, do we prefer larger or smaller values? Why?
Why does BIC often favor simpler models than AIC?
What does it mean for a model to be parsimonious?
How is PRESS different from SSE?
Why is predicted \(R^2\) useful when the goal is prediction?
If one model has the best AIC and another has the best BIC, does that automatically mean one of them is wrong? Explain.
Why should AIC and BIC usually be compared only among models fit to the same response using the same observations?
Why does exhaustive subset selection become difficult as the number of candidate predictors grows?
How is PRESS connected to leave-one-out cross-validation?
Why can adding random noise predictors increase \(R^2\) without improving the model?
If several criteria disagree, why is it usually better to compare a small set of candidate models than to automatically choose the winner from one criterion?

Solutions

Each predictor has two choices. Every candidate predictor can either be included or excluded. With \(p-1\) predictors, this gives \(2^{p-1}\) possible combinations.
It always rewards adding predictors. \(R^2\) cannot decrease when predictors are added, even if the new predictors add little useful information. This can push us toward models that are unnecessarily complex.
It accounts for model size. Adjusted \(R^2\) penalizes models for using more parameters. It tries to distinguish meaningful improvement from improvement that happens merely because more predictors were added.
Smaller values are preferred. AIC and BIC both combine a measure of lack of fit with a penalty for complexity. A smaller value indicates a better tradeoff between fit and model size.
BIC has a stronger penalty when \(n \ge 8\). AIC uses a penalty of \(2p\), while BIC uses \((\ln n)p\). Since \(\ln n\) is larger than 2 for sample sizes of at least 8, BIC more strongly discourages extra predictors.
It is simple but effective. A parsimonious model accomplishes the modeling goal with as few predictors as reasonably possible.
PRESS measures leave-one-out prediction error. SSE uses residuals from a model fit to all observations. PRESS predicts each observation from a model fit without that observation, so it more directly reflects prediction performance.
It puts prediction performance on an \(R^2\)-like scale. Predicted \(R^2\) uses PRESS to summarize how well the model predicts omitted observations. Larger values suggest better predictive performance.
No. Different criteria emphasize different tradeoffs. AIC may prefer a somewhat larger model, while BIC may prefer a simpler one. This disagreement is a reason to examine the models carefully, not a sign that one criterion is invalid.
The criteria are relative comparisons. AIC and BIC are designed to compare candidate models fit to the same data and response. If the response or observations change, the numerical values are no longer directly comparable.
The number of models grows exponentially. Each predictor can be included or excluded, so adding one more candidate predictor doubles the number of possible subset models.
PRESS sums leave-one-out squared prediction errors. Each observation is predicted by a model that was fit without that observation. The squared errors from those omitted-observation predictions are added to form PRESS.
Ordinary least squares can chase sample noise. Adding predictors gives the model more flexibility, so in-sample SSE can decrease and \(R^2\) can increase even when the new predictors have no real meaning. Prediction-focused criteria may reveal that the added complexity does not generalize.
Criteria are tools, not commandments. AIC, BIC, adjusted \(R^2\), and PRESS emphasize different tradeoffs. Comparing a small set of candidate models lets us consider prediction, interpretation, diagnostics, and subject-matter relevance together.