Introducción y preparación de r

En primer lugar, cargamos las librerías necesarias

if (!require(faraway)) install.packages('faraway', dependencies = T)
## Loading required package: faraway
if (!require(tidyverse)) install.packages('tidyverse', dependencies = T)
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
if (!require(caret)) install.packages('caret', dependencies = T)
## Loading required package: caret
## Loading required package: lattice
## 
## Attaching package: 'lattice'
## 
## The following object is masked from 'package:faraway':
## 
##     melanoma
## 
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
if (!require(leaps)) install.packages('leaps', dependencies = T)
## Loading required package: leaps
if (!require(glmnet)) install.packages('glmnet', dependencies = T)
## Loading required package: glmnet
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Loaded glmnet 4.1-7
if (!require(pls)) install.packages('pls', dependencies = T)
## Loading required package: pls
## 
## Attaching package: 'pls'
## 
## The following object is masked from 'package:caret':
## 
##     R2
## 
## The following object is masked from 'package:stats':
## 
##     loadings
library(faraway)
library(tidyverse)
library(caret)
library(leaps)
library(glmnet)
library(pls)

Regresión Ridge y Lasso

Estudiamos la base de datos

names(seatpos)
## [1] "Age"       "Weight"    "HtShoes"   "Ht"        "Seated"    "Arm"      
## [7] "Thigh"     "Leg"       "hipcenter"
nrow(seatpos)
## [1] 38
summary(seatpos)
##       Age            Weight         HtShoes            Ht       
##  Min.   :19.00   Min.   :100.0   Min.   :152.8   Min.   :150.2  
##  1st Qu.:22.25   1st Qu.:131.8   1st Qu.:165.7   1st Qu.:163.6  
##  Median :30.00   Median :153.5   Median :171.9   Median :169.5  
##  Mean   :35.26   Mean   :155.6   Mean   :171.4   Mean   :169.1  
##  3rd Qu.:46.75   3rd Qu.:174.0   3rd Qu.:177.6   3rd Qu.:175.7  
##  Max.   :72.00   Max.   :293.0   Max.   :201.2   Max.   :198.4  
##      Seated            Arm            Thigh            Leg       
##  Min.   : 79.40   Min.   :26.00   Min.   :31.00   Min.   :30.20  
##  1st Qu.: 85.20   1st Qu.:29.50   1st Qu.:35.73   1st Qu.:33.80  
##  Median : 89.40   Median :32.00   Median :38.55   Median :36.30  
##  Mean   : 88.95   Mean   :32.22   Mean   :38.66   Mean   :36.26  
##  3rd Qu.: 91.62   3rd Qu.:34.48   3rd Qu.:41.30   3rd Qu.:38.33  
##  Max.   :101.60   Max.   :39.60   Max.   :45.50   Max.   :43.10  
##    hipcenter      
##  Min.   :-279.15  
##  1st Qu.:-203.09  
##  Median :-174.84  
##  Mean   :-164.88  
##  3rd Qu.:-119.92  
##  Max.   : -30.95

A continuación dividimos la base de datos en un conjunto de entrenamiento (80%) y en otro de validación (20%)

set.seed(123) 
training.samples <- seatpos$hipcenter %>% createDataPartition(p = 0.8, list = FALSE)
train.data <- seatpos [training.samples, ] 
test.data <- seatpos [-training.samples, ]

Como paso previo generamos nuestras matrices x (variables independientes eliminando el intercepto) e y (variable dependiente)

x <- model.matrix(hipcenter ~ ., train.data)[, -1]
y <- train.data$hipcenter

Regresión Ridge

En primer lugar ajustamos un modelo de regresión ridge alfa = 0.

grid <- 10^seq(10, -2, length = 100)
ridge.mod <- glmnet(x, y, alpha = 0, lambda = grid)

Observamos que nos realizará 100 modelos con 9 variables

dim(coef(ridge.mod))
## [1]   9 100

A modo de ejemplo, podemos observar el lamda de un modelo en concreto. En este caso el modelo 40

ridge.mod$lambda[40]
## [1] 187381.7

Podemos observar gráficamente el modelo lambda 40 con las diferentes variables.

plot(ridge.mod, "lambda", label = TRUE)
abline(v=log(ridge.mod$lambda[40]), col="blue",lwd= 4,lty =3)

