Why I Stopped Reading Coefficients Of Linear Regression

Linear Regression Coefficients

I am going to assume that if you are reading this then you have at least a basic understanding of linear regression. In general we are trying to find a linear relationship between our input variables \(\underline{x}\) and our out put variable \(y\) of the form

\[ E[y] = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_px_p \]

This model is completely interpretable. We know for sure that when we increase our input variable (\(x_i\)) by a value of 1 then our expected output variable (\(E[y]\)) will increase by the value of the coefficient preceding it (\(\beta_i\)).

So what is my problem? Everything I have said above is completely true, the linear regression model is completely interpretable and I can tell you exactly how it will behave when any of its inputs change. The problem comes when you start trying to apply your knowledge of the model to the real world system. As the saying goes ‘All models are wrong’ and it is this inaccuracy in representation of the system is where our problem lies, but maybe not in the way that you might first think.

The Boston Housing Data

For this example I will be using the Boston Housing dataset. This is available to download online or comes as a part of the mlbench library, which is how I will be accessing it.

# Load the mlbench library containing the data set
library(mlbench)

# Load the Boston Housing Dataset
data("BostonHousing")

# Turn all entries into a number
data <- apply(BostonHousing, 2, as.numeric)
data <- as.data.frame(data)

## # A tibble: 506 × 14
##       crim    zn indus  chas   nox    rm   age   dis   rad   tax ptratio     b
##      <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl>
##  1 0.00632  18    2.31     0 0.538  6.58  65.2  4.09     1   296    15.3  397.
##  2 0.0273    0    7.07     0 0.469  6.42  78.9  4.97     2   242    17.8  397.
##  3 0.0273    0    7.07     0 0.469  7.18  61.1  4.97     2   242    17.8  393.
##  4 0.0324    0    2.18     0 0.458  7.00  45.8  6.06     3   222    18.7  395.
##  5 0.0690    0    2.18     0 0.458  7.15  54.2  6.06     3   222    18.7  397.
##  6 0.0298    0    2.18     0 0.458  6.43  58.7  6.06     3   222    18.7  394.
##  7 0.0883   12.5  7.87     0 0.524  6.01  66.6  5.56     5   311    15.2  396.
##  8 0.145    12.5  7.87     0 0.524  6.17  96.1  5.95     5   311    15.2  397.
##  9 0.211    12.5  7.87     0 0.524  5.63 100    6.08     5   311    15.2  387.
## 10 0.170    12.5  7.87     0 0.524  6.00  85.9  6.59     5   311    15.2  387.
## # … with 496 more rows, and 2 more variables: lstat <dbl>, medv <dbl>

If you would like to kow more about each of the variables in the data set you can use the help function

?BostonHousing

Now we are going to train a linear model to the data

mod <- lm(medv ~ ., data = data)
mod$coefficients

##   (Intercept)          crim            zn         indus          chas 
##  3.645949e+01 -1.080114e-01  4.642046e-02  2.055863e-02  2.686734e+00 
##           nox            rm           age           dis           rad 
## -1.776661e+01  3.809865e+00  6.922246e-04 -1.475567e+00  3.060495e-01 
##           tax       ptratio             b         lstat 
## -1.233459e-02 -9.527472e-01  9.311683e-03 -5.247584e-01

Here we can see rm has a coefficient of about 3.8 which is the largest in size (ignoring the intercept), does this mean that it is the most important predictor in the data set? The obvious answer here is: we can’t say. This is because all of our variables are using different units and this will have a big effect on the coefficients. So before we train a linear model we should be centering and scaling our data

# Centering and scaling the data
data <- scale(data, center = TRUE, scale = TRUE)
data <- as.data.frame(data)

## # A tibble: 506 × 14
##      crim      zn  indus   chas    nox     rm     age   dis    rad    tax
##     <dbl>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <dbl>  <dbl>  <dbl>
##  1 -0.419  0.285  -1.29  -0.272 -0.144  0.413 -0.120  0.140 -0.982 -0.666
##  2 -0.417 -0.487  -0.593 -0.272 -0.740  0.194  0.367  0.557 -0.867 -0.986
##  3 -0.417 -0.487  -0.593 -0.272 -0.740  1.28  -0.266  0.557 -0.867 -0.986
##  4 -0.416 -0.487  -1.31  -0.272 -0.834  1.02  -0.809  1.08  -0.752 -1.11 
##  5 -0.412 -0.487  -1.31  -0.272 -0.834  1.23  -0.511  1.08  -0.752 -1.11 
##  6 -0.417 -0.487  -1.31  -0.272 -0.834  0.207 -0.351  1.08  -0.752 -1.11 
##  7 -0.410  0.0487 -0.476 -0.272 -0.265 -0.388 -0.0702 0.838 -0.522 -0.577
##  8 -0.403  0.0487 -0.476 -0.272 -0.265 -0.160  0.978  1.02  -0.522 -0.577
##  9 -0.396  0.0487 -0.476 -0.272 -0.265 -0.930  1.12   1.09  -0.522 -0.577
## 10 -0.400  0.0487 -0.476 -0.272 -0.265 -0.399  0.615  1.33  -0.522 -0.577
## # … with 496 more rows, and 4 more variables: ptratio <dbl>, b <dbl>,
## #   lstat <dbl>, medv <dbl>

