Aplicación de Programas Informáticos 2026

Diseño de Experimentos Factoriales

Author

Affiliation

Cátedra de Cálculo Estadístico y Biometría

Facultad de Ciencias Agrarias, UNCUYO

library(tidyverse)
library(rstatix)
library(flextable)
library(car)
library(agricolae)

Experimentos polifactoriales

El objetivo del trabajo fue evaluar el efecto de dos variedades de soja: Pioneer y Novartis y tres momentos de corte del grano: lechoso, pastoso y duro. Se conocía que en la parcela en estudio existía un gradiente de fertilidad de Norte a Sur. La variable respuesta fue la biomasa total (TMS/Ha).

Datos

library(readxl)
soja <- read_excel("Soja.xlsx")

soja$var = as.factor(soja$var)
soja$mom_corte =as.factor(soja$mom_corte)
soja$bloque = as.factor(soja$bloque)
str(soja)

tibble [18 × 4] (S3: tbl_df/tbl/data.frame)
 $ var      : Factor w/ 2 levels "Novartis","Pionner": 2 2 2 1 1 1 2 2 2 1 ...
 $ mom_corte: Factor w/ 3 levels "GD","GL","GP": 2 3 1 2 3 1 2 3 1 2 ...
 $ bloque   : Factor w/ 3 levels "I","II","III": 1 1 1 1 1 1 2 2 2 2 ...
 $ rend     : num [1:18] 16.2 20.8 23.8 16 24.1 ...

Análisis Descriptivo

flextable(soja %>% 
  group_by(mom_corte, var) %>%
  get_summary_stats(rend, show = c("n", "mean", "sd", "min", "max"))
)

var	mom_corte	variable	n	mean	sd	min	max
Novartis	GD	rend	3	19.900	1.990	17.91	21.89
Pionner	GD	rend	3	24.513	2.575	22.37	27.37
Novartis	GL	rend	3	15.893	1.762	14.08	17.60
Pionner	GL	rend	3	16.200	1.620	14.58	17.82
Novartis	GP	rend	3	23.777	1.945	21.69	25.54
Pionner	GP	rend	3	20.383	2.521	17.68	22.67

Boxplot

# 3. Crear el Boxplot agrupado
ggplot(soja, aes(x = mom_corte, y = rend, fill = var)) +
  geom_boxplot() +
  # Etiquetas y títulos para hacerlo profesional
  labs(
    title = "Distribución del Rendimiento de Soja",
    subtitle = "Comparación por Momento de Corte y Variedad",
    x = "Momento de Corte",
    y = "Rendimiento",
    fill = "Variedad"
  ) +
  # Tema limpio para mejor visualización
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    legend.position = "top" # Mueve la leyenda abajo
  )

Modelo

El modelo que utilizaremos es de efectos fijos ya que hemos fijado los niveles de los factores, por lo tanto el modelo a comprobar es: \[Y_{ijl} =\mu + \alpha_j +\beta_l+𝛿_i+(\alpha\beta)_{jl}+\varepsilon_{ijl}\] Para el análisis:

m_aov2 <- aov(rend ~ var*mom_corte+bloque,soja)
m_aov2

Call:
   aov(formula = rend ~ var * mom_corte + bloque, data = soja)

Terms:
                      var mom_corte    bloque var:mom_corte Residuals
Sum of Squares    1.16536 148.72551  23.74898      48.17204  29.17409
Deg. of Freedom         1         2         2             2        10

Residual standard error: 1.708042
Estimated effects may be unbalanced

ADEVA

anova(m_aov2)

Analysis of Variance Table

Response: rend
              Df  Sum Sq Mean Sq F value    Pr(>F)    
var            1   1.165   1.165  0.3994 0.5415530    
mom_corte      2 148.726  74.363 25.4893 0.0001186 ***
bloque         2  23.749  11.874  4.0702 0.0509051 .  
var:mom_corte  2  48.172  24.086  8.2560 0.0076347 ** 
Residuals     10  29.174   2.917                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Supuestos

Normalidad

Para probar los supuestos debemos correr ANOVA pero sin los bloques

m_aov3 <- aov(rend~var:mom_corte,soja)
shapiro.test(resid(m_aov3))

    Shapiro-Wilk normality test

data:  resid(m_aov3)
W = 0.93053, p-value = 0.1984

Homoscedasticidad

plot(m_aov3)

leveneTest(m_aov3)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  0.1058 0.9889
      12               

leveneTest(m_aov3)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  0.1058 0.9889
      12               

Prueba de comparaciones múltiples

tuk_1<-HSD.test(m_aov2, trt = "var", console = F)
tuk_1 <- data.frame(metodos = rownames(tuk_1$groups),
                    tuk_1$groups)
flextable(tuk_1)

metodos	rend	groups
Pionner	20.36556	a
Novartis	19.85667	a

tuk_2<-HSD.test(m_aov2, trt = "mom_corte", console = F)
tuk_2 <- data.frame(metodos = rownames(tuk_2$groups),
                    tuk_2$groups)
flextable(tuk_2 )

metodos	rend	groups
GD	22.20667	a
GP	22.08000	a
GL	16.04667	b

tuk_3<- HSD.test(m_aov2, trt = c("var","mom_corte"), console = F)
tuk_3 <- data.frame(metodos = rownames(tuk_3$groups),
                    tuk_3$groups)
flextable(tuk_3 )

metodos	rend	groups
Pionner:GD	24.51333	a
Novartis:GP	23.77667	a
Pionner:GP	20.38333	ab
Novartis:GD	19.90000	ab
Pionner:GL	16.20000	b
Novartis:GL	15.89333	b

Gráficos interacción

ggplot(soja) +
  aes(x = mom_corte, y = rend,
      color = var,
      group = var) +
  stat_summary(fun.data = "mean_se", 
               geom = "pointrange", 
               size = 0.5) +
  geom_line(stat = "summary", 
            fun = mean, 
            linewidth = 1.2)+
   labs(
    x = "Momento de Corte",
    y = "Rendimiento medio de soja (TMS/ha)"
  ) 

ggplot(soja) +
  aes(x = var, y = rend,
      color = mom_corte,
      group = mom_corte) +
   stat_summary(fun.data = "mean_se", 
               geom = "pointrange",
               size=0.5) +
  geom_line(stat = "summary", 
            fun = mean, 
            linewidth=1.2)+
   labs(
    x = "Variedad",
    y = "Rendimiento medio de soja (TMS/ha)"
  )