Diseño en Cuadro Latino
M.I.I. Gabriel Grosskelwing Núñez
Caso de Análisis
Se quiere estudiar el efecto de cinco diferentes catalizadores (A,B,C,D,E) sobre el tiempo de reacción de un proceso químico. Cada lote de material sólo permite cinco corridas y cada corrida requiere aproximadamente 1.5 horas, por lo que solo se pueden realizar cinco corridas diarias. El experimentador decide correr los experimentos con un diseño en cuadro latino (DCL)npara controlar activamente los lotes y días. Los datos obtenidos se muestran en la siguiente matriz:
tratamiento=c("A=8","C=11","B=4","D=6","E=4","B=7","E=2","A=9","C=8","D=3",
"D=1","A=7","C=10","E=6","B=3","C=7","D=3","E=1","B=6","A=8",
"E=3","B=8","D=5","A=10","C=8")
matriz_tratamiento=matrix(tratamiento,nrow = 5)
dimnames(matriz_tratamiento)=list(c("1","2","3","4","5"),c("1","2","3","4","5"))
print(matriz_tratamiento)## 1 2 3 4 5
## 1 "A=8" "B=7" "D=1" "C=7" "E=3"
## 2 "C=11" "E=2" "A=7" "D=3" "B=8"
## 3 "B=4" "A=9" "C=10" "E=1" "D=5"
## 4 "D=6" "C=8" "E=6" "B=6" "A=10"
## 5 "E=4" "D=3" "B=3" "A=8" "C=8"Dado lo anterior, realice el análisis estadístico correspondiente y establezca las conclusiones adecuadas.
Modelo Matemático
En este caso, como lo que se persigue comprender el comportamiento del rendimiento del proceso en función de los diferentes catalizadores, considerando dentro del análisis a dos variables de ruido, como lo son, el Lote (1,2,3,4,5) y el día de la semana en que se tomaron las muestras (1,2,3,4,5), las dos variables consideradas en el análsis no pueden ser aleatorizadas, razón por la cual se consideran bloques. Dado lo anterior, se define el modelo matemática de la siguiente manera:
\[y_{ijk}=\mu+\tau_{i}+\beta_{j}+\delta_k+\varepsilon_{ijk}\] Comprendiendo ya el modelo matemático, procederemos a ejecutar el siguiente código, mismo que nos dará los como resultado los coeficientes de Mínimos Cuadrados para el caso de análisis, considerando los datos, organizado en columnas:
library(readxl)
datos=read_excel("dataset.xlsx")
attach(datos)
return(datos)Una vez incluidos los datos, es importante indicar que las variables son factores, para que sean reconocidos como categorías, esto de logra mediante la siguiente secuencia de comandos:
f_trataiento=factor(Catalizador)
f_bloque1=factor(Lote)
f_bloque2=factor(Dia)Con estos datos, ya considerados como factores, se procede a establecer el modelo matemático con la siguientes líneas de comando:
modelo=lm(Tiempo~(f_trataiento+f_bloque1+f_bloque2))
summary(modelo)##
## Call:
## lm(formula = Tiempo ~ (f_trataiento + f_bloque1 + f_bloque2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.24 -1.24 -0.24 0.96 2.36
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.440 1.275 6.619 2.47e-05 ***
## f_trataientoB -2.800 1.118 -2.504 0.027725 *
## f_trataientoC 0.400 1.118 0.358 0.726798
## f_trataientoD -5.000 1.118 -4.471 0.000764 ***
## f_trataientoE -5.200 1.118 -4.650 0.000561 ***
## f_bloque12 1.000 1.118 0.894 0.388803
## f_bloque13 0.600 1.118 0.537 0.601409
## f_bloque14 2.000 1.118 1.788 0.098971 .
## f_bloque15 -0.200 1.118 -0.179 0.861049
## f_bloque22 -1.000 1.118 -0.894 0.388803
## f_bloque23 -1.200 1.118 -1.073 0.304364
## f_bloque24 -1.600 1.118 -1.431 0.178037
## f_bloque25 0.200 1.118 0.179 0.861049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.768 on 12 degrees of freedom
## Multiple R-squared: 0.8184, Adjusted R-squared: 0.6369
## F-statistic: 4.507 on 12 and 12 DF, p-value: 0.007157De aqui es importante resaltar que, aunque el Coeficiente de Determinación indica que la variablidad del tiempo del proceso está explicada en un 81.84% por los factores considerados en el análisis, alguno de ellos, no está aportando de manera significativa, dado que el Coefciente de Determinación Ajustado reporta que la variabilidad del tiempo del proceso está explicada de manera efectiva sólo en un 63.69% por los factores considerados en el análisis.
Análisis Estadístico
Dado lo anterior, prcederemos a ver si existen diferencias significativas entre los niveles del tratamiento, y a establecer si los factores de bloque tienen alguna participación importante y, en caso de tenerla, cuantificar la participación de los mismod dentro del proceso. Para tal fin, es necesario realizar el análsis de varianza correspondiente,la cual tendrá las siguientes hipótesis a probar: \[ H_0:\tau_i=\tau_j=0 \] \[ H_1:\tau_i\neq\tau_j\neq0 \] Para el caso del factor de tratamiento, para el caso de los factores de bloque, solo se procederá a concluir de manera parcial, en caso de ser significativos, que los esfuerzos por aislar las variables hicieron más sensible el análsis estadístico. La Tabla ANOVA se procesa mediante la siguiente secuencia de comandos:
anova=aov(modelo)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## f_trataiento 4 141.44 35.36 11.309 0.000488 ***
## f_bloque1 4 15.44 3.86 1.235 0.347618
## f_bloque2 4 12.24 3.06 0.979 0.455014
## Residuals 12 37.52 3.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Según los resultados reportados por la Tabla ANOVA, el factor de tratamiento, es decir, los diferentes catalizadores que se probaron, tienen efectos significativos en el tiempo de reacción del proceso químico, mientras que los factores de bloque, es decir el Lote y el dia de en que se toman las muestras no tienen efectos significativos.
Para el caso del factor de tratamiento, se procede a analizar, mediante pruebas de rango múltiple, cuales son los niveles del factor que minimizan el tiempo de reacción del proceso, mediante la siguiente secuencia:
library(agricolae)
prueba_rango=duncan.test(y=Tiempo,trt = Catalizador, MSerror = 3.13,DFerror = 12,alpha = 0.05)
print(prueba_rango)## $statistics
## MSerror Df Mean CV
## 3.13 12 5.88 30.08811
##
## $parameters
## test name.t ntr alpha
## Duncan Catalizador 5 0.05
##
## $duncan
## Table CriticalRange
## 2 3.081307 2.437935
## 3 3.225244 2.551818
## 4 3.312453 2.620818
## 5 3.370172 2.666485
##
## $means
## Tiempo std r Min Max Q25 Q50 Q75
## A 8.4 1.140175 5 7 10 8 8 9
## B 5.6 2.073644 5 3 8 4 6 7
## C 8.8 1.643168 5 7 11 8 8 10
## D 3.4 2.073644 5 1 6 2 3 5
## E 3.2 1.923538 5 1 6 2 3 4
##
## $comparison
## NULL
##
## $groups
## Tiempo groups
## C 8.8 a
## A 8.4 a
## B 5.6 b
## D 3.4 b
## E 3.2 b
##
## attr(,"class")
## [1] "group"plot.group(prueba_rango,main="Prueba de Rango Múltiple para el factor Catalizador")
De aquí se puede observar que en general, los niveles del factor catalizador que más tiempo de reacción del proceso químico producen son el C y el A, los cuales, con 8.8 y 8.4 minutos, respectivamente, no tienen diferencias significativas con un 95% de confianza, sin embargo para el caso de los catalizadores B, D y E, los cuales producen menor tiempo de reacción del proceso químico, no se encuentran diferencias significativas a 95% de confianza, lo cual quiere decir que los tres último son igualemente recomendables si se desea reducri el tiempo de reacción del proceso químico.
Pruebas de Adecuación
En lo general, un modelo estadístico, para ser considerado como adecuado, debe cumplir con dos condicionantes, a saber:
1. Que los residuales del modelo provengan de una población con distribución normal copn media igual a cero y varianza constante.
2. Que los residuales sean homocedásticos.
Prueba de Normalidad
Para el caso específico de la normalidad, procederemos a utilizar la Prueba de Bondad de Ajuste a la Distribución Normal de Shapiro-Wilks, aunque cabe mencionar que puede utilizarse cualquier otra, como Kolmogorov-Smirnov o Anderson-Darlind, la prueba de Shapiro-Wilks está diseñada específicamente para la Distribución Normal.
Para realizar esta prueba debemos plantear las siguientes hipótesis de trabajo:
\[ H_0: {x \in N(\mu=0,\sigma^2=Constante)} \] \[ H_1: {x \not\in N(\mu=0,\sigma^2=Constante)} \] La hipótesis nula de la Prueba de Shapiro-Wilks se rechaza cuando \(Valor_p<0.05\), siempre es recomendable realizar una gráfica de probabilidad normal para confirmar la prueba de hipótesis, para lo cual utilizaremos la siguiente secuencia de comandos:
normalidad=shapiro.test(resid(modelo))
print(normalidad)##
## Shapiro-Wilk normality test
##
## data: resid(modelo)
## W = 0.96606, p-value = 0.5476#Gráfica de Probabiliad Normal
qqnorm(resid(modelo), main= "Gráfica de Probabilidad para los Residuales del Modelo", xlab="Cuantiles Teoricos", ylab = "Cuantiles de muestra")
qqline(resid(modelo))
Como es observarse, se acepta la hipótesis nula de la Prueba de Shapiro Wilks, y se concluye que los residuales provienen de una problación con Distribución de Probabilidad Normal con \(\mu=0\), conclusion que es confirmada por la Gráfica de Probabilidad Normal.
Prueba de Homocedasticidad
La Prueba de Homocedasticidad, o de Varianzas Iguales, se realiza mediante la Prueba de Bartlett, misma que para el caso particular, como sólo es significativo el factor de tratamiento, solo se considerará éste para la ejecución de la prueba.
La Prueba de Bartlett tiene las siguientes hipótesis:
\[ H_0:\sigma^2_i=\sigma^2_j=Constante \] \[ H_1:\sigma^2_i\neq\sigma^2_j\neq Constante \] La Prueba de Bartlett se realiza con la siguiente línea de comandos:
homocedasticidad=bartlett.test(resid(modelo)~Catalizador,data=datos)
print(homocedasticidad)##
## Bartlett test of homogeneity of variances
##
## data: resid(modelo) by Catalizador
## Bartlett's K-squared = 2.6529, df = 4, p-value = 0.6175DAdo que el \(Valor_p>0.05\), se acepta la hipótesis nula y se concluye que las varianzas de los residuales del modelo matemático son constantes, y, dado que se cumplen las dos condiciones, se declara al modelo matemático como adecuado para realizar las estimaciones pertinentes del tiempo de reacción del proceso químico en función del catalizador aplicado.