Se midió el efecto de la cantidad de harina, azúcar y vainilla sobre la esponjosidad de mantecadas (índice de crecimiento) en un diseño experimental factorial 2 a la 3 o 2x2x2. La unidad experimental fue una mantecada con dos repeticiones por tratamiento. El objetivo del experimento fue determinar la mejor combinación de los factores, para lograr la mantecada más esponjosa, los niveles de los factores y las combinaciones de los niveles se muestran a continuación, en total se obtuvieron ocho tratamientos.

# Importación de los datos
DATOS <- read.csv("C:/Users/Zuleica Trujano/Documents/Documentos Académicos/MC IA/1er Cuatrimestre/Diseños experimentales/Datos y Análisis/Factorial 2^3/Datos mantecadas.csv")
DATOS
##    ï..Harina Azucar Vainilla Crecimiento
## 1         85     75        3         2.8
## 2         95     75        3         2.9
## 3         85    100        3         2.4
## 4         95    100        3         2.6
## 5         85     75        3         2.5
## 6         95     75        3         2.3
## 7         85    100        3         2.7
## 8         95    100        3         3.0
## 9         85     75        5         3.6
## 10        95     75        5         3.0
## 11        85    100        5         4.0
## 12        95    100        5         4.2
## 13        85     75        5         3.0
## 14        95     75        5         2.8
## 15        85    100        5         4.2
## 16        95    100        5         3.7

Para establecer las variables factor y la respuesta y así trabajar adecuadamente con ellas, las variables harina, azucar y vainilla se establecen como factores usando la función as.factor(). El índice de crecimiento de las mantecadas es la variable respuesta.

# Ver las variables y fijarlas
attach(DATOS)
names(DATOS)
## [1] "ï..Harina"   "Azucar"      "Vainilla"    "Crecimiento"
str(DATOS)
## 'data.frame':    16 obs. of  4 variables:
##  $ ï..Harina  : int  85 95 85 95 85 95 85 95 85 95 ...
##  $ Azucar     : int  75 75 100 100 75 75 100 100 75 75 ...
##  $ Vainilla   : int  3 3 3 3 3 3 3 3 5 5 ...
##  $ Crecimiento: num  2.8 2.9 2.4 2.6 2.5 2.3 2.7 3 3.6 3 ...
# Cambio de variables a objetos
harina <- as.factor(ï..Harina)
azucar <- as.factor(Azucar)
vainilla <- as.factor(Vainilla)
crecimiento <- DATOS$Crecimiento

ANOVA

Se realiza el análisis de varianza (ANOVA) con la librería agricolae, dandole los elementos al modelo, usando un asterisco entre los factores para indicar que se desea calcular el efecto de la interacción. La única interacción que resultó ser significativa fue azúcar con vainilla.

library(agricolae)
modelo <- lm(crecimiento~harina*azucar*vainilla)
anova <- aov(modelo)
summary(anova)
##                        Df Sum Sq Mean Sq F value   Pr(>F)    
## harina                  1  0.031   0.031   0.353 0.569095    
## azucar                  1  0.951   0.951  10.942 0.010731 *  
## vainilla                1  3.331   3.331  38.338 0.000262 ***
## harina:azucar           1  0.076   0.076   0.871 0.378115    
## harina:vainilla         1  0.141   0.141   1.619 0.239005    
## azucar:vainilla         1  0.766   0.766   8.813 0.017908 *  
## harina:azucar:vainilla  1  0.001   0.001   0.007 0.934490    
## Residuals               8  0.695   0.087                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Para ver las interacciones dobles graficamente se usan los siguientes comandos, indicando en cada uno los factores a comparar. Graficamente se observan interacciones entre los pares de factores, aunque estas no resultaron significativas en el ANOVA, con excepción de azúcar*vainilla.

# Gráficas de interacciones
interaction.plot(azucar, harina, crecimiento, fixed=T, xlab="Gramos de azúcar", ylab="Crecimiento de las mantecadas", type = c("b"), trace.label="Gramos de harina")

interaction.plot(vainilla, harina, crecimiento, fixed=T, xlab="Mílitros de esencia de vainilla", ylab="Crecimiento de las mantecadas", type = c("b"), trace.label="Gramos de harina")

interaction.plot(vainilla, azucar, crecimiento, fixed=T, xlab="Mililítros de esencia de vainilla", ylab="Crecimiento de las mantecadas", type = c("b"), trace.label="Gramos de azúcar")

Comparación de medias: Tukey

