Introducción

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:

  1. estimar todos los modelos potencialmente interesantes;

  2. definir un criterio para seleccionar el mejor de entre todos los modelos ajustados en (i).

Criterios

Para cuantificar la “bondad” de un modelo de regresión múltiple podemos utilizar distintos criterios:

El coeficiente de determinación ajustado \(\bar R^2\)

\[ \bar R^2 = 1-\frac{RSS/(n-k-1)}{TSS/(n-1)} \]

El criterio de información de Akaike (Akaike Information Criterion \(\equiv\) AIC)

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. \]

El criterio de información bayesiana (Bayesian Information Criterion \(\equiv\) BIC)

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. \]

La \(C_p\) de Mallows

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

Todas las posibles regresiones

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

Métodos paso a paso (Stepwise Methods)

Eliminación progresiva (Backwards Elimination \(\equiv\) BE)

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

Criterio para eliminar un regresor:

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

Introducción progresiva (Forwards Selection \(\equiv\) FS)

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)

Regresión paso a paso (Stepwise Regression \(\equiv\) SR)

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

Método exhaustivo o análisis de los mejores subconjuntos (Exhaustive or All-Possible-Subsets Method \(\equiv\) APS)

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

LASSO (Least Absolute Shrinkage and Selection Operator)

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

Referencias

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.