Introducción

En el siguiente notebook se encuentra un ejercicio del libro Montgomery, D. (2017). Design and analysis of experiments. Ninth edition. Wiley. Chapter 3: Experiments with a single factor: the analysis of variance, en especifico el 3.26 de este capítulo.

Ejercicio 3.26

The response time in milliseconds was determined for three types of circuits that could be used in an automatic valve shutoff mechanism. The results from a completely randomized experiment are shown in the following table:

# Valores
circuit_type_1 <- c(9,12,10,8,15)
circuit_type_2 <- c(20,21,23,17,30)
circuit_type_3 <- c(6,5,8,16,7)

# Unidos en tabla
data <- data.frame(circuit_type_1,
                   circuit_type_2,
                   circuit_type_3)

# Presentación de la tabla
DT::datatable(t(data),
              caption = 'Tabla 1: Tiempos de respuesta de los circuitos',
              extensions = c('Buttons'),
          options = list(dom='Bfrtip',
                         buttons=c('csv','excel','pdf')),
          escape = T)

(a) Test the hypothesis that the three circuit types have the same response time. Use \(\alpha = 0.01\).

(b) Use Tukey’s test to compare pairs of treatment means. Use \(\alpha = 0.01\).

(c) Use the graphical procedure in Section 3.5.3 to compare the treatment means. What conclusions can you draw? How do they compare with conclusions from part (b)?

(d) Construct a set of orthogonal contrasts, assuming that at the outset of the experiment you suspected the response time of circuit type 2 to be different from the other two.

(e) If you were the design engineer and you wished to minimize the response time, which circuit type would you select?

(f) Analyze the residuals from this experiment. Are the basic analysis of variance assumptions satisfied?

Solución (a)

(a) Test the hypothesis that the three circuit types have the same response time. Use \(\alpha = 0.01\).

Donde asumimos que el tiempo promedio de respuesta es el mismo como hipótesis nula, contra la hipótesis alternativa de que los tiempos de respuesta de al menos uno de los circuitos difiere.

\[H_0:\mu_1=\mu_2=\mu_3 \\ vs\\ H_1: \mu_i\ne \mu_j \] Para al menos un par \((i,j)\) con \(i\ne j\)

Planteando dichas hipótesis para el modelo

\[y_{ij} = \mu + \tau_i + \varepsilon_{ij} \;\;\; \begin{cases} i=1,2,..., a\\ j=1,2,...,n \end{cases}\]

# Acomodar datos para aplicar las funciones apropiadamente
# Tratamiento se crea un columna con la codificación de estos
tipo_circuito <- c(rep("1", 5),
                   rep("2", 5),
                   rep("3", 5))
# Respuestas colmuna con los datos de los tiempos de respuesta
tiempo_respuesta <- c(data$circuit_type_1,
                      data$circuit_type_2,
                      data$circuit_type_3)
# Creación de DataFrame
df <- data.frame(tipo_circuito,
                 tiempo_respuesta)

# Presentación de datos usando 'datatable'
DT::datatable(df,
              caption = 'Tabla 2: Datos en formato apropiado',
              extensions = c('Buttons'),
          options = list(dom='Bfrtip',
                         buttons=c('csv','excel','pdf')),
          escape = TRUE,
          colnames = c("Numeración" =1,
                       "Tipo de Circuito" = 2,
                       "Tiempo de Respuesta" = 3)
          )

Observaciones previas

# Boxplot ggplot2
ggplot2::ggplot(data = df, ggplot2::aes(x = tipo_circuito, y = tiempo_respuesta)) +
  ggplot2::stat_boxplot(geom = "errorbar", #bigotes
                        width = 0.2) +
  ggplot2::geom_boxplot(alpha = 0.9, outlier.colour = "red") + 
  ggplot2::scale_y_continuous(name = "Tiempo de Respuesta") +
  ggplot2::scale_x_discrete(name = "Tipo de Circuito") + 
  ggplot2::ggtitle("Boxplot")

Donde observamos que pareciera existir diferencias entre estos, ademas de que es importante notar la existencia de datos atípicos en el tipo de circuito \(2\) y \(3\), y ademas observamos que la mayor diferencias están entre los circuitos \(1\) y \(3\) con respecto al circuito \(2\)

