Introducción

Los diseños factoriales permiten estudiar el efecto simultáneo de varios factores sobre una variable de respuesta. Son ampliamente usados en experimentación cuando se desea analizar interacciones entre factores.

Diseño Factorial Completo \(2^k\)

Se utiliza cuando hay k factores, cada uno con 2 niveles. El número total de tratamientos es:

\[ n = 2^k \]

Modelo General del Diseño \(2^k\)

Para dos factores \(A\) y \(B\), con una réplica, el modelo es:

\[ Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \varepsilon_{ijk} \]

Donde:

  • \(Y_{ijk}\): observación del tratamiento \(i,j\), réplica \(k\)
  • \(\mu\): media general
  • \(\alpha_i\): efecto del nivel \(i\) del factor A
  • \(\beta_j\): efecto del nivel \(j\) del factor B
  • \((\alpha\beta)_{ij}\): interacción entre \(A\) y \(B\)
  • \(\varepsilon_{ijk}\): error aleatorio \(\sim N(0, \sigma^2)\)

Análisis de Varianza (ANOVA)

Se descompone la variación total:

\[ \text{SCT} = \text{SCA} + \text{SCB} + \text{SCAB} + \text{SCE} \]

Donde:

Grados de Libertad

Fuente Grados de Libertad
A \(a - 1\)
B \(b - 1\)
A:B (interacción) \((a - 1)(b - 1)\)
Error \(ab(r - 1)\)
Total \(abr - 1\)

Cuadro ANOVA

Fuente SC gl CM = SC/gl F
A SCA \(a - 1\) CMA \(F_A = \frac{CMA}{CME}\)
B SCB \(b - 1\) CMB \(F_B = \frac{CMB}{CME}\)
A:B SCAB \((a-1)(b-1)\) CMAB \(F_{AB} = \frac{CMAB}{CME}\)
Error SCE \(ab(r-1)\) CME
Total SCT \(abr - 1\)

Efectos Estimados

Los efectos de los factores y sus interacciones se estiman como:

\[ \text{Efecto de } A = \frac{\bar{Y}_{A+} - \bar{Y}_{A-}}{2} \]

\[ \text{Interacción } AB = \frac{\bar{Y}_{++} - \bar{Y}_{+-} - \bar{Y}_{-+} + \bar{Y}_{--}}{4} \]

Recomendaciones

Código en R (Ejemplo)

# Cargar datos ejemplo
data <- expand.grid(A = c(-1, 1), B = c(-1, 1), rep = 1:2)
set.seed(123)
data$Y <- rnorm(8, mean = 10 + 2*data$A + 3*data$B + 1*data$A*data$B)

# Ajustar modelo factorial
modelo <- aov(Y ~ A * B, data = data)
summary(modelo)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            1  27.02   27.02  29.819 0.005468 ** 
## B            1  70.63   70.63  77.956 0.000908 ***
## A:B          1   1.03    1.03   1.136 0.346533    
## Residuals    4   3.62    0.91                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En los modelos de efectos aleatorios, los niveles de uno o más factores no son fijos, sino que se consideran extraídos aleatoriamente de una población más amplia. Este enfoque permite hacer inferencias sobre todos los niveles potenciales del factor, no solo los observados.

Ejemplo común

Diseño factorial \(A \times B\) con efectos aleatorios en ambos factores.

Modelo Lineal con Efectos Aleatorios

Supongamos dos factores aleatorios \(A\) y \(B\), cada uno con \(a\) y \(b\) niveles, respectivamente, y \(r\) réplicas por combinación.

El modelo es:

\[ Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \varepsilon_{ijk} \]

Donde:

Todos los efectos son no correlacionados entre sí.

Varianzas Esperadas

La tabla de varianzas esperadas ayuda a construir pruebas F apropiadas.

Fuente Varianza Esperada del CM
A \(\sigma^2 + r\sigma^2_{\alpha\beta} + rb\sigma^2_\alpha\)
B \(\sigma^2 + r\sigma^2_{\alpha\beta} + ra\sigma^2_\beta\)
A:B \(\sigma^2 + r\sigma^2_{\alpha\beta}\)
Error \(\sigma^2\)

Tabla ANOVA

