Diseño de experimentos Taguchi

El Diseño de Experimentos Taguchi es una metodología desarrollada por Genichi Taguchi que busca optimizar procesos y productos mediante el uso de matrices ortogonales, reduciendo la variabilidad y mejorando la robustez del diseño. Se enfoca en identificar los factores de control que minimizan el efecto del ruido (factores no controlables).

Tipos de factores

El Diseño de Experimentos Taguchi considera factores que pueden afectar a:

  • la media del proceso
  • la variación del proceso
  • la media y a la variación del proceso
  • o que no afectan a la media y a la variación del proceso.

Estos pueden ser factores que se pueden controlar dentro de los procesos y factores que no se pueden controlar, a los últimos se les da el nombre de factores de ruido. Vamos a suponer que tenemos un factor llamado Temperatura para un proceso de elaboración de ladrillos, del cual registramos sus resultados en dos niveles (\(600°C\) y \(800°C\)). La característica de calidad que se monitorea es la longitud de los ladrillos. A continuación simulamos resultados de las longitudes de los ladrillos, suponiendo que tienen una distribución normal y la longitud deseada es de \(25\) cm. En este caso se simula un caso en el que la temperatura afecta la media de la longitud de los ladrillos.

Temp_600 <- rnorm(50, mean = 25, sd = 0.2)
Temp_800 <- rnorm(50, mean = 30, sd = 0.2)

La media y desviación estándar del largo de los ladrillos a la temperatura de \(600°C\) es:

mean(Temp_600)
## [1] 25.00781
sd(Temp_600)
## [1] 0.2124092

mientras que la media y desviación estándar de la longitud de los ladrillos a la temperatura de \(800°C\) es:

mean(Temp_800)
## [1] 29.9802
sd(Temp_800)
## [1] 0.1946319

Veamos un diagrama de caja y bigotes para este caso:

# Cambiar el arreglo de los datos
Cat_1 <- rep("600_grados",50)
Cat_2 <- rep("800_grados",50)
Cat <- c(Cat_1,Cat_2)
Num <- c(Temp_600,Temp_800)
Boxplot_data <- data.frame(Temperatura = Cat, Valor = Num)
Boxplot_data$Temperatura <- as.factor(Boxplot_data$Temperatura)
boxplot(Valor ~ Temperatura,
        data = Boxplot_data,
        xlab = "Nivel",
        ylab = "Temperatura (°C)",
        col = c("#F54927","#2A27F5"))

Se observa una dispersión similiar, sin embargo la media cambia de posición.

Ahora vamos a ver un caso que afecta a la dispersión de los datos:

Temp_600 <- rnorm(50, mean = 25, sd = 0.2)
Temp_800 <- rnorm(50, mean = 25, sd = 2.0)

La media y desviación a \(600°C\) y \(800°C\) en este caso es:

mean(Temp_600)
## [1] 25.01791
sd(Temp_600)
## [1] 0.2024354
mean(Temp_800)
## [1] 25.01099
sd(Temp_800)
## [1] 1.999428

Y el boxplot es:

Cat_1 <- rep("600_grados",50)
Cat_2 <- rep("800_grados",50)
Cat <- c(Cat_1,Cat_2)
Num <- c(Temp_600,Temp_800)
Boxplot_data <- data.frame(Temperatura = Cat, Valor = Num)
Boxplot_data$Temperatura <- as.factor(Boxplot_data$Temperatura)
boxplot(Valor ~ Temperatura,
        data = Boxplot_data,
        xlab = "Nivel",
        ylab = "Temperatura (°C)",
        col = c("#F54927","#2A27F5"))

Se observa que la mayoría de los datos están centrados en una localización similar para \(600°C\) y para \(800°C\).

Por último, se muestra el caso en que afecta a la media y la desviación estándar. Este tipo de factores son los que deben tener un mayor control.

Temp_600 <- rnorm(50, mean = 25, sd = 0.2)
Temp_800 <- rnorm(50, mean = 30, sd = 2.0)

