Materia: Modelos Lineales
Profesor: Alejandra León

Punto 1 — Base de datos rrates

Descripción

El dataset rrates contiene observaciones sobre la tasa de oxidación del benceno con un catalizador de pentóxido de vanadio. Las variables son: Run (identificador), Conc.O (concentración de oxígeno por 10000 gmole/litro), Temp (temperatura en grados Kelvin) y Rate (velocidad de reacción por 10⁹ gmole/gramo de catalizador por segundo).

Carga de datos

library(GLMsData)
library(car)
library(DescTools)
library(MASS)

data("rrates")
View(rrates)

Selección de familia: Gamma vs Inversa Gaussiana

mod1   <- glm(Rate ~ Conc.O + Temp, data = rrates, family = inverse.gaussian(link = "log"))
mod1.1 <- glm(Rate ~ Conc.O + Temp, data = rrates, family = Gamma(link = "log"))

sd(rrates$Rate) / abs(mean(rrates$Rate))
## [1] 0.6242183
#Como esto da < 0.4 se debe usar gamma (metodo descriptivo)

log(logLik(mod1) / logLik(mod1.1))
## 'log Lik.' 0.00520875 (df=4)
#como dio positivo, entonces se hace inversa gaussiana

Selección de variables

step(mod1)
## Start:  AIC=505.22
## Rate ~ Conc.O + Temp
## 
##          Df Deviance     AIC
## <none>      0.002161  505.22
## - Conc.O  1 0.013878  742.21
## - Temp    1 0.042155 1318.97
## 
## Call:  glm(formula = Rate ~ Conc.O + Temp, family = inverse.gaussian(link = "log"), 
##     data = rrates)
## 
## Coefficients:
## (Intercept)       Conc.O         Temp  
##  -12.056350     0.006133     0.026970  
## 
## Degrees of Freedom: 47 Total (i.e. Null);  45 Residual
## Null Deviance:       0.05324 
## Residual Deviance: 0.002161  AIC: 505.2

Modelo final

modfinal1 <- glm(formula = Rate ~ Conc.O + Temp, family = inverse.gaussian(link = "log"),
                 data = rrates)
summary(modfinal1)
## 
## Call:
## glm(formula = Rate ~ Conc.O + Temp, family = inverse.gaussian(link = "log"), 
##     data = rrates)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.206e+01  6.041e-01  -19.96   <2e-16 ***
## Conc.O       6.133e-03  4.177e-04   14.68   <2e-16 ***
## Temp         2.697e-02  9.455e-04   28.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for inverse.gaussian family taken to be 4.902728e-05)
## 
##     Null deviance: 0.053237  on 47  degrees of freedom
## Residual deviance: 0.002161  on 45  degrees of freedom
## AIC: 505.22
## 
## Number of Fisher Scoring iterations: 4

Interpretación de coeficientes estimados

round(exp(modfinal1$coefficients), 2)
## (Intercept)      Conc.O        Temp 
##        0.00        1.01        1.03
#Por cada unidad que aumenta la concentracion de oxigeno (Conc.O), la velocidad de
#reaccion esperada del benceno se espera que aumente a razon de su coeficiente exp
#Por cada grado Kelvin que aumenta la temperatura, la velocidad de reaccion esperada
#del benceno aumenta a razon de su coeficiente exp con respecto a otras condiciones

Bondad de ajuste

PseudoR2(modfinal1, "Nagelkerke")
## Nagelkerke 
##   0.959409
#existe una excelente explicacion de la velocidad de reaccion del benceno
#por el modelo de regresion inversa gaussiana
1 - pchisq(modfinal1$null.deviance - modfinal1$deviance,
           modfinal1$df.null - modfinal1$df.residual)
## [1] 0.9747851
#Se rechaza Ho, existe evidencia estadistica de que la concentracion de oxigeno
#y la temperatura influyen sobre la velocidad de reaccion del benceno por el
#modelo de regresion inversa gaussiana

Detección de datos atípicos

par(mfrow = c(1, 2))
plot(abs(residuals(modfinal1)))
abline(h = 2, col = "red")
plot(abs(residuals(modfinal1, type = "pearson")))
abline(h = 2, col = "red")