A continuación, a modo de ejemplo, observamos el valor de cada varialbe para el lambda 40.

coef(ridge.mod)[,40]
##   (Intercept)           Age        Weight       HtShoes            Ht 
## -1.654325e+02  2.846412e-04 -3.639223e-04 -1.435656e-03 -1.449569e-03 
##        Seated           Arm         Thigh           Leg 
## -2.945947e-03 -3.469481e-03 -3.288843e-03 -4.700651e-03

A continuación buscarmeos el modelo con el mejor lambda. Para ello usarmeos la validación cruzada y haremos la observación gráfica. El mejor lambda está entre las dos líneas verticales.

sal.cv <- cv.glmnet(x,y,alpha=0)
plot(sal.cv)

Para poder concretar cuál es el mejor, usamos la siguiente formulación:

mejor.lambda <- sal.cv$lambda.min
mejor.lambda
## [1] 44.83993

En el caso del gráfico, sería el logaritmo. Es decir, el mejor lambda sería en 3.803099

log(mejor.lambda)
## [1] 3.803099

Prediccciones del modelo Ridge con el mejor lambda

Con esta formulación podemos observar el peso de las diferentes variables.

test.data.mejor.lambda <- model.matrix(hipcenter ~.,test.data)[,-1]
coef(ridge.mod)[,which(ridge.mod$lambda==mejor.lambda)]
## 9 x 0 sparse Matrix of class "dgCMatrix"
##            
## (Intercept)
## Age        
## Weight     
## HtShoes    
## Ht         
## Seated     
## Arm        
## Thigh      
## Leg

Finalmente, mostramos una predicción con el mejor lambda y el resultado del RMSE (Error cuadrado medio) y el R2 para poder comparar después los diferentes modelos.

ridge.pred <- predict(ridge.mod, s=mejor.lambda, newx = test.data.mejor.lambda)
data.frame(RMSE = RMSE(ridge.pred, test.data$hipcenter))
caret::R2(ridge.pred, test.data$hipcenter)
##             s1
## [1,] 0.9397662

Regresión Lasso

En primer lugar ajustamos un modelo de regresión Lasso ajustaando el alfa = 1.

grid <- 10^seq(10, -2, length = 100)
lasso.mod <- glmnet(x, y, alpha = 1, lambda = grid)

Observamos que nos realizará 100 modelos con 9 variables

dim(coef(lasso.mod))
## [1]   9 100

A modo de ejemplo, podemos observar el lamda de un modelo en concreto. En este caso el modelo 90.

lasso.mod$lambda[90]
## [1] 0.1629751

Podemos observar gráficamente el modelo lambda 90 con las diferentes variables.

plot(lasso.mod, "lambda", label = TRUE)
abline(v=log(lasso.mod$lambda[90]), col="blue",lwd= 4,lty =3)

A continuación, a modo de ejemplo, observamos el valor de cada varialbe para el lambda 90.

coef(lasso.mod)[,90]
##  (Intercept)          Age       Weight      HtShoes           Ht       Seated 
## 466.72650612   0.85359135   0.06436282  -1.94985426   0.00000000   0.23759144 
##          Arm        Thigh          Leg 
##  -2.28758535  -1.07885777  -6.63336101

A continuación buscarmeos el modelo con el mejor lambda. Para ello usarmeos la validación cruzada y haremos la observación gráfica. El mejor lambda está entre las dos líneas verticales.

sal.cv <- cv.glmnet(x,y,alpha=1)
plot(sal.cv)

Para poder concretar cuál es el mejor, usamos la siguiente formulación:

mejor.lambda <- sal.cv$lambda.min
mejor.lambda
## [1] 10.85184

En el caso del gráfico, sería el logaritmo. Es decir, el mejor lambda sería en 3.710065

log(mejor.lambda)
## [1] 2.384335

Prediccciones del modelo Lasso con el mejor lambda

Con esta formulación podemos observar el peso de las diferentes variables.

test.data1 <- model.matrix(hipcenter ~.,test.data)[,-1]
coef(lasso.mod)[,which(lasso.mod$lambda==mejor.lambda)]
## 9 x 0 sparse Matrix of class "dgCMatrix"
##            
## (Intercept)
## Age        
## Weight     
## HtShoes    
## Ht         
## Seated     
## Arm        
## Thigh      
## Leg

