#--------------------------------------# #--- Response Surface Modeling in R ---# #--------------------------------------# #First, install and load the "rsm" package # install.packages("rsm") library(rsm) # Example generating a Box-Behnken design with three factors and two center points (no) bbd(3, n0 = 2, coding = list(x1 ~ (Force - 20)/3, x2 ~ (Rate - 50)/10, x3 ~ Polish - 4)) # Example data set data = ChemReact plot(data) # The data set was collected in two blocks. # Block1 is a 2-level, two-factor factorial design with three repeated center points. # Block 2 is the Central Composite Design (circomscribed) with 3 center points. # The variables are Time = 85 +/- 5 and Temp = 175 +/- 5, # Thus, the coded variables are x1 = (Time-85)/5 and x2 = (Temp-175)/5 CR <- coded.data(ChemReact, x1 ~ (Time - 85)/5, x2 ~ (Temp - 175)/5) CR[1:7,] # Note: If the data are already coded, use as.coded.data() to convert to the proper coded data object # Let's work as though the first block (full factorial) has been finished, # and we'll fit a linear model, first order (FO), to it (Yield is the response) CR.rsm1 <- rsm(Yield ~ FO(x1, x2), data = CR, subset = (Block == "B1")) summary(CR.rsm1) #The fit is not very good. Let's include the interaction term (TWI) and update the model, or start over with a new model (these two lines do the same thing) CR.rsm1.5 <- update(CR.rsm1, . ~ . + TWI(x1, x2)) CR.rsm1.5 <- rsm(Yield ~ FO(x1, x2)+TWI(x1, x2), data = CR, subset = (Block == "B1")) summary(CR.rsm1.5) #This is no better! The reason is the strong quadratic response, with the peak near the center. # Now let's assume the second block has been collected. We use the SO (second order) function, which includes FO and TWI CR.rsm2 <- rsm(Yield ~ Block + SO(x1, x2), data = CR) summary(CR.rsm2) # The secondary point is a maximum (both eigenvalues are negative) and within the experimental design range (no extrapolation) # Also note that the block is significant, meaning that the processes shifted between the first set of data and the second. This is not good. The coefficient is -4.5, meaning the yield shifted down by 4.5% between the two blocks - a more significant effect than either temperatue or time! This is most easily seen by looking at the repeat center points. # We can plot the fitted response as a contour plot. contour(CR.rsm2, ~ x1 + x2, at = summary(CR.rsm2)$canonical$xs)