Capítulo IX. Seleção de variáveis


A seleção de variáveis é muito importante na modelagem estatística. Frequentemente, não estamos apenas interessados em usar um modelo para previsão, mas também precisamos identificar corretamente as variáveis relevantes, ou seja, recuperar o modelo correto sob determinadas premissas. Sabe-se que sob certas condições, o método dos mínimos quadrados ordinários (OLS) produz resultados de previsão ruins e não produz um modelo parcimonioso causando overfitting. Portanto, o objetivo dos métodos de seleção de variáveis é encontrar as variáveis que são mais relevantes para a predição. Tais métodos são particularmente importantes quando o verdadeiro modelo subjacente tem uma representação esparsa, ou seja, muitos parâmetros próximos de zero. A identificação de variáveis relevantes reduzirá o ruído e, portanto, melhorará o desempenho de previsão do modelo ajustado.

Alguns métodos de regularização populares usados são a regressão ridge, seleção de subconjuntos, penalização da norma \(L_1\) e suas modificações e combinações. A regressão ridge, por exemplo, que minimiza uma soma de quadrados residual penalizada usando a penalidade de norma \(L_2\) ao quadrado, é empregada para melhorar a estimativa OLS por meio de uma compensação de variância de viés. No entanto, a regressão Ridge tem a desvantagem de não produzir um modelo parcimonioso, pois mantém todos os preditores no modelo e, portanto, cria um problema de interpretabilidade. Ele também fornece erros de previsão próximos aos do modelo OLS.

Outra abordagem proposta para a seleção de variáveis é o chamado operador de seleção e contração mínima absoluta (Lasso), que visa combinar os recursos da regressão Ridge e a seleção de sub-conjuntos, retendo e encolhendo os coeficientes ou definindo-os como zero.

Referência: Least Angle Regression
Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani. Statistics Department, Stanford University. January 9, 2003.


IX.1. Regressão ridge


Motivação: muitos preditores

Não é incomum ver o número de variáveis de entrada exceder bastante o número de observações, por exemplo, estudos de poluição ambiental, etc.. Com muitos preditores, o ajuste do modelo completo sem penalização resultará em grandes intervalos de previsão e o estimador de regressão de mínimos quadrados pode não existir exclusivamente.

Motivação: \(X\) mal condicionado

Como as estimativas de mínimos quadrados dependem de \((X^\top X)^{-1}\), teríamos problemas na computação dos \(\beta\) se \(X^\top X\) fosse singular ou quase singular. Nesses casos, pequenas alterações nos elementos de \(X\) levam a grandes alterações em \((X^\top X)^{-1}\). O estimador de mínimos quadrado de \(\beta\) pode fornecer um bom ajuste aos dados de treinamento, mas não se encaixará suficientemente bem aos dados de teste.

Uma saída dessa situação é abandonar a exigência de um estimador não viciado ou imparcial. Assumimos apenas que \(X\) e \(Y\) foram centrados, para que não precisemos de um termo constante na regressão:
(a) \(X\) é uma matriz \(n\times p\) com colunas centradas,
(b) \(Y\) é um vetor \(n\times 1\) centrado.

Hoerl and Kennard (1970) propuseram que a instabilidade potencial no estimador de mínimos quadrados \[ \widehat{\beta}=(X^\top X)^{-1}Xy \]

poderia ser melhorado adicionando um pequeno valor constante \(\lambda\) às entradas diagonais da matriz \(X^\top X\) antes de tomar seu inverso. O resultado é o estimador de regressão ridge \[ \widehat{\beta}_\mbox{ridge}=(X^\top X +\lambda \mbox{I}_p)^{-1}Xy\cdot \]

A regressão ridge coloca uma forma específica de restrição nos parâmetros \(\beta\): \(\widehat{\beta}_\mbox{ridge}\) é escolhido para minimizar a soma penalizada dos quadrados: \[ \sum_{i=1}^n \big(y_i-\sum_{j=1}^p x_{ij}\beta_j\big)^2-\lambda \sum_{j=1}^p \beta_j^2, \] que é equivalente à minimização de \[ \sum_{i=1}^n \big(y_i-\sum_{j=1}^p x_{ij}\beta_j\big)^2, \qquad \mbox{restrito à} \qquad \sum_{j=1}^p \beta_j^2<c, \]

para algum \(c> 0\), isto é, restringindo a soma dos coeficientes ao quadrados.

Portanto, a regressão Ridge coloca outras restrições nos parâmetros \(\beta\), no modelo linear. Nesse caso, o que estamos fazendo é que, em vez de apenas minimizar a soma residual de quadrados, também temos um termo de penalidade nos \(\beta\). Esse termo de penalidade é \(\lambda\), uma constante pré-escolhida, vezes a norma quadrada do vetor \(\beta\). Isso significa que, se os \(\beta\) assumirem grandes valores, a função de otimização será penalizada. Preferimos tomar \(\beta\) menores ou \(\beta\) que estão próximos de zero para conduzir o termo de penalidade pequeno.

Propriedades do estimador ridge

\(\widehat{\beta}\) é um estimador não viciado de \(\beta\); \(\widehat{\beta}_\mbox{ridge}\) é um estimador tendencioso de \(\beta\). Para covariáveis ortogonais, \(X^\top X = n\mbox{I}_p\), \(\widehat{\beta}_{ridge}=\frac{n}{n+\lambda}\widehat{\beta}\). Portanto, neste caso, o estimador ridge sempre produz retração em direção a 0. \(\lambda\) controla a quantidade de encolhimento.

Um conceito importante em encolhimento é o grau de liberdade eficaz associado a um conjunto de parâmetros. Em um cenário de regressão ridge:
Se escolhermos \(\lambda = 0\), temos \(p\) parâmetros, pois não há penalização.
Se \(\lambda\) for grande, os parâmetros são fortemente restritos e os graus de liberdade serão efetivamente menores, tendendo a 0 quando \(\lambda\to\infty\).

Os graus de liberdade efetivos associados a \(\beta_1,\beta_2,\cdots,\beta_p\) é definido como \[ \mbox{df}(\lambda)=\mbox{tr}\Big( X\big(X^\top X+\lambda \mbox{I}_p\big)^{-1}X^\top \Big) = \sum_{j=1}^p \dfrac{d_j^2}{d_j^2+\lambda}, \]

onde \(d_j\) são os valores singulares de \(X\). Observe que \(\lambda=0\), que corresponde a nenhum encolhimento, dá \(\mbox{df}(\lambda)=p\), desde que \(X^\top X\) não seja singular, como seria de esperar.

Existe um mapeamento um a um entre \(\lambda\) e os graus de liberdade, então, na prática, pode-se simplesmente escolher os graus de liberdade efetivos que se gostaria de associar ao ajuste e resolver para \(\lambda\).

Como alternativa a um \(\lambda\) escolhido pelo usuário, a validação cruzada é frequentemente usada na escolha de \(\lambda\): selecionamos \(\lambda\) que produz o menor erro de previsão de validação cruzada. O intercepto \(\beta_0\) foi deixado de fora do termo de penalidade porque \(Y\) foi centralizado. A penalização do intercepto faria com que o procedimento dependesse da origem escolhida para \(Y\).

