6 Selección y Regularización
Librerías utilizadas
6.1 Selección
A pesar de su simplicidad, los modelos lineales tienen ventajas de interpretabilidad y con frecuencia muestran buen rendimiento en predicción
\[ Y = \beta_0+\beta_1X_1+\cdots+\beta_pX_p+\epsilon. \]
Esta sección aborda el cambiar el ajuste de mínimos cuadrados ordinarios (OLS) tradicional por métodos alternativos que se apalancan en OLS.
¿Por que usar alternativas a OLS?
- Predictibilidad: Si la relación entre las variables es aproximadamente lineal y si \(n\gg p\) (el número de observaciones es más grande que el número de variables), entonces OLS tiene poco sesgo y rinde bien en datos test. Sin embargo, si \(n\) no es tan grande respecto a \(p\), puede haber mucha variabilidad en OLS, lo que resulta en sobreajuste y mal rendimiento en test.
Algunos de los problemas que surgen cuando \(n\approx p\):
- Sobreajuste (Overfitting)
El error cuadrático medio (MSE) en el conjunto de entrenamiento se define como:
\[ MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 \]
donde:
- \(y_i\) es el valor observado,
- \(\hat{y}_i\) es la predicción del modelo,
- \(n\) es el número de observaciones.
Cuando el número de variables \(p\) se aproxima a \(n\), el modelo puede ajustarse casi perfectamente, lo que minimiza el MSE en el conjunto de entrenamiento, pero aumenta el error en nuevos datos (error de generalización). Este es el problema del sobreajuste.
- Multicolinealidad
El problema de multicolinealidad se detecta cuando las variables explicativas están altamente correlacionadas entre sí. Un indicador común es el Factor de Inflación de la Varianza (VIF) para cada variable \(j\), que se define como:
\[ VIF_j = \frac{1}{1 - R_j^2} \]
donde \(R_j^2\) es el coeficiente de determinación de la regresión de la variable \(X_j\) sobre todas las demás variables.
Cuando \(VIF_j\) es grande (mayor a 10, por ejemplo), indica una alta colinealidad. Esto causa inestabilidad en las estimaciones de los coeficientes \(\hat{\beta}_j\), lo que puede hacer que el modelo sea sensible a pequeños cambios en los datos.
- Modelo No Identificable
En regresión lineal, los coeficientes de los parámetros \(\boldsymbol{\beta} = (\beta_1, \beta_2, \dots, \beta_p)\) se obtienen resolviendo el sistema de ecuaciones lineales:
\[ \boldsymbol{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon} \]
donde:
- \(\boldsymbol{y}\) es el vector de observaciones (de dimensión \(n \times 1\)),
- \(\mathbf{X}\) es la matriz de diseño (de dimensión \(n \times p\)),
- \(\boldsymbol{\beta}\) es el vector de coeficientes,
- \(\boldsymbol{\epsilon}\) es el término de error.
Para resolver \(\boldsymbol{\beta}\), necesitamos invertir la matriz \(\mathbf{X}^T \mathbf{X}\). Si \(p \geq n\), la matriz \(\mathbf{X}^T \mathbf{X}\) no es invertible, lo que significa que el sistema es no identificable y no existe una única solución para \(\boldsymbol{\beta}\).
\[ \text{Si } \det(\mathbf{X}^T \mathbf{X}) = 0, \text{ entonces no hay solución única.} \]
- Varianza Alta
Cuando hay más variables que observaciones, el error estándar de los coeficientes de regresión se incrementa. La varianza de los coeficientes \(\hat{\beta}_j\) está dada por:
\[ \text{Var}(\hat{\beta}_j) = \sigma^2 \left( (\mathbf{X}^T \mathbf{X})^{-1} \right)_{jj} \]
donde \(\sigma^2\) es la varianza del error residual. Si \(\mathbf{X}^T \mathbf{X}\) es mal condicionada (debido a la colinealidad o alta dimensionalidad), los elementos de la matriz inversa \((\mathbf{X}^T \mathbf{X})^{-1}\) pueden ser muy grandes, lo que resulta en una varianza elevada para \(\hat{\beta}_j\). Esto significa que los coeficientes son altamente sensibles a los cambios en los datos.
- Interpretabilidad: Incluir variables independientes sin asociación con la variable dependiente en un modelo de regresión resulta en incluir complejidades innecesarias en el modelo. Si son removidas (se fijan sus coeficientes en cero), se puede obtener un modelo que es más interpretable. En OLS hay muy poca probabilidad de que se tengan coeficientes iguales a cero. Al eliminar variables estamos haciendo selección de variables o feature selection: selección de variables, contracción o regularización y reducción de dimensionalidad.
6.1.1 Selección de variables
6.1.1.1 Selección de las mejores variables
| Algoritmo: Selección de las mejores variables |
|---|
| 1. Sea \(\mathcal{M_0}\) el modelo nulo (sin predictores). Este modelo solo predice la media muestral para cada observación. |
| 2. Para \(k = 1,2,\ldots,p\): |
a. Ajusta todos \({p \choose k}\) modelos que contienen exactamente \(k\) predictores.
b. Selecciona el mejor de los \({p \choose k}\) modelos, nómbralo \(\mathcal{M_k}\). El mejor modelo se define como aquel que tiene el mejor \(RRS\) (suma de residuos al cuadrado), o de manera equivalente el más alto \(R^2\).
|
| 3. Selecciona el mejor modelo de los \(\mathcal{M_0},\ldots,\mathcal{M_k}\) usando el error de predicción en un conjunto de validación, \(C_p\), AIC, BIC o \(R^2\) ajustado. O usa el método de validación cruzada. |
Nota que este algoritmo no puede ser aplicado cuando \(p\) es grande.
También puede tener problemas de estadísticos cuando \(p\) es grande. Cuanto mayor sea el espacio de búsqueda, mayor será la posibilidad de encontrar modelos que se vean bien en los datos de entrenamiento, aunque podrían no tener ningún poder predictivo sobre datos futuros.
Por lo tanto, un espacio de búsqueda enorme puede dar lugar a un sobreajuste y a una alta varianza de las estimaciones de los coeficientes.
Por los motivos antes listado, los métodos stepwise o paso a paso son más atractivos.
\(C_p\) de Mallow
\[ C_p = \frac{1}{n}(RSS+2d\hat{\sigma}^2) \]
donde \(d\) es el total de parámetros usados y \(\hat{\sigma}^2\) es la estimación de la varianza del error \(\epsilon\).
Criterio AIC
Se puede usar en modelos donde el ajuste se realiza con máxima verosimilitud:
\[ AIC = -2\text{log}L+2d \]
Donde \(L\) es el valor máximo de la función de máxima verosimilitud del modelo.
\[ BIC = \frac{1}{n}(RSS+\text{log}(n)d\hat{\sigma}^2) \]
donde \(n\) es el número de observaciones.
\(R^2\) ajustado
\[ R^2_{adj} = 1-\frac{RSS/(n-d-1)}{TSS/(n-1)} \]
donde \(TSS\) es la suma total de cuadrados.
6.1.1.2 Selección stepwise hacia adelante
| Algoritmo: Selección stepwise hacia adelante |
|---|
| 1. Sea \(\mathcal{M_0}\) el modelo nulo (sin predictores). Este modelo solo predice la media muestral para cada observación. |
| 2. Para \(k = 0,2,\ldots,p-1\): |
a. Considera todos los \(p-k\) modelos que aumentan predictores en \(\mathcal{M_k}\) con un predictor adicional.
b. Elige el mejor de estos \(p-k\) modelos, y nómbralo \(\mathcal{M_{k+1}}\). El mejor modelo se define como aquel que tiene el mejor \(RRS\), o el más alto \(R^2\).
|
| 3. Selecciona el mejor modelo de los \(\mathcal{M_0},\ldots,\mathcal{M_p}\) usando el error de predicción en un conjunto de validación, \(C_p\), AIC, BIC o \(R^2\) ajustado. O usa el método de validación cruzada. |
6.1.1.3 Selección stepwise hacia atrás
| Algoritmo: Selección stepwise hacia atrás |
|---|
| 1. Sea \(\mathcal{M_p}\) el modelo completo, que contiene todos los \(p\) predictores. |
| 2. Para \(k = p,p-1,\ldots,1\): |
a. Considera todos los \(k\) modelos que contienen todos menos uno de los predictores en en \(\mathcal{M_k}\), para un total de \(k-1\) predictores.
b. Elige el mejor de estos \(k\) modelos, y nómbralo \(\mathcal{M_{k-1}}\). El mejor modelo se define como aquel que tiene el mejor \(RRS\), o el más alto \(R^2\).
|
| 3. Selecciona el mejor modelo de los \(\mathcal{M_0},\ldots,\mathcal{M_p}\) usando el error de predicción en un conjunto de validación, \(C_p\), AIC, BIC o \(R^2\) ajustado. O usa el método de validación cruzada. |
6.2 Regularización
La regularización es un método que nos permite restringir el proceso de estimación, se usa para evitar un posible sobre ajuste del modelo (overfitting).
Una forma de hacer regularización es donde los coeficientes de variables en el modelo no sean muy grandes (regresión ridge). Otra forma es contraer los coeficientes, llegando a tener coeficientes iguales a cero (regresión lasso).
Como regla general, el segundo enfoque suele ser mejor que el primero.
Los coeficientes regularizados se obtiene usando una función de penalidad \(p(\boldsymbol{\alpha})\) para restringir el tamaño del vector de coeficientes \(\boldsymbol{\alpha} = (\alpha_1,\cdots,\alpha_M)^T\) del predictor \(f(\mathbf{x}) = \sum_{j = 1}^{M} \alpha_jC_j(\mathbf{x})\). Los coeficientes penalizados son obtenidos como solución al problema de minimización:
\[ \hat{\boldsymbol{\alpha}}_{\lambda}= \underset{\boldsymbol{\alpha}}{\mathrm{argmin}} = \left\{\sum_{i = 1}^{n}L(y_i,f(\mathbf{x}_i))+\lambda p(\boldsymbol{\alpha}) \right\}, \] donde \(L\) es una función de pérdida y \(\lambda>0\) es un parámetro de regularización también conocido como ratio de aprendizaje.
Existen diferentes opciones para funciones de pérdida:
- Exponencial:
\[ L(y,f(\mathbf{x})) = e^{-yf(\mathbf{x})}, y \in \{-1,+1\}. \]
- Logística:
\[ L(y,f_t(\mathbf{x})) = \text{log}\{1+ e^{-2yf_t(\mathbf{x})}\},y \in \{-1,+1\}. \]
- Error cuadrático:
\[ L(y,f_t(\mathbf{x})) = \frac{1}{2}(y-f_t(\mathbf{x}))^2, y \in \mathcal{R}. \]
- Error absoluto:
\[ L(y,f_t(\mathbf{x})) = |y-f_t(\mathbf{x})|, y \in \mathcal{R}. \]
- Huber
Hay dos tipos de funciones de penalidad:
- \(L_2\): esta función de penalidad restringe la suma de cuadrados de los coeficientes,
\[ p_2(\boldsymbol{\alpha})=\sum_{j = 1}^{M} \alpha_j^2. \] Cuando \(L\) combinado es una combinación convexa y usamos la pérdida de error cuadrático, el predictor de regresión penalizado óptimo es el estimador de regresión ridge.
- \(L_1\): Los coeficientes se restringen tal que su suma de valores absolutos,
\[ p_1(\boldsymbol{\alpha})=\sum_{j = 1}^{M} |\alpha_j|. \] sea menor que un valor dado. La evidencia empírica sugiere que la penalización \(L_1\) (lasso) funciona mejor cuando hay un número pequeño o mediano de coeficientes verdaderos de tamaño moderado.
6.2.1 Regresión Ridge
En mínimos cuadrados ordinarios, las estimaciones de \(\beta_0,\beta_1,\ldots,\beta_p\) se obtienen minimizando
\[ \text{RSS} = \sum_{i = 1}^n\left(y_i-\beta_0-\sum_{j=1}^p \beta_jx_{ij}\right)^2. \]
La regresión ridge es muy similar al enfoque de mínimos cuadrados, pero los coeficientes de la regresión ridge \(\hat{\beta}^R\) son los valores que minimizan
\[ \sum_{i = 1}^n\left(y_i-\beta_0-\sum_{j=1}^p \beta_jx_{ij}\right)^2+\lambda\sum_{j = 1}^{p} \beta_j^2 = \text{RSS}+\lambda\sum_{j = 1}^{p} \beta_j^2, \]
donde \(\lambda\) es un parámetro ajustable que se determina de manera separada. La influencia de la regularización se controla con \(\lambda\). Valores altos de \(\lambda\) significa más regularización y modelos más simples. Sin embargo, la regresión ridge siempre generará un modelo que incluya todos los predictores. Incrementar el valor de \(\lambda\) tenderá a reducir las magnitudes de los coeficientes, pero no resultará en la exclusión de ninguna de las variables. Notemos que la penalidad no se aplica a \(\beta_0\).
6.2.2 Regresión Lasso
\[ \sum_{i = 1}^n\left(y_i-\beta_0-\sum_{j=1}^p \beta_jx_{ij}\right)^2+\lambda\sum_{j = 1}^{p} |\beta_j| = \text{RSS}+\lambda\sum_{j = 1}^{p} |\beta_j|, \]
La penalización \(L_1\) tiene el efecto de forzar algunas de las estimaciones de coeficientes a ser exactamente iguales a cero cuando el parámetro \(\lambda\) de ajuste es suficientemente grande. Por lo tanto, lasso realiza selección de variables. Como resultado, los modelos generados a partir de lasso son generalmente mucho más fáciles de interpretar que los producidos por regresión ridge. Decimos que lasso produce modelos dispersos (sparse), es decir, modelos que involucran solo un subconjunto de variables predictoras.
6.2.2.1 Ejemplo: datos de crimen
El conjunto de datos es obtenido de (https://archive.ics.uci.edu/ml/datasets/Communities+and+Crime+Unnormalized) y se trata de comunidades en EE.UU. Son datos combinados de datos socioeconómicos del censo de 1990, datos de aplicación de la ley de la encuesta de 1990 sobre gestión de la aplicación de la ley y estadísticas administrativas y datos sobre delitos de la UCR del FBI de 1995.
A continuación se importan y limpian los datos:
# Limpieza de datos
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/CommViolPredUnnormalizedData.txt"
crime = rio::import(uu, na.strings='?')
names(crime) <- c('communityname', 'state', 'countyCode', 'communityCode', 'fold', 'population', 'householdsize', 'racepctblack', 'racePctWhite', 'racePctAsian', 'racePctHisp', 'agePct12t21', 'agePct12t29', 'agePct16t24', 'agePct65up', 'numbUrban', 'pctUrban', 'medIncome', 'pctWWage', 'pctWFarmSelf', 'pctWInvInc', 'pctWSocSec', 'pctWPubAsst', 'pctWRetire', 'medFamInc', 'perCapInc', 'whitePerCap', 'blackPerCap', 'indianPerCap', 'AsianPerCap', 'OtherPerCap', 'HispPerCap', 'NumUnderPov', 'PctPopUnderPov', 'PctLess9thGrade', 'PctNotHSGrad', 'PctBSorMore', 'PctUnemployed', 'PctEmploy', 'PctEmplManu', 'PctEmplProfServ', 'PctOccupManu', 'PctOccupMgmtProf', 'MalePctDivorce', 'MalePctNevMarr', 'FemalePctDiv', 'TotalPctDiv', 'PersPerFam', 'PctFam2Par', 'PctKids2Par', 'PctYoungKids2Par', 'PctTeen2Par', 'PctWorkMomYoungKids', 'PctWorkMom', 'NumKidsBornNeverMar', 'PctKidsBornNeverMar', 'NumImmig', 'PctImmigRecent', 'PctImmigRec5', 'PctImmigRec8', 'PctImmigRec10', 'PctRecentImmig', 'PctRecImmig5', 'PctRecImmig8', 'PctRecImmig10', 'PctSpeakEnglOnly', 'PctNotSpeakEnglWell', 'PctLargHouseFam', 'PctLargHouseOccup', 'PersPerOccupHous', 'PersPerOwnOccHous', 'PersPerRentOccHous', 'PctPersOwnOccup', 'PctPersDenseHous', 'PctHousLess3BR', 'MedNumBR', 'HousVacant', 'PctHousOccup', 'PctHousOwnOcc', 'PctVacantBoarded', 'PctVacMore6Mos', 'MedYrHousBuilt', 'PctHousNoPhone', 'PctWOFullPlumb', 'OwnOccLowQuart', 'OwnOccMedVal', 'OwnOccHiQuart', 'OwnOccQrange', 'RentLowQ', 'RentMedian', 'RentHighQ', 'RentQrange', 'MedRent', 'MedRentPctHousInc', 'MedOwnCostPctInc', 'MedOwnCostPctIncNoMtg', 'NumInShelters', 'NumStreet', 'PctForeignBorn', 'PctBornSameState', 'PctSameHouse85', 'PctSameCity85', 'PctSameState85', 'LemasSwornFT', 'LemasSwFTPerPop', 'LemasSwFTFieldOps', 'LemasSwFTFieldPerPop', 'LemasTotalReq', 'LemasTotReqPerPop', 'PolicReqPerOffic', 'PolicPerPop', 'RacialMatchCommPol', 'PctPolicWhite', 'PctPolicBlack', 'PctPolicHisp', 'PctPolicAsian', 'PctPolicMinor', 'OfficAssgnDrugUnits', 'NumKindsDrugsSeiz', 'PolicAveOTWorked', 'LandArea', 'PopDens', 'PctUsePubTrans', 'PolicCars', 'PolicOperBudg', 'LemasPctPolicOnPatr', 'LemasGangUnitDeploy', 'LemasPctOfficDrugUn', 'PolicBudgPerPop', 'murders', 'murdPerPop', 'rapes', 'rapesPerPop', 'robberies', 'robbbPerPop', 'assaults', 'assaultPerPop', 'burglaries', 'burglPerPop', 'larcenies', 'larcPerPop', 'autoTheft', 'autoTheftPerPop', 'arsons', 'arsonsPerPop', 'ViolentCrimesPerPop', 'nonViolPerPop')
# quitamos variables con pocos datos o poca relevancia
columns_to_keep <- c(5, 6, 11:26,32:103,145)+1
crime = na.omit(crime[,columns_to_keep])Ahora escalamos los datos y los separamos en entrenamiento y prueba:
Escalamos los datos con el comando preProcess de la librería caret:
Separamos los datos en entrenamiento y prueba, 75% y 25% respectivamente con el comando createDataPartition (el parámetro list indica si se desea una lista como resultado):
set.seed(8519)
index <- createDataPartition(y, p=0.75, list=FALSE)
X_train <- X[ index, ]
X_test <- X[-index, ]
y_train <- y[index]
y_test <- y[-index]La librearía glmnet contiene la función del mismo nombre con la que se ajustan modelos GLM con regularización. El parámetro alpha indica el tipo de regularización. Si alpha = 1 es regularización lasso, si alpha = 0 es ridge.
lasso <- caret::train(y= y_train,
x = X_train,
method = 'glmnet',
tuneGrid = expand.grid(alpha = 1, lambda = 1)
)
ridge <- caret::train(y = y_train,
x = X_train,
method = 'glmnet',
tuneGrid = expand.grid(alpha = 0, lambda = 1)
) Realizamos las predicciones usando los modelos estimados:
predictions_lasso <- lasso %>% predict(X_test)
predictions_ridge <- ridge %>% predict(X_test)
# Se imprimen los R2:
data.frame(
Ridge_R2 = R2(predictions_ridge, y_test),
Lasso_R2 = R2(predictions_lasso, y_test)
)## Ridge_R2 Lasso_R2
## 1 0.6797983 0.5988638Revisamos los errores:
# Se imprime los RMSE
data.frame(
Ridge_RMSE = RMSE(predictions_ridge, y_test) ,
Lasso_RMSE = RMSE(predictions_lasso, y_test)
)## Ridge_RMSE Lasso_RMSE
## 1 376.9781 432.2993Veamos los coeficientes obtenidos:
data.frame(
as.data.frame.matrix(coef(lasso$finalModel, lasso$bestTune$lambda)),
as.data.frame.matrix(coef(ridge$finalModel, ridge$bestTune$lambda))
) %>%
rename(Lasso_coef = s1, Ridge_coef = s1.1)## Lasso_coef Ridge_coef
## (Intercept) 1.130420e+03 1132.6080108
## population 0.000000e+00 -20.1633173
## householdsize 5.666758e+02 63.9252079
## agePct12t21 -1.176333e+02 13.3219464
## agePct12t29 -1.313944e+02 -64.1113641
## agePct16t24 0.000000e+00 3.8574424
## agePct65up 0.000000e+00 -2.2223517
## numbUrban 0.000000e+00 -18.9127556
## pctUrban 9.330450e+00 18.0719286
## medIncome -4.036756e+02 -40.1100835
## pctWWage -3.058135e+02 -39.5405590
## pctWFarmSelf -9.881609e+01 -60.7961486
## pctWInvInc -4.248567e+01 -116.9446940
## pctWSocSec -3.671205e+02 -52.0449366
## pctWPubAsst 3.389220e+01 53.6143467
## pctWRetire 8.119208e+01 16.3133640
## medFamInc 0.000000e+00 -28.1067377
## perCapInc -1.123051e+02 -33.5958631
## whitePerCap 1.251772e+02 109.9367670
## NumUnderPov -6.926150e+02 -52.9742410
## PctPopUnderPov -7.922162e+00 -24.7666068
## PctLess9thGrade -6.074581e+01 -55.7898435
## PctNotHSGrad -1.443651e+02 -36.5525729
## PctBSorMore 2.130984e+02 53.7561098
## PctUnemployed 4.605313e+01 45.0329162
## PctEmploy 1.847667e+01 -34.3921946
## PctEmplManu -1.972805e+02 -98.7474609
## PctEmplProfServ 0.000000e+00 21.6229637
## PctOccupManu 4.413362e+02 104.6704919
## PctOccupMgmtProf 1.561198e+02 39.3165819
## MalePctDivorce 0.000000e+00 60.3490606
## MalePctNevMarr 4.990046e+00 29.3530136
## FemalePctDiv -3.961422e+01 -16.0748567
## TotalPctDiv 0.000000e+00 13.7176081
## PersPerFam -2.291248e+02 -3.2830146
## PctFam2Par 0.000000e+00 -83.9731317
## PctKids2Par -6.569872e+02 -150.1979018
## PctYoungKids2Par 6.963574e+01 -23.3918684
## PctTeen2Par 1.088942e+02 -5.9783711
## PctWorkMomYoungKids -1.005379e+02 0.9467168
## PctWorkMom 9.220252e+01 32.3961499
## NumKidsBornNeverMar 1.990928e+02 -13.7559394
## PctKidsBornNeverMar 5.287453e-04 150.8676448
## NumImmig 2.252002e+02 61.0197149
## PctImmigRecent 4.719160e+01 -16.7025432
## PctImmigRec5 -1.030924e+02 -60.9153121
## PctImmigRec8 -1.485752e+02 -8.1538351
## PctImmigRec10 2.976061e+02 104.4518321
## PctRecentImmig -1.202809e+02 -37.2421608
## PctRecImmig5 0.000000e+00 -3.0547671
## PctRecImmig8 0.000000e+00 2.1847088
## PctRecImmig10 -2.604227e+02 30.9523183
## PctSpeakEnglOnly 0.000000e+00 46.8223834
## PctNotSpeakEnglWell 2.643543e+00 25.5978791
## PctLargHouseFam 5.067664e+02 -7.7171054
## PctLargHouseOccup -8.214177e+02 -60.9426156
## PersPerOccupHous 0.000000e+00 -0.2984062
## PersPerOwnOccHous 1.854213e+02 63.0225693
## PersPerRentOccHous -2.851028e+02 -45.3390185
## PctPersOwnOccup 0.000000e+00 24.9957897
## PctPersDenseHous 4.046975e+02 77.2752138
## PctHousLess3BR 3.997264e+01 49.3260040
## MedNumBR -4.049265e+01 -13.3199287
## HousVacant 3.060022e+02 75.2517723
## PctHousOccup 1.590144e+01 -7.5891021
## PctHousOwnOcc 4.838100e+01 14.5887569
## PctVacantBoarded 3.827256e+01 55.2052686
## PctVacMore6Mos 2.473944e+01 19.6278759
## MedYrHousBuilt 1.522281e+02 93.1131168
## PctHousNoPhone 2.271944e+01 26.8616927
## PctWOFullPlumb -1.001378e+02 -30.6673892
## OwnOccLowQuart 0.000000e+00 -31.6475492
## OwnOccMedVal 0.000000e+00 -7.0279234
## OwnOccHiQuart -1.142313e+02 -31.6718068
## OwnOccQrange -2.594090e+01 -23.6864496
## RentLowQ -4.709927e+02 -58.4116395
## RentMedian 2.964951e+01 36.0966448
## RentHighQ -7.345267e+01 -14.2412512
## RentQrange -3.348517e-01 67.1453388
## MedRent 8.418653e+02 75.6144996
## MedRentPctHousInc -4.415355e+01 -12.4628225
## MedOwnCostPctInc -3.418865e+01 35.6726206
## MedOwnCostPctIncNoMtg -5.692937e+01 -80.2774801
## NumInShelters -2.558319e+01 25.6507248
## NumStreet 7.019079e+01 25.6110485
## PctForeignBorn 3.040272e+02 72.1198686
## PctBornSameState -1.064257e+01 -23.3078839
## PctSameHouse85 5.027221e+01 44.9946815
## PctSameCity85 1.737350e+01 2.5094687
## PctSameState85 4.649409e+01 17.7969544
## LemasSwornFT -1.195124e+01 -31.4780697Ahora, elegimos el parámetro de regularización con la ayuda de tuneGrid. Los modelos con la puntuación \(R^2\) cuadrado más alta nos darán los mejores parámetros.
parameters <- c(seq(0.1, 2, by =0.1) , seq(2, 5, 0.5) , seq(5, 500, 5))
lasso <- caret::train(y = y_train,
x = X_train,
method = 'glmnet',
tuneGrid = expand.grid(alpha = 1, lambda = parameters) ,
metric = "Rsquared"
)
ridge <- caret::train(y = y_train,
x = X_train,
method = 'glmnet',
tuneGrid = expand.grid(alpha = 0, lambda = parameters),
metric = "Rsquared"
)
linear <- caret::train(y = y_train,
x = X_train,
method = 'lm',
metric = "Rsquared"
)
print(paste0('Lasso (mejores parámetros): ' , lasso$finalModel$lambdaOpt))## [1] "Lasso (mejores parámetros): 230"## [1] "Ridge (mejores parámetros): 500"Revisamos el \(R^2\) y el \(\text{RMSE}\):
predictions_lasso <- lasso %>% predict(X_test)
predictions_ridge <- ridge %>% predict(X_test)
predictions_lin <- linear %>% predict(X_test)
data.frame(
Ridge_R2 = R2(predictions_ridge, y_test),
Lasso_R2 = R2(predictions_lasso, y_test),
Linear_R2 = R2(predictions_lin, y_test)
)## Ridge_R2 Lasso_R2 Linear_R2
## 1 0.6737963 0.6189773 0.4857502data.frame(
Ridge_RMSE = RMSE(predictions_ridge, y_test) ,
Lasso_RMSE = RMSE(predictions_lasso, y_test) ,
Linear_RMSE = RMSE(predictions_lin, y_test)
)## Ridge_RMSE Lasso_RMSE Linear_RMSE
## 1 375.0716 431.1954 508.1641Veamos los coeficientes estimados:
## [1] "Los mejores coeficientes estimados"data.frame(
ridge = as.data.frame.matrix(coef(ridge$finalModel, ridge$finalModel$lambdaOpt)),
lasso = as.data.frame.matrix(coef(lasso$finalModel, lasso$finalModel$lambdaOpt)),
linear = (linear$finalModel$coefficients)
) %>% rename(lasso = s1, ridge = s1.1)## lasso ridge linear
## (Intercept) 1134.9431646 1147.6506 1122.658115
## population 0.3606813 0.0000 -10871.577689
## householdsize 14.8157146 0.0000 840.365581
## agePct12t21 3.3151221 0.0000 -194.637075
## agePct12t29 -27.4810827 0.0000 -177.592993
## agePct16t24 -11.2192172 0.0000 -60.593345
## agePct65up -2.5271129 0.0000 -25.036687
## numbUrban 0.5120029 0.0000 11053.501418
## pctUrban 22.0555885 0.0000 -33.934401
## medIncome -7.9746127 0.0000 -1248.079959
## pctWWage -8.4534822 0.0000 -474.962612
## pctWFarmSelf -40.9593232 0.0000 -132.054694
## pctWInvInc -51.7442717 0.0000 -2.825268
## pctWSocSec -13.1436408 0.0000 -612.355509
## pctWPubAsst 29.1270431 0.0000 150.967869
## pctWRetire -7.2210071 0.0000 61.432506
## medFamInc -8.2261429 0.0000 652.356848
## perCapInc -1.3994315 0.0000 -148.710400
## whitePerCap 42.6074709 0.0000 104.460205
## NumUnderPov -3.5031683 0.0000 -1345.830586
## PctPopUnderPov 11.0609017 0.0000 -116.646514
## PctLess9thGrade -21.3138453 0.0000 -214.436395
## PctNotHSGrad -0.3011144 0.0000 -144.805596
## PctBSorMore 5.7726396 0.0000 151.457508
## PctUnemployed 28.2245967 0.0000 -10.779315
## PctEmploy -22.1852667 0.0000 68.897050
## PctEmplManu -55.7591311 0.0000 -159.479930
## PctEmplProfServ 15.6349282 0.0000 -49.320027
## PctOccupManu 5.7162408 0.0000 489.784243
## PctOccupMgmtProf 5.6746190 0.0000 405.405934
## MalePctDivorce 41.3376924 0.0000 162.111945
## MalePctNevMarr 21.0935986 0.0000 61.978109
## FemalePctDiv 29.5533304 0.0000 313.512501
## TotalPctDiv 35.7865206 0.0000 -539.740385
## PersPerFam 12.8434139 0.0000 -892.724350
## PctFam2Par -63.3791097 0.0000 155.309904
## PctKids2Par -77.2658142 -274.6816 -1076.322008
## PctYoungKids2Par -48.7372012 0.0000 63.766829
## PctTeen2Par -52.3234082 0.0000 228.107231
## PctWorkMomYoungKids 15.4966945 0.0000 -139.933827
## PctWorkMom 15.5887797 0.0000 148.001458
## NumKidsBornNeverMar 7.3792925 0.0000 823.136293
## PctKidsBornNeverMar 90.1571272 111.9105 -275.716625
## NumImmig 8.8744033 0.0000 246.463966
## PctImmigRecent -14.4238048 0.0000 386.938247
## PctImmigRec5 -9.2828834 0.0000 -488.396949
## PctImmigRec8 7.8179446 0.0000 -283.578005
## PctImmigRec10 32.7165884 0.0000 553.681794
## PctRecentImmig -5.4882242 0.0000 -1499.140927
## PctRecImmig5 3.5802580 0.0000 2495.947793
## PctRecImmig8 7.7343690 0.0000 218.683864
## PctRecImmig10 15.4764031 0.0000 -2393.680669
## PctSpeakEnglOnly 20.7302810 0.0000 -120.364986
## PctNotSpeakEnglWell -3.5952710 0.0000 28.471431
## PctLargHouseFam 15.1526969 0.0000 1734.747671
## PctLargHouseOccup 2.2432740 0.0000 -2074.484590
## PersPerOccupHous 5.8992397 0.0000 213.898997
## PersPerOwnOccHous 21.7596068 0.0000 595.806448
## PersPerRentOccHous -2.9449743 0.0000 -661.765417
## PctPersOwnOccup 1.5507281 0.0000 -1620.311758
## PctPersDenseHous 30.0871764 0.0000 867.027964
## PctHousLess3BR 21.8513253 0.0000 66.384520
## MedNumBR -2.4237333 0.0000 -47.465555
## HousVacant 23.5587002 0.0000 460.698398
## PctHousOccup -22.2653743 0.0000 35.298779
## PctHousOwnOcc -0.1461844 0.0000 1686.420107
## PctVacantBoarded 53.7638070 0.0000 40.344261
## PctVacMore6Mos 13.4318226 0.0000 62.599055
## MedYrHousBuilt 34.6951265 0.0000 158.543121
## PctHousNoPhone 19.3481570 0.0000 53.051654
## PctWOFullPlumb 0.8501767 0.0000 -148.569984
## OwnOccLowQuart -10.5489552 0.0000 -360.390700
## OwnOccMedVal -6.2428537 0.0000 417.935802
## OwnOccHiQuart -4.8448381 0.0000 -391.086418
## OwnOccQrange 4.7102975 0.0000 NA
## RentLowQ -14.9873332 0.0000 -718.093140
## RentMedian 2.9430115 0.0000 260.558601
## RentHighQ -0.9732736 0.0000 -339.982879
## RentQrange 23.5623906 0.0000 NA
## MedRent 10.6792040 0.0000 1195.115997
## MedRentPctHousInc 12.5207380 0.0000 -39.531086
## MedOwnCostPctInc 22.1233026 0.0000 -94.471559
## MedOwnCostPctIncNoMtg -31.7979290 0.0000 6.990649
## NumInShelters 13.5471871 0.0000 -39.885696
## NumStreet 5.7852204 0.0000 26.024385
## PctForeignBorn 17.5176905 0.0000 935.263522
## PctBornSameState -24.2491275 0.0000 124.249332
## PctSameHouse85 19.7928596 0.0000 144.115665
## PctSameCity85 3.2963396 0.0000 -30.985328
## PctSameState85 6.1610943 0.0000 -37.785119
## LemasSwornFT 0.6214069 0.0000 -271.999351Nuestro puntaje subió un poco (lasso). Construyamos gráficos de coeficientes para ver cómo el valor de lambda influye en los coeficientes de ambos modelos.
Usaremos la función glmnet para entrenar los modelos y luego usaremos la función plot() que produce un diagrama de perfil de coeficientes de las rutas de los coeficientes para un objeto glmnet ajustado.
El parámetro xvar de plot() define lo que está en el eje X, y hay 3 valores posibles que puede tomar: norm, lambda o dev, donde norm es la norma \(L_1\) de los coeficientes. , lambda para la secuencia log-lambda y dev es el porcentaje de desviación explicado. Lo configuraremos en lambda.
Para entrenar glment, necesitamos convertir nuestro objeto X_train en una matriz.
# Fijamos los valores de los coeficientes
paramLasso <- seq(0, 500, 5)
paramRidge <- seq(0, 500, 5)
# Convertimos en matriz para glmnet
X_train_m <- as.matrix(X_train)
# Construimos los modelos para los valores de lambda values of lambda
rridge <- glmnet(
x = X_train_m,
y = y_train,
alpha = 0, #Ridge
lambda = paramRidge
)
llaso <- glmnet(
x = X_train_m,
y = y_train,
alpha = 1, #Lasso
lambda = paramLasso
)
plot(rridge,xvar = c("lambda"),label = TRUE)
Usamos CV para visualizar la ganancia en MSE:
## [1] 1312.253
## [1] 86.33751Revisamos el \(R^2\) y el \(\text{RMSE}\):
# Convertimos en matriz para glmnet
X_test_m <- as.matrix(X_test)
predictions_cv_ridge <- cv.out_ridge %>% predict(X_test_m,s =bestlam_ridge)
predictions_cv_lasso <- cv.out_lasso %>% predict(X_test_m,s =bestlam_lasso)
data.frame(
Ridge_R2 = R2(predictions_ridge, y_test),
Lasso_R2 = R2(predictions_lasso, y_test),
Linear_R2 = R2(predictions_lin, y_test),
CV_ridge_R2 = R2(predictions_cv_ridge, y_test),
CV_lasso_R2 = R2(predictions_cv_lasso, y_test)
)## Ridge_R2 Lasso_R2 Linear_R2 s1 s1.1
## 1 0.6737963 0.6189773 0.4857502 0.671287 0.6630695data.frame(
Ridge_RMSE = RMSE(predictions_ridge, y_test) ,
Lasso_RMSE = RMSE(predictions_lasso, y_test) ,
Linear_RMSE = RMSE(predictions_lin, y_test),
RMSE_ridge_R2 = RMSE(predictions_cv_ridge, y_test),
RMSE_lasso_R2 = RMSE(predictions_cv_lasso, y_test)
)## Ridge_RMSE Lasso_RMSE Linear_RMSE RMSE_ridge_R2 RMSE_lasso_R2
## 1 375.0716 431.1954 508.1641 372.903 384.7345