O LASSO

Este post vai tratar de um método de machine learning muito interessante e relativamente simples: o LASSO. LASSO significa Least Absolute Shrinkage and Select Operator. Como o nome sugere, o LASSO seleciona quais regressores são relevantes e quais não são. Ou seja, suponha que você é um pesquisador que tem 50 variáveis que são possíveis candidatos a variáveis explicativas de uma variável de interesse. O LASSO permite que você dê os 50 regressores para o computador e ele escolha quais são os relevantes. O LASSO está centrado na ideia de esparsidade: poucos coeficientes vão ser diferentes de zero; poucas variáveis vão ser relevantes numa regressão.

Uma pergunta justa é: por que não simplesmente testamos isso “na mão”? Isso é, fazemos uma regressão com todas as variáveis, uma sem a primeira variável da base de dados, uma outra sem a segunda… e no fim vemos qual a melhor regressão? Além de problemas de teste - lembre que todo teste tem uma chance de você estar tomando a decisão errada - isso seria um inferno: com relés 50 variáveis, precisaríamos de $2^{50}$ regressões, ou seja, 1.1258999^{15} regressões: esse é um número com 15 zeros.

A ideia original vem de um paper de 1996 de Robert Tibshirani. Não temos solução fechada para encontrar o estimador de LASSO, ao contrário do MQO, onde podemos expressar o estimador como uma multiplicação de matrizes. Portanto, é um estimador que seria impossível sem um computador. Mas a ideia é incrivelmente simples.

Pegue o modelo usual de regressão, $y_i = X_i \beta + u_i$, onde $\beta$ é o coeficiente de interesse, $u$ é um choque aleatório, as observações são indexadas por $i = 1,...,n$ e os coeficientes são indexados por $k = 1,..,K$ . O objetivo do MQO é minimizar a soma dos quadrados dos resíduos, ou seja $\min_{\beta} \sum_{i=1}^n \hat{u}^2$. O LASSO começa dai, mas adiciona uma pequena alteração: vamos limitar o valor da soma dos valores absolutos dos coeficientes para que ele seja menor que uma constante $c$. Ou seja, o LASSO resolve o problema:

$$\min_{\beta} \sum_{i=1}^n \hat{u}^2 \text{ sujeito a } \sum_{k=1}^K |\beta_{k}| < c$$

Veja que não temos como tirar a derivada da função módulo, e portanto não podemos resolver o problema “no braço”.

Podemos reescrever o problema como um la grangeano: $\min_{\beta} \left( \sum_{i=1}^n \hat{u}^2 \right) - \lambda \left( \sum_{k=1}^K |\beta_{k}| \right)$, onde o $\lambda$ é o multiplicador de lagrange. Existe uma função que leva de $c$ para $\lambda$, então podemos ignorar o valor de $c$ e pensar só em termos de $\lambda$. Muitas implementações fazem isso e eu seguirei este caminho.

Veja que estamos trabalhando a soma do valor absoluto dos coeficientes. Chamamos isso de norma $\ell_1$. Aqueles que já fizeram um curso de algebra linear conhecem a norma $\ell_2$, conhecida como norma euclidiana: $\sum_{k=1}^K \beta{}^2$. Existe um método de estimação que usa a norma $\ell_2$ ao invés da $\ell_1$, conhecido como ridge. Usando a norma $\ell_2$, o problema tem solução analítica. Então, por que norma $\ell_1$, que parece ser tão ruim de trabalhar?

Esta é uma das belezas do LASSO: a norma $\ell_2$ não gera esparsidade. A norma $\ell_1$ gera. O motivo é ilustrado na figura abaixo, tirada de Statistical Learning with Sparsity, de Trevor Hastie, Robert Tibshirani and Martin Wainwright (cujo pdf, 100% legal, você encontra aqui)

