Taller Final - Diseño y Analisis de Experimentos
Jonathan Martinez - Efrain Martelo
15 de mayo 2020
DISEÑO CUADRADO GRECOLATINO
Para este ejercicio se tratará el problema desarrollado en la sección 4.2 del libro de “Diseño y Analisis de Experimentos” de Montgomery.
Para esto, se procederá inicialmente con el desarrollo del Diseño Cuadrado Latino y luego se transformará en un Diseño Grecolatino.
Funciones utlizadas
Las siguientes son un listado de funciones y su breve descripción que se usarán en el desarrollo del modelo en Rstudio:
shapiro.test: Realiza la prueba de normalidad de Shapiro-Wilk. Su argumento es el conjunto de datos que se le quiere aplicar la prueba. Viene por defecto en el paquete stats el cual está “pre-instalado” en R lo que significa que no se debe instalar ni llamar ninguna librería para su uso.
lm: Se utiliza para ajustar modelos lineales. Se puede utilizar para llevar a cabo regresión, análisis de varianza de un solo estrato y análisis de covarianza (aunque aov puede proporcionar una interfaz más conveniente para estos).
- aov: Se usa para el ajuste de modelo de análisis de varianza mediante una llamada a lm para cada estrato.
ggplot: Es una librería especializada en realización de gráficos.
TukeyHSD: Crea un conjunto de intervalos de confianza sobre las diferencias entre las medias de los niveles de un factor. Los intervalos se basan en la estadística de rango Studentizado, el método de “Diferencia significativa honesta” de Tukey.
LSD.test: Viene incorporado en la librería agricolae, se requiere descargarla y activarla para su uso. Se usa para realizar múltiples comparaciones de tratamientos, mediante el LSD de Fisher (diferencia menos significativa) y da los valores p ajustados.
HSD.test: Viene incorporado en la librería agricolae, se requiere descargarla y activarla para su uso. realiza múltiples comparaciones de tratamientos mediante Tukey. El nivel por defecto alfa es 0.05.
SNK.test: Viene incorporado en la librería agricolae, se requiere descargarla y activarla para su uso. SNK se deriva de Tukey, pero es menos conservador (encuentra más diferencias). Tukey controla el error para todas las comparaciones, donde SNK solo controla las comparaciones bajo consideración. El nivel por defecto alfa es 0.05.
scheffe.test: Viene incorporado en la librería agricolae, se requiere descargarla y activarla para su uso. El método es muy general, ya que todos los posibles contrastes se pueden probar para determinar su importancia y se pueden construir intervalos de confianza para el lineal correspondiente. La prueba es conservadora.
duncan.test: Viene incorporado en la librería agricolae, se requiere descargarla y activarla para su uso. Esta prueba está adaptada del método Newman-Keuls. La prueba de Duncan no controla la tasa de error familiar en el nivel alfa especificado. Tiene más potencia que las otras pruebas posteriores, pero solo porque no controla la tasa de error correctamente. La tasa de error experimental en: 1- (1-alfa) ^ (a-1); donde “a” es el número de medias y es la tasa de error por comparación. El procedimiento de Duncan es solo un poco más conservador que el LSD. El nivel por defecto alfa es 0.05.
fitted: es una función genérica que extrae los valores ajustados de los objetos devueltos por las funciones de modelado.
bartlett.test: Realiza la prueba de Bartlett de igualdad de varianzas en cada uno de los grupos (muestras).
Levene.test: Viene incorporado en la librería car, se requiere descargarla y activarla para su uso. La prueba de Levene se utiliza para evaluar la probabilidad de la igualdad varianza (de los residuos) en diferentes grupos.
Ejercicio
Un experimentador estudia los efectos que tienen cinco Formulaciones diferentes de la carga propulsora utilizada en los sistemas de expulsión de la tripulación de un avión basado en la rapidez de combustión. Cada formulación se hace con un lote de materia prima diferente que sólo alcanza para probar cinco formulaciones. Además, las formulaciones son preparadas por varios Operadores, y puede haber diferencias sustanciales en las habilidades y experiencia en cada uno de ellos; por tanto, al parecer, hay dos factores perturbadores que serán “calculados en promedio” en dicho diseño: los lotes de materia prima y los operadores. El diseño apropiado para este problema consiste en probar cada formulación exactamente una vez en cada uno de los operadores.
Para esto, se parte del diseño de un Cuadrado Latino.
Luego de realizar el análisis del experimento de cuadrado latino anterior, se supone que existe un factor adicional, Los montajes de prueba, que podría ser importante. Este factor considera 5 montajes de prueba diferentes denotados por letras griegas. En la siguiente tabla se muestra el diseño cuadrado grecolatino resultante:
Desarrollo del Ejercicio en Rstudio
Diseño Cuadrado Latino
Comenzaremos diseñando el Cuadrado Latino, de este nace el Grecolatino.
Primero creamos los vectores teniendo en cuenta que los factores se deben crear como factor.
Formulacion = factor(c("A","B","C","D","E",
"B","C","D","E","A",
"C","D","E","A","B",
"D","E","A","B","C",
"E","A","B","C","D"))
Rapidez =c(-1,-8,-7,1,-3,
-5,-1,13,6,5,
-6,5,1,1,-5,
-1,2,2,-2,4,
-1,11,-4,-3,6)
Lote =factor( c(rep("1",1),
rep("2",1),
rep("3",1),
rep("4",1),
rep("5",1)))
Operador = factor(c(rep("1",5),
rep("2",5),
rep("3",5),
rep("4",5),
rep("5",5)))Creamos a partir de nuestro diseño de cuadrado grecolatino, una tabla estructurada, o dataframe con el propósito de correr nuestro experimento
data = data.frame(Rapidez,Formulacion,Lote,Operador)
Lote=data$Lote
data## Rapidez Formulacion Lote Operador
## 1 -1 A 1 1
## 2 -8 B 2 1
## 3 -7 C 3 1
## 4 1 D 4 1
## 5 -3 E 5 1
## 6 -5 B 1 2
## 7 -1 C 2 2
## 8 13 D 3 2
## 9 6 E 4 2
## 10 5 A 5 2
## 11 -6 C 1 3
## 12 5 D 2 3
## 13 1 E 3 3
## 14 1 A 4 3
## 15 -5 B 5 3
## 16 -1 D 1 4
## 17 2 E 2 4
## 18 2 A 3 4
## 19 -2 B 4 4
## 20 4 C 5 4
## 21 -1 E 1 5
## 22 11 A 2 5
## 23 -4 B 3 5
## 24 -3 C 4 5
## 25 6 D 5 5Creamos el modelo lineal a traves de la función lm
modeloLatino= lm(Rapidez ~ Formulacion + Lote + Operador, data)Obtenemos la tabla Anova con dos funciones, primero, la función aov aplicada al modelo lineal, y luego obteniendo el summmary de ese resultado
anovaLatino=aov(modeloLatino)
summary(anovaLatino)## Df Sum Sq Mean Sq F value Pr(>F)
## Formulacion 4 330 82.50 7.734 0.00254 **
## Lote 4 68 17.00 1.594 0.23906
## Operador 4 150 37.50 3.516 0.04037 *
## Residuals 12 128 10.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusiones del diseño LATINO:
- Existe una diferencia significativa en la rapidez de combustión media generada por las diferentes formulaciones de la carga propulsora.
- Existe evidencia estadística para inferir que tambien hay diferencia entre los operadores.
- No hay evidencia estadística para inferir que existe diferencia entre los lotes de materia prima.
Diseño Cuadrado Grecolatino
Con un nuevo asesoramiento, el experimentador determina que existe un factor adicional, Los montajes de prueba, que podrían ser importantes, como existen 5 montajes de prueba, se diseña un cuadrado grecolatino y se le asigna a cada montaje una letra griega que en la tabla denotaremos con la letra minuscula de la “a” a la “e”.
El diseño cuadrado grecolatino usa un analisis de varianza similar al diseño cuadrado latino, y las letras griegas se organizan de forma similar a las letras latinas en el modelo (una vez en cada renglón y columna) de la siguiente manera:
Tomado del libro “Diseño y Analisis de Experimentos” de Montgomery.
El modelo estadístico del diseño es el siguiente:
Tomado del libro “Diseño y Analisis de Experimentos” de Montgomery.
Teniendo en cuenta lo anterior, se ingresa ahora el vector con el factor “Montaje de prueba”.
Montajes =factor(c("b","c","d","e","f",
"d","e","f","b","c",
"f","b","c","d","e",
"c","d","e","f","b",
"e","f","b","c","d")) Se obtiene ahora la nueva tabla de datos del modelo
dataGreco = data.frame(Rapidez,Formulacion,Lote,Operador,Montajes)
dataGreco## Rapidez Formulacion Lote Operador Montajes
## 1 -1 A 1 1 b
## 2 -8 B 2 1 c
## 3 -7 C 3 1 d
## 4 1 D 4 1 e
## 5 -3 E 5 1 f
## 6 -5 B 1 2 d
## 7 -1 C 2 2 e
## 8 13 D 3 2 f
## 9 6 E 4 2 b
## 10 5 A 5 2 c
## 11 -6 C 1 3 f
## 12 5 D 2 3 b
## 13 1 E 3 3 c
## 14 1 A 4 3 d
## 15 -5 B 5 3 e
## 16 -1 D 1 4 c
## 17 2 E 2 4 d
## 18 2 A 3 4 e
## 19 -2 B 4 4 f
## 20 4 C 5 4 b
## 21 -1 E 1 5 e
## 22 11 A 2 5 f
## 23 -4 B 3 5 b
## 24 -3 C 4 5 c
## 25 6 D 5 5 dPara el próximo paso debemos descargar los siguientes paquetes (gridExtra) y (ggplot2)
library(gridExtra)
library(ggplot2)Gráficos descriptivos
Los siguientes son gráficos descriptivos de los factores con la variable respuesta
Form <- ggplot(dataGreco, aes(x = Formulacion, y = Rapidez, fill=Formulacion)) +
geom_boxplot() + theme(legend.position = "none")
Mont <- ggplot(dataGreco, aes(x = Montajes, y = Rapidez, fill=Montajes)) +
geom_boxplot() + theme(legend.position = "none")
Lot <- ggplot(dataGreco, aes(x = Lote, y = Rapidez, fill=Lote)) +
geom_boxplot() + theme(legend.position = "none")
Oper <- ggplot(dataGreco, aes(x = Operador, y = Rapidez, fill=Operador)) +
geom_boxplot() + theme(legend.position = "none")
grid.arrange(Form,Mont,Lot,Oper, nrow=2,ncol=2)
Con las funciones vistas anteriormente, se obtiene el nuevo modelo lineal y su respectiva tabla ANOVA
modeloGreco= lm(Rapidez ~ Formulacion + Lote + Operador + Montajes, dataGreco)
anovaGreco=aov(modeloGreco)
summary(anovaGreco)## Df Sum Sq Mean Sq F value Pr(>F)
## Formulacion 4 330 82.50 10.000 0.00334 **
## Lote 4 68 17.00 2.061 0.17831
## Operador 4 150 37.50 4.545 0.03293 *
## Montajes 4 62 15.50 1.879 0.20764
## Residuals 8 66 8.25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusiones GRECOLATINO:
Se puede inferir basado en evidencia estadística que las formulaciones (Letras latinas) son significativas, es decir, existe alguna diferencia entre ellas; El error experimental disminuye al sacar la variabilidad debida a los montajes de prueba cuando se pasa al diseño Grecolatino, disminuyen los grados de libertad.
Existe evidencia estadística para inferir que los montajes de prueba no son significativos, puesto que no existe diferencia entre cada uno de ellos (no se rechaza la hipotesis de que sus efectos medios son iguales)
Comparaciones múltiples
TukeyHSD(anovaGreco)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = modeloGreco)
##
## $Formulacion
## diff lwr upr p adj
## B-A -8.4 -14.675865 -2.12413503 0.0107829
## C-A -6.2 -12.475865 0.07586497 0.0529245
## D-A 1.2 -5.075865 7.47586497 0.9596565
## E-A -2.6 -8.875865 3.67586497 0.6270350
## C-B 2.2 -4.075865 8.47586497 0.7462701
## D-B 9.6 3.324135 15.87586497 0.0048412
## E-B 5.8 -0.475865 12.07586497 0.0714956
## D-C 7.4 1.124135 13.67586497 0.0218416
## E-C 3.6 -2.675865 9.87586497 0.3525524
## E-D -3.8 -10.075865 2.47586497 0.3087315
##
## $Lote
## diff lwr upr p adj
## 2-1 4.6 -1.675865 10.875865 0.1753581
## 3-1 3.8 -2.475865 10.075865 0.3087315
## 4-1 3.4 -2.875865 9.675865 0.4006635
## 5-1 4.2 -2.075865 10.475865 0.2340133
## 3-2 -0.8 -7.075865 5.475865 0.9906712
## 4-2 -1.2 -7.475865 5.075865 0.9596565
## 5-2 -0.4 -6.675865 5.875865 0.9993563
## 4-3 -0.4 -6.675865 5.875865 0.9993563
## 5-3 0.4 -5.875865 6.675865 0.9993563
## 5-4 0.8 -5.475865 7.075865 0.9906712
##
## $Operador
## diff lwr upr p adj
## 2-1 7.2 0.924135 13.475865 0.0252496
## 3-1 2.8 -3.475865 9.075865 0.5669444
## 4-1 4.6 -1.675865 10.875865 0.1753581
## 5-1 5.4 -0.875865 11.675865 0.0966410
## 3-2 -4.4 -10.675865 1.875865 0.2028079
## 4-2 -2.6 -8.875865 3.675865 0.6270350
## 5-2 -1.8 -8.075865 4.475865 0.8524755
## 4-3 1.8 -4.475865 8.075865 0.8524755
## 5-3 2.6 -3.675865 8.875865 0.6270350
## 5-4 0.8 -5.475865 7.075865 0.9906712
##
## $Montajes
## diff lwr upr p adj
## c-b -3.2 -9.475865 3.075865 0.4528137
## d-b -2.6 -8.875865 3.675865 0.6270350
## e-b -2.8 -9.075865 3.475865 0.5669444
## f-b 0.6 -5.675865 6.875865 0.9968743
## d-c 0.6 -5.675865 6.875865 0.9968743
## e-c 0.4 -5.875865 6.675865 0.9993563
## f-c 3.8 -2.475865 10.075865 0.3087315
## e-d -0.2 -6.475865 6.075865 0.9999587
## f-d 3.2 -3.075865 9.475865 0.4528137
## f-e 3.4 -2.875865 9.675865 0.4006635Analizando las comparaciones de las diferentes formulaciones, se evidencia que la formulación B es la de menos tiempo de respuesta entre las demás, por esta razón, si lo que se desea es encontrar la formulación más rapida, esta sería la mejor decisión.
Obtenemos una gráfica de estos resultados
par(mfrow=c(2,2))
plot(TukeyHSD(anovaGreco))
Pruebas de comparación de medias
library(agricolae)- Prueba LSD
LSD.test(anovaGreco,"Formulacion",console=TRUE)##
## Study: anovaGreco ~ "Formulacion"
##
## LSD t Test for Rapidez
##
## Mean Square Error: 8.25
##
## Formulacion, means and individual ( 95 %) CI
##
## Rapidez std r LCL UCL Min Max
## A 3.6 4.669047 5 0.6378841 6.5621159 -1 11
## B -4.8 2.167948 5 -7.7621159 -1.8378841 -8 -2
## C -2.6 4.393177 5 -5.5621159 0.3621159 -7 4
## D 4.8 5.403702 5 1.8378841 7.7621159 -1 13
## E 1.0 3.391165 5 -1.9621159 3.9621159 -3 6
##
## Alpha: 0.05 ; DF Error: 8
## Critical Value of t: 2.306004
##
## least Significant Difference: 4.189065
##
## Treatments with the same letter are not significantly different.
##
## Rapidez groups
## D 4.8 a
## A 3.6 a
## E 1.0 ab
## C -2.6 bc
## B -4.8 c- Prueba de Tukey
HSD.test(anovaGreco, "Formulacion",console=TRUE)##
## Study: anovaGreco ~ "Formulacion"
##
## HSD Test for Rapidez
##
## Mean Square Error: 8.25
##
## Formulacion, means
##
## Rapidez std r Min Max
## A 3.6 4.669047 5 -1 11
## B -4.8 2.167948 5 -8 -2
## C -2.6 4.393177 5 -7 4
## D 4.8 5.403702 5 -1 13
## E 1.0 3.391165 5 -3 6
##
## Alpha: 0.05 ; DF Error: 8
## Critical Value of Studentized Range: 4.885754
##
## Minimun Significant Difference: 6.275865
##
## Treatments with the same letter are not significantly different.
##
## Rapidez groups
## D 4.8 a
## A 3.6 ab
## E 1.0 abc
## C -2.6 bc
## B -4.8 c- Prueba de Student-Newman-Keuls (SNK)
SNK.test(anovaGreco, "Formulacion",console=TRUE)##
## Study: anovaGreco ~ "Formulacion"
##
## Student Newman Keuls Test
## for Rapidez
##
## Mean Square Error: 8.25
##
## Formulacion, means
##
## Rapidez std r Min Max
## A 3.6 4.669047 5 -1 11
## B -4.8 2.167948 5 -8 -2
## C -2.6 4.393177 5 -7 4
## D 4.8 5.403702 5 -1 13
## E 1.0 3.391165 5 -3 6
##
## Alpha: 0.05 ; DF Error: 8
##
## Critical Range
## 2 3 4 5
## 4.189065 5.190805 5.817361 6.275865
##
## Means with the same letter are not significantly different.
##
## Rapidez groups
## D 4.8 a
## A 3.6 a
## E 1.0 ab
## C -2.6 bc
## B -4.8 c- Prueba de Scheffé
scheffe.test(anovaGreco, "Formulacion",console=TRUE)##
## Study: anovaGreco ~ "Formulacion"
##
## Scheffe Test for Rapidez
##
## Mean Square Error : 8.25
##
## Formulacion, means
##
## Rapidez std r Min Max
## A 3.6 4.669047 5 -1 11
## B -4.8 2.167948 5 -8 -2
## C -2.6 4.393177 5 -7 4
## D 4.8 5.403702 5 -1 13
## E 1.0 3.391165 5 -3 6
##
## Alpha: 0.05 ; DF Error: 8
## Critical Value of F: 3.837853
##
## Minimum Significant Difference: 7.11756
##
## Means with the same letter are not significantly different.
##
## Rapidez groups
## D 4.8 a
## A 3.6 ab
## E 1.0 abc
## C -2.6 bc
## B -4.8 c- Prueba de Duncan
duncan.test(anovaGreco, "Formulacion",console=TRUE)##
## Study: anovaGreco ~ "Formulacion"
##
## Duncan's new multiple range test
## for Rapidez
##
## Mean Square Error: 8.25
##
## Formulacion, means
##
## Rapidez std r Min Max
## A 3.6 4.669047 5 -1 11
## B -4.8 2.167948 5 -8 -2
## C -2.6 4.393177 5 -7 4
## D 4.8 5.403702 5 -1 13
## E 1.0 3.391165 5 -3 6
##
## Alpha: 0.05 ; DF Error: 8
##
## Critical Range
## 2 3 4 5
## 4.189065 4.365401 4.463963 4.523055
##
## Means with the same letter are not significantly different.
##
## Rapidez groups
## D 4.8 a
## A 3.6 a
## E 1.0 ab
## C -2.6 bc
## B -4.8 c- Prueba de Bonferroni
LSD.test(anovaGreco, "Formulacion", p.adj= "bon",console=TRUE)##
## Study: anovaGreco ~ "Formulacion"
##
## LSD t Test for Rapidez
## P value adjustment method: bonferroni
##
## Mean Square Error: 8.25
##
## Formulacion, means and individual ( 95 %) CI
##
## Rapidez std r LCL UCL Min Max
## A 3.6 4.669047 5 0.6378841 6.5621159 -1 11
## B -4.8 2.167948 5 -7.7621159 -1.8378841 -8 -2
## C -2.6 4.393177 5 -5.5621159 0.3621159 -7 4
## D 4.8 5.403702 5 1.8378841 7.7621159 -1 13
## E 1.0 3.391165 5 -1.9621159 3.9621159 -3 6
##
## Alpha: 0.05 ; DF Error: 8
## Critical Value of t: 3.832519
##
## Minimum Significant Difference: 6.962116
##
## Treatments with the same letter are not significantly different.
##
## Rapidez groups
## D 4.8 a
## A 3.6 ab
## E 1.0 abc
## C -2.6 bc
## B -4.8 cCOMPROBACIÓN DE SUPUESTOS DEL MODELO
Normalidad de los residuos
Con un valor P de 0,87 NO se rechaza la hipotesis de la normalidad de los residuos:
shapiro.test(anovaGreco$res) ##
## Shapiro-Wilk normality test
##
## data: anovaGreco$res
## W = 0.97967, p-value = 0.8786Se evidencia media= 0, lo que supone la normalidad de los datos:
summary(anovaGreco$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.2 -1.2 0.0 0.0 1.0 3.4Graficamente se evidencia una distribución Normal:
boxplot(anovaGreco$residuals,col="lightblue") 
Se evidencia una distribución normal mediante histograma y su linea de densidad
hist(anovaGreco$residuals, col="lightgreen", main = "Histograma de los Residuos",
freq = F, xlab="Residuos",ylab="Densidad")
lines(density(anovaGreco$residuals), col="red", lwd=3)
En la Prueba de Normalidad con Q-Q Plot, se evidencia que los datos se ajustan a la linea de normalidad
qqnorm(anovaGreco$residuals)
qqline(anovaGreco$residuals)
Supuesto de Independencia
plot(anovaGreco$residuals,main = "Residuos vs Observaciones") 
boxplot(anovaGreco$residuals~Formulacion, xlab="Formulacion",ylab="Residuos",
col = c("yellow", "blue", "white","green", "red")) # Homocedasticidad
Gráfico de predichos contra residuos estandarizados
pred=fitted(anovaGreco)
resid=rstandard(anovaGreco)
plot(pred,resid,xlab="Valores predichos", ylab="Residuos estandarizados",abline(h=0))
Homogenidad de Varianzas
Se maneja la hipotesis nula de la homogenidad de varianza
De acuerdo al valor P en cada prueba, No existe evidencia estadística para rechazar que existe homogenidad de varianzas en todos los factores y el tratamiento
bartlett.test(anovaGreco$residuals ~ Formulacion)##
## Bartlett test of homogeneity of variances
##
## data: anovaGreco$residuals by Formulacion
## Bartlett's K-squared = 0.31314, df = 4, p-value = 0.989bartlett.test(anovaGreco$residuals ~ Montajes)##
## Bartlett test of homogeneity of variances
##
## data: anovaGreco$residuals by Montajes
## Bartlett's K-squared = 7.4545, df = 4, p-value = 0.1137bartlett.test(anovaGreco$residuals ~ Lote)##
## Bartlett test of homogeneity of variances
##
## data: anovaGreco$residuals by Lote
## Bartlett's K-squared = 5.7344, df = 4, p-value = 0.2199bartlett.test(anovaGreco$residuals ~ Operador)##
## Bartlett test of homogeneity of variances
##
## data: anovaGreco$residuals by Operador
## Bartlett's K-squared = 0.33032, df = 4, p-value = 0.9878En la Prueba Levene para Homogenidad de varianzas no se rechaza la hipotesis nula que indica homogeneidad de varianza en el tratamiento
# install.packages("car") quien no tenga instalado este paquete lo debe hacer
library(car)## Loading required package: carDataleveneTest(Rapidez~Formulacion, dataGreco, center = "median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 4 0.5858 0.6766
## 20¡¡¡Gracias!!!