Librerías

library(readr)
library(lmtest)
library(MASS)
library(dplyr)
library(car)
library(olsrr)
library(Countr)
library(lmridge)
library(glmnet)
library(ggplot2)
library(caret)
library(pscl)
library(MLmetrics)
library(rcompanion)
library(readxl)

Punto 1

La base de datos baseball.csv 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. Las variables observadas fueron: • salary: salario (en miles de dólares) - variable respuesta -.

• batting: Promedio de bateo.

• OBP: Porcentaje en base - la frecuencia con la que un bateador llega a la base.

• runs: Número de carreras.

• hits: Números de aciertos al bate.

• doubles: Número de dobles.

• tripes: Número de triples.

• homeruns: Número de jonrones.

• RBI: Número de carreras impulsadas.

• walks: Número de base por bolas.

• strikeouts: Número de ponches (strike-outs).

• stolenbases: Número de bases robadas.

• errors: Número de errores.

• names: Nombre del jugador.

Datos baseball

datosbaseball <- read_csv("C:/Users/ronal/Desktop/MODELOS 2/Parcial final/baseball[1].csv")
datosbaseball<-datosbaseball[,-14]
attach(datosbaseball)

División de los datos en Entrenamiento(DatosE) y Validación (DatosV)

set.seed(1)
aleatorios=sample(1:nrow(datosbaseball),100,replace=FALSE)
DatosE=datosbaseball[aleatorios,]
nrow(DatosE)
## [1] 100
DatosV=datosbaseball[-aleatorios,]
nrow(DatosV)
## [1] 32

Con la submuestra de entrenamiento se realizaran los análisis y ajuste de los modelos. Donde en primer lugar se realiza un análisis descriptivo de los datos:

attach(DatosE)
## The following objects are masked from datosbaseball:
## 
##     batting, dooubles, errors, hits, homeruns, OBP, RBI, runs, salary,
##     stolenbases, strikeouts, tripes, walks
MitadDatosE=data.frame(DatosE$salary,DatosE$batting,DatosE$OBP,DatosE$runs,DatosE$hits,DatosE$dooubles,DatosE$tripes)
pairs(MitadDatosE)

cor(MitadDatosE)
##                 DatosE.salary DatosE.batting DatosE.OBP DatosE.runs DatosE.hits
## DatosE.salary       1.0000000      0.4136783  0.3777788   0.6859982   0.6452561
## DatosE.batting      0.4136783      1.0000000  0.7319815   0.5782514   0.7214408
## DatosE.OBP          0.3777788      0.7319815  1.0000000   0.6035893   0.4849843
## DatosE.runs         0.6859982      0.5782514  0.6035893   1.0000000   0.8736676
## DatosE.hits         0.6452561      0.7214408  0.4849843   0.8736676   1.0000000
## DatosE.dooubles     0.5872832      0.5512843  0.3446466   0.7415263   0.8321962
## DatosE.tripes       0.1580018      0.3059972  0.1932866   0.3815204   0.3432410
##                 DatosE.dooubles DatosE.tripes
## DatosE.salary         0.5872832     0.1580018
## DatosE.batting        0.5512843     0.3059972
## DatosE.OBP            0.3446466     0.1932866
## DatosE.runs           0.7415263     0.3815204
## DatosE.hits           0.8321962     0.3432410
## DatosE.dooubles       1.0000000     0.2106851
## DatosE.tripes         0.2106851     1.0000000
otramitaddatosE=data.frame(DatosE$salary,DatosE$homeruns,DatosE$RBI,DatosE$walks,DatosE$strikeouts,DatosE$stolenbases,DatosE$errors)
pairs(otramitaddatosE)

cor(otramitaddatosE)
##                    DatosE.salary DatosE.homeruns  DatosE.RBI DatosE.walks
## DatosE.salary         1.00000000     0.636419592  0.69939929    0.5040916
## DatosE.homeruns       0.63641959     1.000000000  0.89760064    0.4688041
## DatosE.RBI            0.69939929     0.897600640  1.00000000    0.5234674
## DatosE.walks          0.50409163     0.468804143  0.52346739    1.0000000
## DatosE.strikeouts     0.46974803     0.721407109  0.62926060    0.5337602
## DatosE.stolenbases    0.20353346     0.002121256 -0.00518212    0.3297183
## DatosE.errors         0.09347357     0.175628075  0.23053522    0.1199691
##                    DatosE.strikeouts DatosE.stolenbases DatosE.errors
## DatosE.salary              0.4697480        0.203533456    0.09347357
## DatosE.homeruns            0.7214071        0.002121256    0.17562807
## DatosE.RBI                 0.6292606       -0.005182120    0.23053522
## DatosE.walks               0.5337602        0.329718330    0.11996909
## DatosE.strikeouts          1.0000000        0.152056566    0.13304872
## DatosE.stolenbases         0.1520566        1.000000000    0.09574443
## DatosE.errors              0.1330487        0.095744428    1.00000000

