LASSO Adaptativo e Critérios de Informação para LASSO

Daniel Coutinho

Em um post anterior, eu falei do LASSO (Least Absolute Shrinkage and Select Operator). Como vamos explorar uma variação do LASSO hoje, eu vou repetir o problema que o LASSO resolvia:

$\hat{\beta}_{LASSO} \in \arg \min_{\beta} \displaystyle \sum_{i=1}^n (y_i - x_i \beta)^2 + \lambda \sum_{k=0}^p |\beta_k|$

(Onde $|.|$ é o valor absoluto do termo). E como eu já disse, o LASSO nos fornece uma maneira de selecionar quais variáveis entram no modelo ou não. Vamos fazer um pequeno teste com o LASSO: eu vou gerar 50 variáveis normais, independentes, e dessas dez - as dez primeiras - eu colocarei $\beta=1$ . As outra vão ser irrelevantes para o problema e vão ter $\beta=0$ . O tamanho da amostra vai ser igual a 100.

Como de praxe, nós podemos ter diversos objetivos ao estimar um modelo. Eu vou comparar 3 coisas: a quantidade de vezes que o LASSO coloca as variáveis relevantes, a quantidade de vezes que ele exclui as variáveis irrelevantes e quando ele obtém o modelo certo - o que exige colocar todo mundo que é relevante e excluir todos os irrelevantes. Vou fazer só 500 replicações e usar o Cross Validation (que eu já discuti aqui) para escolher o $\lambda$ :

library(glmnet)

coeficientes <- Matrix(0,ncol=500,nrow=51)

for(i in 1:500){
  x <- matrix(rnorm(50*100),ncol = 50) #gerando x
  betas <- c(rep(1,10),rep(0,40)) #
  y <- x%*%betas + rnorm(100)
  modelo <- cv.glmnet(x,y) #estimando usando LASSO e Cross Validation
  coeficientes[,i] <- coef(modelo)
}

#Fim da simulação

analise <- matrix(0,ncol=3,nrow=500)
colnames(analise)<- c("Não zeros certos","Zeros certos", "Modelo certo?")

for(i in 1:500){
  analise[i,1]<- mean(coeficientes[2:11,i] != 0)
  analise[i,2]<- mean(coeficientes[12:51,i] == 0)
  analise[i,3]<- ifelse(analise[i,1]+analise[i,2] == 2,1,0)
}

tabela_final <- colMeans(analise)*100
knitr::kable(tabela_final,caption = "Os valores estão em porcentagem")

Table 1: Os valores estão em porcentagem
	x
Não zeros certos	100.00
Zeros certos	83.31
Modelo certo?	0.20

O LASSO sempre inclui as variáveis relevantes, e exclui as variáveis irrelevantes em 84% das vezes. Mesmo assim, a proporção de vezes que o LASSO consegue recuperar o modelo correto é por volta de 1%. Isso parece esquisito a primeira vista, mas lembre que recuperar o modelo certo envolve acertar todas as relevantes e todas as irrelevantes. Se tivermos 50 variáveis e 500 replicações e em toda replicação o LASSO colocar apenas umas variável irrelevante no modelo, nós teríamos 98% de zeros certos ( $49/50$ ) e exatamente 0 modelos certos.

Parte do problema é que o LASSO penaliza todos os coeficientes igualmente, usando o $\lambda$ . Nós esperaríamos que algumas variáveis sejam mais importantes que outras - e isso pode vir a priori ou ser dito pelos dados. O LASSO adaptativo (adaLASSO) adiciona pesos ( $\omega$ ) para cada uma das variáveis na penalidade. Logo o novo problema a ser resolvido é:

$\hat{\beta}_{adaLASSO} \in \arg \min_{\beta} \displaystyle \sum_{i=1}^n (y_i - x_i \beta)^2 + \lambda \sum_{k=0}^p \frac{|\beta_k|}{\omega_k}$

A única exigência desses pesos é que eles sejam positivos. Veja que nós temos um novo parâmetro a escolher, os pesos. Eis um algoritmo muito simples que gera os pesos baseado nos dados e usa o LASSO:

Estime o modelo usando LASSO. Guarde os coeficientes, que eu chamarei de $\beta_{LASSO}$
Defina $\omega_k = \frac{1}{|\beta_{LASSO}|}$
Estime o adaLASSO usando os pesos definidos em 2.

A boa notícia é que o glmnet já nos oferece uma opção para colocar o peso, via o argumento penalty.factor. Nosso trabalho é basicamente reduzido a escrever umas duas linhas de código a mais: uma que faz o LASSO de “primeiro estágio” e outra que define os pesos.

