Situación 1: Efecto de niveles de inclusión de subproducto fibroso en novillos

A. Exploración de Datos

1. Cargar la base de datos

Se cargan las bases de datos que se necesitan y se crea el objeto NOVILLOS.

library(tidyverse)
library(skimr)
library(summarytools)
library(car)
library(agricolae)

NOVILLOS <- read_csv("novillos.csv")

2. Explorar su estructura

Se revisa la estructura de los datos para corrobar que la base de datos se haya cargado correctamente.

str(NOVILLOS)
## 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>
glimpse(NOVILLOS)
## 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…
skim(NOVILLOS)
Data summary
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 ▇▁▆▁▅

3. Calcular la ganancia de peso promedio por corral

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.

4. Medidas descriptivas de la ganancia de peso por tratamiento

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)
Data summary
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.


B. Análisis de Varianza (ANOVA)

5. Hipótesis

  • H0: No existen diferencias en la ganancia de peso entre tratamientos.
  • H1: Al menos un tratamiento difiere en la ganancia de peso.

6. Ajustar modelo de ANOVA

modelo_dca <- aov(GANANCIA_PROMEDIO ~ TRATAMIENTO, data = NOVILLOS_DIETA)
summary(modelo_dca)
##             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:

  • Normalidad de los residuales
  • Homogeneidad de varianzas

Por lo que es aconsejable chequear si se cumplen estos supuestos.


C. Verificación de Supuestos

7. Gráficos diagnósticos

par(mfrow=c(2,2))
plot(modelo_dca)

par(mfrow=c(1,1))

Los resultados no parecen demasiado normales. R marca tres corrales que parecen afectar más que el resto a los análisis de normalidad:

  • el corral 2 del tratamiento P15 con valores demasiado altos,
  • y los corrales 4 y 5 del tratamiento P25 con valores demasiado bajos y altos, respectivamente.

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.

8. Normalidad de residuos (Shapiro-Wilk)

shapiro.test(residuals(modelo_dca))
## 
##  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.

9. Homogeneidad de varianzas (Levene)

leveneTest(GANANCIA_PROMEDIO ~ TRATAMIENTO, data = NOVILLOS_DIETA)
## 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.

D. Comparación de Medias

10. Prueba de Tukey (si corresponde)

TUKEY <- HSD.test(modelo_dca, "TRATAMIENTO")
TUKEY
## $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"
plot(TUKEY)

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:

Situación 2: Efecto de variedad y tratamiento de semillas en rendimiento de poroto

A. Exploración de Datos

1. Cargar la base de datos

porotos <- read_csv("poroto.csv")
## 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.

2. Explorar la estructura de los datos

glimpse(porotos)
## 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…
skim(porotos)
Data summary
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 ▆▃▇▃▇
head(porotos)
## # 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

3. Medidas descriptivas por factor y combinación de factores

# 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.

4. Gráfico de interacción

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.

B. Análisis de Varianza Factorial (ANOVA)

5. Planteo de hipótesis

  • Factor Variedad (A):

    • H0: No hay diferencias de rendimiento entre variedades
    • H1: Al menos una variedad difiere
  • Factor Dosis (B):

    • H0: No hay diferencias de rendimiento entre dosis
    • H1: Al menos una dosis difiere
  • Interacción (A x B):

    • H0: No hay efecto de interacción entre Variedad y Dosis
    • H1: Existe interacción significativa

6. Ajustar modelo de ANOVA factorial

modelo_porotos <- aov(Rendimiento ~ Variedad * Dosis, data = porotos)
summary(modelo_porotos)
##                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.

C. Verificación de Supuestos

7. Gráficos diagnósticos

par(mfrow = c(2,2))
plot(modelo_porotos)

par(mfrow = c(1,1))

8. Normalidad de residuos (Shapiro-Wilk)

shapiro.test(residuals(modelo_porotos))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo_porotos)
## W = 0.97151, p-value = 0.7043

9. Homogeneidad de varianzas (Levene)

leveneTest(Rendimiento ~ interaction(Variedad, Dosis), data = porotos)
## 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.

D. Comparación de Medias

10. Prueba post-hoc Tukey (si corresponde)

# Comparación de Variedades
TUKEY_variedad <- HSD.test(modelo_porotos, "Variedad")
TUKEY_variedad
## $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"
plot(TUKEY_variedad)

# Comparación de Dosis
TUKEY_dosis <- HSD.test(modelo_porotos, "Dosis")
TUKEY_dosis
## $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"
plot(TUKEY_dosis)

#Comparación de interacciónes (No grafica nada porque la interacción no es significativa???)
#TUKEY_porotos <- HSD.test(modelo_porotos, "Variedad:Dosis")
#TUKEY_porotos
#plot(TUKEY_porotos)

11. Gráfico de barras de medias por tratamiento

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))


E. Conclusiones

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. —