Ajuste de modelo salario (Y) y todas como covariables(X)

\[Salary_i =\beta_0+batting_i\beta_1+OBP_i\beta_2+runs_i\beta_3+hits_i\beta_4+dooubles_i\beta_5\]

\[ +tripes_i\beta_6+homeruns_i\beta_7+RBI_i\beta_8+walks_i\beta_9+strikeouts_i\beta_10+stolenbases_i\beta_11+errors_i\beta_12+\epsilon_i \]

Modelocom=lm(salary~.,data = DatosE)
summary(Modelocom)
## 
## Call:
## lm(formula = salary ~ ., data = DatosE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1921.64  -531.25   -60.05   491.41  2675.65 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -142.8460  1281.2624  -0.111    0.911
## batting      5812.2171 11636.7768   0.499    0.619
## OBP         -2770.8349  9890.1776  -0.280    0.780
## runs            4.3917    15.3532   0.286    0.776
## hits            1.0016     8.5998   0.116    0.908
## dooubles        3.7813    17.2939   0.219    0.827
## tripes        -15.4985    49.7727  -0.311    0.756
## homeruns       17.3848    28.1126   0.618    0.538
## RBI            15.7052    11.3231   1.387    0.169
## walks           6.9320    14.6473   0.473    0.637
## strikeouts     -0.1812     4.9638  -0.036    0.971
## stolenbases    12.8036    12.0740   1.060    0.292
## errors        -21.4075    16.7176  -1.281    0.204
## 
## Residual standard error: 882.4 on 87 degrees of freedom
## Multiple R-squared:  0.563,  Adjusted R-squared:  0.5027 
## F-statistic:  9.34 on 12 and 87 DF,  p-value: 2.57e-11

Observando los resultados del modelo se tiene que su R cuadrado ajustado es de 0.5282 siendo este la variabilidad explicada por todas las covariables, donde la unica covariable significativa es RBI, estos resultados siendo poco informativos, se realiza la validación de los supuestos.

Análisis de los supuestos modelo completo

par(mfrow=c(1,3))
res.stud.baseballcom = studres(Modelocom)
mod.fit.baseballcom  = Modelocom$fitted.values

plot(mod.fit.baseballcom,res.stud.baseballcom , ylab='residuos estudentizados',
     xlab='valores ajustados',main ="Linealidad")
abline(h=0,lty=2)
lines(lowess(res.stud.baseballcom ~mod.fit.baseballcom), col = "red")
plot(mod.fit.baseballcom,abs(res.stud.baseballcom ), 
     ylab='|residuos estudentizados|',
     xlab='valores ajustados',main='Homocedasticidad')
lines(lowess(abs(res.stud.baseballcom )~mod.fit.baseballcom), col = "red")
car::qqPlot(Modelocom,xlab='cuantiles teoricos',ylab='residuos estudentizados',main="Grafico de normalidad",
            distribution = 'norm',col="red",id=FALSE)

nortest::ad.test(res.stud.baseballcom)
## 
##  Anderson-Darling normality test
## 
## data:  res.stud.baseballcom
## A = 0.36613, p-value = 0.4279
bptest(Modelocom,~batting+OBP+runs+hits+dooubles+tripes+homeruns+RBI+walks+strikeouts+stolenbases+errors+I(batting^2)+I(OBP^2)+I(runs^2)+I(hits^2)+I(dooubles^2)+I(tripes^2)+I(homeruns^2)+I(RBI^2)+I(walks^2)+I(strikeouts^2)+I(stolenbases^2)+I(errors^2),data = DatosE)
## 
##  studentized Breusch-Pagan test
## 
## data:  Modelocom
## BP = 38.587, df = 24, p-value = 0.03016

Dados los resultados obtenidos tanto el análisis grafico como la prueba de normalidad de Anderson Darling con un valor p mayor que el nivel de significancia por lo tanto no hay evidencia suficiente para rechazar de que los residuos se distribuyen normal.

El supuesto de linealidad la recta presenta una pequeña curvatura la cual no parece afectar el supuesto de linealidad, sin embargo observando la grafica de Homocedasticidad parece tener una curvatura que rechazaria que los residuos se distribuyen homocedasticamente para confirmar con la prueba formal de Breusch-Pagan se tiene que el valor p es menor que el nivel de significancia (0.05).Dada la informacion suministrada por los datos se rechaza la hipotesis de homocedasticidad. Por esto se realiza una transformacion.

