library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(agricolae)
#Situación 1: Alimentación de bovinos
# A. Exploración de Datos
## 1)
NOVILLOS <- read.csv("novillos.csv", stringsAsFactors = FALSE)
str(NOVILLOS)
## 'data.frame': 45 obs. of 6 variables:
## $ TRATAMIENTO : chr "P15" "P15" "P15" "P15" ...
## $ CORRAL : int 1 1 1 1 1 2 2 2 2 2 ...
## $ ANIMAL_NUM : int 1 2 3 4 5 6 7 8 9 10 ...
## $ PESO_INICIAL: int 121 210 200 222 209 315 312 311 309 300 ...
## $ PESO_FINAL : int 447 446 445 446 444 432 435 436 432 431 ...
## $ GANANCIA : int 326 236 245 224 235 117 123 125 123 131 ...
head(NOVILLOS)
## TRATAMIENTO CORRAL ANIMAL_NUM PESO_INICIAL PESO_FINAL GANANCIA
## 1 P15 1 1 121 447 326
## 2 P15 1 2 210 446 236
## 3 P15 1 3 200 445 245
## 4 P15 1 4 222 446 224
## 5 P15 1 5 209 444 235
## 6 P15 2 6 315 432 117
## 2) Transformación y agregación por corral
NOVILLOS_DIETA <- NOVILLOS %>%
group_by(TRATAMIENTO, CORRAL) %>%
summarise(GANANCIA_PROMEDIO = mean(GANANCIA)) %>%
ungroup()
## `summarise()` has grouped output by 'TRATAMIENTO'. You can override using the
## `.groups` argument.
#Transformar TRATAMIENTO a factor
NOVILLOS_DIETA <- NOVILLOS_DIETA %>%
mutate(TRATAMIENTO = factor(TRATAMIENTO))
#Mostrar la tabla resultante
NOVILLOS_DIETA
## # A tibble: 9 × 3
## TRATAMIENTO CORRAL GANANCIA_PROMEDIO
## <fct> <int> <dbl>
## 1 P15 1 253.
## 2 P15 2 124.
## 3 P15 3 256.
## 4 P25 1 291.
## 5 P25 2 124.
## 6 P25 3 234.
## 7 P35 1 212.
## 8 P35 2 125.
## 9 P35 3 176.
## 3,4) Estadísticas descriptivas por tratamiento (a nivel corral)
desc <- NOVILLOS_DIETA %>%
group_by(TRATAMIENTO) %>%
summarise(n = n(),
media = mean(GANANCIA_PROMEDIO),
sd = sd(GANANCIA_PROMEDIO),
se = sd/sqrt(n),
minimo = min(GANANCIA_PROMEDIO),
maximo = max(GANANCIA_PROMEDIO))
desc
## # A tibble: 3 × 7
## TRATAMIENTO n media sd se minimo maximo
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 P15 3 211. 75.6 43.7 124. 256.
## 2 P25 3 216. 84.6 48.8 124. 291.
## 3 P35 3 171. 43.9 25.4 125. 212.
#El orden de ganacia promedio de peso de mayor a menor fue: P25-P15-P35
# Cabe destacar que las desviaciones estándar son grandes.
# Por lo expuesto, no se observa un patrón claro de efecto consistente del nivel de inclusión del subproducto sobre la ganancia de peso de los novillos.
# B. ANÁLISIS DE VARIANZA (ANOVA)
## 5) Hipótesis
#- H0: Las medias de ganancia por corral son iguales entre tratamientos.
#- H1: Al menos una media es diferente.
## 6) Ajuste del modelo ANOVA
modelo_aov <- aov(GANANCIA_PROMEDIO ~ TRATAMIENTO, data = NOVILLOS_DIETA)
summary(modelo_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## TRATAMIENTO 2 3668 1834 0.372 0.704
## Residuals 6 29614 4936
#Conclusión: p > 0.05, por lo tanto no se encontraron diferencias significativas entre tratamientos.
# C. Verificación de Supuestos
## 7) Gráficos diagnósticos
par(mfrow = c(1,2))
plot(modelo_aov, which = 1) # residuos vs ajustados
plot(modelo_aov, which = 2) # Q-Q plot
par(mfrow = c(1,1))
## 8) Normalidad de residuos (Shapiro-Wilk)
shap <- shapiro.test(residuals(modelo_aov))
shap
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_aov)
## W = 0.88492, p-value = 0.1768
## 9) Homogeneidad de varianzas (Levene)
leveneTest(GANANCIA_PROMEDIO ~ TRATAMIENTO, data = NOVILLOS_DIETA, center = median)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.1718 0.8462
## 6
# Residuos vs. Valores Ajustados: nos permite verificar la homogeneidad de varianzas, ya que los puntos se distribuyan sin un patrón claro alrededor de la línea horizontal de cero.
# Q-Q Plot de Residuos: nos permite evaluar la normalidad de los residuos. Ya que los puntos caen aproximadamente sobre la línea diagonal de referencia, los residuos pueden considerarse normales.
# Normalidad de residuos (Shapiro-Wilk): p = 0.1768 por lo tanto no se rechaza normalidad.
# Homogeneidad de varianzas (Levene): p = 0.8462 por lo tanto las varianzas son homogéneas.
# D. Comparación de medias
# Dado que la prueba ANOVA resultó no significativa, no es estadísticamente pertinente realizar una prueba de comparaciones múltiples, ya que no hay diferencias entre los tratamientos.
#Conclusión general situación 1:no existe un efecto estadísticamente significativo del nivel de inclusión del subproducto fibroso en la dieta (15%, 25% o 35%) sobre la ganancia de peso promedio de los novillos Braford. La falta de diferencias podría deberse al bajo número de réplicas (n=3 corrales por tratamiento), a la variabilidad natural entre animales o a que el efecto del nivel de inclusión del subproducto no es marcado bajo las condiciones del ensayo. Sería necesario aumentar el número de corrales y repetir el ensayo para mejorar el poder estadístico, a fin de determinar si realmente existe un efecto biológico.
#Situación 2:Efecto de variedad y tratamiento en rendimiento de poroto
library(dplyr)
library(car)
library(readxl)
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
#A Exploración de datos
poroto <- read_excel("poroto.xlsx")
poroto <- poroto %>% rename(Rep = Repeticion, Dosis = Agroquimico, Rendimiento = `Rendimiento (kg/parcela)`)
poroto$Variedad <- factor(poroto$Variedad)
poroto$Dosis <- factor(poroto$Dosis)
str(poroto)
## tibble [24 × 4] (S3: tbl_df/tbl/data.frame)
## $ Rep : num [1:24] 1 2 3 4 5 6 1 2 3 4 ...
## $ Variedad : Factor w/ 2 levels "V1","V2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Dosis : Factor w/ 2 levels "D0","D1": 1 1 1 1 1 1 2 2 2 2 ...
## $ Rendimiento: num [1:24] 0.89 0.94 0.85 1.12 1.35 1.05 1.78 2 1.86 2.3 ...
head(poroto)
## # A tibble: 6 × 4
## Rep Variedad Dosis Rendimiento
## <dbl> <fct> <fct> <dbl>
## 1 1 V1 D0 0.89
## 2 2 V1 D0 0.94
## 3 3 V1 D0 0.85
## 4 4 V1 D0 1.12
## 5 5 V1 D0 1.35
## 6 6 V1 D0 1.05
#Descriptivas
poroto %>%
group_by(Variedad, Dosis) %>%
summarise(n = n(), media = mean(Rendimiento), sd = sd(Rendimiento), se = sd/sqrt(n))
## `summarise()` has grouped output by 'Variedad'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 6
## # Groups: Variedad [2]
## Variedad Dosis n media sd se
## <fct> <fct> <int> <dbl> <dbl> <dbl>
## 1 V1 D0 6 1.03 0.185 0.0755
## 2 V1 D1 6 2.04 0.249 0.102
## 3 V2 D0 6 1.74 0.316 0.129
## 4 V2 D1 6 2.63 0.121 0.0494
#Las medias aumentan claramente cuando se aplica el agroquímico (D1) en ambas variedades.
#V2 parte de una media más alta en D0 (1.74 vs 1.03) y en D1 (2.63 vs 2.05).
#Las desviaciones estándar son moderadas.
#El rendimiento mas alto se observa en V2:D1 (Mung tratada) y el más bajo en V1:D0 (Alubia sin tratar).
#Gráfico de interacción
interaction.plot(poroto$Dosis, poroto$Variedad, poroto$Rendimiento,
ylab="Rendimiento (kg/parcel)", xlab="Dosis", trace.label="Variedad", lwd=2)
#Ambas líneas tienen una pendiente positiva, lo que confirma que el tratamiento con agroquímico (D1) es beneficioso para el rendimiento en ambas variedades, en comparación con el control (D0). La línea V2 (Mung) está por encima de V1 (Alubia), reflejando el efecto principal de la variedad.
#Las dos líneas no son paralelas. La línea V1 tiene una pendiente más pronunciada que la V2. Esto indica que existe un beneficio en el incremento del agroquímico (D1 vs D0) mayor, en términos absolutos, para la variedad Alubia (V1), ya que su rendimiento base (D0) es muy bajo.
#En conclusión, la exploración de datos sugiere que tanto la variedad como el agroquímico son importantes.
#B Análisis de Varianza Factorial (ANOVA)
#Hipótesis Factor A: Variedad de Poroto (V1: Alubia, V2: Mung)
#H0:μV1=μV2
#H1:μV1 no es igual μV2μ
#Hipótesis Factor B:Dosis de Agroquímico (D0: Sin tratar, D1: Tratado)
#H0:μD0=μD1
#H1:μD0 no es igual μD1
#Hipótesis Interacción: El efecto combinado de Variedad y Agroquímico
#H0(A×B): (V1D1−V1D0)=(V2D1−V2D0 )
#H1(A×B): (V1D1−V1D0)no es igual a (V2D1−V2D0)
# Ajustar el modelo de ANOVA Factorial
modelo <- aov(Rendimiento ~ Variedad * Dosis, data = poroto)
summary(modelo)
## Df Sum Sq Mean Sq F value Pr(>F)
## Variedad 1 2.529 2.529 48.018 9.95e-07 ***
## Dosis 1 5.425 5.425 103.015 2.46e-09 ***
## Variedad:Dosis 1 0.022 0.022 0.422 0.523
## Residuals 20 1.053 0.053
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Efecto Variedad: significativo a α = 0.01. Rechazo H0
#Efecto Dosis: significativo a α = 0.01. Rechazo H0
#Interacción Variedad × Dosis: no significativo. Acepto H0
# En base a lo mencionado, ambos factores tienen efecto principal significativo sobre el rendimiento, sin embargo no hay interacción significativa.
#C. Verificación de Supuestos
par(mfrow=c(1,2))
plot(modelo, which=1) # residuos vs ajustados
plot(modelo, which=2) # Q-Q
par(mfrow=c(1,1))
#Residuos vs ajustados:la distribución de los puntos parece ser aleatoria alrededor de cero. No hay una evidencia clara de heterocedasticidad.
#Q-Q Plot de residuos: la mayoría de los puntos se ajustan bien a la línea diagonal, especialmente en el centro. Visualmente, el supuesto de normalidad parece cumplirse.
shapiro.test(residuals(modelo)) # normalidad
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo)
## W = 0.97151, p-value = 0.7043
leveneTest(Rendimiento ~ interaction(Variedad, Dosis), data = poroto, center = median) # homogeneidad
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.2728 0.3107
## 20
#Test de normalidad: H0: los residuos provienen de una distribución normal. Como el valor p es mayor que α=0.01, no se rechaza la hipótesis nula, por lo tanto los residuos siguen una distribución normal.
#Test de homogeneidad de varianzas: H0: las varianzas de los grupos son iguales (homocedasticidad).Como el valor de p es mayor que α=0.01), no se rechaza la hipótesis nula. Se concluye que las varianzas son homogéneas, cumpliéndose el supuesto de Homogeneidad de Varianzas.
#D. Comparación de Medias: se puede realizar una prueba de comparaciones múltiples para cada efecto principal, ya que ambos resultaron significativos en el ANOVA.
##1. Comparamos el Efecto Principal de Variedad (Factor A)
comparacion_variedad <- poroto %>%
t_test(Rendimiento ~ Variedad, p.adjust.method = "none")
print("Comparación de Medias para Variedad:")
## [1] "Comparación de Medias para Variedad:"
print(comparacion_variedad)
## # A tibble: 1 × 8
## .y. group1 group2 n1 n2 statistic df p
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 Rendimiento V1 V2 12 12 -2.93 21.8 0.00788
#Como p es menor a α, la comparación de medias confirma que existe una diferencia significativa en el rendimiento promedio general entre la Variedad V1 (Alubia) y la Variedad V2 (Mung). La Variedad V2 es significativamente superior a V1.
##2. Comparamos el Efecto Principal de Agroquímico (Factor B)
comparacion_agroquimico <- poroto %>%
t_test(Rendimiento ~ Dosis, p.adjust.method = "none")
print("Comparación de Medias para Agroquímico:")
## [1] "Comparación de Medias para Agroquímico:"
print(comparacion_agroquimico)
## # A tibble: 1 × 8
## .y. group1 group2 n1 n2 statistic df p
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 Rendimiento D0 D1 12 12 -5.75 21.1 0.0000102
# Como p es menor que α, existe una diferencia significativa en el rendimiento promedio general entre las semillas no tratadas (D0) y las tratadas (D1). El tratamiento D1 es significativamente superior al control D0.
# Calcular medias y error estándar (SE) por combinación de factores
resumen <- poroto %>%
group_by(Variedad, Dosis) %>%
summarise(
media = mean(Rendimiento),
sd = sd(Rendimiento),
n = n(),
se = sd / sqrt(n)
)
## `summarise()` has grouped output by 'Variedad'. You can override using the
## `.groups` argument.
# Mostrar tabla resumen
print(resumen)
## # A tibble: 4 × 6
## # Groups: Variedad [2]
## Variedad Dosis media sd n se
## <fct> <fct> <dbl> <dbl> <int> <dbl>
## 1 V1 D0 1.03 0.185 6 0.0755
## 2 V1 D1 2.04 0.249 6 0.102
## 3 V2 D0 1.74 0.316 6 0.129
## 4 V2 D1 2.63 0.121 6 0.0494
# Gráfico de barras con error estándar
ggplot(resumen, aes(x = Dosis, y = media, fill = Variedad)) +
geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
geom_errorbar(
aes(ymin = media - se, ymax = media + se),
width = 0.2,
position = position_dodge(width = 0.9)
) +
labs(
title = "Rendimiento del poroto (kg/parcela)",
x = "Dosis de agroquímico",
y = "Rendimiento medio (± SE)",
fill = "Variedad"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "top"
)
#E. Conclusiones:El análisis estadístico de ANOVA a un nivel de significancia del 1% (α=0.01) determinó que ambos efectos principales (Variedad y Dosis de Agroquímico) influyen de manera altamente significativa en el rendimiento del poroto, mientras que el efecto de la interacción (Variedad × Agroquímico) no resultó significativo. Por ello el efecto de los dos factores puede considerarse aditivo. Por lo tanto se puede recomendar a los productores priorizar la Variedad Mung (V2), ya que, en promedio, supera consistentemente a la Variedad Alubia (V1) en rendimiento, y usar el nuevo agroquímico (D1), ya que incrementa significativamente el rendimiento general, independientemente de la variedad.