Uma coisa deve ficar evidente: da maneira que os pesos foram estabelecidos no meu algoritmo, coeficientes que são zerados pelo LASSO serão automaticamente excluídos pelo LASSO adaptativo: teremos como peso $1/0$ (com perdão aos matemáticos), que no limite é infinito. Como o LASSO não excluiu as variáveis relevantes em nenhum caso, não vamos nos preocupar. Mas é possível imaginar situações em que o LASSO poderia ter problemas (pense em coeficientes altissimamente correlacionados em uma amostra relativamente pequena). Outra coisa que deve ficar clara é que precisamos selecionar o $\lambda$ agora duas vezes, uma para o LASSO e outra para a etapa “adaptativa” do LASSO!

Vamos repetir a simulação ali de cima, mas usando o adaLASSO. Em ambos os estágios eu vou usar o Cross Validation:

coeficientes_adalasso <- Matrix(0,ncol=500,nrow=51)

for(i in 1:500){
  x <- matrix(rnorm(50*100),ncol = 50) #gerando x
  betas <- c(rep(1,10),rep(0,40)) #
  y <- x%*%betas + rnorm(100)
  primeiro_estagio <- cv.glmnet(x,y) #estimando usando LASSO e Cross Validation
  pesos <- 1/abs(coef(primeiro_estagio)[-1,]) #veja que eu tenho que jogar fora o intercepto
  adalasso <- cv.glmnet(x,y,penalty.factor = pesos)
  coeficientes_adalasso[,i] <- coef(adalasso)
}

#Fim da simulação

analise <- matrix(0,ncol=3,nrow=500)
colnames(analise)<- c("Não zeros certos","Zeros certos", "Modelo certo?")

for(i in 1:500){
  analise[i,1]<- mean(coeficientes_adalasso[2:11,i] != 0)
  analise[i,2]<- mean(coeficientes_adalasso[12:51,i] == 0)
  analise[i,3]<- ifelse(analise[i,1]+analise[i,2] == 2,1,0)
}

tabela_final <- colMeans(analise)*100
knitr::kable(tabela_final,caption = "Os valores estão em porcentagem")

Table 2: Os valores estão em porcentagem
	x
Não zeros certos	100.00
Zeros certos	96.84
Modelo certo?	58.20

A performance do LASSO adaptativo é muito melhor que a do LASSO: em aproximadamente 60% dos casos agora recuperamos o modelo correto, contra um pouco mais de 1% dos casos para o LASSO. O adaLASSO é uma modificação extremamente simples do algoritmo do LASSO que gera excelente resultados. Veja que a escolha do Cross Validation aqui é apenas pela conveniência do glmnet já trazer o Cross Validation. Veja que implementar um critério de informação é extremamente simples. Eu vou descrever o algoritmo para o LASSO, e para o adaLASSO basta juntar os dois procedimentos:

Estime o modelo usando LASSO via o glmnet (não o cv.glmnet!). Isso vai devolver uma matriz de coeficientes, de dimensão $p \times L$ , onde $L$ é a quantidade de lambdas que a função usou.
Calcule o resíduo para cada modelo estimado $u_l = y - X\beta_l$
Calcule a soma do quadrado dos erros para cada modelo $SSR = \sum_{i=1}^{n} u_i^2$ .
Crie um vetor que tem a quantidade de coeficientes não zeros para cada modelo estimado. Vamos nos referir a cada entrada desse modelo como $s$ .
Calcule $IC = n \ln(SSR) + cs$ para cada modelo, onde c é o critério de informação escolhido
Escolha o lambda que minimiza o valor de IC

Vamos escrever uma função que faz cada um dos passos. Veja que você não precisa fazer isso: existe um pacote excelente que já implementa isso no GitHub chamado HDeconometrics, mantido pelo Gabriel Vasconcelos. Mas seguindo o espírito do blog, eu vou implementar “no braço”. Eu vou permitir que o usuário escolha qualquer um dos 3 critérios, e para isso eu usarei uns if:

ic_glmnet <- function(x,y,penalty.factor,ic){
  require(glmnet)
  n <- length(y)
  modelo <- glmnet(x,y) #passo 1
  #passo 2
  x_aux <- cbind(1,x)
  u <- y - x_aux%*%coef(modelo)
  #passo 3
  ssr <- colSums(u^2)
  #passo 4
  conj_ativo <- colSums(coef(modelo) !=0)
  #vamos permitir o usuário escolha qual critério de informação vai ser usado
  if(ic == "aic"|ic=="AIC"){
    ic_val <- 2} else if(ic=="bic"|ic=="BIC"){
      ic_val <- log(n)
    } else if(ic == "hqc"|ic == "HQC"){
      ic_val <- 2* log(log(n))
    } else{
      stop("IC not implemented")
    }
  #passo 5
  val <- n*log(ssr)+ic_val*conj_ativo
  #passo 6
  selecionado <- which.min(val)
  return(list("modelo_selecionado" = coef(modelo)[,selecionado],"modelo_glmnet" = modelo, "lambda_selecionado" = modelo$lambda[selecionado]))
}