modcomtrans=lm(sqrt(salary)~.,data = DatosE)
summary(modcomtrans)
## 
## Call:
## lm(formula = sqrt(salary) ~ ., data = DatosE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.9117  -7.1513  -0.2382   6.3328  22.3891 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  23.61431   14.29238   1.652    0.102
## batting      78.37534  129.80731   0.604    0.548
## OBP         -59.98713  110.32413  -0.544    0.588
## runs          0.10879    0.17126   0.635    0.527
## hits          0.03721    0.09593   0.388    0.699
## dooubles     -0.03730    0.19291  -0.193    0.847
## tripes       -0.36988    0.55521  -0.666    0.507
## homeruns      0.25215    0.31359   0.804    0.424
## RBI           0.12745    0.12631   1.009    0.316
## walks         0.08218    0.16339   0.503    0.616
## strikeouts   -0.01266    0.05537  -0.229    0.820
## stolenbases   0.10878    0.13468   0.808    0.421
## errors       -0.29693    0.18648  -1.592    0.115
## 
## Residual standard error: 9.843 on 87 degrees of freedom
## Multiple R-squared:  0.5817, Adjusted R-squared:  0.524 
## F-statistic: 10.08 on 12 and 87 DF,  p-value: 4.471e-12
par(mfrow=c(1,3))
res.stud.baseballcomtrans = studres(modcomtrans)
mod.fit.baseballcomtrans  = modcomtrans$fitted.values

plot(mod.fit.baseballcomtrans,res.stud.baseballcomtrans , ylab='residuos estudentizados',
     xlab='valores ajustados',main ="Linealidad")
abline(h=0,lty=2)
lines(lowess(res.stud.baseballcomtrans ~mod.fit.baseballcomtrans), col = "red")
plot(mod.fit.baseballcomtrans,abs(res.stud.baseballcomtrans ), 
     ylab='|residuos estudentizados|',
     xlab='valores ajustados',main='Homocedasticidad')
lines(lowess(abs(res.stud.baseballcomtrans )~mod.fit.baseballcomtrans), col = "red")
car::qqPlot(modcomtrans,xlab='cuantiles teoricos',ylab='residuos estudentizados',main="Grafico de normalidad",
            distribution = 'norm',col="red",id=FALSE)

nortest::ad.test(res.stud.baseballcomtrans)
## 
##  Anderson-Darling normality test
## 
## data:  res.stud.baseballcomtrans
## A = 0.27151, p-value = 0.6657
bptest(modcomtrans,~batting+OBP+runs+hits+dooubles+tripes+homeruns+RBI+walks+strikeouts+stolenbases+errors+I(batting^2)+I(OBP^2)+I(runs^2)+I(hits^2)+I(dooubles^2)+I(tripes^2)+I(homeruns^2)+I(RBI^2)+I(walks^2)+I(strikeouts^2)+I(stolenbases^2)+I(errors^2),data = DatosE)
## 
##  studentized Breusch-Pagan test
## 
## data:  modcomtrans
## BP = 28.406, df = 24, p-value = 0.2433

Con la tranformacion de raiz cuadrada sobre la variable de respuesta se corrige el supuesto de homocedasticidad, sin embargo dado a que se tienen muchas covariables (12) puede presentar problemas de multicolinealidad para esto se realiza en análisis del factor de inflacion de varianza:

car::vif(modcomtrans)
##     batting         OBP        runs        hits    dooubles      tripes 
##   17.717078   16.843580   22.153333   18.481402    4.032356    1.664189 
##    homeruns         RBI       walks  strikeouts stolenbases      errors 
##   11.453236   14.070871   14.593344    3.365905    2.556619    1.155366

Según este indicador muestra que existen problemas de multicolinealidad entre las covariables observando valores mayores a 10 los cuales indican esto puede afectar a los calculos y las predicciones del modelo, por esto se corrige mediante la regresion de Ridge

Regresion de Ridge

plot(lm.ridge(sqrt(salary)~.,data=DatosE,lambda=seq(0,3,0.01)))
title(main="Ridge Plot")

