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.
Se utiliza cuando hay k factores, cada uno con 2 niveles. El número total de tratamientos es:
\[ n = 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:
Se descompone la variación total:
\[ \text{SCT} = \text{SCA} + \text{SCB} + \text{SCAB} + \text{SCE} \]
Donde:
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\) |
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\) |
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} \]
# 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.
Diseño factorial \(A \times B\) con efectos aleatorios en ambos factores.
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í.
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\) |
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.
Los componentes de varianza se estiman como:
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.
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í.
Para una observación individual:
\[ \operatorname{Var}(Y_{ijk}) = \sigma^2 + \sigma^2_{\alpha\beta} + \sigma^2_\beta \]
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 \]
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.
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\) |
Los componentes se estiman a partir de los cuadrados medios:
# 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)