r1 <- abs(residuals(modfinal1))
r2 <- abs(residuals(modfinal1, type = "pearson"))
residuales <- data.frame(r1, r2)
residuales[residuales$r1 > 2 & residuales$r2 > 2, ]
#Teniendo en cuenta las graficas y el dataframe anterior, se identifican
#los datos que pueden considerarse atipicos

Detección de datos influyentes

influence.measures(modfinal1)
## Influence measures of
##   glm(formula = Rate ~ Conc.O + Temp, family = inverse.gaussian(link = "log"),      data = rrates) :
## 
##      dfb.1_  dfb.Cn.O dfb.Temp    dffit cov.r   cook.d    hat inf
## 1  -0.26143 -0.599850  0.27480 -0.71845 0.990 1.33e-01 0.1333    
## 2  -0.26535 -0.410595  0.27179 -0.57186 0.973 8.65e-02 0.0980    
## 3   0.17577  0.080672 -0.17314  0.27266 1.068 2.66e-02 0.0663    
## 4   0.24210 -0.014330 -0.23397  0.35695 1.022 4.57e-02 0.0681    
## 5   0.05059 -0.013613 -0.04851  0.07559 1.149 1.95e-03 0.0737    
## 6   0.02034 -0.012164 -0.01926  0.03217 1.174 3.48e-04 0.0897    
## 7  -0.31247  0.244167  0.29384 -0.51774 1.022 7.59e-02 0.1032    
## 8  -0.18216  0.167915  0.17038 -0.31433 1.150 3.04e-02 0.1160    
## 9  -0.13142  0.121839  0.12290 -0.22714 1.179 1.63e-02 0.1165    
## 10 -0.00844  0.006526  0.00794 -0.01395 1.192 6.48e-05 0.1025    
## 11  0.27993 -0.078323 -0.26831  0.41875 0.999 6.28e-02 0.0741    
## 12  0.49865  0.106739 -0.48680  0.74456 0.696 1.97e-01 0.0653   *
## 13 -0.15409 -0.289037  0.15966 -0.37063 1.121 4.03e-02 0.1128    
## 14 -0.03952 -0.117501  0.04251 -0.13203 1.278 5.62e-03 0.1678   *
## 15  0.08601  0.133717 -0.09613 -0.24608 1.015 1.72e-02 0.0411    
## 16 -0.01156 -0.007342  0.01252  0.02843 1.098 2.75e-04 0.0279    
## 17 -0.00827  0.035160  0.00746  0.03976 1.142 5.40e-04 0.0648    
## 18 -0.02889  0.123544  0.02605  0.13954 1.123 7.02e-03 0.0652    
## 19  0.03939 -0.169178 -0.03549 -0.19090 1.105 1.07e-02 0.0655    
## 20  0.02753 -0.116690 -0.02486 -0.13207 1.124 5.35e-03 0.0646    
## 21 -0.02942  0.098775  0.02753  0.11924 1.107 5.10e-03 0.0506    
## 22  0.01822 -0.046781 -0.01759 -0.06271 1.107 1.25e-03 0.0396    
## 23 -0.11114 -0.070623  0.12044  0.27349 0.926 2.84e-02 0.0279    
## 24  0.09221  0.147274 -0.10321 -0.26605 1.002 1.99e-02 0.0420    
## 25 -0.03458  0.003701  0.03653  0.08161 1.076 2.36e-03 0.0239    
## 26 -0.05447  0.006387  0.05751  0.12855 1.047 6.01e-03 0.0238    
## 27 -0.05409 -0.066367  0.05981  0.14547 1.066 7.60e-03 0.0352    
## 28 -0.04980 -0.059791  0.05502  0.13331 1.072 6.35e-03 0.0348    
## 29  0.04770  0.091958 -0.05398 -0.14728 1.094 6.70e-03 0.0495    
## 30  0.08251  0.158173 -0.09333 -0.25418 1.035 1.86e-02 0.0493    
## 31  0.04771  0.080511 -0.05356 -0.14018 1.087 6.05e-03 0.0439    
## 32  0.11410  0.193769 -0.12814 -0.33599 0.951 3.02e-02 0.0441    
## 33  0.18688  0.095214 -0.19326 -0.23635 1.087 1.61e-02 0.0665    
## 34  0.10127  0.041507 -0.10436 -0.12424 1.119 4.75e-03 0.0599    
## 35 -0.08291  0.084348  0.08111  0.12653 1.101 5.99e-03 0.0490    
## 36 -0.08670  0.060145  0.08585  0.11584 1.096 4.96e-03 0.0435    
## 37 -0.09931  0.018210  0.10019  0.11550 1.093 4.88e-03 0.0414    
## 38 -0.10495 -0.010121  0.10696  0.12135 1.099 5.35e-03 0.0464    
## 39 -0.30917 -0.048975  0.31577  0.35982 0.949 5.27e-02 0.0483    
## 40 -0.28320 -0.091288  0.29094  0.33963 0.992 4.57e-02 0.0551    
## 41  0.36485  0.209206 -0.37817 -0.47159 0.952 5.47e-02 0.0715    
## 42  0.15880  0.078588 -0.16414 -0.19989 1.101 1.18e-02 0.0655    
## 43  0.00296 -0.000438 -0.00299 -0.00343 1.116 3.92e-06 0.0418    
## 44  0.13876 -0.097155 -0.13737 -0.18585 1.063 9.53e-03 0.0436    
## 45  0.02402 -0.021562 -0.02360 -0.03484 1.120 3.92e-04 0.0467    
## 46  0.13471 -0.183940 -0.13009 -0.23903 1.063 1.49e-02 0.0565    
## 47 -0.06882  0.010726  0.06949  0.07979 1.105 2.27e-03 0.0417    
## 48 -0.10065  0.015819  0.10163  0.11671 1.093 4.98e-03 0.0417
influencePlot(modfinal1)
#Se evaluan los criterios de datos influyentes; al tener estadistico de Cook
#menor a 2 no se observan datos que a la luz de todos los criterios puedan
#considerarse como influyentes

