Issues with stepwise regression are known to statisticians.
- It yields R-squared values that are badly biased to be high.
- The F and chi-squared tests quoted next to each variable on the printout do not have the claimed distribution.
- The method yields confidence intervals for effects and predicted values that are falsely narrow; see Altman and Andersen (1989).
- It yields p-values that do not have the proper meaning, and the proper correction for them is a difficult problem.
- It gives biased regression coefficients that need shrinkage (the coefficients for remaining variables are too large; see Tibshirani [1996]).
- It has severe problems in the presence of collinearity.
- It is based on methods (e.g., F tests for nested models) that were intended to be used to test prespecified hypotheses.
- Increasing the sample size does not help very much; see Derksen and Keselman (1992).
- It allows us to not think about the problem.
- It uses a lot of paper.
I want to focus on #5. I cannot access the Tibshirani paper, but I have run simulations that seem to show the parameters to be unbiased. I can believe that. Once we do the stepwise procedure and select the features to include in the model, we fit a linear regression. Under common assumptions, the Gauss-Markov theorem takes over and gives us unbiased estimates of the coefficients.
Why doesn't Gauss-Markov assure us of getting unbiased coefficient estimates?
(I realize that #2, #3, and #4 are huge problems for inference, whether the coefficients are biased or not, but I want to focus on the bias part for now.)
EDIT
Simulation code in R
library(MASS)
library(ggplot2)
set.seed(2021)
N <- 100
B <- 1000
x1 <- rep(c(0, 1), 2 * N)
x2 <- rep(c(0, 0, 1, 1), N)
x12 <- x1 * x2
y_hat <- 1 + x1 + 0 * x2 - x12
b0_hat <- b1_hat <- b2_hat <- b12_hat <- rep(0, B)
# Define function (from Hadley Wickham) to suppress function printouts
# https://stackoverflow.com/a/54136863
#
quiet <- function(x) {
sink(tempfile())
on.exit(sink())
invisible(force(x))
}
for (i in 1:B){
y <- y_hat + rnorm(length(y_hat), 0, 0.25)
d <- data.frame(x1 = x1, x2 = x2, x12 = x12, y = y)
L <- lm(y ~ ., data = d)
aic <- quiet(MASS::stepAIC(L))
b0_hat[i] <- data.frame(t(aic$coefficients))$X.Intercept.
if (length(data.frame(t(aic$coefficients))$x1) > 0){
b1_hat[i] <- data.frame(t(aic$coefficients))$x1
}
if (length(data.frame(t(aic$coefficients))$x2) > 0){
b2_hat[i] <- data.frame(t(aic$coefficients))$x2
}
if (length(data.frame(t(aic$coefficients))$x12) > 0){
b12_hat[i] <- data.frame(t(aic$coefficients))$x12
}
if (i < 6 | i %% 50 == 0){
print(paste(i/B*100, "% complete", sep = ""))
}
}
d0 <- data.frame(estimate = b0_hat, coefficient = "Intercept")
d1 <- data.frame(estimate = b1_hat, coefficient = "beta1")
d2 <- data.frame(estimate = b2_hat, coefficient = "beta2")
d12 <- data.frame(estimate = b12_hat, coefficient = "beta12")
d <- rbind(d0, d1, d2, d12)
ggplot(d, aes(x = estimate, fill = coefficient)) +
geom_density() +
facet_grid(rows = vars(coefficient)) + #
theme_bw()
t.test(b0_hat, mu = 1)
t.test(b1_hat, mu = 1)
t.test(b2_hat, mu = 0)
t.test(b12_hat, mu = -1)
Perhaps a key takeaway is that I would be including interaction terms, so even if the experimental design is orthogonal (I don't think mine will be, but maybe), the interactions will introduce correlation.
Nonetheless, those coefficient estimates sure don't look biased, and their $95\%$ confidence intervals contain the true parameter values (ditto when I do $50000$ iterations). With $1000$ iterations of the simulation, that should give us good power to reject pretty small differences.