lamda<-lm.ridge(sqrt(salary)~.,data=DatosE)
kestimado<-lamda[["kHKB"]]
kestimado
## [1] 18.43537
modri<-lmridge(sqrt(salary)~.,data=DatosE,K=1.2)
modri
## Call:
## lmridge.default(formula = sqrt(salary) ~ ., data = DatosE, K = 1.2)
## 
##   Intercept     batting         OBP        runs        hits    dooubles 
##    16.69987    25.62741    11.08284     0.05645     0.03330     0.11190 
##      tripes    homeruns         RBI       walks  strikeouts stolenbases 
##    -0.01522     0.15794     0.06141     0.03665     0.02648     0.04237 
##      errors 
##    -0.06706
summary(modri)
## 
## Call:
## lmridge.default(formula = sqrt(salary) ~ ., data = DatosE, K = 1.2)
## 
## 
## Coefficients: for Ridge parameter K= 1.2 
##               Estimate Estimate (Sc) StdErr (Sc) t-value (Sc) Pr(>|t|)    
## Intercept      16.6999    -4984.9482    482.9108     -10.3227   <2e-16 ***
## batting        25.6274        8.1795      3.3740       2.4242   0.0172 *  
## OBP            11.0828        4.0581      3.4636       1.1716   0.2443    
## runs            0.0564       15.2690      2.2075       6.9169   <2e-16 ***
## hits            0.0333       14.6890      2.7330       5.3747   <2e-16 ***
## dooubles        0.1119       11.4649      3.5382       3.2403   0.0016 ** 
## tripes         -0.0152       -0.3480      4.3008      -0.0809   0.9357    
## homeruns        0.1579       16.7769      3.2125       5.2223   <2e-16 ***
## RBI             0.0614       17.9512      2.7510       6.5253   <2e-16 ***
## walks           0.0366        8.4336      3.4932       2.4143   0.0177 *  
## strikeouts      0.0265        8.6363      3.7036       2.3319   0.0218 *  
## stolenbases     0.0424        4.9506      4.1804       1.1843   0.2393    
## errors         -0.0671       -3.8046      4.4447      -0.8560   0.3942    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Ridge Summary
##         R2     adj-R2   DF ridge          F        AIC        BIC 
##   0.368200   0.289200   3.685900   9.983204 460.059559 930.178974 
## Ridge minimum MSE= 2587.257 at K= 1.2 
## P-value for F-test ( 3.6859 , 94.50179 ) = 1.772519e-06 
## -------------------------------------------------------------------
vif(modri)
##       batting    OBP   runs    hits dooubles  tripes homeruns     RBI   walks
## k=1.2 0.11634 0.1226 0.0498 0.07633  0.12794 0.18903  0.10547 0.07734 0.12471
##       strikeouts stolenbases  errors
## k=1.2    0.14018      0.1786 0.20189

En el gráfico de Ridge observamos la tendencia de estabilización de los coeficientes del modelo a partir del valor de \(\lambda > 1\) aproximadamente. La estimación que nos proporcionan los procedimientos HKB sobre el valor de la constante \(\lambda\) es de 0,6231. Despues de ajustar la regrsion de Ridge observamos cómo los coeficientes que más varían son los relativos a las variables

observamos cómo los coeficientes que más varían son los relativos a las variables batting y OBP , que son las más correlacionadas entre sí (correlación simple de 0.73198145 y correlación parcial de 0.22364)

Como seleccion de variables dada la regresion de ridge se tiene como modelo:

\[Salary_i =\beta_0+runs_i\beta_1+hits_i\beta_2+dooubles_i\beta_3+homeruns_i\beta_4+RBIs_i\beta_5\]

ridgemod = lmridge(salary ~ runs+hits+dooubles+homeruns+RBI,data = DatosE) 
summary(ridgemod)
## 
## Call:
## lmridge.default(formula = salary ~ runs + hits + dooubles + homeruns + 
##     RBI, data = DatosE)
## 
## 
## Coefficients: for Ridge parameter K= 0 
##              Estimate Estimate (Sc) StdErr (Sc) t-value (Sc) Pr(>|t|)  
## Intercept  2.8482e+02   -5.4306e+05  3.5849e+05      -1.5148   0.1332  
## runs       1.5097e+01    4.0837e+03  1.9895e+03       2.0527   0.0429 *
## hits       1.4650e+00    6.4623e+02  2.5489e+03       0.2535   0.8004  
## dooubles   2.1979e+00    2.2519e+02  1.6071e+03       0.1401   0.8889  
## homeruns   1.3728e+01    1.4582e+03  2.3660e+03       0.6163   0.5392  
## RBI        1.2427e+01    3.6327e+03  2.9033e+03       1.2512   0.2140  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Ridge Summary
##         R2     adj-R2   DF ridge          F        AIC        BIC 
##    0.54770    0.52870    4.99997   23.01173 1356.02720 1829.56999 
## Ridge minimum MSE= 27064469 at K= 0 
## P-value for F-test ( 4.99997 , 95.00005 ) = 4.549457e-15 
## -------------------------------------------------------------------

