(Why) Are stepwise regression coefficients biased?

Question

Issues with stepwise regression are known to statisticians.

1. It yields R-squared values that are badly biased to be high.
2. The F and chi-squared tests quoted next to each variable on the printout do not have the claimed distribution.
3. The method yields confidence intervals for effects and predicted values that are falsely narrow; see Altman and Andersen (1989).
4. It yields p-values that do not have the proper meaning, and the proper correction for them is a difficult problem.
5. It gives biased regression coefficients that need shrinkage (the coefficients for remaining variables are too large; see Tibshirani [1996]).
6. It has severe problems in the presence of collinearity.
7. It is based on methods (e.g., F tests for nested models) that were intended to be used to test prespecified hypotheses.
8. Increasing the sample size does not help very much; see Derksen and Keselman (1992).
9. It allows us to not think about the problem.
10. 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.

Answer 1

TL;DR: Yes, it does return biased coefficients.

I believe this is incorrect, if we are using the word "bias" in the way that people talk about it with regard to statistical and machine learning. (EDIT: Actually it IS correct, just a little confusing, explained below in the edit. Leaving this here so others can track my thought process, since it is a little tricky—for me at least!)

It gives biased regression coefficients that need shrinkage (the coefficients for remaining variables are too large; see Tibshirani [1996]).

I think they are using "biased" in the colloquial sense of "not good." Doing stepwise regression for variable selection is going to harm your inference, as p-values are going to be a mess since you "cherry-picked" the best variables. (There's a whole bunch of literature on this on "post-selection inference" specifically with regards to variable selection methods like stepwise or the lasso).

The problem with stepwise is that the coefficients are not biased—they will have high variance and will be sensitive to changing from different testing sets. In the long-run, like you find in your simulations, they will converge to the "true" coefficient and are thus unbiased.

But in the real world, we don't get to do many many studies, we don't know the "true" coefficient, and our model will always be an oversimplification, never fully explaining the true data generating procedure. So we actually introduce bias.

It's weird, because shrinkage methods are a way of introducing bias in an attempt to decrease variance. I wouldn't read too much into it—I really think they're just using the word "biased" in a colloquial sense instead of a statistical or machine learning sense.

I re-read Tibshirani (1996), and stepwise methods aren't included in their simulations. (Apparently this is attributed to Frank Harrell originally, who is on this site, and, uh, exponentially smarter than I am, so I'll see if I can find anything there.)

---

EDIT: OK, so I think I know where this argument comes from, after reading DOI: 10.1016/j.jclinepi.2011.06.016 and DOI: 10.1016/s0895-4356(99)00103-1.

The problem here is that the p-values are used for selection, and multiple regression controls for other variables, but only the variables in the model. So, if we use these invalid p-values and standard errors to estimate the model, we are letting these define what variables are included in the model. So, we end up overestimating the coefficients (biased away from zero), because the estimation of that coefficient is dependent on it's p-value and standard error.

That first reference I pointed to above (10.1016/j.jclinepi.2011.06.016) makes a pretty good analogy:

An intuitive example includes publication bias, where trials with statistically significant results have a higher chance of being published, leading to too extreme effect estimates in summaries of published reports.

And they do show simulation studies where stepwise biases the coefficients by estimating them to be too large.

Think of it this way: Stepwise methods are only going to include the variables on the occasions when they suck up enough variance to be significant. This is going to bias the coefficients upwards, because, in the long-run, we aren't even considering the non-zero coefficients that do not happen to be significant. It is truncating the possible distribution of coefficients we could get, so it moves the mean away from zero. That's bias.

I was confused, because often the goal is not to estimate unbiased coefficients. It's for good predictions. Indeed, we want to bias our coefficients toward zero (lasso) or toward the mean (mixed effects models) to decrease variance.

---

ANOTHER EDIT: Since I wasn't happy with my flip-flopping, I did my own simulation:

library(tidyverse)
set.seed(1839)             # set seed for randomization
n <- 100                   # sample size
B <- c(0, .25, rep(0, 9))  # intercept 0, beta1 = .25, betas 2 through 9 = 0
eps_sd <- 1                # sd of the error
iter <- 10000              # do it 10000 times

res <- map(seq_len(iter), \(zzz) {
  # make data:
  X <- cbind(1, replicate(length(B) - 1, rnorm(n, 0, 1)))
  y <- X %*% B + rnorm(n, 0, eps_sd)
  dat <- as.data.frame(X[, -1])
  dat <- cbind(dat, y)
  
  # do stepwise:
  m0 <- lm(y ~ 1, dat) # empty model
  m1 <- lm(y ~ ., dat) # saturated model
  fw <- step(m0, direction = "forward", scope = formula(m1), trace = 0)
  
  # return coefficient for V1 if it was estimated, NA otherwise:
  return(ifelse("V1" %in% names(coef(fw)), coef(fw)[["V1"]], NA))
})

# check it
res %>% 
  unlist() %>% 
  summary()

That should show you that the mean coefficient is .28, which is larger than .25.

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
-0.1947  0.2133  0.2699  0.2793  0.3357  0.7255    1435 

The main thing you might be running into with your simulation: If the effect size is big enough (the power is high enough) that you just so happen to estimate the coefficient on every single iteration, then it will be unbiased—because there is no "selection bias," as they were all selected. What you need to do is make the sample size low enough or the effect weak enough so that the model does NOT estimate it every time, meaning there is a selection bias toward coefficients that are sufficiently away from zero so as to be significant.

(Note that I have 1435 NA values where it was not estimated—this is where the selection bias comes from: We should be including those in our calculation of the mean coefficient, but we're not because forward stepwise introduced selection bias.)

Note that my simulation is one coefficient of .25 and the rest of 0—all of them uncorrelated. Things might get funkier when there are correlated variables.