第五章 回归设计

第五章回归设计第1页，课件共78页，创作于2023年2月§5.1回归设计的基本概念回归设计方法是由英国统计学家G.Box在20世纪50年代初针对化工生产提出的。回归设计也称为响应面设计，目的是寻求试验指标与各定量因子间的定量规律，找到工作条件的最优值（最优工艺、最佳配方等）。它是在多元线性回归的基础上用主动收集数据的方法获得具有较好性质的回归方程的一种试验设计方法。第2页，课件共78页，创作于2023年2月5.1.1回归分析——数据处理由被动变主动

古典的回归分析方法只是被动地处理已有的试验数据，对试验点的安排不提任何要求，试验点散乱而不均匀，预测值的标准误很大，且对于回归方程的精度研究也很少。其后果：（1）盲目增加试验次数，这些试验数据还不能提供充分的信息，在许多复因子试验问题中达不到试验目的。（2）对模型的合适性有时无法检验，因为在被动处理数据时在同一试验点上不一定存在重复试验数据。为了适应寻求最佳工艺、最佳配方、建立生产过程的数学模型等的需要，人们就要求以较少的试验次数建立精度较高的回归方程。第3页，课件共78页，创作于2023年2月

为此，要求摆脱古典回归分析的被动局面，主动把试验的安排、数据的处理和回归方程的精度统一起来考虑，即根据试验目的和数据分析的要求来选择试验点，不仅使得在每一个试验点上获得的数据含有最大的信息，从而减少试验次数，而且使数据的统计分析具有一些较好的性质。这就是二十世纪五十年代发展起来的“回归设计”所研究的问题。

回归设计的分类：根据建立的回归方程的次数不同，回归设计通常有一次回归设计、二次回归设计等；根据设计的性质又有正交设计、旋转设计、通用设计和最优设计等。

本章仅介绍二次回归的各种设计方法。第4页，课件共78页，创作于2023年2月5.1.2多项式回归模型

在一些试验中希望建立试验指标y与各个定量因子之间关系的定量表达式，即回归方程，以便通过该回归方程找出使指标满足极值要求的各因子的取值。

可以假定y与间有如下关系：

这里是的一个函数，其图形也称为响应曲面。

是随机误差，通常假定它服从均值为0，方差为

的正态分布。

第5页，课件共78页，创作于2023年2月

试验设计中，我们称为因子或自变量。称

的可能取值的空间为因子空间。我们的任务就是从因子空间中寻找一个最佳工艺条件（最优点）

，使y满足要求。

当f的函数形式已知时，可以通过最优化的方法去寻找

。在许多情况下f的形式并不知道，这时常常用一个多项式去逼近它，即假定：这里各为未知参数，称为回归系数，通常需要通过试验数据对它们进行估计。第6页，课件共78页，创作于2023年2月在实际中常用如下的一次与二次回归方程：若用表示相应的估计，则称为y关于的多项式回归方程。第7页，课件共78页，创作于2023年2月5.1.3多元线性回归

多项式回归模型，在对变量作了变换并重新命名后也可以看成是一个多元线性回归模型。比如对二次回归模型令即变成五元线性回归模型。

1）回归模型

假定回归模型为：

第8页，课件共78页，创作于2023年2月记随机变量的观测向量为回归参数向量为，随机误差向量为结构矩阵上述模型可以表示为矩阵形式：第9页，课件共78页，创作于2023年2月

2）回归系数的最小二乘估计估计回归模型中回归系数的方法是最小二乘法。记回归系数的最小二乘估计为，应满足如下正规方程组：

当存在时，最小二乘估计为：在求得了最小二乘估计后，可以写出回归方程：

第10页，课件共78页，创作于2023年2月

3）对回归方程的显著性检验

对回归方程的显著性检验是指检验如下假设：

H0：

H1：不全为0

则平方和分解式

其中

为残差平方和，自由度为

为回归平方和，自由度为

当H0为真时，有

给定的显著性水平，拒绝域为

第11页，课件共78页，创作于2023年2月4）对回归系数的显著性检验

当回归方程显著时，可进一步检验某个回归系数是否为0，也即检验如下假设：

