第四章方差分量线性回归模型

1、第四章方差分量线性回归模型本章考虑的线性模型不仅有固定效应、随机误差，而且有随机效应。我们先从随机效应角度理解回归概念，导出方差分量模型，然后研究模型三种主要解法。最后本章介绍关于方差分量模型的两个前沿研究成果，是作者近期在应用数学学报与国际数学杂志Communications in Statistics上发表的。第一节随机效应与方差分量模型一、随机效应回归模型前面所介绍的回归模型不仅都是线性的，而且自变量看作是固定效应。我们从资料对出发建立回归模型，过去一直是把Y看作随机的，X1，Xp看作非随机的。但是实际上，自变量也经常是随机的，而并不是我们可以事先设计好的设计矩阵。我们把自变量也是随机变

2、量的回归模型称为随机效应回归模型。究竟一个回归模型的自变量是随机的还是非随机的，要视具体情况而定。比如一般情况下消费函数可写为（）这里X是居民收入，T是税收，C0是生存基本消费，b是待估系数。加上随机扰动项，就是一元线性回归模型（）那么自变量到底是固定效应还是随机效应?那要看你采样情况。如果你是按一定收入的家庭去调查他的消费，那是取设计矩阵，固定效应。如果你是随机抽取一些家庭，不管他收入如何都登记他的收入与消费，那就是随机效应。对于随机效应的回归模型，我们可以从条件期望的角度推导出与最小二乘法则等价的回归函数。我们希望通过X预测Y，也就是要寻找一个函数，当X的观察值为x时，这个预测的误差平均起

3、来应达到最小，即（）这里min是对一切X的可测函数L(X)取极小。由于当（）时，容易证明（）故当时，（）要使上式左边极小，只有取。这个结果告诉我们，预测函数取作条件期望E(Y|X)时，可使预测误差最小。我们还可以证明，此时M(X)=E(Y|X)与Y具有最大相关，即（）这里表示相关系数。这是因为当时，易证,同时,于是等号当且仅当（）时成立，此时L(X)是M(X)的线性函数。(4.1.3)与(4.1.7)表达了的极好性质，我们称（）为Y关于X的回归曲线。上面的L(X)可取一切函数。如果限定L(X)是X的线性函数，即要限定（）这里是对X的一切线性函数取极小，则称满足上式的线性函数为Y关于X的回归直线

4、。我们可以求出的解。记，则（）这里（）（）（）对L(0，)求微分(矩阵微商公式)得：（）解得（）这里当然假定存在，否则使用广义逆。此时的预测误差方差是（）（）为复相关系数。它指出了Y与多元变量之间的线性相关程度，是一元相关系数（）的推广。从条件期望角度我们导出的随机效应回归模型的回归直线表达式，与从最小二乘角度导出的固定效应的回归方程，表达式是等价的，所以从计算角度，我们不怎么区分。二、方差分量模型概念上段我们建立了随机效应概念，将自变量也视作随机变量，这就可以导出方差分量模型。方差分量模型研究工作的奠基人是我国最早的统计学家许宝驭马录先生。还是刚才提到的消费函数回归模型，我们作随机抽样。考虑

5、居民按职业的分类，如工人、教师、医生、律师、店员等等，记为，我们从这些职业中随机抽取了n个样本，则模型可写为（）这里Xi可看作是第i种职业对收入的效应。如果我们事先安排好取哪个职业的，当然Xi是固定效应。可是我们现在对职业选取是随机的，而且我们还想研究职业效应的方差，这就导入了方差分量模型，因为现在Cij的方差由两部分组成：（）为了数学符号统一，我们将经济学中的符号改过来，刚才建立的模型是（）它有一项固定效应，一项随机效应1，一项随机误差。如果还要考虑地区因素对消费的影响，还可以加进第二个随机效应2，于是可得模型（）这次我们省掉了取值的标记，Y的方差由三项组成。一般地，我们建立方差分量模型如下

6、：（）这里有固定效应向量，随机效应向量（）并且将随机误差项也并入了随机效应向量去。设计矩阵X以及（）都是已知的。对于随机效应,合理的假定是（）当然以后有时还可以考虑i是向量的情况，不过这里假定每个i是一维变量。记，（）则方差分量模型可记为（）模型的主要任务是要估计固定效应向量与方差分量。和一般的多元线性回归模型相比，就是待估的方差多了。通过这些介绍，我们就可以方便地将各种经济方面的普通线性回归模型改造成方差分量模型，当然要根据实际。第二节方差分量模型的解法对于方差分量模型（）一般都采用二步估计法，首先估计方差分量，然后再估计固定效应。按照广义最小二乘（）其中（）所以方差分量模型解法的关键是估计