Finalmente, mostramos una predicción con el mejor lambda y el resultado del RMSE (Error cuadrado medio) y el R2 para poder comparar después los diferentes modelos.

pred.lasso <- predict(lasso.mod, s=mejor.lambda, newx = test.data1)
data.frame(RMSE = RMSE(pred.lasso, test.data$hipcenter))
caret::R2(pred.lasso, test.data$hipcenter)
##             s1
## [1,] 0.9108594

Conclusiones de la regresión Ridge y Lasso

Tal como hemos observado, en este caso la regresión Lasso ofrece el modelo más ofrecía el RMSE más bajo, aunque el R2 más alto lo ofrecía la regresión Ridge.

Regresión RMSE R2
Ridge 14.39057 0.9397662
Lasso 14.02363 0.9108594

Componentes Principales

Observamos la base de datos mtcars y el nombre de las diferentes variables.

names(mtcars)
##  [1] "mpg"  "cyl"  "disp" "hp"   "drat" "wt"   "qsec" "vs"   "am"   "gear"
## [11] "carb"
summary(mtcars)
##       mpg             cyl             disp             hp       
##  Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0  
##  1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5  
##  Median :19.20   Median :6.000   Median :196.3   Median :123.0  
##  Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7  
##  3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0  
##  Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0  
##       drat             wt             qsec             vs        
##  Min.   :2.760   Min.   :1.513   Min.   :14.50   Min.   :0.0000  
##  1st Qu.:3.080   1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000  
##  Median :3.695   Median :3.325   Median :17.71   Median :0.0000  
##  Mean   :3.597   Mean   :3.217   Mean   :17.85   Mean   :0.4375  
##  3rd Qu.:3.920   3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000  
##  Max.   :4.930   Max.   :5.424   Max.   :22.90   Max.   :1.0000  
##        am              gear            carb      
##  Min.   :0.0000   Min.   :3.000   Min.   :1.000  
##  1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000  
##  Median :0.0000   Median :4.000   Median :2.000  
##  Mean   :0.4062   Mean   :3.688   Mean   :2.812  
##  3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000  
##  Max.   :1.0000   Max.   :5.000   Max.   :8.000

En primer lugar ajustamos el modelo de regresión con componentes principales y estudiamos dicho modelo.

modelo_Componentes <- pcr (mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb, data = mtcars, scale = TRUE , validation = "CV")
summary(modelo_Componentes)
## Data:    X dimension: 32 10 
##  Y dimension: 32 1
## Fit method: svdpc
## Number of components considered: 10
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV           6.123    2.646    2.722    2.587    2.604    2.655    2.704
## adjCV        6.123    2.635    2.709    2.569    2.585    2.634    2.679
##        7 comps  8 comps  9 comps  10 comps
## CV       2.953    3.020    3.461     3.309
## adjCV    2.916    2.978    3.395     3.243
## 
## TRAINING: % variance explained
##      1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
## X      57.60    84.10    90.07    92.77    94.99    97.09    98.42    99.23
## mpg    82.53    82.63    85.40    85.41    85.47    85.56    85.58    85.85
##      9 comps  10 comps
## X      99.76     100.0
## mpg    86.09      86.9

Para elejir el número de componentes principales, en primer lugar analizamos el error cuadrático medio de la raíz de prueba (prueba RMSE). En este caso vemos que si añadimos 2 componentes el error desciende pero que, al añadir más aumenta. Por otro lado, al analizar la varianza, observamos que hasta la inclusión de dos componentes se eleva considerablemente, pero después no tanto.

Esto lo podemos visibilizar también gráficamente con la siguiente instrucción

validationplot (modelo_Componentes)

validationplot (modelo_Componentes, val.type = "MSEP")

validationplot (modelo_Componentes, val.type = "R2")

Otra forma de estudiarlo sería con la siguiente formulación.

regfit.full.mtcars <- regsubsets(mpg ~., mtcars)
reg.summary.mtcars <- summary(regfit.full.mtcars)
reg.summary.mtcars$rsq
## [1] 0.7528328 0.8302274 0.8496636 0.8578510 0.8637377 0.8667078 0.8680976
## [8] 0.8687064

Tal como se observa, lo optimo sería añadir dos variables (hasta 0.8302274). Con la siguiente formualción observaremos este aspecto graficamente, mostraanso el RSS y el RSq ajustado.

