Metodología de Superficie de Respuesta (Análisis estadístico)
Jorge Ramos-Frutos jorgearmando.rf@jiquilpan.tecnm.mx
2025-08-26
Metodología de Superficie de Respuesta
Paquetes utilizados en este tema
En Metodología de Superficie de Respuesta (MSR) se utilizan los paquetes siguientes:
library(rsm)
library(ggplot2)
library(car)## Loading required package: carDatalibrary(gridExtra)
library(scatterplot3d)
library(plot3D)## Warning in fun(libname, pkgname): couldn't connect to display ":0"library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:gridExtra':
##
## combine## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutDiseños experimentales MSR
Los diseños experimentales más comúnes para MSR con los que se obtienen modelos de segundo orden son los siguientes:
- Diseño Central Compuesto.
- Diseño Centrado en Caras.
- Diseño Box-Behnken
A continuación se muestra el código para construir estos diseños de experimentos:
Diseño Central Compuesto
Este diseño combina puntos factoriales (\(2^k\)), puntos axiales (\(2k\)) y puntos centrales (\(n0\)). La cantidad total de corridas es:
\[ n_{corridas} = 2^{k} + 2k + n_{0} \] y la distancia de los puntos axiales (\(\alpha\)) se obtiene de la siguiente manera:
\[ \alpha = (2^{k})^{1/4} \] para que el diseño sea rotable. Las ventajas de este modelo son:
- Eficiente para estimar modelos cuadráticos.
- Permite la rotabilidad del diseño.
- Buena relación entre número de puntos y información obtenida.
A continuación se presenta un diseño central compuesto para dos factores, cuatro puntos al centro y una réplica:
# Diseño Central Compuesto para 2 factores
design_ccc <- ccd(2,
n0 = 2, # Número de puntos centrales
alpha = "rotatable", # Diseño rotable (alpha = √2)
randomize = FALSE, # Aleatorizar
oneblock = TRUE, # Un bloque
coding = list(x1 ~ (Factor_A - 0)/1,
x2 ~ (Factor_B - 0)/1))
# Ver el diseño generado
print(design_ccc)## run.order std.order Factor_A Factor_B
## 1 1 1 -1.000000 -1.000000
## 2 2 2 1.000000 -1.000000
## 3 3 3 -1.000000 1.000000
## 4 4 4 1.000000 1.000000
## 5 5 5 0.000000 0.000000
## 6 6 6 0.000000 0.000000
## 7 1 1 -1.414214 0.000000
## 8 2 2 1.414214 0.000000
## 9 3 3 0.000000 -1.414214
## 10 4 4 0.000000 1.414214
## 11 5 5 0.000000 0.000000
## 12 6 6 0.000000 0.000000
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Factor_A - 0)/1
## x2 ~ (Factor_B - 0)/1# Mostrar en formato de data frame
design_df <- as.data.frame(design_ccc)
print(design_df)## run.order std.order x1 x2
## 1 1 1 -1.000000 -1.000000
## 2 2 2 1.000000 -1.000000
## 3 3 3 -1.000000 1.000000
## 4 4 4 1.000000 1.000000
## 5 5 5 0.000000 0.000000
## 6 6 6 0.000000 0.000000
## 7 1 1 -1.414214 0.000000
## 8 2 2 1.414214 0.000000
## 9 3 3 0.000000 -1.414214
## 10 4 4 0.000000 1.414214
## 11 5 5 0.000000 0.000000
## 12 6 6 0.000000 0.000000# Puntos factoriales (esquinas)
puntos_factoriales <- expand.grid(
Factor_A = c(-1, 1),
Factor_B = c(-1, 1)
)
# Puntos axiales (alpha = √2 ≈ 1.414 para diseño rotable)
alpha <- (4)**(1/4)
puntos_axiales <- rbind(
data.frame(Factor_A = c(-alpha, alpha, 0, 0),
Factor_B = c(0, 0, -alpha, alpha))
)
# 4 Puntos centrales
puntos_centrales <- data.frame(
Factor_A = rep(0, 4),
Factor_B = rep(0, 4)
)
# Combinar todos los puntos
simple_ccc <- rbind(
cbind(puntos_factoriales, Tipo = "Factorial"),
cbind(puntos_axiales, Tipo = "Axial"),
cbind(puntos_centrales, Tipo = "Central")
)
# Añadir número de punto
simple_ccc$Punto <- 1:nrow(simple_ccc)
# Reordenar columnas
simple_ccc <- simple_ccc[, c("Punto", "Factor_A", "Factor_B", "Tipo")]
# Mostrar diseño
print(simple_ccc)## Punto Factor_A Factor_B Tipo
## 1 1 -1.000000 -1.000000 Factorial
## 2 2 1.000000 -1.000000 Factorial
## 3 3 -1.000000 1.000000 Factorial
## 4 4 1.000000 1.000000 Factorial
## 5 5 -1.414214 0.000000 Axial
## 6 6 1.414214 0.000000 Axial
## 7 7 0.000000 -1.414214 Axial
## 8 8 0.000000 1.414214 Axial
## 9 9 0.000000 0.000000 Central
## 10 10 0.000000 0.000000 Central
## 11 11 0.000000 0.000000 Central
## 12 12 0.000000 0.000000 Central# Visualización del diseño
# Gráfico del Diseño Central Compuesto
plot(simple_ccc$Factor_A, simple_ccc$Factor_B,
pch = ifelse(simple_ccc$Tipo == "Central", 16,
ifelse(simple_ccc$Tipo == "Axial", 17, 15)),
col = ifelse(simple_ccc$Tipo == "Central", "red",
ifelse(simple_ccc$Tipo == "Axial", "blue", "green")),
xlab = "Factor A", ylab = "Factor B",
main = "Diseño Central Compuesto",
xlim = c(-1.5, 1.5), ylim = c(-1.5, 1.5),
cex = 1.5)
abline(h = 0, v = 0, lty = 2)
grid()
legend("topright",
legend = c("Factorial", "Axial", "Central"),
pch = c(15, 17, 16),
col = c("green", "blue", "red"))
Ahora se presenta un DCC para tres factores, una réplica, y cuatro puntos centrales:
design_ccc_3 <- ccd(3,
n0 = 4,
alpha = "rotatable",
randomize = FALSE,
oneblock = TRUE)
# Ver el diseño generado
print(design_ccc_3)## run.order std.order x1.as.is x2.as.is x3.as.is
## 1 1 1 -1.000000 -1.000000 -1.000000
## 2 2 2 1.000000 -1.000000 -1.000000
## 3 3 3 -1.000000 1.000000 -1.000000
## 4 4 4 1.000000 1.000000 -1.000000
## 5 5 5 -1.000000 -1.000000 1.000000
## 6 6 6 1.000000 -1.000000 1.000000
## 7 7 7 -1.000000 1.000000 1.000000
## 8 8 8 1.000000 1.000000 1.000000
## 9 9 9 0.000000 0.000000 0.000000
## 10 10 10 0.000000 0.000000 0.000000
## 11 11 11 0.000000 0.000000 0.000000
## 12 12 12 0.000000 0.000000 0.000000
## 13 1 1 -1.681793 0.000000 0.000000
## 14 2 2 1.681793 0.000000 0.000000
## 15 3 3 0.000000 -1.681793 0.000000
## 16 4 4 0.000000 1.681793 0.000000
## 17 5 5 0.000000 0.000000 -1.681793
## 18 6 6 0.000000 0.000000 1.681793
## 19 7 7 0.000000 0.000000 0.000000
## 20 8 8 0.000000 0.000000 0.000000
## 21 9 9 0.000000 0.000000 0.000000
## 22 10 10 0.000000 0.000000 0.000000
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.is# Mostrar en formato de data frame
design_df <- as.data.frame(design_ccc_3)
print(design_df)## run.order std.order x1 x2 x3
## 1 1 1 -1.000000 -1.000000 -1.000000
## 2 2 2 1.000000 -1.000000 -1.000000
## 3 3 3 -1.000000 1.000000 -1.000000
## 4 4 4 1.000000 1.000000 -1.000000
## 5 5 5 -1.000000 -1.000000 1.000000
## 6 6 6 1.000000 -1.000000 1.000000
## 7 7 7 -1.000000 1.000000 1.000000
## 8 8 8 1.000000 1.000000 1.000000
## 9 9 9 0.000000 0.000000 0.000000
## 10 10 10 0.000000 0.000000 0.000000
## 11 11 11 0.000000 0.000000 0.000000
## 12 12 12 0.000000 0.000000 0.000000
## 13 1 1 -1.681793 0.000000 0.000000
## 14 2 2 1.681793 0.000000 0.000000
## 15 3 3 0.000000 -1.681793 0.000000
## 16 4 4 0.000000 1.681793 0.000000
## 17 5 5 0.000000 0.000000 -1.681793
## 18 6 6 0.000000 0.000000 1.681793
## 19 7 7 0.000000 0.000000 0.000000
## 20 8 8 0.000000 0.000000 0.000000
## 21 9 9 0.000000 0.000000 0.000000
## 22 10 10 0.000000 0.000000 0.000000# Puntos factoriales (esquinas)
puntos_factoriales <- expand.grid(
Factor_A = c(-1, 1),
Factor_B = c(-1, 1),
Factor_C = c(-1, 1)
)
# Puntos axiales
alpha <- (2**3)**(1/4)
puntos_axiales <- rbind(
data.frame(Factor_A = c(-alpha, alpha, 0, 0, 0, 0),
Factor_B = c(0, 0, -alpha, alpha, 0, 0),
Factor_C = c(0, 0, 0, 0, -alpha, alpha))
)
# 4 Puntos centrales
puntos_centrales <- data.frame(
Factor_A = rep(0, 4),
Factor_B = rep(0, 4),
Factor_C = rep(0, 4)
)
# Combinar todos los puntos
simple_ccc <- rbind(
cbind(puntos_factoriales, Tipo = "Factorial"),
cbind(puntos_axiales, Tipo = "Axial"),
cbind(puntos_centrales, Tipo = "Central")
)
# Añadir número de punto
simple_ccc$Punto <- 1:nrow(simple_ccc)
# Reordenar columnas
simple_ccc <- simple_ccc[, c("Punto", "Factor_A", "Factor_B", "Factor_C", "Tipo")]
# Mostrar diseño
print(simple_ccc)## Punto Factor_A Factor_B Factor_C Tipo
## 1 1 -1.000000 -1.000000 -1.000000 Factorial
## 2 2 1.000000 -1.000000 -1.000000 Factorial
## 3 3 -1.000000 1.000000 -1.000000 Factorial
## 4 4 1.000000 1.000000 -1.000000 Factorial
## 5 5 -1.000000 -1.000000 1.000000 Factorial
## 6 6 1.000000 -1.000000 1.000000 Factorial
## 7 7 -1.000000 1.000000 1.000000 Factorial
## 8 8 1.000000 1.000000 1.000000 Factorial
## 9 9 -1.681793 0.000000 0.000000 Axial
## 10 10 1.681793 0.000000 0.000000 Axial
## 11 11 0.000000 -1.681793 0.000000 Axial
## 12 12 0.000000 1.681793 0.000000 Axial
## 13 13 0.000000 0.000000 -1.681793 Axial
## 14 14 0.000000 0.000000 1.681793 Axial
## 15 15 0.000000 0.000000 0.000000 Central
## 16 16 0.000000 0.000000 0.000000 Central
## 17 17 0.000000 0.000000 0.000000 Central
## 18 18 0.000000 0.000000 0.000000 Central# Definir colores y formas por tipo de punto
colores <- ifelse(simple_ccc$Tipo == "Central", "red",
ifelse(simple_ccc$Tipo == "Axial", "blue", "green"))
formas <- ifelse(simple_ccc$Tipo == "Central", "circle",
ifelse(simple_ccc$Tipo == "Axial", "diamond", "square"))
# Gráfico interactivo 3D
fig <- plot_ly(
data = simple_ccc,
x = ~Factor_A, y = ~Factor_B, z = ~Factor_C,
type = "scatter3d",
mode = "markers",
marker = list(size = 6,
color = colores,
symbol = formas),
text = ~paste("Punto:", Punto, "<br>Tipo:", Tipo)
)
# Ajustar layout
fig <- fig %>% layout(
scene = list(
xaxis = list(title = "Factor A"),
yaxis = list(title = "Factor B"),
zaxis = list(title = "Factor C")
),
legend = list(title = list(text = "Tipo de punto"))
)
figDiseño Centrado en Caras
# Diseño Central Compuesto para 2 factores
design_ccf <- ccd(2,
n0 = 2, # Número de puntos centrales
alpha = "faces", # Diseño rotable (alpha = √2)
randomize = FALSE, # Aleatorizar
oneblock = TRUE, # Un bloque
coding = list(x1 ~ (Factor_A - 0)/1,
x2 ~ (Factor_B - 0)/1))
# Ver el diseño generado
print(design_ccf)## run.order std.order Factor_A Factor_B
## 1 1 1 -1 -1
## 2 2 2 1 -1
## 3 3 3 -1 1
## 4 4 4 1 1
## 5 5 5 0 0
## 6 6 6 0 0
## 7 1 1 -1 0
## 8 2 2 1 0
## 9 3 3 0 -1
## 10 4 4 0 1
## 11 5 5 0 0
## 12 6 6 0 0
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Factor_A - 0)/1
## x2 ~ (Factor_B - 0)/1# Mostrar en formato de data frame
design_df <- as.data.frame(design_ccf)
print(design_df)## run.order std.order x1 x2
## 1 1 1 -1 -1
## 2 2 2 1 -1
## 3 3 3 -1 1
## 4 4 4 1 1
## 5 5 5 0 0
## 6 6 6 0 0
## 7 1 1 -1 0
## 8 2 2 1 0
## 9 3 3 0 -1
## 10 4 4 0 1
## 11 5 5 0 0
## 12 6 6 0 0# Puntos factoriales (esquinas)
puntos_factoriales <- expand.grid(
Factor_A = c(-1, 1),
Factor_B = c(-1, 1)
)
# Puntos centrados en cara
alpha <- 1
puntos_axiales <- rbind(
data.frame(Factor_A = c(-alpha, alpha, 0, 0),
Factor_B = c(0, 0, -alpha, alpha))
)
# 4 Puntos centrales
puntos_centrales <- data.frame(
Factor_A = rep(0, 4),
Factor_B = rep(0, 4)
)
# Combinar todos los puntos
simple_ccf <- rbind(
cbind(puntos_factoriales, Tipo = "Factorial"),
cbind(puntos_axiales, Tipo = "Cara"),
cbind(puntos_centrales, Tipo = "Central")
)
# Añadir número de punto
simple_ccf$Punto <- 1:nrow(simple_ccf)
# Reordenar columnas
simple_ccf <- simple_ccf[, c("Punto", "Factor_A", "Factor_B", "Tipo")]
# Mostrar diseño
print(simple_ccf)## Punto Factor_A Factor_B Tipo
## 1 1 -1 -1 Factorial
## 2 2 1 -1 Factorial
## 3 3 -1 1 Factorial
## 4 4 1 1 Factorial
## 5 5 -1 0 Cara
## 6 6 1 0 Cara
## 7 7 0 -1 Cara
## 8 8 0 1 Cara
## 9 9 0 0 Central
## 10 10 0 0 Central
## 11 11 0 0 Central
## 12 12 0 0 Central# Visualización del diseño
# Gráfico del Diseño Centrado en Caras
plot(simple_ccf$Factor_A, simple_ccf$Factor_B,
pch = ifelse(simple_ccf$Tipo == "Central", 16,
ifelse(simple_ccf$Tipo == "Cara", 17, 15)),
col = ifelse(simple_ccc$Tipo == "Central", "red",
ifelse(simple_ccf$Tipo == "Cara", "blue", "green")),
xlab = "Factor A", ylab = "Factor B",
main = "Diseño Centrado en Caras",
xlim = c(-1.5, 1.5), ylim = c(-1.5, 1.5),
cex = 1.5)
abline(h = 0, v = 0, lty = 2)
grid()
legend("topright",
legend = c("Factorial", "Cara", "Central"),
pch = c(15, 17, 16),
col = c("green", "blue", "red"))
Tres factores
design_ccf_3 <- ccd(3,
n0 = 4,
alpha = "faces",
randomize = FALSE,
oneblock = TRUE)
# Ver el diseño generado
print(design_ccf_3)## run.order std.order x1.as.is x2.as.is x3.as.is
## 1 1 1 -1 -1 -1
## 2 2 2 1 -1 -1
## 3 3 3 -1 1 -1
## 4 4 4 1 1 -1
## 5 5 5 -1 -1 1
## 6 6 6 1 -1 1
## 7 7 7 -1 1 1
## 8 8 8 1 1 1
## 9 9 9 0 0 0
## 10 10 10 0 0 0
## 11 11 11 0 0 0
## 12 12 12 0 0 0
## 13 1 1 -1 0 0
## 14 2 2 1 0 0
## 15 3 3 0 -1 0
## 16 4 4 0 1 0
## 17 5 5 0 0 -1
## 18 6 6 0 0 1
## 19 7 7 0 0 0
## 20 8 8 0 0 0
## 21 9 9 0 0 0
## 22 10 10 0 0 0
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.is# Mostrar en formato de data frame
design_df <- as.data.frame(design_ccf_3)
print(design_df)## run.order std.order x1 x2 x3
## 1 1 1 -1 -1 -1
## 2 2 2 1 -1 -1
## 3 3 3 -1 1 -1
## 4 4 4 1 1 -1
## 5 5 5 -1 -1 1
## 6 6 6 1 -1 1
## 7 7 7 -1 1 1
## 8 8 8 1 1 1
## 9 9 9 0 0 0
## 10 10 10 0 0 0
## 11 11 11 0 0 0
## 12 12 12 0 0 0
## 13 1 1 -1 0 0
## 14 2 2 1 0 0
## 15 3 3 0 -1 0
## 16 4 4 0 1 0
## 17 5 5 0 0 -1
## 18 6 6 0 0 1
## 19 7 7 0 0 0
## 20 8 8 0 0 0
## 21 9 9 0 0 0
## 22 10 10 0 0 0# Puntos factoriales (esquinas)
puntos_factoriales <- expand.grid(
Factor_A = c(-1, 1),
Factor_B = c(-1, 1),
Factor_C = c(-1, 1)
)
# Puntos axiales
alpha <- 1
puntos_axiales <- rbind(
data.frame(Factor_A = c(-alpha, alpha, 0, 0, 0, 0),
Factor_B = c(0, 0, -alpha, alpha, 0, 0),
Factor_C = c(0, 0, 0, 0, -alpha, alpha))
)
# 4 Puntos centrales
puntos_centrales <- data.frame(
Factor_A = rep(0, 4),
Factor_B = rep(0, 4),
Factor_C = rep(0, 4)
)
# Combinar todos los puntos
simple_ccf <- rbind(
cbind(puntos_factoriales, Tipo = "Factorial"),
cbind(puntos_axiales, Tipo = "Cara"),
cbind(puntos_centrales, Tipo = "Central")
)
# Añadir número de punto
simple_ccf$Punto <- 1:nrow(simple_ccf)
# Reordenar columnas
simple_ccf <- simple_ccf[, c("Punto", "Factor_A", "Factor_B", "Factor_C", "Tipo")]
# Mostrar diseño
print(simple_ccf)## Punto Factor_A Factor_B Factor_C Tipo
## 1 1 -1 -1 -1 Factorial
## 2 2 1 -1 -1 Factorial
## 3 3 -1 1 -1 Factorial
## 4 4 1 1 -1 Factorial
## 5 5 -1 -1 1 Factorial
## 6 6 1 -1 1 Factorial
## 7 7 -1 1 1 Factorial
## 8 8 1 1 1 Factorial
## 9 9 -1 0 0 Cara
## 10 10 1 0 0 Cara
## 11 11 0 -1 0 Cara
## 12 12 0 1 0 Cara
## 13 13 0 0 -1 Cara
## 14 14 0 0 1 Cara
## 15 15 0 0 0 Central
## 16 16 0 0 0 Central
## 17 17 0 0 0 Central
## 18 18 0 0 0 Central# Definir colores y formas por tipo de punto
colores <- ifelse(simple_ccf$Tipo == "Central", "red",
ifelse(simple_ccf$Tipo == "Cara", "blue", "green"))
formas <- ifelse(simple_ccf$Tipo == "Central", "circle",
ifelse(simple_ccf$Tipo == "Cara", "diamond", "square"))
# Gráfico interactivo 3D
fig <- plot_ly(
data = simple_ccf,
x = ~Factor_A, y = ~Factor_B, z = ~Factor_C,
type = "scatter3d",
mode = "markers",
marker = list(size = 6,
color = colores,
symbol = formas),
text = ~paste("Punto:", Punto, "<br>Tipo:", Tipo)
)
# Ajustar layout
fig <- fig %>% layout(
scene = list(
xaxis = list(title = "Factor A"),
yaxis = list(title = "Factor B"),
zaxis = list(title = "Factor C")
),
legend = list(title = list(text = "Tipo de punto"))
)
figDiseño Box-Behnken
# Diseño Box-Behnken con 3 factores
design_bbd_3 <- bbd(3,
n0 = 4,
randomize = FALSE)
print(design_bbd_3)## run.order std.order x1.as.is x2.as.is x3.as.is
## 1 1 1 -1 -1 0
## 2 2 2 1 -1 0
## 3 3 3 -1 1 0
## 4 4 4 1 1 0
## 5 5 5 -1 0 -1
## 6 6 6 1 0 -1
## 7 7 7 -1 0 1
## 8 8 8 1 0 1
## 9 9 9 0 -1 -1
## 10 10 10 0 1 -1
## 11 11 11 0 -1 1
## 12 12 12 0 1 1
## 13 13 13 0 0 0
## 14 14 14 0 0 0
## 15 15 15 0 0 0
## 16 16 16 0 0 0
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.isdesign_df <- as.data.frame(design_bbd_3)
print(design_df)## run.order std.order x1 x2 x3
## 1 1 1 -1 -1 0
## 2 2 2 1 -1 0
## 3 3 3 -1 1 0
## 4 4 4 1 1 0
## 5 5 5 -1 0 -1
## 6 6 6 1 0 -1
## 7 7 7 -1 0 1
## 8 8 8 1 0 1
## 9 9 9 0 -1 -1
## 10 10 10 0 1 -1
## 11 11 11 0 -1 1
## 12 12 12 0 1 1
## 13 13 13 0 0 0
## 14 14 14 0 0 0
## 15 15 15 0 0 0
## 16 16 16 0 0 0# Visualización 3D
fig <- plot_ly(
data = design_df,
x = ~x1, y = ~x2, z = ~x3,
type = "scatter3d",
mode = "markers",
colors = c("purple", "orange", "red"),
marker = list(size = 6),
hoverinfo = "text"
) %>%
layout(
scene = list(
xaxis = list(title = "Factor A"),
yaxis = list(title = "Factor B"),
zaxis = list(title = "Factor C")
)
)
figAnálisis estadístico
El efecto de cuatro variables independientes tiempo, cantidad de encimas, velocidad de agitación, y concentración inicial de tinte han sido estudiados para obtener la decolorización máxima. A continuación se codifica el diseño de experimentos en un dataframe. Se plantea un diseño central compuesto (DCC).
# Cargar librerías necesarias
library(rsm)
library(ggplot2)
library(car)
library(gridExtra)
# Crear el dataframe con los datos reales del estudio
# (Valores codificados según el diseño experimental utilizado)
decolorizacion <- data.frame(
x1 = c(-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1, -(2**4)**(1/4),(2**4)**(1/4),
0,0,0,0,0,0,0,0,0,0,0),
x2 = c(-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,-1,-1,1,1,0,0,-(2**4)**(1/4),(2**4)**(1/4),
0,0,0,0,0,0,0,0,0),
x3 = c(-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1,
0,0,0,0,-(2**4)**(1/4),(2**4)**(1/4),0,0,0,0,0,0,0),
x4 = c(-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,
-(2**4)**(1/4),(2**4)**(1/4),0,0,0,0,0)
)
# Decodificar variables (valores reales)
decolorizacion$Reaction_time <- 50 + 25 * decolorizacion$x1 # 45-75 min
decolorizacion$Enzyme_amount <- 2.5 + 1.5 * decolorizacion$x2 # 1.5-2.5 mL
decolorizacion$Agitation_speed <- 75 + 25 * decolorizacion$x3 # 25-75 rpm
decolorizacion$dye_concentration <- 100 + 50 * decolorizacion$x4 # 50-100 mg/L
decolorizacion$Y <- c(88.99,90.75,84.33,84.82,81.93,82.43,89.84,89.59,82.68,
83.44,88.74,90.20,91.23,85.15,91.07,85.78,86.58,89.15,
87.54,88.18,87.22,78.40,92.02,87.77,86.00,84.28,85.47,
84.10,83.56) #%
# Ver las primeras filas del dataframe
head(decolorizacion)## x1 x2 x3 x4 Reaction_time Enzyme_amount Agitation_speed dye_concentration
## 1 -1 -1 -1 -1 25 1 50 50
## 2 1 -1 -1 -1 75 1 50 50
## 3 -1 1 -1 -1 25 4 50 50
## 4 1 1 -1 -1 75 4 50 50
## 5 -1 -1 1 -1 25 1 100 50
## 6 1 -1 1 -1 75 1 100 50
## Y
## 1 88.99
## 2 90.75
## 3 84.33
## 4 84.82
## 5 81.93
## 6 82.43Se plantea un modelo de primer orden
model_Y <- rsm(Y ~ FO(x1,x2,x3,x4), data = decolorizacion)
summary(model_Y)##
## Call:
## rsm(formula = Y ~ FO(x1, x2, x3, x4), data = decolorizacion)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.594483 0.638030 135.7217 <2e-16 ***
## x1 -0.062917 0.701349 -0.0897 0.9293
## x2 0.793750 0.701349 1.1317 0.2689
## x3 -0.607083 0.701349 -0.8656 0.3953
## x4 -0.120417 0.701349 -0.1717 0.8651
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.07932, Adjusted R-squared: -0.07413
## F-statistic: 0.5169 on 4 and 24 DF, p-value: 0.724
##
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3, x4) 4 24.409 6.1023 0.5169 0.72403
## Residuals 24 283.329 11.8054
## Lack of fit 20 279.212 13.9606 13.5629 0.01064
## Pure error 4 4.117 1.0293
##
## Direction of steepest ascent (at radius 1):
## x1 x2 x3 x4
## -0.06238713 0.78706941 -0.60197382 -0.11940318
##
## Corresponding increment in original units:
## x1 x2 x3 x4
## -0.06238713 0.78706941 -0.60197382 -0.11940318Segundo orden incompleto
model_Y <- rsm(Y ~ FO(x1,x2,x3,x4) + TWI(x1,x2,x3,x4), data = decolorizacion)
summary(model_Y)##
## Call:
## rsm(formula = Y ~ FO(x1, x2, x3, x4) + TWI(x1, x2, x3, x4), data = decolorizacion)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.594483 0.667685 129.6936 <2e-16 ***
## x1 -0.062917 0.733948 -0.0857 0.9326
## x2 0.793750 0.733948 1.0815 0.2938
## x3 -0.607083 0.733948 -0.8271 0.4190
## x4 -0.120417 0.733948 -0.1641 0.8715
## x1:x2 -0.033125 0.898899 -0.0369 0.9710
## x1:x3 -0.974375 0.898899 -1.0840 0.2927
## x1:x4 -0.728125 0.898899 -0.8100 0.4285
## x2:x3 0.831875 0.898899 0.9254 0.3670
## x2:x4 0.550625 0.898899 0.6126 0.5478
## x3:x4 0.829375 0.898899 0.9227 0.3684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.2438, Adjusted R-squared: -0.1763
## F-statistic: 0.5803 on 10 and 18 DF, p-value: 0.8093
##
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3, x4) 4 24.409 6.1023 0.4720 0.755663
## TWI(x1, x2, x3, x4) 6 50.620 8.4366 0.6526 0.687939
## Residuals 18 232.709 12.9283
## Lack of fit 14 228.592 16.3280 15.8629 0.008165
## Pure error 4 4.117 1.0293
##
## Stationary point of response surface:
## x1 x2 x3 x4
## 11.4984588 15.1846393 0.1587821 -0.9896943
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 1.02049718 -0.01487919 -0.37128198 -0.63433601
##
## $vectors
## [,1] [,2] [,3] [,4]
## x1 0.4751593 0.62298478 0.2777712 -0.55583883
## x2 -0.3899935 0.77992264 -0.1490188 0.46628229
## x3 -0.5958203 0.03098473 -0.4194098 -0.68420282
## x4 -0.5168432 -0.05148384 0.8513119 -0.07409805Segundo orden completo
model_Y <- rsm(Y ~ FO(x1,x2,x3,x4) + TWI(x1,x2,x3,x4) + PQ(x1,x2,x3,x4), data = decolorizacion)
summary(model_Y)##
## Call:
## rsm(formula = Y ~ FO(x1, x2, x3, x4) + TWI(x1, x2, x3, x4) +
## PQ(x1, x2, x3, x4), data = decolorizacion)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.682000 1.492141 56.7520 < 2e-16 ***
## x1 -0.062917 0.681066 -0.0924 0.92771
## x2 0.793750 0.681066 1.1655 0.26331
## x3 -0.607083 0.681066 -0.8914 0.38779
## x4 -0.120417 0.681066 -0.1768 0.86219
## x1:x2 -0.033125 0.834132 -0.0397 0.96888
## x1:x3 -0.974375 0.834132 -1.1681 0.26226
## x1:x4 -0.728125 0.834132 -0.8729 0.39744
## x2:x3 0.831875 0.834132 0.9973 0.33555
## x2:x4 0.550625 0.834132 0.6601 0.51989
## x3:x4 0.829375 0.834132 0.9943 0.33695
## x1^2 0.767104 0.655028 1.1711 0.26111
## x2^2 0.765854 0.655028 1.1692 0.26185
## x3^2 -0.496646 0.655028 -0.7582 0.46090
## x4^2 1.274604 0.655028 1.9459 0.07203 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.4936, Adjusted R-squared: -0.0129
## F-statistic: 0.9745 on 14 and 14 DF, p-value: 0.5189
##
## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3, x4) 4 24.409 6.1023 0.5482 0.703452
## TWI(x1, x2, x3, x4) 6 50.620 8.4366 0.7578 0.614233
## PQ(x1, x2, x3, x4) 4 76.855 19.2138 1.7259 0.200170
## Residuals 14 155.854 11.1324
## Lack of fit 10 151.737 15.1737 14.7415 0.009752
## Pure error 4 4.117 1.0293
##
## Stationary point of response surface:
## x1 x2 x3 x4
## -0.1812208 -0.3273305 -0.5129482 0.2330637
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 1.7629064 0.7516346 0.5817804 -0.7854047
##
## $vectors
## [,1] [,2] [,3] [,4]
## x1 0.4357358 -0.57531620 0.6399602 -0.2638115
## x2 -0.3462799 -0.81569729 -0.4038709 0.2271926
## x3 -0.2998867 -0.02980293 -0.2063566 -0.9309118
## x4 -0.7747855 0.05254514 0.6202872 0.1104095Se obtiene el mejor modelo eliminando los factores que no son significativos, respetando la jerarquía de los factores.
model_Y <- rsm(Y ~ FO(x1,x2,x3,x4) + TWI(x1,x2,x3,x4) + PQ(x1,x2,x3,x4), data = decolorizacion)Supuestos
shapiro.test(model_Y$residuals)##
## Shapiro-Wilk normality test
##
## data: model_Y$residuals
## W = 0.97433, p-value = 0.6813Graficos de contorno y perspectiva
# contour and perspective plots
persp(model_Y, x1~x2, col = "green3", contours = "colors",
zlab = "Specificity",
xlabs=c("Reaction time", "Enzyme amount"),
cex.lab=1.6,
cex.axes=1.5,
theta = 125,
phi = 20)## Warning in persp.default(dat$x, dat$y, dat$z, zlim = dat$zlim, theta = theta, :
## "cex.axes" is not a graphical parameter## Warning in persp.default(dat$x, dat$y, dat$z, xlab = dat$labs[1], ylab =
## dat$labs[2], : "cex.axes" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "cex.axes" is not a graphical
## parameter
persp(model_Y, x1~x3, col = "green3", contours = "colors",
zlab = "Specificity",
xlabs=c("Reaction time", "Agitation speed"),
cex.lab=1.6,
cex.axes=1.5,
theta = 40,
phi = 20)## Warning in persp.default(dat$x, dat$y, dat$z, zlim = dat$zlim, theta = theta, :
## "cex.axes" is not a graphical parameter## Warning in persp.default(dat$x, dat$y, dat$z, xlab = dat$labs[1], ylab =
## dat$labs[2], : "cex.axes" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "cex.axes" is not a graphical
## parameter
persp(model_Y, x2~x3, col = "green3", contours = "colors",
zlab = "Specificity",
xlabs=c("Enzyme amount", "Agitation speed"),
cex.lab=1.6,
cex.axes=1.5,
theta = 40,
phi = 20)## Warning in persp.default(dat$x, dat$y, dat$z, zlim = dat$zlim, theta = theta, :
## "cex.axes" is not a graphical parameter## Warning in persp.default(dat$x, dat$y, dat$z, xlab = dat$labs[1], ylab =
## dat$labs[2], : "cex.axes" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "cex.axes" is not a graphical
## parameter
persp(model_Y, x3~x4, col = "green3", contours = "colors",
zlab = "Yield (%)",
xlabs=c("Agitation speed", "Dye concentration"),
cex.lab=1.6,
cex.axes=1.5,
theta = 40,
phi = 20)## Warning in persp.default(dat$x, dat$y, dat$z, zlim = dat$zlim, theta = theta, :
## "cex.axes" is not a graphical parameter## Warning in persp.default(dat$x, dat$y, dat$z, xlab = dat$labs[1], ylab =
## dat$labs[2], : "cex.axes" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "cex.axes" is not a graphical
## parameter