Introducción

En el siguiente Notebook se encuentra el ejercicio 5.24 del libro Montgomery, D. (2017). Design and analysis of exp eriments. Ninth edition. Wiley. Chapter 5: Intro duction to Factorial Designs

Ejercicio 5.24

The quality control department of a fabric finishing plant is studying the effect of several factors on the dyeing of cotton–synthetic cloth used to manufacture men’s shirts. Three operators, three cycle times, and two temperatures were selected, and three small specimens of cloth were dyed under each set of conditions. The finished cloth was compared to a standard, and a numerical score was assigned. The results are as follows. Analyze the data and draw conclusions. Comment on the model’s adequacy.

# Operador 300 C
y1 <- c(23,24,25,36,35,36,28,24,27)
y2 <- c(27,28,26,34,38,39,35,35,34)
y3 <- c(31,32,29,33,34,35,26,27,25)
# Operador 350 C
x1 <- c(24,23,28,37,39,35,26,29,25)
x2 <- c(38,36,35,34,38,36,36,37,34)
x3 <- c(34,36,39,34,36,31,28,26,24)
# Ciclos de tiempos
tiempo <- c(40,40,40,50,50,50,60,60,60)
# Creación del dataframe
datos <- data.frame(tiempo,y1,y2,y3,x1,x2,x3)
# Presentación
DT::datatable(data = datos,
              colnames = c("Tiempos de Ciclos",
                           "(Operador 1) 300°C","(Operador 2) 300°C","(Operador 3) 300°C",
                           "(Operador 1) 350°C","(Operador 2) 350°C","(Operador 3) 350°C")
              )

Solución

  • Como es expresado el objetivo es determinar el efecto en el teñido de telas de algodón y sintéticas, siendo esta nuestra variable dependiente/respuesta.

  • Los factores a los cuales la variable de respuesta es sometida son, tres operadores a tres ciclos de tiempo distinto, con dos niveles de temperatura/tratamiento.

  • Donde tres muestras de telas fueron seleccionadas para cada condición, esto es, 3 replicas.

Sea

\[A = \text{Operador}\rightarrow\text{Tres Niveles (Operador 1, 2, 3)}\\ B = \text{Ciclos de Tiempo}\rightarrow\text{Tres Niveles (40, 50, 60)}\\ C = \text{Temperatura}\rightarrow\text{Dos Niveles (300°C y 350°C)}\]

# Mejorando la presentación de los datos para el análisis
efecto <- c(y1,y2,y3,x1,x2,x3)
operador <- factor(c(rep(1,9), rep(2,9), rep(3,9)))
tiempo <- factor(c(rep(40,3), rep(50,3), rep(60,3)))
temperatura <- factor(c(rep(300,27), rep(350,27)))
datos <- data.frame(efecto,operador,tiempo, temperatura)
# Presentación
DT::datatable(datos, caption = "Datos")
# Analísis Descriptitivo/Exploratorio
summarytools::dfSummary(datos,
                        plain.ascii = FALSE,
                        style = "grid",
                        graph.magnif = 0.75,
                        valid.col = FALSE,
                        tmp.img.dir  = "/tmp")

Data Frame Summary

datos

Dimensions: 54 x 4
Duplicates: 2

No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 efecto
[numeric]
Mean (sd) : 31.6 (5.1)
min < med < max:
23 < 34 < 39
IQR (CV) : 9 (0.2)
16 distinct values 0
(0.0%)
2 operador
[factor]
1. 1
2. 2
3. 3
18 (33.3%)
18 (33.3%)
18 (33.3%)
0
(0.0%)
3 tiempo
[factor]
1. 40
2. 50
3. 60
18 (33.3%)
18 (33.3%)
18 (33.3%)
0
(0.0%)
4 temperatura
[factor]
1. 300
2. 350
27 (50.0%)
27 (50.0%)
0
(0.0%)
  • Donde observamos en el histograma una posible bimodalidad, por la cantidad de valores en la cola izquierda.

  • Donde el efecto promedio del teñido de telas de algodón y sintéticas es de 31.6 con una desviación estándar de 5.6.

Calculo de la tabla ANOVA

Ajustando el modelo

\[y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_3+\beta_{12}x_1x_2+\beta_{13}x_1x_3+\beta_{23}x_2x_3+\beta_{123}x_1x_2x_3+\varepsilon\]

Donde contrastamos

# Función del Modelo
modelo1 <- lm(efecto ~ operador + tiempo + temperatura + operador*tiempo + operador*temperatura + tiempo*temperatura + operador*tiempo*temperatura, data = datos)

