Diseño factorial completo \(2^{5}\) en \(2^{1}\) bloques en R

Librerías

library(tidyverse)
library(FrF2)
library(broom)
library(flextable)

Introducción

Para realizar un diseño factorial completo \(2^5\) (32 combinaciones) dispuesto en bloques, podemos utilizar la interacción de mayor orden (\(ABCDE\)) como contraste definidor de bloqueo.

Al confundir \(ABCDE\) con los bloques, sacrificamos la estimación de esta interacción de 5 factores a cambio de aislar la variabilidad por bloques. En este ejemplo, replicamos el experimento 2 veces (\(N = 64\) corridas) distribuidas en 4 bloques de tamaño 16.

A <- B <- C <- D <- E <- c(-1, 1)

# Base completa 2^5
diseno_base <- tidyr::expand_grid(A, B, C, D, E)

# Experimento replicado con confusión de ABCDE
diseno <- diseno_base |> 
  rbind(diseno_base) |> 
  dplyr::mutate(ABCDE = A*B*C*D*E) |> 
  dplyr::arrange(ABCDE) |> 
  dplyr::mutate(Bloque = factor(c(rep(1, 16), 
                                  rep(3, 16), 
                                  rep(2, 16), 
                                  rep(4, 16)))) |> 
  dplyr::select(-ABCDE) # Se elimina el contraste tras asignar bloques

diseno

# A tibble: 64 × 6
       A     B     C     D     E Bloque
   <dbl> <dbl> <dbl> <dbl> <dbl> <fct> 
 1    -1    -1    -1    -1    -1 1     
 2    -1    -1    -1     1     1 1     
 3    -1    -1     1    -1     1 1     
 4    -1    -1     1     1    -1 1     
 5    -1     1    -1    -1     1 1     
 6    -1     1    -1     1    -1 1     
 7    -1     1     1    -1    -1 1     
 8    -1     1     1     1     1 1     
 9     1    -1    -1    -1     1 1     
10     1    -1    -1     1    -1 1     
# ℹ 54 more rows

Recuerde que el diseño se ejecuta en orden aleatorio, en este caso, intrabloque.

# Aleatorización intrabloque
set.seed(123)
diseno <- diseno |> 
  dplyr::group_by(Bloque) |> 
  dplyr::slice_sample(prop = 1) |> 
  dplyr::ungroup()

diseno

# A tibble: 64 × 6
       A     B     C     D     E Bloque
   <dbl> <dbl> <dbl> <dbl> <dbl> <fct> 
 1     1     1     1    -1     1 1     
 2     1     1     1     1    -1 1     
 3    -1    -1     1    -1     1 1     
 4     1     1    -1     1     1 1     
 5     1    -1    -1     1    -1 1     
 6    -1    -1    -1     1     1 1     
 7    -1     1    -1     1    -1 1     
 8    -1     1    -1    -1     1 1     
 9    -1    -1     1     1    -1 1     
10     1    -1     1     1     1 1     
# ℹ 54 more rows

Análisis del modelo

Trabajemos con el diseño creado con R base, agreguemos una respuesta ficticia, solo a modo de ejemplo.

R <- rnorm(64, 30, 2)

resultados <- cbind(diseno, R)

Ajustamos el modelo saturado. Al estar Bloque definido como factor (3 GL), se absorbe la variación correspondiente a \(ABCDE\), por lo que la interacción A:B:C:D:E resulta inestimable (aparece como NA / 0 GL), demostrando la confusión deseada.

modelo1 <- lm(R ~ Bloque + A*B*C*D*E, 
              data = resultados)

broom::tidy(modelo1) |> 
  dplyr::mutate(across(where(is.numeric), ~ round(., 2))) |> 
  flextable::flextable()