Multicolinealidad

cor(rrates[, c("Conc.O", "Temp")])
##            Conc.O       Temp
## Conc.O 1.00000000 0.03742216
## Temp   0.03742216 1.00000000
#Como no hay correlaciones mayores a 0.7 se dice que no hay correlaciones
#entre predictores; no hay multicolinealidad

Predicción

Si se somete el catalizador a una temperatura de 623K con una concentración de oxígeno de 10 * 10000 gmole/Litro, ¿cuál es la velocidad de reacción del benceno?

pred1 <- data.frame(Conc.O = 10, Temp = 623)
predict(modfinal1, pred1, type = "response")
##      1 
## 122.42
#Con una temperatura de 623K y una concentracion de oxigeno de 10 (en unidades
#de 10000 gmole/litro) se espera la velocidad de reaccion del benceno indicada
#en la prediccion anterior (en unidades de 10^9 gmole/gramo de catalizador por segundo)

Punto 2 — Base de datos yieldden

Descripción

El dataset yieldden contiene observaciones sobre el rendimiento medio por planta de tres variedades de cebolla. Las variables son: Yield (rendimiento por planta en gramos), Dens (densidad de siembra en plantas por pie cuadrado) y Var (variedad: 1, 2 o 3).

Carga de datos y preparación

data("yieldden")
View(yieldden)

yieldden$var2 <- ifelse(yieldden$Var == 2, 1, 0)
yieldden$var3 <- ifelse(yieldden$Var == 3, 1, 0)
yieldden <- yieldden[, -3]
yieldden

Selección de familia: Gamma vs Inversa Gaussiana

mod2   <- glm(Yield ~ ., data = yieldden, family = inverse.gaussian(link = "log"))
mod2.1 <- glm(Yield ~ ., data = yieldden, family = Gamma(link = "log"))

sd(yieldden$Yield) / abs(mean(yieldden$Yield))
## [1] 0.6086186
#Como esto da < 0.4 se debe usar gamma (metodo descriptivo)

log(logLik(mod2) / logLik(mod2.1))
## 'log Lik.' -0.01701755 (df=5)
#como dio positivo, entonces se hace inversa gaussiana

Selección de variables