Con un R cuadrado ajustado de 0.5237 explicando la variabilidad del modelo y como unica variable significativa runs es una propuesta de modelo con el cumplimiento de todos los supuestos.

Selección de variables

#seleccion de variables 
selecmod=ols_step_best_subset(modcomtrans)
selecmod
##                                        Best Subsets Regression                                       
## -----------------------------------------------------------------------------------------------------
## Model Index    Predictors
## -----------------------------------------------------------------------------------------------------
##      1         runs                                                                                   
##      2         runs RBI                                                                               
##      3         runs RBI errors                                                                        
##      4         runs tripes RBI errors                                                                 
##      5         runs hits homeruns RBI errors                                                          
##      6         runs hits tripes RBI stolenbases errors                                                
##      7         runs hits tripes homeruns RBI stolenbases errors                                       
##      8         batting runs hits tripes homeruns RBI stolenbases errors                               
##      9         batting OBP runs tripes homeruns RBI walks stolenbases errors                          
##     10         batting OBP runs hits tripes homeruns RBI walks stolenbases errors                     
##     11         batting OBP runs hits tripes homeruns RBI walks strikeouts stolenbases errors          
##     12         batting OBP runs hits dooubles tripes homeruns RBI walks strikeouts stolenbases errors 
## -----------------------------------------------------------------------------------------------------
## 
##                                                       Subsets Regression Summary                                                      
## --------------------------------------------------------------------------------------------------------------------------------------
##                        Adj.        Pred                                                                                                
## Model    R-Square    R-Square    R-Square     C(p)        AIC         SBIC        SBC          MSEP         FPE        HSP       APC  
## --------------------------------------------------------------------------------------------------------------------------------------
##   1        0.4927      0.4875      0.4717     9.5098    752.5014    468.4993    760.3169    10430.7516    106.3932    1.0753    0.5280 
##   2        0.5593      0.5503      0.5327    -2.3456    740.4235    457.1655    750.8441     9155.3833     94.2900    0.9536    0.4679 
##   3        0.5676      0.5541      0.5298    -2.0707    740.5234    457.6065    753.5492     9077.6245     94.3866    0.9553    0.4684 
##   4        0.5709      0.5529      0.5247    -0.7569    741.7575    459.1460    757.3885     9104.1967     95.5625    0.9682    0.4742 
##   5        0.5736      0.5509      0.5121     0.6860    743.1312    460.8418    761.3674     9144.6430     96.8903    0.9829    0.4808 
##   6        0.5767      0.5494      0.5061     2.0426    744.4031    462.4799    765.2445     9176.9818     98.1389    0.9969    0.4870 
##   7        0.5790      0.5470      0.4929     3.5640    745.8580    464.3023    769.3045     9227.3895     99.5884    1.0133    0.4942 
##   8        0.5799      0.5429       0.484     5.3832    747.6513    466.4208    773.7030     9310.6515    101.4052    1.0337    0.5032 
##   9        0.5808      0.5389      0.4758     7.1823    749.4211    468.5306    778.0780     9393.6177    103.2347    1.0545    0.5123 
##  10        0.5813      0.5343      0.4606     9.0822    751.3062    470.7380    782.5683     9489.4545    105.2227    1.0772    0.5222 
##  11        0.5815      0.5292      0.4421    11.0374    753.2548    472.9971    787.1220     9593.5920    107.3219    1.1014    0.5326 
##  12        0.5817      0.5240      0.4276    13.0000    755.2118    475.2645    791.6842     9700.9759    109.4777    1.1265    0.5433 
## --------------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria 
##  SBIC: Sawa's Bayesian Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
##  MSEP: Estimated error of prediction, assuming multivariate normality 
##  FPE: Final Prediction Error 
##  HSP: Hocking's Sp 
##  APC: Amemiya Prediction Criteria

4 runs + tripes + RBI + errors

5 runs + hits + homeruns + RBI + errors

Los modelos fueron seleccionado a partir de los criterios AIC y SBIC donde estos presentaron los menores valores

mod4<-lm(salary~runs + tripes + RBI + errors,data = DatosE)
mod5<-lm(salary~runs + hits + tripes + RBI + errors,data = DatosE)

