Se cargan las bases de datos que se necesitan y se crea el objeto
NOVILLOS.
Se revisa la estructura de los datos para corrobar que la base de datos se haya cargado correctamente.
## spc_tbl_ [45 × 6] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ TRATAMIENTO : chr [1:45] "P15" "P15" "P15" "P15" ...
## $ CORRAL : num [1:45] 1 1 1 1 1 2 2 2 2 2 ...
## $ ANIMAL_NUM : num [1:45] 1 2 3 4 5 6 7 8 9 10 ...
## $ PESO_INICIAL: num [1:45] 121 210 200 222 209 315 312 311 309 300 ...
## $ PESO_FINAL : num [1:45] 447 446 445 446 444 432 435 436 432 431 ...
## $ GANANCIA : num [1:45] 326 236 245 224 235 117 123 125 123 131 ...
## - attr(*, "spec")=
## .. cols(
## .. TRATAMIENTO = col_character(),
## .. CORRAL = col_double(),
## .. ANIMAL_NUM = col_double(),
## .. PESO_INICIAL = col_double(),
## .. PESO_FINAL = col_double(),
## .. GANANCIA = col_double()
## .. )
## - attr(*, "problems")=<externalptr>
## Rows: 45
## Columns: 6
## $ TRATAMIENTO <chr> "P15", "P15", "P15", "P15", "P15", "P15", "P15", "P15", "…
## $ CORRAL <dbl> 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 1, 1, 1, 1, …
## $ ANIMAL_NUM <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1, 2, …
## $ PESO_INICIAL <dbl> 121, 210, 200, 222, 209, 315, 312, 311, 309, 300, 120, 12…
## $ PESO_FINAL <dbl> 447, 446, 445, 446, 444, 432, 435, 436, 432, 431, 417, 41…
## $ GANANCIA <dbl> 326, 236, 245, 224, 235, 117, 123, 125, 123, 131, 297, 28…
| Name | NOVILLOS |
| Number of rows | 45 |
| Number of columns | 6 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 5 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| TRATAMIENTO | 0 | 1 | 3 | 3 | 0 | 3 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| CORRAL | 0 | 1 | 2.00 | 0.83 | 1 | 1 | 2 | 3 | 3 | ▇▁▇▁▇ |
| ANIMAL_NUM | 0 | 1 | 8.00 | 4.37 | 1 | 4 | 8 | 12 | 15 | ▇▇▇▇▇ |
| PESO_INICIAL | 0 | 1 | 231.89 | 78.12 | 112 | 200 | 221 | 310 | 321 | ▅▁▆▁▇ |
| PESO_FINAL | 0 | 1 | 431.36 | 12.95 | 411 | 417 | 433 | 444 | 449 | ▇▁▅▃▇ |
| GANANCIA | 0 | 1 | 199.47 | 77.51 | 97 | 123 | 198 | 245 | 335 | ▇▁▆▁▅ |
NOVILLOS_DIETA <- NOVILLOS %>%
group_by(TRATAMIENTO, CORRAL) %>%
summarise(GANANCIA_PROMEDIO = mean(GANANCIA, na.rm = TRUE), .groups="drop") %>%
ungroup() %>%
mutate(TRATAMIENTO = factor(TRATAMIENTO))
NOVILLOS_DIETA## # A tibble: 9 × 3
## TRATAMIENTO CORRAL GANANCIA_PROMEDIO
## <fct> <dbl> <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.
#Transformación de variables y agregación por corral
# Agrupar los datos por 'CORRAL' y 'TRATAMIENTO' y calcular la ganancia de peso promedio.
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))# Calcular promedio de ganancia por corral
promedio_corral <- NOVILLOS %>%
group_by(CORRAL) %>%
summarise(ganancia_promedio = mean(GANANCIA, na.rm = TRUE))
# Gráfico de barras
ggplot(promedio_corral, aes(x = CORRAL, y = ganancia_promedio, fill = CORRAL)) +
geom_col(show.legend = FALSE) +
labs(title = "Ganancia de peso promedio por corral",
x = "Corral",
y = "Ganancia promedio (kg)") +
theme_minimal()Se nota claramente que el corral N°1 registra las mayores ganancias de peso, seguidos en orden decreciente por los corrales 3 y 2.
NOVILLOS_DIETA %>%
group_by(TRATAMIENTO) %>%
summarise(
media = mean(GANANCIA_PROMEDIO),
sd = sd(GANANCIA_PROMEDIO),
n = n()
)## # A tibble: 3 × 4
## TRATAMIENTO media sd n
## <fct> <dbl> <dbl> <int>
## 1 P15 211. 75.6 3
## 2 P25 216. 84.6 3
## 3 P35 171. 43.9 3
# Resumen descriptivo de GANANCIA por TRATAMIENTO
NOVILLOS %>%
group_by(TRATAMIENTO) %>%
skim(GANANCIA)| Name | Piped data |
| Number of rows | 45 |
| Number of columns | 6 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | TRATAMIENTO |
Variable type: numeric
| skim_variable | TRATAMIENTO | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|---|
| GANANCIA | P15 | 0 | 1 | 211.13 | 73.22 | 117 | 128.0 | 224 | 266.5 | 326 | ▇▂▆▂▆ |
| GANANCIA | P25 | 0 | 1 | 216.20 | 91.22 | 97 | 125.0 | 223 | 296.5 | 335 | ▇▁▂▁▇ |
| GANANCIA | P35 | 0 | 1 | 171.07 | 62.36 | 98 | 122.5 | 135 | 229.0 | 293 | ▇▁▂▃▁ |
# Boxplot de ganancia por tratamiento
ggplot(NOVILLOS, aes(x = TRATAMIENTO, y = GANANCIA, fill = TRATAMIENTO)) +
geom_boxplot(alpha = 0.7) +
labs(title = "Distribución de la ganancia de peso por tratamiento",
x = "Tratamiento",
y = "Ganancia de peso (kg)") +
theme_minimal()El primer análisis demuestra que la mayor media (la mayor ganancia de peso) se alcanza en los tratamientos P25 y P15. Esto sugiere que hay un efecto del tratamiento sobre la ganancia de peso. Sin embargo, las dispersiones de los datos, medidas como desviaciones estándar, son grandes. Esto puede indicar que los tratamientos tienen mayores efectos para un grupo particular de animales, o para un corral en particular. Para comprobar esto, se puede hacer un ANOVA.
## Df Sum Sq Mean Sq F value Pr(>F)
## TRATAMIENTO 2 3668 1834 0.372 0.704
## Residuals 6 29614 4936
Los resultados indican que no hay evidencia que permita descartar la H0 (no existen diferencias en la ganancia de peso entre tratamientos). Sin embargo, se trata de un ANOVA clásico, paramétrico, que asume:
Por lo que es aconsejable chequear si se cumplen estos supuestos.
Los resultados no parecen demasiado normales. R marca tres corrales que parecen afectar más que el resto a los análisis de normalidad:
Un gráfico de boxplot por corral y tratamiento permite comprobar esta suposición.
# Crear identificador único de corral
NOVILLOS <- NOVILLOS %>%
mutate(CORRAL_ID = paste(TRATAMIENTO, CORRAL, sep = "_"))
# Boxplot por corral único
ggplot(NOVILLOS, aes(x = CORRAL_ID, y = GANANCIA, fill = TRATAMIENTO)) +
geom_boxplot(alpha = 0.7) +
labs(title = "Ganancia de peso por corral (9 corrales únicos)",
x = "Corral (Tratamiento_Corral)",
y = "Ganancia de peso (kg)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))Sin embargo, lo importante para el ANOVA no es que los datos sean normales, sino que los residuos lo sean. Se corre entonces un test de normalidad.
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_dca)
## W = 0.88492, p-value = 0.1768
Los resultados muestran que los datos no se desvían de la normalidad de manera significativa. El valor de p obtenido (p = 0.1768) no permite descartar la hipótesis nula de que los residuos son normales, por lo que el ANOVA parece ser correcto y aplicable en este caso.
Corremos el test de Levene para revisar la homogeneidad de varianzas.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.1718 0.8462
## 6
Nuevamente, el valor de p no permite descartar la hipótesis nula y asumimos que las varianzas son similares en todos los grupos.
A pesar de que el ANOVA demostró que no hay diferencias significativas, podemos cerciorarnos haciendo test post hoc de comparación de medias.
## $statistics
## MSerror Df Mean CV MSD
## 4935.636 6 199.4667 35.22096 176.003
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey TRATAMIENTO 3 4.339195 0.05
##
## $means
## GANANCIA_PROMEDIO std r se Min Max Q25 Q50 Q75
## P15 211.1333 75.64981 3 40.56121 123.8 256.4 188.5 253.2 254.8
## P25 216.2000 84.58463 3 40.56121 124.2 290.6 179.0 233.8 262.2
## P35 171.0667 43.92554 3 40.56121 124.8 212.2 150.5 176.2 194.2
##
## $comparison
## NULL
##
## $groups
## GANANCIA_PROMEDIO groups
## P25 216.2000 a
## P15 211.1333 a
## P35 171.0667 a
##
## attr(,"class")
## [1] "group"
El test de Tukey indica que no hay diferencias significativas entre los tratamientos. Sin embargo, los errores estándar son relativamente grandes comparados con las diferencias de medias, lo que sugiere que la variabilidad impide detectar diferencias significativas.
Perfecto. Te paso un R Markdown completo para Situación 2, listo para copiar y pegar después de tu Situación 1, sin incluir el YAML:
## Rows: 24 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Variedad, Dosis
## dbl (2): Repeticion, Rendimiento
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 24
## Columns: 4
## $ Variedad <chr> "V1", "V1", "V1", "V1", "V1", "V1", "V1", "V1", "V1", "V1"…
## $ Dosis <chr> "D0", "D0", "D0", "D0", "D0", "D0", "D1", "D1", "D1", "D1"…
## $ Repeticion <dbl> 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2…
## $ Rendimiento <dbl> 89, 94, 85, 112, 135, 105, 178, 200, 186, 230, 240, 193, 1…
| Name | porotos |
| Number of rows | 24 |
| Number of columns | 4 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| numeric | 2 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Variedad | 0 | 1 | 2 | 2 | 0 | 2 | 0 |
| Dosis | 0 | 1 | 2 | 2 | 0 | 2 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Repeticion | 0 | 1 | 3.50 | 1.74 | 1 | 2.0 | 3.5 | 5.0 | 6 | ▇▃▃▃▃ |
| Rendimiento | 0 | 1 | 186.38 | 62.65 | 85 | 134.5 | 190.0 | 242.5 | 280 | ▆▃▇▃▇ |
## # A tibble: 6 × 4
## Variedad Dosis Repeticion Rendimiento
## <chr> <chr> <dbl> <dbl>
## 1 V1 D0 1 89
## 2 V1 D0 2 94
## 3 V1 D0 3 85
## 4 V1 D0 4 112
## 5 V1 D0 5 135
## 6 V1 D0 6 105
# Resumen por Variedad
porotos %>%
group_by(Variedad) %>%
summarise(
media = mean(Rendimiento),
sd = sd(Rendimiento),
n = n()
)## # A tibble: 2 × 4
## Variedad media sd n
## <chr> <dbl> <dbl> <int>
## 1 V1 154. 56.8 12
## 2 V2 219. 51.8 12
# Resumen por Dosis
porotos %>%
group_by(Dosis) %>%
summarise(
media = mean(Rendimiento),
sd = sd(Rendimiento),
n = n()
)## # A tibble: 2 × 4
## Dosis media sd n
## <chr> <dbl> <dbl> <int>
## 1 D0 139. 44.5 12
## 2 D1 234. 36.0 12
# Resumen por combinación Variedad x Dosis
porotos %>%
group_by(Variedad, Dosis) %>%
summarise(
media = mean(Rendimiento),
sd = sd(Rendimiento),
n = n()
)## `summarise()` has grouped output by 'Variedad'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 5
## # Groups: Variedad [2]
## Variedad Dosis media sd n
## <chr> <chr> <dbl> <dbl> <int>
## 1 V1 D0 103. 18.5 6
## 2 V1 D1 204. 24.9 6
## 3 V2 D0 174. 31.6 6
## 4 V2 D1 263. 12.1 6
ggplot(porotos, aes(x=interaction(Variedad, Dosis), y=Rendimiento, fill=Variedad)) +
geom_boxplot(alpha=0.7) +
labs(
title="Rendimiento de Poroto por Variedad y Dosis",
x="Variedad y Dosis",
y="Rendimiento (kg/parcela)"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle=45, hjust=1))Los datos muestran que la Variedad 2 tiene un mayor rendimiento que la Variedad 1, independientemente del tratamiento. Además, el tratamiento con agroquímicos aumenta el rendimiento significativamente, independientemente de la variedad.
Cuando se analizan los datos agrupados, se nota que la combinación V2+D1 tiene el mayor rendimiento, mientras que la combinación V1+D0 tiene el menor rendimiento. Además, se ve que el fertilizante parece tener mayor efecto en V1 que en V2 (diferencia 101 kg vs 89 kg), pero V2+D1 sigue siendo el máximo absoluto.
Podemos concluir que tanta la variedad como el tratamiento parecen influir en el rendimiento, por lo que se debería hacer un grafioco de interacciones.
interaction_plot <- porotos %>%
group_by(Variedad, Dosis) %>%
summarise(media = mean(Rendimiento), .groups="drop")
ggplot(interaction_plot, aes(x = Dosis, y = media, color = Variedad, group = Variedad)) +
geom_line(size = 1.2) +
geom_point(size = 3) +
labs(title = "Gráfico de interacción Variedad x Dosis",
x = "Dosis del agroquímico",
y = "Rendimiento promedio (kg/parcela)",
color = "Variedad") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
El gráfico de interacción muestra que las líneas de Variedad × Dosis se acercan bastante a ser paralelas y no se cruzan en el rango estudiado, lo que sugiere que la interacción podría ser débil y debe confirmarse en el ANOVA.
Factor Variedad (A):
Factor Dosis (B):
Interacción (A x B):
## Df Sum Sq Mean Sq F value Pr(>F)
## Variedad 1 25285 25285 48.018 9.95e-07 ***
## Dosis 1 54245 54245 103.015 2.46e-09 ***
## Variedad:Dosis 1 222 222 0.422 0.523
## Residuals 20 10531 527
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Esto confirma que el efecto de la dosis es aproximadamente igual para ambas variedades, y que las diferencias se pueden explicar principalmente por los efectos principales de Variedad y Dosis (ambas con efectos estadisticamente significativos).
En otras palabras: no hay evidencia de que la respuesta a la dosis dependa de la variedad, por lo que puedes interpretar los efectos de Variedad y Dosis de forma independiente.
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_porotos)
## W = 0.97151, p-value = 0.7043
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.2728 0.3107
## 20
Los graficos y los tests parecen indicar que los residuos son normales y que por lo tanto es correcto usar técnicas parametricas.
## $statistics
## MSerror Df Mean CV MSD
## 526.575 20 186.375 12.31239 19.54165
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Variedad 2 2.949998 0.05
##
## $means
## Rendimiento std r se Min Max Q25 Q50 Q75
## V1 153.9167 56.82582 12 6.624292 85 240 102.25 156.5 194.75
## V2 218.8333 51.76667 12 6.624292 133 280 184.25 230.0 262.50
##
## $comparison
## NULL
##
## $groups
## Rendimiento groups
## V2 218.8333 a
## V1 153.9167 b
##
## attr(,"class")
## [1] "group"
## $statistics
## MSerror Df Mean CV MSD
## 526.575 20 186.375 12.31239 19.54165
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Dosis 2 2.949998 0.05
##
## $means
## Rendimiento std r se Min Max Q25 Q50 Q75
## D0 138.8333 44.53157 12 6.624292 85 210 102.25 134 178.75
## D1 233.9167 35.96073 12 6.624292 178 280 198.25 245 262.50
##
## $comparison
## NULL
##
## $groups
## Rendimiento groups
## D1 233.9167 a
## D0 138.8333 b
##
## attr(,"class")
## [1] "group"
medias_tratamientos <- porotos %>%
group_by(Variedad, Dosis) %>%
summarise(media = mean(Rendimiento), .groups="drop")
ggplot(medias_tratamientos, aes(x = interaction(Variedad, Dosis), y = media, fill = Variedad)) +
geom_col(position = "dodge") +
labs(title = "Rendimiento promedio por combinación de Variedad y Dosis",
x = "Tratamiento (Variedad_Dosis)",
y = "Rendimiento promedio (kg/parcela)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))Los resultados demostraron que: 1) La variedad 2 produce mas rendimientos, independientemente del uso o no de fertilizantes. 2) El agegado de fertilizantes causa un aumento significativo del rendimiento, independientemente de la variedad de porotos ensayada. 3) Los datos son normales y sus varioanzas son homogeneas, por lo que es correcto usar herramientas parametrcas para su uso. 4) La combinación Variedad 2 + fertilizante es la que produjo mayores rendimientos, la que produjo menores rendmiento fue Variedad 1 sin fertilizantes. —