step(mod2)
## Start:  AIC=221.15
## Yield ~ Dens + var2 + var3
## 
##        Df Deviance    AIC
## - var2  1 0.022427 220.48
## <none>    0.021228 221.15
## - var3  1 0.023430 221.60
## - Dens  1 0.238033 460.70
## 
## Step:  AIC=220.8
## Yield ~ Dens + var3
## 
##        Df Deviance    AIC
## - var3  1 0.023672 220.18
## <none>    0.022427 220.80
## - Dens  1 0.241942 462.97
## 
## Step:  AIC=220.42
## Yield ~ Dens
## 
##        Df Deviance    AIC
## <none>    0.023672 220.42
## - Dens  1 0.243071 458.22
## 
## Call:  glm(formula = Yield ~ Dens, family = inverse.gaussian(link = "log"), 
##     data = yieldden)
## 
## Coefficients:
## (Intercept)         Dens  
##     4.54257     -0.05847  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  28 Residual
## Null Deviance:       0.2431 
## Residual Deviance: 0.02367   AIC: 220.4

Modelo final

modfinal2 <- glm(formula = Yield ~ Dens + var2 + var3,
                 family = inverse.gaussian(link = "log"), data = yieldden)
summary(modfinal2)
## 
## Call:
## glm(formula = Yield ~ Dens + var2 + var3, family = inverse.gaussian(link = "log"), 
##     data = yieldden)
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.636361   0.092376  50.190  < 2e-16 ***
## Dens        -0.059688   0.003573 -16.708 2.01e-15 ***
## var2        -0.105770   0.088109  -1.200    0.241    
## var3        -0.128371   0.081381  -1.577    0.127    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for inverse.gaussian family taken to be 0.000897547)
## 
##     Null deviance: 0.243071  on 29  degrees of freedom
## Residual deviance: 0.021228  on 26  degrees of freedom
## AIC: 221.15
## 
## Number of Fisher Scoring iterations: 7

Interpretación de coeficientes estimados

round(exp(modfinal2$coefficients), 2)
## (Intercept)        Dens        var2        var3 
##      103.17        0.94        0.90        0.88
#Por cada planta adicional por pie cuadrado de densidad de siembra, el rendimiento
#esperado por planta cambia a razon de su coeficiente exp con respecto a otras densidades
#Si la variedad de siembra es 2, el rendimiento esperado aumenta a razon de su
#coeficiente exp con respecto a la variedad de referencia (variedad 1)
#Si la variedad de siembra es 3, el rendimiento esperado cambia a razon de su
#coeficiente exp con respecto a la variedad de referencia (variedad 1)

Bondad de ajuste

PseudoR2(modfinal2, "Nagelkerke")
## Nagelkerke 
##   0.912736
#existe una excelente explicacion del rendimiento por planta de las cebollas
#por el modelo de regresion inversa gaussiana
1 - pchisq(modfinal2$null.deviance - modfinal2$deviance,
           modfinal2$df.null - modfinal2$df.residual)
## [1] 0.9739885
#Se rechaza Ho, existe evidencia estadistica de que la densidad de siembra y
#la variedad influyen sobre el rendimiento por planta por el modelo de
#regresion inversa gaussiana

Detección de datos atípicos

par(mfrow = c(1, 2))
plot(abs(residuals(modfinal2)))
abline(h = 2, col = "red")
plot(abs(residuals(modfinal2, type = "pearson")))
abline(h = 2, col = "red")

r1 <- abs(residuals(modfinal2))
r2 <- abs(residuals(modfinal2, type = "pearson"))
residuales <- data.frame(r1, r2)
residuales[residuales$r1 > 2 & residuales$r2 > 2, ]
#Teniendo en cuenta las graficas y el dataframe anterior, se identifican
#los datos que pueden considerarse atipicos

Detección de datos influyentes