Como o estimador ridge é linear, é simples calcular a matriz de variância-covariância \[ \mbox{Var}(\widehat{\beta}_{ridge})=\sigma^2 \big(X^\top X+\lambda \mbox{I}_p\big)^{-1} X^\top X\big(X^\top X+\lambda \mbox{I}_p\big)^{-1}\cdot \] 


library(MASS)
data(Boston)
str(Boston)
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

Veremos agora como ajustar uma regressão ridge usando o pacote de funções Caret. Usaremos esta biblioteca, pois ela nos fornece muitos recursos para modelagem da vida real.

Para isso, usamos o método train. Passamos os parâmetros, ou seja, escolhemos como resposta medv e as restantes variáveis no arquivo de dados como explicativas, além disso, indicamos method = ‘ridge’ para dizer indicar qual modelo ajustar. 

library(caret)
set.seed(1)
modelo <- train(
  medv ~ .,
  data = Boston,
  method = 'ridge'
)
modelo
## Ridge Regression 
## 
## 506 samples
##  13 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 506, 506, 506, 506, 506, 506, ... 
## Resampling results across tuning parameters:
## 
##   lambda  RMSE      Rsquared   MAE     
##   0e+00   4.914881  0.7152900  3.484317
##   1e-04   4.914758  0.7153023  3.484052
##   1e-01   4.917725  0.7162571  3.412941
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 1e-04.

Podemos mostrar o modelo para ver como a regressão ridge reduz diferentes coeficientes nos preditores. 

plot(modelo)

Também podemos criar um gráfico de importância de cada variável para ver uma classificação de quais preditores têm mais impacto. 

plot(varImp(modelo))

Pré-processamento com Caret: um recurso que usamos do Caret é o pré-processamento. Muitas vezes, na ciência de dados da vida real, queremos executar algum pré-processamento antes da modelagem. Vamos centralizar e dimensionar nossos dados passando o seguinte para o método train: 

set.seed(1)
modelo2 <- train(
  medv ~ .,
  data = Boston,
  method = 'ridge',
  preProcess = c("center", "scale")
)
modelo2
## Ridge Regression 
## 
## 506 samples
##  13 predictor
## 
## Pre-processing: centered (13), scaled (13) 
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 506, 506, 506, 506, 506, 506, ... 
## Resampling results across tuning parameters:
## 
##   lambda  RMSE      Rsquared   MAE     
##   0e+00   4.914881  0.7152900  3.484317
##   1e-04   4.914758  0.7153023  3.484052
##   1e-01   4.917725  0.7162571  3.412941
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 1e-04.

Dividindo o conjunto de dados: muitas vezes, quando estamos modelando, queremos dividir nossos dados em um conjunto de treinamento e teste. Desta forma, podemos verificar se há overfitting. Podemos usar o método createDataPartition para fazer isso. Neste exemplo, usamos o medv de destino para dividir em 80/20, p = 0.80.

Essa função retornará índices que contêm 80% dos dados que devemos usar para treinamento. Em seguida, usamos os índices para obter nossos dados de treinamento do conjunto de dados. 

set.seed(1)
inTraining <- createDataPartition(Boston$medv, p = .80, list = FALSE)
training <- Boston[inTraining,]
testing  <- Boston[-inTraining,]

Podemos então ajustar nosso modelo novamente usando apenas os dados de treinamento. 

set.seed(1)
modelo3 <- train(
  medv ~ .,
  data = training,
  method = 'ridge',
  preProcess = c("center", "scale")
)
modelo3
## Ridge Regression 
## 
## 407 samples
##  13 predictor
## 
## Pre-processing: centered (13), scaled (13) 
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 407, 407, 407, 407, 407, 407, ... 
## Resampling results across tuning parameters:
## 
##   lambda  RMSE      Rsquared   MAE     
##   0e+00   4.895637  0.7167782  3.454826
##   1e-04   4.895610  0.7167819  3.454690
##   1e-01   4.949115  0.7139877  3.440440
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 1e-04.

Agora, queremos verificar nossos dados no conjunto de teste. Podemos usar o método subset para obter os recursos e o destino de teste. Em seguida, usamos o método de previsão passando em nosso modelo de cima e os recursos de teste.

test.features = subset(testing, select=-c(medv))
test.target = subset(testing, select=medv)[,1]
predictions = predict(modelo3, newdata = test.features)


Por fim, calculamos o a raíz do erro quadrático médio (RMSE) e o \(R^2\) para comparar com o modelo acima. 

# RMSE
sqrt(mean((test.target - predictions)^2))
## [1] 5.120435
# R2
cor(test.target, predictions)^2
## [1] 0.7126274

Validação cruzada: na prática, normalmente não construímos nossos dados no conjunto de treinamento. É comum usar uma estratégia de particionamento de dados como validação cruzada que reamostra e divide nossos dados muitas vezes. Em seguida, treinamos o modelo nessas amostras e escolhemos o melhor modelo. Caret facilita isso com o método trainControl.

Usaremos validação cruzada de 10 vezes neste. Para fazer isso, precisamos passar três parâmetros method = “cv”, number = 10 (para 10 vezes). Armazenamos esse resultado em uma variável. 

set.seed(1)
ctrl <- trainControl(
  method = "cv",
  number = 10,
)

Agora, podemos treinar novamente nosso modelo e passar a resposta em trainControl para o parâmetro trControl. Observe que nossa chamada adicionou trControl = set.seed. 

set.seed(1)
modelo4 <- train(
  medv ~ .,
  data = training,
  method = 'ridge',
  preProcess = c("center", "scale"),
  trControl = ctrl
)
modelo4
## Ridge Regression 
## 
## 407 samples
##  13 predictor
## 
## Pre-processing: centered (13), scaled (13) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 367, 366, 367, 366, 365, 367, ... 
## Resampling results across tuning parameters:
## 
##   lambda  RMSE      Rsquared   MAE     
##   0e+00   4.738682  0.7399459  3.427828
##   1e-04   4.738608  0.7399540  3.427727
##   1e-01   4.765288  0.7393663  3.414240
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 1e-04.

Esses resultados parecem ter melhorado nossa precisão para nossos dados de treinamento. Vamos verificar isso nos dados de teste para ver os resultados. 

test.features = subset(testing, select=-c(medv))
test.target = subset(testing, select=medv)[,1]
predictions = predict(modelo4, newdata = test.features)
# RMSE
sqrt(mean((test.target - predictions)^2))
## [1] 5.120435
# R2
cor(test.target, predictions)^2
## [1] 0.7126274

Ajustando os Hiperparâmetros: para ajustar um modelo ridge, podemos dar ao modelo diferentes valores de \(\lambda\). Caret irá treinar novamente o modelo usando diferentes lambdas e selecionar a melhor versão. 

