1. Carga de datos y conversión a factores

frijol <- read.csv("Datos_frijol_bioinsumos.csv", sep = ";")
frijol$Bioinsumo <- as.factor(frijol$Bioinsumo)
frijol$Bloque    <- as.factor(frijol$Bloque)
str(frijol)
## 'data.frame':    16 obs. of  4 variables:
##  $ Tratamiento: chr  "T1" "T1" "T1" "T1" ...
##  $ Bioinsumo  : Factor w/ 4 levels "B1","B2","B3",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Bloque     : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ Rendimiento: num  2.2 1.7 2.5 1.9 2.4 2.6 3 2.3 3.1 2.5 ...
head(frijol)
##   Tratamiento Bioinsumo Bloque Rendimiento
## 1          T1        B1      1         2.2
## 2          T1        B1      2         1.7
## 3          T1        B1      3         2.5
## 4          T1        B1      4         1.9
## 5          T2        B2      1         2.4
## 6          T2        B2      2         2.6

Bioinsumo y Bloque quedan como factores (4 niveles cada uno). Si Bloque se dejara como int, R lo trataría como una covariable continua.

2. Análisis exploratorio

medias <- frijol %>%
  group_by(Bioinsumo) %>%
  summarise(n = n(),
            Media = mean(Rendimiento),
            DE = sd(Rendimiento),
            Min = min(Rendimiento),
            Max = max(Rendimiento))
knitr::kable(medias, digits = 3, caption = "Rendimiento (kg/parcela) por bioinsumo")
Rendimiento (kg/parcela) por bioinsumo
Bioinsumo n Media DE Min Max
B1 4 2.075 0.350 1.7 2.5
B2 4 2.575 0.310 2.3 3.0
B3 4 2.875 0.263 2.5 3.1
B4 4 3.175 0.287 3.0 3.6 
tapply(frijol$Rendimiento, frijol$Bloque, mean)
##     1     2     3     4 
## 2.675 2.475 3.025 2.525
ggplot(frijol, aes(x = Bioinsumo, y = Rendimiento, fill = Bioinsumo)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.1, size = 2) +
  labs(title = "Rendimiento de fríjol por bioinsumo",
       x = "Bioinsumo (B1 = testigo)", y = "Rendimiento (kg/parcela)") +
  theme_minimal() + theme(legend.position = "none")

El rendimiento medio aumenta de B1 (testigo, 2.08) a B4 (3.18). La dispersión dentro de cada bioinsumo es parecida.

3. Modelo DBCA y tabla ANOVA

\[Y_{ij} = \mu + \tau_i + \beta_j + \varepsilon_{ij}\]

  • H0: μB1 = μB2 = μB3 = μB4
  • H1: al menos una media de bioinsumo es diferente
modelo <- aov(Rendimiento ~ Bloque + Bioinsumo, data = frijol)
summary(modelo)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Bloque       3   0.74  0.2467    6.00 0.015716 *  
## Bioinsumo    3   2.64  0.8800   21.41 0.000196 ***
## Residuals    9   0.37  0.0411                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cv.model(modelo)   # coeficiente de variación (%)
## [1] 7.579766

Lectura: el efecto del bioinsumo es altamente significativo (F = 21.41; p = 0.0002 < 0.05), por lo que se rechaza H0: al menos un bioinsumo da un rendimiento medio diferente. El bloque también es significativo (F = 6.00; p = 0.016). El CV es de 7.6 %, lo que indica buena precisión experimental.

4. Verificación de supuestos

residuos <- residuals(modelo)

# Normalidad
shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.89968, p-value = 0.07942
# Homogeneidad de varianzas
bartlett.test(Rendimiento ~ Bioinsumo, data = frijol)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Rendimiento by Bioinsumo
## Bartlett's K-squared = 0.2339, df = 3, p-value = 0.9719
leveneTest(Rendimiento ~ Bioinsumo, data = frijol)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   0.234 0.8709
##       12
par(mfrow = c(1, 2))
plot(modelo, which = 1)   # residuales vs ajustados
plot(modelo, which = 2)   # QQ-plot

par(mfrow = c(1, 1))

