Ejemplo 2


Se tiene información sobre el desempeño de los equipos en la liga de Futbol Americano (1976).

library("xlsx")
datos=read.xlsx("futlbol.xlsx",1)
head(datos)
##    y   x1   x2   x3   x4 x5  x6   x7   x8   x9
## 1 10 2113 1985 38.9 64.7  4 868 59.7 2205 1917
## 2 11 2003 2855 38.8 61.3  3 615 55.0 2096 1575
## 3 11 2957 1737 40.1 60.0 14 914 65.6 1847 2175
## 4 13 2285 2905 41.6 45.3 -4 957 61.4 1903 2476
## 5 10 2971 1666 39.2 53.8 15 836 66.1 1457 1866
## 6 11 2309 2927 39.7 74.1  8 786 61.0 1848 2339

\(y:\) Juegos ganados (por temporada de 14juegos)

\(x_1:\) Yardas por tierra (temporada)

\(x_2:\) Yardas por aire (temporada)

\(x_3:\) Promedio de pateo (yardas/patada)

\(x_4:\) Porcentaje de goles de campo (GC hechos/GC intentados, temporada)

\(x_5:\) Diferencia de pérdidas de balón (pérdidas ganadas/pérdidas perdidas)

\(x_6:\) Yardas de castigo (temporada)

\(x_7:\) Porcentaje de carreras (jugadas por tierra/jugadas totales)

\(x_8:\) Yardas por tierra del contrario (temporada)

\(x_9:\) Yardas por aire del contrario (temporada)

Ajuste un modelo de RLM


Solución


Ajustaremos un modelo de RLM considerando lo siguiente:

El método backward (eliminación hacia atrás) consiste en comenzar con un modelo que incluya las \(k\) variables regresoras y luego ir eliminando variable por variable de acuerdo a cierto criterio.

El AIC (Criterio de Información de Akaike) es una medida de la calidad de un modelo estadístico, considera la bondad de ajuste del modelo y su complejidad.

\[ AIC=-2log(L)+2*n_{parametros} \]

El mejor modelo es aquel que tenga un \(AIC\) pequeño.


Primer Paso: Relación entre Variables (Indicios)

Podemos hacer gráficas de dispersión para ver el comportamiento de los datos.

library("tidyverse")
library("GGally")

attach(datos)

ggpairs(datos, upper = list(continuous = "smooth"),
        diag = list(continuous = "bar"),lower=list(continuous="cor"))


Segundo Paso: Generación de Modelos y Elección del Mejor Modelo

Emplearemos comandos de R:

library("MASS")

#generamos el modelo con todas las variables regresoras
modelo_completo=lm(y~.,data=datos)

#hagamos un análisis del modelo completo (nadamás para ver que está pasando)
summary(modelo_completo)
## 
## Call:
## lm(formula = y ~ ., data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0408 -0.6802 -0.1131  0.9835  2.9785 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.292e+00  1.281e+01  -0.569 0.576312    
## x1           8.124e-04  2.006e-03   0.405 0.690329    
## x2           3.631e-03  8.410e-04   4.318 0.000414 ***
## x3           1.222e-01  2.590e-01   0.472 0.642750    
## x4           3.189e-02  4.160e-02   0.767 0.453289    
## x5           1.511e-05  4.684e-02   0.000 0.999746    
## x6           1.590e-03  3.248e-03   0.490 0.630338    
## x7           1.544e-01  1.521e-01   1.015 0.323547    
## x8          -3.895e-03  2.052e-03  -1.898 0.073793 .  
## x9          -1.791e-03  1.417e-03  -1.264 0.222490    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.83 on 18 degrees of freedom
## Multiple R-squared:  0.8156, Adjusted R-squared:  0.7234 
## F-statistic: 8.846 on 9 and 18 DF,  p-value: 5.303e-05


Ahora si, emplearemos el método backward con el criterio AIC.

library("MASS")

