# Cargar librerías necesarias
library(ggplot2)
library(dplyr)
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
# Crear el dataframe
datos <- data.frame(
A = rep(c(-1, 1), each = 4),
B = rep(c(-1, 1), times = 4),
C = c(-1, -1, 1, 1, -1, -1, 1, 1),
Rep1 = c(50, 44, 46, 42, 49, 48, 47, 56),
Rep2 = c(54, 42, 48, 43, 46, 45, 48, 54)
)
# Promedio de réplicas
datos$y <- rowMeans(datos[, c("Rep1", "Rep2")])
modelo <- lm(y ~ A + B + C, data = datos)
step(modelo)
## Start: AIC=28.01
## y ~ A + B + C
##
## Df Sum of Sq RSS AIC
## - C 1 1.125 98.750 26.105
## - B 1 6.125 103.750 26.500
## - A 1 18.000 115.625 27.367
## <none> 97.625 28.014
##
## Step: AIC=26.11
## y ~ A + B
##
## Df Sum of Sq RSS AIC
## - B 1 6.125 104.88 24.587
## - A 1 18.000 116.75 25.445
## <none> 98.75 26.105
##
## Step: AIC=24.59
## y ~ A
##
## Df Sum of Sq RSS AIC
## - A 1 18 122.88 23.854
## <none> 104.88 24.587
##
## Step: AIC=23.85
## y ~ 1
##
## Call:
## lm(formula = y ~ 1, data = datos)
##
## Coefficients:
## (Intercept)
## 47.62
anova(modelo)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 18.000 18.000 0.7375 0.4389
## B 1 6.125 6.125 0.2510 0.6427
## C 1 1.125 1.125 0.0461 0.8405
## Residuals 4 97.625 24.406
summary(modelo)
##
## Call:
## lm(formula = y ~ A + B + C, data = datos)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## 5.375 -1.875 -0.375 -3.125 -2.125 -1.375 -2.875 6.375
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47.625 1.747 27.267 1.08e-05 ***
## A 1.500 1.747 0.859 0.439
## B -0.875 1.747 -0.501 0.643
## C 0.375 1.747 0.215 0.841
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.94 on 4 degrees of freedom
## Multiple R-squared: 0.2055, Adjusted R-squared: -0.3904
## F-statistic: 0.3449 on 3 and 4 DF, p-value: 0.7958
# Reordenar para gráficos
datos_largo <- datos %>%
pivot_longer(cols = c("Rep1", "Rep2"), names_to = "Rep", values_to = "Valor")
# Efectos principales
ggplot(datos_largo, aes(x = factor(A), y = Valor)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line", aes(group = 1)) +
labs(x = "Factor A", y = "Número de órdenes", title = "Efecto principal de A")
ggplot(datos_largo, aes(x = factor(B), y = Valor)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line", aes(group = 1)) +
labs(x = "Factor B", y = "Número de órdenes", title = "Efecto principal de B")
ggplot(datos_largo, aes(x = factor(C), y = Valor)) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line", aes(group = 1)) +
labs(x = "Factor C", y = "Número de órdenes", title = "Efecto principal de C")
ggplot(datos_largo, aes(x = factor(A), y = Valor, color = factor(B), group = factor(B))) +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line") +
labs(x = "Factor A", y = "Número de órdenes", color = "B", title = "Interacción A:B")
library(ggplot2)
library(dplyr)
library(tidyr)
datos <- data.frame(
Run = 1:8,
A = c("-", "+", "-", "+", "-", "+", "-", "+"),
B = c("-", "-", "+", "+", "-", "-", "+", "+"),
C = c("-", "-", "-", "-", "+", "+", "+", "+"),
Rep1 = c(50, 44, 46, 42, 49, 48, 47, 56),
Rep2 = c(54, 42, 48, 43, 46, 45, 48, 54)
)
datos_largo <- datos %>%
mutate(
Clase = ifelse(A == "-", "3rd", "1st"),
Folleto = ifelse(B == "-", "BW", "Color"),
Precio = ifelse(C == "-", "$19.95", "$24.95")
) %>%
select(Run, Clase, Folleto, Precio, Rep1, Rep2) %>%
pivot_longer(cols = c(Rep1, Rep2), names_to = "Replica", values_to = "Ordenes")
modelo <- aov(Ordenes ~ Clase * Folleto * Precio, data = datos_largo)
summary(modelo)
## Df Sum Sq Mean Sq F value Pr(>F)
## Clase 1 12.25 12.25 4.083 0.077971 .
## Folleto 1 2.25 2.25 0.750 0.411694
## Precio 1 36.00 36.00 12.000 0.008516 **
## Clase:Folleto 1 42.25 42.25 14.083 0.005602 **
## Clase:Precio 1 100.00 100.00 33.333 0.000418 ***
## Folleto:Precio 1 49.00 49.00 16.333 0.003728 **
## Clase:Folleto:Precio 1 4.00 4.00 1.333 0.281537
## Residuals 8 24.00 3.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(datos_largo, aes(x = Clase, y = Ordenes, color = Folleto, shape = Precio, group = interaction(Folleto, Precio))) +
stat_summary(fun = mean, geom = "point", size = 3) +
stat_summary(fun = mean, geom = "line") +
theme_minimal() +
labs(title = "Interacciones entre factores", y = "Número de órdenes", x = "Tipo de correo") +
theme(legend.position = "bottom")
# Crear los factores codificados (A y B)
A <- rep(c(-1, -1, 1, 1), each = 6) # Medio de cultivo
B <- rep(c(-1, 1, -1, 1), each = 6) # Tiempo
# Respuesta: crecimiento bacteriano (24 observaciones = 4 condiciones * 6 réplicas)
y <- c(
23, 26, 22, 29, 38, 36, # A = -1, B = -1
28, 29, 25, 27, 38, 30, # A = -1, B = +1
20, 25, 21, 24, 38, 35, # A = +1, B = -1
26, 33, 25, 20, 33, 34 # A = +1, B = +1
)
# Crear data frame
datos <- data.frame(A = factor(A), B = factor(B), y = y)
modelo <- aov(y ~ A * B, data = datos)
summary(modelo)
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 12.0 12.04 0.313 0.582
## B 1 5.0 5.04 0.131 0.721
## A:B 1 1.0 1.04 0.027 0.871
## Residuals 20 769.8 38.49
# Ajustar el modelo
modelo <- lm(y ~ A * B, data = datos)
# Obtener los residuos
res <- resid(modelo)
# 1. Gráfico de residuos vs valores ajustados
plot(fitted(modelo), res,
main = "Residuos vs Valores Ajustados",
xlab = "Valores Ajustados",
ylab = "Residuos")
abline(h = 0, col = "red")
# 2. Gráfico Q-Q para evaluar normalidad
qqnorm(res)
qqline(res, col = "blue")
# 3. Prueba de normalidad (Shapiro-Wilk)
shapiro.test(res)
##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.9301, p-value = 0.09801
# Supongamos que ya tienes el diseño codificado con las columnas A, B, C, D, Bloque y respuesta
# Creamos las variables cuadráticas manualmente para simular curvatura (aunque no hay punto central real)
# 4. Prueba de homocedasticidad (Breusch-Pagan)
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.4.3
## Cargando paquete requerido: zoo
##
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
bptest(modelo)
##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 3.2541, df = 3, p-value = 0.3541
# ------------------------------
# Cargar paquetes necesarios
# ------------------------------
library(FrF2)
## Warning: package 'FrF2' was built under R version 4.4.3
## Cargando paquete requerido: DoE.base
## Warning: package 'DoE.base' was built under R version 4.4.3
## Cargando paquete requerido: grid
## Cargando paquete requerido: conf.design
## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design
##
## Adjuntando el paquete: 'DoE.base'
## The following objects are masked from 'package:stats':
##
## aov, lm
## The following object is masked from 'package:graphics':
##
## plot.design
## The following object is masked from 'package:base':
##
## lengths
library(ggplot2)
library(car)
## Warning: package 'car' was built under R version 4.4.3
## Cargando paquete requerido: carData
##
## Adjuntando el paquete: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
# ------------------------------
# Crear data frame con factores codificados y respuesta
# ------------------------------
factores <- expand.grid(
A = c(-1, 1),
B = c(-1, 1),
C = c(-1, 1),
D = c(-1, 1),
E = c(-1, 1)
)
# Respuesta (rendimiento) del enunciado
respuesta <- c(7,8,8,6,9,10,12,10,34,32,35,30,55,50,52,53,
16,18,15,15,20,21,22,20,40,44,45,41,60,61,65,63)
datos <- cbind(factores, Rendimiento = respuesta)
datos <- as.data.frame(datos)
# Supongamos que ya tienes el diseño codificado con las columnas A, B, C, D, Bloque y respuesta
# Creamos las variables cuadráticas manualmente para simular curvatura (aunque no hay punto central real)
datos$A2 <- as.numeric(as.character(datos$A))^2
datos$B2 <- as.numeric(as.character(datos$B))^2
datos$C2 <- as.numeric(as.character(datos$C))^2
datos$D2 <- as.numeric(as.character(datos$D))^2
# ------------------------------
# (a) Gráfico de probabilidad normal de efectos
# ------------------------------
modelo_full <- lm(Rendimiento ~ (A + B + C + D + E)^5, data = datos)
efectos <- coef(modelo_full)[-1] # quitar el intercepto
qqnorm(efectos, main = "QQ Plot de los efectos")
#qqline(efectos)
df <- data.frame(efectos = efectos)
ggplot(df, aes(x = efectos)) +
geom_histogram(aes(y = ..density..), fill = "blue", color = "black", bins = 15, alpha = 0.6) +
geom_density(color = "red", size = 1) +
labs(title = "Histograma de efectos con densidad ajustada",
x = "Efectos", y = "Densidad") +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# ------------------------------
# (b) ANOVA para confirmar significancia
# ------------------------------
anova_modelo <- aov(Rendimiento ~ (A + B + C + D + E)^2, data = datos)
summary(anova_modelo)
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 5 5 2.126 0.164
## B 1 2 2 0.616 0.444
## C 1 1116 1116 449.321 3.90e-13 ***
## D 1 9214 9214 3708.792 < 2e-16 ***
## E 1 751 751 302.201 8.21e-12 ***
## A:B 1 11 11 4.541 0.049 *
## A:C 1 0 0 0.013 0.912
## A:D 1 4 4 1.522 0.235
## A:E 1 5 5 2.126 0.164
## B:C 1 7 7 2.830 0.112
## B:D 1 3 3 1.019 0.328
## B:E 1 1 1 0.314 0.583
## C:D 1 504 504 202.881 1.66e-10 ***
## C:E 1 2 2 0.616 0.444
## D:E 1 0 0 0.013 0.912
## Residuals 16 40 2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ------------------------------
# (c) Modelo con efectos significativos
# Supongamos que A, B y C son significativos
# Asegurarse de que las variables sean factores
datos
## A B C D E Rendimiento A2 B2 C2 D2
## 1 -1 -1 -1 -1 -1 7 1 1 1 1
## 2 1 -1 -1 -1 -1 8 1 1 1 1
## 3 -1 1 -1 -1 -1 8 1 1 1 1
## 4 1 1 -1 -1 -1 6 1 1 1 1
## 5 -1 -1 1 -1 -1 9 1 1 1 1
## 6 1 -1 1 -1 -1 10 1 1 1 1
## 7 -1 1 1 -1 -1 12 1 1 1 1
## 8 1 1 1 -1 -1 10 1 1 1 1
## 9 -1 -1 -1 1 -1 34 1 1 1 1
## 10 1 -1 -1 1 -1 32 1 1 1 1
## 11 -1 1 -1 1 -1 35 1 1 1 1
## 12 1 1 -1 1 -1 30 1 1 1 1
## 13 -1 -1 1 1 -1 55 1 1 1 1
## 14 1 -1 1 1 -1 50 1 1 1 1
## 15 -1 1 1 1 -1 52 1 1 1 1
## 16 1 1 1 1 -1 53 1 1 1 1
## 17 -1 -1 -1 -1 1 16 1 1 1 1
## 18 1 -1 -1 -1 1 18 1 1 1 1
## 19 -1 1 -1 -1 1 15 1 1 1 1
## 20 1 1 -1 -1 1 15 1 1 1 1
## 21 -1 -1 1 -1 1 20 1 1 1 1
## 22 1 -1 1 -1 1 21 1 1 1 1
## 23 -1 1 1 -1 1 22 1 1 1 1
## 24 1 1 1 -1 1 20 1 1 1 1
## 25 -1 -1 -1 1 1 40 1 1 1 1
## 26 1 -1 -1 1 1 44 1 1 1 1
## 27 -1 1 -1 1 1 45 1 1 1 1
## 28 1 1 -1 1 1 41 1 1 1 1
## 29 -1 -1 1 1 1 60 1 1 1 1
## 30 1 -1 1 1 1 61 1 1 1 1
## 31 -1 1 1 1 1 65 1 1 1 1
## 32 1 1 1 1 1 63 1 1 1 1
datos$A <- factor(datos$A)
datos$B <- factor(datos$B)
datos$C <- factor(datos$C)
datos$D <- factor(datos$D)
datos$E <- factor(datos$E)
# Modelo con efectos principales y solo las interacciones significativas
modelo_sig <- aov(respuesta ~ C + D + E + A:B + C:D, data = datos)
# Ver resumen si deseas
summary(modelo_sig)
## Df Sum Sq Mean Sq F value Pr(>F)
## C 1 1116 1116 441.000 < 2e-16 ***
## D 1 9214 9214 3640.111 < 2e-16 ***
## E 1 751 751 296.605 5.17e-15 ***
## A:B 3 18 6 2.383 0.0944 .
## C:D 1 504 504 199.123 4.07e-13 ***
## Residuals 24 61 3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ------------------------------
# (d) Gráfico normal de residuos
# ------------------------------
# Extraer los residuos
residuos <- residuals(modelo_sig)
# Ver los residuos
print(residuos)
## 1 2 3 4 5
## 6.250000e-01 1.250000e+00 1.115774e-14 -7.993606e-15 -1.250000e+00
## 6 7 8 9 10
## -6.250000e-01 1.250000e-01 1.250000e-01 1.625000e+00 -7.500000e-01
## 11 12 13 14 15
## 1.000000e+00 -2.000000e+00 2.875000e+00 -2.500000e+00 -1.750000e+00
## 16 17 18 19 20
## 1.250000e+00 -6.250000e-02 1.562500e+00 -2.687500e+00 -6.875000e-01
## 21 22 23 24 25
## 6.250000e-02 6.875000e-01 4.375000e-01 4.375000e-01 -2.062500e+00
## 26 27 28 29 30
## 1.562500e+00 1.312500e+00 -6.875000e-01 -1.812500e+00 -1.187500e+00
## 31 32
## 1.562500e+00 1.562500e+00
# Opción base R
qqnorm(residuos, main = "QQ Plot de los residuos")
#qqline(residuos, col = "red", lwd = 2)
# ------------------------------
# (e) Gráficos de residuos
# ------------------------------
# Supongamos que ya ajustaste el modelo significativo:
modelo_sig <- aov(respuesta ~ C + D + E + A:B + C:D, data = datos)
# Extraer residuos y valores ajustados
residuos <- residuals(modelo_sig)
predichos <- fitted(modelo_sig)
# Graficar residuos vs valores predichos
plot(predichos, residuos,
main = "Residuos vs Valores Predichos",
xlab = "Valores Predichos",
ylab = "Residuos",
pch = 19, col = "blue")
abline(h = 0, col = "red", lty = 2)
# Graficar residuos vs cada factor
par(mfrow = c(2, 3)) # Mostrar varios gráficos juntos
for (factor in c("A", "B", "C", "D", "E")) {
plot(datos[[factor]], residuos,
main = paste("Residuos vs", factor),
xlab = factor,
ylab = "Residuos",
pch = 19, col = "darkgreen")
abline(h = 0, col = "red", lty = 2)
}
par(mfrow = c(1, 1)) # Restaurar layout
# ------------------------------
# (f) Interpretación de interacciones (ver summary)
# ------------------------------
# Ya incluido en summary(modelo_reducido)
# ------------------------------
# (g) Recomendaciones operativas
# ------------------------------
# Se pueden usar predicciones para identificar configuraciones óptimas
#datos$Predicho <- predict(modelo_reducido)
#datos[order(-datos$Predicho), ][1:5, ] # top 5 combinaciones
# ------------------------------
# (h) Proyección a diseño 2^k y resumen por combinación significativa
# Supongamos que A, B, y C son importantes
datos$ABC <- interaction(datos$A, datos$B, datos$C)
agregado <- aggregate(Rendimiento ~ ABC, data = datos, FUN = function(x) c(prom = mean(x), rango = max(x) - min(x)))
print(agregado)
## ABC Rendimiento.prom Rendimiento.rango
## 1 -1.-1.-1 24.25 33.00
## 2 1.-1.-1 25.50 36.00
## 3 -1.1.-1 25.75 37.00
## 4 1.1.-1 23.00 35.00
## 5 -1.-1.1 36.00 51.00
## 6 1.-1.1 35.50 51.00
## 7 -1.1.1 37.75 53.00
## 8 1.1.1 36.50 53.00
# -------------------------
# Preparar los datos
# -------------------------
tratamientos <- c("(1)", "a", "b", "ab", "c", "ac", "bc", "abc",
"d", "ad", "bd", "abd", "cd", "acd", "bcd", "abcd")
# Codificación de los factores en niveles -1 / +1
factores <- expand.grid(
A = c(-1,1),
B = c(-1,1),
C = c(-1,1),
D = c(-1,1)
)
# Agregar nombre de tratamientos en el mismo orden
factores$Tratamiento <- tratamientos
# Agregar la respuesta de la primera réplica (como pide el ejercicio)
respuesta_I <- c(90, 74, 81, 83, 77, 81, 88, 73,
98, 72, 87, 85, 99, 79, 87, 80)
# -------------------------
# Crear el data frame completo
# -------------------------
datos <- cbind(factores, Y = respuesta_I)
# -------------------------
# Confundir ABCD con bloques
# -------------------------
# ABCD = producto de A*B*C*D: +1 o -1 → define bloque
datos$Bloque <- ifelse(datos$A * datos$B * datos$C * datos$D == 1, "Bloque1", "Bloque2")
datos$Bloque <- as.factor(datos$Bloque)
# Cargar paquetes necesarios
library(dplyr)
# Crear data frame con réplicas
tratamientos <- c("(1)", "a", "b", "ab", "c", "ac", "bc", "abc",
"d", "ad", "bd", "abd", "cd", "acd", "bcd", "abcd")
rep1 <- c(90, 74, 81, 83, 77, 81, 88, 73, 98, 72, 87, 85, 99, 79, 87, 80)
rep2 <- c(93, 78, 85, 80, 78, 80, 82, 70, 95, 76, 83, 86, 90, 75, 84, 80)
# Calcular promedio de ambas réplicas
y <- (rep1 + rep2) / 2
# Crear diseño codificado (4 factores)
diseño <- expand.grid(A = c(-1, 1), B = c(-1, 1),
C = c(-1, 1), D = c(-1, 1))
# Agregar la respuesta
datos <- cbind(diseño, y = y)
# Crear variable de curvatura: número de factores en +1 (nivel alto)
datos$curvatura <- apply(datos[, 1:4], 1, function(x) sum(x == 1))
# Modelo factorial completo (sin curvatura)
modelo_lineal <- lm(y ~ A * B * C * D, data = datos)
# Modelo incluyendo curvatura como covariable
modelo_con_curva <- lm(y ~ A * B * C * D + curvatura, data = datos)
# Ver resúmenes
summary(modelo_con_curva)
##
## Call:
## lm.default(formula = y ~ A * B * C * D + curvatura, data = datos)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82.78125 NaN NaN NaN
## A -4.53125 NaN NaN NaN
## B -0.65625 NaN NaN NaN
## C -1.34375 NaN NaN NaN
## D 1.96875 NaN NaN NaN
## curvatura NA NA NA NA
## A:B 2.03125 NaN NaN NaN
## A:C 0.34375 NaN NaN NaN
## B:C -0.28125 NaN NaN NaN
## A:D -1.09375 NaN NaN NaN
## B:D -0.09375 NaN NaN NaN
## C:D 0.84375 NaN NaN NaN
## A:B:C -2.59375 NaN NaN NaN
## A:B:D 2.34375 NaN NaN NaN
## A:C:D -0.46875 NaN NaN NaN
## B:C:D -0.46875 NaN NaN NaN
## A:B:C:D 1.21875 NaN NaN NaN
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 15 and 0 DF, p-value: NA
# Comparar modelos
anova(modelo_lineal, modelo_con_curva)
## Analysis of Variance Table
##
## Model 1: y ~ A * B * C * D
## Model 2: y ~ A * B * C * D + curvatura
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 0 0
## 2 0 0 0 0
# Diseño factorial 2^5 sin replicación con tres respuestas simuladas
# Cargar paquetes necesarios
library(FrF2) # Para generar el diseño
library(ggplot2) # Para visualización
library(car) # Para ANOVA
# Factores y niveles
factors <- c("SiH2Cl2", "GeH4", "HCl", "B2H6", "Temp")
levels_list <- list(
SiH2Cl2 = c(8.0, 12.0),
GeH4 = c(7.2, 10.8),
HCl = c(3.2, 4.8),
B2H6 = c(4.4, 6.6),
Temp = c(740, 760)
)
# Crear diseño factorial 2^5
set.seed(123)
design <- FrF2(nruns = 32, nfactors = 5, factor.names = factors)
## creating full factorial with 32 runs ...
# Convertir niveles codificados (-1, 1) a niveles reales
for (f in factors) {
low <- levels_list[[f]][1]
high <- levels_list[[f]][2]
design[[f]] <- ifelse(design[[f]] == -1, low, high)
}
# Simular respuestas (sustituye por tus datos reales si los tienes)
design$y1 <- 50 + 5*design$SiH2Cl2 - 2*design$Temp + rnorm(32, 0, 2)
design$y2 <- 100 - 3*design$GeH4 + 1.5*design$B2H6 + rnorm(32, 0, 2)
design$y3 <- 30 + 2*design$HCl + 3*design$Temp + rnorm(32, 0, 2)
# Ajustar modelos lineales para cada respuesta
mod1 <- lm(y1 ~ . -y2 -y3, data = design)
mod2 <- lm(y2 ~ . -y1 -y3, data = design)
mod3 <- lm(y3 ~ . -y1 -y2, data = design)
# ANOVA
summary(aov(mod1))
## Df Sum Sq Mean Sq F value Pr(>F)
## SiH2Cl2 1 3569 3569 966.630 <2e-16 ***
## GeH4 1 0 0 0.042 0.839
## HCl 1 3 3 0.683 0.416
## B2H6 1 12 12 3.122 0.089 .
## Temp 1 12502 12502 3385.712 <2e-16 ***
## Residuals 26 96 4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(mod2))
## Df Sum Sq Mean Sq F value Pr(>F)
## SiH2Cl2 1 8.7 8.7 2.968 0.0968 .
## GeH4 1 852.4 852.4 291.211 1.22e-15 ***
## HCl 1 13.6 13.6 4.637 0.0407 *
## B2H6 1 100.2 100.2 34.233 3.62e-06 ***
## Temp 1 0.2 0.2 0.072 0.7901
## Residuals 26 76.1 2.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(mod3))
## Df Sum Sq Mean Sq F value Pr(>F)
## SiH2Cl2 1 1 1 0.213 0.649
## GeH4 1 0 0 0.010 0.922
## HCl 1 99 99 35.139 2.97e-06 ***
## B2H6 1 0 0 0.045 0.833
## Temp 1 28864 28864 10230.176 < 2e-16 ***
## Residuals 26 73 3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Opcional: gráfico de efectos para cada respuesta
library(DoE.base)
eff1 <- DanielPlot(mod1, main = "Daniel Plot - y1")
## simulated critical values not available for all requests, used conservative ones
eff2 <- DanielPlot(mod2, main = "Daniel Plot - y2")
## simulated critical values not available for all requests, used conservative ones
eff3 <- DanielPlot(mod3, main = "Daniel Plot - y3")
## simulated critical values not available for all requests, used conservative ones
library(FrF2)
# Crear diseño 2^4 con codificación estándar y orden aleatorio (como en tabla)
design <- FrF2(nruns = 16, nfactors = 4, factor.names = c("A", "B", "C", "D"), randomize = FALSE)
## creating full factorial with 16 runs ...
#actorial completo 2^4 sin aleatorización
design <- FrF2(nruns = 16, nfactors = 4, factor.names = c("A", "B", "C", "D"), randomize = FALSE)
## creating full factorial with 16 runs ...
# Vector de rendimientos (respuesta)
yield <- c(12, 18, 13, 16, 17, 15, 20, 15, 10, 25, 13, 24, 19, 21, 17, 23)
# Agregar la respuesta al diseño
design$Yield <- yield
# Ajustar modelo con interacciones completas hasta el orden 4
modelo <- lm(Yield ~ A*B*C*D, data = design)
# Graficar efectos estimados en gráfico de Daniel (normal probability plot)
par(mfrow=c(1,1))
DanielPlot(modelo, alpha = 0.05, main = "Gráfico de probabilidad normal de los efectos")
library(FrF2)
library(daewr)
## Warning: package 'daewr' was built under R version 4.4.3
design <- FrF2(nruns = 16, nfactors = 4, factor.names = c("A", "B", "C", "D"), randomize = FALSE)
## creating full factorial with 16 runs ...
design$Yield <- c(12, 18, 13, 16, 17, 15, 20, 15, 10, 25, 13, 24, 19, 21, 17, 23)
modelo <- lm(Yield ~ A*B*C*D, data = design)
DanielPlot(modelo, alpha = 0.05)
modelo_reducido <- lm(Yield ~ A + B + D + B:D, data = design)
summary(modelo_reducido)
##
## Call:
## lm.default(formula = Yield ~ A + B + D + B:D, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.500 -2.375 0.000 2.500 6.250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.738e+01 9.756e-01 17.810 1.85e-09 ***
## A1 2.250e+00 9.756e-01 2.306 0.0416 *
## B1 2.500e-01 9.756e-01 0.256 0.8025
## D1 1.625e+00 9.756e-01 1.666 0.1240
## B1:D1 4.441e-16 9.756e-01 0.000 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.902 on 11 degrees of freedom
## Multiple R-squared: 0.4259, Adjusted R-squared: 0.2171
## F-statistic: 2.04 on 4 and 11 DF, p-value: 0.158
anova(modelo_reducido)
## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 81.00 81.000 5.3194 0.04156 *
## B 1 1.00 1.000 0.0657 0.80248
## D 1 42.25 42.250 2.7746 0.12396
## B:D 1 0.00 0.000 0.0000 1.00000
## Residuals 11 167.50 15.227
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1. Gráfico de residuos vs valores ajustados
plot(modelo_reducido$fitted.values, resid(modelo_reducido),
main = "Residuos vs Valores Ajustados",
xlab = "Valores ajustados", ylab = "Residuos")
#abline(h = 0, col = "red", lty = 2)
# 2. Normal Q-Q plot (residuos)
qqnorm(resid(modelo_reducido), main = "Q-Q Plot de los residuos")
#qqline(resid(modelo_reducido), col = "red")
# 3. Histograma de residuos
hist(resid(modelo_reducido), main = "Histograma de residuos", xlab = "Residuos")
# 4. Prueba de normalidad de Shapiro-Wilk
shapiro.test(resid(modelo_reducido))
##
## Shapiro-Wilk normality test
##
## data: resid(modelo_reducido)
## W = 0.99104, p-value = 0.9996
# 5. Prueba de homocedasticidad (Breusch-Pagan)
library(lmtest)
bptest(modelo_reducido)
##
## studentized Breusch-Pagan test
##
## data: modelo_reducido
## BP = 4.1382, df = 4, p-value = 0.3876
#9.11
# Cargar paquetes si es necesario
# install.packages("dplyr")
library(dplyr)
# Definir los niveles de los factores
botella <- rep(c("Plástico", "28mm", "38mm"), times = 3*3*2)
estante <- rep(c("Permanente", "Pasillo", "Enfriador"), each = 3, times = 6)
trabajador <- rep(c("1", "2", "3"), each = 9, times = 2)
replica <- rep(c("I", "II"), each = 27)
# Datos de tiempo (en minutos)
tiempo <- c(
# Réplica I
3.45, 4.14, 5.80, 4.07, 4.38, 5.48, 4.20, 4.26, 5.67,
4.80, 5.22, 6.21, 4.52, 5.15, 6.25, 4.96, 5.17, 6.03,
4.08, 3.94, 5.14, 4.30, 4.53, 4.99, 4.17, 4.86, 4.85,
# Réplica II
3.36, 4.19, 5.23, 3.52, 4.26, 4.85, 3.68, 4.37, 5.58,
4.40, 4.70, 5.88, 4.44, 4.65, 6.20, 4.39, 4.75, 6.38,
3.65, 4.08, 4.49, 4.04, 4.08, 4.59, 3.88, 4.48, 4.90
)
# Crear el data frame
datos <- data.frame(
Replica = as.factor(replica),
Trabajador = as.factor(trabajador),
Botella = as.factor(botella),
Estante = as.factor(estante),
Tiempo = tiempo
)
# Ajustar modelo factorial con bloques (réplica)
modelo <- aov(Tiempo ~ Botella * Estante * Trabajador + Error(Replica), data = datos)
# Resumen del modelo
summary(modelo)
##
## Error: Replica
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 1.07 1.07
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Botella 2 17.751 8.876 241.549 < 2e-16 ***
## Estante 2 0.410 0.205 5.585 0.0096 **
## Trabajador 2 7.663 3.832 104.274 3.82e-13 ***
## Botella:Estante 4 0.116 0.029 0.789 0.5429
## Botella:Trabajador 4 1.679 0.420 11.423 1.74e-05 ***
## Estante:Trabajador 4 0.109 0.027 0.743 0.5712
## Botella:Estante:Trabajador 8 0.553 0.069 1.881 0.1068
## Residuals 26 0.955 0.037
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1