［转载］均匀设计中有重复试验的统计分析

Statistacal Analysis With Repeated Trials
in Uniform Design

SUN Shang-gong
（Beijing Medicai University,100083）

The weighted regression analysis is listed briefly in paper,and the reciprocal of trial variance is taken as the weight.
　　The theoretical optimal coefficient of determination in data and the degree of lack of fit on the regression function could be estimated.Some problems of SASsoftware are pointed out by an example.Author does not think it is good that using average value of the repeated trial results analyse data and the stepwise procedure of regression analysis select variables if size of variables is more than sample size.
Keywords:Weighted regression,coefficient of determination,stepwise procedure of regression analysis,SAS statistical software.

一、加权回归分析
　　设有p个自变量x₁，x₂，…，x_p；因变量为y，则一般的线性回归是：

y=β₀+β₁x₁+β₂x₂+…+β_px_p=X′β+ε　(1)

其中x=(1，x₁，…x_p)′，β=(β₀，β₁，…，β_p)′，为未知参数向量，ε是随机误差，一般叫残差。假设我们有m组的均匀设计，第i组的试验条件为x_(i)=(1，x_i1，…，x_ip)′，它的n_i次重复试验结果为y_i1，…，y_ini，则(1)变成

y_il=x′_(i)β+ε_il　　　　i=1,…,m　　　　l=1,…,n_i(2)

　　记n=n₁+…+n_m，X及Y分别为n×(p+1)及n×1的阵，它们的前n₁行为第1组试验的条件及结果，…，最后n_m行是第m组试验的条件及结果。这时(2)可写成

Y=Xβ+ε　(3)

这里的ε是由ε_il组成的n×1向量。一般的回归中常假定(2)中的ε_il具有独立性及N(0，σ²)的分布，其中σ²与试验号无关，称为方差齐性。而在均匀设计中，由于X布满空间，不同的试验条件往往有很不相同的随机波动，因此方差齐性的假定不易被满足。因此，我们假定(2)中的ε_il满足独立性及

ε_il～N(0，σ²_i)　(4)

在(4)条件下的回归是最重要的也最简单的加权回归分析。令n×n的权阵为

W=diag(w₁,…,w₁,…,w_m,…,w_m)σ²

其中W的对角线对应于(4)中n个方差的倒数，即w_i=1／σ²_i，而σ²是参数。记β的估计为b=(b₀，b₁,…,b_p)′，y_il的估计为_il，这时y的估计即为=Xb，求b的方法是使下述加权残差平方和SS_e为最小

　(5)

令

Z=W^1／2Y,　　　　Q=W^1／2X,　　　　f=W^1／2ε　(6)

这时(3)变为

Z=Qβ+f　(7)

这里f是残差向量，它的每一个分量的分布为N(0，σ²)，也就是说，虽然(3)是不等方差，但经(6)变换后的(7)，则是方差齐性的回归分析，且(3)与(7)有相同的回归系数及相同的残差平方和。这时可得

(X′WX)b=X′WY　(8)

由(8)求得b，即得SS_e为

SS_e=Y′WY-(Y′WX)，df_e=n-p-1　(9)

记加权平均_w为

　(10)

记y的加权离差平方和为SS_T，则

SS_T=Y′WY-SS_bo,　　df_T=n-1　(11)

其中

　(12)

为p个变量引起的加权回归平方和，自由度为p，显然，

SS_T=SS_e+SS_reg

检验p个变量有否线性回归关系的是

　(13)

定义模型(1)的决定系数D(即复相关系数R的平方)为

　(14)

由(8)得b的方差为

V(b)=(Q′Q)^-1σ²=σ²(X′WX)^-1

σ²是(7)中残差f的方差，它的估计为

=SS_e／(n-p-1)　(15)

由公式(6)，可以看出，在y₀=x₀′b+ε₀中，ε₀的残差方差为

σ²₀=σ²／w₀(16)

这w₀是对应于x₀条件下的权，σ由(15)估计。记　

　(17)

这SS_pe是试验的纯误差(pure error)平方和，它的值与回归模型无关。显然，有

(18)

记上式中后一项为SS_lack，自由度为m-p-l，则(18)为

SS_e=SSpe+SS_lack