par (mfrow = c(2, 2))
plot (reg.summary.mtcars$rss, xlab = "Number of Variables", ylab = "RSS", type = "l")
plot (reg.summary.mtcars$adjr2 , xlab = "Number of Variables", ylab = "Adjusted RSq", type = "l")
which.max(reg.summary.mtcars$adjr2 )
## [1] 5
points(2,reg.summary.mtcars$adjr2[2], col = "red", cex = 2, pch = 20)

Para conocer qué variables serían las adecuadas, recurrimos a la función regsubsets que realiza la mejor selección de subconjunto identificando el mejor modelo que contiene un número determinado de predictores

regfit.full.mtcars <- regsubsets(mpg ~., mtcars)
summary(regfit.full.mtcars)
## Subset selection object
## Call: regsubsets.formula(mpg ~ ., mtcars)
## 10 Variables  (and intercept)
##      Forced in Forced out
## cyl      FALSE      FALSE
## disp     FALSE      FALSE
## hp       FALSE      FALSE
## drat     FALSE      FALSE
## wt       FALSE      FALSE
## qsec     FALSE      FALSE
## vs       FALSE      FALSE
## am       FALSE      FALSE
## gear     FALSE      FALSE
## carb     FALSE      FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: exhaustive
##          cyl disp hp  drat wt  qsec vs  am  gear carb
## 1  ( 1 ) " " " "  " " " "  "*" " "  " " " " " "  " " 
## 2  ( 1 ) "*" " "  " " " "  "*" " "  " " " " " "  " " 
## 3  ( 1 ) " " " "  " " " "  "*" "*"  " " "*" " "  " " 
## 4  ( 1 ) " " " "  "*" " "  "*" "*"  " " "*" " "  " " 
## 5  ( 1 ) " " "*"  "*" " "  "*" "*"  " " "*" " "  " " 
## 6  ( 1 ) " " "*"  "*" "*"  "*" "*"  " " "*" " "  " " 
## 7  ( 1 ) " " "*"  "*" "*"  "*" "*"  " " "*" "*"  " " 
## 8  ( 1 ) " " "*"  "*" "*"  "*" "*"  " " "*" "*"  "*"

Elejiríamos un modelo en el que se recojan dos variables, es decir: cyl y wt

Conclusiones

Para comprobar si hemos realizado bien la elección de componentes principales, mostramos el modelo de regresión tanto con todos los componentes como con únicamente las dos variables seleccionadas

lm.fit <- lm(mpg ~ . , mtcars)
summary(lm.fit)
## 
## Call:
## lm(formula = mpg ~ ., data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4506 -1.6044 -0.1196  1.2193  4.6271 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 12.30337   18.71788   0.657   0.5181  
## cyl         -0.11144    1.04502  -0.107   0.9161  
## disp         0.01334    0.01786   0.747   0.4635  
## hp          -0.02148    0.02177  -0.987   0.3350  
## drat         0.78711    1.63537   0.481   0.6353  
## wt          -3.71530    1.89441  -1.961   0.0633 .
## qsec         0.82104    0.73084   1.123   0.2739  
## vs           0.31776    2.10451   0.151   0.8814  
## am           2.52023    2.05665   1.225   0.2340  
## gear         0.65541    1.49326   0.439   0.6652  
## carb        -0.19942    0.82875  -0.241   0.8122  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.65 on 21 degrees of freedom
## Multiple R-squared:  0.869,  Adjusted R-squared:  0.8066 
## F-statistic: 13.93 on 10 and 21 DF,  p-value: 3.793e-07
lm.fit.PCA <- lm(mpg ~ cyl + wt, mtcars)
summary(lm.fit.PCA)
## 
## Call:
## lm(formula = mpg ~ cyl + wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2893 -1.5512 -0.4684  1.5743  6.1004 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.6863     1.7150  23.141  < 2e-16 ***
## cyl          -1.5078     0.4147  -3.636 0.001064 ** 
## wt           -3.1910     0.7569  -4.216 0.000222 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.568 on 29 degrees of freedom
## Multiple R-squared:  0.8302, Adjusted R-squared:  0.8185 
## F-statistic: 70.91 on 2 and 29 DF,  p-value: 6.809e-12

Podemos observar que el R2 con el modelo propuesto es sensiblemente mejor