阶乘 - 四向方差分析 - 如何找到统计上有效的组合
Factorial - Four Way ANOVA - How to find a statistically effective combination
我必须分析实验数据集以找到最有效的分子生物学反应组合。
实验有四个因素：温度、转速、时间、催化activity。我正在测量反应效率 (EE)。我怎样才能找到最高效率 (EE) 的四个因素的有效组合？
- 没有重复测量。所有数据均为独立实验数据
据我了解 - EE 是参数数据，因素是分类数据（固定组合）。
我是否必须进行 Fourway ANOVA？
如果是这样，这个模型对分析是否正确
library(lsmeans)
lm(EE ~ Temperature + RPM + Time+ Catalytic +
Temperature:RPM +
Temperature:Time +
Temperature:Catalytic +
RPM:Time+
RPM+Catalytic+
Time+Catalytic+
Temperature:RPM:Time +
Temperature:RPM:Catalytic+
Temperature:Time:Catalytic+
RPM:Time:Catalytic+
Temperature:RPM:Time:Catalytic, "data")
然后，我怎样才能得到每个成对比较的重要值？
这里是示例数据集。
> dput(df)
structure(list(TEMPERATURE = c(40, 40, 40, 40, 40, 40, 40, 40,
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40,
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40,
42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5,
42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5,
42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5,
42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 45, 45, 45, 45, 45,
45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
45, 45, 45), RPM = c(150, 150, 150, 150, 150, 150, 150, 150,
150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 200,
200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 150, 150, 150, 150, 150, 150, 150,
150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150,
200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 200, 150, 150, 150, 150, 150, 150,
150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150,
150, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 200, 200), TIME = c(24, 24, 24,
24, 24, 48, 48, 48, 48, 48, 72, 72, 72, 72, 72, 96, 96, 96, 96,
96, 24, 24, 24, 24, 24, 48, 48, 48, 48, 48, 72, 72, 72, 72, 72,
96, 96, 96, 96, 96, 24, 24, 24, 24, 24, 48, 48, 48, 48, 48, 72,
72, 72, 72, 72, 96, 96, 96, 96, 96, 24, 24, 24, 24, 24, 48, 48,
48, 48, 48, 72, 72, 72, 72, 72, 96, 96, 96, 96, 96, 24, 24, 24,
24, 24, 48, 48, 48, 48, 48, 72, 72, 72, 72, 72, 96, 96, 96, 96,
96, 24, 24, 24, 24, 24, 48, 48, 48, 48, 48, 72, 72, 72, 72, 72,
96, 96, 96, 96, 96), CAT = c(4, 6, 8, 10, 12, 4, 6, 8, 10, 12,
4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10,
12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6,
8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12,
4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10,
12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6,
8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12),
EE = c(50, 53, 54, 57, 59, 53, 56, 59, 61, 64, 57, 58, 60,
62, 63, 56, 54, 52, 55, 55, 44, 48, 50, 50, 54, 49, 52, 56,
57, 56, 52, 56, 57, 58, 66, 46, 48, 48, 52, 49, 53, 57, 59,
62, 64, 54, 58, 60, 64, 66, 55, 59, 61, 63, 65, 54, 59, 64,
65, 67, 49, 51, 53, 54, 59, 50, 54, 63, 64, 64, 52, 56, 56,
59, 57, 52, 55, 58, 60, 63, 52, 56, 58, 61, 63, 54, 55, 58,
63, 63, 56, 58, 62, 62, 65, 57, 59, 62, 63, 66, 42, 42, 51,
54, 56, 46, 50, 52, 56, 58, 48, 51, 54, 55, 57, 48, 53, 56,
57, 61)), class = c("spec_tbl_df", "tbl_df", "tbl", "data.frame"
), row.names = c(NA, -120L), spec = structure(list(cols = list(
TEMPERATURE = structure(list(), class = c("collector_double",
"collector")), RPM = structure(list(), class = c("collector_double",
"collector")), TIME = structure(list(), class = c("collector_double",
"collector")), CAT = structure(list(), class = c("collector_double",
"collector")), EE = structure(list(), class = c("collector_double",
"collector"))), default = structure(list(), class = c("collector_guess",
"collector")), skip = 1), class = "col_spec"))
对于筛选 DOE，您收集的数据比需要的多。
这是一个起点，欢迎补充。
我会模拟你所有因素的线性组合：
model <-lm(EE ~ TEMPERATURE + RPM + TIME +CAT , data=df)
summary(model)
Call:
lm(formula = EE ~ TEMPERATURE + RPM + TIME + CAT, data = df)
Residuals:
Min 1Q Median 3Q Max
-10.1850 -1.5742 0.3383 1.7767 9.7033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.50833 6.47932 7.950 1.42e-12 ***
TEMPERATURE 0.27000 0.14245 1.895 0.06056 .
RPM -0.10533 0.01163 -9.056 4.10e-15 ***
TIME 0.03639 0.01084 3.358 0.00107 **
CAT 1.20417 0.10281 11.713 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.185 on 115 degrees of freedom
Multiple R-squared: 0.6706, Adjusted R-squared: 0.6591
F-statistic: 58.52 on 4 and 115 DF, p-value: < 2.2e-16
查看斜率估计的标志。虽然这个简单的模型假设线性关系，但斜率的物理意义是将低值的平均值与高值的平均值进行比较。
例如，温度项的斜率是正的。这意味着当温度从低值 (40C) 升高到高值 (45C) 时，效率会提高。
温度、时间和 CAT 作为正斜率，我会取可用的最大值。
RPM 有一个负斜率，所以我会选择可用的最低值。
因此我的实验预测我的预测将产生 EE=66，
虽然实验的最高结果是：
df[which.max(df$EE),]
# A tibble: 1 x 5
TEMPERATURE RPM TIME CAT EE
<dbl> <dbl> <dbl> <dbl> <dbl>
1 42.5 150 96 12 67
现在您可以通过查看此模型的结果来研究非线性关系：
通话：
lm(formula = EE ~ TEMPERATURE * RPM * TIME * CAT, data = df)
此处交互项的斜率比线性项小几个数量级。这可能会产生误导，因为变量没有标准化。
祝你好运。
尝试这样的事情：
library(rsm)
mod = rsm(EE ~ SO(Temperature, RPM, Time, Catalytic), data = data)
summary(mod)
这将拟合二阶曲面（模型方程包括所有预测变量、双向交互作用和正方形）。摘要显示固定点和相关统计数据。如果所有的特征值都是负的，那么它就是一个峰。否则你有某种鞍点。
此模型与 OP 中的模型不同，因为它没有三向或四向交互作用，但包括预测变量的平方，这对于拟合二阶响应曲面非常重要。
更多详细信息
我不得不稍微修改一下以说明 R 区分大小写的事实！
> mod = rsm(EE ~ SO(TEMPERATURE, RPM, TIME, CAT), data = df)
这里有一个问题，RPM 只有两个值，所以我们无法估计纯二次效应。所以只有一个 NA 系数，这会扰乱驻点的计算。但是，我们仍然可以绘制拟合曲面（尽管有一些警告消息）
> par(mfrow = c(2,3))
> contour(mod, ~TEMPERATURE+RPM+TIME+CAT)
看起来我们最好使用较大的 CAT 和较低的 RPM（参见该图），所以再看一遍：
> par(mfrow=c(1,1))
> contour(mod, ~ TEMPERATURE + TIME, at = list(CAT = 12, RPM = 150))
所以在视觉上，我们似乎在温度 43.5、时间 65、催化剂 12 和 rpm 150 左右获得最佳响应。
如果您坚持将这些建模为因子，可以做到，但您需要将所有预测变量转换为因子。这是一个常见的错误；您可以设计一个仅包含定量变量的几个不同值的实验，但 R 不会读懂您的想法并认为它是一个因素；你必须把它转换成一个。在下文中，我选择了最多具有 2 种交互作用的模型。
> facmod = lm(EE ~ (factor(TEMPERATURE) + factor(TIME) + factor(RPM) + factor(CAT))^2, data = df)
> library(emmeans)
> emmip(facmod, TIME ~ TEMPERATURE | CAT*RPM)
最高拟合响应是在催化剂 12、RPM 150、温度 42.5 或更高以及时间 96 时。很明显，150 RPM 更好（左右比较），高 CAT 更好（比较面板）垂直）。这些是不同的模型和有些不同的结果。我更喜欢 rsm 方法，因为它更系统。
我必须分析实验数据集以找到最有效的分子生物学反应组合。
实验有四个因素：温度、转速、时间、催化activity。我正在测量反应效率 (EE)。我怎样才能找到最高效率 (EE) 的四个因素的有效组合？
- 没有重复测量。所有数据均为独立实验数据
据我了解 - EE 是参数数据，因素是分类数据（固定组合）。 我是否必须进行 Fourway ANOVA？
如果是这样，这个模型对分析是否正确
library(lsmeans)
lm(EE ~ Temperature + RPM + Time+ Catalytic +
Temperature:RPM +
Temperature:Time +
Temperature:Catalytic +
RPM:Time+
RPM+Catalytic+
Time+Catalytic+
Temperature:RPM:Time +
Temperature:RPM:Catalytic+
Temperature:Time:Catalytic+
RPM:Time:Catalytic+
Temperature:RPM:Time:Catalytic, "data")
然后，我怎样才能得到每个成对比较的重要值？
这里是示例数据集。
> dput(df)
structure(list(TEMPERATURE = c(40, 40, 40, 40, 40, 40, 40, 40,
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40,
40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40,
42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5,
42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5,
42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5,
42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 42.5, 45, 45, 45, 45, 45,
45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
45, 45, 45), RPM = c(150, 150, 150, 150, 150, 150, 150, 150,
150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 200,
200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 150, 150, 150, 150, 150, 150, 150,
150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150,
200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 200, 150, 150, 150, 150, 150, 150,
150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150,
150, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 200, 200), TIME = c(24, 24, 24,
24, 24, 48, 48, 48, 48, 48, 72, 72, 72, 72, 72, 96, 96, 96, 96,
96, 24, 24, 24, 24, 24, 48, 48, 48, 48, 48, 72, 72, 72, 72, 72,
96, 96, 96, 96, 96, 24, 24, 24, 24, 24, 48, 48, 48, 48, 48, 72,
72, 72, 72, 72, 96, 96, 96, 96, 96, 24, 24, 24, 24, 24, 48, 48,
48, 48, 48, 72, 72, 72, 72, 72, 96, 96, 96, 96, 96, 24, 24, 24,
24, 24, 48, 48, 48, 48, 48, 72, 72, 72, 72, 72, 96, 96, 96, 96,
96, 24, 24, 24, 24, 24, 48, 48, 48, 48, 48, 72, 72, 72, 72, 72,
96, 96, 96, 96, 96), CAT = c(4, 6, 8, 10, 12, 4, 6, 8, 10, 12,
4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10,
12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6,
8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12,
4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10,
12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6,
8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12, 4, 6, 8, 10, 12),
EE = c(50, 53, 54, 57, 59, 53, 56, 59, 61, 64, 57, 58, 60,
62, 63, 56, 54, 52, 55, 55, 44, 48, 50, 50, 54, 49, 52, 56,
57, 56, 52, 56, 57, 58, 66, 46, 48, 48, 52, 49, 53, 57, 59,
62, 64, 54, 58, 60, 64, 66, 55, 59, 61, 63, 65, 54, 59, 64,
65, 67, 49, 51, 53, 54, 59, 50, 54, 63, 64, 64, 52, 56, 56,
59, 57, 52, 55, 58, 60, 63, 52, 56, 58, 61, 63, 54, 55, 58,
63, 63, 56, 58, 62, 62, 65, 57, 59, 62, 63, 66, 42, 42, 51,
54, 56, 46, 50, 52, 56, 58, 48, 51, 54, 55, 57, 48, 53, 56,
57, 61)), class = c("spec_tbl_df", "tbl_df", "tbl", "data.frame"
), row.names = c(NA, -120L), spec = structure(list(cols = list(
TEMPERATURE = structure(list(), class = c("collector_double",
"collector")), RPM = structure(list(), class = c("collector_double",
"collector")), TIME = structure(list(), class = c("collector_double",
"collector")), CAT = structure(list(), class = c("collector_double",
"collector")), EE = structure(list(), class = c("collector_double",
"collector"))), default = structure(list(), class = c("collector_guess",
"collector")), skip = 1), class = "col_spec"))
对于筛选 DOE，您收集的数据比需要的多。 这是一个起点，欢迎补充。
我会模拟你所有因素的线性组合：
model <-lm(EE ~ TEMPERATURE + RPM + TIME +CAT , data=df)
summary(model)
Call:
lm(formula = EE ~ TEMPERATURE + RPM + TIME + CAT, data = df)
Residuals:
Min 1Q Median 3Q Max
-10.1850 -1.5742 0.3383 1.7767 9.7033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.50833 6.47932 7.950 1.42e-12 ***
TEMPERATURE 0.27000 0.14245 1.895 0.06056 .
RPM -0.10533 0.01163 -9.056 4.10e-15 ***
TIME 0.03639 0.01084 3.358 0.00107 **
CAT 1.20417 0.10281 11.713 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.185 on 115 degrees of freedom
Multiple R-squared: 0.6706, Adjusted R-squared: 0.6591
F-statistic: 58.52 on 4 and 115 DF, p-value: < 2.2e-16
查看斜率估计的标志。虽然这个简单的模型假设线性关系，但斜率的物理意义是将低值的平均值与高值的平均值进行比较。 例如，温度项的斜率是正的。这意味着当温度从低值 (40C) 升高到高值 (45C) 时，效率会提高。
温度、时间和 CAT 作为正斜率，我会取可用的最大值。
RPM 有一个负斜率，所以我会选择可用的最低值。
因此我的实验预测我的预测将产生 EE=66， 虽然实验的最高结果是：
df[which.max(df$EE),]
# A tibble: 1 x 5
TEMPERATURE RPM TIME CAT EE
<dbl> <dbl> <dbl> <dbl> <dbl>
1 42.5 150 96 12 67
现在您可以通过查看此模型的结果来研究非线性关系： 通话：
lm(formula = EE ~ TEMPERATURE * RPM * TIME * CAT, data = df)
此处交互项的斜率比线性项小几个数量级。这可能会产生误导，因为变量没有标准化。
祝你好运。
尝试这样的事情：
library(rsm)
mod = rsm(EE ~ SO(Temperature, RPM, Time, Catalytic), data = data)
summary(mod)
这将拟合二阶曲面（模型方程包括所有预测变量、双向交互作用和正方形）。摘要显示固定点和相关统计数据。如果所有的特征值都是负的，那么它就是一个峰。否则你有某种鞍点。
此模型与 OP 中的模型不同，因为它没有三向或四向交互作用，但包括预测变量的平方，这对于拟合二阶响应曲面非常重要。
更多详细信息
我不得不稍微修改一下以说明 R 区分大小写的事实！
> mod = rsm(EE ~ SO(TEMPERATURE, RPM, TIME, CAT), data = df)
这里有一个问题，RPM 只有两个值，所以我们无法估计纯二次效应。所以只有一个 NA 系数，这会扰乱驻点的计算。但是，我们仍然可以绘制拟合曲面（尽管有一些警告消息）
> par(mfrow = c(2,3))
> contour(mod, ~TEMPERATURE+RPM+TIME+CAT)
看起来我们最好使用较大的 CAT 和较低的 RPM（参见该图），所以再看一遍：
> par(mfrow=c(1,1))
> contour(mod, ~ TEMPERATURE + TIME, at = list(CAT = 12, RPM = 150))
所以在视觉上，我们似乎在温度 43.5、时间 65、催化剂 12 和 rpm 150 左右获得最佳响应。
如果您坚持将这些建模为因子，可以做到，但您需要将所有预测变量转换为因子。这是一个常见的错误；您可以设计一个仅包含定量变量的几个不同值的实验，但 R 不会读懂您的想法并认为它是一个因素；你必须把它转换成一个。在下文中，我选择了最多具有 2 种交互作用的模型。
> facmod = lm(EE ~ (factor(TEMPERATURE) + factor(TIME) + factor(RPM) + factor(CAT))^2, data = df)
> library(emmeans)
> emmip(facmod, TIME ~ TEMPERATURE | CAT*RPM)
最高拟合响应是在催化剂 12、RPM 150、温度 42.5 或更高以及时间 96 时。很明显，150 RPM 更好（左右比较），高 CAT 更好（比较面板）垂直）。这些是不同的模型和有些不同的结果。我更喜欢 rsm 方法，因为它更系统。