A imagem supõe que você só tem duas variáveis, por motivos ilustrativos. O $\hat{\beta}$ é o máximo global, o conhecido estimador de MQO. Em vermelho são as curvas de nível do estimador. O que o Ridge (à direita) faz é impor uma restrição no formato de uma bola (azul). Veja que a curva de nível encosta na bola fora do eixo vertical - onde $\beta_1$ seria zero. Já o LASSO é ilustrado na figura à esquerda. Veja que, nesse caso, a curva de nível acerta na quina da restrição, efetivamente zerando o coeficiente $\beta_1$.

Veja que, por enquanto, eu escondi o $\lambda$ (ou o $c$) debaixo do tapete: eu ainda não expliquei como escolher essa variável, que de fato está no centro do problema. O pesquisador pode ter uma vaga ideia de quanto os coeficientes devem somar. Mas uma ideia exata do tamanho é complicada, então gostariamos de deixar os dados nos dizerem o $\lambda$ certo. Existem várias estratégias, e infelizmente o glmnet só implementa uma.

Com alguma ideia do que o LASSO está fazendo, podemos colocar a mão na massa e discutir o pacote padrão para o LASSO no R, o glmnet.

##O pacote glmnet

Obviamente, começamos carregando o pacote:

library(glmnet)

## Carregando pacotes exigidos: Matrix

## Loaded glmnet 4.0-2

O comando para rodar uma regressão usando o LASSO é o glmnet. Ele recebe basicamente duas coisas, uma matriz com regressores e o vetor que é a variável dependente. Até ai, muito parecido com o MQO. Mas a grande diferença em termos de output é que o glmnet devolve uma matriz de coeficientes. O que o glmnet faz é estimar o modelo para um vetor de $\lambda$, de tal forma que você vá de nenhuma variável incluida no modelo até todas as variáveis. Isso não significa que cada $\lambda$ estimado adicione uma variável nova ao modelo: em alguns casos o LASSO não seleciona mais ninguém, mas como ele tem mais “orçamento” para alocar, ele aumenta o valor dos coeficientes das variáveis incluídas. O glmnet resolve o problema de escolher qual $\lambda$ usar simplesmente estimando com vários valores diferentes.

Vamos testar o glmnet gerando uma matriz com 50 variáveis e 2000 observações e criando uma variável dependente que depende só das 10 primeiras variáveis:

x <- matrix(rnorm(2000*50), ncol = 50) #regressores
betas <- c(rep(1,10),rep(0,40)) # as 10 primeiras serão relevantes com um coeficiente 1. As outras são irrelevantes - e portanto tem um coeficiente 0
y <- x%*%betas + rnorm(2000) 

modelo <- glmnet(x,y)
coef(modelo)[,1:5] # só pegando as 5 primeiras colunas. Veja que o comando summary não nos dá o que queremos no glmnet

## 51 x 5 sparse Matrix of class "dgCMatrix"
##                    s0          s1         s2         s3         s4
## (Intercept) 0.0899452 0.088733817 0.08744122 0.07980104 0.07205909
## V1          .         .           0.07186523 0.15461998 0.22925068
## V2          .         0.006269833 0.09268133 0.17163545 0.24390558
## V3          .         .           .          0.04081292 0.12344005
## V4          .         0.023074549 0.10526749 0.18520641 0.26065968
## V5          .         .           .          0.08775326 0.16789139
## V6          .         .           0.05578103 0.13964770 0.21611625
## V7          .         .           0.01023765 0.09799273 0.17949380
## V8          .         0.095583403 0.17473050 0.24967983 0.32036144
## V9          .         .           0.04487896 0.12992889 0.20907527
## V10         .         .           .          0.06677614 0.14574878
## V11         .         .           .          .          .         
## V12         .         .           .          .          .         
## V13         .         .           .          .          .         
## V14         .         .           .          .          .         
## V15         .         .           .          .          .         
## V16         .         .           .          .          .         
## V17         .         .           .          .          .         
## V18         .         .           .          .          .         
## V19         .         .           .          .          .         
## V20         .         .           .          .          .         
## V21         .         .           .          .          .         
## V22         .         .           .          .          .         
## V23         .         .           .          .          .         
## V24         .         .           .          .          .         
## V25         .         .           .          .          .         
## V26         .         .           .          .          .         
## V27         .         .           .          .          .         
## V28         .         .           .          .          .         
## V29         .         .           .          .          .         
## V30         .         .           .          .          .         
## V31         .         .           .          .          .         
## V32         .         .           .          .          .         
## V33         .         .           .          .          .         
## V34         .         .           .          .          .         
## V35         .         .           .          .          .         
## V36         .         .           .          .          .         
## V37         .         .           .          .          .         
## V38         .         .           .          .          .         
## V39         .         .           .          .          .         
## V40         .         .           .          .          .         
## V41         .         .           .          .          .         
## V42         .         .           .          .          .         
## V43         .         .           .          .          .         
## V44         .         .           .          .          .         
## V45         .         .           .          .          .         
## V46         .         .           .          .          .         
## V47         .         .           .          .          .         
## V48         .         .           .          .          .         
## V49         .         .           .          .          .         
## V50         .         .           .          .          .