La media y desviación estándar para ambos niveles se muestran a continuación:

mean(Temp_600)
## [1] 24.99471
sd(Temp_600)
## [1] 0.1968247
mean(Temp_800)
## [1] 29.93926
sd(Temp_800)
## [1] 2.03048

Y el diagrama de caja y bigote muestra lo siguiente:

Cat_1 <- rep("600_grados",50)
Cat_2 <- rep("800_grados",50)
Cat <- c(Cat_1,Cat_2)
Num <- c(Temp_600,Temp_800)
Boxplot_data <- data.frame(Temperatura = Cat, Valor = Num)
Boxplot_data$Temperatura <- as.factor(Boxplot_data$Temperatura)
boxplot(Valor ~ Temperatura,
        data = Boxplot_data,
        xlab = "Nivel",
        ylab = "Temperatura (°C)",
        col = c("#F54927","#2A27F5"))

Aquí se observa un cambio en la posición de la mayoría de los datos, y también se observa un cambio en la dispersión. Cuando un factor tiene este impacto en la variable de respuesta debe ser tratado con mucho cuidado.

Arreglos ortogonales Taguchi

Existen varios arreglos ortogonales propuestos en la metodología Taguchi. Se va a emplear el siguiente paquete:

library(DoE.base) 
## Loading required package: grid
## Loading required package: conf.design
## Registered S3 method overwritten by 'DoE.base':
##   method           from       
##   factorize.factor 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':
## 
##     lengths

A continuación, se explican algunos de ellos:

L4

Este diseño trabaja con máximo tres factores y 2 niveles.

diseno_L4 <- oa.design(nfactors = 3,  # Máximo de factores para L4
                       nlevels = 2,   # Todos los factores con 2 niveles
                       replications = 2, # Recomendable al menos 2 réplicas
                       randomize = TRUE,
                       seed = 123)

print(diseno_L4)
##   run.no run.no.std.rp A B C
## 1      1           3.1 2 1 2
## 2      2           4.1 2 2 1
## 3      3           1.1 1 1 1
## 4      4           2.1 1 2 2
## 5      5           3.2 2 1 2
## 6      6           2.2 1 2 2
## 7      7           4.2 2 2 1
## 8      8           1.2 1 1 1
## class=design, type= oa 
## NOTE: columns run.no and run.no.std.rp  are annotation, 
##  not part of the data frame

L8

El diseño L8 trabaja máximo con 7 factores y tres niveles.

diseno_L8 <- oa.design(nfactors = 7,  # Máximo de factores para L8
                       nlevels = 2,   # Todos los factores con 2 niveles
                       replications = 1, # Número de réplicas
                       randomize = TRUE,
                       seed = 123)

print(diseno_L8)
##   A B C D E F G
## 1 2 2 1 1 2 1 2
## 2 2 2 2 2 1 1 1
## 3 1 2 1 2 1 2 2
## 4 2 1 2 1 1 2 2
## 5 1 1 2 2 2 1 2
## 6 1 2 2 1 2 2 1
## 7 2 1 1 2 2 2 1
## 8 1 1 1 1 1 1 1
## class=design, type= oa

L9

El diseño L9 trabaja con 4 factores como máximo en 3 niveles.

diseno_L9 <- oa.design(nfactors = 4, 
                       nlevels = 3,
                       replications = 2,  # 2 réplicas para mayor confiabilidad
                       randomize = TRUE,
                       seed = 123)