# Construcción tabla ANOVA
summary(aov(modelo1))
##                             Df Sum Sq Mean Sq F value   Pr(>F)    
## operador                     2  261.3  130.67  39.864 7.44e-10 ***
## tiempo                       2  436.0  218.00  66.508 8.14e-13 ***
## temperatura                  1   50.1   50.07  15.277 0.000393 ***
## operador:tiempo              4  355.7   88.92  27.127 1.98e-10 ***
## operador:temperatura         2   11.3    5.63   1.718 0.193895    
## tiempo:temperatura           2   78.8   39.41  12.023 0.000100 ***
## operador:tiempo:temperatura  4   46.2   11.55   3.523 0.015870 *  
## Residuals                   36  118.0    3.28                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Usando un \(\alpha = 0.01\), concluimos que:

  • El operador, el tiempo, y la temperatura son altamente significativos, la combinación de operador-tiempo y tiempo-temperatura, también son influyentes.

Así corriendo donde nuestro mayor interés son las interacciones

\[y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_3+\beta_{12}x_1x_2+\beta_{23}x_2x_3+\varepsilon\]

# Función de regresión
modelo2 <- lm(efecto ~ operador*tiempo + tiempo*temperatura, data = datos)
# Bondad de ajuste/Tabla ANOVA
summary(aov(modelo2))
##                    Df Sum Sq Mean Sq F value   Pr(>F)    
## operador            2  261.3  130.67  31.281 4.80e-09 ***
## tiempo              2  436.0  218.00  52.187 4.11e-12 ***
## temperatura         1   50.1   50.07  11.987 0.001244 ** 
## operador:tiempo     4  355.7   88.92  21.286 1.19e-09 ***
## tiempo:temperatura  2   78.8   39.41   9.434 0.000413 ***
## Residuals          42  175.4    4.18                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En donde ahora se tiene un modelo que todos son factores son significativos, lo que nos dice que podríamos medir el efecto de teñido de telas de algodón y sintéticas, sujetos a la interacción de los tres operadores y los tres niveles de tiempo y estos tres niveles de tiempo y los dos niveles de temperatura.

Gráficos de interes

Gráficos de Efectos

# Gráfico de efectos
plot.design(datos, xlab = "Factor", ylab = "Efecto")
grid(lty = 2)

  • Observamos como con el segundo operador, con el ciclo de tiempo de 50, y con la temperatura a 350 grados celsius se obtiene el efecto mas alto. (En base a la escala que se usa para medir el efecto en la tela.)
# Efectos de interes (interacciones)
par(mfrow = c(1, 2))
plot.design(efecto ~ operador*tiempo, data = datos,
            xlab = "Factor", ylab = "Efecto")
grid(lty = 2)
plot.design(efecto ~ tiempo*temperatura, data = datos,
            xlab = "Factor", ylab = "Efecto")
grid(lty = 2)

  • Dichas observaciones las podemos ver individualmente también al gráficar por separado los efectos con respecto a los operadores y al tiempo.

  • Otra observación interesante es como el máximo ciclo de tiempo, genera la respuesta mas baja junto con el primero operador.

Gráficos de Interacción

par(mfrow = c(1, 2))
# Interacción tiempo-operador
interaction.plot(x.factor = datos$tiempo,
                 trace.factor = datos$operador,
                 response = datos$efecto,
                 fun = mean,
                 xlab = "Tiempo",
                 ylab = "Media del Efecto",
                 col = c("pink","blue", "green"),
                 lty = 1,
                 lwd = 2,
                 trace.label = "Operador")
grid(lty = 2)
# Interacción tiempo-temperatura
interaction.plot(x.factor = datos$tiempo,
                 trace.factor = datos$temperatura,
                 response = datos$efecto,
                 fun = mean,
                 xlab = "Tiempo",
                 ylab = "Media del Efecto",
                 col = c("pink","blue"),
                 lty = 1,
                 lwd = 2,
                 trace.label = "Temperatura")
grid(lty = 2)

  • Donde primero observamos como el primer y el segundo operador con un ciclo de tiempo de 50, tienen efectos de respuesta similares, como el uno y el tres a un ciclo de tiempo de 60.

  • En base al criterio de decisión (si se busca una respuesta alta o baja) se puede empezar a buscar decisiones sobre el mejor operador y ciclo de tiempo.

  • Para el caso de la temperatura, en relación a los ciclos de tiempos observamos como, para ambos niveles de temperatura con ciclos de tiempo de 50 y 60 se obtienen promedios de respuesta similares a excepción del ciclo de tiempo de 40. En donde las respuestas varían.

  • De igual forma en base a lo que se busque, una respuesta alta o baja, se puede proponer una mejor combinación.

  • Asumiendo una respuesta alta, (un alto teñido de la tela), se podría recomendar usar a la metodología del segundo operador con una temperatura de 300 o 350 grados Celsius (la que sea mas barata, dado que el efecto parece decirnos que es similar), con un ciclo de tiempo de 50 (sean minutos/segundos o la unidad de tiempo que usen, que no fue dada en el enunciado).

