Objetivo

En este ejercicio se busca encontrar un modelo de regresion lineal que explique la variable respuesta \(y\) en funcion de las covariables presentadas en la tabla b3 del texto.

Base de datos

A continuacion se muestra el encabezado de la base de datos y la definicion de las variables.

Nota: Type of transmission (1=automatic, 0=manual).

Antes de iniciar es necesario revisar si hay NA's y eliminarlos.

require(MPV)
table.b3[22:26, ] # Can you see the missing values?
##        y    x1  x2  x3  x4   x5 x6 x7    x8   x9  x10 x11
## 22 21.47 360.0 180 290 8.4 2.45  2  3 214.2 76.3 4250   1
## 23 16.59 400.0 185  NA 7.6 3.08  4  3 196.0 73.0 3850   1
## 24 31.90  96.9  75  83 9.0 4.30  2  5 165.2 61.8 2275   0
## 25 29.40 140.0  86  NA 8.0 2.92  2  4 176.4 65.4 2150   0
## 26 13.27 460.0 223 366 8.0 3.00  4  3 228.0 79.8 5430   1
datis <- table.b3[-c(23,25), ]

El objeto datis tiene la base de datos sin las lineas con NA, lo mismo se hubiese podido realizar usando la funcion na.omit.

Diagramas de dispersion

A continuacion se muestran los diagramas de dispersion para las variables de la base de datos.

Aplicacion del metodo backward

Vamos a crear un modelo saturado, es decir, el modelo mayor a considerar.

full.model <- lm(y ~ ., data = datis)
summary(full.model)
## 
## Call:
## lm(formula = y ~ ., data = datis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3441 -1.6711 -0.4486  1.4906  5.2508 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 17.339838  30.355375   0.571   0.5749  
## x1          -0.075588   0.056347  -1.341   0.1964  
## x2          -0.069163   0.087791  -0.788   0.4411  
## x3           0.115117   0.088113   1.306   0.2078  
## x4           1.494737   3.101464   0.482   0.6357  
## x5           5.843495   3.148438   1.856   0.0799 .
## x6           0.317583   1.288967   0.246   0.8082  
## x7          -3.205390   3.109185  -1.031   0.3162  
## x8           0.180811   0.130301   1.388   0.1822  
## x9          -0.397945   0.323456  -1.230   0.2344  
## x10         -0.005115   0.005896  -0.868   0.3971  
## x11          0.638483   3.021680   0.211   0.8350  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.227 on 18 degrees of freedom
## Multiple R-squared:  0.8355, Adjusted R-squared:  0.7349 
## F-statistic:  8.31 on 11 and 18 DF,  p-value: 5.231e-05

De la tabla anterior se puede pensar en que hay un efecto de enmascaramiento entre las variables ya que ninguna parece significativa marginalmente.

Se usa la funcion stepAIC y se elije trace=TRUE para obtener detalles del proceso de seleccion.

