可以根据 StatEase 中的旧结果重新计算 R 中的响应曲面设计实验吗？
Possible to re-calculate a response surface design experiment in R from old results in StatEase?
我正在尝试复制一些几年前使用旧的 StatEase 6 软件生成的结果。我做过约束设计，响应面模型，D-最优设计
参数 1 的值介于 12 到 10.5 之间。参数 2 的范围从 18 到 3。我想优化我的反应产率。
我只想要可以生成相关 3Dplot 和相应显着性统计分析的最终方程式。
当时table软件生成的（以及实验室生成的相应结果）如下table所示：
|区块 1|
Parameter1
Parameter2
Yield
10.5
3.0
5.8
10.5
18.0
6.1
11.25
10.5
7.9
11.25
10.5
8.7
11.25
10.5
11.8
12.00
3.0
4.8
12.00
18.0
8.7
|区块 2|
Parameter1
Parameter2
Yield
10.19
10.5
8.5
11.25
0.0
0.21
11.25
10.5
9.4
11.25
10.5
8.5
11.25
10.5
10.1
12.31
10.5
8.9
12.31
21.1
8.5
即使参考了这个不错的教程，我也无法插入这些值并进行分析：https://bookdown.org/gerhard_krennrich/doe_and_optimization/doe2.html#optimal-designs
非常感谢您的帮助。
我花了一些时间试图更多地了解我所看到的内容。如果这不是您要找的，请告诉我。这应该可以解决您对响应曲面模型的请求。
我也包括了 RSM 的假设 - 假设这是一张 LP。
rsm 是这里使用的主包；它代表响应面模型（方法？）。软件包编写者编写了 响应面分析。
区块 1 的模型较差（可能不是线性的）；街区 2 看起来好多了。 block 1 的所有模型都以“fit”结尾；块 2 以“fit2.”结尾。
library(tidyverse)
library(rsm)
library(car)
library(plotly)
#----------------- first block ---------------------
b1 <- data.frame(
Parameter1 = c(10.5, 10.5, 11.25, 11.25, 11.25, 12, 12),
Parameter2 = c(3, 18, 10.5, 10.5, 10.5, 3, 18),
Yield = c(5.8, 6.1, 7.9, 8.7, 11.8, 4.8, 8.7)
)
names(b1) <- c("x1", "x2", "y")
# [1] "x1" "x2" "y"
#----------- second-order model (SO) ---------------
rsm.fit <- rsm(y ~ SO(x1, x2), data = b1)
summary(rsm.fit)
这在 R 中生成了一条警告，指出“不能使用 'rsm' 方法。”
它 returns 是一个 LM（线性模型）对象。你得到很多相同的信息，所以我把它留在这里。
与 statease 一样，它执行几种不同的模型 - 在摘要中进行了预测。这是摘要的一部分：
# Coefficients: (1 not defined because of singularities)
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -680.3533 354.0967 -1.921 0.195
# FO(x1, x2)x1 123.5200 62.9764 1.961 0.189
# FO(x1, x2)x2 -1.6600 2.0645 -0.804 0.506
# TWI(x1, x2) 0.1600 0.1831 0.874 0.474
# PQ(x1, x2)x1^2 -5.5407 2.7970 -1.981 0.186
# PQ(x1, x2)x2^2 NA NA NA NA
#
# Residual standard error: 2.06 on 2 degrees of freedom
# Multiple R-squared: 0.7461, Adjusted R-squared: 0.2384
# F-statistic: 1.469 on 4 and 2 DF, p-value: 0.4433
我运行其他几个模型，先从一个一阶模型开始。这返回了一个 RSM。
#--------------- first order (FO) model ------------------
# first order model
rsm2.fit <- rsm(y ~ FO(x1, x2), data = b1)
summary(rsm2.fit)
#----- first-order, 2-way interaction model (TWI) --------
rsm3.fit <- rsm(y ~ FO(x1, x2) + TWI(x1, x2), data = b1)
summary(rsm3.fit)
#----- first-order, 2-way NO interaction model -----------
rsm4.fit <- rsm(y ~ FO(x1, x2) + PQ(x1, x2), data = b1)
summary(rsm4.fit)
#-------------------- list the models --------------------
# create a list of the models to make assumption and plotting easier
fit <- list(rsm.fit = rsm.fit,
rsm2.fit = rsm2.fit,
rsm3.fit = rsm3.fit,
rsm4.fit = rsm4.fit)
#--------------------- assumptions -----------------------
# plot residuals
par(mfrow = c(2,2)) # four plots per page
lapply(fit, plot, which = 1:4)
#qqplot
qqPlot(rsm.fit)
qqPlot(rsm2.fit)
qqPlot(rsm3.fit)
qqPlot(rsm4.fit)
lapply(fit, cooks.distance)
#-------------------- visualize RSM -----------------------
# plot model contours
lapply(1:length(fit), function(x) {
contour(fit[[x]], ~ x1 + x2, image = T,
main = paste0("Block 1: Model ", x))
})
# create base R persp and plot_ly 3D plots
# empty lists to store plots of each model
p.fit <- vector(mode = "list", length = length(fit)) # persp plots
plt.fit <- vector(mode = "list", length = length(fit)) # plotly plots
for(i in 1:length(fit)){
p <- persp(fit[[i]],
x2~x1, zlab = "y",
main = paste("Block 1: Model ", i))
p.fit[[i]] <- p # so you can still view persp plots
x1 = p$`x1 ~ x2`$x # for plotly
x2 = p$`x1 ~ x2`$y
y = p$`x1 ~ x2`$z
plt.fit[[i]] <- plotly_build(
plot_ly(x = ~x1, y = ~x2, z = ~y) %>%
add_surface(contours = list(
z = list(usecolormap = T, # create contours
show = T,
highlightcolor = "#ff0000",
# highlight contour in red when hovering
project = list(z = T))),
hovertemplate = paste0('x1: %{x}<br>',
'x2: %{y}<br>',
'y: %{z}<extra></extra>')) %>%
layout(title = list(text = paste0("Block 1: Model ", i),
y = .95)) # pull title down a bit
)
}
plt.fit[[1]]
plt.fit[[2]]
plt.fit[[3]]
plt.fit[[4]]
轮廓和视角
如果你以前没有用过 Plotly，它的交互性很强。
我对第二个块又做了一遍——除了数据之外都是一样的。
#--------------------- 2nd Block ----------------------
b2 <- data.frame(
Parameter1 = c(10.19, 11.25, 11.25, 11.25, 11.25, 12.31, 12.31),
Parameter2 = c(10.5, 0, 10.5, 10.5, 10.5, 10.5, 21.1),
Yield = c(8.5, 0.21, 9.4, 8.5, 10.1, 8.9, 8.5)
)
names(b2) <- c("x1", "x2", "y")
#---------------------- SO model ----------------------
rsm.fit2 <- rsm(y ~ SO(x1, x2), data = b2)
summary(rsm.fit2)
块 2 上的二阶模型没有错误。
#---------------------- FO Model ----------------------
rsm2.fit2 <- rsm(y ~ FO(x1, x2), data = b2)
summary(rsm2.fit2)
#--------------------- FO TWI model -------------------
rsm3.fit2 <- rsm(y ~ FO(x1, x2) + TWI(x1, x2), data = b2)
summary(rsm3.fit2)
#----------------- FO NO interaction model--------------
rsm4.fit2 <- rsm(y ~ FO(x1, x2) + PQ(x1, x2), data = b2)
summary(rsm4.fit2)
#-------------------- list the models ------------------
# create a list of the models to make assumption and plotting easier
fit2 <- list(rsm.fit2 = rsm.fit2,
rsm2.fit2 = rsm2.fit2,
rsm3.fit2 = rsm3.fit2,
rsm4.fit2 = rsm4.fit2)
#-------------------- assumptions ----------------------
# plot residuals
lapply(fit2, plot, which = 1:4)
#qqplot
qqPlot(rsm.fit2)
qqPlot(rsm2.fit2)
qqPlot(rsm3.fit2)
qqPlot(rsm4.fit2)
# cooks distance
lapply(fit2, cooks.distance)
#-------------------- visualize RSM -----------------------
# plot model contours
lapply(1:length(fit2), function(x) {
contour(fit2[[x]], ~ x1 + x2, image = T,
main = paste0("Block 2: Model ", x))
})
# base R contours (scroll through plot pane to view)
# create base R persp and plot_ly 3D plots
# empty lists to store plots of each model
p.fit2 <- vector(mode = "list", length = length(fit2)) # persp plots
plt.fit2 <- vector(mode = "list", length = length(fit2)) # plotly plots
for(i in 1:length(fit2)){
p <- persp(fit2[[i]],
x2~x1, zlab = "y",
main = paste("Block 2: Model ", i))
p.fit2[[i]] <- p # so you can still view persp plots
x1 = p$`x1 ~ x2`$x # for plotly
x2 = p$`x1 ~ x2`$y
y = p$`x1 ~ x2`$z
plt.fit2[[i]] <- plotly_build(
plot_ly(x = ~x1, y = ~x2, z = ~y) %>%
add_surface(contours = list(
z = list(usecolormap = T, # create contours
show = T,
highlightcolor = "#ff0000",
# highlight contour in red when hovering
project = list(z = T))),
hovertemplate = paste0('x1: %{x}<br>',
'x2: %{y}<br>',
'y: %{z}<extra></extra>')) %>%
layout(title = list(text = paste0("Block 2: Model ", i),
y = .95)) # pull title down a bit
)
}
plt.fit2[[1]]
plt.fit2[[2]]
plt.fit2[[3]]
plt.fit2[[4]]
# reset plots per page
par(mfrow = c(1,1))
轮廓和视角
情节观点
我正在尝试复制一些几年前使用旧的 StatEase 6 软件生成的结果。我做过约束设计，响应面模型，D-最优设计
参数 1 的值介于 12 到 10.5 之间。参数 2 的范围从 18 到 3。我想优化我的反应产率。 我只想要可以生成相关 3Dplot 和相应显着性统计分析的最终方程式。
当时table软件生成的（以及实验室生成的相应结果）如下table所示：
|区块 1|
| Parameter1 | Parameter2 | Yield |
|---|---|---|
| 10.5 | 3.0 | 5.8 |
| 10.5 | 18.0 | 6.1 |
| 11.25 | 10.5 | 7.9 |
| 11.25 | 10.5 | 8.7 |
| 11.25 | 10.5 | 11.8 |
| 12.00 | 3.0 | 4.8 |
| 12.00 | 18.0 | 8.7 |
|区块 2|
| Parameter1 | Parameter2 | Yield |
|---|---|---|
| 10.19 | 10.5 | 8.5 |
| 11.25 | 0.0 | 0.21 |
| 11.25 | 10.5 | 9.4 |
| 11.25 | 10.5 | 8.5 |
| 11.25 | 10.5 | 10.1 |
| 12.31 | 10.5 | 8.9 |
| 12.31 | 21.1 | 8.5 |
即使参考了这个不错的教程，我也无法插入这些值并进行分析：https://bookdown.org/gerhard_krennrich/doe_and_optimization/doe2.html#optimal-designs
非常感谢您的帮助。
我花了一些时间试图更多地了解我所看到的内容。如果这不是您要找的，请告诉我。这应该可以解决您对响应曲面模型的请求。
我也包括了 RSM 的假设 - 假设这是一张 LP。
rsm 是这里使用的主包；它代表响应面模型（方法？）。软件包编写者编写了 响应面分析。
区块 1 的模型较差（可能不是线性的）；街区 2 看起来好多了。 block 1 的所有模型都以“fit”结尾；块 2 以“fit2.”结尾。
library(tidyverse)
library(rsm)
library(car)
library(plotly)
#----------------- first block ---------------------
b1 <- data.frame(
Parameter1 = c(10.5, 10.5, 11.25, 11.25, 11.25, 12, 12),
Parameter2 = c(3, 18, 10.5, 10.5, 10.5, 3, 18),
Yield = c(5.8, 6.1, 7.9, 8.7, 11.8, 4.8, 8.7)
)
names(b1) <- c("x1", "x2", "y")
# [1] "x1" "x2" "y"
#----------- second-order model (SO) ---------------
rsm.fit <- rsm(y ~ SO(x1, x2), data = b1)
summary(rsm.fit)
这在 R 中生成了一条警告，指出“不能使用 'rsm' 方法。” 它 returns 是一个 LM（线性模型）对象。你得到很多相同的信息，所以我把它留在这里。
与 statease 一样，它执行几种不同的模型 - 在摘要中进行了预测。这是摘要的一部分：
# Coefficients: (1 not defined because of singularities)
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -680.3533 354.0967 -1.921 0.195
# FO(x1, x2)x1 123.5200 62.9764 1.961 0.189
# FO(x1, x2)x2 -1.6600 2.0645 -0.804 0.506
# TWI(x1, x2) 0.1600 0.1831 0.874 0.474
# PQ(x1, x2)x1^2 -5.5407 2.7970 -1.981 0.186
# PQ(x1, x2)x2^2 NA NA NA NA
#
# Residual standard error: 2.06 on 2 degrees of freedom
# Multiple R-squared: 0.7461, Adjusted R-squared: 0.2384
# F-statistic: 1.469 on 4 and 2 DF, p-value: 0.4433
我运行其他几个模型，先从一个一阶模型开始。这返回了一个 RSM。
#--------------- first order (FO) model ------------------
# first order model
rsm2.fit <- rsm(y ~ FO(x1, x2), data = b1)
summary(rsm2.fit)
#----- first-order, 2-way interaction model (TWI) --------
rsm3.fit <- rsm(y ~ FO(x1, x2) + TWI(x1, x2), data = b1)
summary(rsm3.fit)
#----- first-order, 2-way NO interaction model -----------
rsm4.fit <- rsm(y ~ FO(x1, x2) + PQ(x1, x2), data = b1)
summary(rsm4.fit)
#-------------------- list the models --------------------
# create a list of the models to make assumption and plotting easier
fit <- list(rsm.fit = rsm.fit,
rsm2.fit = rsm2.fit,
rsm3.fit = rsm3.fit,
rsm4.fit = rsm4.fit)
#--------------------- assumptions -----------------------
# plot residuals
par(mfrow = c(2,2)) # four plots per page
lapply(fit, plot, which = 1:4)
#qqplot
qqPlot(rsm.fit)
qqPlot(rsm2.fit)
qqPlot(rsm3.fit)
qqPlot(rsm4.fit)
lapply(fit, cooks.distance)
#-------------------- visualize RSM -----------------------
# plot model contours
lapply(1:length(fit), function(x) {
contour(fit[[x]], ~ x1 + x2, image = T,
main = paste0("Block 1: Model ", x))
})
# create base R persp and plot_ly 3D plots
# empty lists to store plots of each model
p.fit <- vector(mode = "list", length = length(fit)) # persp plots
plt.fit <- vector(mode = "list", length = length(fit)) # plotly plots
for(i in 1:length(fit)){
p <- persp(fit[[i]],
x2~x1, zlab = "y",
main = paste("Block 1: Model ", i))
p.fit[[i]] <- p # so you can still view persp plots
x1 = p$`x1 ~ x2`$x # for plotly
x2 = p$`x1 ~ x2`$y
y = p$`x1 ~ x2`$z
plt.fit[[i]] <- plotly_build(
plot_ly(x = ~x1, y = ~x2, z = ~y) %>%
add_surface(contours = list(
z = list(usecolormap = T, # create contours
show = T,
highlightcolor = "#ff0000",
# highlight contour in red when hovering
project = list(z = T))),
hovertemplate = paste0('x1: %{x}<br>',
'x2: %{y}<br>',
'y: %{z}<extra></extra>')) %>%
layout(title = list(text = paste0("Block 1: Model ", i),
y = .95)) # pull title down a bit
)
}
plt.fit[[1]]
plt.fit[[2]]
plt.fit[[3]]
plt.fit[[4]]
轮廓和视角
如果你以前没有用过 Plotly，它的交互性很强。
我对第二个块又做了一遍——除了数据之外都是一样的。
#--------------------- 2nd Block ----------------------
b2 <- data.frame(
Parameter1 = c(10.19, 11.25, 11.25, 11.25, 11.25, 12.31, 12.31),
Parameter2 = c(10.5, 0, 10.5, 10.5, 10.5, 10.5, 21.1),
Yield = c(8.5, 0.21, 9.4, 8.5, 10.1, 8.9, 8.5)
)
names(b2) <- c("x1", "x2", "y")
#---------------------- SO model ----------------------
rsm.fit2 <- rsm(y ~ SO(x1, x2), data = b2)
summary(rsm.fit2)
块 2 上的二阶模型没有错误。
#---------------------- FO Model ----------------------
rsm2.fit2 <- rsm(y ~ FO(x1, x2), data = b2)
summary(rsm2.fit2)
#--------------------- FO TWI model -------------------
rsm3.fit2 <- rsm(y ~ FO(x1, x2) + TWI(x1, x2), data = b2)
summary(rsm3.fit2)
#----------------- FO NO interaction model--------------
rsm4.fit2 <- rsm(y ~ FO(x1, x2) + PQ(x1, x2), data = b2)
summary(rsm4.fit2)
#-------------------- list the models ------------------
# create a list of the models to make assumption and plotting easier
fit2 <- list(rsm.fit2 = rsm.fit2,
rsm2.fit2 = rsm2.fit2,
rsm3.fit2 = rsm3.fit2,
rsm4.fit2 = rsm4.fit2)
#-------------------- assumptions ----------------------
# plot residuals
lapply(fit2, plot, which = 1:4)
#qqplot
qqPlot(rsm.fit2)
qqPlot(rsm2.fit2)
qqPlot(rsm3.fit2)
qqPlot(rsm4.fit2)
# cooks distance
lapply(fit2, cooks.distance)
#-------------------- visualize RSM -----------------------
# plot model contours
lapply(1:length(fit2), function(x) {
contour(fit2[[x]], ~ x1 + x2, image = T,
main = paste0("Block 2: Model ", x))
})
# base R contours (scroll through plot pane to view)
# create base R persp and plot_ly 3D plots
# empty lists to store plots of each model
p.fit2 <- vector(mode = "list", length = length(fit2)) # persp plots
plt.fit2 <- vector(mode = "list", length = length(fit2)) # plotly plots
for(i in 1:length(fit2)){
p <- persp(fit2[[i]],
x2~x1, zlab = "y",
main = paste("Block 2: Model ", i))
p.fit2[[i]] <- p # so you can still view persp plots
x1 = p$`x1 ~ x2`$x # for plotly
x2 = p$`x1 ~ x2`$y
y = p$`x1 ~ x2`$z
plt.fit2[[i]] <- plotly_build(
plot_ly(x = ~x1, y = ~x2, z = ~y) %>%
add_surface(contours = list(
z = list(usecolormap = T, # create contours
show = T,
highlightcolor = "#ff0000",
# highlight contour in red when hovering
project = list(z = T))),
hovertemplate = paste0('x1: %{x}<br>',
'x2: %{y}<br>',
'y: %{z}<extra></extra>')) %>%
layout(title = list(text = paste0("Block 2: Model ", i),
y = .95)) # pull title down a bit
)
}
plt.fit2[[1]]
plt.fit2[[2]]
plt.fit2[[3]]
plt.fit2[[4]]
# reset plots per page
par(mfrow = c(1,1))
轮廓和视角
情节观点