Para terminar de confirmar dicha intuición gráfica realizamos el análisis de varianzas.

# Función del ANOVA
circuitos.aov <- aov(formula = tiempo_respuesta ~ tipo_circuito,
                     data = df)
# Datos del ajuste
summary(circuitos.aov)

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## tipo_circuito  2  543.6   271.8   16.08 0.000402 ***
## Residuals     12  202.8    16.9                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Donde primero observamos que la suma de cuadrados de los tratamientos (tipo de circuitos) es superior la suma de cuadrados de los errores, con lo cual concluimos que si existe una variación entre los grupos.

También los cuadrados del tratamiento (tipos de circuito) es superior que los cuadrados medios de los errores, y estas difieren entre ellas y considerando que estos son estimadores de la varianza nos da indicativos de que sus varianzas (variación entre los grupos) son diferentes lo que nos puede llevar a un posible rechazo de la \(H_0\)

Por ultimo observamos que el valor \(F\) se aleja de 1, lo cual como observamos con el p-valor nos conduce a rechazar la hipótesis nula con un nivel de significancia \(\alpha = 0.01\) esto es, que con un nivel de confianza del 99% concluimos que si hay diferencia en el tiempo de respuesta promedio de los circuitos del mecanismo de cierre de válvulas.

Otra visualización de este con ggstatsplot es

# ggstatsplot
ggstatsplot::ggbetweenstats(
  data = df,
  x = tipo_circuito,
  y = tiempo_respuesta,
  type = "parametric",
  pairwise.display = "significant",
  conf.level = 0.99,
  ggtheme = ggplot2::theme_gray()
  ) +
  ggplot2::ylab("Tiempos de Respuesta") +
  ggplot2::xlab("Tipo de Circuito")

El cual usa un valor \(F_{Welch}\) el cual ajusta el denominador del estadístico \(F\) de tal forma que este tenga la misma expectativa que el numerador cuando la hipótesis nula es cierta, el cual sigue bastante alejado de 1, y con un p-valor de \(0.000466\) con un nivel de significancia del \(\alpha=0.01\) se rechaza la hipótesis nula, esto es que con un nivel de confianza del 99% se concluye que existen diferencias entre los tiempos de respuesta promedios de los circuitos.

Con respecto al tamaño del efecto donde usando el estudio de Field(2013) siendo mayor este a 0.14 es lo suficientemente grande para concluir que el efecto de los tipos de circuito en los tiempo de respuesta de cierre de las válvulas es considerable.

Observando las estadísticas en la parte inferior del gráfico (Bayesianas) nos indica que el posterior que es equivalente a un coeficiente de determinación es de 0.66 y un estudio de Cohen(1998) considera que si este es superior a 0.26 es substancial, lo que nos indica que el ajuste es bueno (desde la perspectiva Bayesiana) y el Bias Factor, que es el equivalente al p-value Frecuentista nos indica -4.44, lo cual es considerado como fuerte evidencia para rechazar a \(H_0\).

(Se que la ULA no es una escuela Bayesiana, solo se mencionan estas métricas como comentarios)

Una ventaja del paquete ggstatsplot es la cantidad de información condensada que nos da en un gráfico publicable como indica la parte derecha del gráfico la función corre automáticamente la prueba Games-Howell y corrige los p-valores del problema de comparaciones múltiples con el método Holm-Bonferroni, que se ven en el gráfico (0,02 y 0,01) respectivamente.

Donde muestra en particular que los circuitos \(1-2\) y \(3-2\) son los que muestran ser significativamente diferentes.

Solución (b)

(b) Use Tukey’s test to compare pairs of treatment means. Use \(\alpha = 0.01\).

Donde planteamos que

\[H_0 : \mu_1=\mu_2=\mu_3\\ vs\\ H_1:\mu_i\ne\mu_j \\ i\ne j\]

# Función del test
post_hoc <- TukeyHSD(circuitos.aov,
                    conf.level = 0.99)
post_hoc

##   Tukey multiple comparisons of means
##     99% family-wise confidence level
## 
## Fit: aov(formula = tiempo_respuesta ~ tipo_circuito, data = df)
## 
## $tipo_circuito
##      diff        lwr       upr     p adj
## 2-1  11.4   2.123163 20.676837 0.0023656
## 3-1  -2.4 -11.676837  6.876837 0.6367043
## 3-2 -13.8 -23.076837 -4.523163 0.0005042

