Diseños Factoriales
Jaime Yahir Garrido Del Valle
21 de mayo de 2022
Diseños \(2^{k}\)
Para el primer desarrollo se trabaja con un diseño completo \(2^{k}\) donde \(k = 2\). La base de datos será la propuesta por Gutierrez(1998). Para este ejemplo, interesa conocer los efectos en la respuesta media de la vibración de una ranura. Tenemos dos factores \(x_{1}\) el cual es tamaño de la broca y \(x_{2}\) la velocidad. Para \(x_{1}: -1 = 1/16 \text{pulgadas}, +1 = 1/8 \text{pulgadas}\). Para el fector 2 \(x_{2}: -1 = 40, +1 = 90\).
library(dplyr)
y = c( 18.2, 27.2, 15.9, 41.0, 18.9, 24.0, 14.5, 43.9, 12.9, 22.4, 15.1,
36.3, 14.4, 22.5, 14.2, 39.9)
x1 = c(-1, +1, -1, +1, -1, +1, -1, +1, -1, +1, -1, +1, -1, +1, -1, +1 )
x2 = c(-1, -1, +1, +1, -1, -1, +1, +1, -1, -1, +1, +1, -1, -1, +1, +1 )
datos = data.frame(y, x1, x2)
head(datos) y x1 x2
1 18.2 -1 -1
2 27.2 1 -1
3 15.9 -1 1
4 41.0 1 1
5 18.9 -1 -1
6 24.0 1 -1Para empezar se realiza una inspección gráfica con gráficos de cajas y alambres:
par(mfrow = c(1,2))
boxplot(y ~ x1, main = "Vibraciones por tamaño")
boxplot(y ~ x2, main = "Vibraciones por velocidad")
par(mfrow = c(1,1))Como se observa, parece que existen diferencias entre los dos niveles de velocidad de de la broca y el tamaño de la broca para explicar las vibraciones. Para conocer si existe evidencia de efectos en las interacciones, se calcula el mapa de contorno de la regresión lineal buscando curvaturas sobre las curvas de nivel del plano \(x_{1}, x_{2}\):
library(rsm)
modelo1 <- rsm(y ~ FO(x1, x2, x1*x2), data = datos)
par(mfrow = c(1,2))
persp(modelo1, x2 ~ x1,
zlab = "y",
contours = list(z = "bottom", col = "colors"), # posicion y color
at = c(summary(modelo1$canonical$xs)),
theta = 15, # coordenadas graficas
phi = 20)
# Grafico de contornos
contour(modelo1, ~ x1 + x2, image = TRUE)
Como se observa en la curva de nivel del hiperplano de la regresión muestra cierto nivel de curvatura por lo que quizá se detecten diferencias significativas en el ANOVA. Se calcula ahora el análisis de varianza:
anova(lm(y ~ x1 + x2 + x1*x2))## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 1107.23 1107.23 185.252 1.175e-08 ***
## x2 1 227.26 227.26 38.023 4.826e-05 ***
## x1:x2 1 303.63 303.63 50.801 1.201e-05 ***
## Residuals 12 71.72 5.98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(x.factor = x1, trace.factor = x2, response = y)
orden = seq(1,16,1)
modelo1 = lm(y ~ x1 + x2 + x1*x2)
par(mfrow=c(2,2))
residuos = modelo1$residuals
estimado = modelo1$fitted.values
plot(y = residuos, x = estimado, main = "A")
plot(x1, residuos, main = "B")
qqnorm(residuos, main = "C")
qqline(residuos)
plot(residuos, orden, main = "D")
El siguiente ejemplo procede de Monstgomery(), en el se presenta un estudio sobre el efecto del porcentaje de carbonatación (A), la preparación de la operación (B) yla velocidad de linea sobre la altura de llenado de una bebida carbonatada. La variable de respuesta (y) es la desviación de la altura de llenado. Este es un ejemplo aplicado al control de calidad. Los datos son los siguientes:
A = c(-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1)
B = c(-1,-1,+1,+1,-1,-1,+1,+1,-1,-1,+1,+1,-1,-1,+1,+1)
C = c(-1,-1,-1,-1,+1,+1,+1,+1,-1,-1,-1,-1,+1,+1,+1,+1)
y = c(-3,0,-1,2,-1,2,1,6, -1,1,0,3,0,1,1,5)
datos = data.frame(A,B,C,y)
head(datos)## A B C y
## 1 -1 -1 -1 -3
## 2 1 -1 -1 0
## 3 -1 1 -1 -1
## 4 1 1 -1 2
## 5 -1 -1 1 -1
## 6 1 -1 1 2Realizamos una inspección con un gráfico de cajas y alambres:
par(mfrow = c(1,3))
boxplot(y ~ A)
boxplot(y ~ B)
boxplot(y ~ C)
Realizamos el ANOVA para analizar tanto si existen diferencias significativas entre factores como para analizar la significancia de las interacciones:
modelo1 = lm(y ~ A*B*C, data = datos)
anova(modelo1)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 36.00 36.000 57.6 6.368e-05 ***
## B 1 20.25 20.250 32.4 0.0004585 ***
## C 1 12.25 12.250 19.6 0.0022053 **
## A:B 1 2.25 2.250 3.6 0.0943498 .
## A:C 1 0.25 0.250 0.4 0.5447373
## B:C 1 1.00 1.000 1.6 0.2415040
## A:B:C 1 1.00 1.000 1.6 0.2415040
## Residuals 8 5.00 0.625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Siguiendo a Montgomery() solo los efectos principales son significativos y la interacción de A y B pero solo con una significacia estadística del 90%. El modelo de regresión propuesto y el gráfico de contorno es el siguiente:
modelo2 = lm(y ~ A + B + C + A*B, data = datos)
summary(modelo2)##
## Call:
## lm(formula = y ~ A + B + C + A * B, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1250 -0.4375 -0.1250 0.4375 1.1250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.000 0.203 4.927 0.000452 ***
## A 1.500 0.203 7.391 1.38e-05 ***
## B 1.125 0.203 5.543 0.000175 ***
## C 0.875 0.203 4.311 0.001233 **
## A:B 0.375 0.203 1.848 0.091700 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8118 on 11 degrees of freedom
## Multiple R-squared: 0.9071, Adjusted R-squared: 0.8733
## F-statistic: 26.84 on 4 and 11 DF, p-value: 1.267e-05library(rsm)
modelo2 <- rsm(y ~ FO(A, B, A*B), data = datos)
par(mfrow = c(1,2))
persp(modelo2, A ~ B,
zlab = "y",
contours = list(z = "bottom", col = "colors"), # posicion y color
#at = c(summary(modelo1$canonical$xs)),
theta = 15, # coordenadas graficas
phi = 20)
contour(modelo2, ~ A + B, image = TRUE)
Se calcula ahora el gráfico de interacción entre variables A y B:
interaction.plot(x.factor = A, trace.factor = B, response = y)
Finalmente comprobamos los supuestos
orden = seq(1,16,1)
modelo1 = lm(y ~ A + B + C + A*B)
par(mfrow=c(2,3))
residuos = modelo1$residuals
estimado = modelo1$fitted.values
plot(y = residuos, x = estimado)
qqnorm(residuos)
qqline(residuos)
plot(residuos, orden)
plot(A, residuos)
plot(B, residuos)
plot(C, residuos)
Diseños Fraccionados
Para este caso seguiremos a Montgomery(). Se investigaron cinco factores de un proceso de manufactura de un circuito integrado en un diseño \(2^{5-1}\). Los factores fueron: A = Ajuste de apertura (pequeña y grande), B = tiempo de exposición (20% abajo del nominal, 20% arriba del nominal), C = tiempo de desarrollo (30s, 40s), D = Tamaño de la máscara (pequeña, grande) y E = Tiempo de grabado (14.5min, 15.5min) y el y = rendimiento.
A = c(-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1,-1,+1)
B = c(-1,-1,+1,+1,-1,-1,+1,+1,-1,-1,+1,+1,-1,-1,+1,+1)
C = c(-1,-1,-1,-1,+1,+1,+1,+1,-1,-1,-1,-1,+1,+1,+1,+1)
D = c(-1,-1,-1,-1,-1,-1,-1,-1,+1,+1,+1,+1,+1,+1,+1,+1)
E = c(+1,-1,-1,+1,-1,+1,+1,-1,-1,+1,+1,-1,+1,-1,-1,+1)
y = c(8,9,34,52,16,22,45,60,6,10,30,50,15,21,44,63)
datos = data.frame(A,B,C,D,E,y)
head(datos)## A B C D E y
## 1 -1 -1 -1 -1 1 8
## 2 1 -1 -1 -1 -1 9
## 3 -1 1 -1 -1 -1 34
## 4 1 1 -1 -1 1 52
## 5 -1 -1 1 -1 -1 16
## 6 1 -1 1 -1 1 22Realizamos la revisión de los dos primeros efectos del diseño experimental:
modelo3 = lm(y ~ (.)^2, data = datos)
summary(modelo3)##
## Call:
## lm(formula = y ~ (.)^2, data = datos)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.3125 NA NA NA
## A 5.5625 NA NA NA
## B 16.9375 NA NA NA
## C 5.4375 NA NA NA
## D -0.4375 NA NA NA
## E 0.3125 NA NA NA
## A:B 3.4375 NA NA NA
## A:C 0.1875 NA NA NA
## A:D 0.5625 NA NA NA
## A:E 0.5625 NA NA NA
## B:C 0.3125 NA NA NA
## B:D -0.0625 NA NA NA
## B:E -0.0625 NA NA NA
## C:D 0.4375 NA NA NA
## C:E 0.1875 NA NA NA
## D:E -0.6875 NA NA NA
##
## 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: NAComo se observa no es posible estimar la significacia estadística de los coeficientes de la regresión debido a que se requiere más información de la disponible. Para aplicar el diseño fraccionado se realiza una inspección gráfica.
Se realiza el gráfico de probabilidad normal o gráfico de Daniel. Este gráfico muestra los efectos estandarizados relacionada con una línea de ajuste de distribución para el caso en que todos los efectos son cero. Los efectos estandarizados son estadísticos t y prueban la hipótesis nula de que el efecto es igual a 0. Este gráfico ayuda a definir la posible significancia estadística de los efectos de un modelo factorial donde no es posible inspeccionar mediante regresión la signifcancia estadística de los coeficientes.
library(daewr)
LGB(coef(modelo3)[-1], rpt = FALSE)
Una vez seleccionados los efectos solo tomaremos la fracción de efectos de A, B, C y AB.
modelo4 = lm(y ~ A + B + C + A*B, data = datos)
anova(modelo4)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 495.1 495.1 193.195 2.535e-08 ***
## B 1 4590.1 4590.1 1791.244 1.560e-13 ***
## C 1 473.1 473.1 184.610 3.214e-08 ***
## A:B 1 189.1 189.1 73.781 3.302e-06 ***
## Residuals 11 28.2 2.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se muestra el gráfico de la interacción de los factores AB
interaction.plot(x.factor = A, trace.factor = B, response = y)
Finalmente revisamos los supuestos del modelo
orden = seq(1,16,1)
modelo1 = lm(y ~ A + B + C + A*B)
par(mfrow=c(2,3))
residuos = modelo1$residuals
estimado = modelo1$fitted.values
plot(y = residuos, x = estimado)
qqnorm(residuos)
qqline(residuos)
plot(residuos, orden)
plot(A, residuos)
plot(B, residuos)
plot(C, residuos)
El en siguiente ejemplo se utiliza la base de datos de Carreón (2022) donde se buscan los factores e interacciones importantes en un proceso químico. Para este caso se busca resolver un análisis fraccionario en un diseño \(2^{5-1}\).
datos <- read.csv("https://raw.githubusercontent.com/jpch26/Fractional-Designs-with-R/master/data/yield_data.csv")
head(datos)## FeedRate Catalyst AgitationSpeed Temperature Concentration Yield
## 1 -1 -1 -1 -1 1 56
## 2 1 -1 -1 -1 -1 53
## 3 -1 1 -1 -1 -1 63
## 4 1 1 -1 -1 1 65
## 5 -1 -1 1 -1 -1 53
## 6 1 -1 1 -1 1 55Si se realiza la regresión por mínimos cuadrados ya que solo se tienen 16 observaciones.
modelo4 <- lm(Yield ~ (.)^2, data = datos)
summary(modelo4)##
## Call:
## lm(formula = Yield ~ (.)^2, data = datos)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.525e+01 NA NA NA
## FeedRate -1.000e+00 NA NA NA
## Catalyst 1.025e+01 NA NA NA
## AgitationSpeed 1.544e-16 NA NA NA
## Temperature 6.125e+00 NA NA NA
## Concentration -3.125e+00 NA NA NA
## FeedRate:Catalyst 7.500e-01 NA NA NA
## FeedRate:AgitationSpeed 2.500e-01 NA NA NA
## FeedRate:Temperature -3.750e-01 NA NA NA
## FeedRate:Concentration 6.250e-01 NA NA NA
## Catalyst:AgitationSpeed 7.500e-01 NA NA NA
## Catalyst:Temperature 5.375e+00 NA NA NA
## Catalyst:Concentration 6.250e-01 NA NA NA
## AgitationSpeed:Temperature 1.250e-01 NA NA NA
## AgitationSpeed:Concentration 1.125e+00 NA NA NA
## Temperature:Concentration -4.750e+00 NA NA NA
##
## 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: NAPara realizar el diseño fraccionado se realiza el gráfico de Daniel, ahora con la librería DanielPlot del paquete FrF2
library("FrF2")## Loading required package: DoE.base## Loading required package: grid## Loading required package: conf.design##
## Attaching package: 'DoE.base'## The following objects are masked from 'package:stats':
##
## aov, lm## The following object is masked from 'package:graphics':
##
## plot.design## The following object is masked from 'package:base':
##
## lengthsDanielPlot(modelo4)
Como se observa el número de efectos principales que son significativos es reducido, se pueden ampliar nuestra visión realizando un gráfico de pareto:
library(ggplot2)
coeficientes <- coefficients(modelo4)[-1] # We discard the intercept
efectos <- data.frame(
Efecto = names(coeficientes),
Valor = unname(coeficientes),
ValorAbsoluto = abs(unname(coeficientes)),
Signo = as.character(unname(sign(coeficientes)))
)
ggplot(efectos, aes(ValorAbsoluto, reorder(Efecto, -ValorAbsoluto, abs))) +
geom_col(aes(fill = Signo)) +
xlab("Magnitud del Efecto") +
ylab("Efecto") +
theme_minimal()
Se toma solo los factores con los efectos más grandes para realizar el análisis de varianza.
modelo5 <- lm(Yield ~ Catalyst*Temperature*Concentration, data = datos)
anova(modelo5)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Catalyst 1 1681.00 1681.00 213.4603 4.725e-07 ***
## Temperature 1 600.25 600.25 76.2222 2.316e-05 ***
## Concentration 1 156.25 156.25 19.8413 0.002127 **
## Catalyst:Temperature 1 462.25 462.25 58.6984 5.953e-05 ***
## Catalyst:Concentration 1 6.25 6.25 0.7937 0.398998
## Temperature:Concentration 1 361.00 361.00 45.8413 0.000142 ***
## Catalyst:Temperature:Concentration 1 1.00 1.00 0.1270 0.730795
## Residuals 8 63.00 7.88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comos se puede observar de los efectos seleccionados solo Catalyst:Concentration no fue significativo, así como el efecto de tercer orden Catalyst:Temperature:Concentration.