每一项回归系数j=1,2,…,p逐一进行。常用的检验方法是t检验或等价的F检验，F统计量为：其中是中的第j+1个对角元素。记分子为，它是因子的回归平方和。分母是模型中的无偏估计。第12页，课件共78页，创作于2023年2月

当H0j为真时，有。给定的显著性水平，当时拒绝原假设H0j，即认为显著不为零，回归关系显著；否则人为回归关系不显著，可以将对应的变量从回归方程中删除。注：当有不显著的系数时，一般情况下一次只能删除一个F值最小的变量，重新计算回归系数，再重新检验。通常要到余下的系数都显著时为止。

第13页，课件共78页，创作于2023年2月

5）失拟检验当在某些点有重复试验数据，便可以对试验指标y

的期望是否是的函数进行检验，这种检验称为失拟检验，它检验如下假设：当在有些试验点上有mi重复试验时，试验点为n，总试验次数为N，残差平方和可进一步分解为组内平方和与组间平方和，其中组内平方和就是纯误差平方和，记为，组间平方和称为失拟平方和，记为，即：

第14页，课件共78页，创作于2023年2月，，，，检验统计量为

在H0为真时，，对于给定的显著性水平下，拒绝域为：

当拒绝H0时，需要寻找原因，改变模型，否则接受线性回归模型合适，可以将Se与SLf合并作为SE检验方程是否显著。其中第15页，课件共78页，创作于2023年2月5.1.4因子水平的编码在回归问题中各因子的量纲不同，其取值的范围也不同，为了数据处理的方便，对所有的因子作一个线性变换，使所有因子的取值范围都转化为中心在原点的一个单位“立方体”中，这一变换称为对因子水平的编码。方法如下：设因子的取值范围为：与分别为因子的下水平()与上水平()其中心也称为零水平：

因子的变化半径为令编码值，而实际值此变换式就称为“编码变换”

第16页，课件共78页，创作于2023年2月例5.1.1因子z的取值范围为：10~30，，对其作编码：

编码后，10对应-1，30对应1，20对应0。变换后，正交点在编码空间为中心在原点的立方体，其边长为2。

编码变换后，zm对应的编码为

，zM对应的编码为

，z0对应的编码为0。这样不管什么取值范围，都转化为值域[-1,1]或[-,]。见示意图。第17页，课件共78页，创作于2023年2月§5.2Box－Benhken设计Box－Behnken设计是由统计学家Box和Behnken提出的一种比较常用的回归设计方法，适用于2至5个因子的优化实验。Box－Behnken设计首先假定实验范围内因子存在二次项，其基试验点的选取为编码立方体的每条棱的中点，即任意两因子做22交互，而其它因子固定在0水平。再加上中心点。第18页，课件共78页，创作于2023年2月三因子Box-Behnken设计试验点示意图Box－Benhken设计第19页，课件共78页，创作于2023年2月例题：对超高压杀灭枯草芽孢杆菌效果Y的研究发现：温度、压力、保压时间是灭活枯草芽孢杆菌显著影响因子。研究结果表明杀灭6个数量级的枯草芽孢杆菌的杀菌条件：温度为T=30～60℃，压力为P=200～600MPa，保压时间为M=10～20min，试分析最优杀菌工艺参数。Box－Benhken设计第20页，课件共78页，创作于2023年2月Box－Benhken设计题解：本试验采用Box-Behnken设计，以温度T，压力P，保压时间M三个外界因子为自变量，并以+1、0、-1分别代表自变量的高、中、低水平，对自变量进行编码。超高压杀灭菌的数量级Y为响应值（Y=-log10Nt/N0

，即经超高压作用后枯草芽孢杆菌死亡的数量级，Nt为超高压处理后1ml菌液中的活菌数，N0为对照1ml菌液中的活菌数）第21页，课件共78页，创作于2023年2月Box－Benhken设计因子代号水平编码-101温度(℃)T304560压力(MPa)P200400600保压时间(min)M101520试验因子的水平及编码表第22页，课件共78页，创作于2023年2月Box－Benhken设计试验设计与试验结果列表TrialNo.TPMY1-10-14.27210-15.443-1015.1141015.795-1-102.1161-103.217-1106.0481106.8790-1-12.70100-113.441101-16.23120116.43130005.45140005.32150005.67160005.43170005.23第23页，课件共78页，创作于2023年2月Box－Benhken设计分析结果FactorDFSSMSFPT42.0412470.51031213.670.0020P426.7978746.699469179.46<.0001M40.7164850.1791214.800.0352在0.05水平下，只有温度（T）压力（P）和保压时间（M）与灭菌效果都存在显著的回归关系。在T=60.37,P=663.87,M=13.51时，灭菌效果最大，达到6.79。需要进行试验验证。第24页，课件共78页，创作于2023年2月Box－Benhken设计T=60.37,P=663.87,M=13.51时，极大值Y=6.79第25页，课件共78页，创作于2023年2月§5.3二次回归的中心组合设计

