Automatic Model Selection (Optional)

NoteOptional topic

The main course compares a small set of scientifically justified models and checks predictive performance. Automatic searches and additional selection criteria are supplementary. A search-selected model does not inherit the usual fixed-model interpretation of its p-values and confidence intervals.

Forward Selection

In forward selection, we begin with a small model, often the intercept-only model, and add predictors one at a time. At each step, we add the variable that gives the best improvement according to some criterion, such as:

  • partial \(F\) test;
  • smallest \(p\)-value;
  • largest drop in AIC;
  • largest increase in adjusted \(R^2\).

This continues until no candidate addition meets the chosen rule. Backward Elimination In backward elimination, we begin with a larger model and remove predictors one at a time. At each step, we remove the least useful variable according to a chosen rule. This continues until all remaining variables satisfy the stopping criterion. Backward elimination requires that the initial model be estimable, and it may be sensitive to collinearity and hierarchical structure. Stepwise Selection

Stepwise selection combines forward and backward ideas. A variable may enter at one stage and later be removed if its contribution becomes weak after other variables enter the model. This procedure is popular in software because it is automated, but it should be used cautiously. Automatic procedures can be unstable and may encourage overly mechanical modelling. Problems With Automatic Selection

Automatic variable selection can create several difficulties:

  • it ignores model uncertainty;
  • repeated searching inflates the chance of false discoveries;
  • selected coefficients and \(p\)-values may look more certain than they really are;
  • different but similar datasets may lead to different selected models;
  • scientific structure may be lost if the procedure is used blindly.

So automated methods should be treated as exploratory tools, not as final arbiters of truth.

Additional comparison criteria

Criteria for Comparing Models

Several criteria are commonly used in model comparison.

Adjusted R Squared

Adjusted \(R^2\) is

\[ R^2_{\mathrm{adj}} = 1 - \frac{\mathrm{SSE}/(n-p)}{\mathrm{SST}/(n-1)}. \]

Unlike ordinary \(R^2\), adjusted \(R^2\) can decrease when unnecessary predictors are added. It rewards fit but penalizes unnecessary complexity in a simple way.

AIC

The Akaike information criterion is typically written as

\[ \mathrm{AIC} = -2\log L + 2k, \]

where \(L\) is the fitted likelihood and \(k\) is the number of estimated parameters. Smaller AIC indicates a better tradeoff between fit and complexity. AIC is often more prediction-oriented than strict parsimony-oriented.

BIC

The Bayesian information criterion is

\[ \mathrm{BIC} = -2\log L + k\log n. \]

Because the penalty term is stronger than that of AIC when \(n\) is moderate or large, BIC tends to favor smaller models.

Mallows’ Cp

Another classical criterion is Mallows’ \(C_p\), which compares model bias and variance relative to a fuller model. It is less frequently emphasized in introductory software workflows today, but it remains conceptually important in regression theory.

An executable selection example

Generate data with multicollinearity

Show R code 
set.seed(8561)
n <- 80
x1 <- rnorm(n)
x2 <- 0.92 * x1 + rnorm(n, sd = 0.25)
x3 <- rnorm(n)
y <- 3 + 2 * x1 - 1.5 * x2 + 1.2 * x3 + rnorm(n, sd = 1)

dat <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3)

fit_full <- lm(y ~ x1 + x2 + x3, data = dat)
summary(fit_full)

