knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
install.packages("car")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
install.packages("ggplot2")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
install.packages("GLMsData")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
install.packages("lmtest")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)
library(car)
## Loading required package: carData
library(ggplot2)
library(GLMsData)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Desarrollo del Taller Punto 1

data(paper)
head(paper)
##   Strength Hardwood
## 1      6.3      1.0
## 2     11.1      1.5
## 3     20.0      2.0
## 4     24.0      3.0
## 5     26.1      4.0
## 6     30.0      4.5
  1. Ajuste del MRLS
mod1 <- lm(Strength ~ Hardwood, data = paper)
summary(mod1)
## 
## Call:
## lm(formula = Strength ~ Hardwood, data = paper)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.986  -3.749   2.938   7.675  15.840 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  21.3213     5.4302   3.926  0.00109 **
## Hardwood      1.7710     0.6478   2.734  0.01414 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.82 on 17 degrees of freedom
## Multiple R-squared:  0.3054, Adjusted R-squared:  0.2645 
## F-statistic: 7.474 on 1 and 17 DF,  p-value: 0.01414
R2 <- summary(mod1)$r.squared
R2_adj <- summary(mod1)$adj.r.squared
cat("R2:", R2, "\n")
## R2: 0.3053739
cat("R2 ajustado:", R2_adj, "\n")
## R2 ajustado: 0.2645135

El modelo lineal muestra una relación positiva y significativa entre la concentración de madera (Hardwood) y la resistencia del papel (Strength). Además, el R² indica que un porcentaje alto de la variabilidad en la resistencia del papel puede ser explicado por la concentración de madera en este modelo simple.

  1. atipicos
par(mfrow = c(2, 2))
plot(mod1)

par(mfrow = c(1, 1))

cooks <- cooks.distance(mod1)
plot(cooks, type = "h", main = "Distancia de Cook - Paper", col = "blue")
abline(h = 4/nrow(paper), col = "red", lty = 2)

which(cooks > 4/nrow(paper))
##  1 18 19 
##  1 18 19
  1. Estimacines
nuevo <- data.frame(Hardwood = 6.5)
predict(mod1, newdata = nuevo, interval = "confidence")
##        fit      lwr      upr
## 1 32.83267 27.01915 38.64619
predict(mod1, newdata = nuevo, interval = "prediction")
##        fit      lwr      upr
## 1 32.83267 7.234445 58.43089
  1. Supuestos
shapiro.test(residuals(mod1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod1)
## W = 0.93232, p-value = 0.1911
bptest(mod1)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 4.1288, df = 1, p-value = 0.04216
dwtest(mod1)
## 
##  Durbin-Watson test
## 
## data:  mod1
## DW = 0.24689, p-value = 9.336e-10
## alternative hypothesis: true autocorrelation is greater than 0

Punto 2: “butterfat”

data("butterfat")
head(butterfat)
##   Butterfat    Breed    Age
## 1      3.74 Ayrshire Mature
## 2      4.01 Ayrshire  2year
## 3      3.77 Ayrshire Mature
## 4      3.78 Ayrshire  2year
## 5      4.10 Ayrshire Mature
## 6      4.06 Ayrshire  2year
  1. Ajute del modelo
mod2 <- lm(Butterfat ~ Breed + Age, data = butterfat)
summary(mod2)
## 
## Call:
## lm(formula = Butterfat ~ Breed + Age, data = butterfat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0202 -0.2373 -0.0640  0.2617  1.2098 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.00770    0.10135  39.541  < 2e-16 ***
## BreedCanadian          0.37850    0.13085   2.893  0.00475 ** 
## BreedGuernsey          0.89000    0.13085   6.802 9.48e-10 ***
## BreedHolstein-Fresian -0.39050    0.13085  -2.984  0.00362 ** 
## BreedJersey            1.23250    0.13085   9.419 3.16e-15 ***
## AgeMature              0.10460    0.08276   1.264  0.20937    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4138 on 94 degrees of freedom
## Multiple R-squared:  0.6825, Adjusted R-squared:  0.6656 
## F-statistic: 40.41 on 5 and 94 DF,  p-value: < 2.2e-16
R2_2 <- summary(mod2)$r.squared
R2_adj2 <- summary(mod2)$adj.r.squared
cat("R2:", R2_2, "\n")
## R2: 0.6824944
cat("R2 ajustado:", R2_adj2, "\n")
## R2 ajustado: 0.6656058

Se ajustó un modelo de regresión múltiple incluyendo las variables categóricas raza (Breed), con los niveles Ayrshire, Canadian, Guernsey, Holstein-Fresian y Jersey, y edad (Age), con los niveles Mature y 2year, para explicar la variable respuesta butterfat. Los coeficientes permiten analizar el efecto de cada raza y categoría de edad sobre la variable respuesta. Además, el R² ajustado permite evaluar qué tan bien se ajusta el modelo en general.

  1. Datos atípicos
par(mfrow = c(2, 2))
plot(mod2)

par(mfrow = c(1, 1))

cooks2 <- cooks.distance(mod2)
plot(cooks2, type = "h", main = "Distancia de Cook - Butterfat", col = "darkgreen")
abline(h = 4/nrow(butterfat), col = "red", lty = 2)

which(cooks2 > 4/nrow(butterfat))
## 32 48 60 82 96 99 
## 32 48 60 82 96 99

Mediante la distancia de Cook y los gráficos de diagnóstico aplicados a la base butterfat, se identificaron las observaciones que superan el umbral de influencia (4/n). Esto permite determinar si existen datos atípicos o influyentes que puedan afectar los coeficientes estimados del modelo múltiple.

  1. Estimación para “animal”
nuevo_vaca <- data.frame(Breed = "Holstein-Fresian", Age = "2year")
predict(mod2, newdata = nuevo_vaca, interval = "confidence")
##      fit      lwr      upr
## 1 3.6172 3.415958 3.818442
  1. Supuestos del modelo
shapiro.test(residuals(mod2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod2)
## W = 0.96347, p-value = 0.007168
bptest(mod2)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod2
## BP = 14.739, df = 5, p-value = 0.01154
vif(mod2)
##       GVIF Df GVIF^(1/(2*Df))
## Breed    1  4               1
## Age      1  1               1