它充分地说明：如果所选模型不正确，或者讲，模型的回归拟合不充分，则SS_lack显著地不为零，这时残差方差的估计永远有。另一方面，我们也可以看到，如果我们全部的重复试验y_il都用_i取代，这时(18)中第一项不存在，而后一项中由于_i是_i的估计，所以(18)将变得很小，这显然很易造成有很高的复相关系数。但这个残差方差并不能反映一次试验的波动，也不能做下面的是否有“拟合不足”的检验：

　(19)

可以看出，n=m或pm-l时，不能作“拟合不足”的检验。另一个反映模型好坏的指标为

　(20)

这是数据中的理论决定系数，即回旭模型的决定系数的上界是D_max。我们可以用D／D_max作为衡量已有回归模型与“最佳”模型的差异程度。

二、例及问题
　　例［1］(P.32)：为了研究环境污染对人体健康的危害，考察6种金属含量：x₁(Cd),x₂(Cu),x₃(Zn),x₄(Ni),x₅(Cr),x₆(Pb),y是老鼠身上某种细胞的死亡率，先用U₁₇(17¹⁶)均匀表，每一试验条件下重复3次，这里m=17,n_i=3,n=51,数据见［1］中32页。这17组的y值的样本方差s²_i的最小值为0.0837，最大值为8.0882，相差近100倍，方差齐性检验时P=0.001，因此，该试验应当使用加权回归分析。如何估计理论方差σ²_i?由于s²_i随试验条件不同而有差异，因此很自然的想法以试验条件x₁，…，x₆为自变量以s²_i为因变量做回归分析。我们采用二次项法，但这时未知参数共28个，远远超过试验条件组数17。对28个未知参数用逐步回归是不能允许的。原因将在讨论中叙述。此处我们先用单变量法选部分变量：计算s_i及s²_i对27个变量(线性及二次项)的相关系数，从中选出相关最大的15个变量：x₁～x₆，x²₂，x²₃，x²₅，x²₆，x₁x₅，x₂x₃，x₂x₅，x₃x₆，x₄x₆。再在这15个变量上对y做后退法回归，得σ_i的预报公式为

=0.3346+0.0456x₁-0.1565x₂+0.0357x₄+0.2530x₆+0.0090x²₂
-0.0149x²₆+0.0069x₂x₅-0.0046x₃x₆-0.0046x₄x₆　(21)

入选的变量全部p0.05，上公式R=0.9728，残差标准差0.2247。用上述公式的²_i估计σ²_i及求对应的17组数w_i,结果如下：

10.12562，3.89230，1.31365，0.63019，0.85739，2.95520，6.5250，0.77804，
22.21846，7.74011，5.52281，0.12801，0.44826，0.89733，0.76944，
2.59852，2.58384

由公式(21)，我们可以求出在任何试验条件下的权w。
　　下面做加权回归分析：
　　1.使用SAS统计软件中的RS-regression analysis.
不考虑二次型时，SAS软件中RS回归的语法如下：
　　Proc rsreg;
　　Model y=xl-x6／lackfit／covar=6;
　　Weight w;
　　输出结果见Table l.由Table l可见，6个变量全部显著。
由公式(20)得

D_max=0.99933,R²／D_max=82.89%

后一值说明线性模型拟合y不是很充分，而由Table l中的“Lack ofFit”的F值及prob值也可看出线性回归的拟合是显著地“拟合不足”。因此考察二次型回归是必要的。语法是：
　　Proc rsreg;
　　Modely=xl-x6／lackfit／nocode;
Weight w;
输出结果中，R²=0.9993，已达D_max值。但我们发现“Lack of fit”的平方和是负值，且自由度df₁=-1!原因是公式(8)左边的系数阵明显是“病态”(rank(X)=17)。也说是说，从理论上看，(8)的解有无限多组。还有一个怪现象：选中的17个变量的p值全为1.0，即没有一个变量是有统计意义的，这结果是很难被接受的。看来RS回归的软件编写中存在有不足之处。
　　2.二次型回归。
　　为了使(8)中不出现“病态”，我们必须先删去一部分变量，确保(8)中的未知参数不超过17个。一个简单方法是计算y与27个x的加权简单相关系数r_iy，选出绝对值最大的16个变量，结果是：

