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.
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.
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.
\(\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)
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.
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
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.
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.
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.
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 α entreset.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
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:
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