influence.measures(modfinal2)
## Influence measures of
##   glm(formula = Yield ~ Dens + var2 + var3, family = inverse.gaussian(link = "log"),      data = yieldden) :
## 
##      dfb.1_ dfb.Dens dfb.var2 dfb.var3   dffit cov.r   cook.d    hat inf
## 1   0.25798  -0.1881 -0.16076 -0.14126  0.2588 1.160 0.018979 0.0917    
## 2   0.06735  -0.0488 -0.04214 -0.03712  0.0676 1.279 0.001140 0.0917    
## 3  -0.02292   0.0154  0.01507  0.01364 -0.0233 1.286 0.000126 0.0912    
## 4  -0.16289   0.1053  0.10937  0.10006 -0.1666 1.232 0.005759 0.0911    
## 5  -0.36928   0.2226  0.25733  0.23980 -0.3839 1.024 0.025482 0.0911    
## 6  -0.31972   0.1409  0.25304  0.24940 -0.3639 1.059 0.024215 0.0942    
## 7  -0.24968   0.0486  0.23343  0.24426 -0.3421 1.123 0.022721 0.1068    
## 8  -0.24938  -0.0384  0.28387  0.31403 -0.4485 1.095 0.038390 0.1315    
## 9   0.12015   0.3031 -0.30276 -0.38061  0.6460 1.305 0.105738 0.2573    
## 10  0.08711   0.5926 -0.43694 -0.57635  1.0601 1.250 0.279714 0.3350    
## 11  0.27115  -0.3474  0.27880 -0.02785  0.6107 0.679 0.121464 0.0806    
## 12  0.17109  -0.2192  0.18581 -0.01757  0.3968 0.966 0.047845 0.0809    
## 13  0.12022  -0.1540  0.14868 -0.01235  0.3002 1.087 0.025939 0.0817    
## 14  0.02895  -0.0371  0.04282 -0.00297  0.0808 1.261 0.001639 0.0831    
## 15 -0.03604   0.0462 -0.07865  0.00370 -0.1330 1.244 0.003791 0.0871    
## 16 -0.03769   0.0483 -0.09696  0.00387 -0.1585 1.233 0.005297 0.0890    
## 17 -0.09192   0.1178 -0.32239  0.00944 -0.5018 0.886 0.040377 0.0929    
## 18  0.01142  -0.0146 -0.11305 -0.00117 -0.1534 1.305 0.005182 0.1265    
## 19  0.06815  -0.0873 -0.18279 -0.00700 -0.2424 1.400 0.012892 0.1957    
## 20  0.23817  -0.3052 -0.41398 -0.02446 -0.5662 1.526 0.067466 0.3139   *
## 21  0.30451  -0.3902 -0.09170  0.24301  0.5530 0.816 0.096331 0.0920    
## 22  0.15917  -0.2040 -0.04793  0.13718  0.2995 1.120 0.025537 0.0917    
## 23  0.03091  -0.0396 -0.00931  0.02861  0.0603 1.280 0.000896 0.0915    
## 24 -0.11706   0.1500  0.03525 -0.12603 -0.2478 1.169 0.011964 0.0912    
## 25 -0.08340   0.1069  0.02511 -0.10922 -0.1994 1.209 0.008133 0.0913    
## 26 -0.16911   0.2167  0.05092 -0.25853 -0.4498 0.947 0.033848 0.0919    
## 27 -0.01454   0.0186  0.00438 -0.13860 -0.2010 1.248 0.008624 0.1100    
## 28  0.00883  -0.0113 -0.00266 -0.37649 -0.5368 0.967 0.051190 0.1216    
## 29 -0.03762   0.0482  0.01133  0.07669  0.1147 1.498 0.003180 0.2249   *
## 30 -0.39834   0.5104  0.11995  0.51649  0.8499 1.573 0.179154 0.3809   *
influencePlot(modfinal2)
#Se evaluan los criterios de datos influyentes; al tener estadistico de Cook
#menor a 2 no se observan datos que a la luz de todos los criterios puedan
#considerarse como influyentes

Multicolinealidad

cor(yieldden[, c("Dens", "var2", "var3")])
##             Dens       var2        var3
## Dens  1.00000000 -0.1550058  0.02098874
## var2 -0.15500577  1.0000000 -0.50000000
## var3  0.02098874 -0.5000000  1.00000000
#Como no hay correlaciones mayores a 0.7 se dice que no hay correlaciones
#entre predictores; no hay multicolinealidad

Predicción

Con el tipo de siembra 2 y una densidad promedio de 10 plantas por pie cuadrado, ¿cuánto es el rendimiento de la planta?

pred2 <- data.frame(Dens = 10, var2 = 1, var3 = 0)
predict(modfinal2, pred2, type = "response")
##        1 
## 51.09626
#Con una densidad de 10 plantas por pie cuadrado y variedad de siembra 2
#se espera el rendimiento por planta (en gramos) indicado en la prediccion anterior