一、中心组合设计方案

中心组合设计中的试验点由三部分组成：（1）将编码值-1与1看成每个因子的两个水平，采用二水平正交表安排试验，可以是全因子试验，也可以是其1/2实施，1/4实施等，称这种试验点为正交点。这样的点有mc=2p个，选取正交表的p个基本列构成。（2）在每一因子的坐标轴上取两个试验点，该因子的编码值分别为-与，其它因子的编码值为0。由于有p个因子，因此这部分试验点共有2p个。称这种试验点为星号点或主轴点。（3）在试验区域的中心进行m0次重复试验，这时每个因子的编码值均为0。称这种试验点为中心点。第26页，课件共78页，创作于2023年2月

如p=2的中心组合设计方案是：第27页，课件共78页，创作于2023年2月

如p=2的中心组合设计试验点的分布图第28页，课件共78页，创作于2023年2月二、中心组合设计方案的特点该方案总试验次数n为：

每个因子都取5个水平，故该方案所布的试验点范围较广。该方案还有较大的灵活性，因为在方案中留有两个待定参数m0（中心点的试验次数）和（主轴点的位置），使二次回归设计具有正交性、旋转性、通用性和D最优性优良性质等。中心点处的m0次重复，可以准确估计试验误差，从而使对方程与系数的检验有了可靠依据。第29页，课件共78页，创作于2023年2月§5.4二次回归正交设计

二次回归正交设计是二次回归中心组合设计的一种常用设计方法。如果一个设计具有正交性，则数据分析将是十分方便的，由于所得的回归系数的估计值之间互不相关，因此删除某些因子时不会影响其它的回归系数的估计，从而很容易写出所有回归系数为显著的回归方程。我们可以适当选择m0与使二次回归中心组合设计具有正交性。第30页，课件共78页，创作于2023年2月5.4.1二次中心组合设计的结构矩阵X与系数矩阵

p=2的中心组合设计回归模型的结构式为

结构矩阵如下：x0x1x2x1x2(x1)2(x2)2第31页，课件共78页，创作于2023年2月这里mc=4，2p=4，则n=mc+2p+m0=8+m0，再记那么各平方项所对应的行、列在非对角线都有非0元素第32页，课件共78页，创作于2023年2月5.4.2正交性的实现

要使中心组合设计具有正交性，就要求为对角阵。首先利用“中心化”变换使各平方项列的和为0，为此把列的元素减去该列的均值，而此时的阵为：

这里GG是p阶对称方阵：

第33页，课件共78页，创作于2023年2月为使设计成为正交的只要设法使g=0因为n=mc+2p+m0设定m0后，因此可以适当选取使g=0。对不同的因子个数p与中心点重复次数m0，对应的值见表5.4.1。

第34页，课件共78页，创作于2023年2月表5.4.1二次回归正交设计的参数值表

第35页，课件共78页，创作于2023年2月

1．回归系数的估计

在对列作了中心化变换后，我们可以首先建立y

关于诸的回归方程：

再记，其中

则各项回归系数：5.4.3统计分析第36页，课件共78页，创作于2023年2月2．对回归方程与回归系数的检验

由于是正交设计，各项的回归平方和为

总回归平方和为仍然用表示总平方和，其自由度为，则残差平方和为回归方程的检验：第37页，课件共78页，创作于2023年2月

若在中心点上有重复试验的话，还可以进一步对进行分解：则失拟项：可用对二次回归模型的合适性进行检验。各项回归系数可用进行检验。第38页，课件共78页，创作于2023年2月

例5.4.1为提高钻头的寿命，在数控机床上进行试验，考察钻头的寿命与钻头轴向振动频率F及振幅A的关系。

