Punto 1 (30pts) La base de datos baseball.csv 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. Las variables observadas fueron:

Comenzamos leyendo la base de datos baseball

## # A tibble: 6 × 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
## # … with 4 more variables: strikeouts <dbl>, stolenbases <dbl>, errors <dbl>,
## #   name <chr>

Se procede a realizar un analisis descriptivo con la variable respuesta:

Observando el histograma de la variable salario, se puede decir que la mayoria de los jugadores ganan entre 0 a 4500 dolares; sin embargo algunos jugadores ganan incluso hasta mas de 5000 dolares.

Se quiere un modelo que prediga el sueldo del jugador con base a su rendimiento. Divida la muestra en dos partes de forma aleatoria: una para estimar el modelo - entrenamiento - (100 observaciones) y otra para hacer una evaluación del modelo ajustado - validación - (32 observaciones).

Se procede a dividir la muestra en dos partes de forma aleatoria:

set.seed(3542)
entrenamiento<-sample_n(datos, size=100)
validacion<- datos[!(datos$name %in% entrenamiento$name),]

Ajuste un modelo de regresión lineal considerando la variable salary como respuesta y todas las demás como covariables (excepto names, obviamente). Si es necesario, utilice transformaciones. Comente los resultados más relevantes.

Ahora se realiza el ajuste del modelo:

modelo<-lm(salary~.,data = entrenamiento[,-14])
summary(modelo)
## 
## Call:
## lm(formula = salary ~ ., data = entrenamiento[, -14])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1688.6  -577.4  -145.3   572.1  2101.1 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  9.704e+02  7.908e+02   1.227   0.2231  
## batting      5.548e+03  9.583e+03   0.579   0.5641  
## OBP         -5.400e+03  7.170e+03  -0.753   0.4534  
## runs         2.220e+01  1.462e+01   1.519   0.1325  
## hits        -6.531e-01  8.219e+00  -0.079   0.9368  
## dooubles    -1.858e+01  1.867e+01  -0.995   0.3225  
## tripes      -5.684e+01  4.466e+01  -1.273   0.2066  
## homeruns     1.487e+01  2.465e+01   0.603   0.5480  
## RBI          1.570e+01  9.920e+00   1.583   0.1171  
## walks        8.343e-02  1.207e+01   0.007   0.9945  
## strikeouts  -2.680e+00  4.800e+00  -0.558   0.5780  
## stolenbases  1.088e+01  1.101e+01   0.989   0.3255  
## errors      -3.223e+01  1.544e+01  -2.087   0.0398 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 833.9 on 87 degrees of freedom
## Multiple R-squared:  0.554,  Adjusted R-squared:  0.4924 
## F-statistic: 9.004 on 12 and 87 DF,  p-value: 5.784e-11

Como se puede observar, los resultados del modelo nos arroja que ninguna de las variables se muestran significativas para el modelo a excepcion de errors. Además el \(R^{2}\) fue de 55.4%, por lo que se dice que el modelo solamente explica poco mas de la mitad de la variación de la variable respuesta (es poco) y además, el \(R^{2}\) ajustado fue de 49.24%.

Miramos los VIF, donde se aprecia que varios VIF son mayores a 10, lo que indica problemas de multicolinealidad que pueden ser corregidos con la seleccion de variables que se realizara posteriormente.

vif(modelo)
##     batting         OBP        runs        hits    dooubles      tripes 
##   20.006490   16.240929   23.614127   20.693109    4.761134    1.791253 
##    homeruns         RBI       walks  strikeouts stolenbases      errors 
##    8.957574   11.214526   11.318402    3.243120    2.313928    1.181935

Supuestos

Homocedasticidad

###RESIDUOS
res.stud.base = studres(modelo)
modelo.fit.base= modelo$fitted.values
par(mfrow=c(1,2))
plot(modelo.fit.base,res.stud.base, ylab='residuos estudentizados',
     xlab='valores ajustados',main='(a)', col = "BLUE")
abline(h=0,lty=2)
lines(lowess(res.stud.base~modelo.fit.base), col = 2)

plot(modelo.fit.base,abs(res.stud.base), 
     ylab='valor absoluto de los residuos estudentizados',
     xlab='valores ajustados',main='(b)', col = "BLUE")
lines(lowess(abs(res.stud.base)~modelo.fit.base), col = 2)

Con las graficas de los residuos estudentizados vemos que no se nota un patron muy marcado o alguna tendencia, además de que los residuos se distribuyen alrededor de cero. Entonces parece que no hay problemas de heterocedasticidad y lo comprobamos con la prueba de Breusch-Pagan.

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 13.234, df = 12, p-value = 0.3523