Donde primero observamos que la diferencia de las medias de los circuitos donde observamos que la mas baja es la diferencia de la media del circuito \(3\) y el circuito \(1\) es la mas baja y las diferencias de medias de los circuitos \(2-1\) y \(3-2\) son las mas altas, algo que ya se observaba en los Boxplots.

Ahora si observamos los p-valores directamente vemos como precisamente las comparaciones de \(2-1\) y \(3-2\) que son precisamente las que se veían mas claras en los Boxplots son las que presentan un p-valor inferior al \(0.01\) por lo tanto con un nivel de confianza del 99% concluimos que entre los circuitos \(2-1\) y \(3-2\) existen diferencias significativas en los tiempos de respuesta del mecanismo de cierre de las válvulas, donde la diferencia de las medias entre \(3-2\) parece ser superior.

Para los intervalos de confianza estos se interpretan mejor con el siguiente gráfico

plot(post_hoc)
grid(lty = 3,
     col = "grey")

Donde se ve claramente que es el grupo \(3-2\) y \(2-1\) las comparaciones que presentan diferencias entre sus medias, ya que el cero no esta en estos.

Solución (c)

(c) Use the graphical procedure in Section 3.5.3 to compare the treatment means. What conclusions can you draw? How do they compare with conclusions from part (b)?

La lógica de este método es que gracias al teorema central del limite conocemos que la distribución de los estadísticos \(\bar{y}_{i}\) siendo esta de \(\bar{y}_{i} \sim N(\bar{y}_{..},\sigma/\sqrt{n})\)

Donde se aproxima la varianza con \(\sqrt{MS_E/n}\)

Siendo esta de

desv <- sqrt(16.9/5) 
desv

## [1] 1.838478

Así si los tratamientos son iguales las distribuciones deberían de estar en una posición donde se vea que estas provienen de la misma distribución.

# Semilla
set.seed(1234)
# Simulación de 500 observaciones
treat1 <- rnorm(500, mean(circuit_type_1), desv)
treat2 <- rnorm(500, mean(circuit_type_2), desv)
treat3 <- rnorm(500, mean(circuit_type_3), desv)

# dataframe de simulaciones
dfs <- data.frame(treat1,
                  treat2,
                  treat3)
# Uso de reshape para tenerlos en un solo dataframe de dos columnas, (es mas rapido con ese paquete, para no repetir el proceso de como el dataframe del análisis de varianza de la pregunta 'a')
dfs.plot <- reshape::melt(dfs)
#Observación
DT::datatable(dfs.plot)

Así con los datos en un formato mas cómodo para ggplot2, gráficamos

# Gráfico usando ggplot2
ggplot2::ggplot(dfs.plot, ggplot2::aes(x = value, fill = variable)) + 
  ggplot2::geom_density() +
  ggplot2::ggtitle("Densidades") + 
  ggplot2::ylab("Densidad") +
  ggplot2::xlab("Observaciones")

Donde observamos como los tratamientos \(1-3\) al gráficar sus estadísticos no tienen una diferencia visual significativa lo que nos da a entender que provienen de la misma población pero observamos como las densidades de los tratamientos \(2-1\) y \(3-2\) difieren.

Siendo estas las mismas conclusiones que se obtuvieron al realizar el ANOVA mas la prueba de Tukey.

Solución (d)

(d) Construct a set of orthogonal contrasts, assuming that at the outset of the experiment you suspected the response time of circuit type 2 to be different from the other two.

Si sospechamos que el tiempo de respuesta 2 es el que presenta diferencias de los otros dos, el contraste ortogonal vendría siendo de la siguiente forma:

# Declaramos los tratamientos/tipo de circuitos como factores
fact <- factor(df$tipo_circuito)
contrasts(fact)

##   2 3
## 1 0 0
## 2 1 0
## 3 0 1

Planteamos el contraste

\[H_0 : \mu_1 -2\mu_2 + \mu_3 =0\\ H_0: 2\mu_2 = \mu_1 + \mu_3 \\ H_1 : 2\mu_2 \ne \mu_1 + \mu_3\\ \Rightarrow\\ C_1 = y_{1.}-2y_{2.}+y_{3.}\]