print(diseno_L9)
##    run.no run.no.std.rp A B C D
## 1       1           3.1 1 3 2 3
## 2       2           6.1 2 3 1 2
## 3       3           9.1 3 3 3 1
## 4       4           2.1 1 2 3 2
## 5       5           8.1 3 2 1 3
## 6       6           5.1 2 2 2 1
## 7       7           7.1 3 1 2 2
## 8       8           1.1 1 1 1 1
## 9       9           4.1 2 1 3 3
## 10     10           6.2 2 3 1 2
## 11     11           1.2 1 1 1 1
## 12     12           2.2 1 2 3 2
## 13     13           3.2 1 3 2 3
## 14     14           5.2 2 2 2 1
## 15     15           9.2 3 3 3 1
## 16     16           4.2 2 1 3 3
## 17     17           8.2 3 2 1 3
## 18     18           7.2 3 1 2 2
## class=design, type= oa 
## NOTE: columns run.no and run.no.std.rp  are annotation, 
##  not part of the data frame

###L12

El diseño L12 trabaja máximo con 11 factores con 2 niveles.

diseno_L12 <- oa.design(nfactors = 11,  # Máximo de factores para L12
                        nlevels = 2,    # Todos los factores con 2 niveles
                        replications = 1,
                        randomize = TRUE,
                        seed = 123)

print(diseno_L12)
##    A B C D E F G H J K L
## 1  1 1 2 1 2 2 2 1 1 1 2
## 2  2 2 2 2 2 2 2 2 2 2 2
## 3  2 2 1 1 1 2 1 1 2 1 2
## 4  1 1 2 1 1 2 1 2 2 2 1
## 5  1 2 2 2 1 1 1 2 1 1 2
## 6  2 2 2 1 1 1 2 1 1 2 1
## 7  1 2 1 2 2 2 1 1 1 2 1
## 8  1 2 1 1 2 1 2 2 2 1 1
## 9  2 1 2 2 2 1 1 1 2 1 1
## 10 2 1 1 2 1 2 2 2 1 1 1
## 11 1 1 1 2 1 1 2 1 2 2 2
## 12 2 1 1 1 2 1 1 2 1 2 2
## class=design, type= oa

Análisis estadístico

El análisis estadístico del Diseño de Experimentos Taguchi se hace para la media y la razón señal ruido (\(S/N\hspace{0.2cm} \text{ratio}\)). Para un mejor entendimiento vamos a desarrollar un ejemplo:

En una empresa que elabora harinas preparadas para pasteles se quiere introducir una nueva mezcla en el mercado. Interesa diseñar la mezcla de tal manera que sea robusta al hecho de que el cliente no se apegue del todo a las instrucciones de la caja. Los factores de diseño son la cantidad de harina (\(H\)), de azúcar (\(A\)), y de huevo (\(E\)). La variable de respuesta es una evaluación subjetiva de los pasteles mediante un panel de jueces usando una escala hedórica de uno a siete. El diseño robusto empleado y los resultados se obtienen a continuación. Utilice la razón señal-ruido del tipo nominal es mejor tipo 2.

H <- c(-1,1,-1,1,-1,1,-1,1)
A <- c(-1,-1,1,1,-1,-1,1,1)
E <- c(-1,-1,-1,-1,1,1,1,1)
replica_1 <- c(3.1,3.2,5.3,4.1,5.9,6.9,3.0,4.5)
replica_2 <- c(1.1,3.8,3.7,4.5,4.2,5.0,3.1,3.9)
replica_3 <- c(5.7,4.9,5.1,6.4,6.8,6.0,6.3,5.5)
replica_4 <- c(6.4,4.3,6.7,5.8,6.5,5.9,6.4,5.0)
replica_5 <- c(1.3,2.1,2.9,5.2,3.5,5.7,3.0,5.4)

Datos <- data.frame(H,A,E,replica_1,replica_2,replica_3,replica_4,replica_5)

Para el análisis estadístico se cargan los siguientes paquetes:

library(DoE.base)
library(FrF2)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)

Se debe obtener la media de las cuatro réplicas para este caso y la razón señal-ruido. La media se obtiene para cada corrida con la siguiente ecuación.

\[\bar{x} = \frac{1}{n}\sum_{i=1}^{n}x_{i}\]

Mientras que la razón señal-ruido se obtiene con la ecuación siguiente:

\[\text{S/N} = -10 \log_{10}(S^{2})\] Para ello se van a generar dos columnas más en el dataframe que son la media y la razón señal-ruido:

Datos <- Datos %>% 
  rowwise() %>%  # Para aplicar las funciones por fila
  mutate(Media = mean(c(replica_1,replica_2,replica_3,replica_4,replica_5)),
         Varianza = var(c(replica_1,replica_2,replica_3,replica_4,replica_5)),
         SN = -10*log10(Varianza)) %>% 
  ungroup()  # Remover el agrupamiento por fila

print(Datos[,c(1,2,3,9,11)])
## # A tibble: 8 × 5
##       H     A     E Media     SN
##   <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1    -1    -1    -1  3.52 -7.78 
## 2     1    -1    -1  3.66 -0.618
## 3    -1     1    -1  4.74 -3.40 
## 4     1     1    -1  5.2   0.580
## 5    -1    -1     1  5.38 -3.26 
## 6     1    -1     1  5.9   3.33 
## 7    -1     1     1  4.36 -5.19 
## 8     1     1     1  4.86  3.54

Se debe cambiar a factores las variables numéricas \(H\), \(A\), y \(E\).

Datos$A <- as.factor(Datos$A)
Datos$H <- as.factor(Datos$H)
Datos$E <- as.factor(Datos$E)

Se realiza un Análisis de Varianza (ANOVA, por sus siglas en inglés) para la media y la razón señal ruido.

ANOVA para la media

ANOVA_media <- aov(Media ~ A + H + E, data = Datos)
summary(ANOVA_media)
##             Df Sum Sq Mean Sq F value Pr(>F)
## A            1 0.0613  0.0613   0.083  0.788
## H            1 0.3281  0.3281   0.445  0.541
## E            1 1.4281  1.4281   1.935  0.237
## Residuals    4 2.9518  0.7380

Se observa un valor-\(p\) alto para los factores \(A\), y \(H\). Se eliminan del modelo para ver si el valor-\(p\) del factor \(E\) es significativo.

ANOVA_media <- aov(Media ~ E, data = Datos)
summary(ANOVA_media)
##             Df Sum Sq Mean Sq F value Pr(>F)
## E            1  1.428  1.4281   2.565   0.16
## Residuals    6  3.341  0.5569

Aunque se eliminen los otros dos factores el factor \(E\) no es significativo

ANOVA para la razón señal-ruido

ANOVA_SN <- aov(SN ~ A + H + E, data = Datos)
summary(ANOVA_SN)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## A            1   1.86    1.86   0.595 0.48348   
## H            1  87.47   87.47  27.960 0.00614 **
## E            1  11.61   11.61   3.711 0.12634   
## Residuals    4  12.51    3.13                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El valor-\(p\) más alto es el del factor \(A\) y el único valor significativo en este modelo es el \(H\). Se elimina del modelo los factores \(A\) y \(E\).

ANOVA_SN <- aov(SN ~ H, data = Datos)
summary(ANOVA_SN)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## H            1  87.47   87.47    20.2 0.00413 **
## Residuals    6  25.98    4.33                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Por lo que, sólo se tiene un factor significativo en el ANOVA para la razón señal-ruido. Este factor es el \(H\).

Gráficos de efectos principales para la media

Aunque el ANOVA para la media no muestra factores significativos, se realizan los gráficos de efectos principales con la Media como variable de respuesta.

Factor A