set.seed(1)
tuneGrid <- expand.grid(
  .lambda = seq(0, .1, by = 0.01)
)
modelo5 <- train(
  medv ~ .,
  data = training,
  method = 'ridge',
  preProcess = c("center", "scale"),
  trControl = ctrl,
  tuneGrid = tuneGrid
)
modelo5
## Ridge Regression 
## 
## 407 samples
##  13 predictor
## 
## Pre-processing: centered (13), scaled (13) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 367, 366, 367, 366, 365, 367, ... 
## Resampling results across tuning parameters:
## 
##   lambda  RMSE      Rsquared   MAE     
##   0.00    4.738682  0.7399459  3.427828
##   0.01    4.733617  0.7405339  3.419158
##   0.02    4.731912  0.7408108  3.413409
##   0.03    4.732410  0.7408942  3.409016
##   0.04    4.734467  0.7408497  3.405893
##   0.05    4.737699  0.7407168  3.403919
##   0.06    4.741863  0.7405205  3.403654
##   0.07    4.746800  0.7402774  3.405022
##   0.08    4.752398  0.7399990  3.407386
##   0.09    4.758580  0.7396933  3.410556
##   0.10    4.765288  0.7393663  3.414240
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 0.02.

Finalmente, podemos mostrar novamente o modelo para ver como ele se comporta em diferentes parâmetros de ajuste. 

plot(modelo5)


IX.2. Lasso


Tibshirani (1996) introduziu o Lasso para modelos lineares generalizados, onde a variável de resposta \(y\) é contínua e não categórica. Lasso tem duas características importantes. Primeiro, possui um termo de \(L_1\)-penalidade que realiza retração em coeficientes de maneira semelhante à regressão ridge, onde é usada uma \(L_2\)-penalidade.

Segundo, diferentemente da regressão ridge, o Lasso realiza a seleção de subconjuntos variável, impulsionando alguns coeficientes exatamente zero devido à natureza da restrição, onde a função objetivo pode tocar a área de restrição quadrática em um canto. Por esse motivo, o Lasso é capaz de produzir soluções esparsas e, portanto, é capaz de combinar boas características do procedimento de regressão e seleção de subconjuntos de regressão ridge. Ele produz modelos interpretáveis e tem a estabilidade da regressão ridge.


IX.2.1. Lasso no modelo de regressão linear


O modelo de regressão linear pode ser escrito como \[ y=X\beta +\epsilon, \] onde \(y\) é um vetor \(n\times 1\) das observações da variável resposta, \(X=(x_1^\top,\cdots,x_n^\top)^\top\), \(x_i\in\mathbb{R}^p\), \(i=1,\cdots,n\) a matriz de planejamento contendo as informações de \(p\) variáveis explicativas da resposta e \(\epsilon=(\epsilon_1,\cdots,\epsilon_n)^\top\) o vetor de ero aleatório com \(\mbox{E}(\epsilon_i)=0\) e \(\mbox{Var}(\epsilon_i)=\sigma^2\).

Nesta estrutura \(\mbox{E}(y|X)=X\beta\) com \(\beta=(\beta_1,\cdots,\beta_p)^\top\). Suponha ainda que as colunas de \(X\) sejam padronizadas de modo que \(\frac{1}{n}\sum_{i=1}^n x_{ij}=0\) e \(\frac{1}{n}\sum_{i=1}^n x_{ij}^2=1\). O estimador Lasso \(\widehat{\beta}\) pode ser definido da seguinte maneira \[ \widehat{\beta}=\arg \max_{\beta} \left( \sum_{i=1}^n \big(y_1-x_i^\top\beta \big)^2\right), \qquad \mbox{restrito à} \qquad \sum_{j=1}^p |\beta_j|\leq s, \]

sendo \(s\geq 0\) é o parâmetro de ajuste que controla a quantidade de encolhimento.

Para o estimador OLS \(\widehat{\beta}^0 = (X^\top X)^{-1} X^\top y\) Uma escolha de parâmetros de ajuste \(s<s_0\), onde \(s_0=\sum_{j=1}^p |\widehat{\beta}_j^0|\), causará encolhimento das soluções para 0 e, finalmente, alguns coeficientes são exatamente iguais a 0. Para Valores \(s\geq s_0\) os coeficientes Lasso são iguais aos coeficientes OLS não penalizados.

Uma representação aleternativa é \[ \widehat{\beta}=\arg \min_\beta \left( \sum_{i=1}^n \big(y_i-x_i^\top \beta\big)^2 +\lambda \sum_{j=1}^p |\beta_j|\right), \] com um parâmetro de ajuste \(\lambda\geq 0\).

À medida que \(\lambda\) aumenta, as estimativas da regressão Lasso são continuamente reduzidas para zero. Então, se \(\lambda\) for muito grande, alguns coeficientes são exatamente zero. Para \(\lambda=0\) os coeficientes Lasso coincidem com a estimativa OLS.

Vejamos o procedimento no seguint exemplo. 

library(glmnet)
data(swiss) ## Lendo o arquivo de dados
head(swiss)
##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6

Dados de fertilidade Suíça e indicadores socioeconômicos em 1888. Medidas de fecundidade padronizadas e indicadores socioeconômicos para cada uma das 47 províncias francófonas da Suíça por volta de 1888. Arquivo de dados com 47 observações em 6 variáveis, cada uma delas em porcentagem, ou seja, no intervalo [0, 100].
Fertility: medida de fertilidade padronizada.
Agriculture: % de homens envolvidos na agricultura como ocupação.
Examination: % de recrutas que receberam a nota mais alta no exame do exército.
Education: % educação além da escola primária para recrutas.
Catholic: % católica da população, em oposição a protestante.
Infant.Motality: % de nascidos vivos que vivem menos de 1 ano.
Todas as variáveis, exceto Fertilidty, fornecem proporções da população.

A Suíça, em 1888, entrava em um período conhecido como transição demográfica; ou seja, sua fecundidade começava a cair do alto nível típico dos países subdesenvolvidos. Os dados coletados são para 47 províncias de língua francesa por volta de 1888.

Observação: o arquivos para todos os 182 distritos em 1888 e outros anos estão disponíveis em https://opr.princeton.edu/archive/pefp/switz.aspx. 

x.vars <- model.matrix(Fertility ~ . , swiss)[,-1]
y.var <- swiss$Fertility
lambda.seq <- 10^seq(2, -2, by = -.1)
lambda.seq
##  [1] 100.00000000  79.43282347  63.09573445  50.11872336  39.81071706
##  [6]  31.62277660  25.11886432  19.95262315  15.84893192  12.58925412
## [11]  10.00000000   7.94328235   6.30957344   5.01187234   3.98107171
## [16]   3.16227766   2.51188643   1.99526231   1.58489319   1.25892541
## [21]   1.00000000   0.79432823   0.63095734   0.50118723   0.39810717
## [26]   0.31622777   0.25118864   0.19952623   0.15848932   0.12589254
## [31]   0.10000000   0.07943282   0.06309573   0.05011872   0.03981072
## [36]   0.03162278   0.02511886   0.01995262   0.01584893   0.01258925
## [41]   0.01000000