summary(mod4)
## 
## Call:
## lm(formula = salary ~ runs + tripes + RBI + errors, data = DatosE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1794.79  -535.59   -38.81   472.64  2496.41 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  350.180    221.440   1.581 0.117116    
## runs          18.453      5.418   3.406 0.000969 ***
## tripes       -23.725     41.980  -0.565 0.573304    
## RBI           17.864      4.697   3.803 0.000253 ***
## errors       -18.186     15.472  -1.175 0.242781    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 853.3 on 95 degrees of freedom
## Multiple R-squared:  0.5538, Adjusted R-squared:  0.535 
## F-statistic: 29.47 on 4 and 95 DF,  p-value: 6.182e-16
summary(mod5)
## 
## Call:
## lm(formula = salary ~ runs + hits + tripes + RBI + errors, data = DatosE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1828.51  -515.43   -60.55   485.23  2511.95 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  319.637    244.503   1.307 0.194303    
## runs          17.106      7.044   2.428 0.017065 *  
## hits           1.316      4.367   0.301 0.763814    
## tripes       -25.419     42.555  -0.597 0.551736    
## RBI           17.334      5.037   3.441 0.000866 ***
## errors       -19.080     15.828  -1.205 0.231036    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 857.4 on 94 degrees of freedom
## Multiple R-squared:  0.5542, Adjusted R-squared:  0.5305 
## F-statistic: 23.37 on 5 and 94 DF,  p-value: 3.455e-15
car::vif(mod4)
##     runs   tripes      RBI   errors 
## 2.950523 1.265979 2.589548 1.058274
car::vif(mod5)
##     runs     hits   tripes      RBI   errors 
## 4.938620 5.046987 1.288451 2.949387 1.096846
Modelo N° covariables N° covariables significativas \[ R^2 \] \[R^2_{adj}\] AIC
Completo tranformado 12 1 0.5817 0.524 1658.57
Ridge 5 1 0.5437 0.5237 1356.02720
Modelo 4 (selección) 4 2 0.5709 0.5529 741.7575
Modelo 5 (selección) 5 2 0.5736 0.5509 743.1312

Dado a los ajustes realizados en todos los modelos los cuales cumplen con los supuestos en cada uno se observa que el R cuadrado no varía mucho entre ellos, pero los modelos con el R cuadrado adj más altos fueron los modelos 4 y 5 que también cuentan con el menor AIC por lo tanto por principio de parsimonia y los resultado obtenidos fueron los modelos 4 y 5 los que mejor explicar la variable de respuesta salario para los jugadores de baseball con respecto a las variables de evaluación de rendimiento, con variables significativas en el salario son las runs y el RBI siendo estas 2 las que más influyen en el salario de los jugadores.

Predicciones

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

attach(DatosV)
#Predicción Modelo completo tranformado
x.test1=model.matrix(sqrt(salary)~.,data=DatosV)
pred1=predict(modcomtrans,newx = x.test1)[1:32]
RMSE1=RMSE(pred1,sqrt(DatosV$salary))
RMSE1

#Predicción Modelo Ridge 
x.test2=model.matrix(sqrt(salary)~.,data=DatosV)
pred2=predict(modri , s=1,newx = x.test2)[1:32]
RMSE2=RMSE(pred2,sqrt(DatosV$salary))
RMSE2

#Prediccion Modelo 4
x.test3=model.matrix(sqrt(salary)~runs + tripes + RBI + errors,data=DatosV)
pred3=predict(mod4,newx = x.test3)[1:32]
RMSE3=RMSE(pred3,sqrt(DatosV$salary))
RMSE3

#Prediccion Modelo 5
x.test4=model.matrix(sqrt(salary)~runs + hits + tripes + RBI + errors,data=DatosV)
pred4=predict(mod5,newx = x.test4)[1:32]
RMSE4=RMSE(pred4,sqrt(DatosV$salary))
RMSE4
Modelo RMSEP
Completo transformado 17.86306
Ridge 10.18142
Modelo 4 (selección) 17.25489
Modelo 5(selección) 17.7434

EL modelo que obtuvo al error cuadrático medio de la predicción (RMSEP) con menor valor, siendo así el que mejor predice el salario de los jugadores para la submuestra de validación dado su rendimiento (covaraibles) es el modelo de Ridge, en comparación con los demás modelos. Por lo tanto si la finalidad del estudio es determinar el salario de los jugadores el modelo que mejor predicciones tendrá es el de Ridge.

Punto 2

Datos de fertilidad

Datos=fertility
attach(Datos)

Datos=Datos %>% select(-"german",-"voc_train")

La base de dator fertility (de la librería Countr) de R, contiene una muestra de 1243 mujeres casadas (de 44 años o más) de Alemania en 1985. Se quiere determinar que factores influyeron en el número de hijos que tuvieron estas mujeres. La data cuenta con nueve variables, pero se trabajará con las siguientes:

  • children: número de hijos (variable de respuesta).
  • years_school: años de escolaridad (educación media).
  • university: educación universitaria (si/no).
  • religion: religión que profesa (católica, protestante, musulmana u otras).
  • rural: su vivienda se encuentra en un área rural (si/no).
  • year_birth: año de nacimiento (2 últimos dígitos).
  • age_marriage: edad al momento de contraer matrimonio.