7、方差分量。以下介绍的方法，也都是针对方差分量估计方法而言的。一、方差分析法先从一个简单的模型结合数据结构形象地说明方法。考虑模型（）0为总平均，是固定效应，1,，m是随机效应，。对于随机误差。这个模型如果记作方差分量模型的标准形式是（）其中设计阵X=(1,1,，1)，随机效应矩阵为（）我们手中资料只有我们采用(4.2.4)记法方便一些，将资料Y排成表ji12k组内平均1Y11Y12Y1k2Y21Y22Y2kmYm1Ym2Ymk方差分析主要掌握三点，一是计算组内差、组间差，二是作平方和分解，三是计算各自的自由度。先计算总平均：（）总变差(全体资料与总平均的偏差平方和)：（）各组平均(各组资料横向

8、相加并平均)（）组间差(各组平均数与总平均数的偏差平方和)（）组内差(各组数据与本组平均数的偏差平方和)（）则必有平方和分解（）将各平方和除以各自的自由度。ST有一个约束 (4.2.7)，自由度为；SA有m组差，1个约束，自由度为m1；Se有mk组差，m个约束，自由度为mk-m。注意有自由度分解：（）于是算出均方：（）（）（）因为假定为随机效应，可以算出各均方的均值：（）（）以代者，代替，得方程组：（）解得（）这样就作好了方差分量的估计，然后可以按(4.2.2)作出的估计。因为这里的方差分量是由方差分析法作出的，故称为方差分析法。推广到一般的方差分量模型时，基本原则是类似的。我们不妨考虑方差分

9、量模型（）先对总平方和YY作平方和分解（）其中S是在模型Y=X+中，的回归平方和：（）是在模型中，消去影响后1的平方和（）类似地，是在模型中消去和1影响后，2的平方和：（）最后的S为残差平方和（）可以验证（）（）（）（）这里（）（）（）这里P*表示关于*的投影阵。下面计算各平方和的均值。（）因为，所以上式第一项为0。在第三项中，（）在第六项中（）所以最后有（）其中（）（）（）（）类似还可以求得（）（）（）于是我们得到方程组（）解此方程组，就可以得到的估计。然后进入二步估计的第二步，就可以得到关于固定效应的估计。算例4.2.1 市场收益率与股利和换手率的关系考虑一个随机效应的多元线性模型U的形式

10、如同(4.2.6)。问题的实际背景是，观测对象被分成了m组，可能存在一个随机效应向量对各组资料有不同的作用。模型也可以写作数据结构及具体数值如下表所示，m=6,k=6。这些资料采自96上海股票市场资料总汇。我们研究目的一是看过去一年的股利收入与当年换手率对当年市场收益率有何影响，二是想知道是否存在一个潜在的尚未观测到的随机效应，对行业有明显影响。当然这种情况采用方差分量模型比较合适。要注意本例是两个方差量，上一章第二节模型(3.2.10)也是两个待估的方差量。它们的随机效应作用范围不一样，不是一回事。表4.2.1 1996年股市资料类别股号股名1996年收益率%1995年股利%1996年日换手

11、率商业类628新世界64.769203.12631中百一店46.84511.81.68632华联商厦41.95811.31.81655豫园商城16.19511.21.10682新百公司79.9115.23.36694大连商场91.3885.84.26电子类602真空电子33.112103.52651飞乐音响8.10801.95800天津磁卡271.76353.74839四川长虹381.686604.41850华东计算机14.63813.27870厦华电子68.5793.94.20化工类617联合化纤-21.8710.51.53618卤碱化工22.37022.63672广东化纤11.86004.

12、65688上海石化179.8173.54.38886湖北兴化236.32852.34.45889南京化纤-33.1222.44.644.11医药类664哈医药191.66654.32671天目药业111.135164.11812华北制药152.0158.44.11842中西药业13.8212.61.71849四药股份17.8922.51.68779四川制药-24.744012.28钢铁类608异型钢管8.389102.24665沪昌特钢75.391.874.01674四川峨铁35.93233.64808马钢股份86.5280.54.52845钢管股份-25.17001.61894广钢股份51.

13、3712.77.08机械类604二纺机14.4102.31605轻工机械6.12203.50610中纺机0.70102.27806昆明机床41.8521.14.22841上柴股份66.981202.20862南通股份41.093201.36首先我们作普通最小二乘回归，得到，然后计算。此时的已消除固定效应影响，我们将它排成平面表，以作方差分析，计算与。计算过程从(4.2.7)至(4.2.20)。从下面计算过程可以看到，平方和分解式是满足的：即(ST=SA+S)。对于本例资料，随机误差4055.6远大于随机效应方差，组内差远大于组间差，可以认为随机效应不明显，即行业差别不明显。对于选定的方差分量模

