RSM
denniskiprotichsang
2023-07-13
####import the data into r
library("rsm")
ChemReact## Time Temp Block Yield
## 1 80.00 170.00 B1 80.5
## 2 80.00 180.00 B1 81.5
## 3 90.00 170.00 B1 82.0
## 4 90.00 180.00 B1 83.5
## 5 85.00 175.00 B1 83.9
## 6 85.00 175.00 B1 84.3
## 7 85.00 175.00 B1 84.0
## 8 85.00 175.00 B2 79.7
## 9 85.00 175.00 B2 79.8
## 10 85.00 175.00 B2 79.5
## 11 92.07 175.00 B2 78.4
## 12 77.93 175.00 B2 75.6
## 13 85.00 182.07 B2 78.5
## 14 85.00 167.93 B2 77.0the coded data for the two blocks , block 1 is chmreact 1 and chemreact2
CR1<-coded.data(ChemReact1,x1~(Time-85)/5,x2~(Temp-175)/5)
CR1## Time Temp Yield
## 1 80 170 80.5
## 2 80 180 81.5
## 3 90 170 82.0
## 4 90 180 83.5
## 5 85 175 83.9
## 6 85 175 84.3
## 7 85 175 84.0
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Time - 85)/5
## x2 ~ (Temp - 175)/5data frame the data
as.data.frame(CR1)## x1 x2 Yield
## 1 -1 -1 80.5
## 2 -1 1 81.5
## 3 1 -1 82.0
## 4 1 1 83.5
## 5 0 0 83.9
## 6 0 0 84.3
## 7 0 0 84.0create a Box Behnken of 3 faactors and two center poinnts
bbd(3, n0 = 2, coding =list(x1 ~ (Force - 20)/3, x2 ~ (Rate - 50)/10, x3 ~ Polish - 4))## run.order std.order Force Rate Polish
## 1 1 9 20 40 3
## 2 2 10 20 60 3
## 3 3 6 23 50 3
## 4 4 14 20 50 4
## 5 5 2 23 40 4
## 6 6 11 20 40 5
## 7 7 7 17 50 5
## 8 8 5 17 50 3
## 9 9 1 17 40 4
## 10 10 3 17 60 4
## 11 11 13 20 50 4
## 12 12 4 23 60 4
## 13 13 8 23 50 5
## 14 14 12 20 60 5
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Force - 20)/3
## x2 ~ (Rate - 50)/10
## x3 ~ Polish - 4fit the first order response surface model to the first block
CR1.rsm <- rsm(Yield ~ FO(x1, x2), data = CR1)
summary(CR1.rsm)##
## Call:
## rsm(formula = Yield ~ FO(x1, x2), data = CR1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82.81429 0.54719 151.3456 1.143e-08 ***
## x1 0.87500 0.72386 1.2088 0.2933
## x2 0.62500 0.72386 0.8634 0.4366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.3555, Adjusted R-squared: 0.0333
## F-statistic: 1.103 on 2 and 4 DF, p-value: 0.4153
##
## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 4.6250 2.3125 1.1033 0.41534
## Residuals 4 8.3836 2.0959
## Lack of fit 2 8.2969 4.1485 95.7335 0.01034
## Pure error 2 0.0867 0.0433
##
## Direction of steepest ascent (at radius 1):
## x1 x2
## 0.8137335 0.5812382
##
## Corresponding increment in original units:
## Time Temp
## 4.068667 2.906191###The summary of this model indicates that there is a significant lack of fit (p = 0.01034) ###This suggests that it may be a good idea to test a higher order model
###lets model with two way interaction
CR1.rsmi <- update(CR1.rsm, . ~ . + TWI(x1, x2))
summary(CR1.rsmi)##
## Call:
## rsm(formula = Yield ~ FO(x1, x2) + TWI(x1, x2), data = CR1)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82.81429 0.62948 131.5604 9.683e-07 ***
## x1 0.87500 0.83272 1.0508 0.3705
## x2 0.62500 0.83272 0.7506 0.5074
## x1:x2 0.12500 0.83272 0.1501 0.8902
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.3603, Adjusted R-squared: -0.2793
## F-statistic: 0.5633 on 3 and 3 DF, p-value: 0.6755
##
## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2) 2 4.6250 2.3125 0.8337 0.515302
## TWI(x1, x2) 1 0.0625 0.0625 0.0225 0.890202
## Residuals 3 8.3211 2.7737
## Lack of fit 1 8.2344 8.2344 190.0247 0.005221
## Pure error 2 0.0867 0.0433
##
## Stationary point of response surface:
## x1 x2
## -5 -7
##
## Stationary point in original units:
## Time Temp
## 60 140
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.0625 -0.0625
##
## $vectors
## [,1] [,2]
## x1 0.7071068 -0.7071068
## x2 0.7071068 0.7071068###The p-value for the lack of fit test is still very small (p = 0.005221). To investigate further, more data is needed so we can combine two blocks of the greater ChemReact experiment (ChemReact1 + ChemReact2):
###combining the two blocks
( CR2 <- djoin(CR1, ChemReact2) )## Time Temp Yield Block
## 1 80.00 170.00 80.5 1
## 2 80.00 180.00 81.5 1
## 3 90.00 170.00 82.0 1
## 4 90.00 180.00 83.5 1
## 5 85.00 175.00 83.9 1
## 6 85.00 175.00 84.3 1
## 7 85.00 175.00 84.0 1
## 8 85.00 175.00 79.7 2
## 9 85.00 175.00 79.8 2
## 10 85.00 175.00 79.5 2
## 11 92.07 175.00 78.4 2
## 12 77.93 175.00 75.6 2
## 13 85.00 182.07 78.5 2
## 14 85.00 167.93 77.0 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Time - 85)/5
## x2 ~ (Temp - 175)/5###The larger dataset allows us to fit a full second-order model to the data. This can be accomplished by using the following code:
CR2.rsm <- rsm(Yield ~ Block + SO(x1, x2), data = CR2)
summary(CR2.rsm)##
## Call:
## rsm(formula = Yield ~ Block + SO(x1, x2), data = CR2)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.095427 0.079631 1056.067 < 2.2e-16 ***
## Block2 -4.457530 0.087226 -51.103 2.877e-10 ***
## x1 0.932541 0.057699 16.162 8.444e-07 ***
## x2 0.577712 0.057699 10.012 2.122e-05 ***
## x1:x2 0.125000 0.081592 1.532 0.1694
## x1^2 -1.308555 0.060064 -21.786 1.083e-07 ***
## x2^2 -0.933442 0.060064 -15.541 1.104e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9981, Adjusted R-squared: 0.9964
## F-statistic: 607.2 on 6 and 7 DF, p-value: 3.811e-09
##
## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Block 1 69.531 69.531 2611.0950 2.879e-10
## FO(x1, x2) 2 9.626 4.813 180.7341 9.450e-07
## TWI(x1, x2) 1 0.063 0.063 2.3470 0.1694
## PQ(x1, x2) 2 17.791 8.896 334.0539 1.135e-07
## Residuals 7 0.186 0.027
## Lack of fit 3 0.053 0.018 0.5307 0.6851
## Pure error 4 0.133 0.033
##
## Stationary point of response surface:
## x1 x2
## 0.3722954 0.3343802
##
## Stationary point in original units:
## Time Temp
## 86.86148 176.67190
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.9233027 -1.3186949
##
## $vectors
## [,1] [,2]
## x1 -0.1601375 -0.9870947
## x2 -0.9870947 0.1601375#####In this model, the lack of fit is not significant (p = 0.6851) ###The stationary point values are like the coordinates for the center of the plot ###In the Eigenanalysis, both eigenvalues are negative ###This is indicative that the stationary point is a maximum ##This is an ideal situation in which the optimal conditions for this system are in close range of the stationary point ##The next step would be to collect more data around this estimated optimum point when ##Time = 86.86148 min ##Temp = 176.67190 C
lets visualize the contour plot
par(mfrow = c(2,3))
contour(CR2.rsm, ~ x1 + x2, image = TRUE, at = summary(CR2.rsm)$canonical$xs)
#### import the data in r of an expriment of an paper helicopter
flight
heli## block A R W L ave logSD
## 1 1 11.8 2.26 1.00 1.5 367 72
## 2 1 13.0 2.26 1.00 1.5 369 72
## 3 1 11.8 2.78 1.00 1.5 374 74
## 4 1 13.0 2.78 1.00 1.5 370 79
## 5 1 11.8 2.26 1.50 1.5 372 72
## 6 1 13.0 2.26 1.50 1.5 355 81
## 7 1 11.8 2.78 1.50 1.5 397 72
## 8 1 13.0 2.78 1.50 1.5 377 99
## 9 1 11.8 2.26 1.00 2.5 350 90
## 10 1 13.0 2.26 1.00 2.5 373 86
## 11 1 11.8 2.78 1.00 2.5 358 92
## 12 1 13.0 2.78 1.00 2.5 363 112
## 13 1 11.8 2.26 1.50 2.5 344 76
## 14 1 13.0 2.26 1.50 2.5 355 69
## 15 1 11.8 2.78 1.50 2.5 370 91
## 16 1 13.0 2.78 1.50 2.5 362 71
## 17 1 12.4 2.52 1.25 2.0 377 51
## 18 1 12.4 2.52 1.25 2.0 375 74
## 19 2 11.2 2.52 1.25 2.0 361 111
## 20 2 13.6 2.52 1.25 2.0 364 93
## 21 2 12.4 2.00 1.25 2.0 355 100
## 22 2 12.4 3.04 1.25 2.0 373 80
## 23 2 12.4 2.52 0.75 2.0 361 71
## 24 2 12.4 2.52 1.75 2.0 360 98
## 25 2 12.4 2.52 1.25 1.0 380 69
## 26 2 12.4 2.52 1.25 3.0 360 74
## 27 2 12.4 2.52 1.25 2.0 370 86
## 28 2 12.4 2.52 1.25 2.0 368 74
## 29 2 12.4 2.52 1.25 2.0 369 89
## 30 2 12.4 2.52 1.25 2.0 366 76
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (A - 12.4)/0.6
## x2 ~ (R - 2.52)/0.26
## x3 ~ (W - 1.25)/0.25
## x4 ~ (L - 2)/0.5heli.rsm <- rsm(ave ~ block + SO(x1, x2, x3, x4), data = heli)
summary(heli.rsm)##
## Call:
## rsm(formula = ave ~ block + SO(x1, x2, x3, x4), data = heli)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 372.800000 1.506375 247.4815 < 2.2e-16 ***
## block2 -2.950000 1.207787 -2.4425 0.0284522 *
## x1 -0.083333 0.636560 -0.1309 0.8977075
## x2 5.083333 0.636560 7.9856 1.398e-06 ***
## x3 0.250000 0.636560 0.3927 0.7004292
## x4 -6.083333 0.636560 -9.5566 1.633e-07 ***
## x1:x2 -2.875000 0.779623 -3.6877 0.0024360 **
## x1:x3 -3.750000 0.779623 -4.8100 0.0002773 ***
## x1:x4 4.375000 0.779623 5.6117 6.412e-05 ***
## x2:x3 4.625000 0.779623 5.9324 3.657e-05 ***
## x2:x4 -1.500000 0.779623 -1.9240 0.0749257 .
## x3:x4 -2.125000 0.779623 -2.7257 0.0164099 *
## x1^2 -2.037500 0.603894 -3.3739 0.0045424 **
## x2^2 -1.662500 0.603894 -2.7530 0.0155541 *
## x3^2 -2.537500 0.603894 -4.2019 0.0008873 ***
## x4^2 -0.162500 0.603894 -0.2691 0.7917877
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9555, Adjusted R-squared: 0.9078
## F-statistic: 20.04 on 15 and 14 DF, p-value: 6.54e-07
##
## Analysis of Variance Table
##
## Response: ave
## Df Sum Sq Mean Sq F value Pr(>F)
## block 1 16.81 16.81 1.7281 0.209786
## FO(x1, x2, x3, x4) 4 1510.00 377.50 38.8175 1.965e-07
## TWI(x1, x2, x3, x4) 6 1114.00 185.67 19.0917 5.355e-06
## PQ(x1, x2, x3, x4) 4 282.54 70.64 7.2634 0.002201
## Residuals 14 136.15 9.72
## Lack of fit 10 125.40 12.54 4.6660 0.075500
## Pure error 4 10.75 2.69
##
## Stationary point of response surface:
## x1 x2 x3 x4
## 0.8607107 -0.3307115 -0.8394866 -0.1161465
##
## Stationary point in original units:
## A R W L
## 12.916426 2.434015 1.040128 1.941927
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 3.258222 -1.198324 -3.807935 -4.651963
##
## $vectors
## [,1] [,2] [,3] [,4]
## x1 0.5177048 0.04099358 -0.7608371 -0.38913772
## x2 -0.4504231 0.58176202 -0.5056034 0.45059647
## x3 -0.4517232 0.37582195 0.1219894 -0.79988915
## x4 0.5701289 0.72015994 0.3880860 0.07557783Predict the response variable at the stationary point
#predict the response variable using this model at these parameters
max <- data.frame(x1 = 0.8607107, x2 = -0.3307115, x3 = -0.8394866, x4 = -0.1161465, block = "1") #use the coded values
predict(heli.rsm, newdata = max) ## 1
## 372.1719many second order have a significant p values , means that they contribute significantly to the model and that is canonical analysis is relevant
stationarity points seems to be close to optimized regionn but since the eigenvalues ARE of mixed sign so this stationary points is a saddle points (neither min or max ) .thus further analysis is required.
###Displaying a response surface
##The canonical analyses and the breakdown of their tables are important for understanding the behavior of a second order response surface ##However, effective graphs tend to be easier to interpret ##rsm includes functions for data visualization ##In the paper helicopter example, since there are four predictor variables, there are six pairs of predictors
##Visualize the behavior of the fitted surface around the stationary point:
par(mfrow = c(2, 3))
contour(heli.rsm, ~ x1 + x2 + x3 + x4,
image = TRUE, #fills spaces between the contour lines with color; color levels are consistent across plots
at = summary(heli.rsm)$canonical$xs) #location of stationary point in the coded units 
##direction for further experimentation
##In many first-order and some second-order cases, we find that the saddle point or stationary point is far from the optimal experimental region ##If this is the case then the most practical the next is to determine where in which direction to investigate further ##This is a facet of RSM called direction of steepest ascent ##The rsm summary table provides information about this concept and we can zoom into that using the steepest function:
steepest(CR1.rsm, #dataframe to access
dist = c(0,0.5,1)) #distance ## Path of steepest ascent from ridge analysis:## dist x1 x2 | Time Temp | yhat
## 1 0.0 0.000 0.000 | 85.000 175.000 | 82.814
## 2 0.5 0.407 0.291 | 87.035 176.455 | 83.352
## 3 1.0 0.814 0.581 | 89.070 177.905 | 83.890##yhat is the predicted value of the response variable in the model ##any set of distances can be factored in using the dist argument ##However, while the fitted values are displayed, remember that these are only predictions ##As the distance along the path increases, the predictions become less reliable ##The real use case of the path of steepest ascent is to determine the next set of experiments to get close to the optimal experimental conditions ##For a second-order model, the steepest function works, but it employs the ‘ridge analysis method’ which is a type of analog of the steepest ascent ##In this method, for a specified distance, d, a point at which the predicted response is a maximum among all predictor combinations at radius d can be determined ##It is sensible to use this method when the stationary point is somewhat far away
steepest(heli.rsm, #dataframe to access
dist = c(0,0.5,1)) #distance ## Path of steepest ascent from ridge analysis:## dist x1 x2 x3 x4 | A R W L | yhat
## 1 0.0 0.000 0.000 0.000 0.000 | 12.4000 2.52000 1.250 2.0000 | 372.800
## 2 0.5 -0.127 0.288 0.116 -0.371 | 12.3238 2.59488 1.279 1.8145 | 377.106
## 3 1.0 -0.351 0.538 0.312 -0.700 | 12.1894 2.65988 1.328 1.6500 | 382.675