Parcial MLG
Victoria Perdomo-Ana Palechor
12/8/2022
Modelo de predicción del salario del jugador dado su rendimiento.
La base de datos baseball contiene información del salario y medidas del rendimiento de 128 jugadores de béisbol de las Grandes Ligas de los EEUU en las temporadas de 1991 y 1992.
library(readxl)
baseball<- read_excel(file.choose(), sheet = 1)
baseball## # A tibble: 132 x 14
## salary batting OBP runs hits dooubles tripes homeruns RBI walks
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3300 0.272 0.302 69 153 21 4 31 104 22
## 2 2600 0.269 0.335 58 111 17 2 18 66 39
## 3 2500 0.249 0.337 54 115 15 1 17 73 63
## 4 2175 0.291 0.379 104 170 32 2 26 100 87
## 5 600 0.258 0.37 34 86 14 1 14 38 15
## 6 2600 0.3 0.368 69 141 22 3 19 75 53
## 7 1907 0.225 0.292 60 130 22 1 13 73 50
## 8 990 0.29 0.349 59 141 30 2 16 64 42
## 9 925 0.246 0.323 22 81 14 0 6 26 22
## 10 6100 0.302 0.391 102 174 44 6 18 100 90
## # ... with 122 more rows, and 4 more variables: strikeouts <dbl>,
## # stolenbases <dbl>, errors <dbl>, name <chr>A continuación presentamos cómo se divide la muestra para el ajuste del modelo y para la validación. Para conformar la base final eliminamos la columna “name” dado que no se tomará como covariable.
set.seed(12)
subm<-sample(132,100)
entrenamiento<-baseball[subm, ]
validacion<-baseball[-subm,]
e_final<-entrenamiento[,-c(14)]Modelo de regresión lineal múltiple
Se presenta entonces la estimación de un modelo de regresión multiple cuya variable dependiente será el salario del jugador:
##
## Call:
## lm(formula = salary ~ ., data = e_final)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1706.1 -610.9 -135.2 471.1 2833.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1151.207 839.173 1.372 0.1736
## batting 2043.649 10733.388 0.190 0.8494
## OBP -3422.738 7887.432 -0.434 0.6654
## runs 23.164 15.490 1.495 0.1384
## hits -4.292 8.898 -0.482 0.6307
## dooubles -8.242 18.678 -0.441 0.6601
## tripes -49.147 55.485 -0.886 0.3782
## homeruns -3.972 27.139 -0.146 0.8840
## RBI 25.661 11.376 2.256 0.0266 *
## walks -1.240 12.295 -0.101 0.9199
## strikeouts -3.444 4.671 -0.737 0.4630
## stolenbases 6.504 11.730 0.554 0.5807
## errors -19.125 17.791 -1.075 0.2853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 892.9 on 87 degrees of freedom
## Multiple R-squared: 0.5309, Adjusted R-squared: 0.4662
## F-statistic: 8.204 on 12 and 87 DF, p-value: 4.246e-10Validación de supuestos
Como se sabe, el modelo de regresión lineal ha de cumplir una serie de supuestos que garanticen su correcta aplicación. A continuación se validan estos supuestos para el modelo anteriormente ajustado:
Supuesto de independencia
El supuesto de la independencia de los residuos implica que los errores sean independientes entre sí y por lo tanto constituyan una variable aleatoria. Para verificarlo usaremos la prueba de Durbin-Watson. Donde el estadístico es de 1.928832 con un valor p asociado de 0.716. Por lo tanto, hay evidencia suficiente para no rechazar la hipótesis nula es decir que los residuos no están autocorrelacionados.
## Loading required package: carData## lag Autocorrelation D-W Statistic p-value
## 1 0.03535436 1.928832 0.716
## Alternative hypothesis: rho != 0Supuesto de varianza constante
El gráfico para comprobar la homocedasticidad es el ya conocido de Residuos frente a Valores ajustados. 
En el gráfico de los residuos estudentizados contra los valores ajustados, los residuos parecen no tener una tendencia, por lo que se podría considerar que los errores si tienen varianza constante. Para verificar lo anterior, usaremos la prueba de Breusch-Pagan, donde el valor p obtenido es de 0.2666, el cual es mayor al nivel de significancia = 0.05, por lo tanto hay evidencia suficiente para decir que se cumple la homocedasticidad en los errores.
## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## studentized Breusch-Pagan test
##
## data: m_bsb
## BP = 14.556, df = 12, p-value = 0.2666Supuesto de normalidad
Para hacer uso de la estimación por intervalo del modelo de regresión es necesario la normalidad en la distribución de los errores. Para comprobar esto, haremos uso de la siguiente gráfica: 
## [1] 24 92##
## One-sample Kolmogorov-Smirnov test
##
## data: m_bsb$residuals
## D = 0.072272, p-value = 0.6731
## alternative hypothesis: two-sidedAquí vemos que para los datos normales, los cuantiles muestrales y teóricos siguen aproximadamente la línea recta de referencia y la mayoría de estos puntos se encuentran dentro de la zona de confianza, por lo tanto los datos presentan un comportamiento normal. Para comprobar esto usaremos el test de kolmogorov smirnov, donde se obtuvo un valor p de 0.6731, entonces hay evidencia suficente para aceptar la hipotesis nula, por lo cual se cumple el supuesto de normalidad.
Ya comprobamos que los supuestos se cumplen para el modelo, es decir que realizando el modelo lineal podemos concluir sobre la bondad del ajuste y la independencia o dependencia de los factores sin problema.
Proceso de selección de variables comparando todos los modelos
Para obtener todos los modelos usamos la función step:
step(m_bsb)## Start: AIC=1370.96
## salary ~ batting + OBP + runs + hits + dooubles + tripes + homeruns +
## RBI + walks + strikeouts + stolenbases + errors
##
## Df Sum of Sq RSS AIC
## - walks 1 8103 69364310 1369.0
## - homeruns 1 17076 69373283 1369.0
## - batting 1 28900 69385108 1369.0
## - OBP 1 150121 69506328 1369.2
## - dooubles 1 155224 69511431 1369.2
## - hits 1 185519 69541726 1369.2
## - stolenbases 1 245089 69601296 1369.3
## - strikeouts 1 433313 69789520 1369.6
## - tripes 1 625487 69981694 1369.9
## - errors 1 921280 70277487 1370.3
## <none> 69356207 1371.0
## - runs 1 1782862 71139069 1371.5
## - RBI 1 4056594 73412801 1374.6
##
## Step: AIC=1368.97
## salary ~ batting + OBP + runs + hits + dooubles + tripes + homeruns +
## RBI + strikeouts + stolenbases + errors
##
## Df Sum of Sq RSS AIC
## - homeruns 1 11015 69375325 1367.0
## - batting 1 102857 69467167 1367.1
## - dooubles 1 150277 69514587 1367.2
## - hits 1 203813 69568123 1367.3
## - stolenbases 1 273577 69637887 1367.4
## - strikeouts 1 473292 69837602 1367.7
## - OBP 1 503549 69867859 1367.7
## - tripes 1 637398 70001707 1367.9
## - errors 1 920662 70284971 1368.3
## <none> 69364310 1369.0
## - runs 1 2184899 71549209 1370.1
## - RBI 1 4184289 73548599 1372.8
##
## Step: AIC=1366.99
## salary ~ batting + OBP + runs + hits + dooubles + tripes + RBI +
## strikeouts + stolenbases + errors
##
## Df Sum of Sq RSS AIC
## - batting 1 93785 69469110 1365.1
## - dooubles 1 146395 69521720 1365.2
## - hits 1 201664 69576989 1365.3
## - stolenbases 1 285041 69660366 1365.4
## - OBP 1 495074 69870399 1365.7
## - strikeouts 1 606234 69981558 1365.9
## - tripes 1 628239 70003563 1365.9
## - errors 1 919810 70295135 1366.3
## <none> 69375325 1367.0
## - runs 1 2374961 71750285 1368.3
## - RBI 1 8581690 77957015 1376.7
##
## Step: AIC=1365.12
## salary ~ OBP + runs + hits + dooubles + tripes + RBI + strikeouts +
## stolenbases + errors
##
## Df Sum of Sq RSS AIC
## - hits 1 107986 69577096 1363.3
## - dooubles 1 143868 69612978 1363.3
## - stolenbases 1 358404 69827515 1363.6
## - tripes 1 544284 70013394 1363.9
## - OBP 1 734934 70204044 1364.2
## - errors 1 895246 70364356 1364.4
## - strikeouts 1 941832 70410942 1364.5
## <none> 69469110 1365.1
## - runs 1 2700350 72169460 1366.9
## - RBI 1 9007985 78477095 1375.3
##
## Step: AIC=1363.28
## salary ~ OBP + runs + dooubles + tripes + RBI + strikeouts +
## stolenbases + errors
##
## Df Sum of Sq RSS AIC
## - stolenbases 1 359884 69936980 1361.8
## - dooubles 1 385037 69962134 1361.8
## - tripes 1 722960 70300056 1362.3
## - OBP 1 808196 70385292 1362.4
## - strikeouts 1 845561 70422657 1362.5
## - errors 1 1115506 70692602 1362.9
## <none> 69577096 1363.3
## - runs 1 2645761 72222857 1365.0
## - RBI 1 9277962 78855058 1373.8
##
## Step: AIC=1361.79
## salary ~ OBP + runs + dooubles + tripes + RBI + strikeouts +
## errors
##
## Df Sum of Sq RSS AIC
## - dooubles 1 592042 70529022 1360.6
## - tripes 1 743196 70680176 1360.8
## - strikeouts 1 855552 70792533 1361.0
## - OBP 1 968311 70905291 1361.2
## - errors 1 1011411 70948391 1361.2
## <none> 69936980 1361.8
## - runs 1 5950553 75887533 1368.0
## - RBI 1 9370566 79307546 1372.4
##
## Step: AIC=1360.64
## salary ~ OBP + runs + tripes + RBI + strikeouts + errors
##
## Df Sum of Sq RSS AIC
## - strikeouts 1 533880 71062902 1359.4
## - tripes 1 657452 71186474 1359.6
## - OBP 1 901702 71430724 1359.9
## - errors 1 993526 71522548 1360.0
## <none> 70529022 1360.6
## - runs 1 5434127 75963149 1366.1
## - RBI 1 8999426 79528448 1370.7
##
## Step: AIC=1359.39
## salary ~ OBP + runs + tripes + RBI + errors
##
## Df Sum of Sq RSS AIC
## - tripes 1 579462 71642364 1358.2
## - OBP 1 681754 71744655 1358.3
## - errors 1 952385 72015286 1358.7
## <none> 71062902 1359.4
## - runs 1 4980802 76043703 1364.2
## - RBI 1 8617036 79679938 1368.8
##
## Step: AIC=1358.2
## salary ~ OBP + runs + RBI + errors
##
## Df Sum of Sq RSS AIC
## - OBP 1 519206 72161570 1356.9
## - errors 1 907452 72549817 1357.5
## <none> 71642364 1358.2
## - runs 1 4651012 76293376 1362.5
## - RBI 1 11467611 83109975 1371.0
##
## Step: AIC=1356.92
## salary ~ runs + RBI + errors
##
## Df Sum of Sq RSS AIC
## - errors 1 918846 73080416 1356.2
## <none> 72161570 1356.9
## - runs 1 4307976 76469547 1360.7
## - RBI 1 12341568 84503138 1370.7
##
## Step: AIC=1356.19
## salary ~ runs + RBI
##
## Df Sum of Sq RSS AIC
## <none> 73080416 1356.2
## - runs 1 3951944 77032360 1359.5
## - RBI 1 12178286 85258702 1369.6##
## Call:
## lm(formula = salary ~ runs + RBI, data = e_final)
##
## Coefficients:
## (Intercept) runs RBI
## 346.12 11.69 20.13Por lo anterior, los 3 modelos que planteamos teniendo en cuenta los valores AIC son:
mod1<-lm(salary~OBP + runs + RBI + errors, data=baseball)
summary(mod1)##
## Call:
## lm(formula = salary ~ OBP + runs + RBI + errors, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1716.6 -586.4 -74.7 467.8 2577.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 462.515 619.660 0.746 0.45680
## OBP -469.588 2116.066 -0.222 0.82474
## runs 15.185 4.891 3.105 0.00235 **
## RBI 19.272 4.033 4.778 4.81e-06 ***
## errors -15.504 13.871 -1.118 0.26579
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 850.1 on 127 degrees of freedom
## Multiple R-squared: 0.5336, Adjusted R-squared: 0.5189
## F-statistic: 36.33 on 4 and 127 DF, p-value: < 2.2e-16mod2<-lm(salary ~ runs + RBI + errors, data=baseball)
summary(mod2)##
## Call:
## lm(formula = salary ~ runs + RBI + errors, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1719.02 -585.19 -85.25 482.62 2571.01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 331.311 184.838 1.792 0.075425 .
## runs 14.619 4.157 3.517 0.000605 ***
## RBI 19.401 3.976 4.879 3.1e-06 ***
## errors -15.568 13.816 -1.127 0.261932
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 846.9 on 128 degrees of freedom
## Multiple R-squared: 0.5334, Adjusted R-squared: 0.5225
## F-statistic: 48.78 on 3 and 128 DF, p-value: < 2.2e-16mod3<-lm(salary~runs+RBI, data= baseball)
summary(mod3)##
## Call:
## lm(formula = salary ~ runs + RBI, data = baseball)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1755.52 -586.09 -70.49 499.07 2467.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 262.807 174.736 1.504 0.135020
## runs 14.369 4.155 3.458 0.000738 ***
## RBI 19.037 3.967 4.799 4.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 847.8 on 129 degrees of freedom
## Multiple R-squared: 0.5288, Adjusted R-squared: 0.5215
## F-statistic: 72.39 on 2 and 129 DF, p-value: < 2.2e-16Notamos que la diferencia entre el coeficiente de determinación \(R^2\) entre los cuatro modelos es pequeña. El modelo cuyo \(R^2\) es más grande es el modelo 1 (mod1): salary ~ OBP + runs + RBI + errors, es decir que las covariables explican el 53.36% de la variabilidad de los datos del salario de los jugadores. Sin embargo, se aclara que en este modelo no todas las variables son significativas.
Sólo en el modelo 3 (mod 3) todas las covariables son significativas y estas son runs y RBI. Es evidente la mejoría del modelo ya que en el modelo completo sólo la variable RBI era significativa según el valor p asociado.
Predicciones
Usando la muestra definida anteriormente para la validación se predicen valores del salario de los jugadores usando el modelo completo y los 3 modelos propuestos. Adicionalmente, se calcula el error cuadrado medio de predicción que usaremos para comparar los modelos.
predict(m_bsb, validacion)## 1 2 3 4 5 6 7 8
## 1038.8842 1190.8400 1126.4241 918.7851 3120.6843 2175.2735 861.1170 1678.5151
## 9 10 11 12 13 14 15 16
## 2693.1903 2604.3916 1538.2151 2652.3358 1703.6038 848.0763 3233.1653 2203.1436
## 17 18 19 20 21 22 23 24
## 2149.4648 2011.3190 898.4217 3003.2141 2746.3406 1877.1911 1037.9492 2018.0379
## 25 26 27 28 29 30 31 32
## 2525.2252 1891.1971 1029.9674 3175.5191 3213.7015 3221.3159 2822.0596 880.0949predict(mod1,validacion)## 1 2 3 4 5 6 7 8
## 1068.4524 1263.2828 1177.0305 821.0273 3200.9489 2032.9705 1218.2463 1749.1098
## 9 10 11 12 13 14 15 16
## 3004.7167 1972.8337 1589.1759 2707.2516 1616.7917 772.6433 3377.8408 2491.9070
## 17 18 19 20 21 22 23 24
## 2651.2195 2213.6290 1137.3692 3279.1573 2498.3966 1802.8877 1164.0994 1926.2514
## 25 26 27 28 29 30 31 32
## 2515.4066 1848.0859 926.2793 2971.8255 3228.5336 3155.6673 2820.7552 837.1018predict(mod2,validacion)## 1 2 3 4 5 6 7 8
## 1079.5078 1276.7716 1203.5420 805.0085 3197.8469 2039.6616 1236.5141 1763.8319
## 9 10 11 12 13 14 15 16
## 3000.3182 1973.1369 1577.3534 2705.8283 1619.6881 763.1231 3397.8341 2496.1544
## 17 18 19 20 21 22 23 24
## 2627.9919 2205.2747 1147.4101 3302.6916 2491.9522 1803.4496 1191.6378 1925.2161
## 25 26 27 28 29 30 31 32
## 2504.6689 1852.3943 922.7723 2986.0003 3222.0306 3172.8509 2808.6435 855.0752predict(mod3,validacion)## 1 2 3 4 5 6 7 8
## 1073.8785 1237.3343 1180.9544 773.9581 3092.9511 2184.5050 1289.4116 1700.3454
## 9 10 11 12 13 14 15 16
## 3250.6435 1921.6419 1532.2217 2595.5182 1833.9680 855.1382 3487.6902 2434.5392
## 17 18 19 20 21 22 23 24
## 2564.2242 2133.5233 1156.1542 3287.2564 2475.8200 1708.2205 1107.6493 1904.3517
## 25 26 27 28 29 30 31 32
## 2550.2205 1832.8725 858.7106 2991.6387 3116.6558 3068.1509 2771.4377 945.2890predicciones<-data.frame(real=validacion$salary,pm_bsb=predict(m_bsb, validacion),pmod1=predict(mod1, validacion),pmod2=predict(mod2, validacion),pmod3=predict(mod3, validacion))
mean ((predicciones$real - predicciones$pm_bsb)^2)## [1] 625218.7mean ((predicciones$real - predicciones$pmod1)^2)## [1] 610640.7mean ((predicciones$real - predicciones$pmod2)^2)## [1] 604816.6mean ((predicciones$real - predicciones$pmod3)^2)## [1] 604315.6El MSE mide la precisión de predicción del modelo. Por su fórmula, cuánto más pequeño sea este valor, habrá más precisión en la predicción de valores, en este caso del salario de los jugadores.
Los MSE menores corresponden al modelo 2 (mod2) y modelo 3 (mod3), aunque la diferencia entre ellos es muy pequeña, el modelo salary~runs+RBI es el que se considera tiene mayor precisión en la predicción del salario de los jugadores.
Sobre base de datos fertility
library(Countr)## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod carfertilidad<-fertility
fert_final<-fertilidad[,c(1,3,5,6,7,8,9)]
fert_final$university<-as.factor(fert_final$university)
fert_final$religion<-as.factor(fert_final$religion)
fert_final$rural<-as.factor(fert_final$rural)Ajuste modelo Poisson
Se ajusta un modelo Poisson con variable dependiente “children” y las variables years school, university, religion, year birthday, rural, age married como covariables.
md_fert<-glm(children~.,family="poisson",data=fert_final)
summary(md_fert)##
## Call:
## glm(formula = children ~ ., family = "poisson", data = fert_final)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1942 -0.7075 -0.1098 0.4416 4.2936
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.670323 0.280905 5.946 2.74e-09 ***
## years_school -0.044555 0.028134 -1.584 0.1133
## universityyes 0.135115 0.146449 0.923 0.3562
## religionMuslim 0.150217 0.067106 2.239 0.0252 *
## religionOther 0.621705 0.083801 7.419 1.18e-13 ***
## religionProtestant 0.011385 0.069016 0.165 0.8690
## year_birth 0.002884 0.002292 1.258 0.2083
## ruralyes 0.078748 0.037221 2.116 0.0344 *
## age_marriage -0.030950 0.006505 -4.758 1.96e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1204.3 on 1242 degrees of freedom
## Residual deviance: 1057.4 on 1234 degrees of freedom
## AIC: 4244.7
##
## Number of Fisher Scoring iterations: 5Primero, aclaramos que la variable religión tiene 4 categorias que son: catholic, muslim, other y protestant. El modelo toma como base católica, esto es, cuando la mujer profese esta religión las covariables religionMuslim,religionOther y religionProtestant toman el valor de 0. Cuando no es así, por ejemplo, religionProtestant toma el valor de 1 cuando la mujer pertenece a dicha religión y 0 cuando no. Análogamente, se definen las covariables religionMuslim y religionOther.
La variable university toma como base la opción “yes”, por lo que, universityyes toma el valor de 1 cuando la mujer estudió en la universidad y 0 cuando no. Y la variable rural tiene como base la opción “yes” por tal razón ruralyes toma el valor 1 cuando la mujer reside en zona rural y 0 cuando no es el caso.
Del resumen del modelo, podemos afirmar que la covariable religionOther y age_marriage son significativas. Mientras que, years_school, universityyes y year_birth se concluye que no son significativas observando los valores-p asociados.
Sobre los coeficientes podemos deducir que si la mujer pertenece a otra religión y tiene educación universitaria se estima que tenga un mayor número de hijos que una mujer que pertenezca a la religión católica, musulman o protestante y que no tenga educación superior.
Sobredispersión
Es necesario descartar que el modelo presenta sobredispersión ya que así sabremos si el valor de la varianza y esperanza no difieren. Cuando estos valores no coinciden se dice hay sobredispersión.
Para probar si hay sobredispersión en los datos, usaremos la función dispersiontest que prueba la hipótesis nula de equidispersión en un modelo Poisson frente a la alternativa de sobredispersión y/o subdispersión.
library(AER)## Loading required package: sandwich## Loading required package: survivaldispersiontest(md_fert)##
## Overdispersion test
##
## data: md_fert
## z = -3.3515, p-value = 0.9996
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 0.8249273Ya que el valor p= 0.9996 concluimos que hay evidencia suficiente para pensar que no hay sobredispersión. No es necesario ajustar un modelo binomial negativo.