Dividindo os dados em teste e treinamento. 

set.seed(86)
train = sample(1:nrow(x.vars), nrow(x.vars)/2)
train
##  [1]  3 29  1 46 27 23 37  2 38 26 44 22 33 15 21 43 10 28 25 19  9 32 35
x.test = (-train)
y.test = y.var[x.test]
swiss[train,]
##              Fertility Agriculture Examination Education Catholic
## Franches-Mnt      92.5        39.7           5         5    93.40
## Vevey             58.3        26.8          25        19    18.46
## Courtelary        80.2        17.0          15        12     9.96
## Rive Droite       44.7        46.6          16        29    50.43
## Paysd'enhaut      72.0        63.5           6         3     2.56
## Nyone             56.6        50.9          22        12    15.14
## Sierre            92.2        84.6           3         3    99.46
## Delemont          83.1        45.1           6         9    84.84
## Sion              79.3        63.1          13        13    96.83
## Payerne           74.2        58.1          14         8     5.23
## ValdeTravers      67.6        18.7          25         7     8.65
## Moudon            65.0        55.1          14         3     4.52
## Herens            77.3        89.7           5         2   100.00
## Cossonay          61.7        69.3          22         5     2.82
## Morges            65.5        59.8          22        10     5.23
## Val de Ruz        77.6        37.6          15         7     4.97
## Sarine            82.9        45.2          16        13    91.38
## Rolle             60.5        60.8          16        10     7.72
## Oron              72.5        71.2          12         1     2.40
## La Vallee         54.3        15.2          31        20     2.15
## Gruyere           82.4        53.3          12         7    97.67
## Entremont         69.3        84.9           7         6    99.68
## Monthey           79.4        64.9           7         3    98.22
##              Infant.Mortality
## Franches-Mnt             20.2
## Vevey                    20.9
## Courtelary               22.2
## Rive Droite              18.2
## Paysd'enhaut             18.0
## Nyone                    16.7
## Sierre                   16.3
## Delemont                 22.2
## Sion                     18.1
## Payerne                  23.8
## ValdeTravers             19.5
## Moudon                   22.4
## Herens                   18.3
## Cossonay                 18.7
## Morges                   18.0
## Val de Ruz               20.0
## Sarine                   24.4
## Rolle                    16.3
## Oron                     21.0
## La Vallee                10.8
## Gruyere                  21.0
## Entremont                19.8
## Monthey                  20.2
cv.output <- cv.glmnet(x.vars[train,], y.var[train],
                       alpha = 1, lambda = lambda.seq, 
                       nfolds = 5, type.measure = "mse")
best.lam <- cv.output$lambda.min
best.lam
## [1] 0.01


Usando este valor, vamos treinar o modelo Lasso novamente. 

# Rebuilding the model with best lamda value identified
lasso.best <- glmnet(x.vars[train,], y.var[train], alpha = 1, lambda = best.lam)
pred <- predict(lasso.best, s = best.lam, newx = x.vars[x.test,])

Por fim, combinamos os valores previstos e os valores reais para ver os dois valores lado a lado. 


final <- cbind(y.test, pred)
# Checking the first six obs
head(final)
##            y.test       s1
## Moutier      85.8 79.06310
## Neuveville   76.9 61.76096
## Porrentruy   76.1 90.75805
## Broye        83.8 75.86269
## Glane        92.4 77.58111
## Veveyse      87.1 80.81532

Para obter a lista de variáveis importantes, basta investigar os coeficientes beta do melhor modelo final. 

coef(lasso.best)
## 6 x 1 sparse Matrix of class "dgCMatrix"
##                          s0
## (Intercept)      92.8500268
## Agriculture      -0.2976592
## Examination      -0.5447304
## Education        -1.0955085
## Catholic          0.1338946
## Infant.Mortality  0.3330710

IX.3. Rede elástica


Embora o Lasso seja amplamente utilizado na seleção de variáveis, ele apresenta várias desvantagens. Zou and Hastie (2005) afirmaram que:
1. se \(p > n\), o Lasso seleciona no máximo \(n\) variáveis antes de saturar;
2. se houver um grupo de variáveis com correlação muito alta, então o Lasso tende a selecionar apenas uma variável desse grupo;
3. para a condição usual \(n > p\), se houver altas correlações entre os preditores, o desempenho de predição do Lasso é dominado pela regressão Ridge, ver Tibshirani (1996).

Zou and Hastie (2005) introduziram a rede elástica (elastic net) que combina boas características das penalidades da norma \(L_1\) e da norma \(L_2\). A rede elástica é um método de regressão regularizado que supera as limitações do Lasso. Este método é muito útil quando \(p\gg n\) ou existem muitas variáveis correlacionadas. As vantagens são:
(1) um grupo de variáveis correlacionadas pode ser selecionado sem omissões arbitrárias,
(2) o número de variáveis selecionadas não é mais limitado pelo tamanho da amostra.


IX.3.1. Rede elástica no modelo de regressão linear


Descrevemos a rede elástica no modelo de regressão linear. Por simplicidade, assumimos que os \(x_{ij}\) são padronizados de tal forma que \(\sum_{i=1}^n x_{ij}=0\) e \(\frac{1}{n}\sum_{i=1}^n x_{ij}^2=1\). A penalidade da rede elástica \(P_\alpha(\beta)\) leva à seguinte modificação do problema para obter o estimador \(\widehat{\beta}\) \[ \arg \min_\beta \left( \dfrac{1}{2n} \sum_{i=1}^n \big( y_i-x_i^\top\beta\big)^2+\lambda P_\alpha(\beta)\right), \]

onde \[ P_\alpha(\beta)=\frac{1}{2}(1-\alpha)|| \beta||_2^2+\alpha || \beta||_1 = \sum_{j=1}^p \left(\dfrac{1}{2}(1-\alpha)\beta_j^2+\alpha|\beta_j| \right)\cdot \]

A penalidade \(P_\alpha(\beta)\) é um compromisso entre a regressão Ridge e o Lasso. Se \(\alpha=0\) então o critério é a regressão Ridge e se \(\alpha=1\) o método será o Lasso. Praticamente, para pequenos \(\epsilon > 0\), a rede elástica com \(\alpha=1-\epsilon\) funciona como o Lasso, mas remove degenerescências e comportamentos erráticos de seleção de variável causado por correlação extrema. Dado um específico \(\lambda\), à medida que \(\alpha\)˛ aumenta de 0 para 1, a esparsidade das soluções de rede elástica aumenta monotonicamente de 0 para a esparsidade das soluções de Lasso.

O problema de otimização de rede elástica pode ser representado como o problema Lasso usual, usando vetores \(X\) e \(y\) modificados, conforme mostrado no exemplo a seguir.


Exemplo IX.5