Veja que eu faço a função retorna os coeficientes escolhidos, o modelo inteiro escolhido e o $\lambda$ escolhido ( o motivo para isso vai ficar claro). Vamos fazer um pequeno teste da nossa função:

## Carregando pacotes exigidos: glmnet

## Carregando pacotes exigidos: Matrix

## Loaded glmnet 4.0-2

## (Intercept)          V1          V2          V3          V4          V5 
## -0.07163931  0.91825550  0.82485362  0.80345212  0.81528969  0.68988950 
##          V6          V7          V8          V9         V10         V11 
##  0.89619273  0.82549050  0.81241610  0.90796060  0.77718892  0.00000000 
##         V12         V13         V14         V15         V16         V17 
##  0.00000000  0.00000000 -0.11836918  0.00000000  0.00000000  0.00000000 
##         V18         V19         V20         V21         V22         V23 
## -0.01061692  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 
##         V24         V25         V26         V27         V28         V29 
##  0.00000000  0.00000000  0.02245120  0.00000000  0.00000000  0.00000000 
##         V30         V31         V32         V33         V34         V35 
##  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 
##         V36         V37         V38         V39         V40         V41 
##  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 
##         V42         V43         V44         V45         V46         V47 
##  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000  0.08652614 
##         V48         V49         V50 
##  0.00000000  0.00000000  0.00000000

E uma simulação, com uma mudança: nós não vamos escolher o $\lambda$ duas vezes. Nós vamos repetir o $\lambda$ escolhido no primeiro estágio para o segundo estágio. Isso tem dois objetivos:

Isola o efeito da seleção do lambda dos efeitos da pessagem do LASSO adaptativo
Ao invés de termos de selecionar os paramêtros duas vezes, selecionamos apenas uma vez.

Vamos a simulação:

coeficientes_adalasso_bic <- Matrix(0,ncol=500,nrow=51)

for(i in 1:500){
  x <- matrix(rnorm(50*100),ncol = 50) #gerando x
  betas <- c(rep(1,10),rep(0,40)) #
  y <- x%*%betas + rnorm(100)
  primeiro_estagio <- ic_glmnet(x,y,ic = "BIC") #estimando usando LASSO e BIC
  pesos <- 1/abs(primeiro_estagio$modelo_selecionado[-1]) #veja que eu tenho que jogar fora o intercepto
  adalasso <- glmnet(x,y,lambda = primeiro_estagio$lambda_selecionado,penalty.factor = pesos)
  coeficientes_adalasso_bic[,i] <- coef(adalasso)
}

#Fim da simulação

analise <- matrix(0,ncol=3,nrow=500)
colnames(analise)<- c("Não zeros certos","Zeros certos", "Modelo certo?")

for(i in 1:500){
  analise[i,1]<- mean(coeficientes_adalasso_bic[2:11,i] != 0)
  analise[i,2]<- mean(coeficientes_adalasso_bic[12:51,i] == 0)
  analise[i,3]<- ifelse(analise[i,1]+analise[i,2] == 2,1,0)
}

tabela_final <- colMeans(analise)*100
knitr::kable(tabela_final,caption = "Os valores estão em porcentagem")

Table 3: Os valores estão em porcentagem
	x
Não zeros certos	99.980
Zeros certos	95.305
Modelo certo?	44.000

Veja que o BIC é bastante bem sucedido em recuperar o modelo certo. Veja que como mantivemos a penalidade parada, todo o ganho vem de acertar os pesos melhor, ilustrando bem qual é o segredo do adaLASSO. Apesar disso, ele é pior do que o Cross Validation. Isso vem com duas observações, que já foram discutidas no post de Cross Validation mas eu repito aqui:

CV é um processo computacionalmente intensivo: quebre em blocos os seus dados, estime o modelo repetidas vezes.
Mais importante, Cross Validation supõe independência entre as observações, ou seja, não podemos usar ele para dados de séries temporais

O BIC não sofre de nenhum desses problemas, o que faz dele uma alternativa atraente para selecionar modelos com LASSO adaptativo.

Esse é um post que aborda duas coisas: LASSO adaptativo e seleção do parâmetro de penalidade via critérios de informação. Ambos tem aplicações interessantes e não precisam ser empregadas em conjunto.

Mostramos como o adaLASSO é bastante superior ao LASSO em seleção de modelos. Uma pergunta justa é se a diferença de previsão dos modelos são tão diferentes assim…

LASSO

LASSO Adaptativo e Critérios de Informação para LASSO

Daniel Coutinho

O LASSO