require(MASS)  # Para poder usar la funcion stepAIC
modback <- stepAIC(full.model, trace=TRUE, direction="backward")
## Start:  AIC=78.96
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11
## 
##        Df Sum of Sq    RSS    AIC
## - x11   1     0.465 187.87 77.036
## - x6    1     0.632 188.03 77.063
## - x4    1     2.418 189.82 77.346
## - x2    1     6.462 193.86 77.979
## - x10   1     7.836 195.24 78.190
## - x7    1    11.065 198.47 78.683
## <none>              187.40 78.962
## - x9    1    15.758 203.16 79.384
## - x3    1    17.770 205.17 79.679
## - x1    1    18.736 206.14 79.820
## - x8    1    20.047 207.45 80.011
## - x5    1    35.864 223.26 82.215
## 
## Step:  AIC=77.04
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10
## 
##        Df Sum of Sq    RSS    AIC
## - x6    1     0.536 188.40 75.121
## - x4    1     2.363 190.23 75.411
## - x2    1     6.642 194.51 76.078
## - x10   1     7.985 195.85 76.285
## <none>              187.87 77.036
## - x7    1    14.124 201.99 77.211
## - x9    1    16.914 204.78 77.622
## - x3    1    17.815 205.68 77.754
## - x1    1    18.280 206.15 77.822
## - x8    1    20.301 208.17 78.114
## - x5    1    36.370 224.24 80.345
## 
## Step:  AIC=75.12
## y ~ x1 + x2 + x3 + x4 + x5 + x7 + x8 + x9 + x10
## 
##        Df Sum of Sq    RSS    AIC
## - x4    1     3.451 191.85 73.666
## - x2    1     6.932 195.33 74.205
## - x10   1     9.351 197.75 74.574
## <none>              188.40 75.121
## - x7    1    14.473 202.87 75.342
## - x3    1    17.802 206.20 75.830
## - x9    1    18.146 206.55 75.880
## - x1    1    18.780 207.18 75.972
## - x8    1    21.244 209.65 76.326
## - x5    1    39.332 227.73 78.809
## 
## Step:  AIC=73.67
## y ~ x1 + x2 + x3 + x5 + x7 + x8 + x9 + x10
## 
##        Df Sum of Sq    RSS    AIC
## - x2    1    10.780 202.63 73.306
## - x7    1    11.113 202.97 73.355
## <none>              191.85 73.666
## - x10   1    14.988 206.84 73.923
## - x1    1    16.602 208.46 74.156
## - x9    1    18.072 209.92 74.366
## - x3    1    21.314 213.17 74.826
## - x8    1    28.835 220.69 75.867
## - x5    1    40.323 232.18 77.389
## 
## Step:  AIC=73.31
## y ~ x1 + x3 + x5 + x7 + x8 + x9 + x10
## 
##        Df Sum of Sq    RSS    AIC
## - x7    1    10.457 213.09 72.815
## - x3    1    10.595 213.23 72.835
## - x1    1    11.998 214.63 73.032
## - x9    1    12.643 215.28 73.122
## - x10   1    13.887 216.52 73.295
## <none>              202.63 73.306
## - x8    1    27.665 230.30 75.145
## - x5    1    30.191 232.82 75.472
## 
## Step:  AIC=72.82
## y ~ x1 + x3 + x5 + x8 + x9 + x10
## 
##        Df Sum of Sq    RSS    AIC
## - x3    1    4.8720 217.96 71.494
## - x9    1    5.2049 218.29 71.539
## - x1    1    5.3212 218.41 71.555
## <none>              213.09 72.815
## - x10   1   18.3677 231.46 73.296
## - x5    1   23.3458 236.44 73.934
## - x8    1   26.0316 239.12 74.273
## 
## Step:  AIC=71.49
## y ~ x1 + x5 + x8 + x9 + x10
## 
##        Df Sum of Sq    RSS    AIC
## - x1    1     0.765 218.73 69.599
## - x9    1     5.863 223.82 70.290
## <none>              217.96 71.494
## - x10   1    20.291 238.25 72.164
## - x5    1    23.020 240.98 72.506
## - x8    1    31.634 249.59 73.559
## 
## Step:  AIC=69.6
## y ~ x5 + x8 + x9 + x10
## 
##        Df Sum of Sq    RSS    AIC
## - x9    1     5.097 223.82 68.290
## <none>              218.73 69.599
## - x5    1    40.404 259.13 72.684
## - x8    1    57.407 276.13 74.591
## - x10   1   135.105 353.83 82.029
## 
## Step:  AIC=68.29
## y ~ x5 + x8 + x10
## 
##        Df Sum of Sq    RSS    AIC
## <none>              223.82 68.290
## - x5    1    36.314 260.14 70.800
## - x8    1    52.960 276.78 72.661
## - x10   1   194.838 418.66 85.076

Para obtener un resumen del proceso se usa:

modback$anova
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11
## 
## Final Model:
## y ~ x5 + x8 + x10
## 
## 
##    Step Df   Deviance Resid. Df Resid. Dev      AIC
## 1                            18   187.4007 78.96155
## 2 - x11  1  0.4648362        19   187.8655 77.03587
## 3  - x6  1  0.5356445        20   188.4012 75.12128
## 4  - x4  1  3.4514854        21   191.8526 73.66591
## 5  - x2  1 10.7796848        22   202.6323 73.30587
## 6  - x7  1 10.4571693        23   213.0895 72.81545
## 7  - x3  1  4.8720101        24   217.9615 71.49363
## 8  - x1  1  0.7654631        25   218.7270 69.59881
## 9  - x9  1  5.0970905        26   223.8241 68.28989

Para ver la tabla de resultados del modelo modback.

