第十讲(1)回归正交设计

第七章回归的正交设计

古典回归是被动的，对试验不做任何要求，回归方程精度也很少研究。这不仅盲目增加了试验次数，且数据信息不充分，达不到试验目的。

寻求最佳农艺措施和模式化载培及建立数学模型，主动地把试验安排、数据处理和回归方程的精度统一起来考虑。即选择试验点使数据含有最大信息、减少试验次数、数据分析具有较好的性质。这就是“回归设计与分析”研究的问题。

§7.1什么是回归设计1

如粮食生产农艺最优化,在不完全了解生产过程物理、化学和生物原理的情况下,回归分析是一个有效方法。

类似“黑箱”理论，输入和输出已知，但内部结构不清楚。输入是一些因子z1,z2,…，zp(如N肥、P肥、密度等)，输出就是要最优化的农艺指标y(如产量)，当然可以有二个或更多的农艺指标。

y总可以表示为因子z1,z2,…，zp的一个函数，即

y=f(z1,z2,…,zp)并且可以用一个多项式来近似表达这个函数。2

以因子z1,z2,…，zp

为坐标的空间称为因子空间，寻找最优化农艺措施，也就是要在因子空间中寻找y为最大值或最小值的点。

如何在因子空间中选择适当的试验点，使得能以最快的速度建立一个有效的多项式，从而解决最优化问题呢？这就要用到回归的设计与分析这一工具。

回归设计内容相当丰富，有回归正交、旋转、Ｄ-最优设计等。主要介绍这三种设计。按照回归模型的次数，回归设计又分为一次回归设计、二次回归设计等，首先讨论一次回归正交设计。

3

一次回归的正交设计主要是运用二水平正交表（如L4(23),L8(27),L12(211)等）进行设计，设计与分析的主要步骤和性质如下：1.确定因子的变化范围§7.2一次正交回归设计为因子zj的零水平，

如要研究p个因子z1,z2,…，zp与某项指标y的数量关系，那么首先要确定每个因子zj变化的下界和上界，假如试验就在水平z1j和z2j上进行，那么分别称z1j和z2j为因子zj的下水平和上水平，并称它们的算术平均41.确定因子的变化范围为因子zj的变化区间。

如要研究p个因子z1,z2,…，zp与某项指标y的数量关系，那么首先要确定每个因子zj变化的下界和上界，假如试验就在水平z1j和z2j上进行，那么分别称z1j和z2j为因子zj的下水平和上水平，并称它们的算术平均为因子zj的零水平，它们差的一半5例7.1

为了研究小麦高产栽培技术，选择影响小麦产量三个主要因素水分状况z1（全生育期土壤湿度占田间持水量％,75~95）、肥料z2(公斤/亩，20~40)、密度z3(万株/亩，45~65)进行回归正交试验。

z1:z11=75,z12=95,z01=(75+95)/2=85,1=(95-75)/2=10z2:z12=20,z22=40,z02=(20+40)/2=30,2=(40-20)/2=10z3:z13=45,z23=65,z03=(45+65)/2=55,3=(65-45)/2=1062.对每个因子zj的水平进行编码编码就是对因子的取值作如下线性变换：这样就建立了因子zj与xj取值的一一对应关系:下水平z1j←→-1零水平z0j←→0上水平

z2j←→+1

具体编码工作可在表上进行。7

从上面知道，当zj在区间［z1j,z2j］内变化时，它的编码值xj就在区间［-1,+1］内变化。因子水平编码表因子

z1

z2……zp

下水平(-1)零水平(0)上水平(+1)

z11

z12……z1pz01

z02……z0pz21

z22……z2p变动区间(j)

1

2

……p

对因子zj编码后。y对z1,z2,…,zp的回归问题，就转化为y对x1,x2,…,xp的回归问题。因此可以在x1,x2,…,xp为坐标轴的编码空间中选择试验点进行回归设计。8对例7.1的因子水平编码如下表。因子水平编码表因子

z1

z2

z3

下水平(-1)零水平(0)上水平(+1)

752045853055954065变动区间(j)

10101093.选择适当的二水平正交表