x₁-x₃，x²₁-x²₄，x₁x₂，x₁x₃，x₁x₄，x₁x₆，x₂x₃，x₂x₄，x₂x₆，x₃x₄，x₃x₆

由这16个变量采用一般加权回归的后退法。这时选中的16个变量全部显著，且R²=0.9993，也达到D_max值。RS回归的优点是可自动计算二次型的极值，而一般回归就无此功能，因此我们只好用另外法(比如MATLAB)计算二次型方阵A的特征根，其结果为-2.84，1.54，-0.20，-2.50，-3.60，4.50。可见，特征根有负有正，说明这二次型y没有极值。它的最大(小)值在边界上。

三、讨论
　　1.试验有重复很有好处：它可以估计试验误差，也可以估计模型是否有“拟合不足”。在重复试验中，重复数不宜多，我以为以3次为宜。如想多做试验，意见以增加试验组数(m)为佳。
　　2.用重复试验的均数取代试验结果的做法不妥；另外，在做变量筛选中，当“自变量总数超过试验组数”时，不应采用经典的逐步回归选变量法。这虽然是国内外很流行的方法，但不妥。［5］中介绍，国外早就用Monte-Carlo法证明：上述方法很易把毫无意义的变量(比如白噪声)引进方程，在［4］(p.139)也介绍，国外用专业知识或理论上早就判定对y无意义的自变量，在人为地改造增大它的数目后，也可以被逐步回归引进公式。这一类研究已经不少了，我们不应再用了。
　　正如本文例子一样，在R大体不变的情形下，二次型(或其它回归)的公式不是唯一的。我们应当尽可能地结合专业知识，找一个在专业上可取的公式。另外，我们不应过份地追求R²要近似于1。本文例中，交叉积平方和为1152.84，它只占回归平方和56078的2%，虽然它的统计检验是p=0.0000，但这么小的平方和在理解公式及求最大(小)值时会带来很大的不便。我们为什么不可以去掉一切交叉积项，而只保留线性及平方项呢?去除交叉积项后，上例的结果，可求得
　　y=53.89-0.23(x₁-5.6)²-0.19(x₂-12)²-0.06(x₃-9.92)²+0.08(x₄-3.69)²+
　　0.15(x₅-5.6)²-0.05(x₆-9.7)²由此式求最大(小)值自然非常方便，也易于解析每一个变量的作用。
　　3.权的影响问题。已经有不少文章指出，在统计检验中，回归分析中正态分布的条件是不重要的，相对说来，等方差性的要求却重要得多［6］。在本文例子中，权(方差的倒数)的影响是很大的：当不使用权时(即w_i=l)，它与y的线性回归中仅能选中3个变量：x₁，x₂，x₄，R²=0.7026；它与加权的Table l的结果差别很大。而笔者用另法求权，即权w_i改成下列值时：
　　8.6100，2.3411，3.1404，0.4615，0.8770，3.1086，1.0589，1.2117，5.7569，3.3398，
　　8.2797，0.1384，0.4518，0.9379，1.1971，2.2259，2.5600
则6个线性项中选中的是x₁～x₅，R²=0.8106。这说明权的大小对回归分析的结果影响很大。求权的方法不是唯一的，但本文方法应是合理的。精细的求法，可用“迭代法”：从目前已求得的回归公式可以求出51个_il，由此可以求值的组内方差s²_i(1)，用s²₁(1)去代替前面的σ²_i，再求权，记为W(1)从而又得y新的回归公式；一直往下做，直到σ²_i的估计值(或W)相当稳定为止。当然这个过程是否收敛是有待证明的。这是方开泰教授建议的方法。
　　4.从本文的例子可见，SAS软件中的RS回归程序是有缺点的，但它可以自动求极值，为保持优点，改变缺点，改造RS回归看来是必要的。
　　Table l.RS regression for xl-x6(linear Model)in SAS software

　　proc rsreg;
　　model y=xl-x6／lackfit／covar=6;weight w;
　　Output:
　　　　　　Response Surtace for Variable Y
　　　　　　Response Mean　　　　35.4534466
　　　　　　Root MSE　　　　　　14.795262
　　　　　　R-Square　　　　　　0.8284
　　　　　　Coef.of Variation　　　41.7315