Análisis de supuestos

Gráficos de diagnóstico

# Gráficos de diagnósticos
par(mfrow = c(2,2))
plot(aov(modelo2))

  • Donde observamos que no parece haber un patrón marcado con respecto a los residuos y los valores ajustados.

  • De igual forma con los valores estandarizados.

  • Aún así se encuentran algunos posibles valores atípicos.

  • El Q-Q Plot muestra ciertas desviaciones altas de la posible normalidad en relación a sus colas.

  • Residuos contra los niveles de factores no muestran un patrón marcado exceptuando a los valores atípicos que se observan.

Chequeo Analítico de Normalidad

patchwork::wrap_plots(ggplot2::ggplot(data = datos, ggplot2::aes(x = residuals(aov(modelo2)))) +
  ggplot2::geom_histogram(ggplot2::aes(y=..density..),
                          binwidth = 1) +
  ggplot2::geom_density() +
  ggplot2::ggtitle("Histograma de los residuos")+
  ggplot2::xlab("Residuos")+
  ggplot2::ylab("Frecuencia"),
  #Q-Q Plot
  ggpubr::ggqqplot(residuals(aov(modelo2)),
                 conf.int.level = 0.99,
                 ggtheme = ggplot2::theme_gray(),
                 title = "QQ-Plot",
                 ylab = "Valores Teoricos",
                 xlab = "Valores de Muestra"),
  ncol = 2
  )

\[H_0: \text{Normalidad}\\ vs\\ H_1: \text{No Normalidad}\]

shapiro.test(residuals(aov(modelo2)))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(modelo2))
## W = 0.98673, p-value = 0.8104

Donde con un \(\alpha=0.01\), y un p-value del 0.8104, no tenemos suficiente evidencia para rechazar \(H_0\) y decir que los residuos no son normales.

Igualdad de varianzas/Homocedasticidad

Efecto y Tiempo
bartlett.test(efecto ~ tiempo, data = datos)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  efecto by tiempo
## Bartlett's K-squared = 12.926, df = 2, p-value = 0.00156
Efecto y Operador
bartlett.test(efecto ~ operador, data = datos)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  efecto by operador
## Bartlett's K-squared = 2.5971, df = 2, p-value = 0.2729
Efecto y Temperatura
bartlett.test(efecto ~ temperatura, data = datos)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  efecto by temperatura
## Bartlett's K-squared = 0.1218, df = 1, p-value = 0.7271
Residuos
lmtest::bptest(aov(modelo2))
## 
##  studentized Breusch-Pagan test
## 
## data:  aov(modelo2)
## BP = 16.75, df = 11, p-value = 0.1155

Conclusión de la homocedasticidad

A excepción del efecto de respuesta y del tiempo no tenemos evidencia para rechazar el que haya homocedasticidad.

Comparaciones Multiples

Tiempo

plot(TukeyHSD(aov(modelo2), which="tiempo", conf.level = .99))
grid(lty = 2)

Operador

plot(TukeyHSD(aov(modelo2), which="operador", conf.level = .99))
grid(lty = 2)

Conclusión

  • Donde determinamos que en el caso de el tiempo, para la combinación de ciclos de 60-50 y 50-40 estas son significativamente distintas de 60-40

  • Para los operadores las combinaciones 3-2 y 2-1 son significativamente distintas de 3-1

Conclusión

  • Se encontró que la interacción de los operadores y el ciclo de tiempo del proceso y ademas, el ciclo de tiempo y las temperaturas, si afectan a la respuesta del teñido de la tela. Otra forma de verlo es como si esta de verdad depende de las interacciones de esta.

  • En particular los ciclos de tiempo de 50 con el segundo y tercer operador y cualquier temperatura generaban las respuestas de teñido mas altas. Con respecto a las comparaciones múltiples dependiendo del criterio de decisión del cliente se decidiría sobre cual combinación es mas optima.

  • Considerando que tan adecuado es el modelo, se puede decir que cumple con el supuesto de normalidad en sus residuos, y sus residuos son homocedasticos, exceptuando el problema de la homocedasticidad con respecto a la variable de respuesta y los ciclos de tiempo, lo cual es de considerarse con precaución, podría de manera conservadora decirse que es un modelo adecuado, dado que se detecta que efectivamente que el efecto del teñido es dependiente de los factores.