Obviamente, ainda temos que selecionar qual dos modelos nós queremos. O pacote nos permite selecionar por Cross Validation (CV). Em linhas gerais, o que o Cross Validation faz é separar os dados em $k$ blocos - 5,10, ou até mesmo blocos de uma única observação (conhecido como leave one out). Depois, escolha $k-1$ blocos e estime o modelo usando o glmnet. Para avaliar qual $\lambda$ devemos usar, use o bloco que não foi usado para estimar para ver qual modelo gera o menor erro segundo alguma métrica (erro quadrático médio, por exemplo). O CV vai permitir com que cada bloco tenha uma chance de ser “fora da amostra”.

Isso pode parecer muito computacionalmente intensivo - e de fato é - mas a implementação do glmnet é tão eficiente que normalmente num piscar de olhos o modelo é estimado - no meu Dell Vostro de 2012 com i5, 6GB e Windows 10, o tempo é de apenas 0.19s. O comando que faz a seleção do $\lambda$ é o cv.glmnet. Vamos testar no exemplo acima:

cross_val <- cv.glmnet(x,y)

coef(cross_val)

## 51 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept) -0.001739634
## V1           0.937859461
## V2           0.930282416
## V3           0.908357779
## V4           0.977442350
## V5           0.929193703
## V6           0.942473364
## V7           0.953166419
## V8           0.991665138
## V9           0.960682012
## V10          0.896159947
## V11          .          
## V12          .          
## V13          .          
## V14          .          
## V15          .          
## V16          .          
## V17          .          
## V18          .          
## V19          .          
## V20         -0.010307812
## V21          .          
## V22          .          
## V23          .          
## V24          .          
## V25          .          
## V26          .          
## V27          .          
## V28          .          
## V29          .          
## V30          .          
## V31          .          
## V32          .          
## V33          .          
## V34          .          
## V35          .          
## V36          .          
## V37          .          
## V38          .          
## V39          .          
## V40          .          
## V41          .          
## V42          .          
## V43          .          
## V44          .          
## V45          .          
## V46          .          
## V47          .          
## V48          .          
## V49          .          
## V50          .

O LASSO com CV põe apenas um coeficiente a mais. O cv.glmnet também fornece um plot com o número de coeficientes e o erro médio do Cross Validation, o que ajuda a ilustrar o ponto:

plot(cross_val)

Uma das linhas tracejadas marca o número de coeficientes escolhidos: 11, número em cima. Depois dessa linha tracejada, o Mean Square Error não flutua muito. E também não há muita variação para imediatamente a direita, onde temos o modelo verdadeiro com 10 coeficientes.

(O autor agradece aos Professores Marcelo Medeiros e Pedro Souza pelas longas explicações sobre o LASSO e métodos correlatos. Pedro Cava, o outro autor deste blog, foi primordial em incentivar esse texto. Todos os erros neste texto são, como de praxe, de minha exclusiva culpa)