# ================================
# EJERCICIO 1 – DBCA
# ================================
# Datos del ejercicio
tratamiento <- factor(rep(c("A","B","C"), each = 4))
bloque <- factor(rep(c("1","2","3","4"), times = 3))
respuesta <- c(3,4,2,6, # A
7,9,3,10, # B
4,6,3,7) # C
datos <- data.frame(tratamiento, bloque, respuesta)
datos
## tratamiento bloque respuesta
## 1 A 1 3
## 2 A 2 4
## 3 A 3 2
## 4 A 4 6
## 5 B 1 7
## 6 B 2 9
## 7 B 3 3
## 8 B 4 10
## 9 C 1 4
## 10 C 2 6
## 11 C 3 3
## 12 C 4 7
# ================================
# MODELO ANOVA DBCA
# ================================
modelo <- aov(respuesta ~ tratamiento + bloque, datos)
summary(modelo)
## Df Sum Sq Mean Sq F value Pr(>F)
## tratamiento 2 25.17 12.583 13.73 0.00577 **
## bloque 3 42.00 14.000 15.27 0.00324 **
## Residuals 6 5.50 0.917
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ================================
# GRÁFICAS
# ================================
# Gráfica de cajas por tratamiento
boxplot(respuesta ~ tratamiento,
main = "Diagrama de cajas por tratamiento",
xlab = "Tratamiento", ylab = "Respuesta",
col = c("skyblue","salmon","lightgreen"))
# Gráfica de cajas por bloque
boxplot(respuesta ~ bloque,
main = "Diagrama de cajas por bloque",
xlab = "Bloque", ylab = "Respuesta",
col = rainbow(4))
# Gráfica de interacción
interaction.plot(bloque, tratamiento, respuesta,
main = "Gráfica de interacción",
xlab = "Bloque", ylab = "Respuesta",
col = 1:3, lwd = 2)
# ================================
# COMPARACIONES MÚLTIPLES
# ================================
TukeyHSD(modelo, "tratamiento")
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = respuesta ~ tratamiento + bloque, data = datos)
##
## $tratamiento
## diff lwr upr p adj
## B-A 3.50 1.4227684 5.5772316 0.0049834
## C-A 1.25 -0.8272316 3.3272316 0.2339241
## C-B -2.25 -4.3272316 -0.1727684 0.0366516
# ================================
# SUPUESTOS DEL MODELO
# ================================
# Normalidad de residuos
res <- residuals(modelo)
shapiro.test(res)
##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.93679, p-value = 0.4576
qqnorm(res)
qqline(res)
# Homogeneidad de varianzas
plot(modelo, which = 1)
# ==========================================
# EJERCICIO 2 – DISEÑO COMPLETAMENTE AL AZAR
# ==========================================
# Datos del ejercicio 2
tratamiento <- factor(rep(c("A","B","C","D","E"), each = 5))
respuesta <- c(
4.0, 4.0, 5.0, 0.5, 3.0, # A
5.0, 6.0, 2.0, 4.0, 4.0, # B
4.5, 4.0, 3.5, 2.0, 3.0, # C
2.5, 4.0, 6.5, 4.5, 4.0, # D
4.0, 4.0, 3.5, 2.0, 4.0 # E
)
datos2 <- data.frame(tratamiento, respuesta)
datos2
## tratamiento respuesta
## 1 A 4.0
## 2 A 4.0
## 3 A 5.0
## 4 A 0.5
## 5 A 3.0
## 6 B 5.0
## 7 B 6.0
## 8 B 2.0
## 9 B 4.0
## 10 B 4.0
## 11 C 4.5
## 12 C 4.0
## 13 C 3.5
## 14 C 2.0
## 15 C 3.0
## 16 D 2.5
## 17 D 4.0
## 18 D 6.5
## 19 D 4.5
## 20 D 4.0
## 21 E 4.0
## 22 E 4.0
## 23 E 3.5
## 24 E 2.0
## 25 E 4.0
# ==========================================
# ANOVA – DCA
# ==========================================
modelo2 <- aov(respuesta ~ tratamiento, data = datos2)
summary(modelo2)
## Df Sum Sq Mean Sq F value Pr(>F)
## tratamiento 4 4.46 1.115 0.626 0.649
## Residuals 20 35.60 1.780
# ==========================================
# GRÁFICAS
# ==========================================
# Diagrama de cajas por tratamiento
boxplot(respuesta ~ tratamiento,
main = "Boxplot por tratamiento (Ejercicio 2)",
xlab = "Tratamiento", ylab = "Respuesta",
col = rainbow(5))
# Medias por tratamiento
plot(tapply(respuesta, tratamiento, mean),
type = "b", pch = 19, lwd = 2,
main = "Medias por tratamiento",
xlab = "Tratamiento", ylab = "Media",
xaxt = "n")
axis(1, at = 1:5, labels = c("A","B","C","D","E"))
# ==========================================
# COMPARACIONES MÚLTIPLES (Tukey)
# ==========================================
TukeyHSD(modelo2, "tratamiento")
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = respuesta ~ tratamiento, data = datos2)
##
## $tratamiento
## diff lwr upr p adj
## B-A 0.9 -1.624968 3.424968 0.8211096
## C-A 0.1 -2.424968 2.624968 0.9999512
## D-A 1.0 -1.524968 3.524968 0.7595518
## E-A 0.2 -2.324968 2.724968 0.9992377
## C-B -0.8 -3.324968 1.724968 0.8745474
## D-B 0.1 -2.424968 2.624968 0.9999512
## E-B -0.7 -3.224968 1.824968 0.9182606
## D-C 0.9 -1.624968 3.424968 0.8211096
## E-C 0.1 -2.424968 2.624968 0.9999512
## E-D -0.8 -3.324968 1.724968 0.8745474
# ==========================================
# SUPUESTOS
# ==========================================
# Normalidad
res2 <- residuals(modelo2)
shapiro.test(res2)
##
## Shapiro-Wilk normality test
##
## data: res2
## W = 0.95716, p-value = 0.3608
qqnorm(res2)
qqline(res2)
# Homogeneidad de varianzas
plot(modelo2, which = 1)
# ==========================================
# EJERCICIO 3 – DISEÑO FACTORIAL 3x2 CON RÉPLICAS
# ==========================================
# Datos del Ejercicio 3
operador <- factor(rep(c(1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
2,2,2,2,2,2,2,2,
2,2,2,2,2,2,2,2,
3,3,3,3,3,3,3,3,
3,3,3,3,3,3,3,3), times = 1))
equipo <- factor(rep(c(1,1,1,1,1,1,1,1,
2,2,2,2,2,2,2,2), times = 3))
respuesta <- c(
# Operador 1 – Equipo 1
1.328, 0.985, 1.316, 1.553, 1.310, 1.113, 1.057, 1.144,
# Operador 1 – Equipo 2
1.273, 0.985, 1.134, 1.412, 0.917, 0.789, 0.671, 0.554,
# Operador 2 – Equipo 1
1.269, 1.268, 1.091, 1.195, 1.380, 1.093, 0.984, 1.087,
# Operador 2 – Equipo 2
1.036, 0.783, 1.108, 1.129, 1.132, 0.201, 0.900, 0.916,
# Operador 3 – Equipo 1
1.440, 1.079, 1.389, 1.611, 1.445, 1.150, 1.190, 1.247,
# Operador 3 – Equipo 2
1.454, 1.063, 1.219, 1.602, 1.583, 1.018, 1.050, 0.997
)
datos3 <- data.frame(operador, equipo, respuesta)
datos3
## operador equipo respuesta
## 1 1 1 1.328
## 2 1 1 0.985
## 3 1 1 1.316
## 4 1 1 1.553
## 5 1 1 1.310
## 6 1 1 1.113
## 7 1 1 1.057
## 8 1 1 1.144
## 9 1 2 1.273
## 10 1 2 0.985
## 11 1 2 1.134
## 12 1 2 1.412
## 13 1 2 0.917
## 14 1 2 0.789
## 15 1 2 0.671
## 16 1 2 0.554
## 17 2 1 1.269
## 18 2 1 1.268
## 19 2 1 1.091
## 20 2 1 1.195
## 21 2 1 1.380
## 22 2 1 1.093
## 23 2 1 0.984
## 24 2 1 1.087
## 25 2 2 1.036
## 26 2 2 0.783
## 27 2 2 1.108
## 28 2 2 1.129
## 29 2 2 1.132
## 30 2 2 0.201
## 31 2 2 0.900
## 32 2 2 0.916
## 33 3 1 1.440
## 34 3 1 1.079
## 35 3 1 1.389
## 36 3 1 1.611
## 37 3 1 1.445
## 38 3 1 1.150
## 39 3 1 1.190
## 40 3 1 1.247
## 41 3 2 1.454
## 42 3 2 1.063
## 43 3 2 1.219
## 44 3 2 1.602
## 45 3 2 1.583
## 46 3 2 1.018
## 47 3 2 1.050
## 48 3 2 0.997
# ==========================================
# MODELO FACTORIAL 3×2
# ==========================================
modelo3 <- aov(respuesta ~ operador * equipo, data = datos3)
summary(modelo3)
## Df Sum Sq Mean Sq F value Pr(>F)
## operador 2 0.5341 0.2670 4.797 0.01329 *
## equipo 1 0.4796 0.4796 8.615 0.00538 **
## operador:equipo 2 0.1006 0.0503 0.903 0.41300
## Residuals 42 2.3380 0.0557
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ==========================================
# GRÁFICAS
# ==========================================
# Boxplot por operador
boxplot(respuesta ~ operador,
main = "Boxplot por Operador",
xlab = "Operador", ylab = "Respuesta",
col = c("skyblue","salmon","lightgreen"))
# Boxplot por equipo
boxplot(respuesta ~ equipo,
main = "Boxplot por Equipo",
xlab = "Equipo", ylab = "Respuesta",
col = c("orange","lightblue"))
# Gráfica de interacción
interaction.plot(operador, equipo, respuesta,
main = "Interacción Operador × Equipo",
xlab = "Operador", ylab = "Respuesta",
trace.label = "Equipo",
col = c("red","blue"), lwd = 2)
# ==========================================
# SUPUESTOS DEL MODELO
# ==========================================
# Normalidad de residuos
res3 <- residuals(modelo3)
shapiro.test(res3)
##
## Shapiro-Wilk normality test
##
## data: res3
## W = 0.97452, p-value = 0.3759
qqnorm(res3)
qqline(res3)
# Homogeneidad de varianzas
plot(modelo3, which = 1)
# ==========================================
# COMPARACIONES MÚLTIPLES
# ==========================================
# Si los factores principales son significativos:
TukeyHSD(modelo3, "operador")
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = respuesta ~ operador * equipo, data = datos3)
##
## $operador
## diff lwr upr p adj
## 2-1 -0.0605625 -0.26322386 0.1420989 0.7495455
## 3-1 0.1872500 -0.01541136 0.3899114 0.0752937
## 3-2 0.2478125 0.04515114 0.4504739 0.0133037
TukeyHSD(modelo3, "equipo")
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = respuesta ~ operador * equipo, data = datos3)
##
## $equipo
## diff lwr upr p adj
## 2-1 -0.1999167 -0.3373678 -0.06246555 0.005385
Datos para termómetros Rtd y Mer
dia <- factor(rep(1:5, 5))
silo <- factor(rep(LETTERS[1:5], each = 5))
Rtd <- c(4.1, 4.3, 4.5, 4.2, 4.4, # A
4.3, 4.5, 4.7, 4.4, 4.6, # B
4.0, 4.2, 4.4, 4.1, 4.3, # C
4.2, 4.4, 4.6, 4.3, 4.5, # D
4.4, 4.6, 4.8, 4.5, 4.7) # E
Mer <- c(4.0, 4.2, 4.4, 4.1, 4.3, # A
4.2, 4.4, 4.6, 4.3, 4.5, # B
3.9, 4.1, 4.3, 4.0, 4.2, # C
4.1, 4.3, 4.5, 4.2, 4.4, # D
4.3, 4.5, 4.7, 4.4, 4.6) # E
datos4 <- data.frame(silo, dia, Rtd, Mer)
# a) Conjetura inicial
cat("a) Conjetura inicial:\n")
## a) Conjetura inicial:
cat(" Las mediciones Rtd parecen sistemáticamente mayores (0.1°C aprox.)\n")
## Las mediciones Rtd parecen sistemáticamente mayores (0.1°C aprox.)
cat(" que las mediciones Mer, lo que sugiere un posible sesgo.\n\n")
## que las mediciones Mer, lo que sugiere un posible sesgo.
# b) ANOVA solo para Rtd
cat("b) Análisis para termómetro Rtd:\n")
## b) Análisis para termómetro Rtd:
modelo_rtd <- aov(Rtd ~ silo + dia, data = datos4)
anova_rtd <- anova(modelo_rtd)
print(anova_rtd)
## Analysis of Variance Table
##
## Response: Rtd
## Df Sum Sq Mean Sq F value Pr(>F)
## silo 4 0.5 0.125 6.2397e+29 < 2.2e-16 ***
## dia 4 0.5 0.125 6.2397e+29 < 2.2e-16 ***
## Residuals 16 0.0 0.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nConclusiones para Rtd:\n")
##
## Conclusiones para Rtd:
if(anova_rtd["silo", "Pr(>F)"] < 0.05) {
cat("- Diferencias significativas entre silos\n")
} else {
cat("- No hay diferencias significativas entre silos\n")
}
## - Diferencias significativas entre silos
# c) ANOVA solo para Mer
cat("\nc) Análisis para termómetro Mer:\n")
##
## c) Análisis para termómetro Mer:
modelo_mer <- aov(Mer ~ silo + dia, data = datos4)
anova_mer <- anova(modelo_mer)
print(anova_mer)
## Analysis of Variance Table
##
## Response: Mer
## Df Sum Sq Mean Sq F value Pr(>F)
## silo 4 0.5 0.125 1.2016e+29 < 2.2e-16 ***
## dia 4 0.5 0.125 1.2016e+29 < 2.2e-16 ***
## Residuals 16 0.0 0.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nConclusiones para Mer:\n")
##
## Conclusiones para Mer:
if(anova_mer["silo", "Pr(>F)"] < 0.05) {
cat("- Diferencias significativas entre silos\n")
} else {
cat("- No hay diferencias significativas entre silos\n")
}
## - Diferencias significativas entre silos
# d) Comparación de conclusiones
cat("\nd) Comparación de conclusiones:\n")
##
## d) Comparación de conclusiones:
cat(" Ambos termómetros muestran el mismo patrón:\n")
## Ambos termómetros muestran el mismo patrón:
if((anova_rtd["silo", "Pr(>F)"] < 0.05) == (anova_mer["silo", "Pr(>F)"] < 0.05)) {
cat(" Conclusiones coincidentes en cuanto a diferencias entre silos\n")
} else {
cat(" Conclusiones diferentes\n")
}
## Conclusiones coincidentes en cuanto a diferencias entre silos
# e) Datos pareados
cat("\ne) Análisis de datos pareados (diferencias):\n")
##
## e) Análisis de datos pareados (diferencias):
datos4$diferencia <- datos4$Rtd - datos4$Mer
modelo_diff <- aov(diferencia ~ silo + dia, data = datos4)
anova_diff <- anova(modelo_diff)
print(anova_diff)
## Analysis of Variance Table
##
## Response: diferencia
## Df Sum Sq Mean Sq F value Pr(>F)
## silo 4 3.270e-32 8.166e-33 0.0286 0.9982
## dia 4 2.770e-32 6.933e-33 0.0243 0.9987
## Residuals 16 4.568e-30 2.855e-31
cat("\nConclusiones sobre diferencias:\n")
##
## Conclusiones sobre diferencias:
if(anova_diff["silo", "Pr(>F)"] < 0.05) {
cat(" Hay diferencias significativas en las discrepancias entre silos\n")
} else {
cat(" No hay diferencias significativas en las discrepancias entre silos\n")
}
## No hay diferencias significativas en las discrepancias entre silos
cat("\nMedia de la diferencia:", round(mean(datos4$diferencia), 4))
##
## Media de la diferencia: 0.1
if(mean(datos4$diferencia) > 0) {
cat("El termómetro Rtd mide consistentemente más alto que el Mer\n")
}
## El termómetro Rtd mide consistentemente más alto que el Mer
EJERCICIO 5. Datos del cuadro latino
lote <- factor(rep(1:5, each = 5))
dia <- factor(rep(1:5, 5))
catalizador <- 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")
tiempo <- c(8, 11, 9, 10, 12,
7, 10, 8, 9, 11,
9, 12, 10, 11, 13,
8, 11, 9, 10, 12,
10, 13, 11, 12, 14)
datos5 <- data.frame(lote, dia, catalizador, tiempo)
# a) Aleatorización
cat("a) Aleatorización:\n")
## a) Aleatorización:
cat(" Se aleatorizó el orden de aplicación de catalizadores\n")
## Se aleatorizó el orden de aplicación de catalizadores
cat(" dentro de cada fila y columna del cuadro latino\n\n")
## dentro de cada fila y columna del cuadro latino
# b) Modelo e hipótesis
cat("b) Modelo: Y_ijk = μ + α_i + τ_j + β_k + ε_ijk\n")
## b) Modelo: Y_ijk = μ + α_i + τ_j + β_k + ε_ijk
cat(" donde:\n")
## donde:
cat(" α_i: efecto del lote i\n")
## α_i: efecto del lote i
cat(" τ_j: efecto del catalizador j\n")
## τ_j: efecto del catalizador j
cat(" β_k: efecto del día k\n\n")
## β_k: efecto del día k
cat(" Hipótesis para tratamientos:\n")
## Hipótesis para tratamientos:
cat(" H0: τ_A = τ_B = τ_C = τ_D = τ_E = 0\n")
## H0: τ_A = τ_B = τ_C = τ_D = τ_E = 0
cat(" H1: Al menos un τ_j ≠ 0\n\n")
## H1: Al menos un τ_j ≠ 0
# c) Análisis ANOVA
cat("c) Análisis ANOVA:\n")
## c) Análisis ANOVA:
modelo5 <- aov(tiempo ~ lote + dia + catalizador, data = datos5)
anova5 <- anova(modelo5)
print(anova5)
## Analysis of Variance Table
##
## Response: tiempo
## Df Sum Sq Mean Sq F value Pr(>F)
## lote 4 26 6.5 3.0439e+30 <2e-16 ***
## dia 4 50 12.5 5.8537e+30 <2e-16 ***
## catalizador 4 0 0.0 2.0213e+00 0.1553
## Residuals 12 0 0.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Prueba de diferencias entre catalizadores
cat("\nDiferencias entre catalizadores:\n")
##
## Diferencias entre catalizadores:
if(anova5["catalizador", "Pr(>F)"] < 0.05) {
cat("Existen diferencias significativas entre catalizadores (p =",
round(anova5["catalizador", "Pr(>F)"], 4), ")\n")
# Comparación múltiple
tukey5 <- TukeyHSD(modelo5, "catalizador")
cat("\nComparaciones múltiples (Tukey):\n")
print(tukey5$catalizador)
} else {
cat("No existen diferencias significativas entre catalizadores (p =",
round(anova5["catalizador", "Pr(>F)"], 4), ")\n")
}
## No existen diferencias significativas entre catalizadores (p = 0.1553 )
# d) Factores de ruido
cat("\nd) Efecto de factores de ruido:\n")
##
## d) Efecto de factores de ruido:
if(anova5["lote", "Pr(>F)"] < 0.05) {
cat("Los lotes afectan significativamente el tiempo de reacción (p =",
round(anova5["lote", "Pr(>F)"], 4), ")\n")
} else {
cat("Los lotes no afectan significativamente el tiempo de reacción\n")
}
## Los lotes afectan significativamente el tiempo de reacción (p = 0 )
if(anova5["dia", "Pr(>F)"] < 0.05) {
cat("Los días afectan significativamente el tiempo de reacción (p =",
round(anova5["dia", "Pr(>F)"], 4), ")\n")
} else {
cat("Los días no afectan significativamente el tiempo de reacción\n")
}
## Los días afectan significativamente el tiempo de reacción (p = 0 )
# e) Gráficas de medias
cat("\ne) Gráficas de medias:\n")
##
## e) Gráficas de medias:
# Medias por catalizador
medias_cat <- aggregate(tiempo ~ catalizador, datos5, mean)
p1 <- ggplot(medias_cat, aes(x = catalizador, y = tiempo, fill = catalizador)) +
geom_bar(stat = "identity") +
labs(title = "Tiempo promedio por catalizador", x = "Catalizador", y = "Tiempo") +
theme_minimal()
# Medias por lote
medias_lote <- aggregate(tiempo ~ lote, datos5, mean)
p2 <- ggplot(medias_lote, aes(x = lote, y = tiempo, fill = lote)) +
geom_bar(stat = "identity") +
labs(title = "Tiempo promedio por lote", x = "Lote", y = "Tiempo") +
theme_minimal()
# Medias por día
medias_dia <- aggregate(tiempo ~ dia, datos5, mean)
p3 <- ggplot(medias_dia, aes(x = dia, y = tiempo, fill = dia)) +
geom_bar(stat = "identity") +
labs(title = "Tiempo promedio por día", x = "Día", y = "Tiempo") +
theme_minimal()
grid.arrange(p1, p2, p3, ncol = 3)
# Mejor tratamiento (menor tiempo)
cat("\nMejor catalizador (menor tiempo promedio):",
medias_cat$catalizador[which.min(medias_cat$tiempo)],
"con tiempo =", min(medias_cat$tiempo), "\n")
##
## Mejor catalizador (menor tiempo promedio): A con tiempo = 10.4
# f) Verificación de supuestos
cat("\nf) Verificación de supuestos:\n")
##
## f) Verificación de supuestos:
par(mfrow = c(2, 2))
plot(modelo5)
EJERCICIO 6.Eliminar factor día del ejercicio 5
cat("6. Eliminación del factor día\n")
## 6. Eliminación del factor día
# a) Justificación de eliminación
cat("a) Justificación de eliminación:\n")
## a) Justificación de eliminación:
cat(" Si el factor día no fue significativo (p > 0.05) en el ejercicio 5,\n")
## Si el factor día no fue significativo (p > 0.05) en el ejercicio 5,
cat(" entonces su eliminación está justificada.\n\n")
## entonces su eliminación está justificada.
# b) Diseño y modelo
cat("b) Diseño: Diseño en Bloques Completos Aleatorizados (DBCA)\n")
## b) Diseño: Diseño en Bloques Completos Aleatorizados (DBCA)
cat(" Modelo: Y_ij = μ + τ_i + β_j + ε_ij\n")
## Modelo: Y_ij = μ + τ_i + β_j + ε_ij
cat(" donde τ_i: efecto del catalizador, β_j: efecto del lote\n\n")
## donde τ_i: efecto del catalizador, β_j: efecto del lote
# c) Prueba de hipótesis
cat("c) Análisis ANOVA sin factor día:\n")
## c) Análisis ANOVA sin factor día:
modelo6 <- aov(tiempo ~ catalizador + lote, data = datos5)
anova6 <- anova(modelo6)
print(anova6)
## Analysis of Variance Table
##
## Response: tiempo
## Df Sum Sq Mean Sq F value Pr(>F)
## catalizador 4 0 0.000 0.00 1.0000
## lote 4 26 6.500 2.08 0.1311
## Residuals 16 50 3.125
cat("\nConclusiones:\n")
##
## Conclusiones:
if(anova6["catalizador", "Pr(>F)"] < 0.05) {
cat("Existen diferencias significativas entre catalizadores (p =",
round(anova6["catalizador", "Pr(>F)"], 4), ")\n")
} else {
cat("No existen diferencias significativas entre catalizadores (p =",
round(anova6["catalizador", "Pr(>F)"], 4), ")\n")
}
## No existen diferencias significativas entre catalizadores (p = 1 )
# d) Comparación de cuadro medio del error
cat("\nd) Comparación del cuadro medio del error (MSE):\n")
##
## d) Comparación del cuadro medio del error (MSE):
cat(" Ejercicio 5 (con día): MSE =", round(anova5["Residuals", "Mean Sq"], 4), "\n")
## Ejercicio 5 (con día): MSE = 0
cat(" Ejercicio 6 (sin día): MSE =", round(anova6["Residuals", "Mean Sq"], 4), "\n")
## Ejercicio 6 (sin día): MSE = 3.125
if(anova6["Residuals", "Mean Sq"] > anova5["Residuals", "Mean Sq"]) {
cat(" El MSE aumentó al eliminar el día, lo que indica que el día\n")
cat(" explicaba parte de la variabilidad no explicada.\n")
} else {
cat(" El MSE disminuyó o se mantuvo similar.\n")
}
## El MSE aumentó al eliminar el día, lo que indica que el día
## explicaba parte de la variabilidad no explicada.
# e) ¿Por qué mismas conclusiones?
cat("\ne) Razón de mismas conclusiones en tratamientos:\n")
##
## e) Razón de mismas conclusiones en tratamientos:
cat(" Si el factor día no era significativo, su eliminación no debería\n")
## Si el factor día no era significativo, su eliminación no debería
cat(" afectar considerablemente las conclusiones sobre los catalizadores.\n")
## afectar considerablemente las conclusiones sobre los catalizadores.
Ejercicio 7 - Eliminación adicional del factor lote
Eliminar ambos factores de bloque
cat("7. Eliminación adicional del factor lote\n")
## 7. Eliminación adicional del factor lote
# a) Justificación de eliminación
cat("a) Justificación de doble eliminación:\n")
## a) Justificación de doble eliminación:
cat(" Si ni día ni lote fueron significativos en ejercicios previos,\n")
## Si ni día ni lote fueron significativos en ejercicios previos,
cat(" entonces su eliminación simultánea está justificada.\n\n")
## entonces su eliminación simultánea está justificada.
# b) Diseño y modelo
cat("b) Diseño: Diseño Completamente Aleatorizado (DCA)\n")
## b) Diseño: Diseño Completamente Aleatorizado (DCA)
cat(" Modelo: Y_ij = μ + τ_i + ε_ij\n")
## Modelo: Y_ij = μ + τ_i + ε_ij
cat(" donde τ_i: efecto del catalizador\n\n")
## donde τ_i: efecto del catalizador
# c) Prueba de hipótesis
cat("c) Análisis ANOVA (DCA):\n")
## c) Análisis ANOVA (DCA):
modelo7 <- aov(tiempo ~ catalizador, data = datos5)
anova7 <- anova(modelo7)
print(anova7)
## Analysis of Variance Table
##
## Response: tiempo
## Df Sum Sq Mean Sq F value Pr(>F)
## catalizador 4 0 0.0 0 1
## Residuals 20 76 3.8
cat("\nConclusiones:\n")
##
## Conclusiones:
if(anova7["catalizador", "Pr(>F)"] < 0.05) {
cat("Existen diferencias significativas entre catalizadores (p =",
round(anova7["catalizador", "Pr(>F)"], 4), ")\n")
} else {
cat("No existen diferencias significativas entre catalizadores (p =",
round(anova7["catalizador", "Pr(>F)"], 4), ")\n")
}
## No existen diferencias significativas entre catalizadores (p = 1 )
# d) Comparación de MSE
cat("\nd) Comparación del cuadro medio del error (MSE):\n")
##
## d) Comparación del cuadro medio del error (MSE):
cat(" Ejercicio 5 (cuadro latino): MSE =", round(anova5["Residuals", "Mean Sq"], 4), "\n")
## Ejercicio 5 (cuadro latino): MSE = 0
cat(" Ejercicio 6 (DBCA): MSE =", round(anova6["Residuals", "Mean Sq"], 4), "\n")
## Ejercicio 6 (DBCA): MSE = 3.125
cat(" Ejercicio 7 (DCA): MSE =", round(anova7["Residuals", "Mean Sq"], 4), "\n")
## Ejercicio 7 (DCA): MSE = 3.8
# e) Razón de mismas conclusiones
cat("\ne) Razón de mismas conclusiones:\n")
##
## e) Razón de mismas conclusiones:
cat(" Si los factores de bloque no eran significativos, eliminarlos\n")
## Si los factores de bloque no eran significativos, eliminarlos
cat(" no cambia las conclusiones sobre los tratamientos principales.\n\n")
## no cambia las conclusiones sobre los tratamientos principales.
# f) Efecto si factores de bloque hubieran sido significativos
cat("f) Efecto si factores de bloque fueran significativos:\n")
## f) Efecto si factores de bloque fueran significativos:
cat(" Si los factores de bloque hubieran sido significativos,\n")
## Si los factores de bloque hubieran sido significativos,
cat(" eliminarlos aumentaría el error experimental y podría\n")
## eliminarlos aumentaría el error experimental y podría
cat(" enmascarar diferencias reales entre tratamientos.\n")
## enmascarar diferencias reales entre tratamientos.
Ejercicio 8 - Cuadro Latino con tres factores
# Ejercicio 8 - Cuadro Latino con tres factores
# Datos del cuadro latino con tres factores
proveedor <- factor(c("A", "B", "C", "B", "C", "A", "C", "A", "B"))
inspector <- factor(c("I", "I", "I", "II", "II", "II", "III", "III", "III"))
escala <- factor(c(1, 2, 3, 2, 3, 1, 3, 1, 2))
peso <- c(14.8, 15.1, 14.9, 15.2, 15.0, 14.7, 15.3, 14.8, 15.1)
datos8 <- data.frame(proveedor, inspector, escala, peso)
# Función segura para obtener valores de ANOVA
obtener_valor_anova <- function(anova_tab, factor_nombre) {
if(factor_nombre %in% rownames(anova_tab)) {
return(list(
existe = TRUE,
suma_cuadrados = anova_tab[factor_nombre, "Sum Sq"],
grados_libertad = anova_tab[factor_nombre, "Df"],
cuadro_medio = anova_tab[factor_nombre, "Mean Sq"],
valor_f = anova_tab[factor_nombre, "F value"],
p_valor = anova_tab[factor_nombre, "Pr(>F)"]
))
} else {
return(list(
existe = FALSE,
suma_cuadrados = NA,
grados_libertad = NA,
cuadro_medio = NA,
valor_f = NA,
p_valor = NA
))
}
}
# a) ¿Hay diferencias entre los proveedores?
cat("8a) Diferencias entre proveedores:\n")
## 8a) Diferencias entre proveedores:
modelo8 <- aov(peso ~ proveedor + inspector + escala, data = datos8)
anova8 <- anova(modelo8)
print(anova8)
## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## proveedor 2 0.228889 0.114444 6.4375 0.05619 .
## inspector 2 0.028889 0.014444 0.8125 0.50568
## Residuals 4 0.071111 0.017778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
info_proveedor <- obtener_valor_anova(anova8, "proveedor")
if(info_proveedor$existe) {
if(info_proveedor$p_valor < 0.05) {
cat("\nExisten diferencias significativas entre proveedores (p =",
round(info_proveedor$p_valor, 4), ")\n")
} else {
cat("\nNo existen diferencias significativas entre proveedores (p =",
round(info_proveedor$p_valor, 4), ")\n")
}
}
##
## No existen diferencias significativas entre proveedores (p = 0.0562 )
# b) ¿Hay diferencias entre inspectores y entre escalas?
cat("\n8b) Diferencias entre inspectores y entre escalas:\n")
##
## 8b) Diferencias entre inspectores y entre escalas:
info_inspector <- obtener_valor_anova(anova8, "inspector")
if(info_inspector$existe) {
if(info_inspector$p_valor < 0.05) {
cat("Existen diferencias significativas entre inspectores (p =",
round(info_inspector$p_valor, 4), ")\n")
} else {
cat("No existen diferencias significativas entre inspectores (p =",
round(info_inspector$p_valor, 4), ")\n")
}
} else {
cat("Factor 'inspector' no encontrado en ANOVA\n")
}
## No existen diferencias significativas entre inspectores (p = 0.5057 )
info_escala <- obtener_valor_anova(anova8, "escala")
if(info_escala$existe) {
if(info_escala$p_valor < 0.05) {
cat("Existen diferencias significativas entre escalas (p =",
round(info_escala$p_valor, 4), ")\n")
} else {
cat("No existen diferencias significativas entre escalas (p =",
round(info_escala$p_valor, 4), ")\n")
}
} else {
cat("Factor 'escala' no encontrado en ANOVA\n")
}
## Factor 'escala' no encontrado en ANOVA
# c) Si el peso debe ser 15 g, ¿cuál proveedor es mejor?
cat("\n8c) Mejor proveedor (más cercano a 15g):\n")
##
## 8c) Mejor proveedor (más cercano a 15g):
medias_prov <- aggregate(peso ~ proveedor, datos8, mean)
medias_prov$diferencia <- abs(medias_prov$peso - 15)
cat("Medias por proveedor:\n")
## Medias por proveedor:
print(medias_prov)
## proveedor peso diferencia
## 1 A 14.76667 0.23333333
## 2 B 15.13333 0.13333333
## 3 C 15.06667 0.06666667
mejor <- medias_prov[which.min(medias_prov$diferencia), ]
cat("\nMejor proveedor: ", mejor$proveedor,
" con peso promedio = ", round(mejor$peso, 4),
" (diferencia = ", round(mejor$diferencia, 4), ")\n", sep = "")
##
## Mejor proveedor: 3 con peso promedio = 15.0667 (diferencia = 0.0667)
# d) Si algún factor de bloque es no significativo, elimínelo y haga el análisis adecuado
cat("\n8d) Eliminación de factores no significativos:\n")
##
## 8d) Eliminación de factores no significativos:
# Determinar qué factores mantener
factores_a_mantener <- c("proveedor") # Factor de tratamiento siempre se mantiene
# Verificar inspector
if(info_inspector$existe) {
if(info_inspector$p_valor < 0.05) {
factores_a_mantener <- c(factores_a_mantener, "inspector")
cat("Inspector: MANTENER (significativo, p =",
round(info_inspector$p_valor, 4), ")\n")
} else {
cat("Inspector: ELIMINAR (no significativo, p =",
round(info_inspector$p_valor, 4), ")\n")
}
}
## Inspector: ELIMINAR (no significativo, p = 0.5057 )
# Verificar escala
if(info_escala$existe) {
if(info_escala$p_valor < 0.05) {
factores_a_mantener <- c(factores_a_mantener, "escala")
cat("Escala: MANTENER (significativo, p =",
round(info_escala$p_valor, 4), ")\n")
} else {
cat("Escala: ELIMINAR (no significativo, p =",
round(info_escala$p_valor, 4), ")\n")
}
}
# Crear fórmula para el modelo final
if(length(factores_a_mantener) == 1) {
formula_final <- "peso ~ proveedor"
nombre_diseno <- "Diseño Completamente Aleatorizado (DCA)"
} else if(length(factores_a_mantener) == 2) {
formula_final <- paste("peso ~", paste(factores_a_mantener, collapse = " + "))
nombre_diseno <- "Diseño en Bloques"
} else {
formula_final <- "peso ~ proveedor + inspector + escala"
nombre_diseno <- "Cuadro Latino"
}
cat("\nModelo final: ", nombre_diseno, "\n", sep = "")
##
## Modelo final: Diseño Completamente Aleatorizado (DCA)
cat("Fórmula: ", formula_final, "\n", sep = "")
## Fórmula: peso ~ proveedor
# Ajustar modelo final
modelo_final <- aov(as.formula(formula_final), data = datos8)
cat("\nAnálisis ANOVA del modelo final:\n")
##
## Análisis ANOVA del modelo final:
anova_final <- anova(modelo_final)
print(anova_final)
## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## proveedor 2 0.22889 0.114444 6.8667 0.02811 *
## Residuals 6 0.10000 0.016667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Comparación de modelos
cat("\nComparación de modelos:\n")
##
## Comparación de modelos:
cat("Modelo completo (Cuadro Latino):\n")
## Modelo completo (Cuadro Latino):
cat(" R² ajustado =", round(summary.lm(modelo8)$adj.r.squared, 4), "\n")
## R² ajustado = 0.5676
cat(" MSE =", round(anova8["Residuals", "Mean Sq"], 4), "\n\n")
## MSE = 0.0178
cat("Modelo final (", nombre_diseno, "):\n", sep = "")
## Modelo final (Diseño Completamente Aleatorizado (DCA)):
cat(" R² ajustado =", round(summary.lm(modelo_final)$adj.r.squared, 4), "\n")
## R² ajustado = 0.5946
cat(" MSE =", round(anova_final["Residuals", "Mean Sq"], 4), "\n")
## MSE = 0.0167
# Gráficas de diagnóstico del modelo final
cat("\nDiagnóstico del modelo final:\n")
##
## Diagnóstico del modelo final:
par(mfrow = c(2, 2))
plot(modelo_final, main = paste("Diagnóstico -", nombre_diseno))
# Si solo queda el factor proveedor, hacer comparaciones múltiples
if(length(factores_a_mantener) == 1 ||
(length(factores_a_mantener) == 2 &&
!"inspector" %in% factores_a_mantener &&
!"escala" %in% factores_a_mantener)) {
cat("\nComparaciones múltiples de proveedores (Tukey):\n")
tukey_proveedores <- TukeyHSD(modelo_final, "proveedor")
print(tukey_proveedores$proveedor)
# Gráfica de medias
medias_plot <- ggplot(medias_prov, aes(x = proveedor, y = peso, fill = proveedor)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 15, linetype = "dashed", color = "red", size = 1) +
geom_text(aes(label = round(peso, 2)), vjust = -0.5, size = 4) +
labs(title = "Peso promedio por proveedor",
subtitle = "Línea roja = objetivo de 15g",
x = "Proveedor", y = "Peso promedio (g)") +
ylim(0, max(medias_prov$peso) * 1.2) +
theme_minimal() +
theme(legend.position = "none")
print(medias_plot)
}
##
## Comparaciones múltiples de proveedores (Tukey):
## diff lwr upr p adj
## B-A 0.36666667 0.04324214 0.6900912 0.03046498
## C-A 0.30000000 -0.02342452 0.6234245 0.06589643
## C-B -0.06666667 -0.39009119 0.2567579 0.80848888
# Verificación de supuestos del modelo final
cat("\nVerificación de supuestos del modelo final:\n")
##
## Verificación de supuestos del modelo final:
residuos_final <- residuals(modelo_final)
# Normalidad
cat("1. Prueba de normalidad (Shapiro-Wilk):\n")
## 1. Prueba de normalidad (Shapiro-Wilk):
shapiro_final <- shapiro.test(residuos_final)
cat(" p-value =", round(shapiro_final$p.value, 4))
## p-value = 0.4756
if(shapiro_final$p.value > 0.05) {
cat(" → No se rechaza normalidad ✓\n")
} else {
cat(" → Se rechaza normalidad ✗\n")
}
## → No se rechaza normalidad ✓
# Homogeneidad de varianzas (si hay más de un proveedor)
if(length(unique(datos8$proveedor)) > 1) {
cat("\n2. Prueba de homogeneidad de varianzas (Levene):\n")
levene_final <- tryCatch({
leveneTest(peso ~ proveedor, data = datos8)
}, error = function(e) {
return(NULL)
})
if(!is.null(levene_final)) {
cat(" p-value =", round(levene_final$`Pr(>F)`[1], 4))
if(levene_final$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas ✓\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas ✗\n")
}
}
}
##
## 2. Prueba de homogeneidad de varianzas (Levene):
## p-value = 0.4219 → No se rechaza homogeneidad de varianzas ✓
# Conclusión final
cat("\n" , strrep("=", 50), "\n", sep = "")
##
## ==================================================
cat("CONCLUSIÓN FINAL EJERCICIO 8:\n")
## CONCLUSIÓN FINAL EJERCICIO 8:
cat(strrep("=", 50), "\n")
## ==================================================
cat("1. Proveedor recomendado: ", mejor$proveedor, "\n", sep = "")
## 1. Proveedor recomendado: 3
cat(" - Peso promedio: ", round(mejor$peso, 3), "g\n", sep = "")
## - Peso promedio: 15.067g
cat(" - Diferencia respecto al objetivo (15g): ",
round(mejor$diferencia, 3), "g\n", sep = "")
## - Diferencia respecto al objetivo (15g): 0.067g
cat(" - Desempeño: ",
ifelse(mejor$diferencia < 0.1, "Excelente",
ifelse(mejor$diferencia < 0.2, "Bueno", "Aceptable")), "\n\n", sep = "")
## - Desempeño: Excelente
cat("2. Diseño experimental seleccionado: ", nombre_diseno, "\n", sep = "")
## 2. Diseño experimental seleccionado: Diseño Completamente Aleatorizado (DCA)
cat(" - Factores incluidos: ", paste(factores_a_mantener, collapse = ", "), "\n", sep = "")
## - Factores incluidos: proveedor
cat(" - R² ajustado del modelo: ",
round(summary.lm(modelo_final)$adj.r.squared, 4), "\n\n", sep = "")
## - R² ajustado del modelo: 0.5946
cat("3. Recomendaciones:\n")
## 3. Recomendaciones:
cat(" a) Seleccionar al proveedor ", mejor$proveedor, " como principal\n", sep = "")
## a) Seleccionar al proveedor 3 como principal
cat(" b) ",
ifelse(info_proveedor$existe && info_proveedor$p_valor < 0.05,
"Las diferencias entre proveedores son estadísticamente significativas",
"No hay evidencia estadística de diferencias entre proveedores"), "\n", sep = "")
## b) No hay evidencia estadística de diferencias entre proveedores
cat(" c) Considerar otros factores como costo y disponibilidad\n")
## c) Considerar otros factores como costo y disponibilidad
cat(strrep("=", 50), "\n")
## ==================================================
Ejercicio 9 - Cuadro Grecolatino
# Datos del cuadro grecolatino
semana <- factor(rep(1:2, each = 16))
dia_semana <- factor(rep(rep(1:4, each = 4), 2))
chofer <- factor(rep(c("C1", "C2", "C3", "C4"), 8))
marca <- factor(c("M1", "M2", "M3", "M4",
"M2", "M3", "M4", "M1",
"M3", "M4", "M1", "M2",
"M4", "M1", "M2", "M3",
"M1", "M2", "M3", "M4",
"M2", "M3", "M4", "M1",
"M3", "M4", "M1", "M2",
"M4", "M1", "M2", "M3"))
ruta <- factor(c("A", "B", "C", "D",
"B", "C", "D", "A",
"C", "D", "A", "B",
"D", "A", "B", "C",
"A", "B", "C", "D",
"B", "C", "D", "A",
"C", "D", "A", "B",
"D", "A", "B", "C"))
costo <- c(150, 155, 160, 165,
152, 157, 162, 167,
153, 158, 163, 168,
151, 156, 161, 166,
149, 154, 159, 164,
150, 155, 160, 165,
152, 157, 162, 167,
148, 153, 158, 163)
datos9 <- data.frame(semana, dia_semana, chofer, marca, ruta, costo)
# a) Análisis de varianza
cat("9a) Análisis de varianza del cuadro grecolatino:\n")
## 9a) Análisis de varianza del cuadro grecolatino:
modelo9 <- aov(costo ~ semana + dia_semana + chofer + marca + ruta, data = datos9)
anova9 <- anova(modelo9)
print(anova9)
## Analysis of Variance Table
##
## Response: costo
## Df Sum Sq Mean Sq F value Pr(>F)
## semana 1 24.5 24.50 93.546 3.464e-09 ***
## dia_semana 3 49.5 16.50 63.000 1.141e-10 ***
## chofer 3 1000.0 333.33 1272.727 < 2.2e-16 ***
## marca 3 0.0 0.00 0.000 1
## Residuals 21 5.5 0.26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# b) Comparaciones múltiples para factores significativos
cat("\n9b) Pruebas de comparaciones múltiples:\n")
##
## 9b) Pruebas de comparaciones múltiples:
# Función para extraer p-value de forma segura
get_p_value <- function(anova_table, factor_name) {
if(factor_name %in% rownames(anova_table)) {
return(anova_table[factor_name, "Pr(>F)"])
} else {
return(NA)
}
}
# Solo para factores significativos
p_ruta <- get_p_value(anova9, "ruta")
if(!is.na(p_ruta) && p_ruta < 0.05) {
cat("\nComparaciones para rutas (Tukey):\n")
tukey_ruta <- TukeyHSD(modelo9, "ruta")
print(tukey_ruta$ruta)
}
p_chofer <- get_p_value(anova9, "chofer")
if(!is.na(p_chofer) && p_chofer < 0.05) {
cat("\nComparaciones para choferes (Tukey):\n")
tukey_chofer <- TukeyHSD(modelo9, "chofer")
print(tukey_chofer$chofer)
}
##
## Comparaciones para choferes (Tukey):
## diff lwr upr p adj
## C2-C1 5 4.28677 5.71323 7.382983e-14
## C3-C1 10 9.28677 10.71323 4.363176e-14
## C4-C1 15 14.28677 15.71323 4.363176e-14
## C3-C2 5 4.28677 5.71323 7.382983e-14
## C4-C2 10 9.28677 10.71323 4.363176e-14
## C4-C3 5 4.28677 5.71323 7.382983e-14
p_marca <- get_p_value(anova9, "marca")
if(!is.na(p_marca) && p_marca < 0.05) {
cat("\nComparaciones para marcas (Tukey):\n")
tukey_marca <- TukeyHSD(modelo9, "marca")
print(tukey_marca$marca)
}
# c) Gráficas de medias y diagramas de dispersión
cat("\n9c) Representación gráfica:\n")
##
## 9c) Representación gráfica:
# Gráficas de medias
p1 <- ggplot(datos9, aes(x = ruta, y = costo, fill = ruta)) +
geom_boxplot() +
labs(title = "Costos por Ruta", x = "Ruta", y = "Costo") +
theme_minimal()
p2 <- ggplot(datos9, aes(x = chofer, y = costo, fill = chofer)) +
geom_boxplot() +
labs(title = "Costos por Chofer", x = "Chofer", y = "Costo") +
theme_minimal()
p3 <- ggplot(datos9, aes(x = marca, y = costo, fill = marca)) +
geom_boxplot() +
labs(title = "Costos por Marca", x = "Marca", y = "Costo") +
theme_minimal()
p4 <- ggplot(datos9, aes(x = dia_semana, y = costo, fill = dia_semana)) +
geom_boxplot() +
labs(title = "Costos por Día", x = "Día de la semana", y = "Costo") +
theme_minimal()
grid.arrange(p1, p2, p3, p4, ncol = 2)
# Diagramas de dispersión
p5 <- ggplot(datos9, aes(x = ruta, y = costo)) +
geom_jitter(width = 0.2, aes(color = chofer), size = 3) +
labs(title = "Dispersión: Costo vs Ruta (color por chofer)") +
theme_minimal()
p6 <- ggplot(datos9, aes(x = ruta, y = costo)) +
geom_jitter(width = 0.2, aes(color = marca), size = 3) +
labs(title = "Dispersión: Costo vs Ruta (color por marca)") +
theme_minimal()
grid.arrange(p5, p6, ncol = 2)
# d) Mejor y peor ruta
cat("\n9d) Mejor y peor ruta:\n")
##
## 9d) Mejor y peor ruta:
medias_rutas <- aggregate(costo ~ ruta, datos9, mean)
mejor_ruta <- medias_rutas[which.min(medias_rutas$costo), ]
peor_ruta <- medias_rutas[which.max(medias_rutas$costo), ]
cat("Mejor ruta:", mejor_ruta$ruta, "(costo promedio =", round(mejor_ruta$costo, 2), ")\n")
## Mejor ruta: 1 (costo promedio = 158.12 )
cat("Peor ruta:", peor_ruta$ruta, "(costo promedio =", round(peor_ruta$costo, 2), ")\n")
## Peor ruta: 1 (costo promedio = 158.12 )
# e) Diferencias entre choferes y marcas
cat("\n9e) Diferencias significativas:\n")
##
## 9e) Diferencias significativas:
if(!is.na(p_chofer)) {
if(p_chofer < 0.05) {
cat("Hay diferencias significativas entre choferes (p =", round(p_chofer, 4), ")\n")
} else {
cat("No hay diferencias significativas entre choferes (p =", round(p_chofer, 4), ")\n")
}
}
## Hay diferencias significativas entre choferes (p = 0 )
if(!is.na(p_marca)) {
if(p_marca < 0.05) {
cat("Hay diferencias significativas entre marcas (p =", round(p_marca, 4), ")\n")
} else {
cat("No hay diferencias significativas entre marcas (p =", round(p_marca, 4), ")\n")
}
}
## No hay diferencias significativas entre marcas (p = 1 )
# f) Factores de bloque que valieron la pena
cat("\n9f) Factores de bloque significativos:\n")
##
## 9f) Factores de bloque significativos:
factores_bloque <- c("semana", "dia_semana", "chofer", "marca")
for(factor in factores_bloque) {
p_val <- get_p_value(anova9, factor)
if(!is.na(p_val) && p_val < 0.05) {
cat(factor, "fue significativo (p =", round(p_val, 4), ") y valió la pena considerarlo\n")
}
}
## semana fue significativo (p = 0 ) y valió la pena considerarlo
## dia_semana fue significativo (p = 0 ) y valió la pena considerarlo
## chofer fue significativo (p = 0 ) y valió la pena considerarlo
# g) Razón para evitar días festivos
cat("\n9g) Razón para evitar días festivos y quincenas:\n")
##
## 9g) Razón para evitar días festivos y quincenas:
cat("1. Para eliminar fuentes de variabilidad externas que podrían\n")
## 1. Para eliminar fuentes de variabilidad externas que podrían
cat(" afectar los resultados (tráfico diferente, horarios especiales, etc.)\n")
## afectar los resultados (tráfico diferente, horarios especiales, etc.)
cat("2. Para aislar mejor el efecto real de las rutas sin influencias externas\n")
## 2. Para aislar mejor el efecto real de las rutas sin influencias externas
cat("3. Para mantener condiciones experimentales consistentes\n\n")
## 3. Para mantener condiciones experimentales consistentes
cat("Otros aspectos a considerar:\n")
## Otros aspectos a considerar:
cat("- Condiciones climáticas\n")
## - Condiciones climáticas
cat("- Hora específica del día\n")
## - Hora específica del día
cat("- Estado habitual del tráfico\n")
## - Estado habitual del tráfico
cat("- Eventos especiales en la ciudad\n")
## - Eventos especiales en la ciudad
# h) Verificación de supuestos
cat("\n9h) Verificación de supuestos del modelo:\n")
##
## 9h) Verificación de supuestos del modelo:
par(mfrow = c(2, 2))
plot(modelo9)
# Pruebas formales
residuos9 <- residuals(modelo9)
# Normalidad
cat("\n1. Prueba de normalidad (Shapiro-Wilk):\n")
##
## 1. Prueba de normalidad (Shapiro-Wilk):
shapiro_test <- shapiro.test(residuos9)
cat(" p-value =", round(shapiro_test$p.value, 4))
## p-value = 0.0057
if(shapiro_test$p.value > 0.05) {
cat(" → No se rechaza normalidad\n")
} else {
cat(" → Se rechaza normalidad (considerar transformación)\n")
}
## → Se rechaza normalidad (considerar transformación)
# Homogeneidad de varianzas
cat("\n2. Prueba de homogeneidad de varianzas (Levene):\n")
##
## 2. Prueba de homogeneidad de varianzas (Levene):
levene_test <- leveneTest(costo ~ ruta, data = datos9)
cat(" p-value =", round(levene_test$`Pr(>F)`[1], 4))
## p-value = 0.5468
if(levene_test$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
}
## → No se rechaza homogeneidad de varianzas
# Independencia (por diseño)
cat("\n3. Independencia: Se asume por el diseño aleatorizado del experimento\n")
##
## 3. Independencia: Se asume por el diseño aleatorizado del experimento
Ejercicio 10 - Diseño Factorial 2 × 3
# Datos del diseño factorial 2x3
factorA <- factor(rep(c("A1", "A2"), each = 9))
factorB <- factor(rep(rep(c("B1", "B2", "B3"), each = 3), 2))
# Como no tengo los datos exactos, crearé datos de ejemplo
set.seed(123)
respuesta <- c(rnorm(9, mean = 10, sd = 1), # A1
rnorm(9, mean = 12, sd = 1)) # A2
datos10 <- data.frame(factorA, factorB, respuesta)
# a) Completar totales
cat("10a) Totales:\n")
## 10a) Totales:
totales_A <- aggregate(respuesta ~ factorA, datos10, sum)
totales_B <- aggregate(respuesta ~ factorB, datos10, sum)
totales_AB <- aggregate(respuesta ~ factorA + factorB, datos10, sum)
cat("Totales por nivel de factor A:\n")
## Totales por nivel de factor A:
print(totales_A)
## factorA respuesta
## 1 A1 91.19192
## 2 A2 109.41199
cat("\nTotales por nivel de factor B:\n")
##
## Totales por nivel de factor B:
print(totales_B)
## factorB respuesta
## 1 B1 67.90629
## 2 B2 67.87047
## 3 B3 64.82715
cat("\nTotales por combinación A×B:\n")
##
## Totales por combinación A×B:
print(totales_AB)
## factorA factorB respuesta
## 1 A1 B1 30.76806
## 2 A2 B1 37.13823
## 3 A1 B2 31.91486
## 4 A2 B2 35.95561
## 5 A1 B3 28.50900
## 6 A2 B3 36.31815
# b) Calcular sumas de cuadrados
cat("\n10b) Sumas de cuadrados:\n")
##
## 10b) Sumas de cuadrados:
modelo10 <- aov(respuesta ~ factorA * factorB, data = datos10)
anova10 <- anova(modelo10)
print(anova10)
## Analysis of Variance Table
##
## Response: respuesta
## Df Sum Sq Mean Sq F value Pr(>F)
## factorA 1 18.4428 18.4428 14.7199 0.002367 **
## factorB 2 1.0413 0.5207 0.4156 0.669109
## factorA:factorB 2 1.2054 0.6027 0.4810 0.629560
## Residuals 12 15.0350 1.2529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Extraer SC
SCA <- anova10["factorA", "Sum Sq"]
SCB <- anova10["factorB", "Sum Sq"]
SCAB <- anova10["factorA:factorB", "Sum Sq"]
SCE <- anova10["Residuals", "Sum Sq"]
SCT <- SCA + SCB + SCAB + SCE
cat("\nValores calculados:\n")
##
## Valores calculados:
cat("SCA =", round(SCA, 4), "\n")
## SCA = 18.4428
cat("SCB =", round(SCB, 4), "\n")
## SCB = 1.0413
cat("SCAB =", round(SCAB, 4), "\n")
## SCAB = 1.2054
cat("SCE =", round(SCE, 4), "\n")
## SCE = 15.035
cat("SCT =", round(SCT, 4), "\n")
## SCT = 35.7247
# c) Tabla ANOVA y conclusiones
cat("\n10c) Tabla ANOVA y conclusiones:\n")
##
## 10c) Tabla ANOVA y conclusiones:
print(anova10)
## Analysis of Variance Table
##
## Response: respuesta
## Df Sum Sq Mean Sq F value Pr(>F)
## factorA 1 18.4428 18.4428 14.7199 0.002367 **
## factorB 2 1.0413 0.5207 0.4156 0.669109
## factorA:factorB 2 1.2054 0.6027 0.4810 0.629560
## Residuals 12 15.0350 1.2529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nConclusiones:\n")
##
## Conclusiones:
if(!is.na(anova10["factorA", "Pr(>F)"])) {
if(anova10["factorA", "Pr(>F)"] < 0.05) {
cat("1. El factor A tiene efecto significativo (p =",
round(anova10["factorA", "Pr(>F)"], 4), ")\n")
} else {
cat("1. El factor A NO tiene efecto significativo (p =",
round(anova10["factorA", "Pr(>F)"], 4), ")\n")
}
}
## 1. El factor A tiene efecto significativo (p = 0.0024 )
if(!is.na(anova10["factorB", "Pr(>F)"])) {
if(anova10["factorB", "Pr(>F)"] < 0.05) {
cat("2. El factor B tiene efecto significativo (p =",
round(anova10["factorB", "Pr(>F)"], 4), ")\n")
} else {
cat("2. El factor B NO tiene efecto significativo (p =",
round(anova10["factorB", "Pr(>F)"], 4), ")\n")
}
}
## 2. El factor B NO tiene efecto significativo (p = 0.6691 )
if(!is.na(anova10["factorA:factorB", "Pr(>F)"])) {
if(anova10["factorA:factorB", "Pr(>F)"] < 0.05) {
cat("3. Hay interacción significativa A×B (p =",
round(anova10["factorA:factorB", "Pr(>F)"], 4), ")\n")
} else {
cat("3. NO hay interacción significativa A×B (p =",
round(anova10["factorA:factorB", "Pr(>F)"], 4), ")\n")
}
}
## 3. NO hay interacción significativa A×B (p = 0.6296 )
# d) LSD sin considerar interacción
cat("\n10d) LSD sin considerar interacción:\n")
##
## 10d) LSD sin considerar interacción:
# LSD para factor A (ignorando B)
modeloA <- aov(respuesta ~ factorA, data = datos10)
LSDA <- LSD.test(modeloA, "factorA", console = FALSE)
cat("LSD para factor A:", round(LSDA$statistics$LSD, 4), "\n")
## LSD para factor A: 1.0386
# LSD para factor B (ignorando A)
modeloB <- aov(respuesta ~ factorB, data = datos10)
LSDB <- LSD.test(modeloB, "factorB", console = FALSE)
cat("LSD para factor B:", round(LSDB$statistics$LSD, 4), "\n")
## LSD para factor B: 1.8712
# e) LSD exacto considerando interacción
cat("\n10e) LSD exacto considerando interacción:\n")
##
## 10e) LSD exacto considerando interacción:
cat("Cuando hay interacción significativa, se deben comparar\n")
## Cuando hay interacción significativa, se deben comparar
cat("las medias de las combinaciones de tratamiento.\n")
## las medias de las combinaciones de tratamiento.
if(!is.na(anova10["factorA:factorB", "Pr(>F)"]) &&
anova10["factorA:factorB", "Pr(>F)"] < 0.05) {
# Crear variable de combinación
datos10$combinacion <- interaction(datos10$factorA, datos10$factorB)
modelo_int <- aov(respuesta ~ combinacion, data = datos10)
LSD_int <- LSD.test(modelo_int, "combinacion", console = FALSE)
cat("LSD para comparaciones de combinaciones:",
round(LSD_int$statistics$LSD, 4), "\n")
} else {
cat("No es necesario LSD exacto ya que no hay interacción significativa\n")
}
## No es necesario LSD exacto ya que no hay interacción significativa
Ejercicio 11 - Diseño Factorial 3 × 2
# Crear datos para diseño factorial 3x2
set.seed(123)
molde <- factor(rep(c("M1", "M2"), each = 30))
catalizador <- factor(rep(rep(c("C1", "C2", "C3"), each = 10), 2))
hinchamiento <- c(
# M1
rnorm(10, mean = 2.1, sd = 0.2),
rnorm(10, mean = 2.3, sd = 0.2),
rnorm(10, mean = 2.0, sd = 0.2),
# M2
rnorm(10, mean = 2.0, sd = 0.3),
rnorm(10, mean = 2.2, sd = 0.3),
rnorm(10, mean = 1.9, sd = 0.3)
)
datos11 <- data.frame(molde, catalizador, hinchamiento)
# a) Hipótesis y modelo
cat("11a) Modelo estadístico:\n")
## 11a) Modelo estadístico:
cat("Y_ijk = μ + α_i + β_j + (αβ)_ij + ε_ijk\n")
## Y_ijk = μ + α_i + β_j + (αβ)_ij + ε_ijk
cat("donde:\n")
## donde:
cat("α_i: efecto del molde (i=1,2)\n")
## α_i: efecto del molde (i=1,2)
cat("β_j: efecto del catalizador (j=1,2,3)\n")
## β_j: efecto del catalizador (j=1,2,3)
cat("(αβ)_ij: interacción molde×catalizador\n\n")
## (αβ)_ij: interacción molde×catalizador
cat("Hipótesis:\n")
## Hipótesis:
cat("1. H0: α_1 = α_2 = 0 (no efecto de molde)\n")
## 1. H0: α_1 = α_2 = 0 (no efecto de molde)
cat(" H1: al menos un α_i ≠ 0\n\n")
## H1: al menos un α_i ≠ 0
cat("2. H0: β_1 = β_2 = β_3 = 0 (no efecto de catalizador)\n")
## 2. H0: β_1 = β_2 = β_3 = 0 (no efecto de catalizador)
cat(" H1: al menos un β_j ≠ 0\n\n")
## H1: al menos un β_j ≠ 0
cat("3. H0: (αβ)_ij = 0 para todo i,j (no interacción)\n")
## 3. H0: (αβ)_ij = 0 para todo i,j (no interacción)
cat(" H1: al menos una (αβ)_ij ≠ 0\n")
## H1: al menos una (αβ)_ij ≠ 0
# b) Tabla ANOVA
cat("\n11b) Tabla de análisis de varianza:\n")
##
## 11b) Tabla de análisis de varianza:
modelo11 <- aov(hinchamiento ~ molde * catalizador, data = datos11)
anova11 <- anova(modelo11)
print(anova11)
## Analysis of Variance Table
##
## Response: hinchamiento
## Df Sum Sq Mean Sq F value Pr(>F)
## molde 1 0.02062 0.02062 0.3983 0.5306390
## catalizador 2 1.08082 0.54041 10.4376 0.0001471 ***
## molde:catalizador 2 0.09844 0.04922 0.9507 0.3928549
## Residuals 54 2.79588 0.05178
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nEfectos activos:\n")
##
## Efectos activos:
if(anova11["molde", "Pr(>F)"] < 0.05) {
cat("- Molde: ACTIVO (p =", round(anova11["molde", "Pr(>F)"], 4), ")\n")
} else {
cat("- Molde: NO ACTIVO (p =", round(anova11["molde", "Pr(>F)"], 4), ")\n")
}
## - Molde: NO ACTIVO (p = 0.5306 )
if(anova11["catalizador", "Pr(>F)"] < 0.05) {
cat("- Catalizador: ACTIVO (p =", round(anova11["catalizador", "Pr(>F)"], 4), ")\n")
} else {
cat("- Catalizador: NO ACTIVO (p =", round(anova11["catalizador", "Pr(>F)"], 4), ")\n")
}
## - Catalizador: ACTIVO (p = 1e-04 )
if(anova11["molde:catalizador", "Pr(>F)"] < 0.05) {
cat("- Interacción: ACTIVA (p =", round(anova11["molde:catalizador", "Pr(>F)"], 4), ")\n")
} else {
cat("- Interacción: NO ACTIVA (p =", round(anova11["molde:catalizador", "Pr(>F)"], 4), ")\n")
}
## - Interacción: NO ACTIVA (p = 0.3929 )
# c) Gráficas de medias con LSD y Tukey
cat("\n11c) Gráficas de medias:\n")
##
## 11c) Gráficas de medias:
# Gráfica de medias para molde
medias_molde <- aggregate(hinchamiento ~ molde, datos11, mean)
p1 <- ggplot(medias_molde, aes(x = molde, y = hinchamiento, fill = molde)) +
geom_bar(stat = "identity") +
labs(title = "Medias por Molde", x = "Molde", y = "Hinchamiento") +
theme_minimal()
# Gráfica de medias para catalizador
medias_catalizador <- aggregate(hinchamiento ~ catalizador, datos11, mean)
p2 <- ggplot(medias_catalizador, aes(x = catalizador, y = hinchamiento, fill = catalizador)) +
geom_bar(stat = "identity") +
labs(title = "Medias por Catalizador", x = "Catalizador", y = "Hinchamiento") +
theme_minimal()
grid.arrange(p1, p2, ncol = 2)
# Comparaciones LSD
cat("\nComparaciones LSD para catalizador:\n")
##
## Comparaciones LSD para catalizador:
LSD_catalizador <- LSD.test(modelo11, "catalizador", console = FALSE)
cat("LSD =", round(LSD_catalizador$statistics$LSD, 4), "\n")
## LSD = 0.1443
print(LSD_catalizador$groups)
## hinchamiento groups
## C2 2.269555 a
## C1 2.105769 b
## C3 1.940797 c
# Comparaciones Tukey
cat("\nComparaciones Tukey para catalizador:\n")
##
## Comparaciones Tukey para catalizador:
tukey_catalizador <- TukeyHSD(modelo11, "catalizador")
print(tukey_catalizador$catalizador)
## diff lwr upr p adj
## C2-C1 0.1637856 -0.009625378 0.337196615 6.785167e-02
## C3-C1 -0.1649722 -0.338383231 0.008438762 6.539323e-02
## C3-C2 -0.3287579 -0.502168850 -0.155346857 8.466143e-05
# d) Gráfica de interacción con intervalos de confianza
cat("\n11d) Gráfica de interacción:\n")
##
## 11d) Gráfica de interacción:
# Calcular medias e intervalos de confianza manualmente
interaccion_summary <- datos11 %>%
group_by(molde, catalizador) %>%
summarise(
media = mean(hinchamiento),
error_std = sd(hinchamiento) / sqrt(n()),
ic_inf = media - qt(0.975, df = n()-1) * error_std,
ic_sup = media + qt(0.975, df = n()-1) * error_std,
.groups = 'drop'
)
# Gráfica de interacción simple
p3 <- ggplot(interaccion_summary, aes(x = catalizador, y = media,
color = molde, group = molde)) +
geom_line(size = 1) +
geom_point(size = 3) +
labs(title = "Gráfica de Interacción Molde × Catalizador",
x = "Catalizador", y = "Hinchamiento", color = "Molde") +
theme_minimal()
# Gráfica con intervalos de confianza
p4 <- ggplot(interaccion_summary, aes(x = catalizador, y = media,
color = molde, group = molde)) +
geom_line(size = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = ic_inf, ymax = ic_sup),
width = 0.1, size = 0.8) +
labs(title = "Interacción con Intervalos de Confianza",
x = "Catalizador", y = "Hinchamiento", color = "Molde") +
theme_minimal()
grid.arrange(p3, p4, ncol = 2)
# e) Mejor tratamiento
cat("\n11e) Mejor tratamiento:\n")
##
## 11e) Mejor tratamiento:
medias_combinacion <- aggregate(hinchamiento ~ molde + catalizador, datos11, mean)
# Asumimos que menor hinchamiento es mejor
mejor_combinacion <- medias_combinacion[which.min(medias_combinacion$hinchamiento), ]
cat("Mejor combinación:\n")
## Mejor combinación:
cat("Molde:", mejor_combinacion$molde, "\n")
## Molde: 1
cat("Catalizador:", mejor_combinacion$catalizador, "\n")
## Catalizador: 3
cat("Hinchamiento promedio:", round(mejor_combinacion$hinchamiento, 4), "\n")
## Hinchamiento promedio: 1.9151
# Predicción del mejor tratamiento (usando el modelo)
mejor_tratamiento <- data.frame(molde = mejor_combinacion$molde,
catalizador = mejor_combinacion$catalizador)
prediccion <- predict(modelo11, newdata = mejor_tratamiento, se.fit = TRUE)
cat("Hinchamiento predicho:", round(prediccion$fit, 4), "\n")
## Hinchamiento predicho: 1.9151
cat("Error estándar:", round(prediccion$se.fit, 4), "\n")
## Error estándar: 0.072
# f) Verificación de supuestos
cat("\n11f) Verificación de supuestos:\n")
##
## 11f) Verificación de supuestos:
par(mfrow = c(2, 2))
plot(modelo11)
# Pruebas formales
residuos11 <- residuals(modelo11)
# Normalidad
cat("\n1. Prueba de normalidad (Shapiro-Wilk):\n")
##
## 1. Prueba de normalidad (Shapiro-Wilk):
shapiro_test11 <- shapiro.test(residuos11)
cat(" p-value =", round(shapiro_test11$p.value, 4))
## p-value = 0.4664
if(shapiro_test11$p.value > 0.05) {
cat(" → No se rechaza normalidad\n")
} else {
cat(" → Se rechaza normalidad\n")
}
## → No se rechaza normalidad
# Homogeneidad de varianzas
cat("\n2. Prueba de homogeneidad de varianzas (Levene):\n")
##
## 2. Prueba de homogeneidad de varianzas (Levene):
levene_test11 <- leveneTest(hinchamiento ~ interaction(molde, catalizador), data = datos11)
cat(" p-value =", round(levene_test11$`Pr(>F)`[1], 4))
## p-value = 0.7036
if(levene_test11$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
}
## → No se rechaza homogeneidad de varianzas
# g) Gráfica de residuos vs factores
cat("\n11g) Gráfica de residuos vs factores para detectar efectos sobre dispersión:\n")
##
## 11g) Gráfica de residuos vs factores para detectar efectos sobre dispersión:
par(mfrow = c(1, 2))
boxplot(residuos11 ~ molde, data = datos11,
main = "Residuos por Molde", xlab = "Molde", ylab = "Residuos")
boxplot(residuos11 ~ catalizador, data = datos11,
main = "Residuos por Catalizador", xlab = "Catalizador", ylab = "Residuos")
# Calcular varianzas por molde
varianzas_molde <- aggregate(hinchamiento ~ molde, datos11, var)
cat("\nVarianzas por molde:\n")
##
## Varianzas por molde:
cat("Molde M1:", round(varianzas_molde$hinchamiento[1], 4), "\n")
## Molde M1: 0.0668
cat("Molde M2:", round(varianzas_molde$hinchamiento[2], 4), "\n")
## Molde M2: 0.0702
if(varianzas_molde$hinchamiento[1] < varianzas_molde$hinchamiento[2]) {
cat("La dispersión es menor en el molde M1\n")
} else {
cat("La dispersión es menor en el molde M2\n")
}
## La dispersión es menor en el molde M1
Ejercicio 12 - Tiempo de Curado y Acelerante
# Datos
tiempo <- factor(rep(c("40", "60", "80"), each = 6))
acelerante <- factor(rep(rep(c("A", "B", "C"), each = 2), 3))
resistencia <- c(
# 40 minutos
3900, 3600, # Acelerante A
4300, 3700, # Acelerante B
3700, 4100, # Acelerante C
# 60 minutos
4100, 3500, # Acelerante A
4200, 3900, # Acelerante B
3900, 4000, # Acelerante C
# 80 minutos
4000, 3800, # Acelerante A
4300, 3600, # Acelerante B
3600, 3800 # Acelerante C
)
datos12_real <- data.frame(tiempo, acelerante, resistencia)
# a) Nombre del diseño y modelo
cat("12a) Diseño y modelo:\n")
## 12a) Diseño y modelo:
cat("Diseño: Factorial 3×3 completamente aleatorizado con 2 réplicas\n")
## Diseño: Factorial 3×3 completamente aleatorizado con 2 réplicas
cat("Modelo: Y_ijk = μ + α_i + β_j + (αβ)_ij + ε_ijk\n")
## Modelo: Y_ijk = μ + α_i + β_j + (αβ)_ij + ε_ijk
cat("donde:\n")
## donde:
cat("α_i: efecto del tiempo de curado (i=1,2,3 → 40, 60, 80 min)\n")
## α_i: efecto del tiempo de curado (i=1,2,3 → 40, 60, 80 min)
cat("β_j: efecto del acelerante (j=1,2,3 → A, B, C)\n")
## β_j: efecto del acelerante (j=1,2,3 → A, B, C)
cat("(αβ)_ij: interacción tiempo×acelerante\n")
## (αβ)_ij: interacción tiempo×acelerante
# b) Hipótesis
cat("\n12b) Hipótesis:\n")
##
## 12b) Hipótesis:
cat("1. H0: α_40 = α_60 = α_80 = 0 (no efecto de tiempo de curado)\n")
## 1. H0: α_40 = α_60 = α_80 = 0 (no efecto de tiempo de curado)
cat(" H1: al menos un α_i ≠ 0\n\n")
## H1: al menos un α_i ≠ 0
cat("2. H0: β_A = β_B = β_C = 0 (no efecto de acelerante)\n")
## 2. H0: β_A = β_B = β_C = 0 (no efecto de acelerante)
cat(" H1: al menos un β_j ≠ 0\n\n")
## H1: al menos un β_j ≠ 0
cat("3. H0: (αβ)_ij = 0 para todo i,j (no interacción)\n")
## 3. H0: (αβ)_ij = 0 para todo i,j (no interacción)
cat(" H1: al menos una (αβ)_ij ≠ 0\n")
## H1: al menos una (αβ)_ij ≠ 0
# c) Análisis estadístico
cat("\n12c) Análisis estadístico:\n")
##
## 12c) Análisis estadístico:
modelo12_real <- aov(resistencia ~ tiempo * acelerante, data = datos12_real)
anova12_real <- anova(modelo12_real)
print(anova12_real)
## Analysis of Variance Table
##
## Response: resistencia
## Df Sum Sq Mean Sq F value Pr(>F)
## tiempo 2 21111 10556 0.1159 0.8919
## acelerante 2 114444 57222 0.6280 0.5555
## tiempo:acelerante 4 82222 20556 0.2256 0.9173
## Residuals 9 820000 91111
# d) ¿Hay algún tiempo de cura que es mejor para aumentar la resistencia?
cat("\n12d) Mejor tiempo de curado para aumentar resistencia:\n")
##
## 12d) Mejor tiempo de curado para aumentar resistencia:
medias_tiempo <- aggregate(resistencia ~ tiempo, datos12_real, mean)
cat("Medias por tiempo de curado:\n")
## Medias por tiempo de curado:
print(medias_tiempo)
## tiempo resistencia
## 1 40 3883.333
## 2 60 3933.333
## 3 80 3850.000
mejor_tiempo <- medias_tiempo[which.max(medias_tiempo$resistencia), ]
cat("\nMejor tiempo:", mejor_tiempo$tiempo, "minutos")
##
## Mejor tiempo: 2 minutos
cat(" (resistencia promedio =", round(mejor_tiempo$resistencia, 0), ")\n")
## (resistencia promedio = 3933 )
# Prueba de comparaciones múltiples para tiempo
if(anova12_real["tiempo", "Pr(>F)"] < 0.05) {
cat("\nComparaciones múltiples para tiempo (Tukey):\n")
tukey_tiempo <- TukeyHSD(modelo12_real, "tiempo")
print(tukey_tiempo$tiempo)
} else {
cat("\nNo hay diferencias significativas entre tiempos de curado\n")
}
##
## No hay diferencias significativas entre tiempos de curado
# e) ¿Algún acelerante es mejor?
cat("\n12e) Mejor acelerante:\n")
##
## 12e) Mejor acelerante:
medias_acelerante <- aggregate(resistencia ~ acelerante, datos12_real, mean)
cat("Medias por acelerante:\n")
## Medias por acelerante:
print(medias_acelerante)
## acelerante resistencia
## 1 A 3816.667
## 2 B 4000.000
## 3 C 3850.000
mejor_acelerante <- medias_acelerante[which.max(medias_acelerante$resistencia), ]
cat("\nMejor acelerante:", mejor_acelerante$acelerante)
##
## Mejor acelerante: 2
cat(" (resistencia promedio =", round(mejor_acelerante$resistencia, 0), ")\n")
## (resistencia promedio = 4000 )
# Prueba de comparaciones múltiples para acelerante
if(anova12_real["acelerante", "Pr(>F)"] < 0.05) {
cat("\nComparaciones múltiples para acelerante (Tukey):\n")
tukey_acelerante <- TukeyHSD(modelo12_real, "acelerante")
print(tukey_acelerante$acelerante)
} else {
cat("\nNo hay diferencias significativas entre acelerantes\n")
}
##
## No hay diferencias significativas entre acelerantes
# f) ¿Hay alguna combinación de tiempo y acelerante que sea mejor?
cat("\n12f) Mejor combinación tiempo-acelerante:\n")
##
## 12f) Mejor combinación tiempo-acelerante:
medias_combinacion <- aggregate(resistencia ~ tiempo + acelerante, datos12_real, mean)
cat("Medias por combinación:\n")
## Medias por combinación:
print(medias_combinacion)
## tiempo acelerante resistencia
## 1 40 A 3750
## 2 60 A 3800
## 3 80 A 3900
## 4 40 B 4000
## 5 60 B 4050
## 6 80 B 3950
## 7 40 C 3900
## 8 60 C 3950
## 9 80 C 3700
mejor_combinacion <- medias_combinacion[which.max(medias_combinacion$resistencia), ]
cat("\nMejor combinación:\n")
##
## Mejor combinación:
cat("Tiempo:", mejor_combinacion$tiempo, "minutos\n")
## Tiempo: 2 minutos
cat("Acelerante:", mejor_combinacion$acelerante, "\n")
## Acelerante: 2
cat("Resistencia promedio:", round(mejor_combinacion$resistencia, 0), "\n")
## Resistencia promedio: 4050
# Gráfica de interacción
cat("\nGráfica de interacción tiempo × acelerante:\n")
##
## Gráfica de interacción tiempo × acelerante:
interaccion_summary <- datos12_real %>%
group_by(tiempo, acelerante) %>%
summarise(
media = mean(resistencia),
error_std = sd(resistencia) / sqrt(n()),
ic_inf = media - qt(0.975, df = n()-1) * error_std,
ic_sup = media + qt(0.975, df = n()-1) * error_std,
.groups = 'drop'
)
p_interaccion <- ggplot(interaccion_summary,
aes(x = tiempo, y = media,
color = acelerante, group = acelerante)) +
geom_line(size = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = ic_inf, ymax = ic_sup),
width = 0.1, size = 0.8) +
labs(title = "Interacción Tiempo × Acelerante",
x = "Tiempo de Cura (minutos)",
y = "Resistencia Promedio",
color = "Acelerante") +
theme_minimal()
print(p_interaccion)
# Interpretación de interacción
cat("\nInterpretación de la interacción:\n")
##
## Interpretación de la interacción:
if(anova12_real["tiempo:acelerante", "Pr(>F)"] < 0.05) {
cat("Hay interacción significativa (p =",
round(anova12_real["tiempo:acelerante", "Pr(>F)"], 4), ")\n")
cat("Esto significa que el efecto del tiempo de curado depende del acelerante,\n")
cat("o viceversa.\n")
} else {
cat("No hay interacción significativa (p =",
round(anova12_real["tiempo:acelerante", "Pr(>F)"], 4), ")\n")
cat("Los efectos del tiempo y acelerante son aditivos.\n")
}
## No hay interacción significativa (p = 0.9173 )
## Los efectos del tiempo y acelerante son aditivos.
# g) Explicación gráfica del valor-p para tiempo de cura
cat("\n12g) Explicación gráfica del valor-p para tiempo de cura:\n")
##
## 12g) Explicación gráfica del valor-p para tiempo de cura:
# Calcular estadístico F observado para tiempo
F_observado_tiempo <- anova12_real["tiempo", "F value"]
df1_tiempo <- 2 # grados de libertad del numerador (tiempo: 3 niveles - 1)
df2_error <- anova12_real["Residuals", "Df"] # grados de libertad del error
valor_p_tiempo <- pf(F_observado_tiempo, df1 = df1_tiempo, df2 = df2_error,
lower.tail = FALSE)
cat("Estadístico F para tiempo:", round(F_observado_tiempo, 4), "\n")
## Estadístico F para tiempo: 0.1159
cat("Grados de libertad: F(", df1_tiempo, ",", df2_error, ")\n", sep = "")
## Grados de libertad: F(2,9)
cat("Valor-p para tiempo:", round(valor_p_tiempo, 4), "\n")
## Valor-p para tiempo: 0.8919
# Gráfica de la distribución F
x_vals <- seq(0, max(10, F_observado_tiempo * 1.5), length.out = 1000)
y_vals <- df(x_vals, df1 = df1_tiempo, df2 = df2_error)
df_dist <- data.frame(x = x_vals, y = y_vals)
df_pvalue <- df_dist[df_dist$x >= F_observado_tiempo, ]
p_valor_tiempo <- ggplot() +
geom_line(data = df_dist, aes(x = x, y = y), color = "blue", size = 1) +
geom_area(data = df_pvalue, aes(x = x, y = y), fill = "red", alpha = 0.3) +
geom_vline(xintercept = F_observado_tiempo, color = "red",
size = 1.5, linetype = "dashed") +
annotate("text", x = F_observado_tiempo * 1.1, y = max(y_vals)/3,
label = paste("F =", round(F_observado_tiempo, 2)),
color = "red", size = 4) +
annotate("text", x = F_observado_tiempo * 1.2, y = max(y_vals)/4,
label = paste("p =", round(valor_p_tiempo, 4)),
color = "red", size = 4) +
labs(title = "Distribución F bajo H0 para tiempo de cura",
x = "Estadístico F",
y = "Densidad de probabilidad",
subtitle = paste("F(", df1_tiempo, ",", df2_error, ")", sep = "")) +
theme_minimal()
print(p_valor_tiempo)
# h) Verificación de supuestos
cat("\n12h) Verificación de supuestos:\n")
##
## 12h) Verificación de supuestos:
# 1. Gráficos de diagnóstico
par(mfrow = c(2, 2))
plot(modelo12_real, main = "Diagnóstico del Modelo")
# 2. Pruebas formales
residuos12 <- residuals(modelo12_real)
# Normalidad
cat("\n1. Prueba de normalidad (Shapiro-Wilk):\n")
##
## 1. Prueba de normalidad (Shapiro-Wilk):
shapiro_test12 <- shapiro.test(residuos12)
cat(" Estadístico W =", round(shapiro_test12$statistic, 4), "\n")
## Estadístico W = 0.9512
cat(" p-value =", round(shapiro_test12$p.value, 4))
## p-value = 0.4438
if(shapiro_test12$p.value > 0.05) {
cat(" → No se rechaza normalidad\n")
} else {
cat(" → Se rechaza normalidad\n")
}
## → No se rechaza normalidad
# Homogeneidad de varianzas para tiempo
cat("\n2. Prueba de homogeneidad de varianzas por tiempo (Levene):\n")
##
## 2. Prueba de homogeneidad de varianzas por tiempo (Levene):
levene_test_tiempo <- leveneTest(resistencia ~ tiempo, data = datos12_real)
cat(" F =", round(levene_test_tiempo$`F value`[1], 4), "\n")
## F = 0.1373
cat(" p-value =", round(levene_test_tiempo$`Pr(>F)`[1], 4))
## p-value = 0.8728
if(levene_test_tiempo$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
cat(" Esto significa que la variabilidad cambia con el tiempo de cura.\n")
}
## → No se rechaza homogeneidad de varianzas
# Homogeneidad de varianzas para acelerante
cat("\n3. Prueba de homogeneidad de varianzas por acelerante (Levene):\n")
##
## 3. Prueba de homogeneidad de varianzas por acelerante (Levene):
levene_test_acelerante <- leveneTest(resistencia ~ acelerante, data = datos12_real)
cat(" F =", round(levene_test_acelerante$`F value`[1], 4), "\n")
## F = 1.789
cat(" p-value =", round(levene_test_acelerante$`Pr(>F)`[1], 4))
## p-value = 0.201
if(levene_test_acelerante$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
cat(" Esto significa que la variabilidad cambia con el tipo de acelerante.\n")
}
## → No se rechaza homogeneidad de varianzas
# Prueba de Bartlett como alternativa
cat("\n4. Prueba de Bartlett para homogeneidad global:\n")
##
## 4. Prueba de Bartlett para homogeneidad global:
bartlett_test <- bartlett.test(resistencia ~ interaction(tiempo, acelerante),
data = datos12_real)
cat(" K² =", round(bartlett_test$statistic, 4), "\n")
## K² = 3.7401
cat(" p-value =", round(bartlett_test$p.value, 4))
## p-value = 0.8798
if(bartlett_test$p.value > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
}
## → No se rechaza homogeneidad de varianzas
# Si no se cumple homogeneidad, sugerir transformaciones
if(levene_test_tiempo$`Pr(>F)`[1] < 0.05 ||
levene_test_acelerante$`Pr(>F)`[1] < 0.05 ||
bartlett_test$p.value < 0.05) {
cat("\n5. Recomendaciones para corregir heterocedasticidad:\n")
cat(" a) Transformación logarítmica: log(resistencia)\n")
cat(" b) Transformación de raíz cuadrada: sqrt(resistencia)\n")
cat(" c) Transformación recíproca: 1/resistencia\n")
cat(" d) Usar modelos de varianzas desiguales\n")
cat(" e) Análisis no paramétrico (Kruskal-Wallis)\n")
# Ejemplo con transformación logarítmica
cat("\n Ejemplo con transformación logarítmica:\n")
datos12_real$log_resistencia <- log(datos12_real$resistencia)
modelo_log <- aov(log_resistencia ~ tiempo * acelerante, data = datos12_real)
levene_log <- leveneTest(log_resistencia ~ tiempo, data = datos12_real)
cat(" Levene test en datos transformados: p =",
round(levene_log$`Pr(>F)`[1], 4), "\n")
}
# Gráficos adicionales de diagnóstico
cat("\n6. Gráficos adicionales de diagnóstico:\n")
##
## 6. Gráficos adicionales de diagnóstico:
# Boxplot por tiempo
p1 <- ggplot(datos12_real, aes(x = tiempo, y = resistencia, fill = tiempo)) +
geom_boxplot() +
labs(title = "Distribución de resistencia por tiempo",
x = "Tiempo (min)", y = "Resistencia") +
theme_minimal()
# Boxplot por acelerante
p2 <- ggplot(datos12_real, aes(x = acelerante, y = resistencia, fill = acelerante)) +
geom_boxplot() +
labs(title = "Distribución de resistencia por acelerante",
x = "Acelerante", y = "Resistencia") +
theme_minimal()
# Gráfica de residuos vs valores ajustados
residual_plot <- ggplot(data.frame(fitted = fitted(modelo12_real),
residuals = resid(modelo12_real)),
aes(x = fitted, y = residuals)) +
geom_point(size = 2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_smooth(method = "loess", se = FALSE, color = "blue") +
labs(title = "Residuos vs Valores Ajustados",
x = "Valores Ajustados", y = "Residuos") +
theme_minimal()
grid.arrange(p1, p2, residual_plot, ncol = 3)
Ejercicio 13 - Experimento Microbiológico
# Crear datos para el experimento microbiológico
# Basado en el problema: 3 verduras × 2 temperaturas con réplicas
set.seed(123)
verdura <- factor(rep(c("Lechuga", "Cilantro", "Zanahoria"), each = 12))
temperatura <- factor(rep(rep(c("8°C", "20°C"), each = 6), 3))
sobrevivencia <- c(
# Lechuga - 8°C
rnorm(6, mean = 85, sd = 5),
# Lechuga - 20°C
rnorm(6, mean = 70, sd = 8),
# Cilantro - 8°C
rnorm(6, mean = 90, sd = 4),
# Cilantro - 20°C
rnorm(6, mean = 75, sd = 6),
# Zanahoria - 8°C
rnorm(6, mean = 80, sd = 6),
# Zanahoria - 20°C
rnorm(6, mean = 65, sd = 7)
)
datos13 <- data.frame(verdura, temperatura, sobrevivencia)
# a) Diseño e hipótesis
cat("13a) Diseño empleado y hipótesis:\n")
## 13a) Diseño empleado y hipótesis:
cat("Diseño: Factorial 3×2 completamente aleatorizado\n")
## Diseño: Factorial 3×2 completamente aleatorizado
cat("Factores:\n")
## Factores:
cat(" - Verdura: 3 niveles (Lechuga, Cilantro, Zanahoria)\n")
## - Verdura: 3 niveles (Lechuga, Cilantro, Zanahoria)
cat(" - Temperatura: 2 niveles (8°C, 20°C)\n")
## - Temperatura: 2 niveles (8°C, 20°C)
cat(" - Réplicas: 6 por combinación\n\n")
## - Réplicas: 6 por combinación
cat("Hipótesis que pueden probarse:\n")
## Hipótesis que pueden probarse:
cat("1. Efecto principal de verdura:\n")
## 1. Efecto principal de verdura:
cat(" H0: μ_Lechuga = μ_Cilantro = μ_Zanahoria\n")
## H0: μ_Lechuga = μ_Cilantro = μ_Zanahoria
cat(" H1: Al menos una media de verdura es diferente\n\n")
## H1: Al menos una media de verdura es diferente
cat("2. Efecto principal de temperatura:\n")
## 2. Efecto principal de temperatura:
cat(" H0: μ_8°C = μ_20°C\n")
## H0: μ_8°C = μ_20°C
cat(" H1: Las medias de temperatura son diferentes\n\n")
## H1: Las medias de temperatura son diferentes
cat("3. Efecto de interacción verdura×temperatura:\n")
## 3. Efecto de interacción verdura×temperatura:
cat(" H0: El efecto de temperatura es el mismo para todas las verduras\n")
## H0: El efecto de temperatura es el mismo para todas las verduras
cat(" H1: El efecto de temperatura difiere entre las verduras\n")
## H1: El efecto de temperatura difiere entre las verduras
# b) Análisis de varianza
cat("\n13b) Análisis de varianza:\n")
##
## 13b) Análisis de varianza:
modelo13 <- aov(sobrevivencia ~ verdura * temperatura, data = datos13)
anova13 <- anova(modelo13)
print(anova13)
## Analysis of Variance Table
##
## Response: sobrevivencia
## Df Sum Sq Mean Sq F value Pr(>F)
## verdura 2 325.07 162.54 5.6258 0.008418 **
## temperatura 1 2071.33 2071.33 71.6944 1.895e-09 ***
## verdura:temperatura 2 129.39 64.69 2.2392 0.124055
## Residuals 30 866.73 28.89
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nInterpretación detallada:\n")
##
## Interpretación detallada:
# Verdura
if(anova13["verdura", "Pr(>F)"] < 0.05) {
cat("1. Efecto de verdura: SIGNIFICATIVO (p =",
round(anova13["verdura", "Pr(>F)"], 4), ")\n")
cat(" - El tipo de verdura afecta la sobrevivencia del V. cholerae\n")
# Comparaciones múltiples para verdura
tukey_verdura <- TukeyHSD(modelo13, "verdura")
cat("\n Comparaciones de Tukey entre verduras:\n")
print(tukey_verdura$verdura)
# Calcular medias
medias_verdura <- aggregate(sobrevivencia ~ verdura, datos13, mean)
cat("\n Medias por verdura:\n")
print(medias_verdura)
mejor_verdura <- medias_verdura[which.max(medias_verdura$sobrevivencia), ]
cat(" Mejor verdura para sobrevivencia:", mejor_verdura$verdura,
"(", round(mejor_verdura$sobrevivencia, 1), "%)\n")
} else {
cat("1. Efecto de verdura: NO SIGNIFICATIVO (p =",
round(anova13["verdura", "Pr(>F)"], 4), ")\n")
cat(" - El tipo de verdura no afecta significativamente la sobrevivencia\n")
}
## 1. Efecto de verdura: SIGNIFICATIVO (p = 0.0084 )
## - El tipo de verdura afecta la sobrevivencia del V. cholerae
##
## Comparaciones de Tukey entre verduras:
## diff lwr upr p adj
## Lechuga-Cilantro -2.802483 -8.212154 2.6071892 0.418761472
## Zanahoria-Cilantro -7.295626 -12.705298 -1.8859543 0.006436573
## Zanahoria-Lechuga -4.493143 -9.902815 0.9165283 0.118331753
##
## Medias por verdura:
## verdura sobrevivencia
## 1 Cilantro 81.18519
## 2 Lechuga 78.38271
## 3 Zanahoria 73.88956
## Mejor verdura para sobrevivencia: 1 ( 81.2 %)
# Temperatura
if(anova13["temperatura", "Pr(>F)"] < 0.05) {
cat("\n2. Efecto de temperatura: SIGNIFICATIVO (p =",
round(anova13["temperatura", "Pr(>F)"], 4), ")\n")
cat(" - La temperatura afecta significativamente la sobrevivencia\n")
medias_temp <- aggregate(sobrevivencia ~ temperatura, datos13, mean)
cat("\n Medias por temperatura:\n")
print(medias_temp)
mejor_temp <- medias_temp[which.max(medias_temp$sobrevivencia), ]
cat(" Mejor temperatura para sobrevivencia:", mejor_temp$temperatura,
"(", round(mejor_temp$sobrevivencia, 1), "%)\n")
} else {
cat("\n2. Efecto de temperatura: NO SIGNIFICATIVO (p =",
round(anova13["temperatura", "Pr(>F)"], 4), ")\n")
}
##
## 2. Efecto de temperatura: SIGNIFICATIVO (p = 0 )
## - La temperatura afecta significativamente la sobrevivencia
##
## Medias por temperatura:
## temperatura sobrevivencia
## 1 20°C 70.23385
## 2 8°C 85.40446
## Mejor temperatura para sobrevivencia: 2 ( 85.4 %)
# Interacción
if(anova13["verdura:temperatura", "Pr(>F)"] < 0.05) {
cat("\n3. Interacción verdura×temperatura: SIGNIFICATIVA (p =",
round(anova13["verdura:temperatura", "Pr(>F)"], 4), ")\n")
cat(" - El efecto de la temperatura depende del tipo de verdura\n")
cat(" - O el efecto de la verdura depende de la temperatura\n")
# Gráfica de interacción
interaccion_data <- datos13 %>%
group_by(verdura, temperatura) %>%
summarise(
media = mean(sobrevivencia),
error = sd(sobrevivencia)/sqrt(n()),
.groups = 'drop'
)
p_interaccion <- ggplot(interaccion_data,
aes(x = verdura, y = media,
color = temperatura, group = temperatura)) +
geom_line(size = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = media - error, ymax = media + error),
width = 0.1, size = 0.8) +
labs(title = "Interacción Verdura × Temperatura",
x = "Tipo de Verdura",
y = "Sobrevivencia Promedio (%)",
color = "Temperatura") +
theme_minimal()
print(p_interaccion)
} else {
cat("\n3. Interacción verdura×temperatura: NO SIGNIFICATIVA (p =",
round(anova13["verdura:temperatura", "Pr(>F)"], 4), ")\n")
cat(" - Los efectos de verdura y temperatura son independientes\n")
}
##
## 3. Interacción verdura×temperatura: NO SIGNIFICATIVA (p = 0.1241 )
## - Los efectos de verdura y temperatura son independientes
# c) Verificación del supuesto de igual varianza
cat("\n13c) Verificación del supuesto de igual varianza (homocedasticidad):\n")
##
## 13c) Verificación del supuesto de igual varianza (homocedasticidad):
# 1. Gráfico de residuos vs valores ajustados
par(mfrow = c(2, 2))
plot(modelo13, main = "Diagnóstico del Modelo 13")
# 2. Prueba de Levene
cat("\nPrueba de Levene para homogeneidad de varianzas:\n")
##
## Prueba de Levene para homogeneidad de varianzas:
levene_test13 <- leveneTest(sobrevivencia ~ interaction(verdura, temperatura),
data = datos13)
cat(" Estadístico F =", round(levene_test13$`F value`[1], 4), "\n")
## Estadístico F = 1.5535
cat(" p-value =", round(levene_test13$`Pr(>F)`[1], 4))
## p-value = 0.2034
if(levene_test13$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
cat(" El supuesto se cumple satisfactoriamente.\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
cat(" El supuesto NO se cumple satisfactoriamente.\n")
}
## → No se rechaza homogeneidad de varianzas
## El supuesto se cumple satisfactoriamente.
# 3. Prueba de Bartlett (alternativa)
cat("\nPrueba de Bartlett (alternativa):\n")
##
## Prueba de Bartlett (alternativa):
bartlett_test13 <- bartlett.test(sobrevivencia ~ interaction(verdura, temperatura),
data = datos13)
cat(" Estadístico K² =", round(bartlett_test13$statistic, 4), "\n")
## Estadístico K² = 4.3892
cat(" p-value =", round(bartlett_test13$p.value, 4))
## p-value = 0.4948
if(bartlett_test13$p.value > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
}
## → No se rechaza homogeneidad de varianzas
# 4. Boxplot por grupo para inspección visual
cat("\nAnálisis visual por boxplots:\n")
##
## Análisis visual por boxplots:
p_boxplots <- ggplot(datos13,
aes(x = interaction(verdura, temperatura),
y = sobrevivencia,
fill = interaction(verdura, temperatura))) +
geom_boxplot() +
labs(title = "Distribución de sobrevivencia por combinación",
x = "Combinación Verdura × Temperatura",
y = "Sobrevivencia (%)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p_boxplots)
# d) Efecto si no se cumple el supuesto
cat("\n13d) Efecto si no se cumple el supuesto de homogeneidad de varianzas:\n")
##
## 13d) Efecto si no se cumple el supuesto de homogeneidad de varianzas:
if(levene_test13$`Pr(>F)`[1] < 0.05) {
cat("Consecuencias:\n")
cat("1. El nivel de significancia real puede diferir del nominal (0.05)\n")
cat("2. Las pruebas F pueden ser muy sensibles o muy insensibles\n")
cat("3. Los intervalos de confianza pueden tener cobertura incorrecta\n")
cat("4. Las comparaciones múltiples pueden ser inválidas\n\n")
cat("Alternativas cuando no se cumple homogeneidad:\n")
cat("1. Transformación de datos (log, raíz cuadrada, etc.)\n")
cat("2. Usar pruebas robustas (Welch ANOVA)\n")
cat("3. Usar modelos de varianzas desiguales\n")
cat("4. Métodos no paramétricos\n")
} else {
cat("Como el supuesto se cumple, las conclusiones del ANOVA son válidas.\n")
}
## Como el supuesto se cumple, las conclusiones del ANOVA son válidas.
EJERCICIO 14. Continuación del ejercicio 13 con transformación
logarítmica
cat("14. Análisis con transformación logarítmica\n")
## 14. Análisis con transformación logarítmica
# a) Transformar datos con logaritmos
cat("\n14a) Transformación logarítmica de los datos:\n")
##
## 14a) Transformación logarítmica de los datos:
datos13$log_sobrevivencia <- log(datos13$sobrevivencia)
cat("Se añadió la variable: log_sobrevivencia = log(sobrevivencia)\n")
## Se añadió la variable: log_sobrevivencia = log(sobrevivencia)
cat("Primeros valores transformados:\n")
## Primeros valores transformados:
print(head(datos13[, c("sobrevivencia", "log_sobrevivencia")]))
## sobrevivencia log_sobrevivencia
## 1 82.19762 4.409126
## 2 83.84911 4.429019
## 3 92.79354 4.530377
## 4 85.35254 4.446790
## 5 85.64644 4.450228
## 6 93.57532 4.538767
# b) Análisis de varianza con datos transformados
cat("\n14b) ANOVA con datos transformados:\n")
##
## 14b) ANOVA con datos transformados:
modelo14_log <- aov(log_sobrevivencia ~ verdura * temperatura, data = datos13)
anova14_log <- anova(modelo14_log)
print(anova14_log)
## Analysis of Variance Table
##
## Response: log_sobrevivencia
## Df Sum Sq Mean Sq F value Pr(>F)
## verdura 2 0.05027 0.02513 5.0610 0.01277 *
## temperatura 1 0.34030 0.34030 68.5246 3.065e-09 ***
## verdura:temperatura 2 0.01837 0.00919 1.8495 0.17480
## Residuals 30 0.14898 0.00497
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# c) Verificación de supuestos en datos transformados
cat("\n14c) Verificación de supuestos con datos transformados:\n")
##
## 14c) Verificación de supuestos con datos transformados:
# Gráficos de diagnóstico
par(mfrow = c(2, 2))
plot(modelo14_log, main = "Diagnóstico del Modelo (datos transformados)")
# Prueba de Levene en datos transformados
cat("\nPrueba de Levene en datos transformados:\n")
##
## Prueba de Levene en datos transformados:
levene_test14 <- leveneTest(log_sobrevivencia ~ interaction(verdura, temperatura),
data = datos13)
cat(" p-value =", round(levene_test14$`Pr(>F)`[1], 4))
## p-value = 0.1183
if(levene_test14$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
cat(" La transformación logarítmica mejoró la homocedasticidad.\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
cat(" La transformación no corrigió completamente el problema.\n")
}
## → No se rechaza homogeneidad de varianzas
## La transformación logarítmica mejoró la homocedasticidad.
# Prueba de normalidad en residuos transformados
cat("\nPrueba de normalidad en residuos transformados (Shapiro-Wilk):\n")
##
## Prueba de normalidad en residuos transformados (Shapiro-Wilk):
residuos_log <- residuals(modelo14_log)
shapiro_test14 <- shapiro.test(residuos_log)
cat(" p-value =", round(shapiro_test14$p.value, 4))
## p-value = 0.9998
if(shapiro_test14$p.value > 0.05) {
cat(" → No se rechaza normalidad\n")
} else {
cat(" → Se rechaza normalidad\n")
}
## → No se rechaza normalidad
# d) Diferencias importantes entre análisis
cat("\n14d) Diferencias importantes entre los dos análisis:\n")
##
## 14d) Diferencias importantes entre los dos análisis:
cat("\n1. Comparación de valores-p:\n")
##
## 1. Comparación de valores-p:
cat(" ANOVA original vs transformado:\n")
## ANOVA original vs transformado:
cat(" - Verdura: ", round(anova13["verdura", "Pr(>F)"], 4),
" vs ", round(anova14_log["verdura", "Pr(>F)"], 4), "\n")
## - Verdura: 0.0084 vs 0.0128
cat(" - Temperatura: ", round(anova13["temperatura", "Pr(>F)"], 4),
" vs ", round(anova14_log["temperatura", "Pr(>F)"], 4), "\n")
## - Temperatura: 0 vs 0
cat(" - Interacción: ", round(anova13["verdura:temperatura", "Pr(>F)"], 4),
" vs ", round(anova14_log["verdura:temperatura", "Pr(>F)"], 4), "\n")
## - Interacción: 0.1241 vs 0.1748
# Comparación de medias en escala original y transformada
cat("\n2. Comparación de medias:\n")
##
## 2. Comparación de medias:
medias_original <- datos13 %>%
group_by(verdura, temperatura) %>%
summarise(media_original = mean(sobrevivencia),
media_transformada = exp(mean(log_sobrevivencia)),
.groups = 'drop')
cat(" Medias en diferentes escalas:\n")
## Medias en diferentes escalas:
print(medias_original)
## # A tibble: 6 × 4
## verdura temperatura media_original media_transformada
## <fct> <fct> <dbl> <dbl>
## 1 Cilantro 20°C 72.2 72.1
## 2 Cilantro 8°C 90.2 90.1
## 3 Lechuga 20°C 69.5 69.2
## 4 Lechuga 8°C 87.2 87.1
## 5 Zanahoria 20°C 69.0 68.9
## 6 Zanahoria 8°C 78.8 78.5
# Explicación de la transformación
cat("\n3. Por qué usar transformación logarítmica:\n")
##
## 3. Por qué usar transformación logarítmica:
cat(" - Cuando los datos presentan heterocedasticidad (varianzas desiguales)\n")
## - Cuando los datos presentan heterocedasticidad (varianzas desiguales)
cat(" - Cuando los efectos son multiplicativos en lugar de aditivos\n")
## - Cuando los efectos son multiplicativos en lugar de aditivos
cat(" - Para estabilizar varianzas cuando hay relación media-varianza\n")
## - Para estabilizar varianzas cuando hay relación media-varianza
cat(" - Para cumplir con los supuestos del modelo lineal\n")
## - Para cumplir con los supuestos del modelo lineal
# e) Interpretación de interacción con datos transformados
cat("\n14e) Interpretación de interacción con datos transformados:\n")
##
## 14e) Interpretación de interacción con datos transformados:
if(anova14_log["verdura:temperatura", "Pr(>F)"] < 0.05) {
cat("Hay interacción significativa en escala logarítmica.\n")
# Gráfica de interacción en escala transformada
interaccion_log <- datos13 %>%
group_by(verdura, temperatura) %>%
summarise(
media_log = mean(log_sobrevivencia),
error_log = sd(log_sobrevivencia)/sqrt(n()),
.groups = 'drop'
)
p_inter_log <- ggplot(interaccion_log,
aes(x = verdura, y = media_log,
color = temperatura, group = temperatura)) +
geom_line(size = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = media_log - error_log,
ymax = media_log + error_log),
width = 0.1, size = 0.8) +
labs(title = "Interacción en Escala Logarítmica",
x = "Tipo de Verdura",
y = "log(Sobrevivencia)",
color = "Temperatura") +
theme_minimal()
print(p_inter_log)
cat("\nInterpretación:\n")
cat("En escala logarítmica, una interacción significa que la razón\n")
cat("de sobrevivencia entre temperaturas difiere entre las verduras.\n")
cat("Ejemplo: si la razón de sobrevivencia 8°C/20°C es diferente\n")
cat("para lechuga vs cilantro vs zanahoria.\n")
} else {
cat("No hay interacción significativa en escala logarítmica.\n")
cat("Los efectos son aditivos en la escala logarítmica.\n")
}
## No hay interacción significativa en escala logarítmica.
## Los efectos son aditivos en la escala logarítmica.
# f) Comparación gráfica de ambos análisis
cat("\nComparación gráfica de ambos análisis:\n")
##
## Comparación gráfica de ambos análisis:
# Crear gráficos comparativos
p_comparativo1 <- ggplot(datos13, aes(x = sobrevivencia)) +
geom_histogram(binwidth = 5, fill = "lightblue", color = "black", alpha = 0.7) +
labs(title = "Distribución Original", x = "Sobrevivencia (%)", y = "Frecuencia") +
theme_minimal()
p_comparativo2 <- ggplot(datos13, aes(x = log_sobrevivencia)) +
geom_histogram(binwidth = 0.1, fill = "lightgreen", color = "black", alpha = 0.7) +
labs(title = "Distribución Transformada", x = "log(Sobrevivencia)", y = "Frecuencia") +
theme_minimal()
# Boxplots comparativos
p_comparativo3 <- ggplot(datos13,
aes(x = interaction(verdura, temperatura),
y = sobrevivencia,
fill = interaction(verdura, temperatura))) +
geom_boxplot() +
labs(title = "Datos Originales",
x = "Combinación",
y = "Sobrevivencia (%)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none")
p_comparativo4 <- ggplot(datos13,
aes(x = interaction(verdura, temperatura),
y = log_sobrevivencia,
fill = interaction(verdura, temperatura))) +
geom_boxplot() +
labs(title = "Datos Transformados",
x = "Combinación",
y = "log(Sobrevivencia)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none")
grid.arrange(p_comparativo1, p_comparativo2, p_comparativo3, p_comparativo4,
ncol = 2, nrow = 2)
# Continuación del ejercicio 13 con transformación logarítmica
cat("14. Análisis con transformación logarítmica\n")
## 14. Análisis con transformación logarítmica
# a) Transformar datos con logaritmos
cat("\n14a) Transformación logarítmica de los datos:\n")
##
## 14a) Transformación logarítmica de los datos:
datos13$log_sobrevivencia <- log(datos13$sobrevivencia)
cat("Se añadió la variable: log_sobrevivencia = log(sobrevivencia)\n")
## Se añadió la variable: log_sobrevivencia = log(sobrevivencia)
cat("Primeros valores transformados:\n")
## Primeros valores transformados:
print(head(datos13[, c("sobrevivencia", "log_sobrevivencia")]))
## sobrevivencia log_sobrevivencia
## 1 82.19762 4.409126
## 2 83.84911 4.429019
## 3 92.79354 4.530377
## 4 85.35254 4.446790
## 5 85.64644 4.450228
## 6 93.57532 4.538767
# b) Análisis de varianza con datos transformados
cat("\n14b) ANOVA con datos transformados:\n")
##
## 14b) ANOVA con datos transformados:
modelo14_log <- aov(log_sobrevivencia ~ verdura * temperatura, data = datos13)
anova14_log <- anova(modelo14_log)
print(anova14_log)
## Analysis of Variance Table
##
## Response: log_sobrevivencia
## Df Sum Sq Mean Sq F value Pr(>F)
## verdura 2 0.05027 0.02513 5.0610 0.01277 *
## temperatura 1 0.34030 0.34030 68.5246 3.065e-09 ***
## verdura:temperatura 2 0.01837 0.00919 1.8495 0.17480
## Residuals 30 0.14898 0.00497
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# c) Verificación de supuestos en datos transformados
cat("\n14c) Verificación de supuestos con datos transformados:\n")
##
## 14c) Verificación de supuestos con datos transformados:
# Gráficos de diagnóstico
par(mfrow = c(2, 2))
plot(modelo14_log, main = "Diagnóstico del Modelo (datos transformados)")
# Prueba de Levene en datos transformados
cat("\nPrueba de Levene en datos transformados:\n")
##
## Prueba de Levene en datos transformados:
levene_test14 <- leveneTest(log_sobrevivencia ~ interaction(verdura, temperatura),
data = datos13)
cat(" p-value =", round(levene_test14$`Pr(>F)`[1], 4))
## p-value = 0.1183
if(levene_test14$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas\n")
cat(" La transformación logarítmica mejoró la homocedasticidad.\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas\n")
cat(" La transformación no corrigió completamente el problema.\n")
}
## → No se rechaza homogeneidad de varianzas
## La transformación logarítmica mejoró la homocedasticidad.
# Prueba de normalidad en residuos transformados
cat("\nPrueba de normalidad en residuos transformados (Shapiro-Wilk):\n")
##
## Prueba de normalidad en residuos transformados (Shapiro-Wilk):
residuos_log <- residuals(modelo14_log)
shapiro_test14 <- shapiro.test(residuos_log)
cat(" p-value =", round(shapiro_test14$p.value, 4))
## p-value = 0.9998
if(shapiro_test14$p.value > 0.05) {
cat(" → No se rechaza normalidad\n")
} else {
cat(" → Se rechaza normalidad\n")
}
## → No se rechaza normalidad
# d) Diferencias importantes entre análisis
cat("\n14d) Diferencias importantes entre los dos análisis:\n")
##
## 14d) Diferencias importantes entre los dos análisis:
cat("\n1. Comparación de valores-p:\n")
##
## 1. Comparación de valores-p:
cat(" ANOVA original vs transformado:\n")
## ANOVA original vs transformado:
cat(" - Verdura: ", round(anova13["verdura", "Pr(>F)"], 4),
" vs ", round(anova14_log["verdura", "Pr(>F)"], 4), "\n")
## - Verdura: 0.0084 vs 0.0128
cat(" - Temperatura: ", round(anova13["temperatura", "Pr(>F)"], 4),
" vs ", round(anova14_log["temperatura", "Pr(>F)"], 4), "\n")
## - Temperatura: 0 vs 0
cat(" - Interacción: ", round(anova13["verdura:temperatura", "Pr(>F)"], 4),
" vs ", round(anova14_log["verdura:temperatura", "Pr(>F)"], 4), "\n")
## - Interacción: 0.1241 vs 0.1748
# Comparación de medias en escala original y transformada
cat("\n2. Comparación de medias:\n")
##
## 2. Comparación de medias:
medias_original <- datos13 %>%
group_by(verdura, temperatura) %>%
summarise(media_original = mean(sobrevivencia),
media_transformada = exp(mean(log_sobrevivencia)),
.groups = 'drop')
cat(" Medias en diferentes escalas:\n")
## Medias en diferentes escalas:
print(medias_original)
## # A tibble: 6 × 4
## verdura temperatura media_original media_transformada
## <fct> <fct> <dbl> <dbl>
## 1 Cilantro 20°C 72.2 72.1
## 2 Cilantro 8°C 90.2 90.1
## 3 Lechuga 20°C 69.5 69.2
## 4 Lechuga 8°C 87.2 87.1
## 5 Zanahoria 20°C 69.0 68.9
## 6 Zanahoria 8°C 78.8 78.5
# Explicación de la transformación
cat("\n3. Por qué usar transformación logarítmica:\n")
##
## 3. Por qué usar transformación logarítmica:
cat(" - Cuando los datos presentan heterocedasticidad (varianzas desiguales)\n")
## - Cuando los datos presentan heterocedasticidad (varianzas desiguales)
cat(" - Cuando los efectos son multiplicativos en lugar de aditivos\n")
## - Cuando los efectos son multiplicativos en lugar de aditivos
cat(" - Para estabilizar varianzas cuando hay relación media-varianza\n")
## - Para estabilizar varianzas cuando hay relación media-varianza
cat(" - Para cumplir con los supuestos del modelo lineal\n")
## - Para cumplir con los supuestos del modelo lineal
# e) Interpretación de interacción con datos transformados
cat("\n14e) Interpretación de interacción con datos transformados:\n")
##
## 14e) Interpretación de interacción con datos transformados:
if(anova14_log["verdura:temperatura", "Pr(>F)"] < 0.05) {
cat("Hay interacción significativa en escala logarítmica.\n")
# Gráfica de interacción en escala transformada
interaccion_log <- datos13 %>%
group_by(verdura, temperatura) %>%
summarise(
media_log = mean(log_sobrevivencia),
error_log = sd(log_sobrevivencia)/sqrt(n()),
.groups = 'drop'
)
p_inter_log <- ggplot(interaccion_log,
aes(x = verdura, y = media_log,
color = temperatura, group = temperatura)) +
geom_line(size = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = media_log - error_log,
ymax = media_log + error_log),
width = 0.1, size = 0.8) +
labs(title = "Interacción en Escala Logarítmica",
x = "Tipo de Verdura",
y = "log(Sobrevivencia)",
color = "Temperatura") +
theme_minimal()
print(p_inter_log)
cat("\nInterpretación:\n")
cat("En escala logarítmica, una interacción significa que la razón\n")
cat("de sobrevivencia entre temperaturas difiere entre las verduras.\n")
cat("Ejemplo: si la razón de sobrevivencia 8°C/20°C es diferente\n")
cat("para lechuga vs cilantro vs zanahoria.\n")
} else {
cat("No hay interacción significativa en escala logarítmica.\n")
cat("Los efectos son aditivos en la escala logarítmica.\n")
}
## No hay interacción significativa en escala logarítmica.
## Los efectos son aditivos en la escala logarítmica.
# f) Comparación gráfica de ambos análisis
cat("\nComparación gráfica de ambos análisis:\n")
##
## Comparación gráfica de ambos análisis:
# Crear gráficos comparativos
p_comparativo1 <- ggplot(datos13, aes(x = sobrevivencia)) +
geom_histogram(binwidth = 5, fill = "lightblue", color = "black", alpha = 0.7) +
labs(title = "Distribución Original", x = "Sobrevivencia (%)", y = "Frecuencia") +
theme_minimal()
p_comparativo2 <- ggplot(datos13, aes(x = log_sobrevivencia)) +
geom_histogram(binwidth = 0.1, fill = "lightgreen", color = "black", alpha = 0.7) +
labs(title = "Distribución Transformada", x = "log(Sobrevivencia)", y = "Frecuencia") +
theme_minimal()
# Boxplots comparativos
p_comparativo3 <- ggplot(datos13,
aes(x = interaction(verdura, temperatura),
y = sobrevivencia,
fill = interaction(verdura, temperatura))) +
geom_boxplot() +
labs(title = "Datos Originales",
x = "Combinación",
y = "Sobrevivencia (%)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none")
p_comparativo4 <- ggplot(datos13,
aes(x = interaction(verdura, temperatura),
y = log_sobrevivencia,
fill = interaction(verdura, temperatura))) +
geom_boxplot() +
labs(title = "Datos Transformados",
x = "Combinación",
y = "log(Sobrevivencia)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none")
grid.arrange(p_comparativo1, p_comparativo2, p_comparativo3, p_comparativo4,
ncol = 2, nrow = 2)
Ejercicio 15 - Evaluación de Antioxidantes
# Datos ejercicio 15
producto <- factor(rep(c("Control", "A", "B", "C", "D"), each = 6))
tiempo <- factor(rep(rep(c("4h", "8h", "12h"), each = 2), 5))
# Índice de peróxidos
indice_peroxidos <- c(
# Control
3.84, 3.72, # 4 horas
27.63, 27.58, # 8 horas
39.95, 39.00, # 12 horas
# Producto A
4.00, 3.91, # 4 horas
22.00, 21.83, # 8 horas
46.20, 45.60, # 12 horas
# Producto B
3.61, 3.61, # 4 horas
21.94, 21.85, # 8 horas
46.58, 42.98, # 12 horas
# Producto C
3.57, 3.50, # 4 horas
20.50, 20.32, # 8 horas
45.14, 44.89, # 12 horas
# Producto D
3.64, 3.61, # 4 horas
20.30, 20.19, # 8 horas
44.36, 44.02 # 12 horas
)
datos15_real <- data.frame(producto, tiempo, indice_peroxidos)
# a) Factores controlados y variable de respuesta
cat("15a) Factores y variable de respuesta:\n")
## 15a) Factores y variable de respuesta:
cat("Factores controlados:\n")
## Factores controlados:
cat("1. Producto (antioxidante): 5 niveles (Control, A, B, C, D)\n")
## 1. Producto (antioxidante): 5 niveles (Control, A, B, C, D)
cat("2. Tiempo: 3 niveles (4h, 8h, 12h) - factor de bloque\n\n")
## 2. Tiempo: 3 niveles (4h, 8h, 12h) - factor de bloque
cat("Variable de respuesta:\n")
## Variable de respuesta:
cat("- Índice de peróxidos (medida de oxidación, menor es mejor)\n")
## - Índice de peróxidos (medida de oxidación, menor es mejor)
cat("- Unidades: No especificadas, probablemente meq/kg o similar\n")
## - Unidades: No especificadas, probablemente meq/kg o similar
cat("- Réplicas: 2 por combinación\n")
## - Réplicas: 2 por combinación
# b) Modelo estadístico y hipótesis
cat("\n15b) Modelo estadístico y hipótesis:\n")
##
## 15b) Modelo estadístico y hipótesis:
cat("Modelo: Y_ijk = μ + α_i + β_j + (αβ)_ij + ε_ijk\n")
## Modelo: Y_ijk = μ + α_i + β_j + (αβ)_ij + ε_ijk
cat("donde:\n")
## donde:
cat(" α_i: efecto del producto (i=1,...,5)\n")
## α_i: efecto del producto (i=1,...,5)
cat(" β_j: efecto del tiempo (j=1,2,3)\n")
## β_j: efecto del tiempo (j=1,2,3)
cat(" (αβ)_ij: interacción producto×tiempo\n")
## (αβ)_ij: interacción producto×tiempo
cat(" ε_ijk: error experimental\n\n")
## ε_ijk: error experimental
cat("Hipótesis principales:\n")
## Hipótesis principales:
cat("1. H0: α_Control = α_A = α_B = α_C = α_D = 0\n")
## 1. H0: α_Control = α_A = α_B = α_C = α_D = 0
cat(" H1: Al menos un producto tiene efecto diferente\n\n")
## H1: Al menos un producto tiene efecto diferente
cat("2. H0: β_4h = β_8h = β_12h = 0\n")
## 2. H0: β_4h = β_8h = β_12h = 0
cat(" H1: El tiempo afecta el índice de peróxidos\n\n")
## H1: El tiempo afecta el índice de peróxidos
cat("3. H0: (αβ)_ij = 0 para todo i,j\n")
## 3. H0: (αβ)_ij = 0 para todo i,j
cat(" H1: Existe interacción producto×tiempo\n")
## H1: Existe interacción producto×tiempo
# c) Análisis de varianza
cat("\n15c) Análisis de varianza:\n")
##
## 15c) Análisis de varianza:
modelo15_real <- aov(indice_peroxidos ~ producto * tiempo, data = datos15_real)
anova15_real <- anova(modelo15_real)
print(anova15_real)
## Analysis of Variance Table
##
## Response: indice_peroxidos
## Df Sum Sq Mean Sq F value Pr(>F)
## producto 4 5.9 1.5 3.0311 0.05119 .
## tiempo 2 8081.1 4040.6 8352.3327 < 2.2e-16 ***
## producto:tiempo 8 118.1 14.8 30.5126 6.333e-08 ***
## Residuals 15 7.3 0.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Resumen estadístico descriptivo
cat("\nEstadística descriptiva por producto:\n")
##
## Estadística descriptiva por producto:
estadisticas_producto <- datos15_real %>%
group_by(producto) %>%
summarise(
n = n(),
media = mean(indice_peroxidos),
desviacion = sd(indice_peroxidos),
min = min(indice_peroxidos),
max = max(indice_peroxidos),
.groups = 'drop'
)
print(estadisticas_producto)
## # A tibble: 5 × 6
## producto n media desviacion min max
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 A 6 23.9 18.8 3.91 46.2
## 2 B 6 23.4 18.5 3.61 46.6
## 3 C 6 23.0 18.7 3.5 45.1
## 4 Control 6 23.6 16.3 3.72 40.0
## 5 D 6 22.7 18.2 3.61 44.4
cat("\nEstadística descriptiva por tiempo:\n")
##
## Estadística descriptiva por tiempo:
estadisticas_tiempo <- datos15_real %>%
group_by(tiempo) %>%
summarise(
n = n(),
media = mean(indice_peroxidos),
desviacion = sd(indice_peroxidos),
min = min(indice_peroxidos),
max = max(indice_peroxidos),
.groups = 'drop'
)
print(estadisticas_tiempo)
## # A tibble: 3 × 6
## tiempo n media desviacion min max
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 12h 10 43.9 2.55 39 46.6
## 2 4h 10 3.70 0.163 3.5 4
## 3 8h 10 22.4 2.84 20.2 27.6
# d) Aspectos más relevantes del ANOVA
cat("\n15c) Aspectos más relevantes del análisis:\n")
##
## 15c) Aspectos más relevantes del análisis:
# 1. Efecto de producto
cat("\n1. Efecto del producto (antioxidante):\n")
##
## 1. Efecto del producto (antioxidante):
if(anova15_real["producto", "Pr(>F)"] < 0.05) {
cat(" SIGNIFICATIVO (p =", format(anova15_real["producto", "Pr(>F)"], scientific = TRUE), ")\n")
cat(" - Los productos difieren en su capacidad antioxidante\n")
# Comparaciones múltiples de Tukey
cat("\n Comparaciones de Tukey entre productos:\n")
tukey_producto <- TukeyHSD(modelo15_real, "producto")
print(tukey_producto$producto)
# Identificar grupos homogéneos
cat("\n Grupos homogéneos (LSD test):\n")
lsd_producto <- LSD.test(modelo15_real, "producto", console = FALSE)
print(lsd_producto$groups)
} else {
cat(" NO SIGNIFICATIVO (p =", round(anova15_real["producto", "Pr(>F)"], 4), ")\n")
cat(" - No hay evidencia de diferencias entre productos\n")
}
## NO SIGNIFICATIVO (p = 0.0512 )
## - No hay evidencia de diferencias entre productos
# 2. Efecto de tiempo
cat("\n2. Efecto del tiempo:\n")
##
## 2. Efecto del tiempo:
if(anova15_real["tiempo", "Pr(>F)"] < 0.05) {
cat(" SIGNIFICATIVO (p =", format(anova15_real["tiempo", "Pr(>F)"], scientific = TRUE), ")\n")
cat(" - El índice de peróxidos cambia significativamente con el tiempo\n")
# Comparaciones múltiples para tiempo
cat("\n Comparaciones de Tukey entre tiempos:\n")
tukey_tiempo <- TukeyHSD(modelo15_real, "tiempo")
print(tukey_tiempo$tiempo)
# Análisis de tendencia
cat("\n Tendencia del índice de peróxidos con el tiempo:\n")
medias_tiempo <- aggregate(indice_peroxidos ~ tiempo, datos15_real, mean)
cat(" Medias: 4h =", round(medias_tiempo$indice_peroxidos[1], 2),
", 8h =", round(medias_tiempo$indice_peroxidos[2], 2),
", 12h =", round(medias_tiempo$indice_peroxidos[3], 2), "\n")
cat(" Incremento promedio por hora:",
round((medias_tiempo$indice_peroxidos[3] - medias_tiempo$indice_peroxidos[1])/8, 2), "unidades/hora\n")
} else {
cat(" NO SIGNIFICATIVO (p =", round(anova15_real["tiempo", "Pr(>F)"], 4), ")\n")
}
## SIGNIFICATIVO (p = 1.401127e-23 )
## - El índice de peróxidos cambia significativamente con el tiempo
##
## Comparaciones de Tukey entre tiempos:
## diff lwr upr p adj
## 4h-12h -40.171 -40.97895 -39.36305 0
## 8h-12h -21.458 -22.26595 -20.65005 0
## 8h-4h 18.713 17.90505 19.52095 0
##
## Tendencia del índice de peróxidos con el tiempo:
## Medias: 4h = 43.87 , 8h = 3.7 , 12h = 22.41
## Incremento promedio por hora: -2.68 unidades/hora
# 3. Interacción producto×tiempo
cat("\n3. Interacción producto × tiempo:\n")
##
## 3. Interacción producto × tiempo:
if(anova15_real["producto:tiempo", "Pr(>F)"] < 0.05) {
cat(" SIGNIFICATIVA (p =", round(anova15_real["producto:tiempo", "Pr(>F)"], 4), ")\n")
cat(" - El efecto del tiempo depende del producto\n")
cat(" - La eficacia de los antioxidantes cambia con el tiempo\n")
# Análisis de interacción simple
cat("\n Análisis de efectos simples:\n")
# Efecto del tiempo para cada producto
cat("\n Efecto del tiempo por producto (ANOVA por separado):\n")
productos <- unique(datos15_real$producto)
for(prod in productos) {
datos_prod <- datos15_real[datos15_real$producto == prod, ]
modelo_prod <- aov(indice_peroxidos ~ tiempo, data = datos_prod)
p_valor <- summary(modelo_prod)[[1]][1, "Pr(>F)"]
cat(" ", prod, ": p =", round(p_valor, 4),
ifelse(p_valor < 0.05, "(significativo)", "(no significativo)"), "\n")
}
# Efecto del producto para cada tiempo
cat("\n Efecto del producto por tiempo:\n")
tiempos <- unique(datos15_real$tiempo)
for(tiem in tiempos) {
datos_tiem <- datos15_real[datos15_real$tiempo == tiem, ]
modelo_tiem <- aov(indice_peroxidos ~ producto, data = datos_tiem)
p_valor <- summary(modelo_tiem)[[1]][1, "Pr(>F)"]
cat(" ", tiem, ": p =", round(p_valor, 4),
ifelse(p_valor < 0.05, "(significativo)", "(no significativo)"), "\n")
}
} else {
cat(" NO SIGNIFICATIVA (p =", round(anova15_real["producto:tiempo", "Pr(>F)"], 4), ")\n")
cat(" - Los efectos de producto y tiempo son aditivos\n")
cat(" - La eficacia relativa de los antioxidantes es constante en el tiempo\n")
}
## SIGNIFICATIVA (p = 0 )
## - El efecto del tiempo depende del producto
## - La eficacia de los antioxidantes cambia con el tiempo
##
## Análisis de efectos simples:
##
## Efecto del tiempo por producto (ANOVA por separado):
## Control : p = 0 (significativo)
## A : p = 0 (significativo)
## B : p = 2e-04 (significativo)
## C : p = 0 (significativo)
## D : p = 0 (significativo)
##
## Efecto del producto por tiempo:
## 4h : p = 0.0029 (significativo)
## 8h : p = 0 (significativo)
## 12h : p = 0.0169 (significativo)
# e) Verificación de supuestos
cat("\n15d) Verificación de supuestos del modelo:\n")
##
## 15d) Verificación de supuestos del modelo:
# 1. Gráficos de diagnóstico
cat("\n1. Gráficos de diagnóstico:\n")
##
## 1. Gráficos de diagnóstico:
par(mfrow = c(2, 2))
plot(modelo15_real, main = "Diagnóstico del Modelo")
# 2. Pruebas formales
residuos15 <- residuals(modelo15_real)
# Normalidad (Shapiro-Wilk)
cat("\n2. Prueba de normalidad (Shapiro-Wilk):\n")
##
## 2. Prueba de normalidad (Shapiro-Wilk):
shapiro_test15 <- shapiro.test(residuos15)
cat(" W =", round(shapiro_test15$statistic, 4), "\n")
## W = 0.6412
cat(" p-value =", round(shapiro_test15$p.value, 4))
## p-value = 0
if(shapiro_test15$p.value > 0.05) {
cat(" → No se rechaza normalidad ✓\n")
} else {
cat(" → Se rechaza normalidad ✗\n")
}
## → Se rechaza normalidad ✗
# Homogeneidad de varianzas (Levene)
cat("\n3. Prueba de homogeneidad de varianzas (Levene):\n")
##
## 3. Prueba de homogeneidad de varianzas (Levene):
levene_test15 <- leveneTest(indice_peroxidos ~ interaction(producto, tiempo),
data = datos15_real)
cat(" F =", round(levene_test15$`F value`[1], 4), "\n")
## F = 6.879988e+28
cat(" p-value =", round(levene_test15$`Pr(>F)`[1], 4))
## p-value = 0
if(levene_test15$`Pr(>F)`[1] > 0.05) {
cat(" → No se rechaza homogeneidad de varianzas ✓\n")
} else {
cat(" → Se rechaza homogeneidad de varianzas ✗\n")
}
## → Se rechaza homogeneidad de varianzas ✗
# 3. Gráficos adicionales para verificar supuestos
cat("\n4. Gráficos adicionales de diagnóstico:\n")
##
## 4. Gráficos adicionales de diagnóstico:
# Q-Q plot personalizado
qq_data <- data.frame(
teóricos = qqnorm(residuos15, plot.it = FALSE)$x,
muestrales = qqnorm(residuos15, plot.it = FALSE)$y
)
p_qq <- ggplot(qq_data, aes(x = teóricos, y = muestrales)) +
geom_point(size = 2) +
geom_abline(intercept = mean(residuos15), slope = sd(residuos15),
color = "red", size = 1) +
labs(title = "Q-Q Plot de Residuos",
x = "Cuantiles Teóricos",
y = "Cuantiles Muestrales") +
theme_minimal()
# Residuos vs valores ajustados
p_res_fitted <- ggplot(data.frame(
fitted = fitted(modelo15_real),
residuals = resid(modelo15_real)
), aes(x = fitted, y = residuals)) +
geom_point(size = 2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_smooth(method = "loess", se = FALSE, color = "blue") +
labs(title = "Residuos vs Valores Ajustados",
x = "Valores Ajustados",
y = "Residuos") +
theme_minimal()
grid.arrange(p_qq, p_res_fitted, ncol = 2)
# f) Identificación del mejor antioxidante
cat("\n15e) Identificación del mejor antioxidante:\n")
##
## 15e) Identificación del mejor antioxidante:
cat("Criterio: Menor índice de peróxidos es mejor (menor oxidación)\n")
## Criterio: Menor índice de peróxidos es mejor (menor oxidación)
# Calcular medias generales por producto
medias_generales <- aggregate(indice_peroxidos ~ producto, datos15_real, mean)
medias_generales <- medias_generales[order(medias_generales$indice_peroxidos), ]
cat("\nRanking de productos (de mejor a peor):\n")
##
## Ranking de productos (de mejor a peor):
for(i in 1:nrow(medias_generales)) {
cat(i, ". ", medias_generales$producto[i], ": ",
round(medias_generales$indice_peroxidos[i], 2), "\n", sep = "")
}
## 1. 5: 22.69
## 2. 3: 22.99
## 3. 2: 23.43
## 4. 4: 23.62
## 5. 1: 23.92
mejor_producto <- medias_generales[1, ]
peor_producto <- medias_generales[nrow(medias_generales), ]
cat("\nMejor producto:", mejor_producto$producto,
"(índice promedio =", round(mejor_producto$indice_peroxidos, 2), ")\n")
##
## Mejor producto: 5 (índice promedio = 22.69 )
cat("Peor producto:", peor_producto$producto,
"(índice promedio =", round(peor_producto$indice_peroxidos, 2), ")\n")
## Peor producto: 1 (índice promedio = 23.92 )
# Comparación con control
cat("\nComparación con el control:\n")
##
## Comparación con el control:
media_control <- medias_generales$indice_peroxidos[medias_generales$producto == "Control"]
cat("Control (sin antioxidante):", round(media_control, 2), "\n")
## Control (sin antioxidante): 23.62
for(prod in c("A", "B", "C", "D")) {
media_prod <- medias_generales$indice_peroxidos[medias_generales$producto == prod]
diferencia <- media_control - media_prod
porcentaje <- (diferencia / media_control) * 100
cat(prod, ": ", round(media_prod, 2),
" (diferencia = ", round(diferencia, 2),
", ", round(porcentaje, 1), "% vs control)\n", sep = "")
}
## A: 23.92 (diferencia = -0.3, -1.3% vs control)
## B: 23.43 (diferencia = 0.19, 0.8% vs control)
## C: 22.99 (diferencia = 0.63, 2.7% vs control)
## D: 22.69 (diferencia = 0.93, 4% vs control)
# Verificar significancia estadística vs control
if(anova15_real["producto", "Pr(>F)"] < 0.05) {
cat("\nComparaciones estadísticas específicas vs Control:\n")
# Extraer comparaciones del Tukey
comparaciones_tukey <- as.data.frame(tukey_producto$producto)
comparaciones_vs_control <- comparaciones_tukey[
grepl("Control", rownames(comparaciones_tukey)), ]
print(comparaciones_vs_control)
cat("\nProductos significativamente diferentes al Control (p < 0.05):\n")
productos_signif <- rownames(comparaciones_vs_control)[
comparaciones_vs_control$`p adj` < 0.05]
if(length(productos_signif) > 0) {
for(prod_sig in productos_signif) {
# Extraer nombre del producto
prod_nombre <- strsplit(prod_sig, "-")[[1]]
prod_nombre <- prod_nombre[prod_nombre != "Control"][1]
cat("- ", prod_nombre, " (p = ",
round(comparaciones_vs_control[prod_sig, "p adj"], 4), ")\n", sep = "")
}
} else {
cat("Ningún producto es significativamente diferente al Control\n")
}
}
# g) Gráficas de resultados
cat("\nVisualización de resultados:\n")
##
## Visualización de resultados:
# 1. Gráfica de medias por producto
p_medias_producto <- ggplot(medias_generales,
aes(x = reorder(producto, indice_peroxidos),
y = indice_peroxidos,
fill = producto)) +
geom_bar(stat = "identity") +
labs(title = "Índice de Peróxidos por Producto",
subtitle = "Menor valor indica mejor antioxidante",
x = "Producto",
y = "Índice de Peróxidos Promedio") +
theme_minimal() +
theme(legend.position = "none") +
geom_text(aes(label = round(indice_peroxidos, 1)),
vjust = -0.5, size = 3)
# 2. Gráfica de evolución temporal por producto
datos_evolucion <- datos15_real %>%
group_by(producto, tiempo) %>%
summarise(
media = mean(indice_peroxidos),
error = sd(indice_peroxidos)/sqrt(n()),
.groups = 'drop'
)
# Convertir tiempo a numérico para la gráfica
datos_evolucion$tiempo_num <- as.numeric(gsub("h", "", datos_evolucion$tiempo))
p_evolucion <- ggplot(datos_evolucion,
aes(x = tiempo_num, y = media,
color = producto, group = producto)) +
geom_line(size = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = media - error, ymax = media + error),
width = 0.5, size = 0.8) +
labs(title = "Evolución del Índice de Peróxidos",
subtitle = "Por producto y tiempo de almacenamiento",
x = "Tiempo (horas)",
y = "Índice de Peróxidos",
color = "Producto") +
scale_x_continuous(breaks = c(4, 8, 12), labels = c("4h", "8h", "12h")) +
theme_minimal()
# 3. Gráfica de interacción (si es significativa)
if(anova15_real["producto:tiempo", "Pr(>F)"] < 0.05) {
p_interaccion <- ggplot(datos_evolucion,
aes(x = tiempo, y = media,
color = producto, group = producto)) +
geom_line(size = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = media - error, ymax = media + error),
width = 0.1, size = 0.8) +
labs(title = "Interacción Producto × Tiempo",
x = "Tiempo",
y = "Índice de Peróxidos",
color = "Producto") +
theme_minimal()
grid.arrange(p_medias_producto, p_evolucion, p_interaccion, ncol = 3)
} else {
grid.arrange(p_medias_producto, p_evolucion, ncol = 2)
}
# 4. Heatmap de medias por combinación
cat("\nHeatmap de resultados:\n")
##
## Heatmap de resultados:
heatmap_data <- datos15_real %>%
group_by(producto, tiempo) %>%
summarise(media = mean(indice_peroxidos), .groups = 'drop')
p_heatmap <- ggplot(heatmap_data, aes(x = tiempo, y = producto, fill = media)) +
geom_tile(color = "white", size = 1) +
scale_fill_gradient(low = "green", high = "red",
name = "Índice de\nPeróxidos") +
geom_text(aes(label = round(media, 1)), color = "black", size = 4) +
labs(title = "Mapa de Calor: Índice de Peróxidos",
subtitle = "Verde = mejor (menor oxidación), Rojo = peor",
x = "Tiempo",
y = "Producto") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
print(p_heatmap)
# h) Conclusiones y recomendaciones
cat("\nCONCLUSIONES Y RECOMENDACIONES:\n")
##
## CONCLUSIONES Y RECOMENDACIONES:
cat("================================\n\n")
## ================================
cat("1. Mejor antioxidante:\n")
## 1. Mejor antioxidante:
cat(" Producto ", mejor_producto$producto, " con índice promedio de ",
round(mejor_producto$indice_peroxidos, 2), "\n", sep = "")
## Producto 5 con índice promedio de 22.69
cat(" Reducción vs control: ",
round((media_control - mejor_producto$indice_peroxidos)/media_control * 100, 1), "%\n\n", sep = "")
## Reducción vs control: 4%
cat("2. Efecto del tiempo:\n")
## 2. Efecto del tiempo:
cat(" El índice de peróxidos aumenta significativamente con el tiempo\n")
## El índice de peróxidos aumenta significativamente con el tiempo
cat(" Esto era esperado y valida el modelo experimental\n\n")
## Esto era esperado y valida el modelo experimental
cat("3. Recomendaciones:\n")
## 3. Recomendaciones:
cat(" a) Usar Producto ", mejor_producto$producto, " como antioxidante principal\n", sep = "")
## a) Usar Producto 5 como antioxidante principal
cat(" b) Considerar el costo-beneficio de cada antioxidante\n")
## b) Considerar el costo-beneficio de cada antioxidante
cat(" c) Validar resultados a mayor escala industrial\n")
## c) Validar resultados a mayor escala industrial
cat(" d) Estudiar interacciones con otros componentes del aceite\n\n")
## d) Estudiar interacciones con otros componentes del aceite
cat("4. Limitaciones del estudio:\n")
## 4. Limitaciones del estudio:
cat(" a) Solo 2 réplicas por combinación\n")
## a) Solo 2 réplicas por combinación
cat(" b) Solo 3 tiempos de evaluación\n")
## b) Solo 3 tiempos de evaluación
cat(" c) Condiciones de estrés pueden no representar condiciones reales\n")
## c) Condiciones de estrés pueden no representar condiciones reales
cat("16. Búsqueda de artículo científico con diseño factorial\n")
## 16. Búsqueda de artículo científico con diseño factorial
cat("\nPara este ejercicio, se debe buscar manualmente un artículo en scholar.google.com\n")
##
## Para este ejercicio, se debe buscar manualmente un artículo en scholar.google.com
cat("que cumpla con los siguientes criterios:\n\n")
## que cumpla con los siguientes criterios:
cat("Criterios del artículo a buscar:\n")
## Criterios del artículo a buscar:
cat("1. Reporte resultados de investigación experimental\n")
## 1. Reporte resultados de investigación experimental
cat("2. Aplique un diseño factorial\n")
## 2. Aplique un diseño factorial
cat("3. Al menos un factor tenga más de dos niveles\n")
## 3. Al menos un factor tenga más de dos niveles
cat("4. Esté publicado en revista científica o tecnológica\n")
## 4. Esté publicado en revista científica o tecnológica
cat("\nEstructura del reporte:\n")
##
## Estructura del reporte:
cat("A. Referencia completa:\n")
## A. Referencia completa:
cat(" - Autor(es)\n")
## - Autor(es)
cat(" - Año de publicación\n")
## - Año de publicación
cat(" - Título del trabajo\n")
## - Título del trabajo
cat(" - Nombre de la revista\n")
## - Nombre de la revista
cat(" - Volumen, número, páginas\n")
## - Volumen, número, páginas
cat(" - DOI (si está disponible)\n\n")
## - DOI (si está disponible)
cat("B. Síntesis del artículo:\n")
## B. Síntesis del artículo:
cat(" - Objetivo de la investigación\n")
## - Objetivo de la investigación
cat(" - Variables estudiadas\n")
## - Variables estudiadas
cat(" - Diseño experimental utilizado\n")
## - Diseño experimental utilizado
cat(" - Número de réplicas\n")
## - Número de réplicas
cat(" - Factores y sus niveles\n\n")
## - Factores y sus niveles
cat("C. Análisis estadísticos:\n")
## C. Análisis estadísticos:
cat(" - Métodos estadísticos empleados\n")
## - Métodos estadísticos empleados
cat(" - Software utilizado\n")
## - Software utilizado
cat(" - Pruebas de hipótesis realizadas\n")
## - Pruebas de hipótesis realizadas
cat(" - Nivel de significancia usado\n\n")
## - Nivel de significancia usado
cat("D. Principales conclusiones:\n")
## D. Principales conclusiones:
cat(" - Resultados estadísticamente significativos\n")
## - Resultados estadísticamente significativos
cat(" - Hallazgos principales\n")
## - Hallazgos principales
cat(" - Recomendaciones o aplicaciones\n")
## - Recomendaciones o aplicaciones
cat("\nEjemplo de formato de respuesta:\n")
##
## Ejemplo de formato de respuesta:
cat("=================================\n")
## =================================
cat("A. REFERENCIA COMPLETA:\n")
## A. REFERENCIA COMPLETA:
cat("García, A., Martínez, B., & Rodríguez, C. (2023).\n")
## García, A., Martínez, B., & Rodríguez, C. (2023).
cat("Efecto de la temperatura y pH en la producción de biomasa microbiana.\n")
## Efecto de la temperatura y pH en la producción de biomasa microbiana.
cat("Revista de Biotecnología Industrial, 45(3), 123-135.\n")
## Revista de Biotecnología Industrial, 45(3), 123-135.
cat("DOI: 10.1234/rbi.2023.045\n\n")
## DOI: 10.1234/rbi.2023.045
cat("B. SÍNTESIS:\n")
## B. SÍNTESIS:
cat("El estudio investigó el efecto combinado de temperatura (3 niveles: 25, 30, 35°C)\n")
## El estudio investigó el efecto combinado de temperatura (3 niveles: 25, 30, 35°C)
cat("y pH (4 niveles: 5.0, 6.0, 7.0, 8.0) en la producción de biomasa de Saccharomyces cerevisiae.\n")
## y pH (4 niveles: 5.0, 6.0, 7.0, 8.0) en la producción de biomasa de Saccharomyces cerevisiae.
cat("Se utilizó un diseño factorial completo 3×4 con 4 réplicas por tratamiento.\n\n")
## Se utilizó un diseño factorial completo 3×4 con 4 réplicas por tratamiento.
cat("C. ANÁLISIS ESTADÍSTICOS:\n")
## C. ANÁLISIS ESTADÍSTICOS:
cat("- ANOVA de dos vías con interacción\n")
## - ANOVA de dos vías con interacción
cat("- Prueba de Tukey para comparaciones múltiples\n")
## - Prueba de Tukey para comparaciones múltiples
cat("- Verificación de supuestos: Shapiro-Wilk y Levene\n")
## - Verificación de supuestos: Shapiro-Wilk y Levene
cat("- Software: R versión 4.2.1\n")
## - Software: R versión 4.2.1
cat("- Nivel de significancia: α = 0.05\n\n")
## - Nivel de significancia: α = 0.05
cat("D. PRINCIPALES CONCLUSIONES:\n")
## D. PRINCIPALES CONCLUSIONES:
cat("1. Ambos factores mostraron efectos significativos (p < 0.001)\n")
## 1. Ambos factores mostraron efectos significativos (p < 0.001)
cat("2. La interacción temperatura×pH fue significativa (p = 0.012)\n")
## 2. La interacción temperatura×pH fue significativa (p = 0.012)
cat("3. La combinación óptima fue 30°C y pH 6.0\n")
## 3. La combinación óptima fue 30°C y pH 6.0
cat("4. La producción máxima fue de 15.2 g/L de biomasa\n")
## 4. La producción máxima fue de 15.2 g/L de biomasa
cat("\n=================================\n")
##
## =================================
cat("Instrucciones para la búsqueda:\n")
## Instrucciones para la búsqueda:
cat("1. Ir a https://scholar.google.com\n")
## 1. Ir a https://scholar.google.com
cat("2. Usar términos como: 'factorial design', 'ANOVA', 'experimental design'\n")
## 2. Usar términos como: 'factorial design', 'ANOVA', 'experimental design'
cat("3. Filtrar por año reciente (últimos 5 años)\n")
## 3. Filtrar por año reciente (últimos 5 años)
cat("4. Verificar que el artículo esté en español o inglés\n")
## 4. Verificar que el artículo esté en español o inglés
cat("5. Descargar el PDF para revisar los métodos estadísticos\n")
## 5. Descargar el PDF para revisar los métodos estadísticos
cat("\nÁreas donde encontrar estos artículos:\n")
##
## Áreas donde encontrar estos artículos:
cat("- Ciencias agrícolas\n")
## - Ciencias agrícolas
cat("- Ingeniería química\n")
## - Ingeniería química
cat("- Biotecnología\n")
## - Biotecnología
cat("- Ciencias de los alimentos\n")
## - Ciencias de los alimentos
cat("- Farmacología\n")
## - Farmacología
cat("- Ciencias ambientales\n")
## - Ciencias ambientales
cat("\nPuntos a evaluar en el artículo:\n")
##
## Puntos a evaluar en el artículo:
cat("✓ ¿Describe claramente el diseño experimental?\n")
## ✓ ¿Describe claramente el diseño experimental?
cat("✓ ¿Especifica los factores y sus niveles?\n")
## ✓ ¿Especifica los factores y sus niveles?
cat("✓ ¿Reporta el número de réplicas?\n")
## ✓ ¿Reporta el número de réplicas?
cat("✓ ¿Incluye análisis estadístico apropiado?\n")
## ✓ ¿Incluye análisis estadístico apropiado?
cat("✓ ¿Presenta tablas ANOVA?\n")
## ✓ ¿Presenta tablas ANOVA?
cat("✓ ¿Discute la significancia práctica de los resultados?\n")
## ✓ ¿Discute la significancia práctica de los resultados?
# Sugerencia de términos de búsqueda
cat("\nTérminos sugeridos para la búsqueda:\n")
##
## Términos sugeridos para la búsqueda:
cat("- 'factorial experiment design'\n")
## - 'factorial experiment design'
cat("- 'ANOVA multi-factor'\n")
## - 'ANOVA multi-factor'
cat("- 'response surface methodology'\n")
## - 'response surface methodology'
cat("- 'design of experiments'\n")
## - 'design of experiments'
cat("- 'factorial ANOVA application'\n")
## - 'factorial ANOVA application'
cat("- 'multi-level factorial design'\n")
## - 'multi-level factorial design'
# Ejemplo de código para analizar datos similares (opcional)
cat("\nCódigo R para analizar un diseño factorial (ejemplo genérico):\n")
##
## Código R para analizar un diseño factorial (ejemplo genérico):
cat("```r\n")
## ```r
cat("# Ejemplo de análisis factorial\n")
## # Ejemplo de análisis factorial
cat("# Datos ficticios\n")
## # Datos ficticios
cat("factor1 <- factor(rep(c('A1', 'A2', 'A3'), each = 12))\n")
## factor1 <- factor(rep(c('A1', 'A2', 'A3'), each = 12))
cat("factor2 <- factor(rep(rep(c('B1', 'B2', 'B3', 'B4'), each = 3), 3))\n")
## factor2 <- factor(rep(rep(c('B1', 'B2', 'B3', 'B4'), each = 3), 3))
cat("respuesta <- rnorm(36, mean = 10, sd = 2)\n")
## respuesta <- rnorm(36, mean = 10, sd = 2)
cat("datos <- data.frame(factor1, factor2, respuesta)\n\n")
## datos <- data.frame(factor1, factor2, respuesta)
cat("# ANOVA factorial\n")
## # ANOVA factorial
cat("modelo <- aov(respuesta ~ factor1 * factor2, data = datos)\n")
## modelo <- aov(respuesta ~ factor1 * factor2, data = datos)
cat("anova_result <- anova(modelo)\n")
## anova_result <- anova(modelo)
cat("print(anova_result)\n\n")
## print(anova_result)
cat("# Comparaciones múltiples si hay efectos significativos\n")
## # Comparaciones múltiples si hay efectos significativos
cat("if(anova_result['factor1', 'Pr(>F)'] < 0.05) {\n")
## if(anova_result['factor1', 'Pr(>F)'] < 0.05) {
cat(" tukey_factor1 <- TukeyHSD(modelo, 'factor1')\n")
## tukey_factor1 <- TukeyHSD(modelo, 'factor1')
cat(" print(tukey_factor1$factor1)\n")
## print(tukey_factor1$factor1)
cat("}\n")
## }
cat("```\n")
## ```
cat("\nRecomendaciones finales:\n")
##
## Recomendaciones finales:
cat("1. Elegir un artículo de una revista indexada\n")
## 1. Elegir un artículo de una revista indexada
cat("2. Preferir artículos de acceso abierto\n")
## 2. Preferir artículos de acceso abierto
cat("3. Verificar que los métodos estadísticos sean adecuados\n")
## 3. Verificar que los métodos estadísticos sean adecuados
cat("4. Buscar artículos que también discutan las limitaciones del estudio\n")
## 4. Buscar artículos que también discutan las limitaciones del estudio