A continuación se realiza un análisis descriptivo y exploratorio de las variables.

summary(Datos)
##     children       years_school    university       religion     year_birth   
##  Min.   : 0.000   Min.   : 8.000   no :1207   Catholic  :130   Min.   :40.00  
##  1st Qu.: 1.000   1st Qu.: 9.000   yes:  36   Muslim    :502   1st Qu.:45.00  
##  Median : 2.000   Median : 9.000              Other     : 75   Median :50.00  
##  Mean   : 2.384   Mean   : 9.097              Protestant:536   Mean   :51.99  
##  3rd Qu.: 3.000   3rd Qu.: 9.000                               3rd Qu.:58.00  
##  Max.   :11.000   Max.   :13.000                               Max.   :83.00  
##  rural      age_marriage  
##  no :613   Min.   :17.00  
##  yes:630   1st Qu.:21.00  
##            Median :23.00  
##            Mean   :23.11  
##            3rd Qu.:25.00  
##            Max.   :30.00

Del resúmen descriptivo, vale la pena resaltar el hecho de que un poco menos del 3% de las mujeres manifestó haber accedido a la educación superior. Además, el 50.68% de la población encuestada, vivia en una zona rural. También, llama la atención que el 75% de las mujeres se casó con máximo 25 años y no tuvo más de tres hijos.

El siguiente gráfico relaciona todas la variables entre si.

pairs(Datos,col="purple")

Inicialmente, se observa que las mujeres que alcanzaron un mayor nivel de escolaridad, tuvieron menos de seis hijos, aunque la frecuencia de aquellas que tuvieron más de esos, no llega al 25%. Además, las mujeres que fueron a la Universidad también se reprodujeron en menor medida que aquellas que no lo hicieron. adicionalmente, no se observa ningún tipo de relación entre el número de hijos que tuvieron las mujeres con su edad de matrimonio, su año de nacimiento o si la zona en que vivian era rural o no. En cuanto a la religión, se ve que en todas hubo tendencia a tener mínimo cuatro hijos.

Ajuste del Modelo

A continuación, se ajusta un modelo Poisson para la base de datos mencionada anteriormente, esto para posteriormente comentar los resultados más relevantes de este modelo, incluyendo la interpretación de sus coeficientes.

Se toma el siguiente modelo: \[ ln[E(children_i)]= \beta_0 + years\_school_i \beta_1 + university_i \beta_2 + religionMuslim_i \beta_3 + religionOther_i \beta_4 \] \[ + religionProtestan_i \beta_5 + religionCatho-i \beta_6 + rural_i \beta_7 + year\_birth_i \beta_8 + age\_marriage_i \beta_9 \]

Para este modelo, rural y university se usarán como variables dicotómicas, donde el “No.=0” será tendrá como referencia, y las variables religionMuslim, religionOther, religionProtestant y religionCatho son variables dicotómicas donde son 0 si la mujer no pertenece a la religión y 1 si pertenece (religionCatho es tomado como referencia).

rural=ifelse(rural=="yes",1,0)
university=ifelse(university=="yes",1,0)

Ajustamos el modelo anterior en R, obteniendo el siguiente resumen:

Modelo=glm(children~years_school+university+religion+rural+year_birth+age_marriage,
           data=Datos,family=poisson)
summary(Modelo)
## 
## Call:
## glm(formula = children ~ years_school + university + religion + 
##     rural + year_birth + age_marriage, family = poisson, data = Datos)
## 
## 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    
## ruralyes            0.078748   0.037221   2.116   0.0344 *  
## year_birth          0.002884   0.002292   1.258   0.2083    
## 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

Al tener un modelo Poisson, es conveniente realizar la transformación \(e^{ln[E(y_i)]}\) para poder interpretar de manera sencilla los coeficientes.

Coeficientes=Modelo$coef
IntCoef=exp(Coeficientes); IntCoef
##        (Intercept)       years_school      universityyes     religionMuslim 
##          5.3138849          0.9564231          1.1446679          1.1620863 
##      religionOther religionProtestant           ruralyes         year_birth 
##          1.8621009          1.0114505          1.0819322          1.0028879 
##       age_marriage 
##          0.9695245