Fuente SC gl CM = SC/gl F
A SCA \(a - 1\) CMA \(F_A = \frac{CMA}{CM_{AB}}\)
B SCB \(b - 1\) CMB \(F_B = \frac{CMB}{CM_{AB}}\)
A:B SCAB \((a-1)(b-1)\) CMAB \(F_{AB} = \frac{CMAB}{CME}\)
Error SCE \(abr - ab\) CME
Total SCT \(abr - 1\)

Nota: La prueba F para factores aleatorios no usa siempre el CM del error como denominador. Se debe usar el denominador correcto según la varianza esperada.

Componentes de Varianza

Los componentes de varianza se estiman como:

Código en R (Ejemplo con lme4)

# Instalar si es necesario
# install.packages("lme4")

library(lme4)
## Warning: package 'lme4' was built under R version 4.4.3
## Cargando paquete requerido: Matrix
# Simular datos
set.seed(123)
a <- 3
b <- 4
r <- 2
data <- expand.grid(A = factor(1:a), B = factor(1:b), Rep = 1:r)
n <- nrow(data)

# Efectos aleatorios
alpha <- rnorm(a, 0, sqrt(2))
beta <- rnorm(b, 0, sqrt(3))
alpha_beta <- matrix(rnorm(a*b, 0, sqrt(1)), nrow = a, ncol = b)
epsilon <- rnorm(n, 0, sqrt(4))

data$Y <- with(data,
               10 +
               alpha[as.numeric(A)] +
               beta[as.numeric(B)] +
               mapply(function(i, j) alpha_beta[i,j], as.numeric(A), as.numeric(B)) +
               epsilon)

# Modelo con lmer
modelo <- lmer(Y ~ 1 + (1|A) + (1|B), data = data)

summary(modelo)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Y ~ 1 + (1 | A) + (1 | B)
##    Data: data
## 
## REML criterion at convergence: 103
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.32716 -0.69099 -0.05525  0.74246  1.72242 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  B        (Intercept) 2.078    1.442   
##  A        (Intercept) 5.115    2.262   
##  Residual             2.840    1.685   
## Number of obs: 24, groups:  B, 4; A, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   11.099      1.531   7.251

Los modelos mixtos combinan factores fijos y aleatorios en el mismo análisis. Son apropiados cuando se desea inferir sobre niveles específicos de un factor (efecto fijo) y también sobre una población de niveles (efecto aleatorio).

Este documento presenta la teoría del modelo mixto \(A \times B\), incluyendo estructura de varianza-covarianza, esperanzas de cuadrados medios y análisis de componentes de varianza.

2. Formulación del Modelo Mixto

Consideremos:

El modelo lineal es:

\[ Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \varepsilon_{ijk} \]

Donde:

Todos los efectos aleatorios son independientes entre sí.

3. Varianzas y Covarianzas

3.1 Varianza de una observación

Para una observación individual:

\[ \operatorname{Var}(Y_{ijk}) = \sigma^2 + \sigma^2_{\alpha\beta} + \sigma^2_\beta \]