Lets see what is the “most important” variable according to linear regression now?

mod <- lm(medv ~ ., data = data)
mod$coefficients

##   (Intercept)          crim            zn         indus          chas 
## -1.627634e-16 -1.010171e-01  1.177152e-01  1.533520e-02  7.419883e-02 
##           nox            rm           age           dis           rad 
## -2.238480e-01  2.910565e-01  2.118638e-03 -3.378363e-01  2.897491e-01 
##           tax       ptratio             b         lstat 
## -2.260317e-01 -2.242712e-01  9.243223e-02 -4.074469e-01

Here we we can see lstat is now considered to be the most important. Now we can try and bootstrap these coefficients to try and get some idea of a confidence interval on them. To perform a bootstrap, we will resample our data (randomly select rows) with replacement (meaning we can select the same row more than once) such that until we have a new data set that is the same size as our original one. If we calculate a new set of coefficients on this new data set, they should be similar but slightly different. The amount by which the coefficients vary can then be used to infer the uncertainty related with that coefficient. For this experiment I am going to use 25,000 resamples, for this I am going to create a variable N_SAMPLE to save this. This code can take a long time to run and so if you are recreating this you may wish to try with fewer resamples.

# Setting the seed so that results can be recreated
set.seed(123)
# Setting the number of resamples
N_SAMPLE = 25000

I am then going to create a matrix to store the value of each coefficient, I am also going to create a matrix that will store the rank of each variable according to how important it is.

coefs <- matrix(0, ncol = ncol(data), nrow = N_SAMPLE)
ranks <- matrix(0, ncol = ncol(data), nrow = N_SAMPLE)

I then perform the bootstrapping

for (i in 1:N_SAMPLE) {
  data_samples <- data[sample(1:nrow(data), nrow(data), replace = TRUE),]
  mod <- lm(medv ~ ., data = data_samples)
  coefs[i,] <- mod$coefficients
  ranks[i,] <- -abs(mod$coefficients) %>% rank(ties.method = "random")
}

colnames(coefs) <- names(mod$coefficients)
colnames(ranks) <- names(mod$coefficients)

We can then plot the distribution of these coefficients

By doing this, you can see just how wide the range of these coefficients are under perturbation. Not only this, but some variables that are quite clearly o one side of the 0 (e.g. crim) have multiple coefficients that cross 0 changing the way the variable will be interpretted. This is not to be confused with the ceoffecients centered around 0, these too have their own problems though, with many non-zero coefficients seemingly implying relationships that are not present. This might not be a problem in and of itself, but we can see that the ranges overlap quite substantially, so if we look at the rank of each of these coefficients under resampling we can see that varies wildly too.

This shows that which variables are more important varies substantially under small perturbations of the underlying data. This alone was evidence enough for me that I should stop using the coefficients of linear regression as it can’t even consistently rank the importance of variables.

Why is this happening?

There are multiple reasons why this is happening linked to the resampling of the data. But the main reason is correlation in the underlying data.

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

As we can see from this correlation plot, the variables covary a lot. This can make it hard for linear regression to determine which variable is contributing to a change as they are all varying together. This isn’t necessarily a problem for prediction ut is a big problem for inference. It means that if a variable is truly important, it might not have a large coefficient but something it is correlated with does.

This correlation leads to another problem too, as the way in which we are interpretting the system is by definition at odds with the way in which we are interpretting the model. By reading coefficients we are asking how will the prediction change if this variable changes and nothing else does?. But within the system we know that variables never change independently and a lot of the time the idea of one variable changing whilst another stays the same is nonsensical.

Conclusions

This is a very rough overview of some of the many issues reading coefficients of linear models. It is still something I might do as part of exploratory data analysis or a quick sanity check, but it is now something that I am trying to move away from. In any case where I have to do it I try and keep these ideas in mind and I am very careful about any inferential statements that I make.