Para transformar o problema de otimização rede elástica no Lasso usual, deve-se primeiro aumentar \(y\) com \(p\) zeros adicionais para obter \(\widetilde{y}=(y,0)^\top\). Em seguida, aumente \(X\) com o múltiplo da matriz \(p\times p\) identidade \(\sqrt{\lambda \alpha} \, \mbox{I}\) para obter \(\widetilde{X}=\big(X^\top, \sqrt{\lambda\alpha}\, \mbox{I} \big)^\top\).

Em seguida, defina \(\widetilde{\lambda}=\lambda(1-\alpha)\) e resolva o problema original de minimização do Lasso em termos dos novos \(\widetilde{y}\), \(\widetilde{X}\) e \(\widetilde{\lambda}\). Este novo problema é equivalente ao problema rede elástico original: \[ ||\widetilde{y}-\widetilde{X}\beta||_2^2+\widetilde{\lambda}|| \beta||_1 = \left\lVert\begin{pmatrix} y \\ 0 \end{pmatrix}-\begin{pmatrix} X\beta \\ \sqrt{\lambda\alpha}\, \mbox{I}\beta \end{pmatrix}\right\rVert_2^2 +\lambda(1-\alpha)\left\lVert\beta\right\rVert_1 \\ \qquad \qquad \qquad = \left\lVert y-X\beta\right\rVert_2^2-\lambda\alpha\left\lVert\beta\right\rVert_2^2+\lambda \left\lVert\beta\right\rVert_1-\lambda\alpha \left\lVert\beta\right\rVert_1 \\ \qquad \qquad \qquad = \left\lVert y-X\beta\right\rVert_2^2+\lambda\left( \alpha\left\lVert\beta\right\rVert_2^2+(1-\alpha) \left\lVert\beta\right\rVert_1\right), \] que é equivalente ao problema de rede elástica original.

Seguimos a idéia de Friedman et al. (2010) que usaram um algoritmo de descida de coordenadas para resolver o problema de otimização inicial.

