Diseño de Experimentos 202560
Profe: Julio Hurtado Marquez
31/01/2023
Ejemplos de Diseño de Experimentos con Un Factor - 2025
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
Tematica 1: Diseños en Cuadrados Latinos
1. Problema Ilustrativo: el efecto de cinco diferentes catalizadores (A, B, C, D y E)
Problema: Se quiere estudiar el efecto de cinco diferentes catalizadores (A, B, C, D y 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 sólo se pueden realizar cinco corridas diarias. El experimentador decide correr los experimentos con un diseño en cuadro latino para controlar activamente los lotes y días. Los datos obtenidos son:
| Lote / Día | 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 = 2 | B = 3 | A = 8 | C = 8 |
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
2. Las Hipotesis y el modelo
A.las hipotesis del problema son las siguientes:
H0: Los catalizadores no tienen efecto sobre el tiempo de reacción.
H1: Al menos un catalizador tiene efecto sobre el tiempo de reacción.
las hipotesis del problema con los catalizadores en simbolos son las siguientes:
H0: \(\mu_A = \mu_B = \mu_C = \mu_D = \mu_E\)
H1: Al menos un catalizador tiene efecto sobre el tiempo de reacción.
B.las hipotesis del problema con los lotes son las siguientes:
H0: Los lotes no tienen efecto sobre el tiempo de reacción.
H1: Al menos un lote tiene efecto sobre el tiempo de reacción.
C. las hipotesis del problema con los dias son las siguientes:
H0: Los días no tienen efecto sobre el tiempo de reacción.
H1: Al menos un día tiene efecto sobre el tiempo de reacción.
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
3. Los datos y el DataFrame
# Datos
Lote = factor(rep(1:5, each=5))
Día = factor(rep(1:5, times=5))
Catalizador = factor(c("A","B","D","C","E",
"C","E","A","D","B",
"B","A","C","E","D",
"D","C","E","B","A",
"E","D","B","A","C"))
Tiempo = c(8,7,1,7,3,
11,2,7,3,8,
4,9,10,1,5,
6,8,6,6,10,
4,2,3,8,8)
datos = data.frame(Lote, Día, Catalizador, Tiempo)
datos🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
4. Análisis de varianza
# Análisis de varianza
modelo = aov(Tiempo ~ Lote + Día + Catalizador, data=datos)
summary(modelo)## Df Sum Sq Mean Sq F value Pr(>F)
## Lote 4 15.44 3.86 1.235 0.347618
## Día 4 12.24 3.06 0.979 0.455014
## Catalizador 4 141.44 35.36 11.309 0.000488 ***
## Residuals 12 37.52 3.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14. Interpretacion del Anovas
El valor p del catalizador es 0.0005, por lo que se rechaza la hipótesis nula y se concluye que al menos un catalizador tiene efecto sobre el tiempo de reacción.
El valor p del lote es 0.3476, por lo que no se rechaza la hipótesis nula y se concluye que no evidencia en la muestra de que los lotes tengan efecto sobre el tiempo de reacción.
El valor p del lote es 0.455, por lo que no se rechaza la hipótesis nula y se concluye que no evidencia en la muestra de el Día tengan efecto sobre el tiempo de reacción.
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
5. Boxplot para los Factores
A. Boxplot para Catalizador
# Boxplot para Catalizador con colores
boxplot(Tiempo ~ Catalizador, data=datos, main="Boxplot de Tiempo por Catalizador", xlab="Catalizador", ylab="Tiempo", col=c("red2", "blue2", "yellow2", "orange2","green2"))
B. Boxplot para Lote
# Boxplot para Catalizador con colores
boxplot(Tiempo ~ Lote, data=datos, main="Boxplot de Tiempo por Catalizador", xlab="Lote", ylab="Tiempo", col=c("red2", "blue2", "yellow2", "orange2","green2"))
C. Boxplot para Día
# Boxplot para Catalizador con colores
boxplot(Tiempo ~ Día, data=datos, main="Boxplot de Tiempo por Catalizador", xlab="Día", ylab="Tiempo", col=c("red2", "blue2", "yellow2", "orange2","green2"))
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
6. Prueba de Rangos Multiples - HSD
A. Prueba de Rangos Multiples para Catalizador
TukeyHSD(modelo, "Catalizador")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Tiempo ~ Lote + Día + Catalizador, data = datos)
##
## $Catalizador
## diff lwr upr p adj
## B-A -2.8 -6.3646078 0.7646078 0.1539433
## C-A 0.4 -3.1646078 3.9646078 0.9960012
## D-A -5.0 -8.5646078 -1.4353922 0.0055862
## E-A -5.2 -8.7646078 -1.6353922 0.0041431
## C-B 3.2 -0.3646078 6.7646078 0.0864353
## D-B -2.2 -5.7646078 1.3646078 0.3365811
## E-B -2.4 -5.9646078 1.1646078 0.2631551
## D-C -5.4 -8.9646078 -1.8353922 0.0030822
## E-C -5.6 -9.1646078 -2.0353922 0.0023007
## E-D -0.2 -3.7646078 3.3646078 0.9997349B. Prueba de Rangos Multiples para Lote
TukeyHSD(modelo, "Lote")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Tiempo ~ Lote + Día + Catalizador, data = datos)
##
## $Lote
## diff lwr upr p adj
## 2-1 1.0 -2.564608 4.564608 0.8936609
## 3-1 0.6 -2.964608 4.164608 0.9816047
## 4-1 2.0 -1.564608 5.564608 0.4225127
## 5-1 -0.2 -3.764608 3.364608 0.9997349
## 3-2 -0.4 -3.964608 3.164608 0.9960012
## 4-2 1.0 -2.564608 4.564608 0.8936609
## 5-2 -1.2 -4.764608 2.364608 0.8166339
## 4-3 1.4 -2.164608 4.964608 0.7232162
## 5-3 -0.8 -4.364608 2.764608 0.9489243
## 5-4 -2.2 -5.764608 1.364608 0.3365811C. Prueba de Rangos Multiples para Día
TukeyHSD(modelo, "Día")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Tiempo ~ Lote + Día + Catalizador, data = datos)
##
## $Día
## diff lwr upr p adj
## 2-1 -1.0 -4.564608 2.564608 0.8936609
## 3-1 -1.2 -4.764608 2.364608 0.8166339
## 4-1 -1.6 -5.164608 1.964608 0.6212723
## 5-1 0.2 -3.364608 3.764608 0.9997349
## 3-2 -0.2 -3.764608 3.364608 0.9997349
## 4-2 -0.6 -4.164608 2.964608 0.9816047
## 5-2 1.2 -2.364608 4.764608 0.8166339
## 4-3 -0.4 -3.964608 3.164608 0.9960012
## 5-3 1.4 -2.164608 4.964608 0.7232162
## 5-4 1.8 -1.764608 5.364608 0.5188508🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
7. Verificación de supuestos del ANOVA
A. Grafica para supuestos
par( mfrow = c(2,2) )
plot(modelo, which=5)
plot(modelo, which=1)
plot(modelo, which=2)
plot(residuals(modelo) ~ Tiempo , main="Residuals vs corridas", font.main=1, xlab="Tiempo", ylab="Residuals")
abline(h = 0, lty = 2)
B. Normalidad - Shapiro Wilks
# Prueba de normalidad de Shapiro-Wilk
shapiro.test(modelo$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.96606, p-value = 0.5476C. Homocedasticidad - Bartlett
1. Catalizador - Bartlett
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
bartlett.test(modelo$residuals ~ Catalizador, data=datos)##
## Bartlett test of homogeneity of variances
##
## data: modelo$residuals by Catalizador
## Bartlett's K-squared = 2.6529, df = 4, p-value = 0.61752. Lote - Bartlett
\[H_0:\sigma^2_{L_1}=\sigma^2_{L_2}=\sigma^2_{L_3}=\sigma^2_{L_3}=\sigma^2_{L_5}=\sigma^2\]
bartlett.test(modelo$residuals ~ Lote, data=datos)##
## Bartlett test of homogeneity of variances
##
## data: modelo$residuals by Lote
## Bartlett's K-squared = 0.69132, df = 4, p-value = 0.95243. Día - Bartlett
\[H_0:\sigma^2_{D_1}=\sigma^2_{D_2}=\sigma^2_{D_3}=\sigma^2_{D_3}=\sigma^2_{D_5}=\sigma^2\]
bartlett.test(modelo$residuals ~ Día, data=datos)##
## Bartlett test of homogeneity of variances
##
## data: modelo$residuals by Día
## Bartlett's K-squared = 2.6163, df = 4, p-value = 0.6239D. Homocedasticidad - Levene
1. Catalizador - Levene
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
#Levene
library(car)## Loading required package: carDataleveneTest(modelo$residuals ~ Catalizador, data=datos)2. Lote - Levene
\[H_0:\sigma^2_{L_1}=\sigma^2_{L_2}=\sigma^2_{L_3}=\sigma^2_{L_3}=\sigma^2_{L_5}=\sigma^2\]
#Levene
library(car)
leveneTest(modelo$residuals ~ Lote, data=datos)3. Día - Levene
\[H_0:\sigma^2_{D_1}=\sigma^2_{D_2}=\sigma^2_{D_3}=\sigma^2_{D_3}=\sigma^2_{D_5}=\sigma^2\]
#Levene
library(car)
leveneTest(modelo$residuals ~ Día, data=datos)🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
Tematica 2: Diseños en Cuadrados Greco Latinos
1. Problema Ilustrativo: el efecto de cinco diferentes Formulaciones (A, B, C, D y E)
Este experimento se diseñó para investigar el efecto de cinco Formulaciones diferentes (A, B, C, D, E) de carga propulsora sobre la rapidez de combustión en sistemas de expulsión de tripulación de aviones. El objetivo es aislar el efecto de cada formulación minimizando la influencia de factores perturbadores.
Factores y Niveles
- Factores Controlados:
- Factor de Interés: Formulaciones (A, B, C, D, E)
- Factores de Bloqueo (Factores Perturbadores):
- Lotes de Materia Prima (1, 2, 3, 4, 5)
- Operadores (1, 2, 3, 4, 5)
- Montajes de Prueba (α, β, γ, δ, ε)
La estructura del diseño se presenta en la siguiente tabla:
Tabla 6. Diseño en cuadro grecolatino
| (Operador) | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| 1 | Aα=24 | Bγ=20 | Cε=19 | Dβ=24 | Eδ=24 |
| 2 | Bβ=17 | Cδ=24 | Dα=30 | Eγ=27 | Aε=36 |
| 3 | Cγ=18 | Dε=38 | Eβ=26 | Aδ=27 | Bα=21 |
| 4 | Dδ=26 | Eα=31 | Aγ=26 | Bε=23 | Cβ=22 |
| 5 | Eε=22 | Aβ=30 | Bδ=20 | Cα=29 | Dγ=31 |
Interpretación de la Tabla:
- Cada celda de la tabla representa una combinación única de Operador, Lote, Formulación y Montaje de Prueba, junto con el resultado observado (rapidez de combustión). Por ejemplo, la celda en la fila 1 y columna 1 indica que el Operador 1 utilizó el Lote 1 con la Formulación A y el Montaje de Prueba α, obteniendo una rapidez de combustión de 24.
Análisis
- El diseño de Cuadrado Grecolatino permite realizar un análisis de varianza (ANOVA) para determinar si existen diferencias estadísticamente significativas en la rapidez de combustión entre las diferentes formulaciones, controlando al mismo tiempo el efecto de los lotes, operadores y montajes de prueba.
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
2. Las Hipotesis y el modelo
A.las hipotesis del problema son las siguientes:
H0: Las formulaciones no tienen efecto sobre la rápidez de combustion.
H1: Al menos una formulación tiene efecto sobre la rápidez de combustion.
las hipotesis del problema con las formulaciones en simbolos son las siguientes:
H0: \(\mu_A = \mu_B = \mu_C = \mu_D = \mu_E\)
H1: Al menos una formulación tiene efecto sobre la rápidez de combustion.
B.las hipotesis del problema con los lotes son las siguientes:
H0: Los lotes no tienen efecto sobre la rápidez de combustion.
H1: Al menos un lote tiene efecto la rápidez de combustion.
C. las hipotesis del problema con los operadores son las siguientes:
H0: Los operadores no tienen efecto sobre la rápidez de combustion..
H1: Al menos un operador tiene efecto sobre la rápidez de combustion.
D. las hipotesis del problema con los Montajes son las siguientes:
H0: Los Montajes no tienen efecto sobre la rápidez de combustion..
H1: Al menos un Montajes tiene efecto sobre la rápidez de combustion.
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
3. Los datos y el DataFrame
# Datos rapidez de combustion
Lote = factor(rep(1:5, each=5))
Operador = factor(rep(1:5, times=5))
Montaje = factor(c("α","γ","ε","β","δ",
"β","δ","α","γ","ε",
"γ","ε","β","δ","α",
"δ","α","γ","ε","β",
"ε","β","δ","α","γ"))
Formulación = 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(24,20,19,24,24,
17,24,30,27,36,
18,38,26,27,21,
26,31,26,23,22,
22,30,20,29,31)
datos2 = data.frame(Lote, Operador, Montaje, Formulación, Rapidez)
datos2🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
4. Análisis de varianza
# Análisis de varianza
modelo2 = aov(Rapidez ~ Formulación + Montaje + Operador + Lote, data=datos2)
summary(modelo2)## Df Sum Sq Mean Sq F value Pr(>F)
## Formulación 4 330 82.50 10.000 0.00334 **
## Montaje 4 62 15.50 1.879 0.20764
## Operador 4 150 37.50 4.545 0.03293 *
## Lote 4 68 17.00 2.061 0.17831
## Residuals 8 66 8.25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14. Interpretacion del Anovas
El valor p del catalizador es 0.0005, por lo que se rechaza la hipótesis nula y se concluye que al menos un catalizador tiene efecto sobre el tiempo de reacción.
El valor p del lote es 0.3476, por lo que no se rechaza la hipótesis nula y se concluye que no evidencia en la muestra de que los lotes tengan efecto sobre el tiempo de reacción.
El valor p del lote es 0.455, por lo que no se rechaza la hipótesis nula y se concluye que no evidencia en la muestra de el Día tengan efecto sobre el tiempo de reacción.
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
5. Boxplot para los Factores
A. Boxplot para Catalizador
# Boxplot para Catalizador con colores
boxplot(Rapidez ~ Formulación, data=datos2, main="Boxplot de Rapidez por Formulación", xlab="Formulación", ylab="Rapidez", col=c("red2", "blue2", "yellow2", "orange2","green2"))
B. Boxplot para Lote
# Boxplot para Catalizador con colores
boxplot(Rapidez ~ Montaje, data=datos2, main="Boxplot de Tiempo por Montaje", xlab="Montaje", ylab="Rapidez", col=c("red2", "blue2", "yellow2", "orange2","green2"))
C. Boxplot para Operador
Rapidez ~ Formulación + Montaje + Operador + Lote
# Boxplot para Catalizador con colores
boxplot(Rapidez ~ Operador, data=datos2, main="Boxplot de Rapidez por Operador", xlab="Rapidez", ylab="Rapidez", col=c("red2", "blue2", "yellow2", "orange2","green2"))
C. Boxplot para Lote
# Boxplot para Catalizador con colores
boxplot(Rapidez ~ Lote, data=datos2, main="Boxplot de Rapidez por Lote", xlab="Lote", ylab="Rapidez", col=c("red2", "blue2", "yellow2", "orange2","green2"))
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
6. Prueba de Rangos Multiples - HSD
A. Prueba de Rangos Multiples para Formulación
TukeyHSD(modelo2, "Formulación")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rapidez ~ Formulación + Montaje + Operador + Lote, data = datos2)
##
## $Formulación
## 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.3087315B. Prueba de Rangos Multiples para Montaje
TukeyHSD(modelo2, "Montaje")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rapidez ~ Formulación + Montaje + Operador + Lote, data = datos2)
##
## $Montaje
## diff lwr upr p adj
## β-α -3.2 -9.475865 3.075865 0.4528137
## γ-α -2.6 -8.875865 3.675865 0.6270350
## δ-α -2.8 -9.075865 3.475865 0.5669444
## ε-α 0.6 -5.675865 6.875865 0.9968743
## γ-β 0.6 -5.675865 6.875865 0.9968743
## δ-β 0.4 -5.875865 6.675865 0.9993563
## ε-β 3.8 -2.475865 10.075865 0.3087315
## δ-γ -0.2 -6.475865 6.075865 0.9999587
## ε-γ 3.2 -3.075865 9.475865 0.4528137
## ε-δ 3.4 -2.875865 9.675865 0.4006635C. Prueba de Rangos Multiples para Operador
TukeyHSD(modelo2, "Operador")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rapidez ~ Formulación + Montaje + Operador + Lote, data = datos2)
##
## $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.9906712D. Prueba de Rangos Multiples para Lote
TukeyHSD(modelo2, "Lote")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rapidez ~ Formulación + Montaje + Operador + Lote, data = datos2)
##
## $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🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
7. Verificación de supuestos del ANOVA
A. Grafica para supuestos
par( mfrow = c(2,2) )
plot(modelo, which=5)
plot(modelo, which=1)
plot(modelo, which=2)
plot(residuals(modelo2) ~ Rapidez , main="Residuals vs corridas", font.main=1, xlab="Tiempo", ylab="Residuals")
abline(h = 0, lty = 2)
B. Normalidad - Shapiro Wilks
# Prueba de normalidad de Shapiro-Wilk
shapiro.test(modelo2$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo2$residuals
## W = 0.97967, p-value = 0.8786C. Homocedasticidad - Bartlett
1. Formulación - Bartlett
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
bartlett.test(modelo2$residuals ~ Formulación, data=datos2)##
## Bartlett test of homogeneity of variances
##
## data: modelo2$residuals by Formulación
## Bartlett's K-squared = 0.31314, df = 4, p-value = 0.989Rapidez ~ Formulación + Montaje + Operador + Lote ### 2. Montaje - Bartlett
\[H_0:\sigma^2_{L_1}=\sigma^2_{L_2}=\sigma^2_{L_3}=\sigma^2_{L_3}=\sigma^2_{L_5}=\sigma^2\]
bartlett.test(modelo2$residuals ~ Montaje, data=datos2)##
## Bartlett test of homogeneity of variances
##
## data: modelo2$residuals by Montaje
## Bartlett's K-squared = 7.4545, df = 4, p-value = 0.11373. Operador - Bartlett
\[H_0:\sigma^2_{L_1}=\sigma^2_{L_2}=\sigma^2_{L_3}=\sigma^2_{L_3}=\sigma^2_{L_5}=\sigma^2\]
bartlett.test(modelo2$residuals ~ Operador, data=datos2)##
## Bartlett test of homogeneity of variances
##
## data: modelo2$residuals by Operador
## Bartlett's K-squared = 0.33032, df = 4, p-value = 0.98784. Lote - Bartlett
\[H_0:\sigma^2_{L_1}=\sigma^2_{L_2}=\sigma^2_{L_3}=\sigma^2_{L_3}=\sigma^2_{L_5}=\sigma^2\]
bartlett.test(modelo2$residuals ~ Lote, data=datos2)##
## Bartlett test of homogeneity of variances
##
## data: modelo2$residuals by Lote
## Bartlett's K-squared = 5.7344, df = 4, p-value = 0.2199D. Homocedasticidad - Levene
1. Formulación - Levene
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
#Levene
library(car)
leveneTest(modelo2$residuals ~ Formulación, data=datos2)2. Montaje - Levene
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
#Levene
library(car)
leveneTest(modelo2$residuals ~ Montaje, data=datos2)3. Operador - Levene
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
#Levene
library(car)
leveneTest(modelo2$residuals ~ Operador, data=datos2)Rapidez ~ Formulación + Montaje + Operador + Lote ### 4. Lote - Levene
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
#Levene
library(car)
leveneTest(modelo2$residuals ~ Lote, data=datos2)🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
Tematica 3: Diseños en de Dos Factore
1. Problema Ilustrativo: el efecto de cinco diferentes Formulaciones (A, B, C, D y E)
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
Problema Ilustrativo: Se tiene interés en el rendimiento de un proceso en particular. Para ello se consideran tres factores: Temperatura (A): con niveles 100, 120 y 140. Presión (B): con niveles 400 y 450. Tiempo (C): del lavado del producto en seguida del proceso de enfriamiento, con niveles 30 y 35 minutos. Se realizan tres pruebas en cada combinación de los factores. Los resultados obtenidos se presentan en la tabla, la cual servirá como base para analizar el efecto de los factores y sus interacciones sobre el rendimiento del proceso.
Factores y Niveles
| Temperatura | Presión | Tiempo | Rendimiento |
|---|---|---|---|
| 100 | 400 | 30 | 34 |
| 100 | 400 | 30 | 33 |
| 100 | 400 | 30 | 33 |
| 100 | 450 | 30 | 32 |
| 100 | 450 | 30 | 32 |
| 100 | 450 | 30 | 33 |
| 100 | 400 | 35 | 27 |
| 100 | 400 | 35 | 29 |
| 100 | 400 | 35 | 28 |
| 100 | 450 | 35 | 28 |
| 100 | 450 | 35 | 28 |
| 100 | 450 | 35 | 27 |
| 120 | 400 | 30 | 32 |
| 120 | 400 | 30 | 34 |
| 120 | 400 | 30 | 34 |
| 120 | 450 | 30 | 32 |
| 120 | 450 | 30 | 33 |
| 120 | 450 | 30 | 32 |
| 120 | 400 | 35 | 26 |
| 120 | 400 | 35 | 28 |
| 120 | 400 | 35 | 29 |
| 120 | 450 | 35 | 30 |
| 120 | 450 | 35 | 25 |
| 120 | 450 | 35 | 27 |
| 140 | 400 | 30 | 36 |
| 140 | 400 | 30 | 36 |
| 140 | 400 | 30 | 37 |
| 140 | 450 | 30 | 34 |
| 140 | 450 | 30 | 34 |
| 140 | 450 | 30 | 34 |
| 140 | 400 | 35 | 29 |
| 140 | 400 | 35 | 29 |
| 140 | 400 | 35 | 30 |
| 140 | 450 | 35 | 27 |
| 140 | 450 | 35 | 28 |
| 140 | 450 | 35 | 29 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
2. Las Hipotesis y el modelo
A.las hipotesis del problema son las siguientes -Temperatura
H0: Las Temperaturas no tienen efecto sobre el Rendimiento del proceso.
H1: Al menos una Temperaturas tiene efecto sobre el Rendimiento del proceso.
B.las hipotesis del problema son las siguientes -Preión
H0: Las Presiones no tienen efecto sobre el Rendimiento del proceso.
H1: Al menos una Presión tiene efecto sobre el Rendimiento del proceso.
C.las hipotesis del problema son las siguientes -Tiempo
H0: Las Tiempos no tienen efecto sobre el Rendimiento del proceso.
H1: Al menos una Tiempo tiene efecto sobre el Rendimiento del proceso.
D.las hipotesis del problema son las siguientes -Temperatura vs Tiempo
H0: Las combinaciones de Temperatura vs Tiempo no tienen efecto sobre el Rendimiento del proceso.
H1: Al menos una combinacion de Temperatura vs Tiempo tiene efecto sobre el Rendimiento del proceso.
E.las hipotesis del problema son las siguientes -Temperatura vs Preión
H0: Las combinaciones de Temperatura vs Presión no tienen efecto sobre el Rendimiento del proceso.
H1: Al menos una combinacion de Temperatura vs Presión tiene efecto sobre el Rendimiento del proceso.
F.las hipotesis del problema son las siguientes -Tiempo vs Presión
H0: Las combinaciones de Tiempo vs Presión no tienen efecto sobre el Rendimiento del proceso.
H1: Al menos una combinacion de Tiempo vs Presión tiene efecto sobre el Rendimiento del proceso.
G.las hipotesis del problema son las siguientes -Temperatura vs Tiempo vs Presión
H0: Las combinaciones de Temperatura vs Tiempo vs Presión no tienen efecto sobre el Rendimiento del proceso.
H1: Al menos una combinacion de Temperatura vs Tiempo vs Presión tiene efecto sobre el Rendimiento del proceso.
📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊📊
🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
3. Los datos y el DataFrame
# Datos Rendimiento
datos3 = data.frame(
Temperatura = factor(rep(c(100, 120, 140), each=12)),
Presión = factor(rep(rep(c(400, 450), each=6), times=3)),
Tiempo = factor(rep(c(30, 35), each=3, times=6)),
Rendimiento = c(34, 33, 33, 32, 32, 33, 27, 29, 28, 28, 28, 27,
32, 34, 34, 32, 33, 32, 26, 28, 29, 30, 25, 27,
36, 36, 37, 34, 34, 34, 29, 29, 30, 27, 28, 29)
)
datos3🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯🎯
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
4. Análisis de varianza
Temperatura Presión Tiempo Rendimiento
# Análisis de varianza
modelo3 = aov(Rendimiento ~ Temperatura * Presión * Tiempo, data=datos3)
summary(modelo3)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperatura 2 22.39 11.19 9.595 0.000867 ***
## Presión 1 283.36 283.36 242.881 4.71e-14 ***
## Tiempo 1 10.03 10.03 8.595 0.007294 **
## Temperatura:Presión 2 3.72 1.86 1.595 0.223632
## Temperatura:Tiempo 2 2.72 1.36 1.167 0.328447
## Presión:Tiempo 1 1.36 1.36 1.167 0.290823
## Temperatura:Presión:Tiempo 2 0.06 0.03 0.024 0.976495
## Residuals 24 28.00 1.17
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14. Interpretacion del Anovas
interpretacion del Anova
# Interpretación del ANOVA
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionanova_results <- summary(modelo3)[[1]] %>%
as.data.frame() %>%
mutate(Significant = ifelse(`Pr(>F)` < 0.05, "Yes", "No"))
anova_results🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
5. Boxplot para los Factores Prinicpales
A. Boxplot para Presión
# Boxplot para Catalizador con colores
boxplot(Rendimiento ~ Presión, data=datos3, main="Boxplot de Rendimiento por Presión", xlab="Presión", ylab="Rapidez", col=c("red2", "blue2", "yellow2", "orange2","green2"))
B. Boxplot para Temperatura
# Boxplot para Catalizador con colores
boxplot(Rendimiento ~ Temperatura, data=datos3, main="Boxplot de Rendimiento por Temperatura", xlab="Temperatura", ylab="Rapidez", col=c("red2", "blue2", "yellow2", "orange2","green2"))
C. Boxplot para Tiempo
# Boxplot para Catalizador con colores
boxplot(Rendimiento ~ Tiempo, data=datos3, main="Boxplot de Rendimiento por Tiempo", xlab="Tiempo", ylab="Rapidez", col=c("red2", "blue2", "yellow2", "orange2","green2"))
# 6. Boxplot para los Factores
Interacciones Dobles
A. Boxplot para Temperatura * Presión
# Boxplot bivariante para Temperatura y Presión
boxplot(Rendimiento ~ Temperatura + Presión, data=datos3, main="Boxplot de Rendimiento por Temperatura y Presión", xlab="Temperatura y Presión", ylab="Rendimiento", col=c("red2", "blue2", "yellow2", "orange2","green2"))
## B. Boxplot para Temperatura *
Tiempo
# Boxplot bivariante para Temperatura y Tiempo
boxplot(Rendimiento ~ Temperatura + Tiempo, data=datos3, main="Boxplot de Rendimiento por Temperatura y Tiempo", xlab="Temperatura y Tiempo", ylab="Rendimiento", col=c("red2", "blue2", "yellow2", "orange2","green2"))
## C. Boxplot para Presión *
Tiempo
# Boxplot bivariante para Presión y Tiempo
boxplot(Rendimiento ~ Presión + Tiempo, data=datos3, main="Boxplot de Rendimiento por Presión y Tiempo", xlab="Presión y Tiempo", ylab="Rendimiento", col=c("red2", "blue2", "yellow2", "orange2","green2"))
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
6. Prueba de Rangos Multiples - HSD
A. Prueba de Rangos Multiples para Formulación
TukeyHSD(modelo3, "Temperatura")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rendimiento ~ Temperatura * Presión * Tiempo, data = datos3)
##
## $Temperatura
## diff lwr upr p adj
## 120-100 -0.1666667 -1.2678668 0.9345334 0.9245219
## 140-100 1.5833333 0.4821332 2.6845334 0.0040499
## 140-120 1.7500000 0.6487999 2.8512001 0.0015949B. Prueba de Rangos Multiples para Presión
TukeyHSD(modelo3, "Presión")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rendimiento ~ Temperatura * Presión * Tiempo, data = datos3)
##
## $Presión
## diff lwr upr p adj
## 450-400 -5.611111 -6.354199 -4.868023 0C. Prueba de Rangos Multiples para Operador
TukeyHSD(modelo3, "Tiempo")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rendimiento ~ Temperatura * Presión * Tiempo, data = datos3)
##
## $Tiempo
## diff lwr upr p adj
## 35-30 -1.055556 -1.798644 -0.3124672 0.0072938🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
7. Graficas de Interacción
A. Prueba de Rangos Multiples para Temperatura*Presion
# grafico de interacción Temperatura * Presión
interaction.plot(datos3$Temperatura, datos3$Presión, datos3$Rendimiento, col=c("red2", "blue2", "yellow2", "orange2","green2"), main="Interacción entre Temperatura y Presión", xlab="Temperatura", ylab="Rendimiento")
B. Prueba de Rangos Multiples para Temperatura*Tiempo
# grafico de interacción Temperatura * Tiempo
interaction.plot(datos3$Temperatura, datos3$Tiempo, datos3$Rendimiento, col=c("red2", "blue2", "yellow2", "orange2","green2"), main="Interacción entre Temperatura y Tiempo", xlab="Temperatura", ylab="Rendimiento")
C. Prueba de Rangos Multiples para Presión*Tiempo
# grafico de interacción Presión * Tiempo
interaction.plot(datos3$Presión, datos3$Tiempo, datos3$Rendimiento, col=c("red2", "blue2", "yellow2", "orange2","green2"), main="Interacción entre Presión y Tiempo", xlab="Presión", ylab="Rendimiento")
8. Graficas de Efectos Principales - HSD
A. Gráficos de Medias para Temperatura
# grafico de medias HSD Temperatura
plot(TukeyHSD(modelo3, "Temperatura"), col="blue")
B. Gráficos de Medias para Presión
# grafico de medias HSD Presión
plot(TukeyHSD(modelo3, "Presión"), col="blue")
C. Gráficos de Medias para Tiempo
# grafico de medias HSD Tiempo
plot(TukeyHSD(modelo3, "Tiempo"), col="blue")
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
9. Verificación de supuestos del ANOVA
A. Grafica para supuestos
# supuestos del modelo
par(mfrow=c(2,2))
plot(modelo3)
B. Normalidad - Shapiro Wilks
# Prueba de normalidad de Shapiro-Wilk
shapiro.test(modelo3$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo3$residuals
## W = 0.94271, p-value = 0.0617C. Homocedasticidad - Bartlett
1. Temperatura - Bartlett
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
bartlett.test(modelo3$residuals ~ Temperatura, data=datos3)##
## Bartlett test of homogeneity of variances
##
## data: modelo3$residuals by Temperatura
## Bartlett's K-squared = 11.44, df = 2, p-value = 0.003282. Presión - Bartlett
\[H_0:\sigma^2_{L_1}=\sigma^2_{L_2}=\sigma^2_{L_3}=\sigma^2_{L_3}=\sigma^2_{L_5}=\sigma^2\]
bartlett.test(modelo3$residuals ~ Presión, data=datos3)##
## Bartlett test of homogeneity of variances
##
## data: modelo3$residuals by Presión
## Bartlett's K-squared = 7.9804, df = 1, p-value = 0.0047293. Operador - Bartlett
\[H_0:\sigma^2_{L_1}=\sigma^2_{L_2}=\sigma^2_{L_3}=\sigma^2_{L_3}=\sigma^2_{L_5}=\sigma^2\]
bartlett.test(modelo3$residuals ~ Tiempo, data=datos3)##
## Bartlett test of homogeneity of variances
##
## data: modelo3$residuals by Tiempo
## Bartlett's K-squared = 0.6103, df = 1, p-value = 0.4347D. Homocedasticidad - Levene
1. Temperatura - Levene
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
#Levene
library(car)
leveneTest(modelo3$residuals ~ Temperatura, data=datos3)2. Presión - Levene
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
#Levene
library(car)
leveneTest(modelo3$residuals ~ Presión, data=datos3)3. Tiempo - Levene
\[H_0:\sigma^2_A=\sigma^2_B=\sigma^2_C=\sigma^2_D=\sigma^2_E=\sigma^2\]
#Levene
library(car)
leveneTest(modelo3$residuals ~ Tiempo, data=datos3)🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
Tematica 4. Experimento \(2^5\) no replicado - Ejemplo, pág. 199
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
Ejemplo integrador
En una planta donde se fabrican semiconductores se quiere mejorar el rendimiento del proceso vía diseño de experimentos. De acuerdo con la experiencia del grupo de mejora, los factores que podrían tener mayor influencia sobre la variable de respuesta (rendimiento), así como los niveles de prueba utilizados son los siguientes:
A = Nivel de la abertura (pequeña, grande) = (-1,1).
B = Tiempo de exposición (20% abajo, 20% arriba)= (-1,1)
C = Tiempo de revelado (30 seg, 45 seg)= (-1,1)
D = Dimensión de la máscara (pequeña, grande)= (-1,1)
E = Tiempo de grabado (14.5 min, 15.5 min)= (-1,1)
Se decide correr un experimento 25 con una sola corrida o réplica para estudiar estos cinco factores. Se hacen las 32 corridas a nivel proceso.
Paso 1. Datos
# Supongamos que tenemos un vector “Y” con los rendimientos de las 32 corridas,
# en el orden de la tabla de Yates (orden lexicográfico de factores codificados -1/1).
Y <- c(7, 9, 34, 55, 16, 20, 41, 60,
8, 10, 32, 50, 18, 21, 44, 61,
8, 12, 35, 52, 15, 22, 45, 65,
6, 10, 30, 53, 15, 20, 41, 63)
# Creamos los factores codificados
# expand.grid genera todas las combinaciones
df <- expand.grid(
A = c(-1, 1),
B = c(-1, 1),
C = c(-1, 1),
D = c(-1, 1),
E = c(-1, 1)
)
# Ordenar correctamente al estilo de Yates si es necesario
# (dependiendo del orden en que definiste Y)
# Luego añadir la variable de respuesta
df$rendimiento <- Y
# Revisar
str(df)## 'data.frame': 32 obs. of 6 variables:
## $ A : num -1 1 -1 1 -1 1 -1 1 -1 1 ...
## $ B : num -1 -1 1 1 -1 -1 1 1 -1 -1 ...
## $ C : num -1 -1 -1 -1 1 1 1 1 -1 -1 ...
## $ D : num -1 -1 -1 -1 -1 -1 -1 -1 1 1 ...
## $ E : num -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ rendimiento: num 7 9 34 55 16 20 41 60 8 10 ...
## - attr(*, "out.attrs")=List of 2
## ..$ dim : Named int [1:5] 2 2 2 2 2
## .. ..- attr(*, "names")= chr [1:5] "A" "B" "C" "D" ...
## ..$ dimnames:List of 5
## .. ..$ A: chr [1:2] "A=-1" "A= 1"
## .. ..$ B: chr [1:2] "B=-1" "B= 1"
## .. ..$ C: chr [1:2] "C=-1" "C= 1"
## .. ..$ D: chr [1:2] "D=-1" "D= 1"
## .. ..$ E: chr [1:2] "E=-1" "E= 1"head(df)# Ver el data.frame completo
print(df)## A B C D E rendimiento
## 1 -1 -1 -1 -1 -1 7
## 2 1 -1 -1 -1 -1 9
## 3 -1 1 -1 -1 -1 34
## 4 1 1 -1 -1 -1 55
## 5 -1 -1 1 -1 -1 16
## 6 1 -1 1 -1 -1 20
## 7 -1 1 1 -1 -1 41
## 8 1 1 1 -1 -1 60
## 9 -1 -1 -1 1 -1 8
## 10 1 -1 -1 1 -1 10
## 11 -1 1 -1 1 -1 32
## 12 1 1 -1 1 -1 50
## 13 -1 -1 1 1 -1 18
## 14 1 -1 1 1 -1 21
## 15 -1 1 1 1 -1 44
## 16 1 1 1 1 -1 61
## 17 -1 -1 -1 -1 1 8
## 18 1 -1 -1 -1 1 12
## 19 -1 1 -1 -1 1 35
## 20 1 1 -1 -1 1 52
## 21 -1 -1 1 -1 1 15
## 22 1 -1 1 -1 1 22
## 23 -1 1 1 -1 1 45
## 24 1 1 1 -1 1 65
## 25 -1 -1 -1 1 1 6
## 26 1 -1 -1 1 1 10
## 27 -1 1 -1 1 1 30
## 28 1 1 -1 1 1 53
## 29 -1 -1 1 1 1 15
## 30 1 -1 1 1 1 20
## 31 -1 1 1 1 1 41
## 32 1 1 1 1 1 63🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
Paso 2. Boxplot para efectos principales
# Boxplot para cada factor
par(mfrow=c(2,3)) # Para ver varios gráficos juntos
boxplot(rendimiento ~ A, data=df, main="Boxplot de A", col=c("lightblue", "lightgreen"))
boxplot(rendimiento ~ B, data=df, main="Boxplot de B", col=c("lightblue", "lightgreen"))
boxplot(rendimiento ~ C, data=df, main="Boxplot de C", col=c("lightblue", "lightgreen"))
boxplot(rendimiento ~ D, data=df, main="Boxplot de D", col=c("lightblue", "lightgreen"))
boxplot(rendimiento ~ E, data=df, main="Boxplot de E", col=c("lightblue", "lightgreen"))
par(mfrow=c(1,1)) # Volver a un solo gráfico
Paso 3. Boxplot para efectos de interacción
# Boxplot para interacciones dobles
par(mfrow=c(2,3)) # Para ver varios gráficos juntos
boxplot(rendimiento ~ A + B, data=df, main="Boxplot de A*B", col=c("lightblue", "lightgreen"))
boxplot(rendimiento ~ A + C, data=df, main="Boxplot de A*C", col=c("lightblue", "lightgreen"))
boxplot(rendimiento ~ A + D, data=df, main="Boxplot de A*D", col=c("lightblue", "lightgreen"))
boxplot(rendimiento ~ A + E, data=df, main="Boxplot de A*E", col=c("lightblue", "lightgreen"))
boxplot(rendimiento ~ B + C, data=df, main="Boxplot de B*C", col=c("lightblue", "lightgreen"))
par(mfrow=c(1,1)) # Volver a un solo gráfico
# Efectos principales
library(dplyr)
efectos <- df %>%
summarise(
efecto_A = mean(rendimiento[A == 1]) - mean(rendimiento[A == -1]),
efecto_B = mean(rendimiento[B == 1]) - mean(rendimiento[B == -1]),
efecto_C = mean(rendimiento[C == 1]) - mean(rendimiento[C == -1]),
efecto_D = mean(rendimiento[D == 1]) - mean(rendimiento[D == -1]),
efecto_E = mean(rendimiento[E == 1]) - mean(rendimiento[E == -1])
)
print(efectos)## efecto_A efecto_B efecto_C efecto_D efecto_E
## 1 11.75 34 9.75 -0.875 0.375Efectos de interacción dobles
# Efectos de interacción dobles
interacciones <- df %>%
summarise(
efecto_AB = mean(rendimiento[A == 1 & B == 1]) - mean(rendimiento[A == 1 & B == -1]) -
(mean(rendimiento[A == -1 & B == 1]) - mean(rendimiento[A == -1 & B == -1])),
efecto_AC = mean(rendimiento[A == 1 & C == 1]) - mean(rendimiento[A == 1 & C == -1]) -
(mean(rendimiento[A == -1 & C == 1]) - mean(rendimiento[A == -1 & C == -1])),
efecto_AD = mean(rendimiento[A == 1 & D == 1]) - mean(rendimiento[A == 1 & D == -1]) -
(mean(rendimiento[A == -1 & D == 1]) - mean(rendimiento[A == -1 & D == -1])),
efecto_AE = mean(rendimiento[A == 1 & E == 1]) - mean(rendimiento[A == 1 & E == -1]) -
(mean(rendimiento[A == -1 & E == 1]) - mean(rendimiento[A == -1 & E == -1])),
efecto_BC = mean(rendimiento[B == 1 & C == 1]) - mean(rendimiento[B == 1 & C == -1]) -
(mean(rendimiento[B == -1 & C == 1]) - mean(rendimiento[B == -1 & C == -1]))
)
print(interacciones)## efecto_AB efecto_AC efecto_AD efecto_AE efecto_BC
## 1 15.75 0.75 0 2 0.25Tabla Anova
# Análisis de varianza
modelo_full <- aov(rendimiento ~ (A + B + C + D + E)^2, data=df)
summary(modelo_full)## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 1104 1104 482.512 2.24e-13 ***
## B 1 9248 9248 4040.082 < 2e-16 ***
## C 1 761 761 332.232 3.98e-12 ***
## D 1 6 6 2.676 0.1214
## E 1 1 1 0.491 0.4933
## A:B 1 496 496 216.737 1.01e-10 ***
## A:C 1 1 1 0.491 0.4933
## A:D 1 0 0 0.000 1.0000
## A:E 1 8 8 3.495 0.0800 .
## B:C 1 0 0 0.055 0.8182
## B:D 1 5 5 1.966 0.1800
## B:E 1 2 2 0.874 0.3638
## C:D 1 4 4 1.966 0.1800
## C:E 1 0 0 0.218 0.6465
## D:E 1 10 10 4.423 0.0516 .
## Residuals 16 37 2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# grafico de medias HSD
# Grafico de medias HSD
install.packages("agricolae")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.5'
## (as 'lib' is unspecified)library(agricolae)
hsd <- HSD.test(modelo_full, "A", group = TRUE)
plot(hsd, main = "HSD para A")
hsd <- HSD.test(modelo_full, "B", group = TRUE)
plot(hsd, main = "HSD para B")
hsd <- HSD.test(modelo_full, "C", group = TRUE)
plot(hsd, main = "HSD para C")
hsd <- HSD.test(modelo_full, "D", group = TRUE)
plot(hsd, main = "HSD para D")
hsd <- HSD.test(modelo_full, "E", group = TRUE)
plot(hsd, main = "HSD para E")
elegir el mejor tratamiento
# Identificar el mejor tratamiento
best_treatment <- df[which.max(df$rendimiento), ]
print(best_treatment)## A B C D E rendimiento
## 24 1 1 1 -1 1 65Prueba de normalidad de Shapiro-Wilk
# Prueba de normalidad de Shapiro-Wilk
shapiro.test(modelo_full$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo_full$residuals
## W = 0.96465, p-value = 0.366# Vector de rendimientos (orden Yates)
Y <- c(7, 9, 34, 55, 16, 20, 41, 60,
8, 10, 32, 50, 18, 21, 44, 61,
8, 12, 35, 52, 15, 22, 45, 65,
6, 10, 30, 53, 15, 20, 41, 63)
# Crear los factores codificados
df <- expand.grid(
A = c(-1, 1),
B = c(-1, 1),
C = c(-1, 1),
D = c(-1, 1),
E = c(-1, 1)
)
# Añadir la respuesta
df$rendimiento <- Y
# Verificación
head(df)# -----------------------------------------
# Gráfico de efectos principales (versión base R)
# -----------------------------------------
# Calcular medias por nivel de cada factor
mediaA <- tapply(df$rendimiento, df$A, mean)
mediaB <- tapply(df$rendimiento, df$B, mean)
mediaC <- tapply(df$rendimiento, df$C, mean)
mediaD <- tapply(df$rendimiento, df$D, mean)
mediaE <- tapply(df$rendimiento, df$E, mean)
# Crear el gráfico
par(mfrow = c(1, 5), bg = "gray90", mar = c(4, 4, 3, 1))
plot(c(-1, 1), mediaA, type = "b", col = "blue", lwd = 2,
xaxt = "n", xlab = "A", ylab = "Rendimiento",
main = "A")
axis(1, at = c(-1, 1), labels = c("Bajo", "Alto"))
plot(c(-1, 1), mediaB, type = "b", col = "blue", lwd = 2,
xaxt = "n", xlab = "B", ylab = "", main = "B")
axis(1, at = c(-1, 1), labels = c("Bajo", "Alto"))
plot(c(-1, 1), mediaC, type = "b", col = "blue", lwd = 2,
xaxt = "n", xlab = "C", ylab = "", main = "C")
axis(1, at = c(-1, 1), labels = c("Bajo", "Alto"))
plot(c(-1, 1), mediaD, type = "b", col = "blue", lwd = 2,
xaxt = "n", xlab = "D", ylab = "", main = "D")
axis(1, at = c(-1, 1), labels = c("Bajo", "Alto"))
plot(c(-1, 1), mediaE, type = "b", col = "blue", lwd = 2,
xaxt = "n", xlab = "E", ylab = "", main = "E")
axis(1, at = c(-1, 1), labels = c("Bajo", "Alto"))
title("Gráfica de Efectos Principales para Rendimiento", outer = TRUE)
par(mfrow = c(1, 1)) # Resetear la disposición de gráficos# en un mismo plano las figuras de efevtos principales
library(ggplot2)
install.packages("reshape2")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.5'
## (as 'lib' is unspecified)library(reshape2)
# Transformar los datos para ggplot
df_melt <- melt(df, id.vars = "rendimiento")
# Calcular medias
medias <- aggregate(rendimiento ~ variable + value, data = df_melt, mean)
# Graficar
ggplot(medias, aes(x = value, y = rendimiento, group = variable, color = variable)) +
geom_line() +
geom_point() +
facet_wrap(~ variable, scales = "free_x") +
labs(title = "Gráfica de Efectos Principales para Rendimiento",
x = "Nivel del Factor", y = "Rendimiento") +
theme_minimal()
# corrige el codigo anterior
# Gráfica de efectos principales usando ggplot2
library(ggplot2)
library(reshape2)
# Transformar los datos para ggplot
df_melt <- melt(df, id.vars = "rendimiento")
# Calcular medias
medias <- aggregate(rendimiento ~ variable + value, data = df_melt, mean)
# Graficar
ggplot(medias, aes(x = value, y = rendimiento, group = variable, color = variable)) +
geom_line() +
geom_point() +
facet_wrap(~ variable, scales = "free_x") +
labs(title = "Gráfica de Efectos Principales para Rendimiento",
x = "Nivel del Factor", y = "Rendimiento") +
theme_minimal()
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
# Instalar y cargar el paquete FrF2 (solo si no está instalado)
if (!require(FrF2)) {
install.packages("FrF2")
library(FrF2)
}## Loading required package: FrF2## Loading required package: 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# Cargar otros paquetes útiles
if (!require(ggplot2)) {
install.packages("ggplot2")
library(ggplot2)
}
# Datos de rendimiento
Y <- c(7, 9, 34, 55, 16, 20, 41, 60,
8, 10, 32, 50, 18, 21, 44, 61,
8, 12, 35, 52, 15, 22, 45, 65,
6, 10, 30, 53, 15, 20, 41, 63)
# Crear el diseño factorial usando FrF2
design <- FrF2(
nruns = 32,
nfactors = 5,
factor.names = list(
A = c("-1", "1"),
B = c("-1", "1"),
C = c("-1", "1"),
D = c("-1", "1"),
E = c("-1", "1")
),
randomize = FALSE
)## creating full factorial with 32 runs ...# VERIFICAR que el diseño coincide con tu orden de datos
print("Diseño factorial creado:")## [1] "Diseño factorial creado:"print(design)## A B C D E
## 1 -1 -1 -1 -1 -1
## 2 1 -1 -1 -1 -1
## 3 -1 1 -1 -1 -1
## 4 1 1 -1 -1 -1
## 5 -1 -1 1 -1 -1
## 6 1 -1 1 -1 -1
## 7 -1 1 1 -1 -1
## 8 1 1 1 -1 -1
## 9 -1 -1 -1 1 -1
## 10 1 -1 -1 1 -1
## 11 -1 1 -1 1 -1
## 12 1 1 -1 1 -1
## 13 -1 -1 1 1 -1
## 14 1 -1 1 1 -1
## 15 -1 1 1 1 -1
## 16 1 1 1 1 -1
## 17 -1 -1 -1 -1 1
## 18 1 -1 -1 -1 1
## 19 -1 1 -1 -1 1
## 20 1 1 -1 -1 1
## 21 -1 -1 1 -1 1
## 22 1 -1 1 -1 1
## 23 -1 1 1 -1 1
## 24 1 1 1 -1 1
## 25 -1 -1 -1 1 1
## 26 1 -1 -1 1 1
## 27 -1 1 -1 1 1
## 28 1 1 -1 1 1
## 29 -1 -1 1 1 1
## 30 1 -1 1 1 1
## 31 -1 1 1 1 1
## 32 1 1 1 1 1
## class=design, type= full factorial# AGREGAR la variable de respuesta al diseño
design$Y <- Y
# Verificar que se agregó correctamente
print("Diseño con variable de respuesta:")## [1] "Diseño con variable de respuesta:"print(response.names(design))## NULLprint(head(design))## A B C D E Y
## 1 -1 -1 -1 -1 -1 7
## 2 1 -1 -1 -1 -1 9
## 3 -1 1 -1 -1 -1 34
## 4 1 1 -1 -1 -1 55
## 5 -1 -1 1 -1 -1 16
## 6 1 -1 1 -1 -1 20# ANÁLISIS ESTADÍSTICO COMPLETO
cat("\n=== ANÁLISIS ESTADÍSTICO ===\n")##
## === ANÁLISIS ESTADÍSTICO ===# 1. Modelo lineal con interacciones
modelo <- lm(Y ~ (A + B + C + D + E)^2, data = design)
print("Resumen del modelo lineal:")## [1] "Resumen del modelo lineal:"print(summary(modelo))##
## Call:
## lm.default(formula = Y ~ (A + B + C + D + E)^2, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8125 -0.6250 0.1875 0.5000 2.9375
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.056e+01 2.675e-01 114.271 < 2e-16 ***
## A1 5.875e+00 2.675e-01 21.966 2.24e-13 ***
## B1 1.700e+01 2.675e-01 63.562 < 2e-16 ***
## C1 4.875e+00 2.675e-01 18.227 3.98e-12 ***
## D1 -4.375e-01 2.675e-01 -1.636 0.1214
## E1 1.875e-01 2.675e-01 0.701 0.4933
## A1:B1 3.937e+00 2.675e-01 14.722 1.01e-10 ***
## A1:C1 1.875e-01 2.675e-01 0.701 0.4933
## A1:D1 -7.016e-16 2.675e-01 0.000 1.0000
## A1:E1 5.000e-01 2.675e-01 1.869 0.0800 .
## B1:C1 6.250e-02 2.675e-01 0.234 0.8182
## B1:D1 -3.750e-01 2.675e-01 -1.402 0.1800
## B1:E1 2.500e-01 2.675e-01 0.935 0.3638
## C1:D1 3.750e-01 2.675e-01 1.402 0.1800
## C1:E1 1.250e-01 2.675e-01 0.467 0.6465
## D1:E1 -5.625e-01 2.675e-01 -2.103 0.0516 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.513 on 16 degrees of freedom
## Multiple R-squared: 0.9969, Adjusted R-squared: 0.9939
## F-statistic: 339.2 on 15 and 16 DF, p-value: < 2.2e-16# 2. ANOVA
cat("\n--- ANÁLISIS DE VARIANZA (ANOVA) ---\n")##
## --- ANÁLISIS DE VARIANZA (ANOVA) ---anova_result <- anova(modelo)
print(anova_result)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 1104.5 1104.5 482.5119 2.244e-13 ***
## B 1 9248.0 9248.0 4040.0819 < 2.2e-16 ***
## C 1 760.5 760.5 332.2321 3.983e-12 ***
## D 1 6.1 6.1 2.6758 0.12140
## E 1 1.1 1.1 0.4915 0.49334
## A:B 1 496.1 496.1 216.7372 1.011e-10 ***
## A:C 1 1.1 1.1 0.4915 0.49334
## A:D 1 0.0 0.0 0.0000 1.00000
## A:E 1 8.0 8.0 3.4949 0.07997 .
## B:C 1 0.1 0.1 0.0546 0.81820
## B:D 1 4.5 4.5 1.9659 0.17999
## B:E 1 2.0 2.0 0.8737 0.36382
## C:D 1 4.5 4.5 1.9659 0.17999
## C:E 1 0.5 0.5 0.2184 0.64654
## D:E 1 10.1 10.1 4.4232 0.05163 .
## Residuals 16 36.6 2.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Identificar efectos significativos
significativos <- which(anova_result$`Pr(>F)` < 0.05)
if (length(significativos) > 0) {
cat("\nEFECTOS SIGNIFICATIVOS (p < 0.05):\n")
print(rownames(anova_result)[significativos])
} else {
cat("\nNo se encontraron efectos significativos al nivel alpha = 0.05\n")
}##
## EFECTOS SIGNIFICATIVOS (p < 0.05):
## [1] "A" "B" "C" "A:B"# GRÁFICOS PRINCIPALES
cat("\n=== CREANDO GRÁFICOS ===\n")##
## === CREANDO GRÁFICOS ===# Configurar área de gráficos
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
# 1. GRÁFICA DE INTERACCIONES (IAPlot) - ¡ESTA DEBERÍA FUNCIONAR AHORA!
tryCatch({
IAPlot(design,
response = "Y",
main = "Gráfica de Interacciones para Rendimiento",
cex.main = 0.9)
cat("✓ Gráfica de interacciones creada exitosamente\n")
}, error = function(e) {
cat("✗ Error en IAPlot:", e$message, "\n")
})## ✗ Error en IAPlot: The design obj must have at least one response.# 2. GRÁFICA DE EFECTOS PRINCIPALES
tryCatch({
MEPlot(design,
response = "Y",
main = "Efectos Principales del Rendimiento")
cat("✓ Gráfica de efectos principales creada exitosamente\n")
}, error = function(e) {
cat("✗ Error en MEPlot:", e$message, "\n")
})## ✗ Error en MEPlot: The design obj must have at least one response.# 3. GRÁFICA DE PROBABILIDAD NORMAL (Daniel Plot)
tryCatch({
DanielPlot(design,
response = "Y",
main = "Gráfica de Probabilidad Normal")
cat("✓ Gráfica de probabilidad normal creada exitosamente\n")
}, error = function(e) {
cat("✗ Error en DanielPlot:", e$message, "\n")
})## ✗ Error en DanielPlot: The design fit must have at least one response.# 4. GRÁFICO DE INTERACCIÓN ESPECÍFICO (A vs B)
tryCatch({
interaction.plot(design$A, design$B, design$Y,
main = "Interacción A vs B",
xlab = "Factor A",
ylab = "Rendimiento Y",
col = c("blue", "red"),
lwd = 2,
trace.label = "Factor B")
cat("✓ Gráfico de interacción A vs B creado exitosamente\n")
}, error = function(e) {
cat("✗ Error en interaction.plot:", e$message, "\n")
})## ✓ Gráfico de interacción A vs B creado exitosamente# Restaurar configuración gráfica
par(mfrow = c(1, 1))
# GRÁFICOS ALTERNATIVOS CON GGPLOT2
cat("\n=== GRÁFICOS ALTERNATIVOS CON GGPLOT2 ===\n")##
## === GRÁFICOS ALTERNATIVOS CON GGPLOT2 ===# Convertir a data frame para ggplot
design_df <- as.data.frame(design)
# 1. Boxplot de efectos principales
p1 <- ggplot(design_df, aes(x = A, y = Y)) +
geom_boxplot(fill = "lightblue", alpha = 0.7) +
stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "red") +
labs(title = "Efecto del Factor A sobre el Rendimiento",
x = "Nivel del Factor A",
y = "Rendimiento") +
theme_minimal()
print(p1)
# 2. Gráfico de interacción A vs B
p2 <- ggplot(design_df, aes(x = A, y = Y, color = B, group = B)) +
stat_summary(fun = mean, geom = "line", size = 1) +
stat_summary(fun = mean, geom = "point", size = 3) +
labs(title = "Interacción entre Factor A y B",
x = "Factor A",
y = "Rendimiento Promedio",
color = "Factor B") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.print(p2)
# ANÁLISIS ADICIONAL: EFECTOS ESTIMADOS
cat("\n=== EFECTOS ESTIMADOS ===\n")##
## === EFECTOS ESTIMADOS ===# Calcular efectos principales
efectos <- model.matrix(~ (A + B + C + D + E)^2, data = design)
coeficientes <- coef(modelo)
# Mostrar efectos principales ordenados
efectos_principales <- coeficientes[2:6]
names(efectos_principales) <- c("A", "B", "C", "D", "E")
cat("Efectos principales ordenados por magnitud:\n")## Efectos principales ordenados por magnitud:print(sort(abs(efectos_principales), decreasing = TRUE))## B A C D E
## 17.0000 5.8750 4.8750 0.4375 0.1875# RESUMEN FINAL
cat("\n=== RESUMEN DEL ANÁLISIS ===\n")##
## === RESUMEN DEL ANÁLISIS ===cat("Diseño: Factorial 2^5 con 32 corridas\n")## Diseño: Factorial 2^5 con 32 corridascat("Variable respuesta: Rendimiento (Y)\n")## Variable respuesta: Rendimiento (Y)cat("Rango de rendimiento:", range(Y), "\n")## Rango de rendimiento: 6 65cat("Media del rendimiento:", round(mean(Y), 2), "\n")## Media del rendimiento: 30.56cat("Desviación estándar:", round(sd(Y), 2), "\n")## Desviación estándar: 19.41cat("Número de efectos significativos:", length(significativos), "\n")## Número de efectos significativos: 4# Guardar el diseño con los resultados si es necesario
# write.csv(design_df, "diseno_factorial_con_resultados.csv")# Instalar y cargar paquetes necesarios
if (!require(FrF2)) {
install.packages("FrF2")
library(FrF2)
}
if (!require(ggplot2)) {
install.packages("ggplot2")
library(ggplot2)
}
if (!require(gridExtra)) {
install.packages("gridExtra")
library(gridExtra)
}## Loading required package: gridExtra##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine# Datos de rendimiento
Y <- c(7, 9, 34, 55, 16, 20, 41, 60,
8, 10, 32, 50, 18, 21, 44, 61,
8, 12, 35, 52, 15, 22, 45, 65,
6, 10, 30, 53, 15, 20, 41, 63)
# Crear el diseño factorial
design <- FrF2(
nruns = 32,
nfactors = 5,
factor.names = list(
A = c("-1", "1"),
B = c("-1", "1"),
C = c("-1", "1"),
D = c("-1", "1"),
E = c("-1", "1")
),
randomize = FALSE
)## creating full factorial with 32 runs ...# Agregar la variable de respuesta
design$Y <- Y
# Convertir a data frame para ggplot
design_df <- as.data.frame(design)
# MÉTODO 1: GRÁFICAS DE INTERACCIÓN CON FrF2 (Todas las interacciones)
cat("=== GRÁFICAS DE INTERACCIÓN CON FrF2 ===\n")## === GRÁFICAS DE INTERACCIÓN CON FrF2 ===# Configurar área de gráficos para múltiples gráficos
par(mfrow = c(1, 1))
# IAPlot muestra todas las interacciones de 2 factores
tryCatch({
IAPlot(design,
response = "Y",
main = "TODAS las Interacciones de 2 Factores - Rendimiento",
cex.main = 1)
cat("✓ Gráfica completa de interacciones creada\n")
}, error = function(e) {
cat("Error en IAPlot:", e$message, "\n")
})## Error en IAPlot: The design obj must have at least one response.# MÉTODO 2: GRÁFICAS INDIVIDUALES DE INTERACCIÓN CON BASE R
cat("\n=== GRÁFICAS INDIVIDUALES DE INTERACCIÓN (BASE R) ===\n")##
## === GRÁFICAS INDIVIDUALES DE INTERACCIÓN (BASE R) ===# Configurar para mostrar múltiples gráficos
par(mfrow = c(3, 4), mar = c(4, 4, 2, 1), oma = c(0, 0, 2, 0))
# Lista de todas las combinaciones de 2 factores
factores <- c("A", "B", "C", "D", "E")
combinaciones <- combn(factores, 2, simplify = FALSE)
# Crear gráficos de interacción para cada par de factores
for (i in seq_along(combinaciones)) {
factor1 <- combinaciones[[i]][1]
factor2 <- combinaciones[[i]][2]
interaction.plot(design_df[[factor1]], design_df[[factor2]], design_df$Y,
main = paste("Interacción", factor1, "vs", factor2),
xlab = factor1,
ylab = "Rendimiento",
col = 1:2,
lwd = 2,
trace.label = factor2)
}
title("TODAS LAS INTERACCIONES DE 2 FACTORES - Diseño 2^5", outer = TRUE)
# Restaurar configuración
par(mfrow = c(1, 1))
# MÉTODO 3: GRÁFICAS CON GGPLOT2 (MÁS ELEGANTES)
cat("\n=== GRÁFICAS DE INTERACCIÓN CON GGPLOT2 ===\n")##
## === GRÁFICAS DE INTERACCIÓN CON GGPLOT2 ===# Función para crear gráficos de interacción individuales
crear_grafico_interaccion <- function(fac1, fac2) {
ggplot(design_df, aes(x = .data[[fac1]], y = Y, color = .data[[fac2]], group = .data[[fac2]])) +
stat_summary(fun = mean, geom = "line", size = 1) +
stat_summary(fun = mean, geom = "point", size = 2.5) +
stat_summary(fun = mean, geom = "point", size = 2, color = "white") +
labs(title = paste("Interacción", fac1, "vs", fac2),
x = paste("Factor", fac1),
y = "Rendimiento Promedio",
color = paste("Factor", fac2)) +
theme_minimal() +
theme(plot.title = element_text(size = 10, hjust = 0.5))
}
# Crear todas las combinaciones de gráficos
graficos <- list()
for (i in seq_along(combinaciones)) {
fac1 <- combinaciones[[i]][1]
fac2 <- combinaciones[[i]][2]
graficos[[i]] <- crear_grafico_interaccion(fac1, fac2)
}
# Mostrar todos los gráficos en una cuadrícula
cat("Mostrando todas las interacciones en cuadrícula...\n")## Mostrando todas las interacciones en cuadrícula...grid.arrange(grobs = graficos, ncol = 3,
top = "TODAS LAS INTERACCIONES DE 2 FACTORES - Diseño Factorial 2^5")
# MÉTODO 4: MATRIZ DE INTERACCIONES (HEATMAP)
cat("\n=== MATRIZ DE INTERACCIONES (HEATMAP) ===\n")##
## === MATRIZ DE INTERACCIONES (HEATMAP) ===# Calcular la fuerza de cada interacción
calcular_fuerza_interaccion <- function(fac1, fac2) {
formula <- as.formula(paste("Y ~", fac1, "*", fac2))
modelo_interaccion <- lm(formula, data = design_df)
coeficiente_interaccion <- coef(modelo_interaccion)[paste0(fac1, ":", fac2)]
return(abs(coeficiente_interaccion))
}
# Crear matriz de fuerzas de interacción
matriz_interacciones <- matrix(0, nrow = 5, ncol = 5)
rownames(matriz_interacciones) <- factores
colnames(matriz_interacciones) <- factores
for (i in 1:4) {
for (j in (i+1):5) {
fuerza <- calcular_fuerza_interaccion(factores[i], factores[j])
matriz_interacciones[i, j] <- fuerza
matriz_interacciones[j, i] <- fuerza
}
}
# Graficar heatmap de interacciones
if (!require(reshape2)) {
install.packages("reshape2")
library(reshape2)
}
matriz_melt <- melt(matriz_interacciones)
matriz_melt <- matriz_melt[matriz_melt$Var1 != matriz_melt$Var2, ] # Remover diagonal
ggplot(matriz_melt, aes(x = Var1, y = Var2, fill = value)) +
geom_tile(color = "white") +
geom_text(aes(label = round(value, 2)), color = "black", size = 4) +
scale_fill_gradient(low = "white", high = "red", name = "Fuerza\nInteracción") +
labs(title = "Matriz de Fuerza de Interacciones entre Factores",
x = "Factor",
y = "Factor") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5))## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`geom_text()`).
# MÉTODO 5: GRÁFICAS DE INTERACCIÓN POR NIVELES MÚLTIPLES
cat("\n=== GRÁFICAS DE INTERACCIÓN AVANZADAS ===\n")##
## === GRÁFICAS DE INTERACCIÓN AVANZADAS ===# Gráfico de interacción triple (A vs B, facetado por C)
ggplot(design_df, aes(x = A, y = Y, color = B, group = B)) +
stat_summary(fun = mean, geom = "line", size = 1) +
stat_summary(fun = mean, geom = "point", size = 2) +
facet_wrap(~ C, labeller = label_both) +
labs(title = "Interacción A vs B por Niveles de C",
x = "Factor A",
y = "Rendimiento Promedio",
color = "Factor B") +
theme_minimal()
# ANÁLISIS CUANTITATIVO DE INTERACCIONES
cat("\n=== ANÁLISIS CUANTITATIVO DE INTERACCIONES ===\n")##
## === ANÁLISIS CUANTITATIVO DE INTERACCIONES ===# Modelo con todas las interacciones de 2 factores
modelo_completo <- lm(Y ~ (A + B + C + D + E)^2, data = design_df)
summary_interacciones <- summary(modelo_completo)
# Extraer coeficientes de interacción
coeficientes <- coef(modelo_completo)
interacciones <- coeficientes[grep(":", names(coeficientes))]
cat("Coeficientes de interacción (ordenados por magnitud):\n")## Coeficientes de interacción (ordenados por magnitud):interacciones_ordenadas <- sort(abs(interacciones), decreasing = TRUE)
print(interacciones_ordenadas)## A1:B1 D1:E1 A1:E1 B1:D1 C1:D1 B1:E1
## 3.937500e+00 5.625000e-01 5.000000e-01 3.750000e-01 3.750000e-01 2.500000e-01
## A1:C1 C1:E1 B1:C1 A1:D1
## 1.875000e-01 1.250000e-01 6.250000e-02 7.016351e-16# Identificar interacciones significativas
p_valores <- summary_interacciones$coefficients[,4]
interacciones_significativas <- p_valores[grep(":", names(p_valores))] < 0.05
cat("\nInteracciones significativas (p < 0.05):\n")##
## Interacciones significativas (p < 0.05):print(names(interacciones_significativas)[interacciones_significativas])## [1] "A1:B1"# GRÁFICO FINAL RESUMEN
cat("\n=== GRÁFICO RESUMEN DE INTERACCIONES MÁS FUERTES ===\n")##
## === GRÁFICO RESUMEN DE INTERACCIONES MÁS FUERTES ===# Tomar las 4 interacciones más fuertes
interacciones_fuertes <- names(interacciones_ordenadas)[1:4]
graficos_fuertes <- list()
for (i in 1:4) {
factores_inter <- strsplit(interacciones_fuertes[i], ":")[[1]]
graficos_fuertes[[i]] <- crear_grafico_interaccion(factores_inter[1], factores_inter[2]) +
labs(subtitle = paste("Fuerza:", round(interacciones_ordenadas[i], 3)))
}
# RESUMEN FINAL
cat("\n=== RESUMEN EJECUTIVO ===\n")##
## === RESUMEN EJECUTIVO ===cat("Total de interacciones de 2 factores analizadas:", length(combinaciones), "\n")## Total de interacciones de 2 factores analizadas: 10cat("Interacción más fuerte:", names(interacciones_ordenadas)[1],
"(", round(interacciones_ordenadas[1], 3), ")\n")## Interacción más fuerte: A1:B1 ( 3.937 )cat("Interacciones significativas (p < 0.05):", sum(interacciones_significativas), "\n")## Interacciones significativas (p < 0.05): 1cat("Rango de rendimiento:", range(Y), "\n")## Rango de rendimiento: 6 65🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹
🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹🔹