The Problem With Backward Selection

Backward Step-wise selection/elimination algorithm [linear regression case]

1. Fit a full model using all \(N\) \(x_i\)’s
   • \(\hat{y} = \beta_0 + \beta_1 x_1 + ... + \beta_N x_N\)
2. Remove \(x_i\) with the weakest relationship and refit the model.
   • Usually the least significant variable, if one exists, is dropped, and the model refit.
   • \(\hat{y} = \beta_0 + \beta_1 x_1 + \beta_{i-1} x_{i-1} + \beta_{i+1} x_{i+1} + ... + \beta_N x_N\)
3. Repeat 2 until a stopping criteria is reached.
   • Could be when all remaining coefficient p-values are less than the predefined \(\alpha\) value [say, 0.05], or there is only one coefficient remaining.

P-values and their properties

In order to see and understand why step-wise backward selection may not be an ideal method of covariate selection, we must understand a bit more about p-values. We will accomplish this be examining the distribution of p-values under the null hypothesis where \(x\) is independent of \(y\).

First, let’s set our seed and load some dependencies.

# R 3.3.2
library("dplyr")
library("broom")
set.seed(123)

Let’s say we are interested in understanding the relationship between \(x_1\) and \(y\), which under the null hypothesis, \(x_1\) and \(y\) are independent.

Under the null, what would be expect?

Well, lets set our significance/\(\alpha\)-level to 0.05. By definition of \(\alpha\)-level we know that approximately 5% of our results would be significant despite there being no actual relationship present. Furthermore, because there is no real relationship, we would also expect the resulting p-values to be uniform on (0,1).

Let’s take a look at 1,000 simulated situations where the linear relationship between \(x_1\) and \(y\) is compared where both variables are actually independent of one another - I’ll using dplyr for better readability:

x <- vector()
sim1 <- 1000
# lets run this simulation 10,000 times
for (i in 1:sim1){

        data_set <- data.frame(rnorm(1000), rnorm(1000))

        colnames(data_set) <- c("y","x")

        # lets fit our linear model
        p_val <- lm(y ~ x, data = data_set) %>%

        # then lets extract our 'tidy' output
                broom::tidy() %>%               

        # then we will remove the intercept coefficient
                dplyr::slice(-1) %>%            

        # then we will select select our p-values
                dplyr::select(p.value) #%>%

        # storing our generated p-value in a list of p-values
        x <- append(x, p_val)

}

# storing our list of p-values as a vector for easier plotting etc.
x <- as.numeric(x)

# looking at our results
hist(x,
     col = "steelblue",
     main = "Distribution of p-values",
     # sub = "test",
     xlab = "P-value",
     ylab = "Count")

Distribution of p-values in our null, simple linear regression case. Note - the distribution of p-values is what we’d expect, seemingly uniform distribution with \(49\) [\(4.9\) %] of our simulated examples reaching sub-\(\alpha\) level significance. The take away here is that the percentage of false positives [our realized type I error rate] is very close to our theoretical estimate of 5%, our A priori Type I error rate - exactly what we would hope.

Now that we understand some of the behavior of p-values, let’s take a look at backward selection using 2 variables.

Backward Selection

x2 <- vector()

double_sig <- 0

# lets simulate this situation 10000 times.
sim2_num <- 1000
for (i in 1:sim2_num){

        # creating testing data
        data_set <- data.frame(rnorm(1000),
                rnorm(1000),
                rnorm(1000)
                )

        colnames(data_set) <- c("y", "x1", "x2")

        #fitting lm
        p_vals <- lm(y ~ x1 + x2, data = data_set) %>%

        # then obtaining tidy output
                broom::tidy() %>%  

        # then removing intercept coef        
                dplyr::slice(-1) %>%            

        # then selecting pvalues
                dplyr::select(term, p.value)

        # checking stopping criteria
        # IF both at least one coefficients is not significant
        # then we remove the coefficient with a higher p-value
        # and re-run the regression

        if (sum(p_vals$p.value > 0.05) >= 1) {
                var_keep <- which.min(p_vals$p.value) + 1
                data_set <- data_set %>%
                        dplyr::select(c(1,var_keep))

                p_vals <- lm(y ~ ., data = data_set) %>%
                        broom::tidy() %>%
                        dplyr::slice(-1) %>%
                        dplyr::select(p.value) %>%
                        as.numeric()

                x2 <- append(x2, p_vals)

        # otherwise, we will record that both coefficients were significant
        } else {
                double_sig <- double_sig + 1
        }
}

hist(x2,
     col = "steelblue",
     main = "Distribution of p-values",
     xlab = "P-values")

Other Options and Takeaways

As we saw here, even a small alteration to our methodologies can have striking effects on our type I error rates, and in turn our results and interpretations. It is important when doing any sort of selection method to consider how this type I error rate [which is really one of the most important concepts of hypothesis testing] may be changing, and report it. Model and covariate selection are tricky concepts, partly for the same reason outlined above. Even some ‘robust’ methodologies such as the LASSO can be tainted with bias [although in a bit different form]. The best defense to these issues can be to understand what biases may be occurring and accounting for them. For instance, it still may be possible to utilize backward selection by using a p-value correction method.