运用二水平正交表时，需用“-1”代换通常的二水平正交表中“2”。代换后表中的“+1”既表示因子水平的不同状态，也表示因子水平变化的数量大小。经代换后，正交表的交互作用列还可直接由表中相应几列的对应元素相乘而得到，因此交互列也就不必要了。如L4(23)试验号

x1

x2

x1x21234

1111-1-1-11-1-1-1110

究竟用哪一张二水平正交表，要根据因子的个数而定。正交表确定以后，把各变量放入正交表的某些列上，把这些列取出就组成了一张试验计划。如把三个变量放在正交表L8(27)的x1,x2,x3列上，就得到一张如下的全因子试验计划。试验号

x1

x2

x312345678

11111-11-111-1-1-111-11-1-1-11-1-1-111

用二水平正交表编制的试验计划具有正交性。显然，若以xij表示第i次试验中第j个变量的编码值，于是在试验计划中有任一列的和、任两列的内积等于0，即试验号

x1

x2

x312345678

11111-11-111-1-1-111-11-1-1-11-1-1-1

由此称这种设计为正交设计，所得的的试验计划为正交计划。

12对例7.1,选L8(27)正交表,根据因子水平编码表,“因素顺序上列,水平对号入座”,每行为一个试验处理,得试验计划如下:试验号

x1(z1)x2(z2)x3(z3)12345678

1(95)1(40)1(65)11-1(45)1-1(20)11-1-1-1(75)11-11-1-1-11-1-1-1试验计划134.回归系数的计算与统计检验

根据正交设计进行n次试验,试验结果为y1,y2,…,yn。则一次回归的数学模型为

于是参数的最小二乘估计yi=0+1xi1,+2xi2+…+pxip+ei，i=1,2,…,nb=(XTX)-1XTY=A-1B由正交设计的性质知，A-1=diag(1/n,1/n,…,1/n)B=XTY=(B0,B1,…,Bp)T，B0=Σyi，Bj=Σxijyi，j=1,…,p从而b0=Σyi/n，bj=Σxijyi/n，j=1,…,p(1)回归系数的计算

14由于COV(b)=2(XTX)-1=2diag(1/n,1/n,…,1/n),所以回归系数bj之间不相关。若某个回归系数不显著,可将该项直接剔除,不需重新计算其它回归系数。

例7.1进行试验后得结果如下:试验号

x1(z1)x2(z2)x3(z3)试验结果y12345678

1(95)1(40)1(65)2.111-1(45)2.31-1(20)13.31-1-14.0-1(75)115.0-11-15.6-1-116.9-1-1-17.815b0=Σyi/n=(2.1+…+7.8)/8=37/8=4.625试验号

x1(z1)x2(z2)x3(z3)试验结果y12345678

1(95)1(40)1(65)2.111-1(45)2.31-1(20)13.31-1-14.0-1(75)115.0-11-15.6-1-116.9-1-1-17.8b1=Σxi1yi/n=[12.1+…+(-1)7.8]/8=-13.6/8=-1.7同理得b2=-0.875，b3=-0.3，则回归方程为16(2)回归方程的检验

平方和当FF(p,n-p-1)时，回归方程显著。对例7.1，lyy=(2.12+…+7.82)-(2.1+…+7.8)2/8=30.275u=8(1.72+0.8752+0.32)=29.965，Q=30.275-29.965=0.31F=(29.965/3)/(0.31/4)=128.88>16.69=F0.01(3,4)所以，回归方程极显著。17(3)回归系数的检验

当FjF(1,n-p-1)时，回归系数bj显著。对例7.1：>21.2=F0.01(1,4)>21.2=F0.01(1,4)>7.71=F0.01(1,4)回归系数都显著。185.零水平的重复试验

一次回归方程显著，也不能保证一次回归模型是最好的。为此在零水平(z01,z02,…,z0p)安排一些重复试验，如安排m0次重复试验，其试验结果为y01,y02,…,y0m0。

模型检验(失拟检验)，即检验：零水平试验结果的均值与回归方程中常数项b0是否有显著差异。记若则在水平下，试验中心一次回归与实测值拟合较好。19若则在水平下，试验中心一次回归与实测值拟合较好。否则一次回归拟合不好，需要更高次模型。