在试验中，F与A的变动范围分别为：[125Hz，375Hz]与[1.5，5.5]，采用二次回归正交组合设计，并在中心点重复进行三次试验。1．对因子的取值进行编码现在有两个因子，即p=2，现在中心点进行三次试验，即m0=3，则查得此二次回归正交组合设计中的=1.148。因子F与A的零水平分别是250，3.5;它们的变化半径分别是109，1.74.

因子编码值见表5.4.4。第39页，课件共78页，创作于2023年2月5.4.4第40页，课件共78页，创作于2023年2月

2．试验计划与试验结果

本例的试验计划及试验结果见表5.4.5。表5.4.5试验计划与试验结果

第41页，课件共78页，创作于2023年2月5.4.6

3．参数估计第42页，课件共78页，创作于2023年2月4.模型、方程及系数的检验

本例中由于在中心点有3次重复试验，所以在给出所得到的回归方程之前，先对模型的合适性、方程及系数作显著性检验：中心点上3次试验结果的平均值为=206，由此求得纯误差平方和Se=1026，从而失拟平方和为：SLf=1281.53-1026=255.53，失拟检验的统计量为：

在时，，所以认为模型合适。有关方程与系数的检验见表5.4.7。第43页，课件共78页，创作于2023年2月

由于，所以认为方程显著。又，。所以与的系数在显著性水平0.05上是显著的，x2的系数在显著性水平0.10上是显著的。5.4.7第44页，课件共78页，创作于2023年2月5.写出二次回归方程并求最佳条件

我们可以写出在0.10水平上各系数都显著的回归方程为：再将(5.4.16)代入，即可得y关于x1,x2的二次回归方程：

最后再将编码式代入，即可得y关于F，A的二次回归方程：为延长寿命，可以将回归方程对F与A分别求导，并令其为零以解出最佳水平组合为F=291.58，A=3.50，在该水平组合下，平均寿命的估计是211.6。

第45页，课件共78页，创作于2023年2月5.5二次回归旋转设计5.5.1旋转性回归正交设计的最大优点是试验次数较少，计算简便，又消除了回归系数间的相关性。但是其缺点是预测值的方差依赖于试验点在因子空间中的位置。由于误差的干扰，试验者不能根据预测值直接寻找最优区域。

若能使二次设计具有旋转性，即能使与试验中心距离相等(ρ)的点上预测值的方差相等，即Var()=f(ρ)，那就有助于克服上述缺点。所以试验者常常希望牺牲部分的正交性而获得旋转性，计算的工作量可以交由计算机帮助处理。

第46页，课件共78页，创作于2023年2月5.5.2二次旋转设计

一个中心组合设计要成为二次旋转设计应满足旋转性条件和非退化条件。

1.二次设计的旋转性条件

二次设计的旋转性条件为：第47页，课件共78页，创作于2023年2月2.二次旋转设计的非退化条件：

为使设计是使矩阵不退化，就要求试验点的分布满足：在中心组合设计方案中n个试验点分布在三个不同半径的球面上，其中：个点分布在半径为的球面上；

2p个点分布在半径为的球面上；个点分布在半径为的球面上。满足不会使矩阵

退化的条件。第48页，课件共78页，创作于2023年2月3.的选取为使设计满足旋转性条件只要适当选取参数，在中心组合设计中有：

因此，，为使设计具有旋转性，则要求

即只要：或当对中心组合设计提出进一步的要求时，可以确定设计中的另一个参数m0。

第49页，课件共78页，创作于2023年2月5.5.3二次回归正交旋转设计（Orthogonal）

当要求一个回归设计不仅具有旋转性，还要求具有正交性，或至少是近似正交的，称为二次回归正交旋转设计。这需要使的非对角线元素全为0，即只需要：

在g的表达式中，mc是给定的，因为满足旋转性也已确定，式中，所以g只是m0的函数，所以可解出m0。如果解得的m0是整数，则所得设计为正交旋转设计；如果m0不是整数，则取最接近的整数，这时是近似正交的旋转设计。

第50页，课件共78页，创作于2023年2月

二次回归正交（或近似正交）旋转组合设计的参数与m0见表5.5.1。表5.5.1二次回归正交旋转组合设计参数第51页，课件共78页，创作于2023年2月5.5.4二次回归通用旋转设计(UniformPrecision)