term	estimate	std.error	statistic	p.value
(Intercept)	30.10	0.40	75.44	0.00
Bloque2	0.14	0.56	0.25	0.81
Bloque3	-0.45	0.56	-0.80	0.43
Bloque4	0.69	0.56	1.22	0.23
A	-0.11	0.20	-0.53	0.60
B	0.09	0.20	0.46	0.65
C	-0.15	0.20	-0.75	0.46
D	0.35	0.20	1.77	0.09
E	0.00	0.20	0.01	0.99
A:B	-0.13	0.20	-0.63	0.53
A:C	0.29	0.20	1.46	0.15
B:C	-0.21	0.20	-1.06	0.30
A:D	-0.03	0.20	-0.13	0.90
B:D	-0.06	0.20	-0.28	0.78
C:D	-0.14	0.20	-0.69	0.50
A:E	0.47	0.20	2.37	0.02
B:E	-0.43	0.20	-2.15	0.04
C:E	0.19	0.20	0.98	0.34
D:E	-0.19	0.20	-0.95	0.35
A:B:C	-0.18	0.20	-0.90	0.38
A:B:D	0.25	0.20	1.23	0.23
A:C:D	0.14	0.20	0.71	0.48
B:C:D	0.15	0.20	0.73	0.47
A:B:E	0.27	0.20	1.35	0.19
A:C:E	0.12	0.20	0.59	0.56
B:C:E	0.27	0.20	1.35	0.19
A:D:E	-0.66	0.20	-3.32	0.00
B:D:E	0.40	0.20	2.01	0.05
C:D:E	-0.24	0.20	-1.18	0.25
A:B:C:D	0.17	0.20	0.85	0.40
A:B:C:E	0.05	0.20	0.27	0.79
A:B:D:E	0.06	0.20	0.30	0.76
A:C:D:E	-0.10	0.20	-0.49	0.63
B:C:D:E	0.29	0.20	1.46	0.15
A:B:C:D:E

broom::tidy(aov(modelo1)) %>% 
  dplyr::mutate(across(where(is.numeric), ~ round(., 2))) |>  
  flextable::flextable()

term	df	sumsq	meansq	statistic	p.value
Bloque	3	10.60	3.53	1.39	0.27
A	1	0.71	0.71	0.28	0.60
B	1	0.54	0.54	0.21	0.65
C	1	1.43	1.43	0.56	0.46
D	1	8.02	8.02	3.15	0.09
E	1	0.00	0.00	0.00	0.99
A:B	1	1.02	1.02	0.40	0.53
A:C	1	5.44	5.44	2.13	0.15
B:C	1	2.87	2.87	1.13	0.30
A:D	1	0.04	0.04	0.02	0.90
B:D	1	0.20	0.20	0.08	0.78
C:D	1	1.20	1.20	0.47	0.50
A:E	1	14.25	14.25	5.60	0.02
B:E	1	11.75	11.75	4.61	0.04
C:E	1	2.43	2.43	0.95	0.34
D:E	1	2.32	2.32	0.91	0.35
A:B:C	1	2.06	2.06	0.81	0.38
A:B:D	1	3.85	3.85	1.51	0.23
A:C:D	1	1.29	1.29	0.51	0.48
B:C:D	1	1.36	1.36	0.53	0.47
A:B:E	1	4.67	4.67	1.83	0.19
A:C:E	1	0.87	0.87	0.34	0.56
B:C:E	1	4.62	4.62	1.82	0.19
A:D:E	1	28.03	28.03	11.01	0.00
B:D:E	1	10.25	10.25	4.02	0.05
C:D:E	1	3.54	3.54	1.39	0.25
A:B:C:D	1	1.85	1.85	0.72	0.40
A:B:C:E	1	0.18	0.18	0.07	0.79
A:B:D:E	1	0.24	0.24	0.09	0.76
A:C:D:E	1	0.61	0.61	0.24	0.63
B:C:D:E	1	5.45	5.45	2.14	0.15
Residuals	30	76.40	2.55

Algunos resultados gráficos:

# Modelo auxiliar sin el factor Bloque (solo los 5 factores de 2 niveles)
# porque las funciones no dejan más de dos niveles
modelo_graficos <- lm(R ~ A*B*C*D*E, data = resultados)

FrF2::DanielPlot(modelo_graficos) # gráfico de daniel

FrF2::MEPlot(modelo_graficos) # Efectos principales

FrF2::IAPlot(modelo_graficos) # Efectos de interacción

Para un gráfico de Pareto