14、型，回归结果是它的标准差很小，为1.0084，这正是采用方差分量模型广义最小二乘意义所在。拟合效果图(图4.2.1.1)令人满意。-方差分量模型方差分析法计算程序, 例 4.2.1.第一列为 Y, 以后各列为 X例421.D 数据文件中, m=6, k=6, p=2要显示原始资料吗? 0=不显示, 1=显示（0）先作普通最小二乘, 并打印结果:现在作线性回归显著性检验, 计算t,F,R 统计量请输入显著性水平a, 通常取a=0.01, 0.05, 0.10, a=? （0.05）-线性回归分析计算结果 样本总数 36 自变量个数 2-回归方程 Y = b0+b1*X1+.+b2*X2 Y =

15、5.3188 + 4.5656 X1 + 6.0128 X2回归系数 b0, b1, b2, ., b2 5.3188 4.5656 6.0128-残差平方和: 149400.60 回归平方和: 131759.30误差方差的估计 : 4150.0170 标准差 = 64.4206-线性回归显着性检验显著性水平 : .050-回归方程整体显著性F检验, H0:b0=b1=.=b2=0 F统计量: 14.5517 F临界值F(2, 33) 3.285全相关系数 R : .6846-回归系数逐一显著性t检验, H0:bi=0, i=1,.,2 t 临界值 t( 33) 1.6924回归系数b1-b

16、2的t值: 5.6318 1.1073-打印方差分析资料 Y(I,J) -50.6210 -39.0275 -38.6721 137.5441 -56.0541 -50.4640 -22.4490 -8.9357 -7.8936 8.0546 37.4223 -20.2416 -25.8348 221.1285 -21.4183 83.6328 -4.9701 -18.2668 -46.8722 75.9169 132.1827 -13.6501 51.7486 6.1371 30.6483 -14.9082 -34.5268 -8.9422 -40.1694 -42.8772 33.9744

17、 20.2007 -77.2975 -103.9000 -8.8455 -81.7529计算各种平均: 总平均YYba: .0000各组平均Yba: -16.216 -2.341 39.045 34.244 -18.463 -36.270计算各种变差: 总变差SST, 组间差SSA, 组内差SSESST= 149400.6000 SSA= 27731.9200 SSE= 121668.7000打印方差分量的估计 SIGMAA= 248.4600 SIGMAE= 4055.6240下面计算协方差阵及其逆的一块,并作分解计算广义最小二乘估计, 模型转换Y=PY, X=PX, 下面打印矩阵 P 的一

18、块, PP=的逆SIG1= .1239 -.0007 -.0007 -.0007 -.0007 -.0007SIG1= .0000 .1238 -.0008 -.0008 -.0008 -.0008SIG1= .0000 .0000 .1237 -.0008 -.0008 -.0008SIG1= .0000 .0000 .0000 .1237 -.0008 -.0008SIG1= .0000 .0000 .0000 .0000 .1236 -.0009SIG1= .0000 .0000 .0000 .0000 .0000 .1235下面打印的是广义最小二乘的统计结果:现在作线性回归显著性检验,

19、 计算t,F,R 统计量请输入显著性水平a, 通常取a=0.01, 0.05, 0.10, a=? (0.05)-线性回归分析计算结果 样本总数 36 自变量个数 2-回归方程 Y = b0+b1*X1+.+b2*X2 Y = .7309 + 4.5682 X1 + 5.7726 X2回归系数 b0, b1, b2, ., b2 .7309 4.5682 5.7726-残差平方和: 2290.69 回归平方和: 2014.61误差方差的估计 : 63.6303 标准差 = 7.9769-线性回归显着性检验显著性水平 : .050-回归方程整体显著性F检验, H0:b0=b1=.=b2=0 F统

20、计量: 14.5113 F临界值F(2, 33) 3.285全相关系数 R : .6841-回归系数逐一显著性t检验, H0:bi=0, i=1,.,2 t 临界值 t( 33) 1.6924回归系数b1-b 2的t值: 5.6411 1.0683-比较残差平方和: 普通最小二乘的: 149400.6000 : 广义最小二乘的: 2290.6910下面打印的是利用广义最小二乘的回归系数去计算原始资料的回归拟合结果:要作回归预测吗? 键入 0=不预测, 1=要预测 (0)回归系数: .7309 4.5682 5.7726要打印拟合数据吗? 0=不打印, 1=打印 (0)计算结束。 -计算中需要一

