Experimento en pollos de engorde
Edimer David Jaramillo
21/1/2021
Datos
library(readxl)
datos <- read_excel("experimento_pollos.xlsx", skip = 3)
datos## # A tibble: 96 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
## 7 2 2 6 0 37.3 NA NA NA NA
## 8 2 2 11 0 37.3 NA NA NA NA
## 9 3 1 4 0 37.3 NA NA NA NA
## 10 3 1 16 0 37.3 NA NA NA NA
## # ... with 86 more rows, and 2 more variables: digest_fat <dbl>,
## # digest_prot <dbl>- Conversión de las variables breed y diet a factor:
library(tidyverse)
datos2 <- datos %>%
mutate(breed = as.factor(breed),
diet = as.factor(diet))
datos2## # A tibble: 96 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 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
## 7 2 2 6 0 37.3 NA NA NA NA
## 8 2 2 11 0 37.3 NA NA NA NA
## 9 3 1 4 0 37.3 NA NA NA NA
## 10 3 1 16 0 37.3 NA NA NA NA
## # ... with 86 more rows, and 2 more variables: digest_fat <dbl>,
## # digest_prot <dbl>Curva de crecimiento general
datos2 %>%
ggplot(mapping = aes(x = week, y = weight)) +
geom_point() +
geom_smooth(se = TRUE)
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 promedio por dieta
datos2 %>%
group_by(week, diet) %>%
summarise(promedio_peso = mean(weight)) %>%
ggplot(mapping = aes(x = as.factor(week), y = promedio_peso, color = diet,
group = diet)) +
geom_point() +
geom_line()
“Curva” de crecimiento promedio por dieta y raza
datos2 %>%
group_by(week, diet, breed) %>%
summarise(promedio_peso = mean(weight)) %>%
ggplot(mapping = aes(x = as.factor(week), y = promedio_peso, color = diet,
group = interaction(diet, breed), linetype = breed)) +
geom_point() +
geom_line()
Correlaciones
- Vamos a filtrar sólo datos de la semana 5.
library(janitor)
datos_cor <- datos2 %>%
filter(week == 5) %>%
clean_names()
datos_cor## # 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
Gráfico cuantil cuantil
- Vamos a construir el gráfico inicialmente para la variable peso, sin embargo, lo replicaremos para todas las variables numéricas de interés.
datos_cor %>%
ggplot(mapping = aes(sample = weight)) +
geom_qq() +
geom_qq_line()
Todos los gráficos cuantil cuantil
datos_cor %>%
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
- El juego de hipótesis que subyace a la prueba es el siguiente:
\[H_0: Sí\ hay\ normalidad \\ H_1: No\ hay\ normalidad\]
- El mismo juego de hipótesis podría ser declarado de la siguiente manera:
\[H_0: X \sim N(\mu,\sigma^2) \\ H1: X \nsim N(\mu, \sigma^2)\]
- El nivel de significancia alpha (\(\alpha\)) será del 5% (0.05)
Prueba de Shapiro Wilk peso
shapiro.test(datos_cor$weight)##
## Shapiro-Wilk normality test
##
## data: datos_cor$weight
## W = 0.92461, p-value = 0.2- Conclusión: como el valor p (0.2) es mayor que el nivel de significancia, no existe evidencia para rechazar la hipótesis nula, es decir, que la variable peso se distribuye de forma normal.
Prueba de Shapiro Wilk consumo de alimento
shapiro.test(datos_cor$feed_intake)##
## Shapiro-Wilk normality test
##
## data: datos_cor$feed_intake
## W = 0.95252, p-value = 0.5307Prueba de Shapiro Wilk longitud de la pierna
shapiro.test(datos_cor$length_leg)##
## Shapiro-Wilk normality test
##
## data: datos_cor$length_leg
## W = 0.95171, p-value = 0.5173Prueba de Shapiro Wilk digestibilidad de materia seca
shapiro.test(datos_cor$digest_dm)##
## Shapiro-Wilk normality test
##
## data: datos_cor$digest_dm
## W = 0.88936, p-value = 0.05442Prueba de Shapiro Wilk digestibilidad de fibra cruda
shapiro.test(datos_cor$digest_cf)##
## Shapiro-Wilk normality test
##
## data: datos_cor$digest_cf
## W = 0.85702, p-value = 0.01731Prueba de Shapiro Wilk digestibilidad de grasa
shapiro.test(datos_cor$digest_fat)##
## Shapiro-Wilk normality test
##
## data: datos_cor$digest_fat
## W = 0.92361, p-value = 0.1928Prueba de Shapiro Wilk digestibilidad de proteína
shapiro.test(datos_cor$digest_prot)##
## Shapiro-Wilk normality test
##
## data: datos_cor$digest_prot
## W = 0.89507, p-value = 0.06703Correlaciones no paramétricas de Spearman
datos_cor2 <- datos_cor %>%
select(weight:digest_prot)
datos_cor2## # A tibble: 16 x 7
## weight feed_intake length_leg digest_dm digest_cf digest_fat digest_prot
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1852 16.7 9.9 82.4 -126. 84.2 40.7
## 2 1907 17.9 9.4 70.3 -55.9 93.6 51.0
## 3 1503 19.9 8.1 71.0 40.4 81.4 1.89
## 4 1421 18.0 9 66.0 46.1 86.8 14.1
## 5 1461 14 8.4 73.1 -67.9 93.4 62.8
## 6 1460 15.3 8.6 64.5 -55.2 94.1 60.9
## 7 1253 18 8.1 67.3 55.6 88.4 16.9
## 8 1196 16.8 8 63.9 57.3 87.3 28.4
## 9 1046 11.3 8.6 72.0 -74.6 89.0 46.5
## 10 1057 12.3 8.6 72.4 -78.1 92.8 59.0
## 11 1017 15.5 8 67.7 49.8 81.2 56.5
## 12 1049 15.4 8.3 66.8 62.3 86.9 44.2
## 13 705 7.78 7.6 70.9 -20.5 92.0 40.1
## 14 688 8.52 7.8 72.5 -36.6 92.3 48.1
## 15 704 13.2 7.6 70.7 49.1 84.1 45.5
## 16 699 11.1 7.3 68.1 41.5 85.4 55.0Matriz de correlaciones
cor(datos_cor2, 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 por defecto
library(corrplot)
datos_cor2 %>%
cor(method = "spearman") %>%
corrplot()
Correlograma personalizado
- diag = FALSE: quitar la diagonal del gráfico
- type = “lower” (“upper”): panel superior o inferior del gráfico. En este caso dejamos el panel inferior.
- tl.col = “black”: cambiar el color de las etiquetas a negro.
- tl.srt = 25: rotar las etiquetas para mejorar la lectura del gráfico.
- method = “pie”: cambiar el tipo de gráfico. Por defecto está definido en círculos (circle).
- order = “hclust”: agrupar en función de correlaciones positivas o negativas.
datos_cor2 %>%
cor(method = "spearman") %>%
corrplot(diag = FALSE, type = "lower", tl.col = "black", tl.srt = 25,
method = "pie", order = "hclust")
Prueba t-student
- Para ejecutar la prueba t-student necesitamos validar previamente la normalidad en los grupos que queremos comparar. En este caso vamos a comparar la dieta 1 vs la dieta 2.
- También es necesario validar el supuesto de homocedasticidad (varianza constante) entre los dos grupos.
Verficando la normalidad
peso_dieta1 <- datos_cor %>% filter(diet == "1")
peso_dieta2 <- datos_cor %>% filter(diet == "2")
library(ggpubr)
ggqqplot(peso_dieta1$weight) + labs(title = "Dieta 1")
ggqqplot(peso_dieta2$weight) + labs(title = "Dieta2")
- Prueba de shapiro wilk dieta 1:
shapiro.test(peso_dieta1$weight)##
## Shapiro-Wilk normality test
##
## data: peso_dieta1$weight
## W = 0.91048, p-value = 0.3575- Prueba de shapiro wilk dieta 2:
shapiro.test(peso_dieta2$weight)##
## Shapiro-Wilk normality test
##
## data: peso_dieta2$weight
## W = 0.93125, p-value = 0.5275Homocedasticidad
- Para verificar el supuesto de homocedasticidad podemos utilizar las pruebas de Bartlett o Levene. La primera es de gran utilidad cuando la normalidad se cumple, de lo contrario usaremos la segunda.
\[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(datos_cor$weight ~ datos_cor$diet)##
## Bartlett test of homogeneity of variances
##
## data: datos_cor$weight by datos_cor$diet
## Bartlett's K-squared = 1.3405, df = 1, p-value = 0.247- Ejemplo con el test de Levene:
car::leveneTest(datos_cor$weight ~ datos_cor$diet)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.9932 0.1056
## 14Ejectuando prueba t-student
\[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(datos_cor$weight ~ datos_cor$diet,
alternative = "two.sided",
var.equal = TRUE,
conf.level = 0.95)##
## Two Sample t-test
##
## data: datos_cor$weight by datos_cor$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%, 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. Como el intervalo de confianza para la diferencia de medias contiene al cero (-257.9249 a 591.4249), no existe evidencia para rechazar la hipótesis.
Test de Wilcoxon
wilcox.test(datos_cor$weight ~ datos_cor$diet)##
## Wilcoxon rank sum exact test
##
## data: datos_cor$weight by datos_cor$diet
## W = 39, p-value = 0.5054
## alternative hypothesis: true location shift is not equal to 0Análisis de varianza
Ajuste de modelos
- Vamos a comparar los siguientes modelos:
- Peso = dieta
- Peso = raza (breed)
- Peso = dieta + raza
modelo1 <- aov(weight ~ diet, data = datos_cor)
modelo2 <- aov(weight ~ breed, data = datos_cor)
modelo3 <- aov(weight ~ diet + breed, data = datos_cor)Comparación de modelos
- ¿Cómo elegir el mejor modelo? Podemos utilizar el criterio de información de Akaike (AIC) para seleccionar el mejor modelo. Este criterio es una medida relativa del error, en la medida que sea más bajo mejor será el modelo. Nota: es adimensional, es decir, que se mueve entre -Infinito e Infinito.
AIC(modelo1, modelo2, modelo3)## df AIC
## modelo1 3 240.6753
## modelo2 5 209.0980
## modelo3 6 200.9969Mejor modelo
- Resultados del modelo seleccionado (modelo3):
\[H_0: \mu_{dieta1} = \mu_{dieta2} \\ H_1: \mu_{dieta1} \neq \mu_{dieta2}\]
\[H_0: \mu_{raza1} = \mu_{raza2} = \mu_{raza3} = \mu_{raza4} \\ H_1: \mu_i \neq \mu_j\]
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 ' ' 1Análisis de residuales
- Podemos acceder a los residuales del modelo a través de la función residuals()
- Podemos verificar la normalidad de los residuales a través del gráfico cuantil cuantil o con la prueba de shapiro wilk.
ggqqplot(residuals(modelo3))
shapiro.test(residuals(modelo3))##
## Shapiro-Wilk normality test
##
## data: residuals(modelo3)
## W = 0.93565, p-value = 0.2991- Aplicamos la prueba de Bartlett para verificar la homocedasticidad de los residuales en las dietas:
bartlett.test(residuals(modelo3) ~ datos_cor$diet)##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelo3) by datos_cor$diet
## Bartlett's K-squared = 0.005888, df = 1, p-value = 0.9388- Aplicamos la prueba de Bartlett para verificar la homocedasticidad de los residuales en las razas:
bartlett.test(residuals(modelo3) ~ datos_cor$breed)##
## Bartlett test of homogeneity of variances
##
## data: residuals(modelo3) by datos_cor$breed
## Bartlett's K-squared = 3.243, df = 3, p-value = 0.3556Estimación y comparación de medias (Tukey)
- Podemos estimar las medias a través de la biblioteca emmeans:
library(emmeans)
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.95- Comparamos las medias a través de la prueba de Tukey con la función pairs:
pairs(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 agregar letras para describir mejor las diferencias estadísticas entre razas:
library(multcomp)
library(multcompView)
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.05Análisis de Regresión Lineal
datos_reg <- read_excel("ejemplo_reg.xlsx") %>%
clean_names() %>%
rename(peso_canal = wg, peso_vivo = carcass_wt)
datos_reg## # A tibble: 48 x 28
## cage trt protein_g_kg lipid starch s_f energy peso_canal 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>, peso_vivo <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>- Validamos la normalidad en el peso vivo y en el peso de la canal:
ggqqplot(datos_reg$peso_vivo)
ggqqplot(datos_reg$peso_canal)
shapiro.test(datos_reg$peso_vivo)##
## Shapiro-Wilk normality test
##
## data: datos_reg$peso_vivo
## W = 0.97472, p-value = 0.3822shapiro.test(datos_reg$peso_canal)##
## Shapiro-Wilk normality test
##
## data: datos_reg$peso_canal
## W = 0.97466, p-value = 0.3801- Como el supuesto de normalidad en ambas variables se cumple, vamos a utilizar la correlación paramétrica de Pearson:
cor(datos_reg$peso_vivo, datos_reg$peso_canal, method = "pearson")## [1] 0.7993985- Ejecutamos la prueba de significancia sobre la correlación entre el peso vivo y el peso de la canal:
cor.test(datos_reg$peso_vivo, datos_reg$peso_canal,
method = "pearson",
alternative = "two.sided",
conf.level = 0.95)##
## Pearson's product-moment correlation
##
## data: datos_reg$peso_vivo and datos_reg$peso_canal
## 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- Construimos el modelo de regresión lineal simple:
modelo_reg <- lm(peso_canal ~ peso_vivo, data = datos_reg)
summary(modelo_reg)##
## Call:
## lm(formula = peso_canal ~ peso_vivo, data = datos_reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -175.906 -44.012 -0.053 57.447 147.634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 359.68467 136.86398 2.628 0.0116 *
## peso_vivo 0.76781 0.08508 9.024 9.6e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 72.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 grama que aumente el peso vivo se espera un aumento promedio en el peso de la canal de 0.76781 gramos. Este aumento es estadísticamente significativo (valor p = 9.6e-12). El 63.12% de la variabilidad observada en el peso de la canal es explicada por el peso vivo.
Análisis de residuales:
par(mfrow = c(2, 2))
plot(modelo_reg)
- Podemos realizar la prueba de shapiro wilk para validar la normalidad:
shapiro.test(residuals(modelo_reg))##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_reg)
## W = 0.9836, p-value = 0.7326- Sunpongamos que queremos predecir el peso de la canal para un animal que vivo pesa 1870 gramos:
predict(modelo_reg, newdata = data.frame(peso_vivo = 1870))## 1
## 1795.491Instalar bibliotecas
- car
- ggpubr
- emmeans
- multcomp
- ggsci
- multcompView