所谓一个设计具有通用性是指在与编码中心距离小于1的任意点(x1,x2,…,xp)上的预测值的方差近似相等。

一个旋转设计各点预测值的方差仅与该点到中心的距离ρ有关，即Var()=f(ρ)。而通用设计要求当ρ<1时，f(ρ)基本为一个常数。根据这一要求，可以通过数值的方法来确定m0。当一个回归设计既要具有旋转性又要具有通用性时，称为二次回归通用旋转设计，设计的参数与m0见表5.5.2。

第52页，课件共78页，创作于2023年2月5.5.2第53页，课件共78页，创作于2023年2月5.5.5数据分析

由于正交旋转设计的数据分析同前面的回归正交设计一样，所以下面仅对通用旋转组合设计的数据分析作介绍。

1．回归系数的估计如果记X’Y阵中的元素为：

根据不同的p与实施方案，计算中用到的K，E，F，G的值已列成表格供使用。第54页，课件共78页，创作于2023年2月则各项回归系数的估计为：第55页，课件共78页，创作于2023年2月

2．对回归方程的检验总平方和：残差平方和：回归平方和：各自由度分别为：第56页，课件共78页，创作于2023年2月由于在中心点有m0次重复试验，因此还可将SE分解为：

其自由度分别为：

先检验模型的合适性，即失拟检验：

当模型合适时，再检验方程的显著性：第57页，课件共78页，创作于2023年2月

3．对回归系数的显著性检验

各项的回归平方和及σ2的估计，如下：检验各项回归系数的统计量依次为：

第58页，课件共78页，创作于2023年2月

如果有不显著的项，要删去该项，一次只能剔除一项，由于这里不是正交设计，所以回归系数间具有相关性，删除一个变量后，回归系数需要重新计算。

由于求回归系数的正规方程组的系数矩阵阶数较高，求逆矩阵相当麻烦，通常将这项工作交给计算机协助完成。统计分析方法和前述一样，用反应面回归分析。第59页，课件共78页，创作于2023年2月5.5.6二次通用旋转设计实例

例5.5.1

某三因子（x1,x2,x3）优化试验，试验指标Y具有望大特性，采用二次通用旋转设计，查表得到γ=1.682,n=20，各因子的水平编码见下表。1.因子水平编码表因子编码x1x2x3+γ80120300169.8699.73239.1705570150-140.1440.2760.83-γ30200第60页，课件共78页，创作于2023年2月2.试验设计与结果NOX1X2X3Y11-1-1-130.42-1-1145.33-11-115.94-11132.251-1-143.561-1172.0711-124.2811148.59-1.68180028.4101.68180053.6NOX1X2X3Y1110-1.6818056.81201.6818029.11300-1.681814.314001.681852.11500043.31600042.01700043.51800043.61900042.32000044.0第61页，课件共78页，创作于2023年2月3.模型检验Source

DFSSMSFPModel93696.605410.7338371.8610.0001Linear33431.6311143.8771035.6180.0001Quadratic3176.793658.9312153.353810.0001CrossProduct388.1829.3933326.611470.0001Error1011.045361.104536Lackoffit57.9370271.5874052.5534670.1633PureError53.1083330.621667Total193707.65第62页，课件共78页，创作于2023年2月4.参数估计及其检验TermEstimateStdErrtPr>|t|X17.8188740.2843927.493460.0001X2-8.566080.28439-30.12090.0001X310.80570.2843937.996020.0001X1*X1-0.706850.276846-2.553230.028705X1*X2-1.90.371574-5.113390.000455X1*X32.70.3715747.2663920.0001X2*X2-0.017410.276846-0.06290.951088X2*X3-0.350.371574-0.941940.368422X3*X3-3.464570.276846-12.51440.0001第63页，课件共78页，创作于2023年2月5.去掉不显著的项重新回归——模型检验SourceDFSSMSFPr>FModel73695.62527.9458526.6410.0001Error1212.029731.002478(Lackoffit)78.9213971.2744852.050110.223309(PureError)53.1083330.621667Total193707.65第64页，课件共78页，创作于2023年2月6.重新回归后——参数估计TermEstimateStdErrtPr>|t|X17.8188740.27093328.859050.0001X2-8.566080.270933-31.6170.0001X310.80570.27093339.883270.0001X1*X1-0.705120.262443-2.686770.019788X1*X2-1.90.353991-5.367370.000169X1*X32.70.3539917.6273110.0001X3*X3-3.462850.262443-13.19470.0001第65页，课件共78页，创作于2023年2月7.回归方程及其驻点Y1=43.096+7.819*X1

