Cuando el número \(k\) de variables explicativas en el modelo de regresión múltiple es demasiado alto (en el sentido de que se incluyen variables con baja capacidad de explicar la variable respuesta \(Y\)), se estiman demasiados parámetros \(\beta_j\), el modelo estimado se sobreajusta a los datos (con lo que disminuye el sesgo) y tiene una varianza excesiva. Por el contrario, si el modelo de regresión tiene un déficit de regresores, la varianza del modelo estimado será menor, pero tendrá mayor sesgo y su capacidad predictiva será menor (se infraajusta a los datos). Es deseable alcanzar un equilibrio entre estas dos situaciones extremas, de manera que permanezcan en el modelo final los regresores “importantes”, es decir, aquéllos cuya ausencia empobrece sensiblemente la capacidad predictiva del modelo.
Para seleccionar el modelo “óptimo” de regresión necesitamos:
estimar todos los modelos potencialmente interesantes;
definir un criterio para seleccionar el mejor de entre todos los modelos ajustados en (i).
Para cuantificar la “bondad” de un modelo de regresión múltiple podemos utilizar distintos criterios:
\[ \bar R^2 = 1-\frac{RSS/(n-k-1)}{TSS/(n-1)} \]
El criterio de información propuesto por Akaike (1973) está relacionado con el error de predicción, penalizado por la complejidad del modelo \[ \mbox{AIC} = n \, \log\left(\frac{RSS}{n}\right) +2k. \]
Propuesto por Akaike (1978) y Schwarz (1978), penaliza más la complejidad del modelo que el AIC \[ \mbox{BIC} = n \, \log\left(\frac{RSS}{n}\right) + k \log n. \]
Mallows (1973, 1995) definió el siguiente criterio para comparar un modelo más general con un submodelo \(P\) anidado (es decir, que es un caso particular del modelo general) con \(p\) parámetros \[ C_p = \frac{\mbox{RSS}_P}{s_R^2} - (n-2p) \] siendo \(s_R^2\) la varianza residual del modelo general y \(\mbox{RSS}_P\) la suma de cuadrados residual en el submodelo \(P\).
El número de posibles modelos es \(2^k\) (todos los subconjuntos de regresores que podemos formar). Si \(k\) es “grande”, es inviable y poco práctico ajustar todos estos modelos y cuantificar su capacidad predictiva. Por esta razón, se han propuesto estrategias que seleccionan una cantidad razonable de posibles modelos entre los que luego elegir con algún criterio de optimalidad. No hay un método de selección de variables que sea preferible por encima de los demás: depende de cómo midamos la “bondad” del modelo resultante.
El principio general que rige la selección de variables (y, en general, la selección de modelos) es el de la parsimonia, es decir, se busca un modelo lo más simple posible con la mayor capacidad predictiva posible. En este caso “simple” se refiere a un modelo que incorpore el menor número de regresores posible (lo cual además facilita su interpretabilidad).
Con un número pequeño de regresores sí podemos ajustar los \(2^k\) submodelos.
Ejemplo 1: Cargamos los datos fuel2001
sobre consumo de combustible (y otras variables relacionadas) en
EE.UU.
load("combustible.Rdata")
reg = lm(FuelC ~ Drivers+Income+Miles+MPC+Tax,data=fuel2001)
#install.packages("olsrr")
library(olsrr)
ols_step_all_possible(reg)
## Index N Predictors R-Square Adj. R-Square Mallow's Cp
## 1 1 1 Drivers 0.97038133 0.96977686 22.384750
## 3 2 1 Miles 0.52457680 0.51487429 1066.727069
## 2 3 1 Income 0.04439404 0.02489188 2191.603896
## 4 4 1 MPC 0.04203193 0.02248156 2197.137400
## 5 5 1 Tax 0.03995598 0.02036324 2202.000515
## 7 6 2 Drivers Miles 0.97816641 0.97725667 6.147408
## 8 7 2 Drivers MPC 0.97430639 0.97323582 15.189894
## 6 8 2 Drivers Income 0.97223852 0.97108179 20.034093
## 9 9 2 Drivers Tax 0.97178521 0.97060959 21.096008
## 10 10 2 Income Miles 0.62160831 0.60584199 841.420926
## 13 11 2 Miles MPC 0.58130414 0.56385848 935.837524
## 14 12 2 Miles Tax 0.54813561 0.52930793 1013.538175
## 15 13 2 MPC Tax 0.10024582 0.06275606 2062.765446
## 12 14 2 Income Tax 0.08345953 0.04527034 2102.089044
## 11 15 2 Income MPC 0.05553528 0.01618258 2167.504423
## 20 16 3 Drivers Miles Tax 0.98002394 0.97874888 3.795941
## 19 17 3 Drivers Miles MPC 0.97951538 0.97820785 4.987305
## 16 18 3 Drivers Income Miles 0.97820296 0.97681166 8.061780
## 21 19 3 Drivers MPC Tax 0.97484384 0.97323812 15.930865
## 17 20 3 Drivers Income MPC 0.97445140 0.97282064 16.850194
## 18 21 3 Drivers Income Tax 0.97353854 0.97184951 18.988649
## 23 22 3 Income Miles Tax 0.64333593 0.62057014 792.521772
## 22 23 3 Income Miles MPC 0.62769591 0.60393182 829.160110
## 25 24 3 Miles MPC Tax 0.62113034 0.59694717 844.540605
## 24 25 3 Income MPC Tax 0.10727903 0.05029684 2048.289440
## 29 26 4 Drivers Miles MPC Tax 0.98069740 0.97901891 4.218308
## 27 27 4 Drivers Income Miles Tax 0.98003409 0.97829793 5.772158
## 26 28 4 Drivers Income Miles MPC 0.97965783 0.97788895 6.653593
## 28 29 4 Drivers Income MPC Tax 0.97506019 0.97289151 17.424027
## 30 30 4 Income Miles MPC Tax 0.65654643 0.62668090 763.574837
## 31 31 5 Drivers Income Miles MPC Tax 0.98079059 0.97865621 6.000000
Podemos seleccionar el conjunto de regresores que optimiza alguno de los criterios.
Ejemplo 1:
ols_step_best_subset(reg)
## Best Subsets Regression
## -------------------------------------------
## Model Index Predictors
## -------------------------------------------
## 1 Drivers
## 2 Drivers Miles
## 3 Drivers Miles Tax
## 4 Drivers Miles MPC Tax
## 5 Drivers Income Miles MPC Tax
## -------------------------------------------
##
## Subsets Regression Summary
## -------------------------------------------------------------------------------------------------------------------------------------------------------------
## Adj. Pred
## Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
## -------------------------------------------------------------------------------------------------------------------------------------------------------------
## 1 0.9704 0.9698 0.9638 22.3847 1480.6292 1334.6972 1486.4247 1.121512e+13 228520715448.5286 4581193588.0012 0.0320
## 2 0.9782 0.9773 0.9668 6.1474 1467.0765 1322.3274 1474.8038 8.443198e+12 175209604443.8814 3520760372.9858 0.0246
## 3 0.9800 0.9787 0.9674 3.7959 1464.5418 1320.5126 1474.2009 7.892808e+12 166745750889.3132 3361277982.3537 0.0234
## 4 0.9807 0.9790 0.9674 4.2183 1464.7928 1321.3070 1476.3837 7.7962e+12 167620164347.6546 3392312849.8930 0.0235
## 5 0.9808 0.9787 0.9664 6.0000 1466.5460 1323.3787 1480.0688 7.934892e+12 173562786578.7492 3529386808.4195 0.0243
## -------------------------------------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria
## SBIC: Sawa's Bayesian Information Criteria
## SBC: Schwarz Bayesian Criteria
## MSEP: Estimated error of prediction, assuming multivariate normality
## FPE: Final Prediction Error
## HSP: Hocking's Sp
## APC: Amemiya Prediction Criteria
Se comienza con el conjunto completo de \(k\) variables explicativas. En cada paso del algoritmo BE se elimina el regresor menos explicativo individualmente. A continuación se ajusta el modelo con los \(k-1\) regresores restantes y se elimina un nuevo regresor, el menos significativo. Entonces se ajusta el modelo con \(k-2\) regresores y se elimina una variable explicativa más. El procedimiento se itera, eliminando progresivamente un regresor en cada iteración. El modelo reducido siempre tiene que ser un “submodelo” del modelo más general (es decir, tiene que resultar de eliminar regresores del modelo general).
Como en cada paso del algoritmo se plantea el contraste: \[ \begin{array}{ll} H_0: \mbox{Modelo reducido} \\ H_1: \mbox{Modelo de \(H_0\) con 1 regresor más (modelo ampliado)} \end{array} \] se puede elegir aquella variable explicativa cuyo estadístico del contraste \[ F = \frac{(\mbox{RSS}_0-\mbox{RSS}_1)/(\mbox{df}_0-\mbox{df}_1)}{\mbox{RSS}_1/\mbox{df}_1} \] sea el mínimo, siendo \(\mbox{RSS}_0\) la suma de cuadrados residual (con \(\mbox{df}_0\) grados de libertad) del modelo reducido y \(\mbox{RSS}_1\) la suma de cuadrados residual (con \(\mbox{df}_1\) g.l.) del modelo ampliado.
Como se eliminan las variables explicativas una a una, se cumple que \(\mbox{df}_0-\mbox{df}_1=1\) y \(\mbox{df}_1=n-k_1-1\), siendo \(k_1\) el número de regresores del modelo ampliado (en esa iteración). Dado que la distribución \(t\) de Student con \(\nu\) grados de libertad cumple que \(t_\nu^2=F_ {1;\nu}\), este procedimiento BE equivale a eliminar en cada iteración aquella \(X_j\) cuyo estadístico del contraste individual \(t(\beta_j)=\hat\beta_j/\mbox{s.e.}(\hat\beta_j)\) sea el menor de todos (en el modelo ampliado de esa iteración).
Se suele parar el procedimiento iterativo cuando todas las variables
del modelo tienen un ratio \(F\) mayor
que un valor predeterminado, por ejemplo, el cuantil 0.9 de la \(F_{1;n-k_1-1}\). En la librería
olsrr
se para el procedimiento de eliminar variables cuando
todos los p-valores de los contrastes individuales \(H_0:\beta_j=0\) son mayores que 0.3 (este
valor se puede modificar en el argumento penter
de los
comandos de olsrr
).
La eliminación progresiva tiene el inconveniente de ser computacionalmente costosa (si sólo un subconjunto pequeño de las \(k\) variables originales es significativo) y además, en los primeros pasos, necesita ajustar modelos con tal cantidad de regresores que pueden tener problemas de multicolinealidad.
Hay varias librerías de R
que han implementado la
selección paso a paso, por ejemplo, olsrr
,
Rcmdr
, leaps
,…
Ejemplo 1:
El algoritmo de eliminación progresiva resulta en eliminar únicamente
la variable Income
:
ols_step_backward_p(reg)
##
##
## Elimination Summary
## ------------------------------------------------------------------------------
## Variable Adj.
## Step Removed R-Square R-Square C(p) AIC RMSE
## ------------------------------------------------------------------------------
## 1 Income 0.9807 0.979 4.2183 1464.7928 390709.7110
## ------------------------------------------------------------------------------
Para ver detalles concretos sobre el algoritmo:
ols_step_backward_p(reg,details=TRUE)
## Backward Elimination Method
## ---------------------------
##
## Candidate Terms:
##
## 1 . Drivers
## 2 . Income
## 3 . Miles
## 4 . MPC
## 5 . Tax
##
## We are eliminating variables based on p value...
##
## - Income
##
## Backward Elimination: Step 1
##
## Variable Income Removed
##
## Model Summary
## ------------------------------------------------------------------------
## R 0.990 RMSE 390709.711
## R-Squared 0.981 Coef. Var 15.365
## Adj. R-Squared 0.979 MSE 152654078245.185
## Pred R-Squared 0.967 MAE 248064.575
## ------------------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------------------
## Regression 3.567676e+14 4 8.919189e+13 584.275 0.0000
## Residual 7.022088e+12 46 152654078245.185
## Total 3.637897e+14 50
## ------------------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------------------------
## (Intercept) -178848.203 474154.982 -0.377 0.708 -1133272.680 775576.274
## Drivers 0.618 0.021 0.914 29.273 0.000 0.575 0.660
## Miles 5.582 1.495 0.110 3.735 0.001 2.574 8.590
## MPC 39.030 30.809 0.030 1.267 0.212 -22.985 101.046
## Tax -21549.383 12839.610 -0.036 -1.678 0.100 -47394.177 4295.411
## -------------------------------------------------------------------------------------------------------
##
##
##
## No more variables satisfy the condition of p value = 0.3
##
##
## Variables Removed:
##
## - Income
##
##
## Final Model Output
## ------------------
##
## Model Summary
## ------------------------------------------------------------------------
## R 0.990 RMSE 390709.711
## R-Squared 0.981 Coef. Var 15.365
## Adj. R-Squared 0.979 MSE 152654078245.185
## Pred R-Squared 0.967 MAE 248064.575
## ------------------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------------------
## Regression 3.567676e+14 4 8.919189e+13 584.275 0.0000
## Residual 7.022088e+12 46 152654078245.185
## Total 3.637897e+14 50
## ------------------------------------------------------------------------------
##
## Parameter Estimates
## -------------------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -------------------------------------------------------------------------------------------------------
## (Intercept) -178848.203 474154.982 -0.377 0.708 -1133272.680 775576.274
## Drivers 0.618 0.021 0.914 29.273 0.000 0.575 0.660
## Miles 5.582 1.495 0.110 3.735 0.001 2.574 8.590
## MPC 39.030 30.809 0.030 1.267 0.212 -22.985 101.046
## Tax -21549.383 12839.610 -0.036 -1.678 0.100 -47394.177 4295.411
## -------------------------------------------------------------------------------------------------------
##
##
## Elimination Summary
## ------------------------------------------------------------------------------
## Variable Adj.
## Step Removed R-Square R-Square C(p) AIC RMSE
## ------------------------------------------------------------------------------
## 1 Income 0.9807 0.979 4.2183 1464.7928 390709.7110
## ------------------------------------------------------------------------------
Se comienza con un subconjunto de 0 regresores. En cada iteración se
añade la variable \(X_j\) con el mayor
ratio \(F\) (Izenman 2008). Esto es
equivalente a introducir la variable que tiene menor p-valor en los
contrastes individuales \(H_0:\beta_j=0\) del modelo ampliado.
Alternativamente, se puede añadir la que tiene mayor coeficiente de
correlación parcial con la respuesta \(Y\) (ver definición de correlación parcial
en Peña 2002). Se para el proceso de introducir variables en el modelo
cuando el valor de \(F\) de cada
variable que no está incluida en ese momento en el modelo es menor que
un nivel predeterminado, por ejemplo, el cuantil 0.25 de la distribución
\(F_{1;n-k_1-1}\), o cuando los
p-valores de los contrastes individuales en los modelos ampliados
superen una cota (0.3 en la librería olsrr
por
defecto).
Ejemplo 1: En el primer paso del algoritmo ajustamos
todos los modelos de regresión simple posibles y elegimos la variable
Drivers
porque su p-valor es el menor e inferior a 0.3.
reg1Drivers <- lm(FuelC ~ Drivers,data=fuel2001)
summary(reg1Drivers)$coefficients[,"Pr(>|t|)"]
## (Intercept) Drivers
## 6.145706e-01 4.113896e-39
reg1Income <- lm(FuelC ~ Income,data=fuel2001)
summary(reg1Income)$coefficients[,"Pr(>|t|)"]
## (Intercept) Income
## 0.6579428 0.1377809
reg1Miles <- lm(FuelC ~ Miles,data=fuel2001)
summary(reg1Miles)$coefficients[,"Pr(>|t|)"]
## (Intercept) Miles
## 5.091944e-01 1.887021e-09
reg1MPC <- lm(FuelC ~ MPC,data=fuel2001)
summary(reg1MPC)$coefficients[,"Pr(>|t|)"]
## (Intercept) MPC
## 0.00873296 0.14896499
reg1Tax <- lm(FuelC ~ Tax,data=fuel2001)
summary(reg1Tax)$coefficients[,"Pr(>|t|)"]
## (Intercept) Tax
## 0.005948553 0.159620556
A continuación ajustaríamos todos los modelos con 2 regresores (uno
de ellos Drivers
necesariamente). Vemos que el segundo
regresor seleccionado sería Miles
:
reg2DriversIncome <- lm(FuelC ~ Drivers+Income,data=fuel2001)
summary(reg2DriversIncome)$coefficients[,"Pr(>|t|)"]
## (Intercept) Drivers Income
## 6.890148e-02 1.517511e-38 7.944451e-02
reg2DriversMiles <- lm(FuelC ~ Drivers+Miles,data=fuel2001)
summary(reg2DriversMiles)$coefficients[,"Pr(>|t|)"]
## (Intercept) Drivers Miles
## 3.478349e-02 9.076651e-34 1.409020e-04
reg2DriversMPC <- lm(FuelC ~ Drivers+MPC,data=fuel2001)
summary(reg2DriversMPC)$coefficients[,"Pr(>|t|)"]
## (Intercept) Drivers MPC
## 1.579281e-02 2.228956e-39 9.354085e-03
reg2DriversTax <- lm(FuelC ~ Drivers+Tax,data=fuel2001)
summary(reg2DriversTax)$coefficients[,"Pr(>|t|)"]
## (Intercept) Drivers Tax
## 1.101579e-01 2.003220e-38 1.288118e-01
Ejecutamos el algoritmo FS de golpe con la librería
olsrr
:
ols_step_forward_p(reg)
##
## Selection Summary
## -------------------------------------------------------------------------------
## Variable Adj.
## Step Entered R-Square R-Square C(p) AIC RMSE
## -------------------------------------------------------------------------------
## 1 Drivers 0.9704 0.9698 22.3847 1480.6292 468932.0763
## 2 Miles 0.9782 0.9773 6.1474 1467.0765 406787.0912
## 3 Tax 0.9800 0.9787 3.7959 1464.5418 393215.9549
## 4 MPC 0.9807 0.9790 4.2183 1464.7928 390709.7110
## -------------------------------------------------------------------------------
Si queremos resultados muy detallados de cualquiera de los
procedimientos paso a paso cambiamos el argumento details
a
TRUE
.
Ejemplo 1:
ols_step_forward_p(reg, details = TRUE)
Se trata de un procedimiento híbrido entre introducción y eliminación progresivas. Como en la introducción progresiva de predictores, se comienza con 0 variables. En cada paso, al incluir una nueva variable y estimar el nuevo modelo, se evalúa la pertinencia de eliminar cualquiera de los regresores ya presentes (por ejemplo, mediante un contraste individual con la \(t\) de Student). Se puede rechazar, por tanto, alguna de las \(X_j\) ya incluidas en el modelo.
Ejemplo 1:
ols_step_both_p(reg)
##
## Stepwise Selection Summary
## -------------------------------------------------------------------------------------------
## Added/ Adj.
## Step Variable Removed R-Square R-Square C(p) AIC RMSE
## -------------------------------------------------------------------------------------------
## 1 Drivers addition 0.970 0.970 22.3850 1480.6292 468932.0763
## 2 Miles addition 0.978 0.977 6.1470 1467.0765 406787.0912
## 3 Tax addition 0.980 0.979 3.7960 1464.5418 393215.9549
## -------------------------------------------------------------------------------------------
Consiste en examinar todos los posibles subconjuntos de regresores con un cardinal \(r\leq k\) predeterminado. Esto significa que hay que examinar \(2^r-1\) subconjuntos de regresores (ignoramos el conjunto vacío). Dado que este número es enorme incluso para \(r\)’s moderados, se utilizan algoritmos de optimización (como el branch-and-bound) que reducen drásticamente el número de modelos a examinar.
Sea \(p\in\{1,2,\ldots,r+1\}\) el número de parámetros incluidos en cada submodelo \(P\) de regresión lineal múltiple (con \(p-1\) regresores y un término independiente). Los \(r \choose k\) posibles submodelos se pueden comparar entre sí mediante un criterio de selección de modelos. El criterio suele consistir en minimizar \(\mbox{RSS}_P\), la suma de cuadrados residuales del submodelo \(P\), penalizada (para que no haya sobreajuste): \[ \frac{\mbox{RSS}_P}{n} + \lambda \cdot p \cdot \frac{s_R^2}{n}, \] siendo \(\lambda>0\) un coeficiente o parámetro de penalización y \(s_R^2\) la varianza residual en el modelo ajustado con todos los \(r\) regresores. En ciertos contextos al término \(\mbox{RSS}_P/n\) se le denomina error de aprendizaje o de entrenamiento (learning error, training error, apparent error rate, resubstitution error, …). A \(\lambda \, p \, s_R^2/n\) se le llama término de complejidad, pues aumenta a medida que añadimos predictores al modelo.
Hay diferentes versiones o variaciones de la suma de cuadrados penalizada: por ejemplo, el criterio AIC, el criterio BIC o la \(C_p\) de Mallows.
Ejemplo 1: Tomamos un número máximo de regresores \(r=3\).
library(leaps)
APS <- regsubsets(FuelC~., data = fuel2001,
nbest = 1, # un único modelo óptimo para cada número de predictores
nvmax = 3, # nvmax = NULL para que no haya límite en el número de regresores
force.in = NULL, # regresores que deberían estar en todos los modelos
force.out = "Pop", # variables que no deberían estar en ningún modelo
method = "exhaustive")
summary(APS)
## Subset selection object
## Call: regsubsets.formula(FuelC ~ ., data = fuel2001, nbest = 1, nvmax = 3,
## force.in = NULL, force.out = "Pop", method = "exhaustive")
## 6 Variables (and intercept)
## Forced in Forced out
## Drivers FALSE FALSE
## Income FALSE FALSE
## Miles FALSE FALSE
## MPC FALSE FALSE
## Tax FALSE TRUE
## Pop FALSE FALSE
## 1 subsets of each size up to 3
## Selection Algorithm: exhaustive
## Drivers Income Miles MPC Tax Pop
## 1 ( 1 ) "*" " " " " " " " " " "
## 2 ( 1 ) "*" " " "*" " " " " " "
## 3 ( 1 ) "*" " " "*" " " "*" " "
Se puede hacer una tabla gráfica con los resultados:
plot(APS,scale="adjr2",main="R^2 ajustado")
plot(APS,main="Criterio BIC")
Se trata de una técnica que permite hacer selección de variables (entre otras cosas). Tibshirani (1996) la propuso originalmente en el contexto de la regresión lineal, pero se puede extender a problemas estadísticos diferentes. El objetivo del LASSO es resolver el problema de optimización \[ \min_{\boldsymbol\beta\in\mathbb R^k} \frac{\mbox{RSS}}{n} \mbox{ sujeto a } \sum_{j=1}^k \beta_j \leq t, \] para una cota \(t\) prefijada. Si \(k=2\), entonces \(\mbox{RSS} = \sum_{i=1}^n e_i^2 = \sum_{i=1}^n (y_i-\beta_0 -\beta_1 x_{1i}-\beta_2 x_{2i})^2\) es un polinomio de grado 2 en términos de \(\beta_1\) y \(\beta_2\) y la condición \(\beta_1+\beta_2\leq t\) delimita un rombo en el plano de \(\beta_1\) y \(\beta_2\). Ver el gráfico en el tutorial de r-bloggers. Dada la forma del conjunto \(\{ \boldsymbol\beta\in\mathbb R^k: \sum_{j=1}^k \beta_j \leq t \}\) donde se optimiza, la solución \(\hat\beta\) al anterior problema puede estar en uno de los vértices de esa región, es decir, varias componentes de \(\hat\beta\) pueden ser 0, lo cual conduce a una selección de variables.
El anterior problema de optimización se puede expresar de manera equivalente en forma lagrangiana para un parámetro de penalización \(\lambda\): \[ \min_{\boldsymbol\beta\in\mathbb R^k} \left( \frac{\mbox{RSS}}{n} +\lambda \|\boldsymbol\beta\|_1 \right). \] La cota \(t\) y el parámetro \(\lambda\) son uno función (no necesariamente sencilla) del otro. Pero no es necesario conocer esta función. Simplemente basta resolver el problema en su versión penalizada escogiendo adecuadamente el parámetro \(\lambda\).
library(glmnet)
# Matriz de regresores
X = as.matrix(fuel2001[,-c(2,6)])
# Selección de lambda por validación cruzada
cvfit = cv.glmnet(X,fuel2001$FuelC,alpha=1)
# Estimación de los coeficientes y selección de variables (coeficientes 0)
coef(cvfit, s = "lambda.1se")
## 6 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 3.677408e+04
## Drivers 5.853335e-01
## Income .
## Miles 4.013455e+00
## MPC .
## Tax .
El LASSO da lugar a un problema de optimización menos complejo que la regresión ridge que resuelve el problema \[ \min_{\boldsymbol\beta\in\mathbb R^k} \left( \frac{\mbox{RSS}}{n} +\lambda \|\boldsymbol\beta\|_2 \right). \] Se puede ampliar la información sobre este tipo de técnicas con el libro de Hastie et al. (2009).
Akaike, H. (1973). Information theory and an extension of the maximum likelihood principle. In B. N. Petrov and B. F. Csaki (Eds.), Second International Symposium on Information Theory, (pp. 267–281). Academiai Kiado: Budapest.
Akaike, H. (1978). On newer statistical approaches to parameter estimation and structure determination. International Federation of Automatic Control, 3, 1877–1884.
Hastie, T., Tibshirani, R. and Friedman (2009). The Elements of Statistical Learning. Data Mining, Inference, and Prediction. Second Edition. Springer.
Izenman, A.J. (2008). Modern Multivariate Statistical Techniques. Regression, Classification, and Manifold Learning. Springer.
Mallows, C.L. (1973). Some comments on \(C_P\). Technometrics, 15, 661-675.
Mallows, C.L. (1995). More comments on \(C_P\). Technometrics, 37, 362-372.
Peña, D. (2002). Regresión y diseño de experimentos. Alianza.
Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6, 461–464.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B, 58, 267–88.