efectos <- coef(modelo1)[-c(1, 2, 3, 4, 35)]  # quitamos intercepto, bloques y ABCDE
df_efectos <- data.frame(
  factor = names(efectos),
  efecto = abs(efectos))

ggplot(df_efectos, aes(x = reorder(factor, efecto), y = efecto)) +
  geom_bar(stat = "identity", fill =  "darksalmon") +
  coord_flip() +
  labs(x = "Factor/Interacción", y = "Magnitud del efecto",
       title = "Pareto de efectos factoriales") +
  theme_bw()

¿Qué nos falta?

Pues eliminar todo aquello que no aporte al modelo, esta es una forma de hacerlo automático, pero en DdE puede ser más valioso hacerlo paso a paso.

#forzando que el bloque siempre sea parte del modelo
modelo_final <- step(modelo1, 
                     direction = "backward", 
                     scope = list(lower = ~ Bloque, 
                                  upper = formula(modelo1)), 
                     trace = 0)

broom::tidy(modelo_final) %>% 
  dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>% 
  flextable::flextable()

term	estimate	std.error	statistic	p.value
(Intercept)	30.10	0.38	80.14	0.00
Bloque2	0.14	0.53	0.26	0.79
Bloque3	-0.45	0.53	-0.85	0.40
Bloque4	0.69	0.53	1.29	0.20
A	-0.11	0.19	-0.56	0.58
B	0.09	0.19	0.49	0.63
C	-0.15	0.19	-0.79	0.43
D	0.35	0.19	1.88	0.07
E	0.00	0.19	0.01	0.99
A:B	-0.13	0.19	-0.67	0.51
A:C	0.29	0.19	1.55	0.13
B:C	-0.21	0.19	-1.13	0.27
A:D	-0.03	0.19	-0.14	0.89
B:D	-0.06	0.19	-0.30	0.77
C:D	-0.14	0.19	-0.73	0.47
A:E	0.47	0.19	2.51	0.02
B:E	-0.43	0.19	-2.28	0.03
C:E	0.19	0.19	1.04	0.31
D:E	-0.19	0.19	-1.01	0.32
A:B:D	0.25	0.19	1.31	0.20
B:C:D	0.15	0.19	0.78	0.44
A:B:E	0.27	0.19	1.44	0.16
B:C:E	0.27	0.19	1.43	0.16
A:D:E	-0.66	0.19	-3.52	0.00
B:D:E	0.40	0.19	2.13	0.04
C:D:E	-0.24	0.19	-1.25	0.22
B:C:D:E	0.29	0.19	1.55	0.13

broom::tidy(aov(modelo_final)) %>% 
  dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>% 
  flextable::flextable()

term	df	sumsq	meansq	statistic	p.value
Bloque	3	10.60	3.53	1.57	0.21
A	1	0.71	0.71	0.32	0.58
B	1	0.54	0.54	0.24	0.63
C	1	1.43	1.43	0.63	0.43
D	1	8.02	8.02	3.55	0.07
E	1	0.00	0.00	0.00	0.99
A:B	1	1.02	1.02	0.45	0.51
A:C	1	5.44	5.44	2.41	0.13
B:C	1	2.87	2.87	1.27	0.27
A:D	1	0.04	0.04	0.02	0.89
B:D	1	0.20	0.20	0.09	0.77
C:D	1	1.20	1.20	0.53	0.47
A:E	1	14.25	14.25	6.31	0.02
B:E	1	11.75	11.75	5.20	0.03
C:E	1	2.43	2.43	1.08	0.31
D:E	1	2.32	2.32	1.03	0.32
A:B:D	1	3.85	3.85	1.71	0.20
B:C:D	1	1.36	1.36	0.60	0.44
A:B:E	1	4.67	4.67	2.07	0.16
B:C:E	1	4.62	4.62	2.05	0.16
A:D:E	1	28.03	28.03	12.42	0.00
B:D:E	1	10.25	10.25	4.54	0.04
C:D:E	1	3.54	3.54	1.57	0.22
B:C:D:E	1	5.45	5.45	2.41	0.13
Residuals	37	83.50	2.26