Regularized regression vs simple regression or Pearson's correlation

Question

Let's say I have a response variable (let's say brain pathology of Parkinson's disease patients) and p input variables (let's say expression of p genes). I want to find a subset of genes which might be driving the pathology of the patients.

One way is to simply check the Pearson's correlation of each gene with the pathology and choose the genes with the highest absolute correlation, or apply a linear regression and select the genes with highest absolute regression coefficients.

A more sophisticated way could be to do regularized regression (let's say Elastic Net or LASSO) to assign nonzero coefficients to only a subset of genes.

I am wondering whether those genes that are assigned nonzero coefficients by Elastic Net or Lasso will be the same as the genes that have the highest Pearson's correlation, or the genes that have the highest coefficients assigned by a simple linear regression? I can code and see whether they are the same or not, but I am also wondering the intuition behind they being the same or different; that's why I am asking here. Thanks!

Answer 1

No, this is false in general.

There is one case where this is true: when your predictors are all completely uncorrelated. In this case variables will enter into the LASSO model in the order of the magnitude of their coefficient values in an unregularized regression model.

When the variables in your model are correlated, this is no longer true. Here's a minimal example in R

set.seed(154)

x_uncorr <- matrix(runif(30000), nrow=10000)
y_uncorr <- 1 + 2*x_uncorr[,1] - x_uncorr[,2] + .5*x_uncorr[,3]

sigma <- matrix(c(  1,   0.9,    0,
                    0.9,   1,    0,
                    0,     0,    1), nrow=3, byrow=TRUE
)
x_corr <- x_uncorr %*% sqrtm(sigma)

The original data I generated, x_uncorr has all predictors uncorrelated. The weird matrix square root / matrix multiplication I did here is to produce a specific correlation structure in the predictor data

round(cor(x_corr), 2)
     [,1] [,2] [,3]
[1,] 1.00 0.90 0.02
[2,] 0.90 1.00 0.01
[3,] 0.02 0.01 1.00

I landed on the correlation structure through experementation, but I knew I would find one with the properties I wanted if I looked hard enough.

Now let's create a response and fit a LASSO model

y_corr <- 1 + 2*x_corr[,1] - x_corr[,2] + 0.5*x_corr[,3]

gnet_corr <- glmnet(x_corr, y_corr)
plot(gnet_corr)

The plot shows the coefficient paths as we vary the regularization parameter

The green path, associated with the third variable, enters the model before the red path, associated with the second parameter. This is true even though the value of the parameters in the unregularized model is in the other direction, which you can see by looking at the final values of the two paths.

Additionally, the correlation between these predictors and the response are

round(cor(cbind(x_corr, y_corr))[1:3, 4], 2)                 
0.86 0.63 0.40 

The second varaible is more correlated with the response than the third, but the third varaible is the first of the two to enter the model.