Call:
lm(formula = y ~ x1 + x2 + x3, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max
-2.53882 -0.85013 -0.08111  0.85767  2.24530

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0750     0.1281  24.008  < 2e-16 ***
x1            2.0710     0.4768   4.343 4.28e-05 ***
x2           -1.6962     0.5124  -3.310  0.00143 **
x3            1.1336     0.1277   8.880 2.27e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.117 on 76 degrees of freedom
Multiple R-squared:  0.5981,    Adjusted R-squared:  0.5822
F-statistic:  37.7 on 3 and 76 DF,  p-value: 4.966e-15

Compare with simpler models

Show R code 
fit_small <- lm(y ~ x1 + x3, data = dat)
fit_alt <- lm(y ~ x2 + x3, data = dat)

summary(fit_small)

Call:
lm(formula = y ~ x1 + x3, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max
-2.7697 -0.6645 -0.1148  0.9949  2.3101

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.1177     0.1354   23.02  < 2e-16 ***
x1            0.5343     0.1159    4.61 1.57e-05 ***
x3            1.1421     0.1356    8.42 1.58e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.187 on 77 degrees of freedom
Multiple R-squared:  0.5401,    Adjusted R-squared:  0.5282
F-statistic: 45.22 on 2 and 77 DF,  p-value: 1.027e-13
Show R code 
summary(fit_alt)

Call:
lm(formula = y ~ x2 + x3, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max
-2.7455 -0.7881 -0.1158  1.0147  2.5958

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.1185     0.1417  22.003  < 2e-16 ***
x2            0.4702     0.1301   3.615 0.000534 ***
x3            1.1415     0.1417   8.057 7.91e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.24 on 77 degrees of freedom
Multiple R-squared:  0.4983,    Adjusted R-squared:  0.4853
F-statistic: 38.24 on 2 and 77 DF,  p-value: 2.927e-12
Show R code 
anova(fit_small, fit_full)
Analysis of Variance Table

Model 1: y ~ x1 + x3
Model 2: y ~ x1 + x2 + x3
  Res.Df     RSS Df Sum of Sq      F   Pr(>F)
1     77 108.498
2     76  94.824  1    13.673 10.959 0.001427 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Use AIC-based stepwise selection

Show R code 
fit_null <- lm(y ~ 1, data = dat)
fit_scope <- lm(y ~ x1 + x2 + x3, data = dat)

step_forward <- step(fit_null,
                     scope = list(lower = fit_null, upper = fit_scope),
                     direction = "forward",
                     trace = 0)

step_backward <- step(fit_scope, direction = "backward", trace = 0)

summary(step_forward)

Call:
lm(formula = y ~ x3 + x1 + x2, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max
-2.53882 -0.85013 -0.08111  0.85767  2.24530

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0750     0.1281  24.008  < 2e-16 ***
x3            1.1336     0.1277   8.880 2.27e-13 ***
x1            2.0710     0.4768   4.343 4.28e-05 ***
x2           -1.6962     0.5124  -3.310  0.00143 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.117 on 76 degrees of freedom
Multiple R-squared:  0.5981,    Adjusted R-squared:  0.5822
F-statistic:  37.7 on 3 and 76 DF,  p-value: 4.966e-15
Show R code 
summary(step_backward)

Call:
lm(formula = y ~ x1 + x2 + x3, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max
-2.53882 -0.85013 -0.08111  0.85767  2.24530

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0750     0.1281  24.008  < 2e-16 ***
x1            2.0710     0.4768   4.343 4.28e-05 ***
x2           -1.6962     0.5124  -3.310  0.00143 **
x3            1.1336     0.1277   8.880 2.27e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.117 on 76 degrees of freedom
Multiple R-squared:  0.5981,    Adjusted R-squared:  0.5822
F-statistic:  37.7 on 3 and 76 DF,  p-value: 4.966e-15

Compare AIC, BIC, and adjusted R squared

Show R code 
model_summary <- function(model) {
  c(
    AIC = AIC(model),
    BIC = BIC(model),
    adj_R2 = summary(model)$adj.r.squared
  )
}

rbind(
  small = model_summary(fit_small),
  alt = model_summary(fit_alt),
  full = model_summary(fit_full),
  step_forward = model_summary(step_forward),
  step_backward = model_summary(step_backward)
)
                   AIC      BIC    adj_R2
small         259.4064 268.9345 0.5281786
alt           266.3670 275.8951 0.4852878
full          250.6302 262.5403 0.5822139
step_forward  250.6302 262.5403 0.5822139
step_backward 250.6302 262.5403 0.5822139

Interpretation. Compare the selected variable sets, but do not treat this demonstration as evidence that an automatic search is scientifically optimal. If evaluating a selection procedure by resampling, run the entire selection procedure inside each training fold. Ordinary prediction validation does not repair post-selection confidence intervals.