21、些分析与技巧。因为需要用它及(4.2.2)计算*注意UU是一个分块对角阵，共有m个对角块，每块是一kk方阵，元素皆为1。这样的计算就容易了。恰好X也是分成m块，每块为kp阵。于是而，这就大大简化了计算。如果真的要计算3636的矩阵的逆，一般微机是不可能的。这个程序开始时调用了普通多元线性回归程序对原始资料回归。得到的回归方程为拟合效果也很好(见图4.2.1.2)。但是，这个模型的残差向量的标准差为64.4206，远大于方差分量模型的1.0084。这可以对比说明方差分量模型的作用。二、最小范数二次无偏估计方法方差分量模型中方差分量的最小范数二次无偏估计(Minimum Norm Quadrati

22、c Unbiased Estimator, MINQUE)是提出的，思路类似于他处理奇异广义线性模型的分块逆矩阵法，先提出估计应满足的性质，根据这些性质去求解(一般是一个极小值问题)，若能解出，当然就有那些设定的性质。考虑一般的方差分量模型(4.2.1)，我们要估计方差分量及其线性函数（）首先考虑的估计应具有的形式，因为是估计方差，可考虑采用二次型YAY的形式，即（）A为待估对称矩阵。再来考虑YAY应具有的性质：(1)关于参数的平移不变性。若参数有平移：（）则方差分量的估计应该不变。此时原模型变为（）其二次型估计变为，应该有（）这等价于（）即平移不变性(4.2.50)应满足的充要条件是(2)估

23、计量的无偏性。因为（）今要求对一切2成立：（）则其充要条件应为（）（3)最小范数准则。经过研究，满足平移不变性与无偏性，尚不能唯一确定待估对称矩阵A。于是可以再加一个优良性质。如果随机效应向量i,i=1,，m是已知的，则=c2的估计应该为（）这里（）现在用YAY去估计=c2，在满足不变性条件AX=0时，（）我们自然希望(4.2.56)与(4.2.58)之间相差很小，这只要求矩阵与UAU之间相差很小。我们选用矩阵范数UAU-来度量与UAU之间的差异，即应选A使（）范数可选欧氏范数（）若记V=UU，则因是分块对角阵及无偏性要求：（）是已知的，则极小化范数等价于极小化。总结上述三项优良性要求，求=c

24、2的最小范数无偏估计的问题，归结为求下述极值问题：因为目标函数是矩阵的迹，所以称为最小迹问题。要解决极值问题(L1)的解，可以先对其简化。因为，V正定而存在。令（）则模型(L1)可变为定理4.2.1 极值问题L2的解为（）其中为方程组（）的解，为L(Z)的正交补空间L(Z)上的投影阵：（）证明先证方程组(4.2.64)兼容。设B0满足极值问题(L2)的约束条件，即（）则。因是往上的投影阵，故，记，则（）引进拉直算符，表示将Wj按列拉成一个n21的向量，表示将B0按列拉直成一个n21的向量。定义n2m矩阵G， （）L(G)表示按G的各列展成的子空间，是n2维。定义1为向量在L(G)上的投影（）2

25、为向量在L(G)上的投影（）具体形式,则（）这个情形与最小二乘的基本原理图标(图1.2.1)一样，那里是将向量Y向子空间L(X)投影，于是存在常数1，p,使Y=X=X11+Xpp，即将Y表为X列向量的线性组合。现在是必存在常数，可将表示为G的列向量的线性组合，当然也可以将在L(G)的投影1表成线性组合：（）于是（）现在看(4.2.67)，由公式得（）这就证明了是方程组(4.2.64)的一组解。再证明(4.2.63)的B*是问题(L2)的解。由于故B*满足的约束条件是从方程组tr(BWi) =Ci,i=1,，m中解出来的，当然满足这个约束条件。余下看B*是否是trB2的极小值解。设B为任一满足(

26、L2)约束条件的解，记D=BB*，则D对称，DZ=0,tr(DWj)=0,于是（）因此（）这就证明了B*是极小值解。现在回到原问题(L1)。由于(4.2.62)所定义的三个变换都是可逆变换，故(L1)与(L2)等价。于是(L1)的解存在，（）对于最小范数二次无偏估计，我们介绍了它的原理、算法、解的存在性。深入的讨论还可以引出一些问题，如解的非负性，计算的复杂性等，我们这里不再讨论了。三、极大似然法在假定方差分量模型随机效应服从正态的情况下，可以使用极大似然法求出参数估计。设模型为（）（）其它记号仍同以前。则（）Y的条件密度为（）取对数，略去常数，得似然函数方程（）用矩阵微商公式（）（）（）得似

