# Install and load package
#install.packages("rsm") # run once
library(rsm)
# Generate Box-Behnken Design (3 factors, 3 center points)
design <- bbd(k = 3, n0 = 3)
# Convert coded variables to actual values
design$yeast <- 5 + (design$x1 * 3)
design$time <- 2.5 + (design$x2 * 0.5)
design$temp <- 32.5 + (design$x3 * 2.5)
# Add tensile strength data (example dummy values)
design$ethanol <- c(8.5, 6.2, 8.6, 7.1, 5.8, 7.3, 5.9,
9.0, 6.8, 6.5, 7.2, 8.4, 8.8, 8.7, 8.2)
# Fit second-order response surface model - Quadratic model (curved relationship)
# Quadratic terms (x^2) are included to capture curvature
# and allow identification of optimal conditions (max/min).
# Without these terms, the model is linear and cannot locate an optimum.
model <- rsm(ethanol ~ SO(x1, x2, x3), data = design)
# View results
summary(model)
##
## Call:
## rsm(formula = ethanol ~ SO(x1, x2, x3), data = design)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.96667 0.73265 9.5088 0.0002175 ***
## x1 -0.57500 0.44866 -1.2816 0.2561931
## x2 0.20000 0.44866 0.4458 0.6743988
## x3 0.07500 0.44866 0.1672 0.8737898
## x1:x2 -1.07500 0.63450 -1.6943 0.1509896
## x1:x3 0.17500 0.63450 0.2758 0.7937300
## x2:x3 -0.02500 0.63450 -0.0394 0.9700952
## x1^2 0.10417 0.66040 0.1577 0.8808401
## x2^2 0.70417 0.66040 1.0663 0.3350589
## x3^2 0.25417 0.66040 0.3849 0.7161634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.5475, Adjusted R-squared: -0.267
## F-statistic: 0.6722 on 9 and 5 DF, p-value: 0.7155
##
## Analysis of Variance Table
##
## Response: ethanol
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 3.0100 1.00333 0.6231 0.6301
## TWI(x1, x2, x3) 3 4.7475 1.58250 0.9827 0.4712
## PQ(x1, x2, x3) 3 1.9842 0.66139 0.4107 0.7527
## Residuals 5 8.0517 1.61033
## Lack of fit 3 3.8050 1.26833 0.5973 0.6751
## Pure error 2 4.2467 2.12333
##
## Stationary point of response surface:
## x1 x2 x3
## -0.67634543 -0.65733518 0.05297128
##
## Stationary point in original units:
## x1.as.is x2.as.is x3.as.is
## -0.67634543 -0.65733518 0.05297128
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 1.0236825 0.2603032 -0.2214857
##
## $vectors
## [,1] [,2] [,3]
## x1 0.50843184 0.08849098 0.8565433
## x2 -0.85810765 0.13495153 0.4954183
## x3 0.07175178 0.98689282 -0.1445484
# Perform ANOVA to evaluate significance of factors
anova(model)
## Analysis of Variance Table
##
## Response: ethanol
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 3.0100 1.00333 0.6231 0.6301
## TWI(x1, x2, x3) 3 4.7475 1.58250 0.9827 0.4712
## PQ(x1, x2, x3) 3 1.9842 0.66139 0.4107 0.7527
## Residuals 5 8.0517 1.61033
# Anova Interpretation:
# p > 0.05 → factors are not statistically significant
# Lack-of-fit p > 0.05 → model is adequate
# Find stationary point (optimal condition)
canonical(model)
## $xs
## x1 x2 x3
## -0.67634543 -0.65733518 0.05297128
##
## $eigen
## eigen() decomposition
## $values
## [1] 1.0236825 0.2603032 -0.2214857
##
## $vectors
## [,1] [,2] [,3]
## x1 0.50843184 0.08849098 0.8565433
## x2 -0.85810765 0.13495153 0.4954183
## x3 0.07175178 0.98689282 -0.1445484
# Canonical analysis results:
# xs = stationary point of the response surface in coded variables
# Eigenvalues determine whether the point is a maximum, minimum, or saddle
#Now let us do a contour plot of the model
canon <- canonical(model)
contour(model, ~ x1 + x2,
at = list(x3 = canon$xs[3]),
image = TRUE,
col.regions = topo.colors(20), # ✅ correct
xlabs = c("Yeast (coded)", "Time (coded)"),
main = "Contour Plot: Yeast vs Time")
## Warning in plot.window(...): "col.regions" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "col.regions" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in box(...): "col.regions" is not a graphical parameter
## Warning in title(...): "col.regions" is not a graphical parameter
#to determine optimum point for yeast vs time
points(canon$xs[1], canon$xs[2],
pch = 19, col = "red", cex = 1.5)
text(canon$xs[1], canon$xs[2],
labels = "Optimum",
pos = 4, col = "red")
#contour plot of yeast vs temperature
contour(model, ~ x1 + x3,
at = list(x2 = canon$xs[2]), # fix time
image = TRUE,
col.regions = topo.colors(20), # ✅ correct color argument
xlabs = c("Yeast (coded)", "Temperature (coded)"),
main = "Contour Plot: Yeast vs Temperature")
## Warning in plot.window(...): "col.regions" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "col.regions" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in box(...): "col.regions" is not a graphical parameter
## Warning in title(...): "col.regions" is not a graphical parameter
# Add optimum point
points(canon$xs[1], canon$xs[3],
pch = 19, col = "red", cex = 1.5)
text(canon$xs[1], canon$xs[3],
labels = "Optimum",
pos = 4, col = "red")
#contour plot of time vs temp
contour(model, ~ x2 + x3,
at = list(x1 = canon$xs[1]), # fix Yeast
image = TRUE,
col.regions = topo.colors(20),
xlabs = c("Time (coded)", "Temperature (coded)"),
main = "Contour Plot: Time vs Temperature")
## Warning in plot.window(...): "col.regions" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "col.regions" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in box(...): "col.regions" is not a graphical parameter
## Warning in title(...): "col.regions" is not a graphical parameter
# Add optimum point
points(canon$xs[2], canon$xs[3],
pch = 19, col = "red", cex = 1.5)
text(canon$xs[2], canon$xs[3],
labels = "Optimum",
pos = 4, col = "red")
#if you want to visualize the above in one screen
canon <- canonical(model)
# Arrange 3 plots in one row
par(mfrow = c(1, 3))
# 1. Yeast vs Time
contour(model, ~ x1 + x2,
at = list(x3 = canon$xs[3]),
image = TRUE,
col.regions = topo.colors(20),
xlabs = c("Yeast", "Time"),
main = "Yeast vs Time")
## Warning in plot.window(...): "col.regions" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "col.regions" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in box(...): "col.regions" is not a graphical parameter
## Warning in title(...): "col.regions" is not a graphical parameter
points(canon$xs[1], canon$xs[2], pch = 19, col = "red")
# 2. Yeast vs Temperature
contour(model, ~ x1 + x3,
at = list(x2 = canon$xs[2]),
image = TRUE,
col.regions = topo.colors(20),
xlabs = c("Yeast", "Temperature"),
main = "Yeast vs Temperature")
## Warning in plot.window(...): "col.regions" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "col.regions" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in box(...): "col.regions" is not a graphical parameter
## Warning in title(...): "col.regions" is not a graphical parameter
points(canon$xs[1], canon$xs[3], pch = 19, col = "red")
# 3. Time vs Temperature
contour(model, ~ x2 + x3,
at = list(x1 = canon$xs[1]),
image = TRUE,
col.regions = topo.colors(20),
xlabs = c("Time", "Temperature"),
main = "Time vs Temperature")
## Warning in plot.window(...): "col.regions" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "col.regions" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in box(...): "col.regions" is not a graphical parameter
## Warning in title(...): "col.regions" is not a graphical parameter
points(canon$xs[2], canon$xs[3], pch = 19, col = "red")
# Reset layout
par(mfrow = c(1,1))
#3D surface RSM plots
#Yeast vs time (fix temp) RSM 3D plot
canon <- canonical(model)
persp(model, ~ x1 + x2,
at = list(x3 = canon$xs[3]), # fix temp
zlab = "Strength",
col = "lightblue",
theta = 30,
phi = 25,
main = "3D Surface: Yeast vs Time",
ticktype = "detailed")
# Yeast vs Temp (fix Time) 3D RSM plot
persp(model, ~ x1 + x3,
at = list(x2 = canon$xs[2]), # fix time
zlab = "Strength",
col = "lightgreen",
theta = 30,
phi = 25,
main = "3D Surface: Yeast vs Temp",
ticktype = "detailed")
# Time vs Temp (fix Yeast) 3D RSM plot
persp(model, ~ x2 + x3,
at = list(x1 = canon$xs[1]), # fix yeast
zlab = "Strength",
col = "lightpink",
theta = 30,
phi = 25,
main = "3D Surface: Time vs Temp",
ticktype = "detailed")
#All 3D contour plot for RSM in a single plot
par(mfrow = c(1, 3))
# Yeast vs Time
persp(model, ~ x1 + x2,
at = list(x3 = canon$xs[3]),
zlab = "Strength", col = "lightblue",
theta = 30, phi = 25, main = "Yeast vs Time")
# Yeast vs Temp
persp(model, ~ x1 + x3,
at = list(x2 = canon$xs[2]),
zlab = "Strength", col = "lightgreen",
theta = 30, phi = 25, main = "Yeast vs Temp")
# Time vs Temp
persp(model, ~ x2 + x3,
at = list(x1 = canon$xs[1]),
zlab = "Strength", col = "lightpink",
theta = 30, phi = 25, main = "Time vs Temp")
par(mfrow = c(1,1))