<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/when-the-lasso-fails/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://insightr.wordpress.com/2017/06/14/when-the-lasso-fails/"> R – insightR</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<p><strong>By Gabriel Vasconcelos</strong></p>
</div>
<h2>When the LASSO fails?</h2>
<p>The LASSO has two important uses, the first is forecasting and the second is variable selection. We are going to talk about the second. The variable selection objective is to recover the correct set of variables that generate the data or at least the best approximation given the candidate variables. The LASSO has attracted a lot of attention lately because it allows us to estimate a linear regression with thousands of variables and the model select the right ones for us. However, what many people ignore is when the LASSO fails.</p>
<p><span></span></p>
<p>Like any model, the LASSO also rely on assumptions in order to work. The first is sparsity, i.e. only a small number of variables may actually be relevant. If this assumption does not hold there is no hope to use the LASSO for variable selection. Another assumption is that the irrepresentable condition must hold, this condition may look very technical but it only says that the relevant variable may not be very correlated with the irrelevant variables.</p>
<p>Suppose your candidate variables are represented by the matrix <img alt="X" title="X" src="https://s0.wp.com/latex.php?latex=X&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, where each column is a variable and each line is an observation. We can calculate the covariance matrix <img alt="Sigma=n^{-1} X'X" title="Sigma=n^{-1} X'X" src="https://s0.wp.com/latex.php?latex=%5CSigma%3Dn%5E%7B-1%7D+X%27X&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, which is a symmetric matrix. This matrix may be broken into four pieces:</p>
<p><img alt="displaystyle Sigma=left( begin{array}{cc} C_{1,1} &amp; C_{1,2} \ C_{2,1} &amp; C_{2,2} end{array} right) " title="displaystyle Sigma=left( begin{array}{cc} C_{1,1} &amp; C_{1,2} \ C_{2,1} &amp; C_{2,2} end{array} right) " src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle+%5CSigma%3D%5Cleft%28+%5Cbegin%7Barray%7D%7Bcc%7D+C_%7B1%2C1%7D+%26+C_%7B1%2C2%7D+%5C%5C+C_%7B2%2C1%7D+%26+C_%7B2%2C2%7D+%5Cend%7Barray%7D+%5Cright%29+&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p>The first piece, <img alt="C_{1,1}" title="C_{1,1}" src="https://s0.wp.com/latex.php?latex=C_%7B1%2C1%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, shows the covariances between only the important variables, <img alt="C_{2,2}" title="C_{2,2}" src="https://s0.wp.com/latex.php?latex=C_%7B2%2C2%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is the covariance matrix of the irrelevant variables and <img alt="C_{1,2}" title="C_{1,2}" src="https://s0.wp.com/latex.php?latex=C_%7B1%2C2%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and <img alt="C_{2,1}" title="C_{2,1}" src="https://s0.wp.com/latex.php?latex=C_%7B2%2C1%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> shows the covariances between relevant and irrelevant variables. With that in mind, the irrepresentable condition is:</p>
<p><img alt="displaystyle |C_{2,1}C_{1,1}^{-1}sign(beta)|&lt;1" title="displaystyle |C_{2,1}C_{1,1}^{-1}sign(beta)|&lt;1" src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle+%7CC_%7B2%2C1%7DC_%7B1%2C1%7D%5E%7B-1%7Dsign%28%5Cbeta%29%7C%3C1&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p>The inequality above must hold for all elements. <img alt="sign(beta)" title="sign(beta)" src="https://s0.wp.com/latex.php?latex=sign%28%5Cbeta%29&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is <img alt="1" title="1" src="https://s0.wp.com/latex.php?latex=1&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> for positive values of <img alt="beta" title="beta" src="https://s0.wp.com/latex.php?latex=%5Cbeta&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, <img alt="-1" title="-1" src="https://s0.wp.com/latex.php?latex=-1&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> for negative values and <img alt="0" title="0" src="https://s0.wp.com/latex.php?latex=0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> if <img alt="beta=0" title="beta=0" src="https://s0.wp.com/latex.php?latex=%5Cbeta%3D0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">.</p>
<h2>Example</h2>
<p>For this example we are going to generate two covariance matrices, one that satisfies the irrepresentable condition and one that violates it. Our design will be very simple: only 10 candidate variables where five of them are relevant.</p>
<pre>library(mvtnorm)
library(corrplot)
library(glmnet)
library(clusterGeneration)
</pre>
<pre>k=10 # = Number of Candidate Variables
p=5 # = Number of Relevant Variables
N=500 # = Number of observations
betas=(-1)^(1:p) # = Values for beta
set.seed(12345) # = Seed for replication
sigma1=genPositiveDefMat(k,&quot;unifcorrmat&quot;)$Sigma # = Sigma1 violates the irc
sigma2=sigma1 # = Sigma2 satisfies the irc
sigma2[(p+1):k,1:p]=0
sigma2[1:p,(p+1):k]=0