Suponhamos ter estimativas \(\widetilde{\beta}_k\) para \(k\neq j\). Em seguida, otimizamos o problema original parcialmente em relação a \(\beta\) calculando o gradiente em \(\beta_j=\widetilde{\beta}_j\), que só existe se \(\widetilde{\beta}_j\neq 0\). Com o operador \(S(z,\gamma)\) como \[ \mbox{sgn}(z)\big(|z|-\gamma \big)^+ =\left\{ \begin{array}{rcl} z-\gamma & \mbox{se} & z>0 \quad \mbox{e} \quad \gamma<|z|, \\ z+\gamma & \mbox{se} & z<0 \quad \mbox{e} \quad \gamma<|z| \\ 0 & \mbox{se} & \gamma\geq |z| \end{array}\right., \]

pode-se mostrar que a atualização de coordenadas tem a seguinte expressão \[ \widetilde{\beta}_j =\dfrac{S\Big( \frac{1}{n}\sum_{i=1}^n x_{ij}\big(y_i-\widetilde{y}_i^{(j)} \big), \lambda\alpha\Big)}{1+\lambda (1-\alpha)}, \]

onde \(\widetilde{y}_i^{(j)}=\sum_{k\neq j} x_{ik}\widetilde{\beta}_k\) é um valor ajustado que exclui a contribuição \(x_{ij}\), portanto \(y_i-\widetilde{y}_i^{(j)}\) é um resíduo parcial para ajustar \(\beta_j\).

O algoritmo calcula a estimativa de mínimos quadrados para resíduo parcial \(y_i-\widetilde{y}_i^{(j)}\), aplica a regra de limiares suaves para executar a contribuição do Lasso para a penalidade \(P_\alpha(\beta)\). Posteriormente, um encolhimento proporcional é aplicado à penalidade Ridge. Existem vários métodos usados para atualizar a estimativa atual \(\widetilde{\beta}\). Descrevemos o método de atualização mais simples, a chamada atualização ingênua.

O resíduo parcial pode ser reescrito da seguinte maneira: \[ (y_i-\widetilde{y}_i^{(j)} = y_i-\widehat{y}_i+x_{ij}\widetilde{\beta}_j = r_i+x_{ij}\widetilde{\beta}_j, \]

com \(\widehat{y}_i\) sendo o ajuste atual e o \(r_i\) o resíduo atual. Como \(x_j\) é padronizado, portanto \[ \dfrac{1}{n}\sum_{i=1}^n x_{ij}\big( y_i-\widetilde{y}_i^{(j)}\big)=\dfrac{1}{n}\sum_{i=1}^n x_{ij}r_i + \widetilde{\beta}_j\cdot \] Observe que o primeiro termo no lado direito do novo resíduo parcial é o gradiente da perda em relação a \(\beta_j\).


O exemplo anterior é teórico. Vamos mostrar como proceder para aplicar o procedimento de rede elástica para seleção de variáveis no R.

Exemplo IX.6

Utilizaremos o exemplo Hitters no pacote ISLR. Dados da Major League Baseball das temporadas de 1986 e 1987. Hitters são os rebatedores e este conjunto de dados contém 322 observações dos jogadores da Major League nas 20 variáveis a seguir:

AtBat: Número de vezes no BAT em 1986
Hits: Número de acertos em 1986
Hmrun: Número de home runs em 1986
Runs: Número de corridas em 1986
RBI: Número de execuções em 1986
Walks: Número de caminhadas em 1986
Years: Número de anos nas principais ligas
CAtBat: Número de vezes no BAT durante sua carreira
CHits: Número de acertos durante sua carreira
CHmRun: Número de home runs durante sua carreira
CRuns: Número de corridas durante sua carreira
CRBI: Número de corridas bateu durante sua carreira
CWalks: Número de caminhadas durante sua carreira
League: Fator de níveis A e N indicando a liga de jogadores no final de 1986
Division: Fator de níveis E e W indicando a divisão de jogadores no final de 1986
PutOuts: Número de tacadas em 1986
Assists: Número de assistências em 1986
Errors; Número de erros em 1986
Salary: salário anual em 1987 no dia da abertura em milhares de dólares
NewLeague: Fator de níveis A e N indicando a liga de jogadores no início de 1987

Esse conjunto de dados foi retirado da Biblioteca Statlib, que é mantida na Universidade Carnegie Mellon. Isso faz parte dos dados que foram usados na sessão de pôsteres da seção de gráficos ASA de 1988. Os dados salariais eram originalmente da Sports Illustrated, em 20 de abril de 1987. As estatísticas de 1986 e a carreira foram obtidas na atualização da Enciclopédia de Beisebol de 1987, publicada pela Collier Books, Macmillan Publishing Company, Nova York.

Referência: James, G., Witten, D., Hastie, T., and Tibshirani, R. (2013).
An Introduction to Statistical Learning with applications in R.
https://www.statlearning.com, Springer-Verlag, New York. 

data(Hitters, package = "ISLR")

Removemos os dados ausentes, que estavam todos na variável de resposta, Salário (Salary). 

Hitters = na.omit(Hitters)
head(Hitters)
##                   AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun
## -Alan Ashby         315   81     7   24  38    39    14   3449   835     69
## -Alvin Davis        479  130    18   66  72    76     3   1624   457     63
## -Andre Dawson       496  141    20   65  78    37    11   5628  1575    225
## -Andres Galarraga   321   87    10   39  42    30     2    396   101     12
## -Alfredo Griffin    594  169     4   74  51    35    11   4408  1133     19
## -Al Newman          185   37     1   23   8    21     2    214    42      1
##                   CRuns CRBI CWalks League Division PutOuts Assists Errors
## -Alan Ashby         321  414    375      N        W     632      43     10
## -Alvin Davis        224  266    263      A        W     880      82     14
## -Andre Dawson       828  838    354      N        E     200      11      3
## -Andres Galarraga    48   46     33      N        E     805      40      4
## -Alfredo Griffin    501  336    194      A        W     282     421     25
## -Al Newman           30    9     24      N        E      76     127      7
##                   Salary NewLeague
## -Alan Ashby        475.0         N
## -Alvin Davis       480.0         A
## -Andre Dawson      500.0         N
## -Andres Galarraga   91.5         N
## -Alfredo Griffin   750.0         A
## -Al Newman          70.0         A

Como esse conjunto de dados não é particularmente grande, renunciaremos a uma divisão entre uma parte de treinamento e uma outra de teste e simplesmente usaremos todos os dados como dados de treinamento. 

library(caret); library(glmnet)

Utilizamos o pacote de funções glmnet.

Ambos quantificam o quão “grandes” os coeficientes são. Como Lasso e Ridge, a interceptação não é penalizada e o Glment cuida da padronização internamente. Os coeficientes relatados também estão na escala original. A nova penalidade é λ⋅ (1 - α) 2 vezes a penalidade da crista mais λ⋅α vezes a penalidade do laço. (Dividir a penalidade da crista por 2 é uma conveniência matemática para otimização.) Essencialmente, com a escolha correta de λ e α esses dois “coeficientes de penalidade” podem ser números positivos.

Muitas vezes, é mais útil simplesmente pensar em α como controlar a mistura entre as duas penalidades e λ controlando a quantidade de penalização. α leva valores entre 0 e 1. O uso de α = 1 fornece o laço que já vimos antes. Da mesma forma, α = 0 dá crista. Usamos esses dois antes com o GlmNet () para especificar qual método que desejávamos. Agora também permitimos valores de α entre 

set.seed(42)
cv5 = trainControl(method = "cv", number = 5)

Primeiro, configuramos nossa estratégia de validação cruzada, que será 5 vezes. Em seguida, usamos a função train() com method = “glmnet”, que está realmente ajustando a rede elástica. 

hit.elnet = train(
  Salary ~ ., data = Hitters,
  method = "glmnet",
  trControl = cv5
)
hit.elnet
## glmnet 
## 
## 263 samples
##  19 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 211, 210, 210, 211, 210 
## Resampling results across tuning parameters:
## 
##   alpha  lambda      RMSE      Rsquared   MAE     
##   0.10    0.5105642  335.1157  0.4548961  235.1758
##   0.10    5.1056419  332.4145  0.4631952  231.8598
##   0.10   51.0564193  339.4067  0.4486337  231.1383
##   0.55    0.5105642  334.8926  0.4551329  234.5333
##   0.55    5.1056419  332.7370  0.4649807  230.4455
##   0.55   51.0564193  343.5384  0.4440115  235.9464
##   1.00    0.5105642  334.8944  0.4546176  234.1018
##   1.00    5.1056419  336.0219  0.4590460  230.6486
##   1.00   51.0564193  353.2987  0.4231367  244.1034
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0.1 and lambda = 5.105642.

Observe algumas coisas com esses resultados. Primeiro, tentamos três valores de α, 0,10, 0,55 e 1. Não está totalmente claro por que o Caret não usa 0. Provavelmente usa 0,10 para ajustar um modelo próximo à crista, mas com algum potencial de escassez.

Aqui, o melhor resultado usa α = 0,10, então esse resultado está em algum lugar entre Ridge e Lasso, mas mais perto de Ridge. 

hit.elnet.int = train(
  Salary ~ . ^ 2, data = Hitters,
  method = "glmnet",
  trControl = cv5,
  tuneLength = 10
)
hit.elnet.int
## glmnet 
## 
## 263 samples
##  19 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 211, 210, 209, 212, 210 
## Resampling results across tuning parameters:
## 
##   alpha  lambda       RMSE      Rsquared   MAE     
##   0.1      0.1451882  393.2788  0.4323492  248.4650
##   0.1      0.3354037  392.0764  0.4336160  247.6321
##   0.1      0.7748260  374.8427  0.4559960  237.8825
##   0.1      1.7899484  361.0048  0.4770532  229.1920
##   0.1      4.1350131  350.3773  0.4962656  223.3445
##   0.1      9.5524165  340.4218  0.5122856  219.5706
##   0.1     22.0673211  324.6000  0.5305394  212.8962
##   0.1     50.9783738  313.1050  0.5440708  211.0142
##   0.1    117.7666553  312.6688  0.5445168  214.6348
##   0.1    272.0562482  322.5227  0.5105434  224.8678
##   0.2      0.1451882  410.0915  0.4161913  258.8741
##   0.2      0.3354037  392.0071  0.4330980  246.8087
##   0.2      0.7748260  371.4314  0.4609492  236.5669
##   0.2      1.7899484  359.4744  0.4802813  229.1897
##   0.2      4.1350131  346.8774  0.5028213  222.4905
##   0.2      9.5524165  333.6905  0.5211130  216.3838
##   0.2     22.0673211  317.6890  0.5378235  210.7943
##   0.2     50.9783738  311.1098  0.5493266  212.1244
##   0.2    117.7666553  320.7910  0.5153805  222.3877
##   0.2    272.0562482  326.1054  0.5102624  231.1196
##   0.3      0.1451882  414.4120  0.4113875  261.0514
##   0.3      0.3354037  391.2516  0.4350022  245.1933
##   0.3      0.7748260  368.1012  0.4663713  235.2365
##   0.3      1.7899484  356.2588  0.4851301  227.8150
##   0.3      4.1350131  342.3541  0.5094933  220.9939
##   0.3      9.5524165  327.8091  0.5281987  213.4270
##   0.3     22.0673211  312.8352  0.5450725  209.9155
##   0.3     50.9783738  314.3595  0.5404120  216.6930
##   0.3    117.7666553  321.0825  0.5170018  224.1626
##   0.3    272.0562482  332.1429  0.5100078  239.7011
##   0.4      0.1451882  415.8166  0.4143271  260.5632
##   0.4      0.3354037  390.3954  0.4414857  245.9071
##   0.4      0.7748260  367.3895  0.4690742  235.0114
##   0.4      1.7899484  351.6854  0.4914812  226.3313
##   0.4      4.1350131  339.2011  0.5135001  219.6427
##   0.4      9.5524165  323.3357  0.5319821  212.1000
##   0.4     22.0673211  311.1030  0.5490864  210.5599
##   0.4     50.9783738  318.0806  0.5261530  219.6268
##   0.4    117.7666553  322.1464  0.5179398  226.5787
##   0.4    272.0562482  341.4168  0.5046570  249.4220
##   0.5      0.1451882  410.7874  0.4162946  259.0418
##   0.5      0.3354037  384.6668  0.4475332  244.1645
##   0.5      0.7748260  363.9544  0.4734288  233.5289
##   0.5      1.7899484  349.8862  0.4952438  225.7635
##   0.5      4.1350131  336.9159  0.5169513  217.9557
##   0.5      9.5524165  320.5743  0.5335057  211.4521
##   0.5     22.0673211  310.5651  0.5526770  211.5045
##   0.5     50.9783738  319.0419  0.5214776  221.2613
##   0.5    117.7666553  323.7999  0.5183850  229.2370
##   0.5    272.0562482  353.2482  0.4914794  260.1245
##   0.6      0.1451882  412.1480  0.4164477  258.2229
##   0.6      0.3354037  384.8413  0.4475683  243.7498
##   0.6      0.7748260  364.2765  0.4729895  233.6009
##   0.6      1.7899484  347.8491  0.4988701  224.8212
##   0.6      4.1350131  335.7218  0.5181062  216.7491
##   0.6      9.5524165  318.5155  0.5350509  210.3275
##   0.6     22.0673211  312.5815  0.5480526  214.4155
##   0.6     50.9783738  318.3887  0.5245891  221.3194
##   0.6    117.7666553  325.7653  0.5196660  232.2374
##   0.6    272.0562482  364.7567  0.4844326  271.2375
##   0.7      0.1451882  410.1410  0.4194308  257.5552
##   0.7      0.3354037  381.4274  0.4535890  242.2628
##   0.7      0.7748260  364.2690  0.4737456  233.4038
##   0.7      1.7899484  346.5034  0.5019756  223.8508
##   0.7      4.1350131  333.0896  0.5211407  215.0685
##   0.7      9.5524165  315.9654  0.5377333  209.8188
##   0.7     22.0673211  314.5069  0.5425284  217.1827
##   0.7     50.9783738  317.9467  0.5276893  221.7685
##   0.7    117.7666553  328.4649  0.5198781  235.8628
##   0.7    272.0562482  376.4498  0.4886821  282.5303
##   0.8      0.1451882  407.7863  0.4239645  257.6782
##   0.8      0.3354037  378.7055  0.4552740  241.1232
##   0.8      0.7748260  361.6463  0.4764731  233.0779
##   0.8      1.7899484  343.3875  0.5067585  222.3726
##   0.8      4.1350131  329.8225  0.5254846  213.2195
##   0.8      9.5524165  313.4290  0.5419433  209.7046
##   0.8     22.0673211  316.0156  0.5372505  218.8972
##   0.8     50.9783738  317.7341  0.5308009  222.2297
##   0.8    117.7666553  331.9313  0.5189227  239.8751
##   0.8    272.0562482  389.4807  0.4947194  293.7957
##   0.9      0.1451882  406.1127  0.4283644  255.1716
##   0.9      0.3354037  377.6037  0.4593233  239.2795
##   0.9      0.7748260  360.0576  0.4784634  232.7823
##   0.9      1.7899484  340.9995  0.5106736  220.8525
##   0.9      4.1350131  325.6707  0.5318285  211.0865
##   0.9      9.5524165  310.4904  0.5492193  209.4210
##   0.9     22.0673211  316.2362  0.5356849  219.1407
##   0.9     50.9783738  317.6838  0.5344111  222.6642
##   0.9    117.7666553  335.9835  0.5162603  243.8591
##   0.9    272.0562482  404.4440  0.4928218  306.1844
##   1.0      0.1451882  403.5003  0.4289753  254.6877
##   1.0      0.3354037  376.1830  0.4583933  238.7983
##   1.0      0.7748260  357.1230  0.4816841  231.9271
##   1.0      1.7899484  339.8758  0.5123262  219.5750
##   1.0      4.1350131  320.7589  0.5384688  209.4310
##   1.0      9.5524165  308.5865  0.5555343  209.5209
##   1.0     22.0673211  316.8043  0.5323551  219.1355
##   1.0     50.9783738  317.6560  0.5374296  222.8821
##   1.0    117.7666553  341.3464  0.5075528  248.2012
##   1.0    272.0562482  420.5389  0.4861243  319.8683
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 1 and lambda = 9.552416.

Agora tentamos uma pesquisa de modelo muito maior. Primeiro, estamos expandindo o espaço de recursos para incluir todas as interações. Como estamos usando a regressão penalizada, não precisamos nos preocupar tanto com o excesso de ajuste. Se muitas das variáveis adicionadas não forem úteis, provavelmente usaremos um modelo próximo ao LASSO, o que faz com que muitos deles 0.

Também estamos usando uma grade de ajuste maior. Ao definir o tuneLength = 10, pesquisaremos 10 valores α e 10 valores λ para cada um. Devido a essa grade de ajuste maior, os resultados serão muito grandes.

Para lidar com isso, escrevemos uma função rápida ajudante para extrair a linha com os melhores parâmetros de ajuste.

get.best.result = function(caret.fit) {
  best = which(rownames(caret.fit$results) == rownames(caret.fit$bestTune))
  best.result = caret.fit$results[best, ]
  rownames(best.result) = NULL
  best.result
}

Em seguida, chamamos essa função no objeto treinado.

get.best.result(hit.elnet.int)
##   alpha   lambda     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 1     1 9.552416 308.5865 0.5555343 209.5209 55.15614  0.0866908 21.13291



IX.4. Lasso agrupado


O Lasso agrupadofoi introduzido pela primeira vez por Yuan ande Lin (2006) e foi motivado pelo fato de as variáveis preditoras poderem ocorrer em vários grupos e alguém poderia querer um modelo parcimonioso que use apenas alguns desses grupos. Isto é, suponha que haja \(K\) grupos e o vetor de coeficientes é estruturado da seguinte maneira \[ \beta^G = (\beta_1^\top,\cdots,\beta_K^\top)^\top\in\mathbb{R}^{\sum_k p_k}, \]

onde \(p_k\) é a dimensão do vetor de coeficientes do \(k\)-ésimo grupo, \(k=1,\cdots,K\). Um conjunto esparso de grupos é produzido, embora dentro de cada grupo todas as entradas de \(\beta_k\), \(k=1,\cdots,K\), um elemento correspondente de todo o vetor \(\beta^G\) é zero ou todos eles são diferentes de zero. O problema do Lasso agrupado pode ser formulado em geral como \[ \arg \min_{\beta\in\mathbb{R}^{\sum_k p_k}} \dfrac{1}{n}\left\lVert y-\sum_{k=1}^K X_k \beta_k\right\rVert_2^2 +\lambda \sum_{k=1}^K \sqrt{p_k}||\beta_k||_2, \]

onde \(X_k\)é o \(k\)-ésimo componente da martiz de planejamento \(X\) com colunas correspondentes aos preditores no grupo \(k\), \(\beta_k\) é o vetor de coeficiente para esse grupo e \(p_k\) é a cardinalidade do grupo, isto é, o tamanho do vetor do coeficiente que serve como um peso de equilíbrio no caso de tamanhos de grupo amplamente diferentes. É óbvio que, se os grupos consistirem em elementos únicos, ou seja, \(p_k=1\), \(\forall k\), então o problema do Lasso agrupado se reduz ao Lasso habitual.

O cálculo da solução do Lasso agrupado envolve o cálculo das condições necessárias e suficientes para que \(\widehat{\beta}^G = (\widehat{\beta}_1^\top,\cdots,\widehat{\beta}_K^\top)^\top\) seja solução do problema acima, ou seja, de \[ -X_k^\top \left( y-\sum_{k=1}^K X_k\beta_k\right)+\dfrac{\lambda\beta_k\sqrt{p_k}}{||\beta_k||}=0, \] se \(\beta_k\neq 0\); caso contrário, para \(\beta_k=0\), é válido que \[ \left\lVert X_k^\top \left( y-\sum_{l\neq k} X_l \widehat{\beta}_l \right)\right\rVert \leq \lambda\sqrt{p_k}\cdot \]

As expressões anteriores permitem calcular a solução, a chamada etapa de atualização que pode ser usada para implementar um algoritmo iterativo para resolver o problema do Lasso agrupado. A solução resultante das condições é prontamente mostrada como sendo: \[ \widehat{\beta}_k =\left( \Big(\lambda\sqrt{p_k} ||\widehat{\beta}_k ||^{-1} + X_k^\top X_k\Big)^{-1}\right)^+ X_k^\top \widehat{r}_k, \] onde os resíduos \(\widehat{r}_k\) são definidos como \[ \widehat{r}_k=y-\sum_{l\neq k} X_l \widehat{\beta}_l\cdot \]

Como um caso especial (ortonormal), quando \(X_l^\top X_l=\mbox{I}\), a solução é simplificada para \[ \widehat{\beta}_k = \left( \Big(\lambda\sqrt{p_k} ||\widehat{\beta}_k ||^{-1} + 1\Big)^{-1}\right)^+ X_k^\top \widehat{r}_k\cdot \]

Para obter uma solução completa para esse problema, Yuan and Lin (2006) sugerem o uso de um algoritmo de descida de coordenadas em bloco que aplica iterativamente as estimativas acima para \(k=1,\cdots,K\).

Meier, Van de Geer and Bühlmann (2008) estenderam o grupo Lasso ao caso de regressão logística e demonstraram a convergência de vários algoritmos para o cálculo da solução, bem como os resultados de consistência descritas para o estimador de logit de grupo Lasso. A configuração geral para esse modelo envolve uma variável de resposta binária YI 2 F0; Grupos 1G e K Preditora variável xi d .x> i1; ::::. x> ik />, tanto xi quanto yi são i.i.d., i d 1; ::::. n. Em seguida, o modelo de regressão linear logística pode ser escrito como antes:


Exemplo IX.7 

library(lars)
data("diabetes")
ajuste <- lm(diabetes[,2] ~ diabetes[,1])
summary(ajuste)
## 
## Call:
## lm(formula = diabetes[, 2] ~ diabetes[, 1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -155.829  -38.534   -0.227   37.806  151.355 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       152.133      2.576  59.061  < 2e-16 ***
## diabetes[, 1]age  -10.012     59.749  -0.168 0.867000    
## diabetes[, 1]sex -239.819     61.222  -3.917 0.000104 ***
## diabetes[, 1]bmi  519.840     66.534   7.813 4.30e-14 ***
## diabetes[, 1]map  324.390     65.422   4.958 1.02e-06 ***
## diabetes[, 1]tc  -792.184    416.684  -1.901 0.057947 .  
## diabetes[, 1]ldl  476.746    339.035   1.406 0.160389    
## diabetes[, 1]hdl  101.045    212.533   0.475 0.634721    
## diabetes[, 1]tch  177.064    161.476   1.097 0.273456    
## diabetes[, 1]ltg  751.279    171.902   4.370 1.56e-05 ***
## diabetes[, 1]glu   67.625     65.984   1.025 0.305998    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.15 on 431 degrees of freedom
## Multiple R-squared:  0.5177, Adjusted R-squared:  0.5066 
## F-statistic: 46.27 on 10 and 431 DF,  p-value: < 2.2e-16


ajuste1 <- lars(diabetes[,1],diabetes[,2], type = "lasso")
coef(ajuste1)
##              age        sex       bmi       map        tc       ldl       hdl
##  [1,]   0.000000    0.00000   0.00000   0.00000    0.0000   0.00000    0.0000
##  [2,]   0.000000    0.00000  60.11927   0.00000    0.0000   0.00000    0.0000
##  [3,]   0.000000    0.00000 361.89461   0.00000    0.0000   0.00000    0.0000
##  [4,]   0.000000    0.00000 434.75796  79.23645    0.0000   0.00000    0.0000
##  [5,]   0.000000    0.00000 505.65956 191.26988    0.0000   0.00000 -114.1010
##  [6,]   0.000000  -74.91651 511.34807 234.15462    0.0000   0.00000 -169.7114
##  [7,]   0.000000 -111.97855 512.04409 252.52702    0.0000   0.00000 -196.0454
##  [8,]   0.000000 -197.75650 522.26485 297.15974 -103.9462   0.00000 -223.9260
##  [9,]   0.000000 -226.13366 526.88547 314.38927 -195.1058   0.00000 -152.4773
## [10,]   0.000000 -227.17580 526.39059 314.95047 -237.3410  33.62827 -134.5994
## [11,]  -5.718948 -234.39762 522.64879 320.34255 -554.2663 286.73617    0.0000
## [12,]  -7.011245 -237.10079 521.07513 321.54903 -580.4386 313.86213    0.0000
## [13,] -10.012198 -239.81909 519.83979 324.39043 -792.1842 476.74584  101.0446
##            tch      ltg      glu
##  [1,]   0.0000   0.0000  0.00000
##  [2,]   0.0000   0.0000  0.00000
##  [3,]   0.0000 301.7753  0.00000
##  [4,]   0.0000 374.9158  0.00000
##  [5,]   0.0000 439.6649  0.00000
##  [6,]   0.0000 450.6674  0.00000
##  [7,]   0.0000 452.3927 12.07815
##  [8,]   0.0000 514.7495 54.76768
##  [9,] 106.3428 529.9160 64.48742
## [10,] 111.3841 545.4826 64.60667
## [11,] 148.9004 663.0333 66.33096
## [12,] 139.8579 674.9366 67.17940
## [13,] 177.0642 751.2793 67.62539



IX.5. Exercícios


  1. Obtenha o estimador Lasso explícito para o caso de um planejamento ortonormal.

  2. Otimize o valor de \(S\), de modo que o modelo ajustado no Exemplo IX.3 produça o menor residual.

  3. Otimize o valor de \(s\), de modo que o modelo ajustado no Exemplo IX.4 produça o menor residual.