# 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$clay <- 80 + (design$x1 * 20)
design$water <- 20 + (design$x2 * 5)
design$time <- 4 + (design$x3 * 2)
# Add tensile strength data (example dummy values)
design$strength <- 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(strength ~ SO(x1, x2, x3), data = design)
# View results
summary(model)
##
## Call:
## rsm(formula = strength ~ SO(x1, x2, x3), data = design)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.16667 0.38723 21.0901 4.441e-06 ***
## x1 0.25000 0.23713 1.0543 0.340009
## x2 0.48750 0.23713 2.0559 0.094932 .
## x3 0.03750 0.23713 0.1581 0.880533
## x1:x2 0.35000 0.33535 1.0437 0.344438
## x1:x3 -0.45000 0.33535 -1.3419 0.237344
## x2:x3 0.12500 0.33535 0.3727 0.724616
## x1^2 -1.69583 0.34904 -4.8585 0.004639 **
## x2^2 0.37917 0.34904 1.0863 0.326916
## x3^2 0.12917 0.34904 0.3701 0.726495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.8736, Adjusted R-squared: 0.6461
## F-statistic: 3.839 on 9 and 5 DF, p-value: 0.07627
##
## Analysis of Variance Table
##
## Response: strength
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 2.4125 0.8042 1.7877 0.26588
## TWI(x1, x2, x3) 3 1.3625 0.4542 1.0096 0.46139
## PQ(x1, x2, x3) 3 11.7692 3.9231 8.7211 0.01976
## Residuals 5 2.2492 0.4498
## Lack of fit 3 0.8025 0.2675 0.3698 0.78687
## Pure error 2 1.4467 0.7233
##
## Stationary point of response surface:
## x1 x2 x3
## -0.01434532 -0.66090242 0.14964157
##
## Stationary point in original units:
## x1.as.is x2.as.is x3.as.is
## -0.01434532 -0.66090242 0.14964157
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.4014730 0.1497106 -1.7386836
##
## $vectors
## [,1] [,2] [,3]
## x1 0.06343108 -0.1346532 0.98886047
## x2 0.98284680 -0.1635065 -0.08531001
## x3 0.17317236 0.9773097 0.12197202
# Perform ANOVA to evaluate significance of factors
anova(model)
## Analysis of Variance Table
##
## Response: strength
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 2.4125 0.8042 1.7877 0.26588
## TWI(x1, x2, x3) 3 1.3625 0.4542 1.0096 0.46139
## PQ(x1, x2, x3) 3 11.7692 3.9231 8.7211 0.01976 *
## Residuals 5 2.2492 0.4498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 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)
## Near-stationary-ridge situation detected -- stationary point altered
## Change 'threshold' if this is not what you intend
## $xs
## x1 x2 x3
## 0.02015815 -0.61900558 -0.10078389
##
## $eigen
## eigen() decomposition
## $values
## [1] 0.401473 0.000000 -1.738684
##
## $vectors
## [,1] [,2] [,3]
## x1 0.06343108 -0.1346532 0.98886047
## x2 0.98284680 -0.1635065 -0.08531001
## x3 0.17317236 0.9773097 0.12197202
# 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)
## Near-stationary-ridge situation detected -- stationary point altered
## Change 'threshold' if this is not what you intend
contour(model, ~ x1 + x2,
at = list(x3 = canon$xs[3]),
image = TRUE,
col.regions = topo.colors(20), # ✅ correct
xlabs = c("Clay (coded)", "Water (coded)"),
main = "Contour Plot: Clay vs Water")
## 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 clay vs water
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 clay vs time
contour(model, ~ x1 + x3,
at = list(x2 = canon$xs[2]), # fix water
image = TRUE,
col.regions = topo.colors(20), # ✅ correct color argument
xlabs = c("Clay (coded)", "Time (coded)"),
main = "Contour Plot: Clay 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
# 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 water vs time
contour(model, ~ x2 + x3,
at = list(x1 = canon$xs[1]), # fix clay
image = TRUE,
col.regions = topo.colors(20),
xlabs = c("Water (coded)", "Time (coded)"),
main = "Contour Plot: Water 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
# 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)
## Near-stationary-ridge situation detected -- stationary point altered
## Change 'threshold' if this is not what you intend
# Arrange 3 plots in one row
par(mfrow = c(1, 3))
# 1. Clay vs Water
contour(model, ~ x1 + x2,
at = list(x3 = canon$xs[3]),
image = TRUE,
col.regions = topo.colors(20),
xlabs = c("Clay", "Water"),
main = "Clay vs Water")
## 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. Clay vs Time
contour(model, ~ x1 + x3,
at = list(x2 = canon$xs[2]),
image = TRUE,
col.regions = topo.colors(20),
xlabs = c("Clay", "Time"),
main = "Clay 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[3], pch = 19, col = "red")
# 3. Water vs Time
contour(model, ~ x2 + x3,
at = list(x1 = canon$xs[1]),
image = TRUE,
col.regions = topo.colors(20),
xlabs = c("Water", "Time"),
main = "Water 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[2], canon$xs[3], pch = 19, col = "red")
# Reset layout
par(mfrow = c(1,1))
#3D surface RSM plots
#Clay vs water (fix time) RSM 3D plot
canon <- canonical(model)
## Near-stationary-ridge situation detected -- stationary point altered
## Change 'threshold' if this is not what you intend
persp(model, ~ x1 + x2,
at = list(x3 = canon$xs[3]), # fix time
zlab = "Strength",
col = "lightblue",
theta = 30,
phi = 25,
main = "3D Surface: Clay vs Water",
ticktype = "detailed")
# Clay vs Time (fix Water) 3D RSM plot
persp(model, ~ x1 + x3,
at = list(x2 = canon$xs[2]), # fix water
zlab = "Strength",
col = "lightgreen",
theta = 30,
phi = 25,
main = "3D Surface: Clay vs Time",
ticktype = "detailed")
# Water vs Time (fix clay) 3D RSM plot
persp(model, ~ x2 + x3,
at = list(x1 = canon$xs[1]), # fix clay
zlab = "Strength",
col = "lightpink",
theta = 30,
phi = 25,
main = "3D Surface: Water vs Time",
ticktype = "detailed")
#All 3D contour plot for RSM in a single plot
par(mfrow = c(1, 3))
# Clay vs Water
persp(model, ~ x1 + x2,
at = list(x3 = canon$xs[3]),
zlab = "Strength", col = "lightblue",
theta = 30, phi = 25, main = "Clay vs Water")
# Clay vs Time
persp(model, ~ x1 + x3,
at = list(x2 = canon$xs[2]),
zlab = "Strength", col = "lightgreen",
theta = 30, phi = 25, main = "Clay vs Time")
# Water vs Time
persp(model, ~ x2 + x3,
at = list(x1 = canon$xs[1]),
zlab = "Strength", col = "lightpink",
theta = 30, phi = 25, main = "Water vs Time")
par(mfrow = c(1,1))