原创内容，欢迎转载，转载请注明出处
主笔：于浩
导读
1、条件筛选（optimization of experimental conditions）是优化反应体系、工艺参数或者生物过程中的关键步骤。为了高效、系统、准确地确定最优条件，研究者通常采用一系列的实验设计（Design of Experiments, DoE）方法。在科研论文中常见的条件筛选有：
（1）单因素实验（One-Factor-at-a-Time, OFAT）
（2）正交实验设计（Orthogonal Array Design / Taguchi Method）
（3）Plackett-Burman 设计（P-B Design）
（4）响应面法（Response Surface Methodology, RSM）
（5）其他方法：人工神经网络+遗传算法、机器学习辅助优化、全因子试验、单纯形格点法等
2、在上面的所有的DoE中，RSM（尤其是 BBD ）是当前 SCI 论文中条件优化最主流、最被推崇的方法。因此如果的确要进行条件筛选，完全可以用响应面法。
在发表中英文科研论文的时候，最好采用完整的实验条件优化策略，也就是采用 ”PB实验+最陡爬坡实验+响应面“ 或者 ”单因素+PB实验+最陡爬坡实验+响应面实验“ ，尤其适合于初始变量较多但是真正的关键因子较少且不清楚的情况。 虽然说结果可能跟”单因素+响应面“差不多，但是从实验设计严谨性上来说这是一个最优的方案。这样审稿人会认为你的实验遵循了科研规范，更加容易通过
3、PB筛选和响应面设计分析的工具有很多，比如说R语言、python、matlab，但是最常使用的是收费软件Design-Expert软件，这个软件可以免费下载使用14天。出于知识产权的考虑，如果可以用免费的软件来做的话，建议采用免费的软件做。
如果使用了本教程的python脚本，别忘了引用本网站，谢谢！
一、单因素实验
是否需要进行单因素实验
PB设计就是为了替代复杂的单因素实验而存在的，因此实验设计中可以不使用单因素实验。
在使用 Plackett-Burman（PB）设计 之前，通常不需要系统性地进行完整的单因素实验（即“一次只改变一个因素”的实验），但需要做一些必要的前期准备工作，包括对各因素的水平范围进行初步探索或合理设定，当然这些信息我们完全可以从之前的文章中获取，基本上大家根据前人研究就可以筛选出来主要的一些因素，用于PB设计。
比如说我们要进行某个食用菌液体菌种的发酵条件优化，通过文献就可以知道了主要的因素就是 碳源、氮源、天然基质、无机盐、装液量、接种量、转速、温度等。
什么情况下需要进行单因素实验
在实际科研中，在我们发表高水平论文的时候，实际上我们肯定要对每一个因素进行单因素实验筛选，确定每个因素的大致作用范围和影响显著性，不可能完全依赖于PB筛选，基本上所有的PB筛选的因素的范围也都是单因素筛选之后的。
因此 单因素+响应面 也是一种非常合理的方案，并不能说 ”PB+爬坡+响应面“ 就一定比 ”单因素+响应面“ 高级。
如果只是一个简单的本科毕业设计那么也可以不进行单因素实验。
第2种情况是我们的实验想要引入一些新的因素，这些因素在之前的研究中从来没有报道过，我们突发奇想或者灵光乍现，感觉这个因素会有显著影响，这个时候就可以通过单因素实验来进行初步的筛选，看看因素是否对于结果影响显著，因素的大体的作用范围是什么。
此外，还有一种情况是PB筛选的结果有时候并不是我们需要的结果，有些时候我们会根据经验或者生产实际知道某些因素是重要的，但是PB筛选的结果下某一个因素并不在前三名，这个时候我们可以把范围扩大一下，通过添加单因素筛选这一步，把我们想要的因素给进一步筛选出来。
本教程不涉及单因素实验部分内容，就是纯粹的“PB+爬坡+响应面”实验。
二、Plackett-Burman设计（PB筛选）
PB原理
Plackett-Burman（PB）筛选实验是一种两水平的因子筛选实验，当初始阶段筛选的因子太多了，但是又不清楚哪一个因素是主要因素，PB筛选可以在较少实验的次数下快速识别显著影响 响应值 的关键因子。
主要目的就是快速筛选出来对响应值又显著影响的关键因子，用于后续的爬坡实验和响应面实验进行精细化条件筛选。
PB筛选假设高阶交互作用可以忽略，仅评估各因子的主效应。
PB筛选的实验次数 n 为4的倍数（如8、12、16、20等），最多可以一次性实验筛选n-1个因子，例如12次实验最多可以筛选11个因子（也可以筛选8个、9个和10个因子）
实验内容
1、实验目的：筛选液体发酵培养杏鲍菇菌丝生产麦角甾醇的最优条件。
2、初始的影响因素包括：碳源（葡萄糖）、碳源（蛋白胨）、无机盐（KH2PO4）、装液量、转速、接种量、培养时间、培养温度。
PB筛选实验设计和结果分析
1、因素水平设计表
| 编号 | 因素 | 低水平 (-1) | 高水平 (1) |
|---|---|---|---|
| A | 葡萄糖Glucose (碳源, g/L) | 20 | 40 |
| B | 蛋白胨 Peptone (氮源, g/L) | 5 | 15 |
| C | KH2PO4 (无机盐, g/L) | 1 | 2 |
| D | 装液量Loading_Vol (mL/250mL) | 50 | 100 |
| E | 转速Speed (rpm) | 120 | 180 |
| F | 接种量Inoculum (%) | 5 | 15 |
| G | 培养温度 Temp (℃) | 22 | 28 |
| H | 培养时间Time (d) | 5 | 9 |
2、利用Design Expert设计PB筛选设计
点击“杂项/PB筛选试验”，右侧先选择“因子数“，下面可以不填任何因素名字（因为完全不影响实验分析，后面可以看到），也可以像这样输入各因子的名字，输入单位和 low值、High值，输入了之后直接生成的PB筛选表格就可以带着实验参数。否则还需要自己把”1、-1“转化为实验参数。
这一步也可以什么都不该，直接点击“前进”，因为不影响分析。如果有多个响应数，响应数可以选择多个。
如果前面没有输入各因子的参数，这里生成的就是1、-1矩阵，如果输入了就直接生成实验参数矩阵。
不同电脑生成的矩阵不同，但是都可以用，大家也可以直接参考我们的矩阵，放在本部分最后面。后面还有几列（因为我们是8个因素，最多可以验证11个因素，所以后面还有3列，这几列可以删除，也可以不删除，不影响后续的分析）
3、利用Design Expert分析PB筛选结果
文章中主要展示两个结果：1、方差分析结果；2、Pareto 效应图。
把结果填入最后面响应值这一列，点击左侧的“分析”下面的响应值的图标。
下面这一步很关键，进入“效应/半正态”，在右侧的数值窗口中，把要分析的这些因素选中，右键点击“模型“，只有变成模型，才能进行ANOVA方差分析，否则不能分析。
点击“帕累托”这个窗口，就出现了第一张文章中可以用的图，但是这个图不能直接用，大家需要重新自己处理，或者用python生成，后面会介绍。这个图像直观的展示了各个因素的重要性，值越高重要性越大，可以根据这个值排序。
最后进行ANOVA方差分析
这个表的结果整理一下就可以放到文章里面了。
ANOVA分析里面比较重要的就是这个p-value的表格，也是很多中文文章里面直接放的表格，但是实际上下面的这个回归系数（Coefficient Estimate）也很重要，这个值的正和负关系到爬坡实验的方向，因此发表文章的时候建议也跟p-value放到一个表格里面。
到这里PB筛选实验和分析就完成了。
4、下面是本教程使用的模拟数据
| 序号 | Glucose | Peptone | KH2PO4 | Loading_Vol | Speed | Inoculum | Temp | Time | Ergosterol |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | -1 | 1 | 1 | 1 | -1 | -1 | -1 | |
| 1 | -1 | 1 | 1 | 1 | -1 | -1 | -1 | 1 | |
| 2 | 1 | 1 | 1 | -1 | -1 | -1 | 1 | -1 | |
| 3 | 1 | 1 | -1 | -1 | -1 | 1 | -1 | -1 | |
| 4 | 1 | -1 | -1 | -1 | 1 | -1 | -1 | 1 | |
| 5 | -1 | -1 | -1 | 1 | -1 | -1 | 1 | -1 | |
| 6 | -1 | -1 | 1 | -1 | -1 | 1 | -1 | 1 | |
| 7 | -1 | 1 | -1 | -1 | 1 | -1 | 1 | 1 | |
| 8 | 1 | -1 | -1 | 1 | -1 | 1 | 1 | 1 | |
| 9 | -1 | -1 | 1 | -1 | 1 | 1 | 1 | -1 | |
| 10 | -1 | 1 | -1 | 1 | 1 | 1 | -1 | -1 | |
| 11 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| 序号 | Glucose | Peptone | KH2PO4 | Loading_Vol | Speed | Inoculum | Temp | Time | Ergosterol |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 40 | 5 | 2 | 100 | 180 | 5 | 22 | 5 | 3.05 |
| 1 | 20 | 15 | 2 | 100 | 120 | 5 | 22 | 9 | 3.41 |
| 2 | 40 | 15 | 2 | 50 | 120 | 5 | 28 | 5 | 3.42 |
| 3 | 40 | 15 | 1 | 50 | 120 | 15 | 22 | 5 | 3.24 |
| 4 | 40 | 5 | 1 | 50 | 180 | 5 | 22 | 9 | 3.15 |
| 5 | 20 | 5 | 1 | 100 | 120 | 5 | 28 | 5 | 2.44 |
| 6 | 20 | 5 | 2 | 50 | 120 | 15 | 22 | 9 | 3.02 |
| 7 | 20 | 15 | 1 | 50 | 180 | 5 | 28 | 9 | 3.15 |
| 8 | 40 | 5 | 1 | 100 | 120 | 15 | 28 | 9 | 3.58 |
| 9 | 20 | 5 | 2 | 50 | 180 | 15 | 28 | 5 | 2.32 |
| 10 | 20 | 15 | 1 | 100 | 180 | 15 | 22 | 5 | 3.18 |
| 11 | 40 | 15 | 2 | 100 | 180 | 15 | 28 | 9 | 3.95 |
三、最陡爬坡试验
最陡爬坡实验的设计
1、最陡爬坡实验的目的是什么
PB（Plackett-Burman）实验只能筛选出哪些因素显著，但无法确定这些因素的最佳取值范围。爬坡实验的目的是快速逼近最大响应值（产量）所在的区域，为后续的响应面实验（如 BBD）寻找中心点。
基于 PB 实验建立的一阶线性模型：Y = β0 + Σβi Xi
根据回归系数βi（也就是上面的Coefficient Estimate）的正负和大小，确定各因素的改变方向和步长。通过改变各因素的水平，使实验结果沿着梯度最陡的方向（即响应值增加最快的方向）移动。
2、最陡爬坡实验一般选择几个因素，怎么选择这些因素
通常筛选 PB 实验中显著性较高（p < 0.05）的 3-5 个因素。非显著因素在爬坡实验中固定在 PB 实验时的较优水平。首先应该看p-value的值，从小往大选择，还可以参考帕累托图的值，也就是Effect值，从大往小选择。
如果在响应面中计划选用三个因素，那么爬坡就选择三个就行，如果响应面打算用4因素3水平，那么爬坡就选择4个就行。
3、最陡爬坡实验的具体实验参数如何设置
方向如何确定：
如果回归系数βi是正的，那么说明了增加该因素水平可提高产量，实验中应调高取值；反之则调低。
起点和步长如何确定：
确定起点和步长没有绝对的数学标准，可以采用经验值或者单因素的结果，也可以参考PB的结果。
比如说用PB的地点作为起点，跨越PB筛选的范围。
以PB筛选的中点作为起点，进行爬坡。
不同参数的步长如何确定：
不同参数的步长比例可以参考t值（也就是回归系数/回归系数标准差，python脚本中的t值），在设置多个显著性因素的步长的时候，通常步长与t值（或者回归系数的绝对值）的大小成正比。我们的例子中Glucose的t值是4.9，Peptone是4.8，基本一致，那么在步长设置的时候就应该一样的步长。
最陡爬坡实验的例子和分析
1、三因素的爬坡实验
根据上面的PB筛选的结果，Glucose、Peptone、Time这三个的Effect效应值最高，p-value小于0.05，因此选择这三个因素爬坡。因为这三个的Coefficient Estimate都是正的，所以三个爬坡方向都是从小到大。虽然三个的t值都一样，显然Time不可能跟Glucose一样的步长，所以根据经验，Time的步长设置为2，Glucose的步长设置为15，Peptone的步长设置为8。
爬坡实验计划做5组，具体的设计和随机生成的结果如下：
| 实验组 | Glucose（g/L） | Peptone（g/L) | Time （d） | ergosterol（mg/g DW） |
|---|---|---|---|---|
| step 1 | 20 | 5 | 5 | 3.2 |
| step 2 | 35 | 13 | 7 | 3.9 |
| step 3 | 50 | 21 | 9 | 4.4 |
| step 4 | 65 | 29 | 11 | 4.5（最大值，响应面中心点） |
| step 5 | 80 | 37 | 13 | 4.0 |
2、结果分析
从结果可以看出来，随着各因素增加，麦角甾醇含量提高，step 4达到最大值（顶点），进一步提高参数麦角甾醇含量下降。因此选择step 4的各因素水平作为响应面的中心点，即 Glucose 65 g/L、Peptone 29 g/L、Time 11 d。
3、四因素的爬坡实验
假设还是上面的实验，假设Speed的p-value也显著，那么我计划在响应面中做4因素3水平的实验，我想在爬坡的时候把Speed也加上。通过上面发现Speed的Coefficient Estimate是负值，也就是说转速高了不利于麦角甾醇的生成，那么上面的爬坡实验就应该按照下面的方法设计，Speed从大到小。
| 实验组 | Glucose（g/L） | Peptone（g/L) | Time （d） | Speed（rpm） |
|---|---|---|---|---|
| step 1 | 20 | 5 | 5 | 180 |
| step 2 | 35 | 13 | 7 | 160 |
| step 3 | 50 | 21 | 9 | 140 |
| step 4（最大值） | 65 | 29 | 11 | 120 |
| step 5 | 80 | 37 | 13 | 100 |
4、影响爬坡实验结果的因素
PB 结果决定了爬坡的起点、方向和各因素改变的比例。如果 PB 模型失拟或者显著性判断错误，爬坡实验将无法找到最高点。
因此如果爬坡发现最高点出现在起点或者终点，则需要调整起点和步长，重新进行爬坡。比如说步长过大： 极易“跨过”最优点。如果你发现 Step 1 产量最高，Step 2 就大幅下降，说明你已经跨过了峰值，必须缩小步长重新从 Step 1 附近爬坡。
如果步长过小： 实验组别增加，且可能因为处于误差波动范围内而看不出产量梯度，导致爬坡失败。
四、响应面(Box-Behnken Design) 优化
前面的一个教程用Design Expert v8.0介绍了一个3因素3水平的响应面实验应该如何设计并分析结果。这里我们根据本实验，利用Desing Expert v13.0设计一个4因素3水平的响应面。
做12组实验（+5个中心点一共17组实验），最高的因素和水平的组合就是4因素3水平。如果5因素3水平就要做24组实验（+5个中心点就是29组实验，达不到降低实验次数的效果）
BBD响应面只有3个水平：高中低（对应的值是1、0、-1）不存在4水平的BBD响应面。
如何确定BBD响应面的中心点和范围
1、中心点确定
BBD（Box-Behnken Design）响应面实验的设置与 PB 筛选和最陡爬坡实验的结果密切相关。这三者构成一个逻辑连贯、层层递进的优化流程，BBD 并不是孤立设计的，而是建立在前两步结果基础上的精细建模阶段：
PB 选因素 -> 爬坡定中心 -> BBD 算坐标
因此中心点应该就使用爬坡实验确定的最大值对应的值。
2、高低范围的确定
在中心点基础上，向上和向下浮动一个步长。这个步长通常可以跟爬坡步长一样，也可以小一些。
究竟应该设置多少，还是一个经验值。
3、本示例的值的确定
我们这里还是设计3因素3水平的响应面实验，也就是使用上面的Glucose、Peptone、Time这三个因素，中心点采用爬坡实验给出来的65、29、11这三个数值。高低的水平就采用步长就可以了，具体的BBD参数如下：
（1）Glucose (g/L)：50、65、80
（2）Peptone (g/L)：21、29、37
（3）Time (d)：9、11、13
利用Design Expert来设计BBD响应面实验
1、打开Design Expert软件
从左侧菜单栏选择“响应面/ 随机化 / BBD响应面法”，在右侧选择因子的数量，默认是3，跟我们实验一样。下面输入每个因子的信息，如果不输入后面生成的就是“-1、0-1”矩阵，如果输入了后面就会生成实际实验的矩阵。
中心点默认就是5，我们常用的也是5个，不需要修改。
一个实验可以研究不同的响应面，比如说麦角甾醇浓度、几丁质浓度、三萜浓度，有几个响应值就选择几。
下面输入响应值的信息。
下面是生成的实验矩阵，不同电脑生成的可能不一样，大家可以采用本教程的01矩阵，改写上自己的实验参数，也可以根据生成的来做实验。
注意一定要保存好数据，下次打开有可能就不一样了。
根据实验矩阵做实验，把结果填到表格里面。
点击左侧的“分析”下面的响应值，开始进入分析界面。
点击配置下面的 “开始执行分析”
响应面最核心的结果（放到文章里面的），就是ANOVA方差分析，下面列出来了几个重要的注意的点。
点击上面的 “模型图”，就可以查看响应面的图像了，右侧的 “项” 可以选择不同的因子组合。
如果要查看3D图，需要点击 ”等值线“ 右侧的小三角，在下划菜单中选择 ”切换图表 / 3D曲图“，就可以查看并调整3D图像了。
3组的等高线图和3D图都下载下来，组合到一起就是文章中的图。
如果要寻找最优条件和理论最高值，可以按照下面的方法来分析。
2、下面是本教程使用的原始数据（模拟）
| 序号 | A: Glucose | B: Peptone | C: Time | 麦角甾醇 (mg/g DW) | 备注 |
|---|---|---|---|---|---|
| 1 | -1 (50) | -1 (21) | 0 (11) | 3.82 | 边界点 |
| 2 | 1 (80) | -1 (21) | 0 (11) | 4.05 | 边界点 |
| 3 | -1 (50) | 1 (37) | 0 (11) | 3.98 | 边界点 |
| 4 | 1 (80) | 1 (37) | 0 (11) | 4.28 | 边界点 |
| 5 | -1 (50) | 0 (29) | -1 (9) | 3.75 | 边界点 |
| 6 | 1 (80) | 0 (29) | -1 (9) | 4.02 | 边界点 |
| 7 | -1 (50) | 0 (29) | 1 (13) | 3.88 | 边界点 |
| 8 | 1 (80) | 0 (29) | 1 (13) | 4.12 | 边界点 |
| 9 | 0 (65) | -1 (21) | -1 (9) | 3.92 | 边界点 |
| 10 | 0 (65) | 1 (37) | -1 (9) | 4.15 | 边界点 |
| 11 | 0 (65) | -1 (21) | 1 (13) | 4.08 | 边界点 |
| 12 | 0 (65) | 1 (37) | 1 (13) | 4.25 | 边界点 |
| 13 | 0 (65) | 0 (29) | 0 (11) | 4.52 | 中心点 (最高) |
| 14 | 0 (65) | 0 (29) | 0 (11) | 4.48 | 中心点 |
| 15 | 0 (65) | 0 (29) | 0 (11) | 4.51 | 中心点 |
| 16 | 0 (65) | 0 (29) | 0 (11) | 4.49 | 中心点 |
| 17 | 0 (65) | 0 (29) | 0 (11) | 4.5 | 中心点 |
五、利用python脚本进行PB筛选+最陡爬坡+BBD响应面（推荐）
Design Expert还是收费的，如果单位和个人没有购买，直接使用可能会存在问题。其实Python完全可以做同样的任务，而且分析结果完全一致。
不过用python来做分析，需要具备一定的python知识。需要在自己Windows电脑上面安装python和使用的包，也可以在Linux服务器上面进行操作。
另外本教程的python脚本可能要根据自己的实际实验调整，比如说如果PB筛选中因子数量少于8个，可能脚本给出来的就是一个8组实验的设计，如果要做12组实验，只能用空的因子填充。所以脚本要根据自己的试验情况自己调整。
Plackett-Burman设计的Python脚本
1、点击下载python脚本 PB脚本。
下载脚本后，需要把实验的因素，每个因素的范围，实验结果填到python脚本中
2、在powershell中或者linux系统中运行脚本。
一键可以得到所有的结果了，包括实验矩阵，1、-1矩阵（生成excel表格），ANOVA分析结果和帕累托图。比DE简单多了。
BBD响应面的Python脚本
1、点击下载python脚本 BBD响应面脚本。
下载脚本后，需要把实验的因素，每个因素的范围，实验结果填到python脚本中
2、在powershell中或者linux系统中运行脚本。
一键可以得到所有的结果了，包括01矩阵、实验矩阵、ANOVA分析结果和3D图。比DE简单多了。
评论区