Para saber cual efecto factorial de la interacción azúcarvainilla y de la intracción triple (azúcarvainilla*harina) fue el que produjo mayor índice de crecimiento, se realiza la prueba de la comparación de medias. Recordando el objetivo del estudio que fue determinar la mejor combinación de factores que produjeran el mayor crecimiento en el horneo de las mantecadas. De acuerdo con los resultados de comparación de medias con interacción triple, no hay diferencia significativa entre los primeros tres tratamientos. Es decir, se podría usar cualquiera de estos en la elaboración de las mantecadas; sin embargo, se observa que a cualquier nivel de harina la combinación con 100 gramos de azúcar y 5 mL de vainilla produce el mayor crecimiento, en este caso se recomendaría usar el tratamiento con 85 g de harina, dado que con menor insumo de harina se produce una mayor esponjosidad.

# Tratamientos de la interacción azucar*vainilla: comparación de medias de Tukey
azu_vain <- as.factor(azucar:vainilla)# un solo vector para interacción
azu_vain
##  [1] 75:3  75:3  100:3 100:3 75:3  75:3  100:3 100:3 75:5  75:5  100:5 100:5
## [13] 75:5  75:5  100:5 100:5
## Levels: 75:3 75:5 100:3 100:5
anova_av <- aov(crecimiento~azu_vain)
anova_av
## Call:
##    aov(formula = crecimiento ~ azu_vain)
## 
## Terms:
##                 azu_vain Residuals
## Sum of Squares  5.046875  0.942500
## Deg. of Freedom        3        12
## 
## Residual standard error: 0.2802529
## Estimated effects may be unbalanced
mantecadas_av <- TukeyHSD(anova_av)
mantecadas_av
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = crecimiento ~ azu_vain)
## 
## $azu_vain
##               diff        lwr       upr     p adj
## 75:5-75:3    0.475 -0.1133433 1.0633433 0.1306332
## 100:3-75:3   0.050 -0.5383433 0.6383433 0.9940578
## 100:5-75:3   1.400  0.8116567 1.9883433 0.0000672
## 100:3-75:5  -0.425 -1.0133433 0.1633433 0.1942493
## 100:5-75:5   0.925  0.3366567 1.5133433 0.0026377
## 100:5-100:3  1.350  0.7616567 1.9383433 0.0000956
# Tratamientos de la interacción triple: comparación de medias de Tukey
ctrat <- as.factor(harina:azucar:vainilla) # un solo vector para interacción
ctrat
##  [1] 85:75:3  95:75:3  85:100:3 95:100:3 85:75:3  95:75:3  85:100:3 95:100:3
##  [9] 85:75:5  95:75:5  85:100:5 95:100:5 85:75:5  95:75:5  85:100:5 95:100:5
## 8 Levels: 85:75:3 85:75:5 85:100:3 85:100:5 95:75:3 95:75:5 ... 95:100:5
anova_comp <- aov(crecimiento~ctrat) # modelo ANOVA
mantecadas_hav <- TukeyHSD(anova_comp)
mantecadas_hav
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = crecimiento ~ ctrat)
## 
## $ctrat
##                    diff         lwr         upr     p adj
## 85:75:5-85:75:3    0.65 -0.51633715  1.81633715 0.4327920
## 85:100:3-85:75:3  -0.10 -1.26633715  1.06633715 0.9999444
## 85:100:5-85:75:3   1.45  0.28366285  2.61633715 0.0153120
## 95:75:3-85:75:3   -0.05 -1.21633715  1.11633715 0.9999995
## 95:75:5-85:75:3    0.25 -0.91633715  1.41633715 0.9836621
## 95:100:3-85:75:3   0.15 -1.01633715  1.31633715 0.9992110
## 95:100:5-85:75:3   1.30  0.13366285  2.46633715 0.0283462
## 85:100:3-85:75:5  -0.75 -1.91633715  0.41633715 0.2958262
## 85:100:5-85:75:5   0.80 -0.36633715  1.96633715 0.2414757
## 95:75:3-85:75:5   -0.70 -1.86633715  0.46633715 0.3596161
## 95:75:5-85:75:5   -0.40 -1.56633715  0.76633715 0.8537371
## 95:100:3-85:75:5  -0.50 -1.66633715  0.66633715 0.6912179
## 95:100:5-85:75:5   0.65 -0.51633715  1.81633715 0.4327920
## 85:100:5-85:100:3  1.55  0.38366285  2.71633715 0.0103035
## 95:75:3-85:100:3   0.05 -1.11633715  1.21633715 0.9999995
## 95:75:5-85:100:3   0.35 -0.81633715  1.51633715 0.9151376
## 95:100:3-85:100:3  0.25 -0.91633715  1.41633715 0.9836621
## 95:100:5-85:100:3  1.40  0.23366285  2.56633715 0.0187491
## 95:75:3-85:100:5  -1.50 -2.66633715 -0.33366285 0.0125417
## 95:75:5-85:100:5  -1.20 -2.36633715 -0.03366285 0.0432810
## 95:100:3-85:100:5 -1.30 -2.46633715 -0.13366285 0.0283462
## 95:100:5-85:100:5 -0.15 -1.31633715  1.01633715 0.9992110
## 95:75:5-95:75:3    0.30 -0.86633715  1.46633715 0.9582777
## 95:100:3-95:75:3   0.20 -0.96633715  1.36633715 0.9953832
## 95:100:5-95:75:3   1.35  0.18366285  2.51633715 0.0230227
## 95:100:3-95:75:5  -0.10 -1.26633715  1.06633715 0.9999444
## 95:100:5-95:75:5   1.05 -0.11633715  2.21633715 0.0827233
## 95:100:5-95:100:3  1.15 -0.01633715  2.31633715 0.0536421