# = Verify the irrepresentable condition
irc1=sort(abs(sigma1[(p+1):k,1:p]%*%solve(sigma1[1:p,1:p])%*%sign(betas)))
irc2=sort(abs(sigma2[(p+1):k,1:p]%*%solve(sigma2[1:p,1:p])%*%sign(betas)))
c(max(irc1),max(irc2))
</pre>
<pre>## [1] 3.222599 0.000000
</pre>
<pre># = Have a look at the correlation matrices
par(mfrow=c(1,2))
corrplot(cov2cor(sigma1))
corrplot(cov2cor(sigma2))
</pre>
<p><img alt="plot of chunk unnamed-chunk-20" src="https://i0.wp.com/i.imgur.com/Ln567f6.png"></p>
<p>As you can see, <code>irc1</code> violates the irrepresentable condition and <code>irc2</code> does not. The correlation matrix that satisfies the irrepresentable condition is block diagonal and the relevant variables have no correlation with the irrelevant ones. This is an extreme case, you may have a small correlation and still satisfy the condition.</p>
<p>Now let us check how the LASSO works for both covariance matrices. First we need do understand what is the regularization path. The LASSO objective function penalizes the size of the coefficients and this penalization is controlled by a hyper-parameter <img alt="lambda" title="lambda" src="https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. We can find the exact <img alt="lambda_0" title="lambda_0" src="https://s0.wp.com/latex.php?latex=%5Clambda_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> that is sufficiently big to shrink all variables to zero and for any value smaller than <img alt="lambda_0" title="lambda_0" src="https://s0.wp.com/latex.php?latex=%5Clambda_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> some variable will be included. As we decrease the size of <img alt="lambda" title="lambda" src="https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> more variables are included until we have a model with all variables (or the biggest identified model when we have more variables than observations). This path between the model with all variables and the model with no variables is the regularization path. The code below generates data from multivariate normal distributions for the covariance matrix that violates the irrepresentable condition and the covariance matrix that satisfies it. Then I estimate the regularization path for both case and summarize the information in plots.</p>
<pre>X1=rmvnorm(N,sigma = sigma1) # = Variables for the design that violates the IRC = #
X2=rmvnorm(N,sigma = sigma2) # = Variables for the design that satisfies the IRC = #
e=rnorm(N) # = Error = #
y1=X1[,1:p]%*%betas+e # = Generate y for design 1 = #
y2=X2[,1:p]%*%betas+e # = Generate y for design 2 = #

lasso1=glmnet(X1,y1,nlambda = 100) # = Estimation for design 1 = #
lasso2=glmnet(X2,y2,nlambda = 100) # = Estimation for design 2 = #

## == Regularization path == ##
par(mfrow=c(1,2))
l1=log(lasso1$lambda)
matplot(as.matrix(l1),t(coef(lasso1)[-1,]),type=&quot;l&quot;,lty=1,col=c(rep(1,9),2),ylab=&quot;coef&quot;,xlab=&quot;log(lambda)&quot;,main=&quot;Violates IRC&quot;)
l2=log(lasso2$lambda)
matplot(as.matrix(l2),t(coef(lasso2)[-1,]),type=&quot;l&quot;,lty=1,col=c(rep(1,9),2),ylab=&quot;coef&quot;,xlab=&quot;log(lambda)&quot;,main=&quot;Satisfies IRC&quot;)
</pre>
<p><img alt="plot of chunk unnamed-chunk-21" src="https://i0.wp.com/i.imgur.com/dhlujso.png"></p>
<p>The plot on the left shows the results when the irrepresentable condition is violated and the plot on the right is the case when it is satisfied. The five black lines that slowly converge to zero are the five relevant variables and the red line is an irrelevant variable. As you can see, when the IRC is satisfied all relevant variables shrink very fast to zero as we increase lambda. However, when the IRC is violated one irrelevant variable starts with a very small coefficient that slowly increases before decreasing to zero in the very end of the path. This variable is selected through the entire path, it is virtually impossible to recover the correct set of variables in this case unless you apply a different penalty to each variable. This is precisely what the adaptive LASSO does. Does that mean that the adaLASSO is free from the irrepresentable condition??? The answer is: partially. The adaptive LASSO requires a less restrictive condition called weighted irrepresentable condition, which is much easier to satisfy. The two plots below show the regularization path for the LASSO and the adaLASSO in the case the IRC is violated. As you can see, the adaLASSO selects the correct set of variables in all the path.</p>
<pre>lasso1.1=cv.glmnet(X1,y1)
w.=(abs(coef(lasso1.1)[-1])+1/N)^(-1)
adalasso1=glmnet(X1,y1,penalty.factor = w.)

par(mfrow=c(1,2))
l1=log(lasso1$lambda)
matplot(as.matrix(l1),t(coef(lasso1)[-1,]),type=&quot;l&quot;,lty=1,col=c(rep(1,9),2),ylab=&quot;coef&quot;,xlab=&quot;log(lambda)&quot;,main=&quot;LASSO&quot;)
l2=log(adalasso1$lambda)
matplot(as.matrix(l2),t(coef(adalasso1)[-1,]),type=&quot;l&quot;,lty=1,col=c(rep(1,9),2),ylab=&quot;coef&quot;,xlab=&quot;log(lambda)&quot;,main=&quot;adaLASSO&quot;)
</pre>
<p><img alt="plot of chunk unnamed-chunk-22" src="https://i2.wp.com/i.imgur.com/K3Z1Cvm.png"></p>
<p>The biggest problem is that the irrepresentable condition and its less restricted weighted version are not testable in the real world because we need the populational covariance matrix and the true betas that generate the data. The solution is to study your data as much as possible to at least have an idea of the situation.</p>
<h2>Some articles on the topic</h2>
<p>Zhao, Peng, and Bin Yu. “On model selection consistency of Lasso.” Journal of Machine learning research 7.Nov (2006): 2541-2563.</p>
<p>Meinshausen, Nicolai, and Bin Yu. “Lasso-type recovery of sparse representations for high-dimensional data.” The Annals of Statistics (2009): 246-270.</p>
<p> <a rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/insightr.wordpress.com/345/" target="_blank"><img border="0" alt="" src="http://feeds.wordpress.com/1.0/comments/insightr.wordpress.com/345/"></a> <img border="0" alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" data-recalc-dims="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://insightr.wordpress.com/2017/06/14/when-the-lasso-fails/"> R – insightR</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...

	<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/nn8KfWzT4P4" height="1"></a></div>