对例7.1，在零水平重复2次试验，m0=2,y01=4.5,y02=4.3，S0=0.02，Q=0.31，b0=4.625，n=8，p=3

|t|=1.1078<2.571=t/2(5)，所以，一次回归方程与实测值还是拟合得较好的，此回归方程是适合的。

20

对例7.1，在零水平重复2次试验，m0=2,y01=4.5,y02=4.3，S0=0.02，Q=0.31，b0=4.625，n=8，p=3

|t|=1.1078<2.571=t/2(5)，所以，一次回归方程与实测值还是拟合得较好的，此回归方程是适合的。

将编码公式

x1=(z1-85)/10，x2=(z2-30)/10，x3=(z3-55)/10，代入回归方程

得原始变量的回归方程

216.一次回归正交设计的旋转性

设用正交设计所求得的回归方程为那么回归值的方差

而是p维编码空间内的一个球面，球心在原点，半径为。

位于同一球面上的点的预测值y的方差是相等的。这个性质称为旋转性。22

预测值的方差相等，可以直接比较各预测值的好坏，从而找出预测值相对较优的区域。

在古典回归分析中，由于预测值的方差强烈地依赖于点的位置，所以就不能这样做。今后我们还要专门研究具有旋转性的设计——旋转设计。

显然用回归方程作预报时,预报值y的方差愈小,预报就愈准确,方差愈大,预报值的误差就愈大。有了旋转性以后，这种预报的误差程度可以简单地用

来表示,

大(即点离球心远)误差就大,

小误差也就小。

23

二次回归的正交设计要比一次回归的正交设计复杂一些，p个变量的二次回归方程§7.3二次正交回归设计

为计算二次回归方程的系数，每个变量的水平数应大于3，因而做的试验次数往往比较多。如当p=4时，三水平全因子试验次数是81次，它比4个变量的二次回归的系数C24+2=15要多4倍以上，以致剩余自由度过大。

共有回归系数q=1+C1p+C2p+C1p=C2p+2个，为了得到二次回归方程，试验次数n当然应不小于q。

24

采用“组合设计”。所谓组合设计，就是选择几类具有不同特点的点组合起来形成试验计划。以p=2为例，说明组合设计中试验点在因子空间中的分布。

在二个变量x1,x2场合下，组合设计由n=9个点组成:x1x211-1-1(0,-)(-,0)(,0)(0,)(1,1)(1,-1)(-1,-1)(-1,1)(0,0)25试验号

x1

x21(11)2(1-1)3(-11)4(-1-1)5(0)6(-0)7(0)8(0-)9(00)同理，可以写出p=3的试验计划。组成两水平(+1和-1)全因子试验22分布在x1和x2轴上的星号位置由x1和x2的零水平组成的中心试验点

26一般p个变量的组合设计由下列n个点组成：n=mc+2p+m02p—分布在p个坐标轴上的星号位置点，它们与中心的距离称为星号臂，是待定参数。根据一定的要求调节

，可得到各种具有很好性质的设计（如正交设计、旋转设计）。

m0—在各变量都取零水平的中心点的重复试验次数。它可以只做一次，也可以重复二次或多次。

其中：mc—二水平(+1和-1)的全因子试验的试验点个数2p

，或部分实施的试验点个数2p-1，2p-2等。

27组合设计的优点：