#la función stepAIC nos permite generar modelos paso a paso con eliminación hacia atrás y toma el criterio de AIC.
modelo_backward_AIC= stepAIC(modelo_completo,direction = "backward")
## Start:  AIC=41.48
## y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
## 
##        Df Sum of Sq     RSS    AIC
## - x5    1     0.000  60.293 39.476
## - x1    1     0.549  60.842 39.730
## - x3    1     0.746  61.039 39.821
## - x6    1     0.803  61.096 39.847
## - x4    1     1.968  62.261 40.376
## - x7    1     3.451  63.744 41.035
## <none>               60.293 41.476
## - x9    1     5.348  65.642 41.856
## - x8    1    12.072  72.365 44.587
## - x2    1    62.448 122.741 59.380
## 
## Step:  AIC=39.48
## y ~ x1 + x2 + x3 + x4 + x6 + x7 + x8 + x9
## 
##        Df Sum of Sq     RSS    AIC
## - x1    1     0.553  60.846 37.732
## - x3    1     0.750  61.043 37.822
## - x6    1     0.818  61.111 37.854
## - x4    1     2.053  62.346 38.414
## - x7    1     3.859  64.152 39.213
## <none>               60.293 39.476
## - x9    1     5.351  65.644 39.857
## - x8    1    12.086  72.379 42.592
## - x2    1    66.979 127.272 58.395
## 
## Step:  AIC=37.73
## y ~ x2 + x3 + x4 + x6 + x7 + x8 + x9
## 
##        Df Sum of Sq     RSS    AIC
## - x6    1     0.690  61.536 36.048
## - x3    1     1.715  62.561 36.510
## - x4    1     3.051  63.897 37.102
## <none>               60.846 37.732
## - x9    1     4.852  65.698 37.880
## - x7    1     8.961  69.807 39.579
## - x8    1    16.599  77.445 42.486
## - x2    1    67.010 127.856 56.524
## 
## Step:  AIC=36.05
## y ~ x2 + x3 + x4 + x7 + x8 + x9
## 
##        Df Sum of Sq     RSS    AIC
## - x3    1     1.726  63.262 34.822
## - x4    1     2.767  64.303 35.279
## <none>               61.536 36.048
## - x9    1     4.831  66.367 36.164
## - x7    1     9.390  70.926 38.024
## - x8    1    18.314  79.851 41.343
## - x2    1    66.447 127.984 54.552
## 
## Step:  AIC=34.82
## y ~ x2 + x4 + x7 + x8 + x9
## 
##        Df Sum of Sq     RSS    AIC
## - x4    1     1.743  65.004 33.583
## <none>               63.262 34.822
## - x9    1     5.629  68.891 35.209
## - x8    1    17.701  80.962 39.730
## - x7    1    18.583  81.845 40.033
## - x2    1    75.598 138.860 54.835
## 
## Step:  AIC=33.58
## y ~ x2 + x7 + x8 + x9
## 
##        Df Sum of Sq     RSS    AIC
## <none>               65.004 33.583
## - x9    1     4.866  69.870 33.604
## - x7    1    16.908  81.913 38.057
## - x8    1    23.299  88.303 40.160
## - x2    1    82.892 147.897 54.601

Observemos el mejor modelo, de acuerdo a lo arrojado en stepAIC():

summary(modelo_backward_AIC)
## 
## Call:
## lm(formula = y ~ x2 + x7 + x8 + x9, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3519 -0.5612 -0.0856  0.6972  3.2802 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.8217034  7.7847061  -0.234  0.81705    
## x2           0.0038186  0.0007051   5.416 1.67e-05 ***
## x7           0.2168941  0.0886759   2.446  0.02252 *  
## x8          -0.0040149  0.0013983  -2.871  0.00863 ** 
## x9          -0.0016349  0.0012460  -1.312  0.20244    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.681 on 23 degrees of freedom
## Multiple R-squared:  0.8012, Adjusted R-squared:  0.7666 
## F-statistic: 23.17 on 4 and 23 DF,  p-value: 8.735e-08

A través de la tabla ANOVA se tiene:

anova(modelo_backward_AIC)
## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## x2         1  76.193  76.193 26.9589 2.900e-05 ***
## x7         1 139.501 139.501 49.3585 3.693e-07 ***
## x8         1  41.400  41.400 14.6483  0.000863 ***
## x9         1   4.866   4.866  1.7216  0.202435    
## Residuals 23  65.004   2.826                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


De los resultados anteriores podemos refinar el mejor modelo obtenido:

modelo_refinado=lm(y~0+x2+x7+x8)
summary(modelo_refinado)
## 
## Call:
## lm(formula = y ~ 0 + x2 + x7 + x8)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9707 -0.6953 -0.2601  1.0570  3.7442 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## x2  0.0035213  0.0005969   5.899 3.73e-06 ***
## x7  0.1747501  0.0266986   6.545 7.41e-07 ***
## x8 -0.0050642  0.0006581  -7.695 4.74e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.674 on 25 degrees of freedom
## Multiple R-squared:  0.9584, Adjusted R-squared:  0.9535 
## F-statistic: 192.2 on 3 and 25 DF,  p-value: < 2.2e-16


Tercer Paso: Validación de Supuestos

Ya que hemos obtenido el mejor modelo debemos hacer validación de supuestos.


Relación lineal entre regresores y la variable respuesta

library(ggplot2)
library(gridExtra)

plot1 <- ggplot(data = datos, aes(x2, modelo_refinado$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()

plot2 <- ggplot(data = datos, aes(x7, modelo_refinado$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()

plot3 <- ggplot(data = datos, aes(x8, modelo_refinado$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()

grid.arrange(plot1, plot2,plot3)


Todos los supuestos sobre los Residuales

library("tidyverse")

residuales=data.frame(ajustados=modelo_refinado$fitted.values,residuales_estand=rstandard(modelo_refinado))

ggplot(data=residuales,aes(x=ajustados,y=residuales_estand))+
  geom_point(color="darkblue",size=8)+geom_hline(yintercept = 0)+
  geom_hline(yintercept=2,color="red")+
  geom_hline(yintercept = -2,color="red")

También podemos hacer gráficas como QQPLOT e histogramas para ver gráficamente la normalidad de los residuales.


Pruebas Formales
library("nortest")
ad.test(rstandard(modelo_refinado))
## 
##  Anderson-Darling normality test
## 
## data:  rstandard(modelo_refinado)
## A = 0.519, p-value = 0.1715
library("lmtest")
bptest(modelo_refinado)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_refinado
## BP = 1.5914, df = 2, p-value = 0.4513
library("lmtest")
dwtest(modelo_refinado)
## 
##  Durbin-Watson test
## 
## data:  modelo_refinado
## DW = 1.4702, p-value = 0.0683
## alternative hypothesis: true autocorrelation is greater than 0


Multicolinealidad

Podemos calcular el número de condición (coeficiente \(\kappa\)):

X1=scale(datos[,c(-1)])
A=t(X1)%*%X1
kappa=max(eigen(A)$values)/min(eigen(A)$values)
kappa
## [1] 26.90626

La literatura indica que, en general, si el número de condición es menor que 100 entonces no hay problemas graves de multicolinealidad; si el número de condición oscila entre 100 a 1000 implican multicolinealidad de moderada a fuerte; si \(\kappa>1000\) entonces hay indicio de una fuerte multicolinealidad.


También se sugiere utilizar la matriz de correlaciones entre las variables regresoras para ver si tenemos problemas de multicolinealidad.

library("corrplot")

corrplot(cor(dplyr::select(datos,x2,x7,x8)),method="number",tl.col = "black")

Falta ver si hay datos influyentes y atípicos; eso ya les toca a ustedes


Conclusiones Respecto al Problema

¿Cuáles son las conclusiones respecto al problema?


Conclusiones Respecto a la Metodología Empleada

  • El método backward es una buena opción para no descartar nada desde el principio.

  • El \(AIC\) es una buena medida para comparar modelos en términos de ajuste y complejidad.

  • Existen otros criterios de comparación y métodos de generación de modelos.