3.2 Covarianzas

  • Entre réplicas de una misma combinación \((i, j)\): \[ \operatorname{Cov}(Y_{ijk}, Y_{ijl}) = \sigma^2_\beta + \sigma^2_{\alpha\beta} \]

  • Entre combinaciones con mismo nivel de \(B\): \[ \operatorname{Cov}(Y_{ijk}, Y_{i'jk}) = \sigma^2_\beta \]

  • Entre combinaciones diferentes (diferente \(i, j\)): \[ \operatorname{Cov}(Y_{ijk}, Y_{i'j'k'}) = 0 \]

4. Esperanza de Cuadrados Medios (ECM)

Fuente ECM (Esperanza del Cuadrado Medio)
A (fijo) \(\sigma^2 + r\sigma^2_{\alpha\beta} + br \cdot \text{MSA}^*\)
B (aleatorio) \(\sigma^2 + r\sigma^2_{\alpha\beta} + ar\sigma^2_\beta\)
A:B \(\sigma^2 + r\sigma^2_{\alpha\beta}\)
Error \(\sigma^2\)

Donde \(\text{MSA}^*\) representa el efecto fijo de A.

5. Análisis de Varianza (ANOVA)

Fuente SC gl CM = SC/gl F
A (fijo) SCA \(a - 1\) CMA \(F = \frac{CMA}{CM_{AB}}\)
B (aleat.) SCB \(b - 1\) CMB \(F = \frac{CMB}{CM_{AB}}\)
A:B SCAB \((a-1)(b-1)\) CMAB
Error SCE \(abr - ab\) CME
Total SCT \(abr - 1\)

6. Estimación de Componentes de Varianza

Los componentes se estiman a partir de los cuadrados medios:

7. Ejemplo en R

# Librerías
library(lme4)
library(lmerTest)
## Warning: package 'lmerTest' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
# Parámetros del diseño
a <- 3
b <- 4
r <- 2

# Datos simulados
set.seed(123)
data <- expand.grid(A = factor(1:a), B = factor(1:b), Rep = 1:r)
n <- nrow(data)

# Efectos simulados
alpha <- c(0, 2, -2)  # efecto fijo
beta <- rnorm(b, 0, sqrt(1.5))
alpha_beta <- matrix(rnorm(a * b, 0, sqrt(1)), nrow = a)
error <- rnorm(n, 0, sqrt(2))

# Generar respuesta
data$Y <- with(data,
  10 +
    alpha[as.numeric(A)] +
    beta[as.numeric(B)] +
    mapply(function(i, j) alpha_beta[i, j], as.numeric(A), as.numeric(B)) +
    error
)

# Ajustar modelo mixto: A fijo, B aleatorio
modelo <- lmer(Y ~ A + (1|B) + (1|A:B), data = data)
## boundary (singular) fit: see help('isSingular')
summary(modelo)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ A + (1 | B) + (1 | A:B)
##    Data: data
## 
## REML criterion at convergence: 86.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.28139 -0.69343  0.06363  0.52135  1.94534 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  A:B      (Intercept) 0.000    0.000   
##  B        (Intercept) 1.980    1.407   
##  Residual             2.056    1.434   
## Number of obs: 24, groups:  A:B, 12; B, 4
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   9.7123     0.8672  4.9602  11.200 0.000104 ***
## A2            2.6828     0.7170 18.0000   3.742 0.001493 ** 
## A3           -0.6415     0.7170 18.0000  -0.895 0.382782    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##    (Intr) A2    
## A2 -0.413       
## A3 -0.413  0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
grupo <- rep(paste0("Grupo", 1:4), each = 20)

thorax <- c(
  # Grupo 1
  0.64, 0.70, 0.72, 0.72, 0.72, 0.76, 0.84, 0.84, 0.84, 0.84,
  0.84, 0.88, 0.88, 0.92, 0.92, 0.92, 0.92, 0.92, 0.92, 0.94,
  # Grupo 2
  0.64, 0.68, 0.72, 0.76, 0.76, 0.80, 0.82, 0.82, 0.84, 0.84,
  0.84, 0.84, 0.84, 0.88, 0.88, 0.88, 0.88, 0.88, 0.92, 0.92,
  # Grupo 3
  0.68, 0.68, 0.72, 0.76, 0.78, 0.80, 0.84, 0.84, 0.84, 0.84,
  0.84, 0.88, 0.90, 0.90, 0.90, 0.90, 0.90, 0.92, 0.92, 0.92,
  # Grupo 4
  0.64, 0.64, 0.68, 0.72, 0.72, 0.74, 0.76, 0.78, 0.80, 0.80,
  0.82, 0.82, 0.84, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 0.92
)

longevidad <- c(
  # Grupo 1
  40, 37, 44, 47, 47, 47, 54, 61, 71, 75,
  89, 58, 59, 58, 62, 70, 72, 75, 96, 75,
  # Grupo 2
  46, 42, 65, 46, 58, 42, 50, 80, 63, 65,
  70, 70, 72, 70, 70, 72, 76, 90, 76, 92,
  # Grupo 3
  21, 40, 44, 54, 36, 40, 48, 53, 60, 60,
  65, 68, 48, 56, 68, 75, 81, 48, 68, 68,
  # Grupo 4
  16, 19, 19, 32, 33, 33, 42, 33, 26, 30,
  40, 54, 34, 42, 47, 54, 54, 56, 60, 44
)

# Crear data frame
df <- data.frame(
  grupo = factor(grupo),
  thorax = thorax,
  longevidad = longevidad
)

# Ajustar modelo ANCOVA con interacción (pendiente diferente por grupo)
modelo <- lm(longevidad ~ grupo * thorax, data = df)
# Coeficientes estimados
summary(modelo)$coefficients
##                      Estimate Std. Error     t value     Pr(>|t|)
## (Intercept)        -46.213643   20.39130 -2.26634079 2.643664e-02
## grupoGrupo2         -7.448519   31.68939 -0.23504773 8.148389e-01
## grupoGrupo3        -15.862616   31.52251 -0.50321546 6.163488e-01
## grupoGrupo4        -16.891590   29.21626 -0.57815720 5.649620e-01
## thorax             129.572714   24.31053  5.32990089 1.077612e-06
## grupoGrupo2:thorax  15.697557   38.14237  0.41155171 6.818906e-01
## grupoGrupo3:thorax  10.196092   37.51196  0.27180910 7.865467e-01
## grupoGrupo4:thorax  -2.373174   35.65172 -0.06656548 9.471121e-01
# Tabla ANOVA del modelo con interacción
anova(modelo)
## Analysis of Variance Table
## 
## Response: longevidad
##              Df Sum Sq Mean Sq  F value    Pr(>F)    
## grupo         3 8755.4  2918.5  30.8483 6.200e-13 ***
## thorax        1 9487.1  9487.1 100.2782 2.785e-15 ***
## grupo:thorax  3   27.2     9.1   0.0959    0.9621    
## Residuals    72 6811.7    94.6                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ajustar el modelo completo (con interacción)
modelo_completo <- lm(longevidad ~ grupo * thorax, data = df)

# Ajustar el modelo reducido (sin interacción)
modelo_reducido <- lm(longevidad ~ grupo + thorax, data = df)

# Prueba de hipótesis: ¿las pendientes son iguales? (es decir, sin interacción)
anova(modelo_reducido, modelo_completo)
## Analysis of Variance Table
## 
## Model 1: longevidad ~ grupo + thorax
## Model 2: longevidad ~ grupo * thorax
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     75 6839.0                           
## 2     72 6811.7  3    27.223 0.0959 0.9621
moscas=df

modelo1 <- lm(longevidad ~ grupo + thorax, data = moscas)
modelo2 <- lm(longevidad ~ grupo * thorax, data = moscas)
library(ggplot2)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Datos originales
p <- ggplot(moscas, aes(x = thorax, y = longevidad, color = grupo)) +
  geom_point(alpha = 0.6) +
  labs(title = "Rectas ajustadas por grupo",
       x = "Tamaño del tórax",
       y = "Longevidad") +
  theme_minimal()

# Generar grid para predicciones
nuevo <- expand.grid(
  thorax = seq(min(moscas$thorax), max(moscas$thorax), length.out = 100),
  grupo = levels(moscas$grupo)
)

# Predicciones para cada modelo
nuevo$pred_m1 <- predict(modelo1, newdata = nuevo)
nuevo$pred_m2 <- predict(modelo2, newdata = nuevo)

# Agregar líneas al gráfico
p +
  geom_line(data = nuevo, aes(y = pred_m1), linewidth = 1, linetype = "dashed") +
  geom_line(data = nuevo, aes(y = pred_m2), linewidth = 1, linetype = "solid") +
  scale_color_brewer(palette = "Dark2") +
  labs(caption = "Líneas sólidas: modelo con pendientes distintas\nLíneas punteadas: modelo con pendientes iguales")

# Modelo con pendientes iguales
mod_paralelas <- lm(longevidad ~ grupo + thorax, data = moscas)

# ANOVA secuencial para probar ambas hipótesis
anova(mod_paralelas)
## Analysis of Variance Table
## 
## Response: longevidad
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## grupo      3 8755.4  2918.5  32.006 1.989e-13 ***
## thorax     1 9487.1  9487.1 104.041 8.053e-16 ***
## Residuals 75 6839.0    91.2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod_paralelas)
## 
## Call:
## lm(formula = longevidad ~ grupo + thorax, data = moscas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.7916  -6.9527  -0.3121   4.7475  26.3432 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -50.301     11.201  -4.491 2.53e-05 ***
## grupoGrupo2    5.514      3.024   1.823   0.0722 .  
## grupoGrupo3   -7.338      3.020  -2.430   0.0175 *  
## grupoGrupo4  -18.609      3.057  -6.088 4.49e-08 ***
## thorax       134.473     13.184  10.200 8.05e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.549 on 75 degrees of freedom
## Multiple R-squared:  0.7273, Adjusted R-squared:  0.7128 
## F-statistic: 50.01 on 4 and 75 DF,  p-value: < 2.2e-16
# SEGUNDO PUNTO

# Crear los datos
datos <- data.frame(
  Lote = rep(1:5, each = 5),
  Dia  = rep(1:5, times = 5),
  Tratamiento = c("A","B","D","C","E",   # Lote 1
                  "C","E","A","D","B",   # Lote 2
                  "B","A","C","E","D",   # Lote 3
                  "D","C","E","B","A",   # Lote 4
                  "E","D","B","A","C"),  # Lote 5
  Tiempo = c(8,7,1,7,3,    # Lote 1
             11,2,7,3,8,   # Lote 2
             4,9,10,1,5,   # Lote 3
             6,8,6,6,10,   # Lote 4
             4,2,3,8,8)    # Lote 5
)

# Convertir a factores
datos$Tratamiento <- factor(datos$Tratamiento)
datos$Lote <- factor(datos$Lote)
datos$Dia <- factor(datos$Dia)

# Modelo lineal para Cuadro Latino
modelo <- aov(Tiempo ~ Tratamiento + Lote + Dia, data = datos)

# Tabla ANOVA
summary(modelo)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Tratamiento  4 141.44   35.36  11.309 0.000488 ***
## Lote         4  15.44    3.86   1.235 0.347618    
## Dia          4  12.24    3.06   0.979 0.455014    
## Residuals   12  37.52    3.13                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#c
modelo_cl <- aov(Tiempo ~ Tratamiento + Lote + Dia, data = datos)
summary(modelo_cl)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Tratamiento  4 141.44   35.36  11.309 0.000488 ***
## Lote         4  15.44    3.86   1.235 0.347618    
## Dia          4  12.24    3.06   0.979 0.455014    
## Residuals   12  37.52    3.13                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelo_dca <- aov(Tiempo ~ Tratamiento, data = datos)
summary(modelo_dca)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Tratamiento  4  141.4   35.36   10.85 7.67e-05 ***
## Residuals   20   65.2    3.26                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelo_dbca <- aov(Tiempo ~ Tratamiento + Lote, data = datos)
summary(modelo_dbca)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Tratamiento  4 141.44   35.36  11.370 0.000146 ***
## Lote         4  15.44    3.86   1.241 0.333144    
## Residuals   16  49.76    3.11                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Obtener los cuadrados medios del error
cme_cl   <- summary(modelo_cl)[[1]]["Residuals", "Mean Sq"]
cme_dca  <- summary(modelo_dca)[[1]]["Residuals", "Mean Sq"]
cme_dbca <- summary(modelo_dbca)[[1]]["Residuals", "Mean Sq"]

# Calcular eficiencias relativas
eficiencia_dca_vs_cl  <- cme_dca / cme_cl
eficiencia_dbca_vs_cl <- cme_dbca / cme_cl

# Mostrar resultados
# Base original con nombres
datos <- data.frame(
  Lote = rep(1:5, each = 5),
  Dia = rep(1:5, times = 5),
  Tratamiento = c("A", "B", "D", "C", "E",
                  "C", "E", "A", "D", "B",
                  "B", "A", "C", "E", "D",
                  "D", "C", "E", "B", "A",
                  "E", "D", "B", "A", "C"),
  Tiempo = c(8, 7, 1, 7, 3,
             11, 2, 7, 3, 8,
             4, 9, 10, 1, 5,
             6, 8, 6, 6, 10,
             4, 2, 3, 8, 8)
)

# Estimar valor faltante en Lote 3, Día 4, Tratamiento E
# Guardar valor observado para comparar
valor_observado <- datos$Tiempo[datos$Lote == 3 & datos$Dia == 4 & datos$Tratamiento == "E"]

# Eliminar el valor para estimar
datos$Tiempo[datos$Lote == 3 & datos$Dia == 4 & datos$Tratamiento == "E"] <- NA

# Calcular medias necesarias
media_general <- mean(datos$Tiempo, na.rm = TRUE)
sum_dia_trat <- sum(datos$Tiempo[datos$Dia == 4 & datos$Tratamiento == "E"], na.rm = TRUE)
sum_lote_trat <- sum(datos$Tiempo[datos$Lote == 3 & datos$Tratamiento == "E"], na.rm = TRUE)
sum_lote_dia <- sum(datos$Tiempo[datos$Lote == 3 & datos$Dia == 4], na.rm = TRUE)