# Vector ortogonal 
m <- matrix(nrow = 3, ncol = 2,
            c(1, -2, 1,
              0,0,0)
            )
# Nombramos el vector c_1
colnames(m) <- c("C1", "C2")
#reemplazando
contrasts(fact) <- m
contrasts(fact)

##   C1 C2
## 1  1  0
## 2 -2  0
## 3  1  0

# Función del modelo
circuitos.ortogonales.aov <- aov(formula = tiempo_respuesta ~ fact,
                                 data = df)
# Bondad del ajuste
summary(circuitos.ortogonales.aov)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## fact         1  529.2   529.2   31.67 8.23e-05 ***
## Residuals   13  217.2    16.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Así al contrastar que el circuito \(2\) contra el \(1\) y el \(3\) concluimos que la suma de cuadrados del tratamiento es superior a la del residuo mostrando que si existe variación sistema en en estos, ademas de que sus cuadrados medios difieren mostrando diferencia de medias (al ser estimadores de la varianza), con valor \(F\) bastante alejado de \(1\) y un p-valor menor que \(\alpha = 0.01\) concluimos que si hay suficiente evidencia en la muestra para rechazar la hipótesis nula, esto es, que los tiempos promedios de los circuitos \(1\) y \(3\) difieren del circuito \(2\)

Solución (e)

(e) If you were the design engineer and you wished to minimize the response time, which circuit type would you select?

Como observamos en el Boxplot anterior (el codigo esta en la primera pregunta)

Escogería el \(1\) o el \(3\), al ver la gráfica nos sentiríamos tentados a decir que el \(3\) mas el análisis de varianza nos muestra que entre los circuitos \(1\) y \(3\) no existen diferencias significativas.

Como observamos en la prueba de Tukey (código en el ejercicio b)

Por lo cual se puede escoger cualquiera de esos dos.

Solución (f)

(f) Analyze the residuals from this experiment. Are the basic analysis of variance assumptions satisfied?

Homocedasticidad de los residuos

plot(circuitos.aov$residuals, 
     main = "Residuos",
     xlab = "i",
     ylab = "Residuos")
grid(lty = 1)

Donde observamos un patrón ligeramente creciente al final de este, pero no es claro como para decir que estos tengan un patrón marcado.

Normalidad

ggplot2::ggplot(data = df, ggplot2::aes(x = circuitos.aov$residuals)) +
  ggplot2::geom_histogram(ggplot2::aes(y=..density..),
                          binwidth = 1) +
  ggplot2::geom_density() +
  ggplot2::ggtitle("Histograma de los residuos")+
  ggplot2::xlab("Residuos")+
  ggplot2::ylab("Frecuencia")

Donde claramente se observa como no estan distribuidos normalmente algo que confirmamos con el QQ-Plot

ggpubr::ggqqplot(circuitos.aov$residuals,
                 conf.int.level = 0.99,
                 ggtheme = ggplot2::theme_gray(),
                 title = "QQ-Plot",
                 ylab = "Valores Teoricos",
                 xlab = "Valores de Muestra")

Donde observamos que estos no se quedan en la franja de confianza del 99%.

Valores atípicos

Donde

Como observamos con anterioridad existen valores atípicos en el tratamiento \(2\) y \(3\).

Autocorrelación

acf(circuitos.aov$residuals,
    main = "Autocorrelación",
    xlab= "Rezagos",
    ylab = "ACF")
grid()

Donde observamos que no parecen estar autocorrelacionados, dado que ninguno parece ser significativo.

ggplot2::ggplot(data = df, ggplot2::aes(x=circuitos.aov$fitted.values,
                                        y=circuitos.aov$residuals)) +
  ggplot2::geom_point() + 
  ggplot2::ggtitle("Residuos vs Valores ajustados") + 
  ggplot2::xlab("Valores ajustados") + 
  ggplot2::ylab("Residuos")

Donde observamos lo que esperamos, estos tampoco forman ningún tipo de patrón, inesperado.

Conclusión

Se cumplen todos los supuestos, pero se pone en duda la normalidad, aún así como discutimos en clases, aún así la normalidad este en duda, teniendo sus rezagos incorrelacionados y homogéneos puede seguir siendo un buen modelo.

Conclusión

Así se termina el ejercicio en el cual se puso en practica el ANOVA de efectos fijos.