library(tidyverse)
library(FrF2)
library(broom)
library(flextable)Diseño factorial completo \(2^k\) en R
Librerías
Al inicio de todas las clase se colocarán las librerías que se van a utilizar para el desarrollo de los casos.
Introducción
Diseñar un experimento factorial es bastante simple, y hacerlo en R no es la excepción. En particular, a mi me gusta usar el tidyverse para casi todo.
Por lo general se parte de una matriz de diseño y con base en ella se trabaja con variables codificadas.
Supongamos una situación en la que interesa diseñar un experimento factorial \(2^4\)$2^4$.
A <- B <- C <- D <- c(-1, 1)
diseno <- tidyr::expand_grid(A, B, C, D)
diseno# A tibble: 16 × 4
A B C D
<dbl> <dbl> <dbl> <dbl>
1 -1 -1 -1 -1
2 -1 -1 -1 1
3 -1 -1 1 -1
4 -1 -1 1 1
5 -1 1 -1 -1
6 -1 1 -1 1
7 -1 1 1 -1
8 -1 1 1 1
9 1 -1 -1 -1
10 1 -1 -1 1
11 1 -1 1 -1
12 1 -1 1 1
13 1 1 -1 -1
14 1 1 -1 1
15 1 1 1 -1
16 1 1 1 1Lo anterior corresponde a un diseño factorial \(2^4\), pero sin réplicas. Para replicarlo, se puede usar rbind() o alguna función compuesta que lo genere.
En este caso, vamos a replicarlo 2 veces.
diseno <- rbind(diseno, diseno)
diseno# A tibble: 32 × 4
A B C D
<dbl> <dbl> <dbl> <dbl>
1 -1 -1 -1 -1
2 -1 -1 -1 1
3 -1 -1 1 -1
4 -1 -1 1 1
5 -1 1 -1 -1
6 -1 1 -1 1
7 -1 1 1 -1
8 -1 1 1 1
9 1 -1 -1 -1
10 1 -1 -1 1
# ℹ 22 more rowsFinalmente, solo nos falta aleatorizarlo
diseno <- diseno[order(sample(1:32)), ]Con esto, tenemos un diseño factorial generado \(2^4\), con lo cual se puede expandir a otras cantidades de factores y niveles.
Este diseño también se puede relizar con liberías, por ejemplo usando FrF2.
FrF2::FrF2(nruns = 2^4,
nfactors = 4,
replications = 2,
blocks = 1,
randomize = TRUE) run.no run.no.std.rp A B C D Blocks
1 1 7.1 -1 1 1 -1 .1
2 2 11.1 -1 1 -1 1 .1
3 3 5.1 -1 -1 1 -1 .1
4 4 14.1 1 -1 1 1 .1
5 5 15.1 -1 1 1 1 .1
6 6 8.1 1 1 1 -1 .1
7 7 1.1 -1 -1 -1 -1 .1
8 8 2.1 1 -1 -1 -1 .1
9 9 10.1 1 -1 -1 1 .1
10 10 16.1 1 1 1 1 .1
11 11 4.1 1 1 -1 -1 .1
12 12 3.1 -1 1 -1 -1 .1
13 13 9.1 -1 -1 -1 1 .1
14 14 6.1 1 -1 1 -1 .1
15 15 13.1 -1 -1 1 1 .1
16 16 12.1 1 1 -1 1 .1
17 17 8.2 1 1 1 -1 .2
18 18 2.2 1 -1 -1 -1 .2
19 19 14.2 1 -1 1 1 .2
20 20 11.2 -1 1 -1 1 .2
21 21 16.2 1 1 1 1 .2
22 22 10.2 1 -1 -1 1 .2
23 23 15.2 -1 1 1 1 .2
24 24 5.2 -1 -1 1 -1 .2
25 25 1.2 -1 -1 -1 -1 .2
26 26 6.2 1 -1 1 -1 .2
27 27 12.2 1 1 -1 1 .2
28 28 13.2 -1 -1 1 1 .2
29 29 4.2 1 1 -1 -1 .2
30 30 7.2 -1 1 1 -1 .2
31 31 3.2 -1 1 -1 -1 .2
32 32 9.2 -1 -1 -1 1 .2
class=design, type= full factorial
NOTE: columns run.no and run.no.std.rp are annotation,
not part of the data frameNo obstante, este tipo de librarias, a título personal, no me encantan. ¿Pueden deducir por qué?
Trabajemos con el diseño creado con R base, agreguemos una respuesta ficticia, solo a modo de ejemplo.
R <- rnorm(32, 30, 2)
resultados <- cbind(diseno, R)Se procede con el análisis
modelo1 <- lm(R ~ A*B*C*D,
data = resultados)
broom::tidy(modelo1) %>%
dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>%
flextable::flextable()term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
(Intercept) | 29.58 | 0.34 | 86.55 | 0.00 |
A | -0.35 | 0.34 | -1.03 | 0.32 |
B | -0.03 | 0.34 | -0.09 | 0.93 |
C | -0.43 | 0.34 | -1.27 | 0.22 |
D | -0.73 | 0.34 | -2.15 | 0.05 |
A:B | 0.60 | 0.34 | 1.75 | 0.10 |
A:C | 0.21 | 0.34 | 0.63 | 0.54 |
B:C | 0.14 | 0.34 | 0.40 | 0.70 |
A:D | -0.14 | 0.34 | -0.40 | 0.69 |
B:D | -0.23 | 0.34 | -0.69 | 0.50 |
C:D | 0.29 | 0.34 | 0.86 | 0.40 |
A:B:C | -0.32 | 0.34 | -0.95 | 0.36 |
A:B:D | -0.39 | 0.34 | -1.14 | 0.27 |
A:C:D | 0.17 | 0.34 | 0.49 | 0.63 |
B:C:D | 0.35 | 0.34 | 1.02 | 0.32 |
A:B:C:D | 0.24 | 0.34 | 0.71 | 0.49 |
broom::tidy(aov(modelo1)) %>%
dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>%
flextable::flextable()term | df | sumsq | meansq | statistic | p.value |
|---|---|---|---|---|---|
A | 1 | 3.99 | 3.99 | 1.07 | 0.32 |
B | 1 | 0.03 | 0.03 | 0.01 | 0.93 |
C | 1 | 6.04 | 6.04 | 1.61 | 0.22 |
D | 1 | 17.20 | 17.20 | 4.60 | 0.05 |
A:B | 1 | 11.45 | 11.45 | 3.06 | 0.10 |
A:C | 1 | 1.46 | 1.46 | 0.39 | 0.54 |
B:C | 1 | 0.59 | 0.59 | 0.16 | 0.70 |
A:D | 1 | 0.61 | 0.61 | 0.16 | 0.69 |
B:D | 1 | 1.76 | 1.76 | 0.47 | 0.50 |
C:D | 1 | 2.76 | 2.76 | 0.74 | 0.40 |
A:B:C | 1 | 3.36 | 3.36 | 0.90 | 0.36 |
A:B:D | 1 | 4.87 | 4.87 | 1.30 | 0.27 |
A:C:D | 1 | 0.88 | 0.88 | 0.24 | 0.63 |
B:C:D | 1 | 3.90 | 3.90 | 1.04 | 0.32 |
A:B:C:D | 1 | 1.90 | 1.90 | 0.51 | 0.49 |
Residuals | 16 | 59.80 | 3.74 |
Algunos resultados gráficos:
FrF2::DanielPlot(modelo1) # gráfico de daniel
FrF2::MEPlot(modelo1) # Efectos principales
FrF2::IAPlot(modelo1) # Efectos de interacción
Para un gráfico de Pareto
efectos <- coef(modelo1)[-1] # quitamos intercepto
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.
modelo_final <- step(modelo1, direction = "backward")Start: AIC=52.01
R ~ A * B * C * D
Df Sum of Sq RSS AIC
- A:B:C:D 1 1.9048 61.702 51.011
<none> 59.797 52.007
Step: AIC=51.01
R ~ A + B + C + D + A:B + A:C + B:C + A:D + B:D + C:D + A:B:C +
A:B:D + A:C:D + B:C:D
Df Sum of Sq RSS AIC
- A:C:D 1 0.8826 62.585 49.465
- A:B:C 1 3.3577 65.060 50.706
- B:C:D 1 3.9018 65.604 50.973
<none> 61.702 51.011
- A:B:D 1 4.8668 66.569 51.440
Step: AIC=49.47
R ~ A + B + C + D + A:B + A:C + B:C + A:D + B:D + C:D + A:B:C +
A:B:D + B:C:D
Df Sum of Sq RSS AIC
- A:B:C 1 3.3577 65.942 49.137
- B:C:D 1 3.9018 66.487 49.400
<none> 62.585 49.465
- A:B:D 1 4.8668 67.451 49.862
Step: AIC=49.14
R ~ A + B + C + D + A:B + A:C + B:C + A:D + B:D + C:D + A:B:D +
B:C:D
Df Sum of Sq RSS AIC
- A:C 1 1.4611 67.403 47.839
- B:C:D 1 3.9018 69.844 48.977
<none> 65.942 49.137
- A:B:D 1 4.8668 70.809 49.416
Step: AIC=47.84
R ~ A + B + C + D + A:B + B:C + A:D + B:D + C:D + A:B:D + B:C:D
Df Sum of Sq RSS AIC
- B:C:D 1 3.9018 71.305 47.640
<none> 67.403 47.839
- A:B:D 1 4.8668 72.270 48.070
Step: AIC=47.64
R ~ A + B + C + D + A:B + B:C + A:D + B:D + C:D + A:B:D
Df Sum of Sq RSS AIC
- B:C 1 0.5892 71.895 45.903
- C:D 1 2.7622 74.068 46.856
<none> 71.305 47.640
- A:B:D 1 4.8668 76.172 47.752
Step: AIC=45.9
R ~ A + B + C + D + A:B + A:D + B:D + C:D + A:B:D
Df Sum of Sq RSS AIC
- C:D 1 2.7622 74.657 45.109
<none> 71.895 45.903
- A:B:D 1 4.8668 76.761 45.999
Step: AIC=45.11
R ~ A + B + C + D + A:B + A:D + B:D + A:B:D
Df Sum of Sq RSS AIC
<none> 74.657 45.109
- A:B:D 1 4.8668 79.524 45.130
- C 1 6.0355 80.692 45.597broom::tidy(modelo_final) %>%
dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>%
flextable::flextable()term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
(Intercept) | 29.58 | 0.32 | 92.87 | 0.00 |
A | -0.35 | 0.32 | -1.11 | 0.28 |
B | -0.03 | 0.32 | -0.10 | 0.92 |
C | -0.43 | 0.32 | -1.36 | 0.19 |
D | -0.73 | 0.32 | -2.30 | 0.03 |
A:B | 0.60 | 0.32 | 1.88 | 0.07 |
A:D | -0.14 | 0.32 | -0.43 | 0.67 |
B:D | -0.23 | 0.32 | -0.74 | 0.47 |
A:B:D | -0.39 | 0.32 | -1.22 | 0.23 |
broom::tidy(aov(modelo_final)) %>%
dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>%
flextable::flextable()term | df | sumsq | meansq | statistic | p.value |
|---|---|---|---|---|---|
A | 1 | 3.99 | 3.99 | 1.23 | 0.28 |
B | 1 | 0.03 | 0.03 | 0.01 | 0.92 |
C | 1 | 6.04 | 6.04 | 1.86 | 0.19 |
D | 1 | 17.20 | 17.20 | 5.30 | 0.03 |
A:B | 1 | 11.45 | 11.45 | 3.53 | 0.07 |
A:D | 1 | 0.61 | 0.61 | 0.19 | 0.67 |
B:D | 1 | 1.76 | 1.76 | 0.54 | 0.47 |
A:B:D | 1 | 4.87 | 4.87 | 1.50 | 0.23 |
Residuals | 23 | 74.66 | 3.25 |