-8.5669*X2

+10.806*X3

-0.705*X1*X1-1.9*X1*X2

+2.7*X1*X3

-0.463*X3*X3驻点处X1-4.53107鞍点X26.711891Y=-6.296

X3-0.54514特征值0.865802-1.02896-4.02568在驻点处其二次型矩阵特征值有正有负，此驻点为鞍点，不存在极值，需进行岭脊分析。第66页，课件共78页，创作于2023年2月8.岭脊分析半径Y预测值标准误X1X2X3043.1103180.4286350000.145.7716310.4272330.086905-0.0932550.1097060.248.4318250.4238690.180418-0.1906540.2103300.351.1010420.4211250.278880-0.2914100.3030960.453.7871420.4233790.381000-0.3948910.3891690.556.4962010.4364940.485803-0.5006040.4695790.659.2329320.4666900.592567-0.6081610.5452040.762.0010150.5188080.700760-0.7172610.6167770.864.8033500.5951280.809989-0.8276660.6849050.967.6422440.6956180.919966-0.9391880.7500891.070.5195490.8190191.030477-1.0516750.812741第67页，课件共78页，创作于2023年2月8.岭脊分析（续）半径Y预测值标准误X1X2X31.173.436770.9637831.141363-1.1650.8732051.276.395141.1285371.252504-1.279070.9317661.379.395671.3121991.363813-1.393790.9886631.482.439221.5139641.475222-1.509091.0440961.585.526481.733241.586681-1.624911.0982361.688.658071.9695961.698153-1.74121.1512281.791.83452.2227141.809609-1.857921.2031961.895.056192.4923591.921028-1.975021.2542471.998.323552.7783522.032395-2.092471.3044732.0101.63693.080562.143697-2.210241.353954第68页，课件共78页，创作于2023年2月5.6D最优混合设计5.6.1D最优设计回归设计的正交性、旋转性和通用性各有其优点，但是都没有涉及统计意义上的优劣。

最优设计就是从统计意义上来研究不同试验方案的优劣并建立最优方案。在所有的试验方案中，能使回归预测值与观察值拟合最好的那个方案就是最优方案，即最优设计。判断试验方案的优劣有很多标准，如：D优良性、G优良性、A优良性、E优良性、I优良性、U优良性等。其中以D优良性的使用最为普遍，有人证明：事实上G优良性与D优良性等价。第69页，课件共78页，创作于2023年2月5.6.1D最优设计

1943年，Wald提出信息矩阵的行列式极大值判别法。结构矩阵和信息矩阵与试验方案有关，设试验方案为ε，可以证明，试验方案的信息矩阵A(ε)的行列式的值|A(ε)|越大，回归系数的预测值的方差越小，预测值的方差也越小，试验方案也就越优良。以此为依据，从各种试验方案中挑选最优设计方案的方法，称为D最优设计。在同一模型下，两个试验方案ε1和ε2，如果方案ε1的信息矩阵行列式的值大于ε2的值，即：|A(ε1)|>|A(ε2)|，则说明在D优良性下方案ε1比ε2有更好的优良性。第70页，课件共78页，创作于2023年2月5.6.2D最优混合设计

D优良设计一般不满足正交性和旋转性，于是有人以回归组合设计为基础，提出了部分兼有近似D优良性、正交性和旋转性的D最优混合设计。

D最优混合设计的设计方案处理数较少，试验效率较高，回归系数和待估参数之间拟合得较好，统计分析上也较为方便。到目前为止，人们已经找到了若干2至6个因子的D最优混合设计方案。如:R206,R207,R208,R310A,R310B,R311A,R311B,R312,R3414,R415,R416A,R416B,R416C,R521,R628A,R628B。其中R表示回归设计，第一个数字表示因子数，第二三位数字表示处理个数，最后的字母表示相同因子的方案编号。第71页，课件共78页，创作于2023年2月5.6.2D最优混合设计——方案示例R311ANOX1X2X3Y10