Selección de Variables en modelos lineales
Javier Cara
Curso 2018-19
1 Introducción
library(ISLR)
datos = Hitters
str(datos)## 'data.frame': 322 obs. of 20 variables:
## $ AtBat : int 293 315 479 496 321 594 185 298 323 401 ...
## $ Hits : int 66 81 130 141 87 169 37 73 81 92 ...
## $ HmRun : int 1 7 18 20 10 4 1 0 6 17 ...
## $ Runs : int 30 24 66 65 39 74 23 24 26 49 ...
## $ RBI : int 29 38 72 78 42 51 8 24 32 66 ...
## $ Walks : int 14 39 76 37 30 35 21 7 8 65 ...
## $ Years : int 1 14 3 11 2 11 2 3 2 13 ...
## $ CAtBat : int 293 3449 1624 5628 396 4408 214 509 341 5206 ...
## $ CHits : int 66 835 457 1575 101 1133 42 108 86 1332 ...
## $ CHmRun : int 1 69 63 225 12 19 1 0 6 253 ...
## $ CRuns : int 30 321 224 828 48 501 30 41 32 784 ...
## $ CRBI : int 29 414 266 838 46 336 9 37 34 890 ...
## $ CWalks : int 14 375 263 354 33 194 24 12 8 866 ...
## $ League : Factor w/ 2 levels "A","N": 1 2 1 2 2 1 2 1 2 1 ...
## $ Division : Factor w/ 2 levels "E","W": 1 2 2 1 1 2 1 2 2 1 ...
## $ PutOuts : int 446 632 880 200 805 282 76 121 143 0 ...
## $ Assists : int 33 43 82 11 40 421 127 283 290 0 ...
## $ Errors : int 20 10 14 3 4 25 7 9 19 0 ...
## $ Salary : num NA 475 480 500 91.5 750 70 100 75 1100 ...
## $ NewLeague: Factor w/ 2 levels "A","N": 1 2 1 2 2 1 1 1 2 1 ...Comprobamos si hay missing obsevations (NA) en el salario:
sum(is.na(datos$Salary))## [1] 59Eliminamos estos datos:
datos = na.omit(datos)2 Variables no significativas
m1 = lm(Salary ~ ., data = datos)
summary(m1)##
## Call:
## lm(formula = Salary ~ ., data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -907.62 -178.35 -31.11 139.09 1877.04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 163.10359 90.77854 1.797 0.073622 .
## AtBat -1.97987 0.63398 -3.123 0.002008 **
## Hits 7.50077 2.37753 3.155 0.001808 **
## HmRun 4.33088 6.20145 0.698 0.485616
## Runs -2.37621 2.98076 -0.797 0.426122
## RBI -1.04496 2.60088 -0.402 0.688204
## Walks 6.23129 1.82850 3.408 0.000766 ***
## Years -3.48905 12.41219 -0.281 0.778874
## CAtBat -0.17134 0.13524 -1.267 0.206380
## CHits 0.13399 0.67455 0.199 0.842713
## CHmRun -0.17286 1.61724 -0.107 0.914967
## CRuns 1.45430 0.75046 1.938 0.053795 .
## CRBI 0.80771 0.69262 1.166 0.244691
## CWalks -0.81157 0.32808 -2.474 0.014057 *
## LeagueN 62.59942 79.26140 0.790 0.430424
## DivisionW -116.84925 40.36695 -2.895 0.004141 **
## PutOuts 0.28189 0.07744 3.640 0.000333 ***
## Assists 0.37107 0.22120 1.678 0.094723 .
## Errors -3.36076 4.39163 -0.765 0.444857
## NewLeagueN -24.76233 79.00263 -0.313 0.754218
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 315.6 on 243 degrees of freedom
## Multiple R-squared: 0.5461, Adjusted R-squared: 0.5106
## F-statistic: 15.39 on 19 and 243 DF, p-value: < 2.2e-16Este análisis nos qué variables son significativas en el análisis de Salary en función del resto de variables. Pero quedarnos con el modelo que contiene las variables significativas no nos garantiza que sea el mejor modelo para predecir.
3 Best subset selection
Algoritmo:
Para k = 1, 2, …, p:
- Estimar todos los modelos de k regresores (hay \(\binom{p}{k}\) modelos posibles).
- Elegir el que tenga menor RSS o mayor R\(^2\). Este será el modelo M\(_{k}\).
library(leaps)
m2 = regsubsets(Salary ~ ., data = datos)
summary(m2)## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = datos)
## 19 Variables (and intercept)
## Forced in Forced out
## AtBat FALSE FALSE
## Hits FALSE FALSE
## HmRun FALSE FALSE
## Runs FALSE FALSE
## RBI FALSE FALSE
## Walks FALSE FALSE
## Years FALSE FALSE
## CAtBat FALSE FALSE
## CHits FALSE FALSE
## CHmRun FALSE FALSE
## CRuns FALSE FALSE
## CRBI FALSE FALSE
## CWalks FALSE FALSE
## LeagueN FALSE FALSE
## DivisionW FALSE FALSE
## PutOuts FALSE FALSE
## Assists FALSE FALSE
## Errors FALSE FALSE
## NewLeagueN FALSE FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: exhaustive
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 5 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " " "
## 6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 7 ( 1 ) " " "*" " " " " " " "*" " " "*" "*" "*" " "
## 8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " "*" "*"
## CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1 ( 1 ) "*" " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " "*" " " " " " "
## 4 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 7 ( 1 ) " " " " " " "*" "*" " " " " " "
## 8 ( 1 ) " " "*" " " "*" "*" " " " " " "El resultado son los mejores 8 modelos (por defecto):
- la primera línea es el mejor modelo (en términos de R\(^2\)) de una variable. La variable seleccionada es la que aparece con un asterisco, CRBI.
- la segunda linea es el mejor modelo (en términos de R\(^2\)) de dos variables, Hits y CRBI.
- y así sucesivamente.
Podemos seleccionar el numero de modelos que nos devuelve con nvmax:
m_best = regsubsets(Salary ~ ., data = datos, nvmax = 19)
summary(m_best)## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = datos, nvmax = 19)
## 19 Variables (and intercept)
## Forced in Forced out
## AtBat FALSE FALSE
## Hits FALSE FALSE
## HmRun FALSE FALSE
## Runs FALSE FALSE
## RBI FALSE FALSE
## Walks FALSE FALSE
## Years FALSE FALSE
## CAtBat FALSE FALSE
## CHits FALSE FALSE
## CHmRun FALSE FALSE
## CRuns FALSE FALSE
## CRBI FALSE FALSE
## CWalks FALSE FALSE
## LeagueN FALSE FALSE
## DivisionW FALSE FALSE
## PutOuts FALSE FALSE
## Assists FALSE FALSE
## Errors FALSE FALSE
## NewLeagueN FALSE FALSE
## 1 subsets of each size up to 19
## Selection Algorithm: exhaustive
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 5 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " " "
## 6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 7 ( 1 ) " " "*" " " " " " " "*" " " "*" "*" "*" " "
## 8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " "*" "*"
## 9 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 10 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 11 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 12 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 13 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 14 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" " " " " "*"
## 15 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" "*" " " "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" " " "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"
## CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1 ( 1 ) "*" " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " "*" " " " " " "
## 4 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 7 ( 1 ) " " " " " " "*" "*" " " " " " "
## 8 ( 1 ) " " "*" " " "*" "*" " " " " " "
## 9 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 10 ( 1 ) "*" "*" " " "*" "*" "*" " " " "
## 11 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 14 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"Dentro de cada nivel (es decir, manteniendo fijo el número de variables), se elige la combinación de variables que obtiene el mayor R\(^2\). Pero nosotros queremos el modelo que obtiene el menor test error. Podemos calcular el test error de dos maneras:
- Diréctamente, utilizando subconjuntos de validation o cross-validation.
- Indirectamente, mediante las métricas Cp, AIC, BIC or R\(^2\) ajustado.
\[ Cp = \frac{1}{n}(RSS + 2d\hat{\sigma}^2) \]
\[ AIC = \frac{1}{n\hat{\sigma}^2}(RSS + 2d\hat{\sigma}^2) \]
Cp y AIC son proporcionales, \(C_p = AIC * \hat{\sigma}^2\).
\[ BIC = \frac{1}{n}(RSS + \log(n)d\hat{\sigma}^2) \]
\[ R^2-ajustado = 1 - \frac{RSS/(n-d-1)}{TSS/(n-1)} \]
donde:
- d: número de regresores.
- \(\hat{\sigma}^2\): estimación del error del modelo. En regresión simple es \(RSS/(n-2)\).
- RSS: Residual sum of squares
\[ RSS = \sum _{i=1}^n(y_i - \hat{y}_i)^2 \]
- TSS: Total Sum of Squares
\[ TSS = \sum _{i=1}^n(y_i - \bar{y}_i)^2 \]
En R, estas métricas están se pueden obtener de la siguiente manera:
m_best_summary = summary(m_best)
names(m_best_summary)## [1] "which" "rsq" "rss" "adjr2" "cp" "bic" "outmat" "obj"plot(m_best_summary$adjr2, xlab = "Numero de variables", ylab = "R2 ajustado", type = "b")
Buscamos el máximo:
which.max(m_best_summary$adjr2)## [1] 11Si utilizamos el criterio del Cp (que es equivalente al AIC):
plot(m_best_summary$cp, xlab = "Numero de variables", ylab = "R2 ajustado", type = "b")
which.min(m_best_summary$cp)## [1] 10Justificación: el RSS disminuye con el número de regresores d. Por eso penalizamos incluyendo un término que contiene a d.
Los coeficinetes estimados con ese modelo son:
coef(m_best,10)## (Intercept) AtBat Hits Walks CAtBat
## 162.5354420 -2.1686501 6.9180175 5.7732246 -0.1300798
## CRuns CRBI CWalks DivisionW PutOuts
## 1.4082490 0.7743122 -0.8308264 -112.3800575 0.2973726
## Assists
## 0.28316804 Método Forward-Stepwise
Algoritmo:
- M0 es el modelo sin regresores.
- Para k = 0, …, (p-1)
- A partir del modelo con k regresores, M\(_k\), estimar todos los modelos posibles con (k+1) regresores.
- Elegir el que tenga menor RSS o mayor R\(^2\). Este será el modelo M\(_{k+1}\).
- Elegir el mejor modelo de M0, …, M\(_{p}\) utilizando validación cruzada, Cp, AIC, BIC, R2-ajustado.
m_fwd = regsubsets(Salary ~ ., data = datos, nvmax = 19, method = "forward")
summary(m_fwd)## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = datos, nvmax = 19, method = "forward")
## 19 Variables (and intercept)
## Forced in Forced out
## AtBat FALSE FALSE
## Hits FALSE FALSE
## HmRun FALSE FALSE
## Runs FALSE FALSE
## RBI FALSE FALSE
## Walks FALSE FALSE
## Years FALSE FALSE
## CAtBat FALSE FALSE
## CHits FALSE FALSE
## CHmRun FALSE FALSE
## CRuns FALSE FALSE
## CRBI FALSE FALSE
## CWalks FALSE FALSE
## LeagueN FALSE FALSE
## DivisionW FALSE FALSE
## PutOuts FALSE FALSE
## Assists FALSE FALSE
## Errors FALSE FALSE
## NewLeagueN FALSE FALSE
## 1 subsets of each size up to 19
## Selection Algorithm: forward
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 5 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " " "
## 6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 7 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " "*"
## 9 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 10 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 11 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 12 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 13 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 14 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" " " " " "*"
## 15 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" "*" " " "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" " " "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"
## CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1 ( 1 ) "*" " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " "*" " " " " " "
## 4 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 7 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 8 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 9 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 10 ( 1 ) "*" "*" " " "*" "*" "*" " " " "
## 11 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 14 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"m_fwd_summary = summary(m_fwd)
which.min(m_fwd_summary$cp)## [1] 10coef(m_fwd,10)## (Intercept) AtBat Hits Walks CAtBat
## 162.5354420 -2.1686501 6.9180175 5.7732246 -0.1300798
## CRuns CRBI CWalks DivisionW PutOuts
## 1.4082490 0.7743122 -0.8308264 -112.3800575 0.2973726
## Assists
## 0.28316805 Método Backward-Stepwise
Algoritmo:
- Mp es el modelo con todos los regresores.
- Para k = p, …, 1
- A partir del modelo con k regresores, M\(_k\), estimar todos los modelos posibles con (k-1) regresores.
- Elegir el que tenga menor RSS o mayor R\(^2\). Este será el modelo M\(_{k-1}\).
- Elegir el mejor modelo de M0, …, M\(_{p}\) utilizando validación cruzada, Cp, AIC, BIC, R2 ajustado.
m_bwd = regsubsets(Salary ~ ., data = datos, nvmax = 19, method = "backward")
summary(m_bwd)## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = datos, nvmax = 19, method = "backward")
## 19 Variables (and intercept)
## Forced in Forced out
## AtBat FALSE FALSE
## Hits FALSE FALSE
## HmRun FALSE FALSE
## Runs FALSE FALSE
## RBI FALSE FALSE
## Walks FALSE FALSE
## Years FALSE FALSE
## CAtBat FALSE FALSE
## CHits FALSE FALSE
## CHmRun FALSE FALSE
## CRuns FALSE FALSE
## CRBI FALSE FALSE
## CWalks FALSE FALSE
## LeagueN FALSE FALSE
## DivisionW FALSE FALSE
## PutOuts FALSE FALSE
## Assists FALSE FALSE
## Errors FALSE FALSE
## NewLeagueN FALSE FALSE
## 1 subsets of each size up to 19
## Selection Algorithm: backward
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " "*"
## 2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " "*"
## 3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " "*"
## 4 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " "*"
## 5 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " "*"
## 6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " "*"
## 7 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " "*"
## 8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " "*"
## 9 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 10 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 11 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 12 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 13 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 14 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" " " " " "*"
## 15 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" "*" " " "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" " " "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"
## CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1 ( 1 ) " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " "*" " " " " " "
## 4 ( 1 ) " " " " " " " " "*" " " " " " "
## 5 ( 1 ) " " " " " " " " "*" " " " " " "
## 6 ( 1 ) " " " " " " "*" "*" " " " " " "
## 7 ( 1 ) " " "*" " " "*" "*" " " " " " "
## 8 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 9 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 10 ( 1 ) "*" "*" " " "*" "*" "*" " " " "
## 11 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 14 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"m_bwd_summary = summary(m_bwd)
which.min(m_bwd_summary$cp)## [1] 10coef(m_bwd,10)## (Intercept) AtBat Hits Walks CAtBat
## 162.5354420 -2.1686501 6.9180175 5.7732246 -0.1300798
## CRuns CRBI CWalks DivisionW PutOuts
## 1.4082490 0.7743122 -0.8308264 -112.3800575 0.2973726
## Assists
## 0.2831680Es el mismo que antes.
6 Eligiendo el mejor modelo utilizando subconjuntos de validación
set.seed(1)
n = nrow(datos)
n_train = round(n/2)
n_test = n - n_train
v = 1:n
pos_train = sample(v,n_train,replace = F) # muestreo sin reemplazamiento
pos_test = v[-pos_train]
# dividimos los datos en training set y validation set
datos_train = datos[pos_train,]
datos_test = datos[pos_test,]Estimamos todos los modelos posibles con los datos de entrenamiento
m_val = regsubsets(Salary ~ ., data = datos_train, nvmax = 19)Vamos a predecir con los datos de test. Como no existe la funcion predecir para este tipo de modelos, la construimos nosotros:
X_test = model.matrix(Salary ~ ., data = datos_test)source("MSE.R")
#
mse_val = rep(0,19)
for (i in 1:19){
coef_i = coef(m_val,i)
pred_i = X_test[,names(coef_i)] %*% coef_i
mse_val[i] = MSE(datos_test$Salary,pred_i)
}
which.min(mse_val)## [1] 8Para calcular los coeficientes del modelo de regresión finales, es preferible hacerlo con todos los datos:
m_val_final = regsubsets(Salary ~ ., data = datos, nvmax = 19)
coef(m_val_final,8)## (Intercept) AtBat Hits Walks CHmRun
## 130.9691577 -2.1731903 7.3582935 6.0037597 1.2339718
## CRuns CWalks DivisionW PutOuts
## 0.9651349 -0.8323788 -117.9657795 0.2908431La función que predice a partir de modelos estimados con regsubsets() la vamos a utilizar más veces. Por eso definimos una función para predecir:
source("regsubsets_predict.R")
predict.regsubsets## function (objeto_regsubsets, newdata, id)
## {
## formu = as.formula(objeto_regsubsets$call[[2]])
## X_mat = model.matrix(formu, newdata)
## coefi = coef(objeto_regsubsets, id)
## X_col = names(coefi)
## pred = X_mat[, X_col] %*% coefi
## return(pred)
## }Comprobamos que funciona bien:
head(predict.regsubsets(m_val_final,datos,8))## [,1]
## -Alan Ashby 425.2400
## -Alvin Davis 715.8659
## -Andre Dawson 1153.0071
## -Andres Galarraga 521.4537
## -Alfredo Griffin 603.3258
## -Al Newman 159.5798Tambien funciona con predict()
head(predict(m_val_final,datos,8))## [,1]
## -Alan Ashby 425.2400
## -Alvin Davis 715.8659
## -Andre Dawson 1153.0071
## -Andres Galarraga 521.4537
## -Alfredo Griffin 603.3258
## -Al Newman 159.57987 Eligiendo el mejor modelo utilizando Cross-Validation
Función para obtener las posiciones de train y de test:
source("cross_val_pos.R")Datos de los folds:
num_folds = 10
set.seed(1)
pos = cross_val_pos(nrow(datos),num_folds)Calculamos el error cometido en cada fold por cada modelo:
mse_cv = matrix(0, nrow = num_folds, ncol = 19)
for (i in 1:num_folds){
# datos de training y de validation de cada fold
datos_train = datos[pos$train[[i]],]
datos_test = datos[pos$test[[i]],]
m_cv = regsubsets(Salary ~ .,data = datos_train, nvmax = 19)
for (j in 1:19){
pred = predict(m_cv,newdata = datos_test, id = j)
mse_cv[i,j] = MSE(datos_test$Salary,pred)
}
}mse_cv_med = apply(mse_cv, 2, mean)
plot(mse_cv_med, type = "b")
El que tiene menor error es el de 11 variables. Lo aplicamos a todos los datos:
m_cv_final = regsubsets(Salary ~ ., data = datos, nvmax = 19)
coef(m_cv_final,11)## (Intercept) AtBat Hits Walks CAtBat
## 135.7512195 -2.1277482 6.9236994 5.6202755 -0.1389914
## CRuns CRBI CWalks LeagueN DivisionW
## 1.4553310 0.7852528 -0.8228559 43.1116152 -111.1460252
## PutOuts Assists
## 0.2894087 0.2688277