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 10.1 1 -1 -1 1 .1
2 2 13.1 -1 -1 1 1 .1
3 3 7.1 -1 1 1 -1 .1
4 4 2.1 1 -1 -1 -1 .1
5 5 1.1 -1 -1 -1 -1 .1
6 6 9.1 -1 -1 -1 1 .1
7 7 15.1 -1 1 1 1 .1
8 8 5.1 -1 -1 1 -1 .1
9 9 14.1 1 -1 1 1 .1
10 10 6.1 1 -1 1 -1 .1
11 11 16.1 1 1 1 1 .1
12 12 8.1 1 1 1 -1 .1
13 13 12.1 1 1 -1 1 .1
14 14 3.1 -1 1 -1 -1 .1
15 15 11.1 -1 1 -1 1 .1
16 16 4.1 1 1 -1 -1 .1
17 17 6.2 1 -1 1 -1 .2
18 18 1.2 -1 -1 -1 -1 .2
19 19 7.2 -1 1 1 -1 .2
20 20 11.2 -1 1 -1 1 .2
21 21 12.2 1 1 -1 1 .2
22 22 10.2 1 -1 -1 1 .2
23 23 5.2 -1 -1 1 -1 .2
24 24 13.2 -1 -1 1 1 .2
25 25 14.2 1 -1 1 1 .2
26 26 16.2 1 1 1 1 .2
27 27 15.2 -1 1 1 1 .2
28 28 2.2 1 -1 -1 -1 .2
29 29 8.2 1 1 1 -1 .2
30 30 4.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 librerias, 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) | 30.64 | 0.31 | 97.75 | 0.00 |
A | 0.28 | 0.31 | 0.88 | 0.39 |
B | 0.22 | 0.31 | 0.69 | 0.50 |
C | 0.06 | 0.31 | 0.19 | 0.85 |
D | 0.09 | 0.31 | 0.27 | 0.79 |
A:B | -0.10 | 0.31 | -0.33 | 0.74 |
A:C | 0.18 | 0.31 | 0.58 | 0.57 |
B:C | 0.42 | 0.31 | 1.33 | 0.20 |
A:D | 0.78 | 0.31 | 2.50 | 0.02 |
B:D | -0.37 | 0.31 | -1.18 | 0.26 |
C:D | -0.28 | 0.31 | -0.89 | 0.38 |
A:B:C | -0.01 | 0.31 | -0.02 | 0.98 |
A:B:D | -0.26 | 0.31 | -0.82 | 0.43 |
A:C:D | -0.18 | 0.31 | -0.58 | 0.57 |
B:C:D | -0.08 | 0.31 | -0.26 | 0.80 |
A:B:C:D | -0.20 | 0.31 | -0.63 | 0.54 |
broom::tidy(aov(modelo1)) %>%
dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>%
flextable::flextable()term | df | sumsq | meansq | statistic | p.value |
|---|---|---|---|---|---|
A | 1 | 2.46 | 2.46 | 0.78 | 0.39 |
B | 1 | 1.49 | 1.49 | 0.47 | 0.50 |
C | 1 | 0.12 | 0.12 | 0.04 | 0.85 |
D | 1 | 0.23 | 0.23 | 0.07 | 0.79 |
A:B | 1 | 0.35 | 0.35 | 0.11 | 0.74 |
A:C | 1 | 1.07 | 1.07 | 0.34 | 0.57 |
B:C | 1 | 5.55 | 5.55 | 1.76 | 0.20 |
A:D | 1 | 19.67 | 19.67 | 6.26 | 0.02 |
B:D | 1 | 4.36 | 4.36 | 1.39 | 0.26 |
C:D | 1 | 2.51 | 2.51 | 0.80 | 0.38 |
A:B:C | 1 | 0.00 | 0.00 | 0.00 | 0.98 |
A:B:D | 1 | 2.09 | 2.09 | 0.67 | 0.43 |
A:C:D | 1 | 1.07 | 1.07 | 0.34 | 0.57 |
B:C:D | 1 | 0.21 | 0.21 | 0.07 | 0.80 |
A:B:C:D | 1 | 1.25 | 1.25 | 0.40 | 0.54 |
Residuals | 16 | 50.31 | 3.14 |
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=46.48
R ~ A * B * C * D
Df Sum of Sq RSS AIC
- A:B:C:D 1 1.2525 51.566 45.268
<none> 50.313 46.481
Step: AIC=45.27
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:B:C 1 0.00129 51.567 43.269
- B:C:D 1 0.20979 51.776 43.398
- A:C:D 1 1.07337 52.639 43.927
- A:B:D 1 2.09276 53.659 44.541
<none> 51.566 45.268
Step: AIC=43.27
R ~ A + B + C + D + A:B + A:C + B:C + A:D + B:D + C:D + A:B:D +
A:C:D + B:C:D
Df Sum of Sq RSS AIC
- B:C:D 1 0.20979 51.777 41.399
- A:C:D 1 1.07337 52.640 41.928
- A:B:D 1 2.09276 53.660 42.542
<none> 51.567 43.269
Step: AIC=41.4
R ~ A + B + C + D + A:B + A:C + B:C + A:D + B:D + C:D + A:B:D +
A:C:D
Df Sum of Sq RSS AIC
- A:C:D 1 1.0734 52.850 40.055
- A:B:D 1 2.0928 53.870 40.667
<none> 51.777 41.399
- B:C 1 5.5479 57.325 42.656
Step: AIC=40.06
R ~ A + B + C + D + A:B + A:C + B:C + A:D + B:D + C:D + A:B:D
Df Sum of Sq RSS AIC
- A:C 1 1.0717 53.922 38.698
- A:B:D 1 2.0928 54.943 39.298
- C:D 1 2.5107 55.361 39.540
<none> 52.850 40.055
- B:C 1 5.5479 58.398 41.250
Step: AIC=38.7
R ~ A + B + C + D + A:B + B:C + A:D + B:D + C:D + A:B:D
Df Sum of Sq RSS AIC
- A:B:D 1 2.0928 56.015 37.916
- C:D 1 2.5107 56.433 38.154
<none> 53.922 38.698
- B:C 1 5.5479 59.470 39.832
Step: AIC=37.92
R ~ A + B + C + D + A:B + B:C + A:D + B:D + C:D
Df Sum of Sq RSS AIC
- A:B 1 0.3473 56.362 36.114
- C:D 1 2.5107 58.525 37.319
<none> 56.015 37.916
- B:D 1 4.3596 60.374 38.315
- B:C 1 5.5479 61.563 38.938
- A:D 1 19.6708 75.686 45.547
Step: AIC=36.11
R ~ A + B + C + D + B:C + A:D + B:D + C:D
Df Sum of Sq RSS AIC
- C:D 1 2.5107 58.873 35.509
<none> 56.362 36.114
- B:D 1 4.3596 60.722 36.498
- B:C 1 5.5479 61.910 37.118
- A:D 1 19.6708 76.033 43.694
Step: AIC=35.51
R ~ A + B + C + D + B:C + A:D + B:D
Df Sum of Sq RSS AIC
<none> 58.873 35.509
- B:D 1 4.3596 63.232 35.795
- B:C 1 5.5479 64.421 36.390
- A:D 1 19.6708 78.544 42.733broom::tidy(modelo_final) %>%
dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>%
flextable::flextable()term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
(Intercept) | 30.64 | 0.28 | 110.67 | 0.00 |
A | 0.28 | 0.28 | 1.00 | 0.33 |
B | 0.22 | 0.28 | 0.78 | 0.44 |
C | 0.06 | 0.28 | 0.22 | 0.83 |
D | 0.09 | 0.28 | 0.31 | 0.76 |
B:C | 0.42 | 0.28 | 1.50 | 0.15 |
A:D | 0.78 | 0.28 | 2.83 | 0.01 |
B:D | -0.37 | 0.28 | -1.33 | 0.20 |
broom::tidy(aov(modelo_final)) %>%
dplyr::mutate(across(where(is.numeric), ~ round(., 2))) %>%
flextable::flextable()term | df | sumsq | meansq | statistic | p.value |
|---|---|---|---|---|---|
A | 1 | 2.46 | 2.46 | 1.00 | 0.33 |
B | 1 | 1.49 | 1.49 | 0.61 | 0.44 |
C | 1 | 0.12 | 0.12 | 0.05 | 0.83 |
D | 1 | 0.23 | 0.23 | 0.10 | 0.76 |
B:C | 1 | 5.55 | 5.55 | 2.26 | 0.15 |
A:D | 1 | 19.67 | 19.67 | 8.02 | 0.01 |
B:D | 1 | 4.36 | 4.36 | 1.78 | 0.20 |
Residuals | 24 | 58.87 | 2.45 |