Agenda
This lecture will practice running
- Plackett-Burman Design (PBD)
and
- Response surface method (RSM)
Before begin..
Let’s install packages.
install.packages("daewr")##
## The downloaded binary packages are in
## /var/folders/zh/85dfy3gx6rd659vpk8_r_ppr0000gn/T//RtmpJEQrxm/downloaded_packageslibrary(daewr)PBDes function creates the experimental design of the
PBD. randomize = T will generate a matrix with mixed orders.
daewr::PBDes(nfactors = 11, nruns = 12)## A B C D E F G H J K L
## 1 1 1 -1 1 1 1 -1 -1 -1 1 -1
## 2 -1 1 1 -1 1 1 1 -1 -1 -1 1
## 3 1 -1 1 1 -1 1 1 1 -1 -1 -1
## 4 -1 1 -1 1 1 -1 1 1 1 -1 -1
## 5 -1 -1 1 -1 1 1 -1 1 1 1 -1
## 6 -1 -1 -1 1 -1 1 1 -1 1 1 1
## 7 1 -1 -1 -1 1 -1 1 1 -1 1 1
## 8 1 1 -1 -1 -1 1 -1 1 1 -1 1
## 9 1 1 1 -1 -1 -1 1 -1 1 1 -1
## 10 -1 1 1 1 -1 -1 -1 1 -1 1 1
## 11 1 -1 1 1 1 -1 -1 -1 1 -1 1
## 12 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1# Error examples
# daewr::PBDes(nfactors = 12, nruns = 12)
daewr::PBDes(nfactors = 10, nruns = 12)## A B C D E F G H J K
## 1 1 1 -1 1 1 1 -1 -1 -1 1
## 2 -1 1 1 -1 1 1 1 -1 -1 -1
## 3 1 -1 1 1 -1 1 1 1 -1 -1
## 4 -1 1 -1 1 1 -1 1 1 1 -1
## 5 -1 -1 1 -1 1 1 -1 1 1 1
## 6 -1 -1 -1 1 -1 1 1 -1 1 1
## 7 1 -1 -1 -1 1 -1 1 1 -1 1
## 8 1 1 -1 -1 -1 1 -1 1 1 -1
## 9 1 1 1 -1 -1 -1 1 -1 1 1
## 10 -1 1 1 1 -1 -1 -1 1 -1 1
## 11 1 -1 1 1 1 -1 -1 -1 1 -1
## 12 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1#daewr::PBDes(nfactors = 15, nruns = 24)After setting your higher level and
lower level based on your intuition, you can run a
experiment with 12 experimental groups, and get the result out of the
data.
analyiss of PBD outcome
You can run a miultiple linear regression that is containing all the data!
library(tidyverse)
#<<<<<<< HEAD
download.file("https://raw.githubusercontent.com/minsiksudo/BTE3207_Advanced_Biostatistics/refs/heads/main/PBD_example_data.csv", "PBD_example_data.csv")
pbd_data <- read.csv("PBD_example_data.csv") lm(Biomass ~ Glucose + YE + KH2PO4 + K2HPO4 + MgSO4 + CaCl2 + FeSO4 + Fructose + NH4Cl + TM + Vitamin,
#<<<<<<< HEAD
data = pbd_data) %>% summary##
## Call:
## lm(formula = Biomass ~ Glucose + YE + KH2PO4 + K2HPO4 + MgSO4 +
## CaCl2 + FeSO4 + Fructose + NH4Cl + TM + Vitamin, data = pbd_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04 -0.02 0.00 0.02 0.04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.437500 0.007022 347.134 < 2e-16 ***
## Glucose -0.012500 0.007022 -1.780 0.100360
## YE 0.845833 0.007022 120.458 < 2e-16 ***
## KH2PO4 0.037500 0.007022 5.341 0.000176 ***
## K2HPO4 0.572500 0.007022 81.532 < 2e-16 ***
## MgSO4 0.032500 0.007022 4.628 0.000582 ***
## CaCl2 -0.054167 0.007022 -7.714 5.45e-06 ***
## FeSO4 0.055833 0.007022 7.951 4.00e-06 ***
## Fructose 0.080833 0.007022 11.512 7.67e-08 ***
## NH4Cl -0.017500 0.007022 -2.492 0.028315 *
## TM 0.015833 0.007022 2.255 0.043611 *
## Vitamin -0.019167 0.007022 -2.730 0.018280 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0344 on 12 degrees of freedom
## Multiple R-squared: 0.9994, Adjusted R-squared: 0.9989
## F-statistic: 1953 on 11 and 12 DF, p-value: < 2.2e-16This is your Plackett-Burman result.
RSM
install.packages("rsm")##
## The downloaded binary packages are in
## /var/folders/zh/85dfy3gx6rd659vpk8_r_ppr0000gn/T//RtmpJEQrxm/downloaded_packageslibrary(rsm)using ccd function, you can generate a experimental
design that can be used for RSM experiment and analysis
#2 factor experiment
ccd (2)## run.order std.order x1.as.is x2.as.is Block
## 1 1 7 0.000000 0.000000 1
## 2 2 2 1.000000 -1.000000 1
## 3 3 1 -1.000000 -1.000000 1
## 4 4 8 0.000000 0.000000 1
## 5 5 6 0.000000 0.000000 1
## 6 6 4 1.000000 1.000000 1
## 7 7 3 -1.000000 1.000000 1
## 8 8 5 0.000000 0.000000 1
## 9 1 6 0.000000 0.000000 2
## 10 2 7 0.000000 0.000000 2
## 11 3 5 0.000000 0.000000 2
## 12 4 4 0.000000 1.414214 2
## 13 5 2 1.414214 0.000000 2
## 14 6 3 0.000000 -1.414214 2
## 15 7 8 0.000000 0.000000 2
## 16 8 1 -1.414214 0.000000 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is#3 factor experiment
ccd (3)## run.order std.order x1.as.is x2.as.is x3.as.is Block
## 1 1 8 1.000000 1.000000 1.000000 1
## 2 2 2 1.000000 -1.000000 -1.000000 1
## 3 3 1 -1.000000 -1.000000 -1.000000 1
## 4 4 11 0.000000 0.000000 0.000000 1
## 5 5 5 -1.000000 -1.000000 1.000000 1
## 6 6 9 0.000000 0.000000 0.000000 1
## 7 7 10 0.000000 0.000000 0.000000 1
## 8 8 4 1.000000 1.000000 -1.000000 1
## 9 9 7 -1.000000 1.000000 1.000000 1
## 10 10 6 1.000000 -1.000000 1.000000 1
## 11 11 12 0.000000 0.000000 0.000000 1
## 12 12 3 -1.000000 1.000000 -1.000000 1
## 13 1 1 -1.825742 0.000000 0.000000 2
## 14 2 8 0.000000 0.000000 0.000000 2
## 15 3 3 0.000000 -1.825742 0.000000 2
## 16 4 9 0.000000 0.000000 0.000000 2
## 17 5 10 0.000000 0.000000 0.000000 2
## 18 6 6 0.000000 0.000000 1.825742 2
## 19 7 4 0.000000 1.825742 0.000000 2
## 20 8 7 0.000000 0.000000 0.000000 2
## 21 9 2 1.825742 0.000000 0.000000 2
## 22 10 5 0.000000 0.000000 -1.825742 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.is
## x3 ~ x3.as.isSOdes2 <- ccd (2, n0 = c(4,6), alpha = "rotatable", inscribed = F)
SOdes2## run.order std.order x1.as.is x2.as.is Block
## 1 1 8 0.000000 0.000000 1
## 2 2 5 0.000000 0.000000 1
## 3 3 6 0.000000 0.000000 1
## 4 4 4 1.000000 1.000000 1
## 5 5 1 -1.000000 -1.000000 1
## 6 6 2 1.000000 -1.000000 1
## 7 7 7 0.000000 0.000000 1
## 8 8 3 -1.000000 1.000000 1
## 9 1 8 0.000000 0.000000 2
## 10 2 1 -1.414214 0.000000 2
## 11 3 4 0.000000 1.414214 2
## 12 4 2 1.414214 0.000000 2
## 13 5 5 0.000000 0.000000 2
## 14 6 3 0.000000 -1.414214 2
## 15 7 10 0.000000 0.000000 2
## 16 8 9 0.000000 0.000000 2
## 17 9 6 0.000000 0.000000 2
## 18 10 7 0.000000 0.000000 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ x1.as.is
## x2 ~ x2.as.isWith coding argument, you can set the levels of your
data into actual values
SOdes2 <- ccd (3, n0 = c(4,6), alpha = "rotatable", inscribed = F,
coding = list (
x1 ~ (Glucose - 50)/10,
x2 ~ (NaNO3 - 5)/2,
x3 ~ (K2HPO4 - 2)/1)
)
SOdes2## run.order std.order Glucose NaNO3 K2HPO4 Block
## 1 1 3 40.00000 7.000000 1.0000000 1
## 2 2 7 40.00000 7.000000 3.0000000 1
## 3 3 8 60.00000 7.000000 3.0000000 1
## 4 4 11 50.00000 5.000000 2.0000000 1
## 5 5 1 40.00000 3.000000 1.0000000 1
## 6 6 5 40.00000 3.000000 3.0000000 1
## 7 7 12 50.00000 5.000000 2.0000000 1
## 8 8 4 60.00000 7.000000 1.0000000 1
## 9 9 2 60.00000 3.000000 1.0000000 1
## 10 10 10 50.00000 5.000000 2.0000000 1
## 11 11 6 60.00000 3.000000 3.0000000 1
## 12 12 9 50.00000 5.000000 2.0000000 1
## 13 1 2 66.81793 5.000000 2.0000000 2
## 14 2 4 50.00000 8.363586 2.0000000 2
## 15 3 10 50.00000 5.000000 2.0000000 2
## 16 4 6 50.00000 5.000000 3.6817928 2
## 17 5 9 50.00000 5.000000 2.0000000 2
## 18 6 1 33.18207 5.000000 2.0000000 2
## 19 7 12 50.00000 5.000000 2.0000000 2
## 20 8 3 50.00000 1.636414 2.0000000 2
## 21 9 11 50.00000 5.000000 2.0000000 2
## 22 10 8 50.00000 5.000000 2.0000000 2
## 23 11 5 50.00000 5.000000 0.3182072 2
## 24 12 7 50.00000 5.000000 2.0000000 2
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Glucose - 50)/10
## x2 ~ (NaNO3 - 5)/2
## x3 ~ (K2HPO4 - 2)/1Analysis of RSM
CR1.rsm <- rsm(Yield ~ Block + SO(Time, Temp), data = ChemReact)
summary(CR1.rsm)##
## Call:
## rsm(formula = Yield ~ Block + SO(Time, Temp), data = ChemReact)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3992e+03 9.0919e+01 -15.390 1.180e-06 ***
## BlockB2 -4.4575e+00 8.7226e-02 -51.103 2.877e-10 ***
## Time 8.2097e+00 7.0225e-01 11.691 7.576e-06 ***
## Temp 1.2759e+01 8.8554e-01 14.408 1.848e-06 ***
## Time:Temp 5.0000e-03 3.2637e-03 1.532 0.1694
## Time^2 -5.2342e-02 2.4025e-03 -21.786 1.083e-07 ***
## Temp^2 -3.7338e-02 2.4025e-03 -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(Time, Temp) 2 9.626 4.813 180.7341 9.450e-07
## TWI(Time, Temp) 1 0.062 0.062 2.3470 0.1694
## PQ(Time, Temp) 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:
## Time Temp
## 86.86148 176.67190
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] -0.03693211 -0.05274780
##
## $vectors
## [,1] [,2]
## Time -0.1601375 -0.9870947
## Temp -0.9870947 0.1601375Here, besides the interaction term, all the term is significant.
Anova here tests the “lack of fit” which tests the variance not explaned by model but by the error (from the repeated experiments). This term’s p-values should be greater than effect of others!
To draw contour, use contour function
contour(CR1.rsm, ~ Time + Temp, image = TRUE, at = summary(CR1.rsm)$canonical$xs)
To draw 3d plot, use persp function
persp(CR1.rsm,
Time ~ Temp, #plot xs
col=rainbow(50),
xlabs = c(expression("Time (min)", "Temperature (°C)")), # axis labels
theta=30,
phi=30, r = 120, d=120,
border = NULL,
ltheta = 0,
lphi = 0,
shade = 0.75, zlab="Yield", col.axis=37, font.lab=2, col.lab=35,
contour=("colors"))
Changing theta, phi, and r can rotate the graph
persp(CR1.rsm,
Time ~ Temp, #plot xs
col=rainbow(50),
xlabs = c(expression("Time (min)", "Temperature (°C)")), # axis labels
theta=50,
phi=3, r = 180, d=120,
border = NULL,
ltheta = 0,
lphi = 0,
shade = 0.75, zlab="Yield", col.axis=37, font.lab=2, col.lab=35,
contour=("colors"))
Bibliography
## Computing. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>. We have invested a lot of time and effort in creating R, please cite it when using it for data analysis. See also 'citation("pkgname")' for citing R packages.
## Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
## R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman
## J, reikoch, Beasley W, O'Connor B, Warnes GR, Quinn M, Kamvar ZN, Gao C (2024). yaml: Methods to Convert R Data to YAML and Back_. R package version 2.3.10, <https://CRAN.R-project.org/package=yaml>. ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see 'help("citation")'.
## version 0.4.4, <https://CRAN.R-project.org/package=reactable>.