summary(modback)
## 
## Call:
## lm(formula = y ~ x5 + x8 + x10, data = datis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6101 -1.9868 -0.6613  2.0369  5.8811 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.590404  11.771925   0.390   0.6998    
## x5           2.597240   1.264562   2.054   0.0502 .  
## x8           0.217814   0.087817   2.480   0.0199 *  
## x10         -0.009485   0.001994  -4.757 6.38e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.934 on 26 degrees of freedom
## Multiple R-squared:  0.8035, Adjusted R-squared:  0.7808 
## F-statistic: 35.44 on 3 and 26 DF,  p-value: 2.462e-09

Aplicacion del metodo forward

Para aplicar este metodo se debe crear un modelo vacio (empty.model) del cual iniciara el proceso. Es necesario definir un punto final de busqueda, ese punto es una formula que en este caso llamaremos horizonte. A continuacion el codigo.

empty.model <- lm(y ~ 1, data = datis)
horizonte <- formula(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11)

Se usa la funcion stepAIC y se elije trace=FALSE para que NO se muestren los detalles del proceso de seleccion.

modforw <- stepAIC(empty.model, trace=FALSE, direction="forward", scope=horizonte)
modforw$anova
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## y ~ 1
## 
## Final Model:
## y ~ x1 + x4
## 
## 
##   Step Df  Deviance Resid. Df Resid. Dev       AIC
## 1                          29  1139.1050 111.10402
## 2 + x1  1 866.49528        28   272.6097  70.20532
## 3 + x4  1  18.57161        27   254.0381  70.08861

Para ver la tabla de resultados del modelo modforw.

summary(modforw)
## 
## Call:
## lm(formula = y ~ x1 + x4, data = datis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5011 -2.1243 -0.3884  1.9964  6.9582 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.179421  18.787955   0.382    0.705    
## x1          -0.044479   0.005225  -8.513 3.98e-09 ***
## x4           3.077228   2.190294   1.405    0.171    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.067 on 27 degrees of freedom
## Multiple R-squared:  0.777,  Adjusted R-squared:  0.7605 
## F-statistic: 47.03 on 2 and 27 DF,  p-value: 1.594e-09

Como la variable \(x_4\) no es significativa entonces se puede refinar o actualizar el modelo modforw sacando \(x_4\), esto se puede realizar facilmente por medio de la funcion update asi:

modforw <- update(modforw, y ~ x1)
summary(modforw)
## 
## Call:
## lm(formula = y ~ x1, data = datis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6063 -2.0276 -0.0457  1.4531  7.0213 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 33.490010   1.535476  21.811  < 2e-16 ***
## x1          -0.047026   0.004985  -9.434 3.43e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.12 on 28 degrees of freedom
## Multiple R-squared:  0.7607, Adjusted R-squared:  0.7521 
## F-statistic:    89 on 1 and 28 DF,  p-value: 3.429e-10

Aplicacion del metodo both

Para aplicar este metodo se debe crear un modelo vacio del cual iniciara el proceso. Es necesario definir un punto final de busqueda, ese punto es una formula que en este caso llamaremos horizonte. A continuacion el codigo.

modboth <- stepAIC(empty.model, trace=FALSE, direction="both", scope=horizonte)
modboth$anova
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## y ~ 1
## 
## Final Model:
## y ~ x1 + x4
## 
## 
##   Step Df  Deviance Resid. Df Resid. Dev       AIC
## 1                          29  1139.1050 111.10402
## 2 + x1  1 866.49528        28   272.6097  70.20532
## 3 + x4  1  18.57161        27   254.0381  70.08861

El modelo modboth y modforw son el mismo.

Comparacion de los dos modelos

1. Comparando \(R^2_{Adj}\)

Para extraer el \(R^2_{Adj}\) de la tabla de resultados se usa:

summary(modback)$adj.r.squared
## [1] 0.7808368
summary(modforw)$adj.r.squared
## [1] 0.7521337

2. Comparando \(\hat{\sigma}\)

Para extraer el \(\hat{\sigma}\) de la tabla de resultados se usa:

summary(modback)$sigma
## [1] 2.934045
summary(modforw)$sigma
## [1] 3.120266

3. Comparando los residuales

par(mfrow=c(1, 2))
plot(modback, main="Backward", pch=19, cex=1, which=1)
plot(modforw, main="Forward", pch=19, cex=1, which=1)

Tarea

  1. Que patron observa en los graficos?
  2. Para cada uno de los dos modelos incluya terminos cuadraticos con el objetivo de incluir ese patron cuadratico no explicado y mostrado en los graficos de residuales anteriores.