Verificación de supuestos

Normalidad

Primero se obtiene el conjunto de residuos obtenidos en el ANOVA

# Extraer residuales del ANOVA
residuales <- (anova$res)
residuales
##     1     2     3     4     5     6     7     8     9    10    11    12    13 
##  0.15  0.30 -0.15 -0.20 -0.15 -0.30  0.15  0.20  0.30  0.10 -0.10  0.25 -0.30 
##    14    15    16 
## -0.10  0.10 -0.25

Después se realiza la prueba para Normalidad de Shapiro-Wilk para pocas muestras y se grafican. En la gráfica se espera observar que los puntos se encuentren cercanos a la diagonal. El p-value arrojado por la prueba es mayor que 0.05, por lo que no se rechaza la hipótesis nula para la normalidad de los residuales, es decir, estos siguen una distribución normal.

# NORMALIDAD
# Prueba de Shapiro-Wilk
shapiro.test(residuales) 
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.90657, p-value = 0.1025
# Gráfica Q-Q
plot(anova,2)

Homocedasticidad

Se hace la prueba de Bartlett para cada factor y la interacción. En esta prueba los valores están por encima del valor de p 0.05, por lo que existe homogeneidad de varianza u homocedasticidad de los residuales. Se confirma otro supuesto. En la gráfica se busca observar que los puntos se distribuyan aleatoriamente y sin ningun patrón específico, igualmente se espera que no existan valores atípicos o extremos.

# HOMOCEDASTICIDAD 
# Prueba de Bartlett con un conjunto de residuos
bartlett.test(residuales, harina) # Factor A: harina
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuales and harina
## Bartlett's K-squared = 0.18137, df = 1, p-value = 0.6702
bartlett.test(residuales, azucar) # Factor B: azucar
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuales and azucar
## Bartlett's K-squared = 0.33332, df = 1, p-value = 0.5637
bartlett.test(residuales, vainilla) # Factor C: vainilla
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuales and vainilla
## Bartlett's K-squared = 0.00033816, df = 1, p-value = 0.9853
bartlett.test(residuales, harina:azucar:vainilla) # Interacción ABC
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuales and harina:azucar:vainilla
## Bartlett's K-squared = 1.8289, df = 7, p-value = 0.9687
# Gráfica de residuos contra predichos
plot(anova, 1)

Independencia de los residuales

El valor de p es mayor que 0.05, se rechaza hipótesis nula los residuales son independientes,

# INDEPENDENCIA
# Test de independencia Durbin Watson
library("lmtest")
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
dwtest(anova)
## 
##  Durbin-Watson test
## 
## data:  anova
## DW = 1.6906, p-value = 0.3956
## alternative hypothesis: true autocorrelation is greater than 0

Media cero

No se reachaza hipótesis nula para la media cero de los residuales, estos tienen media cero.

# MEDIA CERO
# Obtener la media de los residuales
media_cero <- sample(0:0, 64, replace=TRUE)
# Hacer prueba T de Student
mean_test <- t.test(residuales, media_cero)
mean_test
## 
##  Welch Two Sample t-test
## 
## data:  residuales and media_cero
## t = -1.9342e-16, df = 15, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1146996  0.1146996
## sample estimates:
##     mean of x     mean of y 
## -1.040834e-17  0.000000e+00

Se cumplieron todos los supuestos de normalidad de los residuales, el uso del ANOVA fue pertinente para el análisis de los datos.