Visualización y análisis estadístico de experimentos
Edimer David Jaramillo
19/1/2021
Lectura de datos
library(readxl)
datos <- read_excel("experimento_pollos.xlsx", skip = 3)
head(datos)## # A tibble: 6 x 11
## breed diet pen week weight `feed intake` length_leg digest_dm digest_cf
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 2 0 45.3 NA NA NA NA
## 2 1 1 9 0 45.3 NA NA NA NA
## 3 1 2 7 0 46.7 NA NA NA NA
## 4 1 2 13 0 46.7 NA NA NA NA
## 5 2 1 5 0 38.7 NA NA NA NA
## 6 2 1 15 0 38.7 NA NA NA NA
## # ... with 2 more variables: digest_fat <dbl>, digest_prot <dbl>- Estructura interna de los datos:
str(datos)## tibble [96 x 11] (S3: tbl_df/tbl/data.frame)
## $ breed : num [1:96] 1 1 1 1 2 2 2 2 3 3 ...
## $ diet : num [1:96] 1 1 2 2 1 1 2 2 1 1 ...
## $ pen : num [1:96] 2 9 7 13 5 15 6 11 4 16 ...
## $ week : num [1:96] 0 0 0 0 0 0 0 0 0 0 ...
## $ weight : num [1:96] 45.3 45.3 46.7 46.7 38.7 38.7 37.3 37.3 37.3 37.3 ...
## $ feed intake: num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ length_leg : num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ digest_dm : num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ digest_cf : num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ digest_fat : num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ digest_prot: num [1:96] NA NA NA NA NA NA NA NA NA NA ...- Editando la clase de breed y diet:
library(tidyverse)
datos2 <- datos %>%
mutate(breed = as.factor(breed),
diet = as.factor(diet))- Estructura de la nueva base de datos:
str(datos2)## tibble [96 x 11] (S3: tbl_df/tbl/data.frame)
## $ breed : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 2 2 2 2 3 3 ...
## $ diet : Factor w/ 2 levels "1","2": 1 1 2 2 1 1 2 2 1 1 ...
## $ pen : num [1:96] 2 9 7 13 5 15 6 11 4 16 ...
## $ week : num [1:96] 0 0 0 0 0 0 0 0 0 0 ...
## $ weight : num [1:96] 45.3 45.3 46.7 46.7 38.7 38.7 37.3 37.3 37.3 37.3 ...
## $ feed intake: num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ length_leg : num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ digest_dm : num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ digest_cf : num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ digest_fat : num [1:96] NA NA NA NA NA NA NA NA NA NA ...
## $ digest_prot: num [1:96] NA NA NA NA NA NA NA NA NA NA ...Visualizaciones
Curva de crecimiento general
datos2 %>%
ggplot(mapping = aes(x = week, y = weight)) +
geom_point() +
geom_smooth()
Curva de crecimiento por dieta
datos2 %>%
ggplot(mapping = aes(x = week, y = weight, color = diet)) +
geom_point() +
geom_smooth(se = FALSE)
Curva de crecimiento por raza
datos2 %>%
ggplot(mapping = aes(x = week, y = weight, color = breed)) +
geom_point() +
geom_smooth(se = FALSE)
“Curva” de crecimiento con promedios por dieta
- Para unir los puntos en cada semana por dieta, convertimos la semana a factor.
- Para unir los puntos con la línea es necesario agregar al mapeo estético el argumento group
datos2 %>%
group_by(week, diet) %>%
summarise(promedio_peso = mean(weight, na.rm = TRUE)) %>%
ggplot(mapping = aes(x = as.factor(week), y = promedio_peso,
color = diet, group = diet)) +
geom_point() +
geom_line()
“curva” de crecimiento con promedios por dieta y raza
datos2 %>%
group_by(week, diet, breed) %>%
summarise(promedio_peso = mean(weight, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(mapping = aes(x = as.factor(week), y = promedio_peso,
group = interaction(breed, diet),
color = diet, lty = breed)) +
geom_point() +
geom_line()
Correlaciones
- Vamos a filtrar sólo los datos de la semana 5.
library(janitor)
datos3 <- datos2 %>%
filter(week == 5) %>%
clean_names()
datos3## # A tibble: 16 x 11
## breed diet pen week weight feed_intake length_leg digest_dm digest_cf
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 2 5 1852 16.7 9.9 82.4 -126.
## 2 1 1 9 5 1907 17.9 9.4 70.3 -55.9
## 3 1 2 7 5 1503 19.9 8.1 71.0 40.4
## 4 1 2 13 5 1421 18.0 9 66.0 46.1
## 5 2 1 5 5 1461 14 8.4 73.1 -67.9
## 6 2 1 15 5 1460 15.3 8.6 64.5 -55.2
## 7 2 2 6 5 1253 18 8.1 67.3 55.6
## 8 2 2 11 5 1196 16.8 8 63.9 57.3
## 9 3 1 4 5 1046 11.3 8.6 72.0 -74.6
## 10 3 1 16 5 1057 12.3 8.6 72.4 -78.1
## 11 3 2 3 5 1017 15.5 8 67.7 49.8
## 12 3 2 10 5 1049 15.4 8.3 66.8 62.3
## 13 4 1 8 5 705 7.78 7.6 70.9 -20.5
## 14 4 1 14 5 688 8.52 7.8 72.5 -36.6
## 15 4 2 1 5 704 13.2 7.6 70.7 49.1
## 16 4 2 12 5 699 11.1 7.3 68.1 41.5
## # ... with 2 more variables: digest_fat <dbl>, digest_prot <dbl>Verificación de normalidad
- La normalidad la podemos verificar con el gráfico cuantil cuantil (qqnorm)
datos3 %>%
ggplot(mapping = aes(sample = weight)) +
geom_qq() +
geom_qq_line()
- Construyendo el mismo gráfico anterior para todas las variables numéricas:
datos3 %>%
select(weight:digest_prot) %>%
pivot_longer(cols = everything(), names_to = "variable", values_to = "valores") %>%
ggplot(mapping = aes(sample = valores)) +
facet_wrap(facets = ~variable, scales = "free") +
geom_qq() +
geom_qq_line()
Prueba de Shapiro-Wilk
- Normalidad del peso. El juego de hipótesis (nula y alternativa) es el siguiente:
\[H_0: Sí\ hay\ normalidad \\ H_1: No\ hay\ normalidad\]
- Podríamos escribir el mismo juego de hipótesis de la siguiente manera:
\[H_0: X \sim N(\mu, \sigma^2) \\ H_1: X \nsim N(\mu, \sigma^2)\]
- El nivel de significancia (\(\alpha\)) será igual a 5%.
Shapiro-Wilk sobre el peso
shapiro.test(datos3$weight)##
## Shapiro-Wilk normality test
##
## data: datos3$weight
## W = 0.92461, p-value = 0.2- Conclusión: como el valor p (0.2) es mayor que el nivel de significancia (5%) no existe evidencia para rechazar la hipótesis nula, es decir, que la variable peso se distribuye de forma normal.
Shapiro-Wilk sobre el consumo
shapiro.test(datos3$feed_intake)##
## Shapiro-Wilk normality test
##
## data: datos3$feed_intake
## W = 0.95252, p-value = 0.5307Shapiro-Wilk sobre la longitud de la pierna
shapiro.test(datos3$length_leg)##
## Shapiro-Wilk normality test
##
## data: datos3$length_leg
## W = 0.95171, p-value = 0.5173Shapiro-Wilk sobre la digestibilidad de la materia seca
shapiro.test(datos3$digest_dm)##
## Shapiro-Wilk normality test
##
## data: datos3$digest_dm
## W = 0.88936, p-value = 0.05442Shapiro-Wilk sobre la digestibilidad de la fibra cruda
shapiro.test(datos3$digest_cf)##
## Shapiro-Wilk normality test
##
## data: datos3$digest_cf
## W = 0.85702, p-value = 0.01731- Conclusión: como el valor p (0.01731) es menor que el nivel de significancia (5%) existe evidencia para rechazar la hipótesis, es decir, que la variable digestibilidad de fibra cruda no se distribuye de forma normal.
Shapiro-Wilk sobre la digestibilidad de la grasa
shapiro.test(datos3$digest_fat)##
## Shapiro-Wilk normality test
##
## data: datos3$digest_fat
## W = 0.92361, p-value = 0.1928Shapiro-Wilk sobre la digestibilidad de la proteína
shapiro.test(datos3$digest_prot)##
## Shapiro-Wilk normality test
##
## data: datos3$digest_prot
## W = 0.89507, p-value = 0.06703Correlaciones no paramétricas de Spearman
Matriz de correlaciones
datos3 %>%
select(weight:digest_prot) %>%
cor(method = "spearman")## weight feed_intake length_leg digest_dm digest_cf
## weight 1.00000000 0.7264706 0.77902003 0.01764706 -0.3205882
## feed_intake 0.72647059 1.0000000 0.43755205 -0.37647059 0.3058824
## length_leg 0.77902003 0.4375520 1.00000000 0.12121374 -0.5336361
## digest_dm 0.01764706 -0.3764706 0.12121374 1.00000000 -0.7382353
## digest_cf -0.32058824 0.3058824 -0.53363611 -0.73823529 1.0000000
## digest_fat 0.19705882 -0.2882353 0.31929474 0.07647059 -0.4588235
## digest_prot -0.12352941 -0.5000000 0.05912866 0.22352941 -0.4294118
## digest_fat digest_prot
## weight 0.19705882 -0.12352941
## feed_intake -0.28823529 -0.50000000
## length_leg 0.31929474 0.05912866
## digest_dm 0.07647059 0.22352941
## digest_cf -0.45882353 -0.42941176
## digest_fat 1.00000000 0.47941176
## digest_prot 0.47941176 1.00000000Correlograma 1
library(corrplot)
datos3 %>%
select(weight:digest_prot) %>%
cor(method = "spearman") %>%
corrplot()
Correlograma 2 (editado)
- diag = FALSE: quitar la diagonal del gráfico
- type = “lower” (“upper”): dejar sólo el panel inferior o superior.
- method = “pie”: tipo de visualización. Consultar la ayuda de la función para ver otros métodos.
- tl.col = “black”: cambiando el color de las etiquetas a negro.
- tl.srt = 25: rotación de las etiquetas.
- order = “hclust”: agrupar las correlaciones.
datos3 %>%
select(weight:digest_prot) %>%
cor(method = "spearman") %>%
corrplot(diag = FALSE, type = "lower", meth = "pie", tl.col = "black",
tl.srt = 25, order = "hclust")
Prueba t-student
- Permite comparar medias de dos grupos. En este ejemplo vamos a comparar el peso de la semana 5 entre las dos dietas.
- Necesitamos validar los supuestos de normalidad y homocedasticidad (varianza constante) en los grupos que estamos comparando.
- Las pruebas analíticas para evaluar la homocedasticidad son el test de Bartlett y el test de Levene. La primera es conveniente cuando la normalidad se cumple, de lo contrario usaremos la segunda.
Evaluando la normalidad
- Podemos evaluar la normalidad en cada grupo a través del gráfico cuantil cuantil.
datos3 %>%
ggplot(mapping = aes(sample = weight)) +
facet_wrap(facets = ~diet, scales = "free") +
geom_qq() +
geom_qq_line()
- Vamos a ejecutar la prueba de Shapiro Wilk sobre el peso en cada dieta:
peso_dieta1 <- datos3 %>% filter(diet == "1")
peso_dieta2 <- datos3 %>% filter(diet == "2")
shapiro.test(peso_dieta1$weight)##
## Shapiro-Wilk normality test
##
## data: peso_dieta1$weight
## W = 0.91048, p-value = 0.3575shapiro.test(peso_dieta2$weight)##
## Shapiro-Wilk normality test
##
## data: peso_dieta2$weight
## W = 0.93125, p-value = 0.5275- Vamos a obtener el mismo resultado de las dos pruebas anterior a través del dplyr:
datos3 %>%
group_by(diet) %>%
summarise(valor_p = shapiro.test(weight)$p.value)## # A tibble: 2 x 2
## diet valor_p
## * <fct> <dbl>
## 1 1 0.357
## 2 2 0.528- Podemos visualizar los mismo gráficos cuantil cuantil a través de la función ggqqplot() del paquete ggpubr:
library(ggpubr)
g1 <- ggqqplot(peso_dieta1$weight) + labs(title = "Dieta 1")
g2 <- ggqqplot(peso_dieta2$weight) + labs(title = "Dieta 2")
ggarrange(g1, g2, ncol = 2, nrow = 1)
Evaluando la homocedasticidad
- Como el supuesto de normalidad de satisface vamos a utilizar la prueba de Bartlett:
\[H_0: Sí\ hay\ homocedasticidad \\ H_1: No\ hay\ homocedasticidad\]
\[H_0: \sigma^2_1 = \sigma^2_2 \\ H_1: \sigma^2_1 \neq \sigma^2_2\]
\[H_0: \frac{\sigma^2_1}{\sigma^2_2} = 1 \\ H_1: \frac{\sigma^2_1}{\sigma^2_2} \neq 1\]
bartlett.test(datos3$weight ~ datos3$diet)##
## Bartlett test of homogeneity of variances
##
## data: datos3$weight by datos3$diet
## Bartlett's K-squared = 1.3405, df = 1, p-value = 0.247- El test de Levene lo podemos ejecutar con la función leveneTest() del paquete car:
car::leveneTest(datos3$weight ~ datos3$diet)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.9932 0.1056
## 14Ejecutando prueba t-student
- El juego de hipótesis que subyace a nuestra comparación es el siguiente:
\[H_0: \mu_1 = \mu_2 \\ H_1: \mu_1 \neq \mu_2\]
\[H_0: \mu_1 - \mu_2 = 0 \\ H_1: \mu_1 - \mu_2 \neq 0\]
t.test(datos3$weight ~ datos3$diet,
alternative = "two.sided",
var.equal = TRUE,
conf.level = 0.95)##
## Two Sample t-test
##
## data: datos3$weight by datos3$diet
## t = 0.84216, df = 14, p-value = 0.4139
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -257.9249 591.4249
## sample estimates:
## mean in group 1 mean in group 2
## 1272.00 1105.25- Conclusión: como el valor p (0.4139) es mayor que el nivel de significancia del 5% (\(alpha\)), no existe evidencia para rechazar la hipótesis nula, es decir, que la media del peso de la dieta 1 es igual a la media del peso de la dieta 2. Además, como el intervalo de confianza para la diferencia de medias contiene al cero, también apoya la conclusión anterior, es decir, que no rechazamos la hipótesis nula.
Test de Wilcoxon
wilcox.test(datos3$weight ~ datos3$diet)##
## Wilcoxon rank sum exact test
##
## data: datos3$weight by datos3$diet
## W = 39, p-value = 0.5054
## alternative hypothesis: true location shift is not equal to 0Análisis de varianza
Ajuste de modelos
- En este caso ejecutaremos un análisis varianza bajo tres parametrizaciones distintas. En cada una de ellas nuestra variable respuesta es el peso de los animales, sin embargo, probamos los siguientes modelos:
- Peso = dieta
- Peso = raza
- Peso = dieta + raza
modelo1 <- aov(weight ~ diet, data = datos3)
modelo2 <- aov(weight ~ breed, data = datos3)
modelo3 <- aov(weight ~ diet + breed, data = datos3)- ¿Cómo elegir el mejor? criterio de información de Akaike (AIC). Es una medida relativa del error, adimensional, es decir, que se mueve entre menos infinito e infinito. En la medida que su valor sea más bajo es porque el error del modelo también lo es.
AIC(modelo1, modelo2, modelo3)## df AIC
## modelo1 3 240.6753
## modelo2 5 209.0980
## modelo3 6 200.9969- Como el valor de AIC es más bajo en el modelo 3, optamos por seleccionar este modelo. Vamos a desplegar los resultados del modelo seleccionado:
Mejor modelo
summary(modelo3)## Df Sum Sq Mean Sq F value Pr(>F)
## diet 1 111222 111222 9.681 0.0099 **
## breed 3 2069121 689707 60.033 4.17e-07 ***
## Residuals 11 126376 11489
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Como el valor p (0.0099) de la fuente de variación dieta es menor que el nivel de significancia, diremos entonces que existen diferencias estadísticas entre las dietas.
- Como el valor p (4.17e-07) de la fuente de variación raza es menor que el nivel de significancia, diremos entonces que al menos entre un par de razas existen diferencias estadísticamente significativas para el peso. Esta conclusión está ligada al siguiente juego de hipótesis:
\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4 \\ H_1: \mu_i \neq \mu_j\]
Diagnósticos del modelo
- Podemos acceder a los residuales del modelo a través de la función residuals():
library(ggpubr)
ggqqplot(residuals(modelo3))
- Ejecutamos la prueba de shapiro wilk sobre los residuales del modelo:
shapiro.test(residuals(modelo3))##
## Shapiro-Wilk normality test
##
## data: residuals(modelo3)
## W = 0.93565, p-value = 0.2991- Ejecutamos la prueba de Bartlett sobre los residuales del modelo, esta prueba la hacemos sobre cada grupo, es decir, sobre la dieta y sobre la raza:
bartlett.test(residuals(modelo3) ~ datos3$diet)##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelo3) by datos3$diet
## Bartlett's K-squared = 0.005888, df = 1, p-value = 0.9388bartlett.test(residuals(modelo3) ~ datos3$breed)##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelo3) by datos3$breed
## Bartlett's K-squared = 3.243, df = 3, p-value = 0.3556- Podemos obtener gráficos de residuales con la función plot():
par(mfrow = c(2, 2))
plot(modelo3)
Medias estimadas
- Medias estimadas para la dieta:
library(emmeans)
medias_dieta <- emmeans(modelo3, specs = "diet")
medias_dieta## diet emmean SE df lower.CL upper.CL
## 1 1272 37.9 11 1189 1355
## 2 1105 37.9 11 1022 1189
##
## Results are averaged over the levels of: breed
## Confidence level used: 0.95- Medias estimadas para las razas:
medias_raza <- emmeans(modelo3, specs = "breed")
medias_raza## breed emmean SE df lower.CL upper.CL
## 1 1671 53.6 11 1553 1789
## 2 1342 53.6 11 1225 1460
## 3 1042 53.6 11 924 1160
## 4 699 53.6 11 581 817
##
## Results are averaged over the levels of: diet
## Confidence level used: 0.95Diferencia de medias (Tukey)
- Podemos obtener la diferencia de medias a través de la función pairs():
pairs(medias_dieta)## contrast estimate SE df t.ratio p.value
## 1 - 2 167 53.6 11 3.111 0.0099
##
## Results are averaged over the levels of: breedpairs(medias_raza)## contrast estimate SE df t.ratio p.value
## 1 - 2 328 75.8 11 4.331 0.0056
## 1 - 3 628 75.8 11 8.292 <.0001
## 1 - 4 972 75.8 11 12.821 <.0001
## 2 - 3 300 75.8 11 3.962 0.0102
## 2 - 4 644 75.8 11 8.490 <.0001
## 3 - 4 343 75.8 11 4.529 0.0041
##
## Results are averaged over the levels of: diet
## P value adjustment: tukey method for comparing a family of 4 estimates- Podemos obtener letras para indicar la diferencia estadística a través del paquete multcomp:
library(multcomp)
cld(medias_raza, Letters = letters)## breed emmean SE df lower.CL upper.CL .group
## 4 699 53.6 11 581 817 a
## 3 1042 53.6 11 924 1160 b
## 2 1342 53.6 11 1225 1460 c
## 1 1671 53.6 11 1553 1789 d
##
## Results are averaged over the levels of: diet
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05Regresión Lineal
- Lectura de datos:
datos_reg <- read_excel("datos_reg.xlsx") %>%
clean_names()
datos_reg## # A tibble: 48 x 28
## cage trt protein_g_kg lipid starch s_f energy wg fi fcr
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 8 202 76 296 4.48 12.4 1464. 2127. 1.45
## 2 2 2 204 26 335 14.4 11.2 1468. 2218. 1.51
## 3 3 6 203 31 374 14.4 12.4 1647. 2306. 1.40
## 4 4 9 214 29 414 20.0 13.5 1676. 2222 1.33
## 5 5 4 188 73 276 4.48 11.2 1374. 2198. 1.60
## 6 6 10 221 36 404 14.4 13.5 1530. 2117. 1.38
## 7 7 7 202 52 350 7.54 12.4 1510. 1981. 1.31
## 8 8 11 205 63 376 7.54 13.5 1796. 2314. 1.29
## 9 9 12 219 75 314 4.48 13.5 1725. 2178 1.26
## 10 10 5 194 26 389 20.0 12.4 1737. 2406 1.39
## # ... with 38 more rows, and 18 more variables:
## # protein_digetsibility_jejunum <dbl>, protein_digetsibility_ileum <dbl>,
## # starch_digetsibility_jejunum <dbl>, starch_digetsibility_ileum <dbl>,
## # ame <dbl>, me_ge <dbl>, n_retention <dbl>, excreta_moisture <dbl>,
## # am_en <dbl>, wate_intake_ml_bird_day <dbl>, carcass_wt <dbl>, yield <dbl>,
## # carcass_energy <dbl>, carcass_protein_dm <dbl>,
## # carcass_fat_calculated <dbl>, carcass_dm <dbl>,
## # carcass_protein_as_is <dbl>, carcass_lipid_as_is <dbl>- Verificamos la normalidad en las variables peso vivo y peso de la canal
ggqqplot(datos_reg$wg)
ggqqplot(datos_reg$carcass_wt)
- Calculamos las correlaciones entre el peso vivo y el peso de la canal:
cor(datos_reg$wg, datos_reg$carcass_wt, method = "pearson")## [1] 0.7993985- Realizamos la prueba estadística sobre la correlación de las dos variables:
cor.test(datos_reg$wg, datos_reg$carcass_wt, method = "pearson",
alternative = "two.sided", conf.level = 0.95)##
## Pearson's product-moment correlation
##
## data: datos_reg$wg and datos_reg$carcass_wt
## t = 9.0243, df = 46, p-value = 9.596e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6666947 0.8829766
## sample estimates:
## cor
## 0.7993985- Ajustamos un modelo de regresión lineal simple con la función lm():
modelo_reg <- lm(carcass_wt ~ wg, data = datos_reg)
summary(modelo_reg)##
## Call:
## lm(formula = carcass_wt ~ wg, data = datos_reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -190.699 -51.277 7.246 52.085 221.591
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 279.56981 147.15370 1.900 0.0637 .
## wg 0.83229 0.09223 9.024 9.6e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 75.74 on 46 degrees of freedom
## Multiple R-squared: 0.639, Adjusted R-squared: 0.6312
## F-statistic: 81.44 on 1 and 46 DF, p-value: 9.596e-12Por cada gramo que aumente el peso vivo se espera un aumento promedio de 0.83229 en el peso de la canal. Del 100% de la variabilidad observada en el peso de la canal, el 63.12 (R-Squared) es explicado por la variable predictora peso vivo, el restante quedará en el error aleatorio. Este modelo es estadísticamente significativo, es decir, que la relación del peso de la canal con el peso vivo no es por azar.
Analizamos los residuales del modelo a través de la función plot():
par(mfrow = c(2, 2))
plot(modelo_reg)
- Prueba de normalidad sobre los residuales del modelo:
shapiro.test(residuals(modelo_reg))##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_reg)
## W = 0.9812, p-value = 0.6293- Prueba de Breush-Pagan para validar la homocedasticidad:
library(lmtest)
bptest(modelo_reg)##
## studentized Breusch-Pagan test
##
## data: modelo_reg
## BP = 0.038264, df = 1, p-value = 0.8449- Predicción de nuevos valores: por ejemplo queremos predecir el peso de la canal para un animal que tiene 1487.70 g de peso vivo:
predict(object = modelo_reg, newdata = data.frame(wg = 1487.70))## 1
## 1517.761Instalar bibliotecas
- multcomp
- multcomView
- emmeans
- ggsci
- corrplot