# 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)

# Build RSM model
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)  7.36667    0.95129  7.7439 0.000574 ***
## x1          -0.58750    0.58254 -1.0085 0.359493    
## x2           0.22500    0.58254  0.3862 0.715209    
## x3           0.21250    0.58254  0.3648 0.730196    
## x1:x2       -0.10000    0.82384 -0.1214 0.908115    
## x1:x3        0.22500    0.82384  0.2731 0.795688    
## x2:x3       -0.05000    0.82384 -0.0607 0.953956    
## x1^2        -0.12083    0.85748 -0.1409 0.893437    
## x2^2         0.15417    0.85748  0.1798 0.864376    
## x3^2         0.27917    0.85748  0.3256 0.757935    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.2371, Adjusted R-squared:  -1.136 
## F-statistic: 0.1727 on 9 and 5 DF,  p-value: 0.9885
## 
## Analysis of Variance Table
## 
## Response: strength
##                 Df  Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3)   3  3.5275 1.17583  0.4331 0.7387
## TWI(x1, x2, x3)  3  0.2525 0.08417  0.0310 0.9918
## PQ(x1, x2, x3)   3  0.4392 0.14639  0.0539 0.9817
## Residuals        5 13.5742 2.71483               
## Lack of fit      3  8.7275 2.90917  1.2005 0.4845
## Pure error       2  4.8467 2.42333               
## 
## Stationary point of response surface:
##         x1         x2         x3 
## -1.7236532 -1.2560726  0.2015254 
## 
## Stationary point in original units:
##   x1.as.is   x2.as.is   x3.as.is 
## -1.7236532 -1.2560726  0.2015254 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1]  0.3171666  0.1513856 -0.1560523
## 
## $vectors
##          [,1]        [,2]       [,3]
## x1  0.2664341 -0.07273493  0.9611048
## x2 -0.2254574  0.96478235  0.1355138
## x3  0.9371136  0.25279368 -0.2406523
# Find optimal conditions
canonical(model)
## $xs
##         x1         x2         x3 
## -1.7236532 -1.2560726  0.2015254 
## 
## $eigen
## eigen() decomposition
## $values
## [1]  0.3171666  0.1513856 -0.1560523
## 
## $vectors
##          [,1]        [,2]       [,3]
## x1  0.2664341 -0.07273493  0.9611048
## x2 -0.2254574  0.96478235  0.1355138
## x3  0.9371136  0.25279368 -0.2406523