# Efectos principales (puede utilizarse para media y razon sn)
ggplot(Datos, aes(x = A, y = Media)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", group = 1) +
  labs(title = "Efecto principal del Factor A", 
       x = "Nivel Factor A", y = "Respuesta Promedio (Media)")

Factor H

# Efectos principales (puede utilizarse para media y razon sn)
ggplot(Datos, aes(x = H, y = Media)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", group = 1) +
  labs(title = "Efecto principal del Factor H", 
       x = "Nivel Factor H", y = "Respuesta Promedio (Media)")

Factor E

# Efectos principales (puede utilizarse para media y razon sn)
ggplot(Datos, aes(x = E, y = Media)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", group = 1) +
  labs(title = "Efecto principal del Factor E", 
       x = "Nivel Factor E", y = "Respuesta Promedio (Media)")

Ningún factor afecta la media, por lo que no se pueden utilizar para llevar el proceso a una media deseada.

Gráficos de efectos principales para la razón señal-ruido

En este caso se tiene como significativo el factor \(H\), sin embargo se realizarán los gráficos de efectos principales de todos los factores.

Factor A

# Efectos principales (puede utilizarse para media y razon sn)
ggplot(Datos, aes(x = A, y = SN)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", group = 1) +
  labs(title = "Efecto principal del Factor A", 
       x = "Nivel Factor A", y = "Respuesta Promedio (razón S/N)")

Factor H

# Efectos principales (puede utilizarse para media y razon sn)
ggplot(Datos, aes(x = H, y = SN)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", group = 1) +
  labs(title = "Efecto principal del Factor H", 
       x = "Nivel Factor H", y = "Respuesta Promedio (razón S/N)")

Factor E

# Efectos principales (puede utilizarse para media y razon sn)
ggplot(Datos, aes(x = E, y = SN)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", group = 1) +
  labs(title = "Efecto principal del Factor E", 
       x = "Nivel Factor E", y = "Respuesta Promedio (razón S/N)")

En este caso, el factor \(H\) es significativo con la razón señal-ruido. La optimización se realiza en dos pasos, que se describen a continuación.

Optimización en dos pasos

La optimización en el Diseño de Experimentos Taguchi se realiza en dos pasos que son:

  • Utilizar el nivel de los factores significativos en el ANOVA de la razón señal-ruido que maximicen la misma.
  • Utilizar el nivel de los factores significativos en el ANOVA de la media que lleven al valor deseado la variable de respuesta que se analiza.

Para este ejercicio, sólo se tiene un factor significativo en el ANOVA de la razón señal-ruido que es el factor \(H\). Por lo que, el nivel que maximiza la razón señal-ruido es el alto. De esta manera se tiene un proceso robusto.

Análisis de supuestos

Por último, se realiza un análisis de los supuestos que debe cumplir un Diseño de Experimentos que son normalidad, homocedasticidad e independencia.

Normalidad

# Residuos del modelo ANOVA de la media
residuos_SN <- residuals(ANOVA_SN)

# Prueba de Shapiro-Wilk
shapiro.test(residuos_SN)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_SN
## W = 0.84513, p-value = 0.085
# Q-Q Plot
qqnorm(residuos_SN)
qqline(residuos_SN, col = "red")

Se tiene normalidad en los residuos del modelo porque se tiene un valor-\(p\) mayor que \(0.05\).

Homocedasticidad

# Levene Test para la media
car::leveneTest(SN ~ E, data = Datos)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.8393 0.3949
##        6
# O prueba de Bartlett
bartlett.test(SN ~ E, data = Datos)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  SN by E
## Bartlett's K-squared = 0.09203, df = 1, p-value = 0.7616

Existe homocedasticidad con este modelo porque se tiene un valor-\(p\) mayor que \(0.05\).

Independencia

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# Durbin-Watson
dwtest(ANOVA_SN)
## 
##  Durbin-Watson test
## 
## data:  ANOVA_SN
## DW = 1.4507, p-value = 0.3195
## alternative hypothesis: true autocorrelation is greater than 0
# Gráfico de residuos vs valores ajustados
plot(fitted(ANOVA_SN), residuos_SN,
     xlab = "Valores ajustados", ylab = "Residuos",
     main = "Residuos vs Valores ajustados")
abline(h = 0, col = "red")

Se prueba que no existe correlación entre los residuos porque se tiene un valor-\(p\) mayor que \(0.05\).