Realizando el test de Breusch-pagan obtuvimos un valor p mayor a 0.05 de tal manera que no hay evidencia suficiente para rechazar la hipotesis nula y de esta manera se concluye que se cumple el supuesto de homocedasticidad.

Normalidad

car::qqPlot(modelo,xlab='cuantiles teóricos',ylab='residuos estudentizados', distribution = 'norm')

## [1] 36 77
shapiro.test(modelo$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo$residuals
## W = 0.98613, p-value = 0.3818

Con la grafica de normalidad, se puede observar que gran parte de los datos siguen el patron lineal; por lo cual se puede decir que se cumple el supuesto de normalidad. Para confirmar lo observado en la gráfica, se procede a realizar el test de Shapiro-Wilk, en el cual se obtiene un valor p de 0.3818, el cual es mayor al nivel de significancia 0.05 y por ende se concluye que el modelo cumple el supuesto de normalidad.

Realice un proceso de selección de variables comparando todos los modelos posibles. Proponga al menos tres modelos diferentes. Compare estos modelos con el modelo completo. Comente los resultados más relevantes.

Modelos propuestos

Se calculan todos los modelos posibles con la siguiente función y luego se proponen 3 modelos segun el \(R^2\) ajustado, Cp de Mallows y BIC o SBC.

all.mod=ols_step_all_possible(modelo)

Modelo 1

El mejor modelo segun el \(R^2\) ajustado más alto es aquel que considera las variables runs, dooubles, tripes, RBI, walks y errors.

R2adj.order = order(all.mod$adjr,decreasing = T)
as.data.frame(all.mod)[R2adj.order[1],c(2:8,10)]
##      n                            predictors   rsquare      adjr   predrsq
## 2382 6 runs dooubles tripes RBI walks errors 0.5450555 0.5157042 0.4776754
##            cp      aic      sbc
## 2382 2.738365 1633.052 1653.893

Ajustando el modelo:

mod1=lm(salary~ runs+dooubles+tripes+RBI+walks+errors, data=entrenamiento[,-14])
summary(mod1)
## 
## Call:
## lm(formula = salary ~ runs + dooubles + tripes + RBI + walks + 
##     errors, data = entrenamiento[, -14])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1807.3  -550.5  -123.2   542.7  1997.7 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  624.615    222.937   2.802  0.00618 ** 
## runs          32.512      8.225   3.953  0.00015 ***
## dooubles     -23.024     15.393  -1.496  0.13811    
## tripes       -72.413     41.220  -1.757  0.08225 .  
## RBI           16.722      5.306   3.152  0.00219 ** 
## walks        -10.714      6.295  -1.702  0.09213 .  
## errors       -32.551     14.557  -2.236  0.02774 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 814.5 on 93 degrees of freedom
## Multiple R-squared:  0.5451, Adjusted R-squared:  0.5157 
## F-statistic: 18.57 on 6 and 93 DF,  p-value: 4.412e-14

Entonces vemos que las covariables runs, RBI y errors resultan tener un aporte significativo al modelo, ademas de que el \(R^2\) ajustado fue de 51.57%.

bptest(mod1)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 6.9364, df = 6, p-value = 0.3268
shapiro.test(mod1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.98703, p-value = 0.4393

En cuanto a los supuestos vemos que en las pruebas de Breusch-Pagan y de Shapiro-Wilk se consigue un valor p mayor a 0.05 en ambos casos, por lo que no habrían problemas de heterocedasticidad y ademas, se cumpliria el supuesto de normalidad.

Por último, vemos que los VIF son menores a 10, por lo que no habrían problemas de multicolinealidad.

car::vif(mod1)
##     runs dooubles   tripes      RBI    walks   errors 
## 7.834080 3.390376 1.598954 3.362412 3.227740 1.100735

Modelo 2

El mejor modelo segun el Cp de Mallows mas bajo es aquel que considera las variables runs, RBI, stolenbases y errors.

cp.order= order(all.mod$cp, decreasing=F)
as.data.frame(all.mod)[cp.order[1],c(2:8,10)]
##     n                  predictors  rsquare      adjr   predrsq       cp
## 663 4 runs RBI stolenbases errors 0.532583 0.5129023 0.4820084 1.171162
##          aic      sbc
## 663 1631.756 1647.387

Se ajusta el modelo

mod2=lm(salary~runs+RBI+stolenbases+errors,data=entrenamiento[,-14])
summary(mod2)
## 
## Call:
## lm(formula = salary ~ runs + RBI + stolenbases + errors, data = entrenamiento[, 
##     -14])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1581.81  -571.97   -41.51   559.69  1998.88 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  392.363    199.164   1.970 0.051745 .  
## runs          10.446      6.359   1.643 0.103728    
## RBI           20.193      5.575   3.622 0.000472 ***
## stolenbases   15.751      9.649   1.632 0.105931    
## errors       -28.504     14.236  -2.002 0.048104 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 816.9 on 95 degrees of freedom
## Multiple R-squared:  0.5326, Adjusted R-squared:  0.5129 
## F-statistic: 27.06 on 4 and 95 DF,  p-value: 5.381e-15

Vemos que el \(R^2\) ajustado fue de 51.29% y que RBI y errors tuvieron un aporte significativo en el modelo.

bptest(mod2)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod2
## BP = 7.2784, df = 4, p-value = 0.1219
shapiro.test(mod2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$residuals
## W = 0.98487, p-value = 0.3108

Pasando a los supuestos, vemos que en las pruebas de Breusch-Pagan y de Shapiro-Wilk se consigue un valor p mayor a 0.05 en ambos casos, por lo que no habrían problemas de heterocedasticidad y ademas, se cumpliria el supuesto de normalidad.

Por último, vemos que los VIF son menores a 10, por lo que no habrían problemas de multicolinealidad.

car::vif(mod2)
##        runs         RBI stolenbases      errors 
##    4.655424    3.691237    1.852934    1.046610

Modelo 3

El mejor modelo segun el BIC o SBC mas bajo es aquel que considera las variables runs y RBI.

BIC.order = order(all.mod$sbc,decreasing = F)
as.data.frame(all.mod)[BIC.order[1],c(2:8,10)]
##    n predictors   rsquare     adjr predrsq      cp      aic      sbc
## 38 2   runs RBI 0.5019315 0.491662  0.4704 3.14983 1634.108 1644.529

Se procede a ajustar el modelo

mod3=lm(salary~ runs+RBI, data=entrenamiento[,-14])
summary(mod3)
## 
## Call:
## lm(formula = salary ~ runs + RBI, data = entrenamiento[, -14])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1516.00  -612.02   -83.94   572.47  2175.83 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  290.997    193.495   1.504 0.135857    
## runs          16.503      4.788   3.447 0.000841 ***
## RBI           14.835      4.715   3.147 0.002194 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 834.5 on 97 degrees of freedom
## Multiple R-squared:  0.5019, Adjusted R-squared:  0.4917 
## F-statistic: 48.88 on 2 and 97 DF,  p-value: 2.082e-15

Se observa que todas las covariables tuvieron un aporte significativo en el modelo y ademas, el \(R^2\) ajustado fue de 0.4917.

Se realiza la comprobación de los supuestos

bptest(mod3)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod3
## BP = 3.5006, df = 2, p-value = 0.1737
shapiro.test(mod3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.98313, p-value = 0.2311

Vemos que en las pruebas de Breusch-Pagan y de Shapiro-Wilk se consigue un valor p mayor a 0.05 en ambos casos, por lo que no habrían problemas de heterocedasticidad y ademas, se cumpliria el supuesto de normalidad.

Por último, vemos que los VIF son menores a 10, por lo que no habrían problemas de multicolinealidad.

car::vif(mod3)
##     runs      RBI 
## 2.529422 2.529422

Comparacion

Para comparar los tres modelos propuestos con respecto al modelo completo, podemos usar los \(R^2\) ajustado y los AIC

R2_Ajd<-c(summary(modelo)$adj.r.squared,summary(mod1)$adj.r.squared,summary(mod2)$adj.r.squared,summary(mod3)$adj.r.squared)
AIC<-c(AIC(modelo),AIC(mod1),AIC(mod2),AIC(mod3))
Tabla<- cbind(R2_Ajd,AIC)
row.names(Tabla)<- c("Modelo completo","Modelo R2 adj","Modelo CP","Modelo BIC")
Tabla
##                    R2_Ajd      AIC
## Modelo completo 0.4924460 1643.073
## Modelo R2 adj   0.5157042 1633.052
## Modelo CP       0.5129023 1631.756
## Modelo BIC      0.4916620 1634.108

Entonces se observa en la tabla que el modelo con el \(R^2\) ajustado más alto fue claramente el primer modelo propuesto, ya que fue elegido mediante este criterio del \(R^2\) ajustado más alto. Y el mejor modelo segun el AIC fue el modelo propuesto bajo el criterio del CP. Aunque los 3 modelos propuestos presentan valores o indicadores similares, se podría decir que de alguna manera, el mejor podria ser el Modelo CP por el hecho de tener una cantidad no muy grande de variables, tener el AIC mas bajo y tener un \(R^2\) ajustado muy similar al del mejor modelo con el \(R^2\) ajustado mas alto.

A partir de los modelos propuestos (y el modelo completo) haga predicciones para los jugadores de la sub-muestra de validación. Compare los modelos usando el error cuadrado medio de predicción.

Teniendo el modelo completo y los 3 modelos propuestos, se procede a realizar predicciones a partir de la muestra de validacion para comparar los modelos usando el error cuadrado medio de predicción donde:

Error cuadrado medio de prediccion = \(\frac{\sum_{i=1}^n(y_i-\hat{y_i})^2}{n}\)

valores=rbind(mean((validacion$salary - predict(modelo, newdata=validacion)) ^ 2), 
mean((validacion$salary - predict(mod1, newdata=validacion)) ^ 2), 
mean((validacion$salary - predict(mod2, newdata=validacion)) ^ 2), 
mean((validacion$salary - predict(mod3, newdata=validacion)) ^ 2)) 
colnames(valores)<- "Error cuadrado medio de predicción"
row.names(valores)<- c("Modelo completo","Modelo R2 adj","Modelo CP","Modelo BIC")
valores
##                 Error cuadrado medio de predicción
## Modelo completo                           997087.1
## Modelo R2 adj                            1017840.8
## Modelo CP                                 923420.2
## Modelo BIC                                843673.5

Vemos que el modelo BIC es el que menor error cuadrado medio de predicción tiene, por lo que sería el modelo a utilizar para mejores predicciones, además de que es un modelo sencillo al tener solo 2 variables.

Punto 2 (30pts) Los datos fertility (de la librera Countr) contienen una muestra de 1243 mujeres casadas (de 44 a nos o mas) de Alemania en 1985. Se quiere determina que factores influyen en el numero de hijos que tienen estas mujeres.

Comenzamos leyendo la base de datos fertility

data("fertility", package = "Countr")
fertility2 <- fertility[ , -c(2,4)]
head(fertility2)
##   children years_school university religion year_birth rural age_marriage
## 1        2            8         no Catholic         42   yes           20
## 2        3            8         no Catholic         55   yes           21
## 3        2            8         no Catholic         51   yes           24
## 4        4            8         no Catholic         54    no           26
## 5        2            8         no Catholic         46   yes           22
## 6        2            8         no Catholic         41    no           18

Se realizan una estadistica descriptiva de la variable respuesta niños:

Observando el histograma de la variable niños, se dice que la mayoria de las mujeres en Alemania tienen entre 1 y 2 hijos; sin embargo tambien se observa que un gran grupo de mujeres tienen entre 3 a 6 hijos, y solo unas pocas tienen mas de 6 hijos.

Ajuste un modelo Poisson para este conjunto de datos usando children como covariable y las otras como covariables. Comente los resultados más relevantes. Interprete los coeficientes que considere más relevantes.

modelopois<-glm(children~.,data = fertility2,family = poisson(link = "log"))
summary(modelopois)
## 
## Call:
## glm(formula = children ~ ., family = poisson(link = "log"), data = fertility2)
## 
## 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: 5

Observando el resumen del modelo de regresion poisson, se dice que solamente el intercepto y las covariables religion musulman, otra religion, rural y edad en que contrajo matrimonio son significativas para el modelo. Aunque la covariable años de educación no es significativa para el modelo, se observa que coeficiente es negativo, por lo que dependiendo de los años que la mujer estudió, asi mismo decrece la cantidad de niños que tienda a tener; La edad en la que contrajo matrimonio es significativa como lo habiamos mecionado anteriormente y tambien presenta un coeficiente negativo, por lo que entre mas edad tenga la mujer cuando contraiga matrimonio, menos será la cantidad de hijos que tendra. Por otro lado los demas coeficientes son positivos, sin embargo son bajos (\(coef<1\)), por lo tanto la cantidad que sumen para la cantidad de hijos en las mujeres va a ser muy lenta.

exp(coef(modelopois))
##        (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.9695245

Elevando los coeficientes a la exponencial, se observa que la covariable universidad se puede concluir que la si la mujer si va a la univerdad, el valor esperado del numero de hijos aumentará en 1.14, con respecto al intercepto se dice que el valor esperado del numero de hijos sin tener en cuenta ningun factor es 5.31.

Evalúe si los datos presentan sobredispersión. Si es necesario, ajuste un modelo binomial-negativo para corregir este problema, ¿presenta este modelo un mejor ajuste que el modelo Poisson?

dispersiontest(modelopois)
## 
##  Overdispersion test
## 
## data:  modelopois
## z = -3.3515, p-value = 0.9996
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##  0.8249273

Observando la prueba de sobredispersion, nos arroja con valor de dispersion 0.8249273, el cual es \(<1\) y por ende se concluye que el modelo de regresion poisson no presenta sobredispersion; de igual manera se confirma al observar el valor p obtenido de la prueba, el cual es mayor que 0.05 y por ende se dice que no hay evidencia suficiente para rechazar la hipotesis nula la cual es que el modelo no presenta sobredispersion.