Modelo lineal general-Parcial 2
Karen Arteaga, Franklin Yandy
2022-08-15
Análisis exploratorio
En esta ocasión disponemos la base de datos baseball que contiene información del salario y medidas del rendimiento de 132 jugadores de béisbol de las Grandes Ligas de los EEUU en las temporadas de 1991 y 1992. En primer lugar haremos una análisis descriptivo de algunas variables.
El salario de los jugadores está entre 109 mil y 6100 mil dólares,además el 50% de los jugadores tiene un salario menor a 2179.
Ahora realizamos un grafico de correlación.
En el gráfico se puede observar que el salario está fuertemente relacionado con las variables runs, homeruns y RBI, además, podemos notar que las variables homeruns y RBI presentan el coeficiente de correlación más alto 0.89.
Selección de datos de entrenamiento y datos de prueba
Se selecciona una submuestra de 100 datos con los cuales se realizará el modelo. Los 32 datos restantes serán utilizados para evaluar los resultados del modelo.
set.seed(2831)
seleccion<-(sample(1:132,100,replace = F))
Datos_entrenamiento<-Datos[c(seleccion),-14]
Datos_prueba<-Datos[-seleccion,-14]A continuación se procede a realizar un modelo de regresión lineal tomando como variable respuesta el salario y como predictoras el resto de variables.
Modelo de regresión lineal
mod1<-lm(salary~.,data=Datos_entrenamiento)
summary(mod1)##
## Call:
## lm(formula = salary ~ ., data = Datos_entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1965.16 -453.49 6.74 478.82 2371.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 921.279 778.273 1.184 0.23974
## batting 7460.855 9392.105 0.794 0.42914
## OBP -8034.232 7168.674 -1.121 0.26548
## runs 8.685 14.211 0.611 0.54270
## hits -4.704 7.561 -0.622 0.53547
## dooubles -1.099 16.694 -0.066 0.94765
## tripes -4.517 45.640 -0.099 0.92140
## homeruns -9.473 25.830 -0.367 0.71470
## RBI 32.835 10.361 3.169 0.00211 **
## walks 14.040 11.418 1.230 0.22216
## strikeouts -4.967 4.567 -1.087 0.27983
## stolenbases 14.322 11.484 1.247 0.21572
## errors -24.572 15.613 -1.574 0.11916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 826.1 on 87 degrees of freedom
## Multiple R-squared: 0.6149, Adjusted R-squared: 0.5618
## F-statistic: 11.58 on 12 and 87 DF, p-value: 1.595e-13El 61% de la variabilidad del salario está siendo explicada por el modelo completo, este modelo tiene un error estimado bastante alto de 826.1, por lo que es necesario hacerle una transformación.
Tras la realización del modelo se obtiene que RBI (número de carreras impulsadas) es la única variable que tiene un aporte significativo en la explicación del salario.
7 Evaluación de supuestos
Homocedasticidad
res.stud<-studres(mod1)
mod.fit<-mod1$fitted.values
par(mfrow=c(1,2))
plot(mod.fit,res.stud,ylab="Residuos estudentizados",xlab="Valores ajustados")
abline(h=0,lty=2)
lines(lowess(res.stud~mod.fit), col = 2)
plot(mod.fit,abs(res.stud),ylab="|Residuos estudentizados|",xlab="Valores ajustados")
lines(lowess(abs(res.stud)~mod.fit), col = 2)
En la figura anterior se observa como la varianza aumenta a medida que lo hacen los valores ajustados, decimos así que el modelo no cumple el supuesto de homocedasticidad.
bptest(mod1)##
## studentized Breusch-Pagan test
##
## data: mod1
## BP = 29.569, df = 12, p-value = 0.003241De la salida anterior se observa que el valor-P es menor que el nivel de significancia usual de 5%, por lo tanto, hay evidencias para decir que no se cumple la homocedasticidad de los errores.
Normalidad
car::qqPlot(mod1)
## 10 86
## 7 78La anterior gráfica muestra los cuantiles muestrales contra los cuantiles teóricos que se esperarían bajo normalidad. Se observa que son algunos datos los que no afecta el cumplimiento de los supuestos.
Transformaciones
Gráficamente el incumplimeinto de los supuestos no es muy notorio pero las prueba nos sugiere que hay heterocedasticidad. Realizaremos una transformación por BoxCox para corregir los supuestos.
boxcox = MASS::boxcox(mod1,lambda=seq(-3,3,length.out = 1000),
ylab='log-verosimilitud')
boxcox$x[boxcox$y ==max(boxcox$y)]## [1] 0.5015015El valor al cual se elevará la variable respuesta es aquel que maximiza la logverosimilitud. De modo que realizamos la transformación Y0.5 para un nuevo modelo.
modBox<-lm(salary^0.5~.,data=Datos_entrenamiento)
summary(modBox)##
## Call:
## lm(formula = salary^0.5 ~ ., data = Datos_entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.4795 -6.5264 0.1644 5.2545 21.1349
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.48139 8.58891 3.782 0.000285 ***
## batting 92.15411 103.64997 0.889 0.376407
## OBP -106.53404 79.11249 -1.347 0.181602
## runs 0.11742 0.15683 0.749 0.456059
## hits -0.02447 0.08344 -0.293 0.769999
## dooubles -0.06797 0.18423 -0.369 0.713072
## tripes -0.19195 0.50368 -0.381 0.704061
## homeruns -0.06175 0.28505 -0.217 0.828995
## RBI 0.31987 0.11434 2.798 0.006338 **
## walks 0.16883 0.12601 1.340 0.183791
## strikeouts -0.05114 0.05040 -1.015 0.313117
## stolenbases 0.15754 0.12674 1.243 0.217195
## errors -0.32057 0.17230 -1.861 0.066189 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.116 on 87 degrees of freedom
## Multiple R-squared: 0.6304, Adjusted R-squared: 0.5794
## F-statistic: 12.36 on 12 and 87 DF, p-value: 3.03e-14El modelo transformado presenta una mejora significativa ya que la estimación del error disminuyó , además se aumentó en 2% el R^{2} respecto al modelo completo sin transformación. Para este modelo también realizamos la prueba de breush pagan, comprobamos que sí se cumple el supuesto de homocedasticidad y por otro lado el de normalidad con la prueba Shapiro Wilks.
##
## studentized Breusch-Pagan test
##
## data: modBox
## BP = 20.489, df = 12, p-value = 0.05838##
## Shapiro-Wilk normality test
##
## data: res.stud2
## W = 0.9925, p-value = 0.856Segunda parte
Para realizar la selección de variables en primer lugar haremos uso de la tecnica random forest para medir la importancia de las variables en funcion del salario.
modeloRF <- randomForest(salary~.,data=Datos_entrenamiento,
ntree=500,importance=TRUE,maxnodes=10)
varImpPlot(modeloRF)
Este método sugiere que las variables más importantes para el modelo son RBi, runs y hits, por el contrario las menos importantes son batting errors y tripes. Teniendo en cuenta lo anterior ajustamos un modelo en funcion de estas variables. ### Modelo Random Forest
modRF1<-lm(salary~RBI+runs+hits+homeruns,data=Datos_entrenamiento)
summary(modRF1)##
## Call:
## lm(formula = salary ~ RBI + runs + hits + homeruns, data = Datos_entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1901.30 -478.96 -47.84 504.98 2145.67
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.712 224.891 0.408 0.684335
## RBI 33.265 9.424 3.530 0.000643 ***
## runs 19.935 7.596 2.624 0.010115 *
## hits -5.623 4.901 -1.147 0.254115
## homeruns -28.829 20.880 -1.381 0.170596
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 818.9 on 95 degrees of freedom
## Multiple R-squared: 0.5868, Adjusted R-squared: 0.5694
## F-statistic: 33.73 on 4 and 95 DF, p-value: < 2.2e-16Este modelo explica el 58% dela variabilidad del salario, y disminuyo el error estimado del modelo completo,es ahora de 819.9.
Modelo StepAIC-Backward
modback$anova## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## salary ~ batting + OBP + runs + hits + dooubles + tripes + homeruns +
## RBI + walks + strikeouts + stolenbases + errors
##
## Final Model:
## salary ~ RBI + walks + strikeouts + stolenbases + errors
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 87 59366771 1355.408
## 2 - dooubles 1 2959.176 88 59369730 1353.412
## 3 - tripes 1 5445.909 89 59375176 1351.422
## 4 - homeruns 1 87298.712 90 59462475 1349.569
## 5 - runs 1 201012.804 91 59663488 1347.906
## 6 - hits 1 88532.946 92 59752021 1346.054
## 7 - batting 1 609521.350 93 60361542 1345.069
## 8 - OBP 1 685653.045 94 61047195 1344.199modSAIC<-lm(salary ~ RBI + walks + strikeouts + stolenbases + errors,data = Datos_entrenamiento)
summary(modSAIC)##
## Call:
## lm(formula = salary ~ RBI + walks + strikeouts + stolenbases +
## errors, data = Datos_entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1925.48 -470.45 -22.38 500.53 2534.63
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 295.616 211.867 1.395 0.1662
## RBI 32.693 4.088 7.998 3.22e-12 ***
## walks 7.870 4.387 1.794 0.0760 .
## strikeouts -5.464 3.382 -1.616 0.1095
## stolenbases 19.353 8.549 2.264 0.0259 *
## errors -25.360 14.818 -1.711 0.0903 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 805.9 on 94 degrees of freedom
## Multiple R-squared: 0.604, Adjusted R-squared: 0.583
## F-statistic: 28.68 on 5 and 94 DF, p-value: < 2.2e-16Este modelo fue hallado mediante la estrategia backward en la cual se empieza con un modelo saturado y se van descartando variables. El modelo explica el 60% de la variabilidad del salario y tiene un error de 805.9
Modelo stepwise-forward
Este modelo fue hallado mediante la estrategia forward donde partido de un modelo “vacio” y se va agregando varibles. Este modelo explica el 58% de la variabilidad del salario y tiene un error de 811.8.
modStep$anova## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 NA NA 99 154177219 1426.844
## 2 + RBI -1 83658394 98 70518825 1350.622
## 3 + runs -1 5403154 97 65115671 1344.651
## 4 + errors -1 1842461 96 63273210 1343.780modSt<-lm(salary ~ RBI + runs + errors,data=Datos_entrenamiento)
summary(modSt)##
## Call:
## lm(formula = salary ~ RBI + runs + errors, data = Datos_entrenamiento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1793.17 -537.74 -44.82 459.60 2574.31
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 196.778 204.905 0.960 0.33930
## RBI 22.246 4.593 4.843 4.89e-06 ***
## runs 14.484 4.849 2.987 0.00358 **
## errors -24.872 14.876 -1.672 0.09779 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 811.8 on 96 degrees of freedom
## Multiple R-squared: 0.5896, Adjusted R-squared: 0.5768
## F-statistic: 45.97 on 3 and 96 DF, p-value: < 2.2e-16###Conclusión A partir de los resultados de todos los modelos podriamos consideras como el mejor modelo el propuesto por StepAIC ya que en comparacion con los demás explica el 60% de la variabilidad del salario y es el que en comparacion con los demás tiene un error estimado más bajo.
Predicciones usando datos de prueba
A continuación comprobaremos qué tan efectivos son los modelos generados utilizando los datos de prueba.
Para cada uno de ellos se realizará la predicción de los valores y se calculará el valor de cuadrado medio de las predicciones. Finalmente, se compararan los valores obtenidos.
Modelo completo
dataPredict<-Datos_prueba[,-1]
predictCom<-predict(modBox,newdata = dataPredict)
MSE_Com<-sum((predictCom-Datos_prueba$salary)^2)/32
MSE_Com## [1] 5502515Modelo StepWise
predictST<-predict(modSt,newdata = dataPredict)
MSE_ST<-sum((predictST-Datos_prueba$salary)^2)/32
MSE_ST## [1] 931049.7Modelo Random Forest
predictRF<-predict(modRF1,newdata = dataPredict)
MSE_RF<-sum((predictRF-Datos_prueba$salary)^2)/32
MSE_RF## [1] 957161.6Modelo Step AIC
predictAIC<-predict(modSAIC,newdata = dataPredict)
MSE_SAIC<-sum((predictAIC-Datos_prueba$salary)^2)/32
MSE_SAIC## [1] 949040.5Tras realizar las predicciones se obtiene que de los modelos realizados mediante selección de variables, ha sido más preciso el modelo obtenido mediante la función StepWise, el siguiente modelo ha sido el obtenido mediante AIC, seguido del de Random Forest.
Sin embargo, el modelo completo al cual se le ha realizado transformación por boxcox ha tenido un valor mucho más bajo que el de los demás.
Concluimos que la transformación ha tenido un gran efecto positivo en la mejora del modelo planteado inicialmente.
Base de datos de fertilidad
data("Fertility")
Fertilidad<-fertility[,c(1,3,5,6,7,8,9)]Esta base de datos cuenta con un total de 1243 datos y 7 variables.
Para esta base de datos realizaremos un modelo poisson para predecir el número de hijos en función del resto de covariables.
mean(Fertilidad$children);var(Fertilidad$children)## [1] 2.383749## [1] 2.330074Esta variable es apropiada para este tipo de modelos debido a que representa un conteo.
Además esta variable tiene un comportamiento muy similar a una variable Poisson, esto debido a que su media y varianza toman valores muy similares.
De modo que procedemos a realizar el modelo Poisson.
modPoisson<-glm(Fertilidad$children~.,family = poisson,data=Fertilidad)
summary(modPoisson)##
## Call:
## glm(formula = Fertilidad$children ~ ., family = poisson, data = Fertilidad)
##
## 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: 5Tras realizar el ajuste obtenemos que las covariables religión Muslim, otra religión, la zona rural y la edad en la que contrajo matrimonio son significantes para explicar el número de hijos.
Ahora, dado que el modelo es de la forma \(Log(y)=\alpha+\beta(x)\)
para realizar la interpretación de los coeficientes se le aplica \(exp\) a los coeficientes de la siguiente manera, con ello se obtiene \(y=exp(\alpha)*exp(\beta(x))\)
exp(modPoisson$coefficients)## (Intercept) years_school universityyes religionMuslim
## 5.3138849 0.9564231 1.1446679 1.1620863
## religionOther religionProtestant year_birth ruralyes
## 1.8621009 1.0114505 1.0028879 1.0819322
## age_marriage
## 0.9695245Por cada unidad que incrementa X, se obtiene un efecto de multiplicar por \(exp(\beta)\) a la media de Y.
Entre más cercanos son estos valores a 1, indican que no hay relación con el número de hijos.
De modo que cuando \(exp(\beta)>1\) se tiene un efecto de incremento en la cantidad de hijos.
cuando \(exp(\beta)<1\) se tiene un efecto de decrecimiento en la cantidad de hijos.
Para este caso se tiene que cuando la religión es “otros”, la cantidad de hijos se multiplica por 1.8.
Por otro lado, por cada año en la edad en la que contrajo matrimonio la cantidad de hijos se multiplica por 0.96, es decir, la relación es negativa, a menor edad mayor cantidad de hijos.
Finalmente, se tiene que hay una relación negativa entre la cantidad de hijos y los años de escolaridad.
Sobredispersión
dispersiontest(modPoisson)##
## Overdispersion test
##
## data: modPoisson
## z = -3.3515, p-value = 0.9996
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 0.8249273Tras realizar la prueba de sobredispersión se obtiene un p valor cercano a 1 que nos lleva a no rechazar la hipótesis nula y concluir que no hay sobredispersión.
Por lo tanto no realizaremos un modelo binomial-negativo.