试验点比全因子试验要少得多，但却仍保持足够的剩余自由度(见表7-12）；它在一次回归基础上获得，如果一次回归不显著，只要在一次回归试验的基础上，再在星号点和中心点补充做一些试验，就可求得二次回归方程

要使组合设计成为正交设计，还要确定适当的星号臂。如p=2时，二次回归组合设计的设计矩阵为

28

要使组合设计成为正交设计，还要确定适当的星号臂。如p=2时，二次回归组合设计的设计矩阵为

试验号

x0x1

x2x1x2

x12

x221111111211-1-11131-11-11141-1-111151

002061-0020710002810-0029000000

对一次计划加入星号点后，并不破坏一次变量和交互效应的正交性。只是被x0和xj2破坏了。

29

为使组合设计具有正交性，必须在使得矩阵C=(XX)-1为对角阵的条件下定出的值。可以证明：当mc=2P时（全因子试验情况），有

4+2P2-2P-1(p+0.5m0)=0当mc=2P-1时（部分实施情况），有

4+2P-12-2P-2(p+0.5m0)=0

对于给定的p和m0，就可计算2的值。常用的2值已计算在表7-14上。然后对平方列进行中心化即可。

30

对于给定的p和m0，就可计算2的值。常用的2值已计算在表7-14上。然后对平方列进行中心化即可。

表7-142值表m0p2345(1/2实施)123456789101.0001.4762.0002.391.1601.6502.1982.581.3171.8312.3902.771.4752.0002.5802.951.6062.1642.7703.141.7422.3252.9503.311.8732.4813.1403.492.0002.6333.3103.662.1232.7823.4903.832.2432.9283.6604.0031

运用组合设计进行二次回归正交设计的步骤和统计分析分述如下：§7.4二次正交回归设计的统计分析1.确定因子的变化范围变化区间为

如某研究中有p个因子z1,z2,…，zp，其第j个因子zj变化的下界和上界分别为z1j和z2j(j=1,2,…,p)，零水平为32

例2

玉米高产栽培试验，选择三个主要因子密度z1：上下界z21=4353株/亩，z11=1647株/亩；施N量z2：z22=16.75公斤/亩，z12=3.25公斤/亩；施P量z3：z23=8.4公斤/亩，z13=1.6公斤/亩。采用二次回归正交设计，取m0=3，查表得=1.353，零水平与变动区间分别为：同理得：z02=10，z03=5同理得：2=2，3=2.5332.编制因子水平编码表与一次回归设计类似，对因子作如下线性变换：则有因子水平编码表xj因子z1z2zp10-1-z21

z22……z2pz01+1

z02+2……z0p+pz01

z02……z0pz01-1

z02-2……z0p-pz11

z12……z1p变化区间j1

2……p34xj因子密度z1(株/亩)施N量z2(公斤/亩)施P量z3(公斤/亩)1.35310-1-1.3534353400030002000164716.75151053.258.47.552.51.6变化区间j100052.5对例2

z21=4353，z11=1647；z22=16.75，z12=3.25；z23=8.4，z13=1.6；=1.353，z01=3000，z02=10，

z03=51=1000，2=2，3=2.5因子水平编码表35对例7.2，试验方案与试验结果如下：3.选择相应的组合设计

进行n=mc+2p+m0次试验。如对于p=3,mc=23,m0=3的情况，由x1、x2、x3所占的列组成试验计划。xj因子密度z1(株/亩)施N量z2(公斤/亩)施P量z3(公斤/亩)1.35310-1-1.3534353400030002000164716.75151053.258.47.552.51.6变化区间j100052.5因子水平编码表36x1(z1)x2(z2)x3(z3)y-1(2000)-1(5)-1(2.5)1136-1-11(7.5)1088-11(15)-11128-11111971(4000)-1-112811-11134911-110751111275xj密度z1(株/亩)施N量z2(公斤/亩)施P量z3(公斤/亩)1.35310-1-1.3534353400030002000164716.75151053.258.47.552.51.6mc=23=8个试验点。“对号入座”37x1(z1)x2(z2)x3(z3)y-1.353(1647)0(10)0(5)11801.353(4353)0013440(3000)-1.353(3.25)0120101.353(16.75)0116500-1.353(1.6)1102001.353(8.4)1218000123600012060001221xj密度z1(株/亩)施N量z2(公斤/亩)施P量z3(公斤/亩)1.35310-1-1.3534353400030002000164716.75151053.258.47.552.51.62p+m0=23+3

=9个试验点。38x1(z1)x2(z2)x3(z3)y-1(2000)-1(5)-1(2.5)1136-1-11(7.5)1088-11(15)-11128-11111971(4000)-1-112811-11134911-110751111275-1.353(1647)0(10)0(5)11801.353(4353)0013440(3000)-1.353(3.25)0120101.353(16.75)0116500-1.353(1.6)1102001.353(8.4)