Interpretación:

  • Shapiro-Wilk: W = 0.900, p = 0.079 > 0.05. No se rechaza la normalidad de los residuales.
  • Bartlett: K² = 0.234, p = 0.972 > 0.05. Levene: p = 0.871 > 0.05. No se rechaza la homogeneidad de varianzas.
  • Los supuestos se cumplen razonablemente, así que el ANOVA es confiable. La independencia la garantiza la aleatorización dentro de cada bloque.

5. Comparación de medias: Tukey (HSD)

tukey <- HSD.test(modelo, "Bioinsumo", alpha = 0.05, console = TRUE)
## 
## Study: modelo ~ "Bioinsumo"
## 
## HSD Test for Rendimiento 
## 
## Mean Square Error:  0.04111111 
## 
## Bioinsumo,  means
## 
##    Rendimiento       std r        se Min Max   Q25  Q50   Q75
## B1       2.075 0.3500000 4 0.1013794 1.7 2.5 1.850 2.05 2.275
## B2       2.575 0.3095696 4 0.1013794 2.3 3.0 2.375 2.50 2.700
## B3       2.875 0.2629956 4 0.1013794 2.5 3.1 2.800 2.95 3.025
## B4       3.175 0.2872281 4 0.1013794 3.0 3.6 3.000 3.05 3.225
## 
## Alpha: 0.05 ; DF Error: 9 
## Critical Value of Studentized Range: 4.41489 
## 
## Minimun Significant Difference: 0.4475788 
## 
## Treatments with the same letter are not significantly different.
## 
##    Rendimiento groups
## B4       3.175      a
## B3       2.875     ab
## B2       2.575      b
## B1       2.075      c
tukey$groups
##    Rendimiento groups
## B4       3.175      a
## B3       2.875     ab
## B2       2.575      b
## B1       2.075      c
grupos <- data.frame(Bioinsumo = rownames(tukey$groups),
                     Media = tukey$groups$Rendimiento,
                     Letra = tukey$groups$groups)

ggplot(grupos, aes(x = reorder(Bioinsumo, -Media), y = Media, fill = Bioinsumo)) +
  geom_col(alpha = 0.8) +
  geom_text(aes(label = Letra), vjust = -0.5, size = 6, fontface = "bold") +
  ylim(0, 3.8) +
  labs(title = "Comparación de medias de Tukey (α = 0.05)",
       x = "Bioinsumo", y = "Rendimiento medio (kg/parcela)") +
  theme_minimal() + theme(legend.position = "none")

# Alternativa base de agricolae: bar.group(tukey$groups, ylim = c(0, 4))
Bioinsumo Media (kg/parcela) Grupo
B4 3.175 a
B3 2.875 ab
B2 2.575 b
B1 (testigo) 2.075 c

DMS de Tukey = 0.448 kg/parcela (CMerror = 0.0411; 9 gl).

6. Conclusión agronómica

  • Todos los bioinsumos superan al testigo. B1 queda solo en el grupo “c”. B2, B3 y B4 aumentan el rendimiento entre 24 % y 53 %.
  • Se recomienda B4 (3.18 kg/parcela, grupo “a”). Es el de mayor rendimiento, 1.10 kg/parcela (+53 %) por encima del testigo, y supera significativamente a B2.
  • B3 es una alternativa válida (2.88 kg/parcela, grupo “ab”). No difiere estadísticamente de B4, así que puede elegirse si es más barato o más fácil de conseguir. La elección final entre B3 y B4 debe apoyarse en un análisis de costo-beneficio.
  • El bloqueo fue útil. El efecto bloque fue significativo (p = 0.016), lo que confirma el gradiente de humedad del lote. La eficiencia relativa del DBCA frente a un DCA es de aproximadamente 2.0: un DCA habría necesitado cerca del doble de repeticiones para lograr la misma precisión. Si no se hubiera bloqueado, esa variación habría inflado el error experimental.
# Eficiencia relativa del DBCA respecto a un DCA
tab <- anova(modelo)
CMb <- tab["Bloque", "Mean Sq"]; CMe <- tab["Residuals", "Mean Sq"]
r <- 4; t <- 4
ER <- ((r - 1) * CMb + r * (t - 1) * CMe) / ((r * t - 1) * CMe)
ER
## [1] 2

7. Publicación

Documento tejido (Knit → HTML) y publicado en RPubs (botón Publish → RPubs en el visor de RStudio).