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))