27、然方程组（）由于似然方程组没有显式解，统计学家提出了一些迭代算法。但是在数值计算技术发达的今天似乎没有必要再介绍这些迭代算法，甚至连求导的似然方程组都没有必要。对于样本，以及已知设计阵，代入到似然函数(4.2.82)中，调用本软件所附的计算极值的程序(由Sargent改进的Powell算法)，就可以计算出，2的估计。第三节方差分量模型参数的广义岭估计在这一节我们先将方差分量模型的方差分量化为派生模型的均值参数，分别作出其相对于LSE和BLUE的广义岭估计，再根据二步估计法作出原模型均值参数的广义二乘估计及其进一步的岭估计，证明了这样不仅使方差分量估计的均方误差减少，而且使原模型均值参数估计的均

28、方误差也不增加和进一步减少。我们还找到了岭参数仅仅依赖于样本的估计。这样既将岭估计方法推进至方差分量模型，也改进了方差分量模型参数的离差均值对应方法。一、方差分量岭估计的构造与性质岭估计和Stein估计是减少均方误差的行之有效的方法。已有文献将岭估计推广到多元线性回归模型和设计阵是列降秩的情况。本节则试图将它们推广到方差分量模型，以期改进离差均值对应方法。设有一般方差分量模型（）这里为已知设计阵，为固定效应向量，为随机效应向量。假定记则模型(4.3.1)可记为（）按照两步估计法我们先求方差分量的估计，为此采取离差均值对应。记，于是有模型（）再记，因为= ，可得派生模型（）此时原模型的方差分量成

29、了派生模型的均值参数。对派生模型(4.3.4)用最小二乘法求出均值参数的最小二乘估计(LSE)：（）是的无偏估计，但由于派生模型是广义线性模型，它未必是的最佳线性无偏估计(BLUE)。当F已知时，(4.3.6)只是在一些特殊情形研究证明了。不难举例说明可能出现负值分量，也不难看出DD很容易呈现病态。当X，V1，Vm中有两个接近相等时，DD的列向量就呈复共线。因此需要对加以改进。对亦然。考虑原模型中心化与派生模型(4.3.4)中心化的关系，有引理若模型(4.3.1)设计矩阵X，U1，Um均已中心化，则派生模型(4.3.4)也已中心化。证明已知这里,则故证毕我们不妨设D已经中心化、标准化。当F未知

30、时，对模型(4.3.4)定义的关于的广义岭估计为（）这里K=diag(k1,km)，ki0，i=1,m。P是正交方阵致P(DD)P=diag(1,m)=。当k1=k2=km时，就成为狭义岭估计。我们将看到狭义岭估计虽然能改进的均方误差，但不一定能保证第二步估计时的均方误差不增加。(K)的第i个分量记作。于是（）当F已知时，定义的关于的广义岭估计为 (4.3.9)（）这里P(DF-1D)P=。自然这里K、P、与(4.3.7)不同。引理4.3.2 若、存在，则证明仅证后式，由定义，注意K0,证毕引理说明、分别是对和的压缩。但这种压缩是模长的压缩，有的分量可能还有伸长。能否使各分量都压缩呢?我们有引

31、理设存在，0t0使(4.3.13)成立，同时有（）或使(4.3.14)成立，同时有（）证明因为，故由引理及定理，存在K0，使(4.3.13)成立，且有。于是即(4.3.24)成立。同理可证(4.3.14)、(4.3.25)可以同时成立。证毕如果我们找不到引理的方法，得不到定理的结论，那么我们的两步估计法对岭估计而言就将失去意义。如果发生了，那么再进一步对采用岭估计我们就无法判断最终的估计的均方误差到底是减少还是增大了。现在有了定理，我们就可以对广义线性模型(4.3.2)继续采用广义岭估计：(4.3.26)这里。是正交方阵致。则又存在，使均方误差进一步缩小：（）且存在，（）利用这种形式，可以作出不依赖于未知参数的的估计。对亦然，不再赘述。第四节方差分量模型参数经验Bayes估计在这一节我们将构造误差正态的方差分量模型参数的经验Bayes估计，定义逆拉直运算，利用作者自己曾经证明的多参数指数族EB估计收敛结果，证明所构造的经验Bayes估计的收敛性。要利用的多参数指数族参数经验Bayes估计收敛速度的结果概述如下：设t=(t1,tp)以及t (1),t(n)为来自多元密度f(t)的当前样本和历史样本，它们都是独立同分布的，f(t)和它的偏导数f