De acuerdo a los coeficientes anteriores, tenemos las siguientes interpretaciones:

  • Cuando las demás covariables se mantiene constantes, el número esperados de hijos por cada mujer cambia por un factor de 0.9564231 por cada año adicional de educación (es decir, una disminución aproximada del (1-0.96) x 100% = 4% en el número esperado de hijos asociados de cada mujer con cada año adicional de educación).

  • Cuando las demás covariables se mantiene constantes, el número esperados de hijos de las mujeres con educación universitaria es 1.1446679 mayor que el numero esperado de hijos de las mujeres que no tienen educación universitaria.

  • Cuando las demás covariables se mantiene constantes, el número esperados de hijos de las mujeres pertenecientes a la religión musulmana es 1.1620863 mayor que el numero esperado de hijos de las mujeres católicas.

  • Cuando las demás covariables se mantiene constantes, el número esperados de hijos de las mujeres pertenecientes a la religión protestante es 1.0114505 mayor que el numero esperado de hijos de las mujeres católicas.

  • Cuando las demás covariables se mantiene constantes, el número esperados de hijos de las mujeres que pertenecen a otra religión (religiones diferentes de la católica, protestante, y la musulmana) es 1.8621009 mayor que el numero esperado de hijos de las mujeres católicas.

  • Cuando las demás covariables se mantiene constantes, el número esperados de hijos de las mujeres que tienen una vivienda en área rural es 1.0819322 mayor que el numero esperado de hijos de las mujeres que no tienen una vivienda en área rural.

  • Cuando las demás covariables se mantiene constantes, el número esperados de hijos incrementa 1.0028879 por cada aumento en el año de nacimiento de la madre.

  • Cuando las demás covariables se mantiene constantes, el número esperados de hijos por cada mujer cambia a un factor de 0.9695245 por cada año adicional que no contrae matrimonio (es decir, una disminución aproximada del (1-0.97) x 100% = 3% en el número esperado de hijos asociados de cada mujer por cada año adicional en el que estas no contraen matrimonio).

Evaluación de sobredispersión

Para evaluar la sobredispersión en los datos, usaremos la razón entre el estadístico \(\chi^2\) de Pearson y \((n-p)\).

X2P = sum(residuals(Modelo,type='pearson')^2)
X2_Gl=X2P/Modelo$df.residual ; X2_Gl
## [1] 0.8306867

Teniendo como resultado para el modelo anterior propuesto, \(\frac{\chi^2}{(n-p)}=0.8306867\). Al tener un resultado cercano a 1, se dice que los datos no parecen presentar sobredispersión.

De igual manera, se calculará el mismo modelo pero esta vez con base a la distribución binomial negativa, a fin de comparar ambos modelos con su respectivo AIC. Por lo que tenemos que:

Modelo2=glm.nb(children~years_school+university+rural+year_birth+age_marriage,data=Datos)
AIC(Modelo); AIC(Modelo2)
## [1] 4244.693
## [1] 4314.096

Como se puede apreciar, el modelo Poisson presenta un AIC de \(4244.693\) y el modelo Binomial negativo un AIC de \(4314.096\), aunque ambos modelos no difieran mucho respecto al estadístico AIC, el modelo Poisson parece presentar un mejor ajuste, teniendo un AIC menor al AIC del modelo Binomial negativo.

Conclusión

Se dará respuesta a la pregunta ¿qué factores influyen al valor esperado del número de hijos por madre?. Tenemos que en el summary los que los p-valor resultan significativos en el modelo para las variables de edad en que la mujer contrajo matrimonio, si la mujer tiene vivienda en área rural y la religión a la que pertecene.

Se realizará una prueba de razón de verosimilitud, para comparar el modelo Poisson propuesto completo con este mismo pero sin la variable religión dicotómica a la que la mujer pertenece, esto para observar si esta variable resulta significativa.

Modelo1=glm(children~years_school+university+rural+year_birth+age_marriage,
            data=Datos,family=poisson)
lrtest0<-lrtest(Modelo1, Modelo); lrtest0
## Likelihood ratio test
## 
## Model 1: children ~ years_school + university + rural + year_birth + age_marriage
## Model 2: children ~ years_school + university + religion + rural + year_birth + 
##     age_marriage
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   6 -2150.0                         
## 2   9 -2113.3  3 73.394  8.005e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Observando los resultados anteriores, existe evidencia estadística para pensar que la variable categórica sí resulta significativa en el modelo.

Por lo tanto, se concluye entonces que los factores religión de la mujer(católicos, protestantes, musulmanes u otros), rural (vivienda en área rural, si/no) influyen positivamente a la esperanza de hijos de cada mujer, si bien, las mujeres que pertenezcan a la religión musulmana o protestante tendrán una esperanza en el número de hijos menor que las mujeres que pertenecen a otra religión; finalmente, la edad en la que la mujer contrajo matrimonio también es un factor que influye negativamente a la esperanza de hijos de cada mujer.