Diseños Factoriales \(2^k\)
Cuando hay multiple factores bajo estudio, el uso de diseños con un solo factor no son eficientes y pueden haber situaciones donde lleven a resultados equivocados. Una mejor estrategia es usar un diseño factorial.
En un diseño factorial las celdas consisten de todas las posibles combinaciones de los niveles de los factores bajo estudio. Y el objetivo principal es estudiar el efecto de varios factores sobre una o varias respuestas, cuando se tiene el mismo interés sobre todos los factores y existe la posibilidad de que el efecto de un factor cambie según los niveles de otros factores, esto es, que los factores interactúen, o exista interacción.
El efecto de un factor es definido como el cambio en la respuesta producido por un cambio en el nivel del factor.
Empesamos estudiando el diseño factorial \(2^k\), es decir, el diseño factorial en el que hay \(k\) factores y todos tienen dos niveles, ya que son ampliamente usados en la practica, su interpretación es sencilla, y sirven para explorar el espacio de los factores.
Ejemplo – Reduction of an enamine
(Carlson R. and Carlson J.E., 2005). Enamines can be reduced to the corresponding saturated amine by treatment with formic acid. To determine suitabile experimental conditions for the reduction of enamines derived from camphor (\(\mathrm{C_{10}H_{16}O}\)), a factorial design was used to explore the reaction of the morpholine enamine. The reaction was conducted by adding formic acid drop-wise to the neat enamine at such a rate that the foaming caused by the evolution of carbon dioxide could be kept under control. The reaction is rapid and completed within a few minutes. The main product in the reaction was the desired bornylmorpholine (mixture of endo and exo isomers). There were also varying amounts of unreacted enamine and some camphor formed by hydrolysis of the enamine.
The experimental procedure is very simple and there are only two variables to consider \(X_1\), the amount of formic acid (1,1.5 mol/mol), and \(X_2\), the reaction temperature (25,100\(^{\circ}C\)). When the reaction was complete, the resulting mixture was treated with aqueous sodium hydroxide. Then the composition of the organic layer was analysed by gas chromatography.
The recovery of organic material was quantitative. The measured responses were: \(Y_1:\) the yield of bornylmorpholine (%), \(Y_2:\) the amount of unreacted enamine (%), and \(Y_3:\) the amount of camphor (%).
| Exp. no. | \(X_1\) | \(X_2\) | \(Y_1\) | Etiqueta |
|---|---|---|---|---|
| 1 | - | - | 80.4 | (1) |
| 2 | + | - | 72.4 | a |
| 3 | - | + | 94.4 | b |
| 4 | + | + | 90.6 | ab |
Note que el nivel bajo de los factores se denota usualmente por
- y el alto por +.
Calculemos los efectos para \(Y_1\):
| . | Nivel bajo \(X_2\) | Nivel alto \(X_2\) | Promedio |
|---|---|---|---|
| Efecto de \(X_1\) | \(\frac{a}{n}-\frac{(1)}{n} = 72.4-80.4= -8\) | \(\frac{ab}{n}-\frac{b}{n} =90.6-94.4 = -3.8\) | \((-8-3.8)/2 = -5.9\) |
Luego aumentar la cantidad de ácido fórmico causa una disminución de 5.9 % de la producción de bornilmorfolina.
Similarmente para \(X_2:\)
| . | Nivel bajo \(X_1\) | Nivel alto \(X_1\) | Promedio |
|---|---|---|---|
| Efecto de \(X_2\) | \(94.4-80.4= 14\) | \(90.6-72.4 = 18.2\) | \((14+18.2)/2 = 16.1\) |
Luego aumentar la temperatura causa un aumento del 16.1% de la producción de bornilmorfolina.
Los efectos principales se puede usar el siguiente código en R:
enamine <- data.frame(x1=gl(2,1,4,labels=c("-","+")),
x2=gl(2,2,4,labels=c("-","+")),
y=c(80.4,72.4, 94.4, 90.6))
coded <- function(x,v="-"){
if(missing(v))
out= x==x[1]
else
out= x==v
return(2*out-1)
}
sum(coded(enamine$x1)*enamine$y)/2## [1] 5.9sum(coded(enamine$x2)*enamine$y)/2## [1] -16.1sum(coded(enamine$x1)*coded(enamine$x2)*enamine$y)/2## [1] 2.1Ventajas del diseño factorial
Al examinar todas las posibles combinaciones de los niveles de los factores, el número de replicas de un nivel específico de un factor es incrementado por el producto del número de niveles de todos los otros factores en el diseño, y por lo tanto se puede obtener la misma potencia con menos réplicas.
Otra ventaja que tiene el uso de un diseño factorial sobre el uso de diseños con un solo factor es que permite el estudio de interacciones entre factores, es decir, se puede detectar situaciones en las que el efecto de un factor cambia dependiendo el nivel de otro factor.
Ejemplo – Reduction of an enamine - cont.
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsenamine %>%
ggplot() +
aes(x = x1, y = y, color = x2) +
geom_line(aes(group = x2)) +
geom_point()+
labs(title = "Gráfico de interacción",
x = "Amount of Formic acid",
y = "Yield of bornylmorpholine %",
color = "Temperature") +
scale_x_discrete(labels = c("1 mol/mol","1.5 mol/mol"))
Ejemplo – Proceso químico
(Montgomery, 2013). An investigation into the effect of the concentration of the reactant and the amount of the catalyst on the conversion (yield) in a chemical process. The objective of the experiment was to determine if adjustments to either of these two factors would increase the yield. Let the reactant concentration be factor A and let the two levels of interest be 15 and 25 percent. The catalyst is factor B, with the high level denoting the use of 2 pounds of the catalyst and the low level denoting the use of only 1 pound. The experiment is replicated three times, so there are 12 runs.
| Factor A | Factor B | Replicas | Total |
|---|---|---|---|
| - | - | 28 25 27 | 80 |
| + | - | 36 32 32 | 100 |
| - | + | 18 19 23 | 60 |
| + | + | 31 30 29 | 90 |
conversion <- data.frame(A=gl(2,3,12,c("-","+")),
B=gl(2,6,12,c("-","+")),
Yield=c(28,25,27,36,32,32,18,19,23,31,30,29) )
n <- 3
sum(coded(conversion$A)*conversion$Yield)/(2*n)## [1] -8.333333sum(coded(conversion$B)*conversion$Yield)/(2*n)## [1] 5sum(coded(conversion$A)*coded(conversion$B)*conversion$Yield)/(2*n)## [1] 1.666667Diseño completamente aleatorizado con dos factores
El modelo estadístico para un diseño completamente aleatorizado con dos factores puede ser escrito así:
\[y_{ijk} = \mu_{ij} + \epsilon_{ijk}\]
donde \(i\) representa el nivel del primer factor, \(j\) representa el nivel del segundo factor, y \(k\) representa el número de réplica. Éste modelo es llamado el modelo de celdas medias y \(\mu_{ij}\) representa la respuesta esperada en la \(ij\)-ésima celda. Otra manera de representar este modelo es con el modelo de efectos:
\[y_{ijk} = \mu + \alpha_i + \beta_j + \alpha\beta_{ij} + \epsilon_{ijk}\]
En este modelo \(\alpha_i\) y \(\beta_j\) son los efectos principales y representan la diferencia entre el promedio marginal de todos los experimentos en el \(i\)-ésimo nivel del primer factor y la media global, y la diferencia entre el promedio marginal de todos los experimentos en el \(j\)-ésimo nivel del segundo factor y la media global, respectivamente. El efecto de la interacción, \(\alpha\beta_{ij}\), representa la diferencia entre la media de la celda \(\mu_{ij}\) y \(\mu + \alpha_i + \beta_j\).
Análisis de varianza
Los supuestos de este modelo son: los errores experimentales son independientes y \(\epsilon_{ijk}\sim N(0,\sigma^2)\). El supuesto de independencia se garantiza si la combinación de tratamientos es asignada aleatoriamente a las unidades experimentales, y los supuestos de igualdad de varianza y normalidad pueden ser verificados como se ha hecho antes.
El análisis de varianza es usado para particonar la variabilidad de la variable respuesta debida a distintas fuentes. En el caso de un modelo de dos factores, la variabilidad se particiona en cuatro fuentes de variación: factor A, factor B, la interacción entre A y B, y el error. Es decir,
\[SST = SSA + SSB + SSAB + SSE,\]
donde el termino de interacción mide la no aditividad de los tratamientos.
| Fuente de variación | Grados de libertad | Suma de Cuadrados | MS | F |
|---|---|---|---|---|
| A | \(a-1\) | \(r b\displaystyle\sum_{i=1}^a (\overline{y}_{i\cdot\cdot} - \overline{\overline{y}})^2\) | \(\frac{SSA}{a-1}\) | \(\frac{MSA}{MSE}\) |
| B | \(b-1\) | \(r a\displaystyle\sum_{j=1}^b (\overline{y}_{\cdot j\cdot} - \overline{\overline{y}})^2\) | \(\frac{SSB}{b-1}\) | \(\frac{MSB}{MSE}\) |
| AB | \((a-1)(b-1)\) | \(r\displaystyle\sum_{i=1}^a\sum_{j=1}^b (\overline{y}_{ij\cdot} -\overline{y}_{i\cdot\cdot}-\overline{y}_{\cdot j\cdot}+ \overline{\overline{y}})^2\) | \(\frac{SSAB}{(a-1)(b-1)}\) | \(\frac{MSAB}{MSE}\) |
| Error | \(ab(r-1)\) | \(\displaystyle\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^{r} (y_{ijk} - \overline{y}_{ij\cdot})^2\) | \(\frac{SSE}{ab(r-1)}\) |
En R se puede estimar el modelo usando las funciones aov
o lm
conversion.fit <- lm(Yield ~ A * B, conversion)
anova(conversion.fit)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 208.333 208.333 53.1915 8.444e-05 ***
## B 1 75.000 75.000 19.1489 0.002362 **
## A:B 1 8.333 8.333 2.1277 0.182776
## Residuals 8 31.333 3.917
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1with(conversion, interaction.plot(A, B, Yield))
library(phia)## Loading required package: car## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## someplot(interactionMeans(conversion.fit))
Modelo de regresión
En un análisis factorial \(2^k\), es fácil expresar los resultados de un experimento en terminos de un modelo de regresión.
Ejemplo - Proceso químico - cont.
El modelo de regresión para el ejemplo del proceso químico (Montgomery, 2013) es
\[y = \beta_0 + \beta_1x_1 + \beta_2x_2+ \epsilon\]
donde \(x_1 = \frac{\mbox{conc} - (\mbox{conc}_{+}+\mbox{conc}_{-})/2}{(\mbox{conc}_{+}-\mbox{conc}_{-})/2}\) y \(x_2 = \frac{\mbox{cat} - (\mbox{cat}_{+}+\mbox{cat}_{-})/2}{(\mbox{cat}_{+}-\mbox{cat}_{-})/2}.\)
conversion2 <- data.frame(conc = rep(c(15,25),4, each=3),
cat=rep(c(1,2),each=6),
Yield=c(28,25,27,36,32,32,18,19,23,31,30,29) )
conversion2$A <- (conversion2$conc - 20)/5
conversion2$B <- (conversion2$cat - 1.5)/.5
conversion2.fit<-lm(Yield~A+B,data=conversion2) # Note que no se incluye interacción
summary(conversion2.fit)##
## Call:
## lm(formula = Yield ~ A + B, data = conversion2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8333 -1.9167 0.3333 1.8333 2.1667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.5000 0.3967 69.314 < 2e-16 ***
## A 4.1667 0.3967 10.502 8.16e-10 ***
## B -2.5000 0.3967 -6.301 3.00e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.944 on 21 degrees of freedom
## Multiple R-squared: 0.8772, Adjusted R-squared: 0.8655
## F-statistic: 75 on 2 and 21 DF, p-value: 2.734e-10conversion2.fit2<-lm(Yield~conc+cat,data=conversion2)
summary(conversion2.fit2)##
## Call:
## lm(formula = Yield ~ conc + cat, data = conversion2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8333 -1.9167 0.3333 1.8333 2.1667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.33333 2.02302 9.062 1.06e-08 ***
## conc 0.83333 0.07935 10.502 8.16e-10 ***
## cat -5.00000 0.79349 -6.301 3.00e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.944 on 21 degrees of freedom
## Multiple R-squared: 0.8772, Adjusted R-squared: 0.8655
## F-statistic: 75 on 2 and 21 DF, p-value: 2.734e-10El modelo ajustado en terminos de las variables codificadas es
\[\widehat{y} = 27.5 + \frac{8.33}{2}x_1 + \frac{-5.00}{2}x_2\]
y en terminos de las variables originales es
\[\widehat{y} = 18.333 + 0.833\cdot\mbox{conc} -5\cdot\mbox{cat}\]
Superficie de respuesta
library(lattice)
tmp <- list(conc=seq(15,25,by=.5),cat=seq(1,2,by=.1))
new.data <- expand.grid(tmp)
new.data$fit <- predict(conversion2.fit2,new.data)
contourplot(fit~conc+cat,new.data,xlab="Reactivo", ylab="Catalisador")
Ejemplo de un diseño factorial \(2^3\), “Fly ash”
(Lawson J. y Erjavec J., 2017). An experiment was set up to study the compressive strength of concrete made from fly ash from two different sources. The experiment factors chosen for study were:
- Fly Ash Source: The composition of the ash in different coals can be quite different. Hence, it is conceivable that the strength of the concrete could vary greatly, depending on the source of the ash. Therefore, it was important to investigate this as one of the factors.
- Water/Cement Ratio: It is important to add the appropriate amount of water, recognizing that what is appropriate may depend on the particular ash used (and its composition).
- Set-up Temperature: Since concrete is generally poured outdoors, it is important to find out if the strength of the concrete depends on the temperature at which it hardens.
The response is the 14 Day Compressive strength, mPa.
unos <- c(-1,1)
ceniza <- data.frame(Fuente=rep(unos,8, each=1), Prop = rep(unos,4,each=2),
Temp = rep(unos,2,each=4), mPa = c(43.4,50.1,27.3,47,
39.4,43.8,25,40.7,42.8,49.5,29.1,48.4,38.4,45,24.2,42.9))
ceniza2 <- data.frame(Fuente= gl(2,1,16,c("-","+")), Prop = gl(2,2,16,c("-","+")),
Temp = gl(2,4,16,c("-","+")), mPa = c(43.4,50.1,27.3,47,
39.4,43.8,25,40.7,42.8,49.5,29.1,48.4,38.4,45,24.2,42.9))
ceniza.fit2<-aov(mPa~Fuente*Prop*Temp,data=ceniza2)
summary(ceniza.fit2)## Df Sum Sq Mean Sq F value Pr(>F)
## Fuente 1 597.8 597.8 691.101 4.71e-09 ***
## Prop 1 287.3 287.3 332.142 8.45e-08 ***
## Temp 1 91.2 91.2 105.436 6.96e-06 ***
## Fuente:Prop 1 150.1 150.1 173.483 1.05e-06 ***
## Fuente:Temp 1 3.1 3.1 3.540 0.0967 .
## Prop:Temp 1 0.0 0.0 0.003 0.9584
## Fuente:Prop:Temp 1 0.3 0.3 0.350 0.5706
## Residuals 8 6.9 0.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Efectos:
ceniza.fit<-lm(mPa~Fuente*Prop*Temp,data=ceniza)
sum(ceniza$mPa*ceniza$Fuente)/8## [1] 12.225sum(ceniza$mPa*ceniza$Prop)/8## [1] -8.475sum(ceniza$mPa*ceniza$Prop*ceniza$Fuente)/8## [1] 6.125summary(ceniza.fit)##
## Call:
## lm(formula = mPa ~ Fuente * Prop * Temp, data = ceniza)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.100 -0.525 0.000 0.525 1.100
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.8125 0.2325 171.227 1.51e-15 ***
## Fuente 6.1125 0.2325 26.289 4.71e-09 ***
## Prop -4.2375 0.2325 -18.225 8.45e-08 ***
## Temp -2.3875 0.2325 -10.268 6.96e-06 ***
## Fuente:Prop 3.0625 0.2325 13.171 1.05e-06 ***
## Fuente:Temp -0.4375 0.2325 -1.882 0.0967 .
## Prop:Temp 0.0125 0.2325 0.054 0.9584
## Fuente:Prop:Temp -0.1375 0.2325 -0.591 0.5706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9301 on 8 degrees of freedom
## Multiple R-squared: 0.9939, Adjusted R-squared: 0.9886
## F-statistic: 186.6 on 7 and 8 DF, p-value: 3.184e-08Interacciones
plot(interactionMeans(ceniza.fit2))
with(ceniza, interaction.plot(Fuente, Prop, mPa)) # Std. Error 0.9301/sqrt(2)
Procedimiento para el análisis de un diseño factorial
- Estimar los efectos de los factores
- Formular el modelo. Con replicación se puede usar el modelo completo. Sin replicación se deben seleccionar un modelo reducido
- Prueba Estadística (ANOVA)
- Refinar el modelo
- Analizar los residuales
- Interpretar los resultados
Diseño \(2^4\): Hidrogenación catalítica
(Godawa,1984; Carlson R. and Carlson J.E., 2005). The experiments were run to determine how four experimental variables influence the yield of tetrahydrofuran in catalytic hydrogenation of furan over a palladium catalyst. The variables and the experimental domain are \(x_1:\) amount of catalyst/sustrate (0.7, 1.0 g/mol), \(x_2:\) hydrogen pressure (45, 55 bar), \(x_3:\) reaction temperature (75, 100 \(^{\circ}C\)), and \(x_4:\) stirring rate (340, 475 rpm). The experimental design and yields are:
hidro<- data.frame(x1=rep(unos,8,each=1),x2=rep(unos,4,each=2),x3=rep(unos,2,each=4),
x4=rep(unos,each=8), y = c(77.5,83.8,87.8,92.9,77.8,83.3,90.0,
94.1,94.1,98.3,94.2,98.3,97,99.3,98,100))# Análisis exploratorio
boxplot(y~x1, data=hidro)
boxplot(y~x4, data=hidro)
# No hay replicaciones, luego no se pueden incluir todas las interacciones
hidro.fit <- lm(y~x1*x2*x3*x4,data=hidro)
summary(hidro.fit)##
## Call:
## lm(formula = y ~ x1 * x2 * x3 * x4, data = hidro)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.6500 NaN NaN NaN
## x1 2.1000 NaN NaN NaN
## x2 2.7625 NaN NaN NaN
## x3 0.7875 NaN NaN NaN
## x4 5.7500 NaN NaN NaN
## x1:x2 -0.1875 NaN NaN NaN
## x1:x3 -0.3625 NaN NaN NaN
## x2:x3 0.3250 NaN NaN NaN
## x1:x4 -0.5250 NaN NaN NaN
## x2:x4 -2.5375 NaN NaN NaN
## x3:x4 0.3875 NaN NaN NaN
## x1:x2:x3 -0.0250 NaN NaN NaN
## x1:x2:x4 0.1375 NaN NaN NaN
## x1:x3:x4 -0.1375 NaN NaN NaN
## x2:x3:x4 -0.1250 NaN NaN NaN
## x1:x2:x3:x4 0.0000 NaN NaN NaN
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 15 and 0 DF, p-value: NA# Normal QQ-plot
plot(qnorm(ppoints(15)), sort(hidro.fit$effects)[-1], xlab="Cuantiles teoréticos",
pch=16,ylim = c(-11,25), ylab = "")
qqline(hidro.fit$effects[-1])
text(qnorm(ppoints(15)),sort(hidro.fit$effects[-1])+2,names(sort(hidro.fit$effects[-1])),cex=.7)
## Lenth Plot
library(BsMD)
LenthPlot(hidro.fit) # 
## alpha PSE ME SME
## 0.050000 0.562500 1.445952 2.935491# No se incluyen interacciones de tercer y cuarto orden
hidro.fit2 <- lm(y~x1+x2+x3+x4+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=hidro)
summary(hidro.fit2) # Efectos/2##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x1:x2 + x1:x3 + x1:x4 +
## x2:x3 + x2:x4 + x3:x4, data = hidro)
##
## Residuals:
## 1 2 3 4 5 6 7 8 9 10 11
## 0.150 0.100 0.125 -0.375 -0.425 0.175 0.150 0.100 -0.100 -0.150 -0.175
## 12 13 14 15 16
## 0.425 0.375 -0.125 -0.100 -0.150
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.6500 0.1040 881.393 3.57e-14 ***
## x1 2.1000 0.1040 20.196 5.50e-06 ***
## x2 2.7625 0.1040 26.567 1.41e-06 ***
## x3 0.7875 0.1040 7.573 0.000637 ***
## x4 5.7500 0.1040 55.297 3.66e-08 ***
## x1:x2 -0.1875 0.1040 -1.803 0.131220
## x1:x3 -0.3625 0.1040 -3.486 0.017543 *
## x1:x4 -0.5250 0.1040 -5.049 0.003937 **
## x2:x3 0.3250 0.1040 3.126 0.026089 *
## x2:x4 -2.5375 0.1040 -24.403 2.15e-06 ***
## x3:x4 0.3875 0.1040 3.727 0.013619 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4159 on 5 degrees of freedom
## Multiple R-squared: 0.999, Adjusted R-squared: 0.9969
## F-statistic: 488.9 on 10 and 5 DF, p-value: 7.805e-07qt(.975,5)*0.4159/sqrt(16) # Intervalo de confianza## [1] 0.2672762anova(hidro.fit2) ## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 70.56 70.56 407.8613 5.504e-06 ***
## x2 1 122.10 122.10 705.7948 1.413e-06 ***
## x3 1 9.92 9.92 57.3555 0.0006368 ***
## x4 1 529.00 529.00 3057.8035 3.658e-08 ***
## x1:x2 1 0.56 0.56 3.2514 0.1312199
## x1:x3 1 2.10 2.10 12.1532 0.0175428 *
## x1:x4 1 4.41 4.41 25.4913 0.0039370 **
## x2:x3 1 1.69 1.69 9.7688 0.0260895 *
## x2:x4 1 103.02 103.02 595.5058 2.154e-06 ***
## x3:x4 1 2.40 2.40 13.8873 0.0136193 *
## Residuals 5 0.86 0.17
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1with(hidro, interaction.plot(x2, x4, y)) # Gráfico de interacción
Superficie de respuesta
hidro2<- data.frame(x1=rep(c(.7,1),8,each=1),x2=rep(c(45,55),4,each=2),
x3=rep(c(75,100),2,each=4), x4=rep(c(340,475),each=8),
y = c(77.5,83.8,87.8,92.9,77.8,83.3,90.0,
94.1,94.1,98.3,94.2,98.3,97,99.3,98,100))
hidro2.fit <- lm(y~x1+x2+x3+x4+x1:x2+x1:x3+x1:x4+x2:x3+x2:x4+x3:x4,data=hidro2)
tmp <- list(x1=seq(.7,1,length=15),x2=seq(45,55,length=15),
x3 = seq(75,100,length=15), x4 = seq(340,475,length=15))
new.data <- expand.grid(tmp)
new.data$fit <- predict(hidro2.fit,new.data)
contourplot(fit~x2+x4,new.data,xlab="x2", ylab="x4",cuts=20)
contourplot(fit~x1+x2,new.data,xlab="x1", ylab="x2",cuts=20)
Número de replicas
Sea \(n_F\) el número total de experimentos en el diseño factorial \(2^k\), es decir, \(n_F = r\cdot 2^k\). Para el diseño factorial \(2^k\) se puede mostrar que el número neceario de corridas del experimento para obtener una potencia de 0.95, cuando el nivel de significancia \(\alpha=0.05\), es
\[n_F^* = \bigg(\frac{8 \sigma}{\Delta}\bigg)^2,\]
donde \(\sigma\) es la desviación estándar del error experimental y \(\Delta\) es el tamaño del efecto práctico.
Si el número de puntos factoriales, \(2^k\), es menor que el número de corridas necesarias para la precisión deseada, entonces todo el diseño factorial debe ser replicado \(r\) veces tal que \(r\cdot 2^k \geq n_F^*\).
Curvatura
Note que en situaciones cuando no hay replicaciones en todos los puntos de diseño, la única manera en que se puede estimar el error es asumiendo que las interacciones de orden alto son cero. Sin embargo, si las interacciones no son cero, el error estimado estará inflado. Luego el estimado de \(s_E\) es también inflado y el valor de la estadística \(t\) será pequeño. Y por lo tanto se puede dar que concluyamos que algunos efectos principales o interacciones no son significativos, cuando en en realidad sí lo son.
El análisis de datos provenientes de un diseño factorial asume que el modelo lineal describe adecuadamente los cambios en la respuesta dado los cambios en los factores. Sin embargo, con solo dos niveles en cada factor, no es posible chequear la validez de este supuesto. Además, si existen no linealidad, el diseño factorial producirá predicciones confiables solamente cerca a los puntos de diseño. En otras regiones, las predicciones no serán confiables. Esto es un gran problema, pues lo se quiere es hacer predicciones confiables sobre toda la región de experimentación.
Experimentos que incluyen replicas en el punto central son bastante comunes en la practica ya que permiten estimar el error y chequear que tan adecuado es el modelo lineal. Ya que al tener puntos centrales, se obtienen replicaciones que pueden ser usadas para hacer una estimación pura del error. Además estos puntos permiten chequear si hay curvatura comparando la respuesta observada en los puntos y el valor predecido en las condiciones centrales.
Ejemplo
En un proceso químico un reactivo se mantiene en un recipiente para enfriarlo. Cuando alcanza la temperatura deseada, un segundo reactivo es rápidamente cargado en una mezcla acuasa/orgánica. Se llevaron a cabo experimentos para optimizar la producción. Se estudiarion dos factores: el primero fue la temperatura en el recipiente y el segundo fue la razón en la mezcla acuosa/orgánica. La respuesta fue el porcentaje de producción. Se realizaron dos replicas en las cuatro condiciones experimentales. Finalmente, dado que los factores son cuantitativos, se hicieron cuatro replicas más en el punto central.
PQ <- data.frame(Temp = factor(c(gl(2,2,8),rep("0",4))),
Razon = factor(c(gl(2,4,8),rep("0",4))),
Prod=c(68.5,68.3, 72,72, 73,70,72.4, 73,74,75,74,73.5))Se puede usar \(C=\bar{Y}_{\mbox{centro}}-\bar{Y}_{\mbox{factorial}}\) para probar si existe alguna curvatura. El estadístio de prueba es
\[ t_C = \frac{C}{S_C},\]
donde
\[S_C = \sqrt{s_P^2\left(\frac{1}{n_c}+\frac{1}{n_F}\right)}.\]
\(n_C\) y \(n_F\) son el número de observaciones en el punto central y en los puntos factoriales, respectivamente. \(S_P\) es el promedio ponderado de los cinco estimados de la varianza, donde los pesos son el número de observaciones menos uno.
library(tidyverse)
mediasVar <- PQ %>% group_by(Temp,Razon) %>% summarize(mu = mean(Prod),v = var(Prod))## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by Temp and Razon.
## ℹ Output is grouped by Temp.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(Temp, Razon))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.sp = sum(mediasVar$v*c(3,1,1,1,1))/7
Curv = mediasVar$mu[1] - mean(mediasVar$mu[-1])
sc = sqrt(sp*(1/4+1/8))
tc = Curv/sc
2*(1-pt(tc,7))## [1] 0.001126661Cuando la prueba de curvatura es significativa, como en este ejemplo, las predicciones lineales no son válidas dentro de la región experimental. Para solucionar este problema, es necesario hacer más experimentos y ajustar un modelo cuadrático.
Diseños fraccionados
Como se mencionó anteriormente los diseños factoriales son eficientes, pero tienen el problema de que cuando el número de factores es grande, el mínimo número de experimentos es muy grande. Por esta razón, cuando se busca seleccionar, de un grupo de \(k\) factores, solamente los factores importantes, se puede limitar el número de tratamientos que se exploran. Es decir, en lugar de estudiar todos las \(2^k\) tratamientos, se estudia una fracción de ellos. Por ejemplo, si se tienen 6 factores, un \(\frac12\) diseño fraccionado consiste de \(2^{6-1}=32\) corridas.
La selección de que diseños eliminar no se puede hacer al azar, debido a que se pueden crear problemas de confusión. Es decir, no se puede distinguir a que factor corresponde el efecto.
El método de diseños fraccionados inicia con un diseño factorial completo con el número de corridas que se quiere usar y luego se le adicionan otros factores.
Diseños fraccionados \(2^{k-1}\)
- Escriba el diseño factorial \(2^{k-1}\)
- Adicione el \(k\)-ésimo factor al diseño asignando usando los niveles de la columna con mayor orden de interacción
- Use las \(2^{k-1}\) columnas más la columna con el \(k\)-ésimo factor para definir el diseño
unos<- c(-1,1)
d41 <- data.frame(x1=rep(unos,4),x2=rep(unos,2,each=2),x3=rep(unos,each=4))
d41$x12 <- d41$x1*d41$x2
d41$x13 <- d41$x1*d41$x3
d41$x23 <- d41$x2*d41$x3
d41$x4 <- d41$x1*d41$x2*d41$x3
d41[,c("x1","x2","x3","x4")]## x1 x2 x3 x4
## 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 1Patrón de confusión
El generador del diseño que se definió anteriormente es \[4 = 123\] Usando este generador se puede encontrar el patrón de confusión: \[4^2=123\cdot 4 \Rightarrow I = 1234. \qquad 24\cdot I = 24\cdot 1234 \Rightarrow 24 = 12^234^2 = 13\]
Ejemplo – Pressure Vessel (Montgomery, 2013)
A chemical product is produced in a pressure vessel. A factorial experiment is carried out in the pilot plant to study the factors thought to influence the filtration rate of this product. The four factors are temperature (\(A\)), pressure (\(B\)), concentration of formaldehyde (\(C\)), and stirring rate (\(D\)). Each factor is present at two levels. The process engineer is interested in maximizing the filtration rate. Current process conditions give filtration rates of around 75 gal/h. The process also currently uses the concentration of formaldehyde, factor C, at the high level. The engineer would like to reduce the formaldehyde concentration as much as possible but has been unable to do so because it always results in lower filtration rates.
unos <- c(-1,1)
filtracion<- data.frame(x1=rep(unos,8,each=1),x2=rep(unos,4,each=2),
x3=rep(unos,2,each=4),x4=rep(unos,each=8),
y = c(45,71,48,65,68,60,80,65,43,100,45,
104,75,86,70,96))
# efectos
sum(filtracion$y*filtracion$x1)/8## [1] 21.625sum(filtracion$y*filtracion$x2)/8## [1] 3.125filtracion.fit <- lm(y~x1*x2*x3*x4,data=filtracion)
summary(filtracion.fit)##
## Call:
## lm(formula = y ~ x1 * x2 * x3 * x4, data = filtracion)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.0625 NaN NaN NaN
## x1 10.8125 NaN NaN NaN
## x2 1.5625 NaN NaN NaN
## x3 4.9375 NaN NaN NaN
## x4 7.3125 NaN NaN NaN
## x1:x2 0.0625 NaN NaN NaN
## x1:x3 -9.0625 NaN NaN NaN
## x2:x3 1.1875 NaN NaN NaN
## x1:x4 8.3125 NaN NaN NaN
## x2:x4 -0.1875 NaN NaN NaN
## x3:x4 -0.5625 NaN NaN NaN
## x1:x2:x3 0.9375 NaN NaN NaN
## x1:x2:x4 2.0625 NaN NaN NaN
## x1:x3:x4 -0.8125 NaN NaN NaN
## x2:x3:x4 -1.3125 NaN NaN NaN
## x1:x2:x3:x4 0.6875 NaN NaN NaN
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 15 and 0 DF, p-value: NAanova(filtracion.fit)## Warning in anova.lm(filtracion.fit): ANOVA F-tests on an essentially perfect
## fit are unreliable## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 1870.56 1870.56 NaN NaN
## x2 1 39.06 39.06 NaN NaN
## x3 1 390.06 390.06 NaN NaN
## x4 1 855.56 855.56 NaN NaN
## x1:x2 1 0.06 0.06 NaN NaN
## x1:x3 1 1314.06 1314.06 NaN NaN
## x2:x3 1 22.56 22.56 NaN NaN
## x1:x4 1 1105.56 1105.56 NaN NaN
## x2:x4 1 0.56 0.56 NaN NaN
## x3:x4 1 5.06 5.06 NaN NaN
## x1:x2:x3 1 14.06 14.06 NaN NaN
## x1:x2:x4 1 68.06 68.06 NaN NaN
## x1:x3:x4 1 10.56 10.56 NaN NaN
## x2:x3:x4 1 27.56 27.56 NaN NaN
## x1:x2:x3:x4 1 7.56 7.56 NaN NaN
## Residuals 0 0.00 NaNplot(qnorm(ppoints(15)), sort(filtracion.fit$effects)[-1], xlab="Cuantiles teoréticos",
pch=16, ylab = "")
qqline(filtracion.fit$effects[-1])
text(qnorm(ppoints(15)),sort(filtracion.fit$effects[-1])+2,names(sort(filtracion.fit$effects[-1])),cex=.7)
LenthPlot(filtracion.fit)
## alpha PSE ME SME
## 0.050000 2.625000 6.747777 13.698960with(filtracion, interaction.plot(x1, x3, y))
with(filtracion, interaction.plot(x1, x4, y))
Let’s see what have happened if a half-fraction of the \(2^4\) design had been run instead of the full factorial.
filtro <- filtracion[c(1,10,11,4,13,6,7,16),]
sum(filtro$y*filtro$x1)/4## [1] 19sum(filtro$y*filtro$x4)/4## [1] 16.5filtro.fit0 <- lm(y~x1*x2*x3*x4,data=filtro)
summary(filtro.fit0)##
## Call:
## lm(formula = y ~ x1 * x2 * x3 * x4, data = filtro)
##
## Residuals:
## ALL 8 residuals are 0: no residual degrees of freedom!
##
## Coefficients: (8 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.75 NaN NaN NaN
## x1 9.50 NaN NaN NaN
## x2 0.75 NaN NaN NaN
## x3 7.00 NaN NaN NaN
## x4 8.25 NaN NaN NaN
## x1:x2 -0.50 NaN NaN NaN
## x1:x3 -9.25 NaN NaN NaN
## x2:x3 9.50 NaN NaN NaN
## x1:x4 NA NA NA NA
## x2:x4 NA NA NA NA
## x3:x4 NA NA NA NA
## x1:x2:x3 NA NA NA NA
## x1:x2:x4 NA NA NA NA
## x1:x3:x4 NA NA NA NA
## x2:x3:x4 NA NA NA NA
## x1:x2:x3:x4 NA NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 7 and 0 DF, p-value: NAfiltro.fit1 <- lm(y~x1*(x2+x3+x4),data=filtro)
summary(filtro.fit1)##
## Call:
## lm(formula = y ~ x1 * (x2 + x3 + x4), data = filtro)
##
## Residuals:
## ALL 8 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.75 NaN NaN NaN
## x1 9.50 NaN NaN NaN
## x2 0.75 NaN NaN NaN
## x3 7.00 NaN NaN NaN
## x4 8.25 NaN NaN NaN
## x1:x2 -0.50 NaN NaN NaN
## x1:x3 -9.25 NaN NaN NaN
## x1:x4 9.50 NaN NaN NaN
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 7 and 0 DF, p-value: NAfiltro.fit2 <- lm(y~x1*(x3+x4),data=filtro)
summary(filtro.fit2)##
## Call:
## lm(formula = y ~ x1 * (x3 + x4), data = filtro)
##
## Residuals:
## 1 10 11 4 13 6 7 16
## -1.25 -0.25 1.25 0.25 -1.25 -0.25 1.25 0.25
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.7500 0.6374 111.00 8.11e-05 ***
## x1 9.5000 0.6374 14.90 0.00447 **
## x3 7.0000 0.6374 10.98 0.00819 **
## x4 8.2500 0.6374 12.94 0.00592 **
## x1:x3 -9.2500 0.6374 -14.51 0.00471 **
## x1:x4 9.5000 0.6374 14.90 0.00447 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.803 on 2 degrees of freedom
## Multiple R-squared: 0.9979, Adjusted R-squared: 0.9926
## F-statistic: 188.6 on 5 and 2 DF, p-value: 0.005282Diseños fraccionados \(2^{k-2}\)
unos<- c(-1,1)
d52 <- data.frame(x1=rep(unos,4),x2=rep(unos,2,each=2),x3=rep(unos,each=4))
d52$x4 <- d52$x1*d52$x2
d52$x5 <- d52$x1*d52$x3
d52$x23 <- d52$x2*d52$x3
d52$x123 <- d52$x1*d52$x2*d52$x3
d52[,c("x1","x2","x3","x4","x5")]## x1 x2 x3 x4 x5
## 1 -1 -1 -1 1 1
## 2 1 -1 -1 -1 -1
## 3 -1 1 -1 -1 1
## 4 1 1 -1 1 -1
## 5 -1 -1 1 1 -1
## 6 1 -1 1 -1 1
## 7 -1 1 1 -1 -1
## 8 1 1 1 1 1\[4 = 12 \qquad 5 = 13 \Rightarrow I = 124 = 135 = 2345\]
Ejemplo - Impurezas (\(2^{7-3}\))
unos<- c(-1,1)
d74 <- data.frame(x1=rep(unos,8),x2=rep(unos,4,each=2),x3=rep(unos,2,each=4),
x4=rep(unos,each=8))
d74$x5 <- d74$x1*d74$x2*d74$x3
d74$x6 <- d74$x2*d74$x3*d74$x4
d74$x7 <- d74$x1*d74$x3*d74$x4
d74[,c("x1","x2","x3","x4","x5","x6","x7")]## x1 x2 x3 x4 x5 x6 x7
## 1 -1 -1 -1 -1 -1 -1 -1
## 2 1 -1 -1 -1 1 -1 1
## 3 -1 1 -1 -1 1 1 -1
## 4 1 1 -1 -1 -1 1 1
## 5 -1 -1 1 -1 1 1 1
## 6 1 -1 1 -1 -1 1 -1
## 7 -1 1 1 -1 -1 -1 1
## 8 1 1 1 -1 1 -1 -1
## 9 -1 -1 -1 1 -1 1 1
## 10 1 -1 -1 1 1 1 -1
## 11 -1 1 -1 1 1 -1 1
## 12 1 1 -1 1 -1 -1 -1
## 13 -1 -1 1 1 1 -1 -1
## 14 1 -1 1 1 -1 -1 1
## 15 -1 1 1 1 -1 1 -1
## 16 1 1 1 1 1 1 1d74$y <- c(2.6,.89,1.82,.49,1.83,1.26,2.23,1.98,0.67,.07,.51,.72,1.61,.08,.53,.16)
d74.fit <- lm(y~x1*(x2+x3+x4+x5+x6+x7),d74)