Como exemplo, usaremos o conjunto de dados interno do R chamado mtcars. Usaremos hp como variável de resposta e as seguintes variáveis como preditores: mpg, wt (peso), drat, e qsec.
Para realizar a regressão Ridge, usaremos funções do pacote glmnet. Este pacote requer que a variável de resposta seja um vetor e o conjunto de variáveis preditoras seja da classe data.matrix.
O código a seguir mostra como definir nossos dados:
#define a variável resposta
y <- mtcars$hp
str(y)
## num [1:32] 110 110 93 110 175 105 245 62 95 123 ...
#define a matriz de variáveis preditoras
x <- data.matrix(mtcars[, c('mpg', 'wt', 'drat', 'qsec')])
head(x)
## mpg wt drat qsec
## Mazda RX4 21.0 2.620 3.90 16.46
## Mazda RX4 Wag 21.0 2.875 3.90 17.02
## Datsun 710 22.8 2.320 3.85 18.61
## Hornet 4 Drive 21.4 3.215 3.08 19.44
## Hornet Sportabout 18.7 3.440 3.15 17.02
## Valiant 18.1 3.460 2.76 20.22
Observe também que a regressão de Ridge requer que os dados sejam padronizados, de modo que cada variável preditora tenha uma média de 0 e um desvio padrão de 1.
Felizmente a função glmnet() executa automaticamente esta padronização para você. Se você já padronizou as variáveis, você pode especificar standardize=False.
library(glmnet)
#ajuste do modelo de regresão ridge
model <- glmnet(x, y, alpha = 0)
summary(model)
## Length Class Mode
## a0 100 -none- numeric
## beta 400 dgCMatrix S4
## df 100 -none- numeric
## dim 2 -none- numeric
## lambda 100 -none- numeric
## dev.ratio 100 -none- numeric
## nulldev 1 -none- numeric
## npasses 1 -none- numeric
## jerr 1 -none- numeric
## offset 1 -none- logical
## call 4 -none- call
## nobs 1 -none- numeric
Em seguida, identificaremos o valor lambda que produz o menor erro quadrático médio de teste (MSE) usando validação cruzada k-fold. Felizmente, o pacote glmnet tem a função cv.glmnet() que executa automaticamente a validação cruzada k-fold usando k = 10 folds.
#executa a validação cruzada k-fold para encontrar o valor lambda ideal
cv_model <- cv.glmnet(x, y, alpha = 0)
#encontra o valor lambda ideal que minimiza o MSE de teste
best_lambda <- cv_model$lambda.min
best_lambda
## [1] 15.99555
#produz gráfico de teste MSE por valor lambda
plot(cv_model)
O valor lambda que minimiza o MSE de teste acaba sendo 10.04567. Por fim, podemos analisar o modelo final produzido pelo valor ótimo de lambda. Podemos usar o seguinte código para obter as estimativas de coeficiente para este modelo:
#encontrando os coeficientes do melhor modelo
best_model <- glmnet(x, y, alpha = 0, lambda = best_lambda)
coef(best_model)
## 5 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 466.378811
## mpg -3.252871
## wt 17.876523
## drat -3.348604
## qsec -16.797192
Também podemos produzir um gráfico para visualizar como as estimativas de coeficiente mudaram como resultado do aumento de lambda:
#produze gráfico Ridge trace plot
plot(model, xvar = "lambda")
grid()
Por fim, podemos calcular o R-quadrado do modelo nos dados de treinamento:
#utilizado para ajustar o melhor modelo para fazer predições
y_predicted <- predict(model, s = best_lambda, newx = x)
#encontrando SST e SSE
sst <- sum((y - mean(y))^2)
sse <- sum((y_predicted - y)^2)
#encontrando o R-quadrado
rsq <- 1 - sse/sst
rsq
## [1] 0.7921595
O R-quadrado acaba sendo 0.8047349. Ou seja, o melhor modelo foi capaz de explicar 80% da variação nos valores de resposta dos dados de treinamento.
library(caret)
set.seed(1)
modelo1 <- train(
hp ~ .,
data = mtcars,
method = 'glmnet',
lambda = best_lambda
)
#obtendo os coeficientes domodelo final
coeficientes <- coef(modelo1$finalModel, modelo1$bestTune$lambda)
coeficientes
## 11 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 181.4098554
## mpg -1.4548172
## cyl 6.8958394
## disp 0.1386153
## drat .
## wt 2.3973643
## qsec -7.4003727
## vs .
## am .
## gear 3.7273806
## carb 10.8394679
#lista com os nomes dos coeficientes selecionados
coeficientes[which(coeficientes != 0),0]
## 8 x 0 sparse Matrix of class "dgCMatrix"
##
## (Intercept)
## mpg
## cyl
## disp
## wt
## qsec
## gear
## carb
plot(varImp(modelo1))