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

第四章方差分量线性回归模型本章考虑的线性模型不仅有固定效应、随机误差，而且有随机效应。我们先从随机效应角度理解回归概念，导出方差分量模型，然后研究模型三种主要解法。最后本章介绍关于方差分量模型的两个前沿研究成果，是作者近期在应用数学学报与国际数学杂志CommunicationsinStatistics上发表的。第一节随机效应与方差分量模型一、随机效应回归模型前面所介绍的回归模型不仅都是线性的，而且自变量看作是固定效应。我们从资料对Yi,Xii,Xpi；出发建立回归模型，过去一直是把Y看作随机的，Xi，Xp看作非随机的。但是实际上，自变量也经常是随机的，而并不是我们可以事先设计好的设计矩阵。我们把自变量也是随机变量的回归模型称为随机效应回归模型。究竟一个回归模型的自变量是随机的还是非随机的，要视具体情况而定。比如一般情况下消费函数可写为CC0b(XT)()这里X是居民收入，T是税收，C0是生存基本消费，b是待估系数。加上随机扰动项，就是元线性回归模型CC0b(XT)()那么自变量到底是固定效应还是随机效应?那要看你采样情况。如果你是按一定收入的家庭去调查他的消费，那是取设计矩阵，固定效应。如果你是随机抽取一些家庭，不管他收入如何都登记他的收入与消费，那就是随机效应。对于随机效应的回归模型，我们可以从条件期望的角度推导出与最小二乘法则等价的回归函数。我们希望通过X预测Y,也就是要寻找一个函数YM(X)M(Xi,Xp)，当X的观察值为x时，这个预测的误差平均起来应达到最小，即EYM(X)2minEYL(X)2()这里min是对一切X的可测函数L(X)取极小。由于当M(X)E(Y|X)()时，容易证明EYM(X)M(X)L(X)0()故当M(X)E(Y|X)时，222EYL(X)2EYM(X)2EM(X)L(X)2()要使上式左边极小，只有取L(X)M(X)E(Y|X)。这个结果告诉我们，预测函数取作条件期望E(Y|X)时，可使预测误差最小。我们还可以证明，此时M(X)=E(Y|X)与Y具有最大相关，即(Y,M(X)mLax(Y,L(X)()这里P表示相关系数。这是因为当M(X)E(Y|X)时，易证Cov(Y,L(X)Cov(M(X),L(X),同时Cov(Y,M(X)Cov(M(X),M(X)，于是222(YL(X)Cov(丫丄(X)Cov(M,L(X)，D(Y)DL(X)D(Y)DL(X)Cov2(M(X),L(X)DM(X)DM(X)DM(X)DL(X)D(Y)DM(X)22(M(X),L(X)(Y,M(X)2(Y,M(X)等号当且仅当|(M(X),L(X)|1()时成立，此时L(X)是M(X)的线性函数。(4.1.3)与(4.1.7)表达了M(X)E(Y|X)的极好性质，我们称YM(X)E(Y|X)()为Y关于X的回归曲线。上面的L(X)可取一切函数。如果限定L(X)是X的线性函数，即要限定E|Y(01X1mXm)|2mjn()L我们可以求出0,1,m的解。记(1：,m)，则L(0,)E|Y(01X1mXm)|2b2Rxx2RxyD(Y)()这里这里min是对X的一切线性函数取极小，则称满足上式的线性函数为Y关于X的回归直线。bE(Y)(01EX1mEXm)RxxE(XEX)(XEX)DX1Cov(X1,X2)COV(Xi,Xm)COV(Xm,Xi)C0V(Xm,X2)D(Xm)Rxy(Cov(Y,XJ,Cov(Y,Xm)对L（卩0,3）求微分（矩阵微商公式（XAX）2AX）得:XRxxRxy)解得EY?E(X)RxxRxy)这里当然假定RxX存在，否则使用广义逆。此时的预测误差方差是L(?0,?)E|Y(0,mXm)|2XY2?RxyD(Y)RxyRxxRxy)1(RxyRxxRxy)2/y)为复相关系数。它指出了Y与多元变量XX1,Xm之间的线性相关程度，是一元相关系数Cov(X,Y)、DXDY的推广。)与从最小二乘角度导出从条件期望角度我们导出的随机效应回归模型的回归直线表达式,的固定效应的回归方程，表达式是等价的，所以从计算角度，我们不怎么区分。二、方差分量模型概念上段我们建立了随机效应概念，将自变量也视作随机变量，这就可以导出方差分量模型。方差分量模型研究工作的奠基人是我国最早的统计学家许宝驭马先录生。还是刚才提到的消费函数回归模型，我们作随机抽样。考虑居民按职业的分类，如工人、教师、医生、律师、店员等等，记为Xi,i1,m，我们从这些职业中随机抽取了n个样本，则模型可写为CijC0b（XiTi）ij,j1,n,i1,m（）这里Xi可看作是第i种职业对收入的效应。如果我们事先安排好取哪个职业的，当然Xi是固定效应。可是我们现在对职业选取是随机的，而且我们还想研究职业效应的方差，这就导入了方差分量模型，因为现在Cij的方差由两部分组成：2222Var（Cij）02b2X22（）为了数学符号统一，我们将经济学中的符号改过来，刚才建立的模型是YijU11iij,i1,m,j1,n（）它有一项固定效应口，一项随机效应E1，一项随机误差&。如果还要考虑地区因素对消费的影响，还可以加进第二个随机效应E2，于是可得模型YU11U22（）这次我们省掉了取值的标记，Y的方差由三项组成。一般地，我们建立方差分量模型如下：YXU11Umm（）这里有固定效应向量卩，随机效应向量（1,2,m）（）并且将随机误差项&也并入了随机效应向量去。设计矩阵X以及U(U1J2,Um)()都是已知的。对于随机效应i,i1,m，合理的假定是E(i)0,Cov(i,j)0,ij()D(i)2/i,i1,m当然以后有时还可以考虑Ei是向量的情况，不过这里假疋每个Ei是一维变量。记ViUiUi,i1,2*2X/,m,1V1mVm,()则方差分量模型可记为E(Y)X,Var(Y)()222模型的主要任务是要估计固定效应向量卩与方差分量1,2,m。和一般的多元线性回归模型相比，就是待估的方差多了。通过这些介绍，我们就可以方便地将各种经济方面的普通线性回归模型改造成方差分量模型，当然要根据实际。第二节方差分量模型的解法对于方差分量模型Y1XU11Ummnpp1nP1P11nPmPm1m2E(Y)X,Var(Y)2口_一般都采用二步估计法，首先估计方差分量i2,m，然后再估计固定效应卩。按照广义最小二乘*（X?1X）X?1Y（）其中m?i2UiUi（）i1所以方差分量模型解法的关键是估计方差分量。以下介绍的方法，也都是针对方差分量估计方法而言的。、方差分析法先从一个简单的模型结合数据结构形象地说明方法。考虑模型)Yij0iij,i1,m,j1,n卩0为总平均，是固定效应,E1,，Em是随机效应,Ei0,Cov(i,j)0,ij，Var(i)A,i1,m。对于随机误差ij,Var(ij)22。这个模型如果记作方差分量模型的标准形式是)其中设计阵X=（1,1,，1），随机效应矩阵为010k01(426)U01101mkm我们手中资料只有Y（Y11,，丫1k,Y21,，丫2k,，丫mk）我们采用（4.2.4）记法方便一些，将资料Y排成表ji12k组内平均1Y11Y12Y1kY12Y21Y22Y2k丫2mYm1Ym2YmkYm方差分析主要掌握三点，一是计算组内差、组间差，二是作平方和分解，三是计算各自的自由度。先计算总平均：mk丫1Y-IJmki1j1()总变差（全体资料与总平均的偏差平方和）：St各组平均（各组资料横向相加并平均）mk(Yiji1j1Y)2-1kY|j11,mY,ikj组间差（各组平均数与总平均数的偏差平方和)kmSaj1(Yii1Y)组内差（各组数据与本组平均数的偏差平方和)mkSi1(Yjj1Yi)2则必有平方和分解StSaS(429)(4210)(4211)(4212)将各平方和除以各自的自由度。St有一个约束Y.（427），自由度为n1mk1；Sa有m组差，1个约束，自由度为m1；Se有mk组差，m个约束，自由度为I、一、mk-m。注意有自由度分解:fAmk1(m1)(mkm)(4213)于是算出均方:QtQaSTmk1七Sa1Smkm(4214)(4215)(4216)因为假定为随机效应，可以算出各均方的均值:E(QA)kE(Q)以Qa代者E(Qa),Q代替E(Q)，得方程组:)kA22QA解得Qe,(QAQ)/k)这样就作好了方差分量的估计，然后可以按(422)作出卩的估计。因为这里的方差分量是由方差分析法作出的，故称为方差分析法。推广到一般的方差分量模型时，基本原则是类似的。我们不妨考虑方差分量模型YXCov(Y)U11U2212U1U122222U2U22I)先对总平方和YY作平方和分解YYS2)其中Sb是在模型Y=X卩+&中,卩的回归平方和:SES()YX(XX)XY)S1是在模型Y中，消去卩影响后Ei的平方和S1SES(,1)SES()类似地，S2是在模型YU11U22中消去卩和Ei影响后,E2的平方和:EES(,1,2)SES(1)最后的S&为残差平方和)SYYSES(,1,2)可以验证SY(ID)YS1Y(DD1)YS2Y(D1D12)YSYD12Y这里DIX(XX)XIPXD1DDU1(U1DU1)U1DDPDUD12D1D1U2(U2D1U2)U2D1D1PD1U2这里P*表示关于*的投影阵。下面计算各平方和的均值。E(S1)X(DD1)Xtr(DD1)2U1U112U2U2222IX(DD1)Xtr(U1DU1)122tr(U1D1U1)122tr(U2DU2)22tr(U2D1U2)22tr(DD1)2因为DX0,D1X0，所以上式第一项为0。在第三项中，tr(U1D1U1)trU1DU1U1DU1(U1DU1)U1DU10)在第六项中所以最后有其中类似还可以求得于是我们得到方程组解此方程组，就可以得到固定效应的估计。tr(DD1)trD1U1(U1D1U1)U1D1tr(U1D1U1)U1D1U1rk(U1D1U1)rk(U1D1)rk(U1X)rk(X)rk(U1X)rk(X)()E(S1)c112(2c3)2r2()c1tr(U1DU1)()c2tr(U2DU2)()c3tr(U2D1U2)()r1rk(X),r1r2rk(U1|X)()E(S222)c22r3()E(S)(n2r1r2r3)()r3rk(XU1U2)r1r2()S12c11(c222c3)2r2S222c22r3()S(nr1r2r3)2221,22,2的估计。然后进入二步估计的第二步，就可以得到关于算例市场收益率与股利和换手率的关系考虑一个随机效应的多元线性模型YXUn1nPp1nmmlU的形式如同(426)。问题的实际背景是，观测对象被分成了不同的作用。模型也可以写作YjXji数据结构及具体数值如下表所示，我们研究目的一是看过去想知道是否存在一个潜在的尚未观测到的随机效应,差分量模型比较合适。要注意本例是两个方差量，上一章第二节模型机效应作用范围不一样，不是一回事。m组，可能存在一个随机效应向量对各组资料有i1,m,j1,km=6,k=6。这些资料采自96上海股票市场资料总汇年的股利收入与当年换手率对当年市场收益率有何影响，二是对行业有明显影响。当然这种情况采用方(3210)也是两个待估的方差量。它们的随表4211996年股市资料类别股号股名1996年收益率%1995年股利%1996年日换手率商业类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.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.3712.77.08机械类604二纺机14.4102.31605轻工机械6.12203.50610中纺机0.70102.27806昆明机床41.8521.14.22841上柴股份66.981202.20862南通股份41.093201.36首先我们作普通最小二乘回归，得到(xx)1xy，然后计算Y*YXj?。此时的Yj已消除固定效应影响，我们将它排成平面表，以作方差分析，计算A与2。计算过程从(4.2.7)至(4.2.20)。从下面计算过程可以看到，平方和分解式是满足的：SSTSSASSE即(St=Sa+Se)。对于本例资料，随机误差4055.6远大于随机效应方差，组内差远大于组间差，可以认为随机效应不明显，即行业差别不明显。对于选定的方差分量模型，回归结果是Y0.13734.7053X13.3258X2它的标准差很小，为1.0084，这正是采用方差分量模型广义最小二乘意义所在。拟合效果图(图4.2.1.1)令人满意。方差分量模型方差分析法计算程序，例421.第一列为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*X2Y=5.3188+4.5656X1+6.0128X2回归系数b0,b1,b2,.,b2残差平方和:149400.60回归平方和:131759.30误差方差的估计:4150.0170标准差=64.4206线性回归显着性检验显著性水平:.050回归方程整体显著性F检验,H0:b0=b1=.=b2=0F统计量:14.5517F临界值F(2,33)3.285全相关系数R:.6846回归系数逐一显著性t检验,H0:bi=0,i=1,.,2t临界值t(33)1.6924回归系数b1-b2的t值打印方差分析资料Y(I,J)-50.6210-39.0275-38.6721137.5441-56.0541-50.4640-22.4490-8.9357-7.89368.054637.4223-20.2416-25.8348221.1285-21.418383.6328-4.9701-18.2668-46.872275.9169132.1827-13.650151.74866.137130.6483-14.9082-34.5268-8.9422-40.1694-42.877233.974420.2007-77.2975-103.9000-8.8455-81.7529计算各种平均:总平均YYba:.0000各组平均Yba:-16.216-2.34139.04534.244计算各种变差:总变差SST,组间差SSA,组内差SSESST=149400.6000SSA=27731.9200SSE=121668.7000打印方差分量的估计SIGMAA=248.4600SIGMAE=4055.6240面计算协方差阵及其逆的一块,并作分解计算广义最小二乘估计,模型转换Y=PY,X=PX,下面打印矩阵P的一块,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-.0009.1235SIG1=.0000.0000.0000.0000.0000下面打印的是广义最小二乘的统计结果:现在作线性回归显著性检验,计算t,F,R统计量请输入显著性水平a,通常取a=0.01,0.05,0.10,a=?(0.05)线性回归分析计算结果样本总数36自变量个数2回归方程Y=b0+b1*X1+.+b2*X2Y=.7309+4.5682X1+5.7726X2回归系数b0,b1,b2,.,b2残差平方和:2290.69回归平方和:2014.61误差方差的估计:63.6303标准差=7.9769线性回归显着性检验显著性水平:.050回归方程整体显著性F检验,H0:b0=b1=.=b2=0F统计量:14.5113F临界值F(2,33)3.285全相关系数R:.6841回归系数逐一显著性t检验,H0:bi=0,i=1,.,2t临界值t(33)1.6924回归系数b1-b2的t值:比较残差平方和:普通最小二乘的:149400.6000:广义最小二乘的:2290.6910下面打印的是利用广义最小二乘的回归系数去计算原始资料的回归拟合结果要作回归预测吗?键入0=不预测,1=要预测(0)回归系数要打印拟合数据吗?0=不打印,1=打印(0)计算结束。圖421.1T-原始数据-拟合数据计算中需要一些分析与技巧。因为Auu2In需要用它及（422）计算卩*(X?1X)1X?1Ykxk方阵，元素皆为1。这样工注意UU是一个分块对角阵，共有m个对角块，每块是的计算就容易了。恰好X也是分成m块，每块为kxp阵。而?1X?1X(X1,Xm)diag(?1,X1?1X1mA11X1m)Xm2Im，这就大大简化了计算。如果真的要计算36X36的矩阵工的逆，一般微机是不可能的。这个程序开始时调用了普通多元线性回归程序对原始资料回归。得到的回归方程为16.0128X264.4206,远大于方差拟合效果也很好（见图4.2.1.2）。但是，这个模型的残差向量的标准差为分量模型的1.0084。这可以对比说明方差分量模型的作用。圖t一原始数据-拟合数据、最小范数二次无偏估计方法方差分量模型中方差分量的最小范数二次无偏估计（MinimumNormQuadraticUnbiasedEstimator,MINQUE）是提出的，思路类似于他处理奇异广义线性模型的分块逆矩阵法，先提出估计应满足的性质,根据这些性质去求解（一般是一个极小值问题），若能解出，当然就有那些设定的性质。考虑一般的方差分量模型（4.2.1），我们要估计方差分量2m及其线性函数(C1,cm),(12,m)YAY的形式，即c?2yAy首先考虑0的估计应具有的形式，因为是估计方差，可考虑采用二次型)A为待估对称矩阵。再来考虑YAY应具有的性质:(1)关于参数卩的平移不变性。若参数卩有平移:0则方差分量的估计应该不变。此时原模型变为YX0XU其二次型估计变为(YX0)A(YX0)，应该有(YX0)A(YX0)YAY这等价于2YAX00XAX00即平移不变性(4.2.50)应满足的充要条件是AX0(2)估计量的无偏性。因为E(YAY)tr(A)XAXm2i2tr(AUiUi)i1今要求对一切(T2成立：m2E(YAY)c2ciii1)tr(AUiUJCi,i1,m()(3)最小范数准则。经过研究，满足平移不变性与无偏性，尚不能唯一确定待估对称矩阵A。于是可以再加一个优良性质。如果随机效应向量Ei,i=1.,m是已知的，则0=Cd2的估计应该为C1(1-)C2(22)Cm(mr-)_()P1P2pm这里diag(-ClIpi，,CmIPm)()Pipm现在用YAY去估计0=C(T2，在满足不变性条件AX=0时，YAY(XU)A(XU)UAU()我们自然希望(4.2.56)与(4.2.58)之间相差很小，这只要求矩阵厶与UAU之间相差很小。我们选用矩阵范数”UAU-AII来度量A与UAU之间的差异，即应选A使UAU|min()范数可选欧氏范数Btr(BB)()若记V=UU，则因A是分块对角阵及无偏性要求：UAUtrUAUUAU2UAU2tr(AUUAUU)2tr(AUU)trtr(AV2)m2i1tr(AV2)n2i1tr(AV2)2tr2tr(AV)-tr2Ci2tr(AUiUi)tr2PiCi+2citrPi2.2tr(4261)A是已知的，则极小化范数|UAU等价于极小化tr(AV)2。总结上述三项优良性要求，求0=c厅2的最小范数无偏估计的问题,归结为求下述极值问题:2(L1)mintr(AV)AX0sttr(AVi)Ci,i1,m因为目标函数是矩阵的迹，所以称为最小迹问题。m要解决极值问题(Li)的解，可以先对其简化。因为VVi,ViUiUii1亠1,V正疋而V存在。BV2AV2,ZV2X,WV2ViV2(4262)则模型(L1)可变为(L2)mintr(B2)BZ0sttr(BW)Ci,i1,m定理421极值问题L2的解为B*Pz(iWi)Pzi1(4263)其中1,m为方程组m(4264)(4265)tr(PzWjPzWj)iCj,j1,mi1的解，Pz为L(Z)的正交补空间L丄(Z)上的投影阵：PzIZ(ZZ)Z证明先证方程组(4.2.64)兼容。设B0满足极值问题(L2)的约束条件，即B0Z0,tr(B0Wi)G,i1,m()则L(B。)L(Z)o因是Pz往L(Z)上的投影阵，故B。BPzPZBFZ，记W；PzWjPz,j1,m，则Cjtr(BoWj)tr(PzBoPzWj)tr(BPZWjPz)tr(BoWj*)()引进拉直算符，Wj表示将Wj按列拉成一个n2X1的向量，Bo表示将Bo按列拉直成一个n2X1的向量。定义n2xm矩阵G,G(W；w2wm)()L(G)表示按G的各列展成的子空间，是n2维。定义口1为向量Bo在L(G)上的投影1PgBo()口2为向量Bo在L丄(G)上的投影2FGBo,2Wio,i1,m()具体形式PgG(GG)G,PgIPg,则12B0)这个情形与最小二乘的基本原理图标（图121）样，那里是将向量Y向子空间L（X）投影，于是存在常数卩1,，卩p,使Y=X3=Xi3i+Xp卩p，即将Y表为X列向量的线性组合。现在是必存在常数0,i1,m,可将B。表示为G的列向量的线性组合，当然也可以将Bo在L(G)的投影口1表成线性组合：m0*1i0Wi*()i1于是m0*B012iWi2()i1现在看(4267)，由公式tr(AB)AB得Cjtr(B0Wj*)B0Wj*m0*(i0Wj*2)Wj*i1mn*i0Wj*Wj*2Wj*i1mi0tr(Wj*Wj*)i1mi0tr(PZWiPZPZWjPZ)i1mi0tr(PZWiPZWj),j1,m()i1这就证明了10,0m是方程组(4.2.64)的一组解。再证明（4263）的B*是问题（L2）的解。由于PzZ0故B*满足的约束条件B*Z0。i,i1,m是从方程组tr(BWi)=C,i=1,，m中解出来的，当然满足这个约束条件。余下看B*是否是trB2的极小值解。设B为任一满足(L2)约束条件的解，记D=BB*，则D对称，DZ=0,tr(DWj)=0,PzDD,于是mtr(B*D)jtr(PzWjPzD)jimmjtr(PZWjD)jtr(WjD)0()jiji因此22tr(B2)tr(BD)(BD)tr(B)tr(D2)tr(B)()这就证明了B*是极小值解。现在回到原问题(L1)。由于(4.2.62)所定义的三个变换都是可逆变换，故(L1)与(L2)等价。于是(L1)的解存在，11*AV2BV2()对于最小范数二次无偏估计，我们介绍了它的原理、算法、解的存在性。深入的讨论还可以引出一些问题，如解的非负性，计算的复杂性等，我们这里不再讨论了。、极大似然法在假定方差分量模型随机效应服从正态的情况下，可以使用极大似然法求出参数估计。设模型为YXU11Umm()2iN(0,iIpi),i1,m()其它记号仍同以前。则YN(X,),2iUiUi(4280)Y的条件密度为f(Y|2)(2严n2exp1(Y1(YX)(4281)取对数，略去常数，得似然函数方程L(,2)21n|1(Y(4282)用矩阵微商公式|A|tXAX|A|tr(A(AA)X,A7)aX(4283)(4284)A1A1t(4285)得似然方程组1YX1XL2i1Vi)1(YX21Vi1(Y0,i1,m(4286)但是在数值计算技术发达由于似然方程组没有显式解，统计学家提出了一些迭代算法。的今天似乎没有必要再介绍这些迭代算法，甚至连求导的似然方程组都没有必要。对于样本Y,Yn),X(X1i,XPi),i1,n，以及已知设计阵VUU，代入到似然函数(4.2.82)中，调用本软件所附的计算极值的程序(由Sargent改进的Powell算法)，就可以计算出B2的估计。第三节方差分量模型参数的广义岭估计在这一节我们先将方差分量模型的方差分量化为派生模型的均值参数，分别作出其相对于LSE和BLUE的广义岭估计，再根据二步估计法作出原模型均值参数的广义二乘估计及其进一步的岭估计，证明了这样不仅使方差分量估计的均方误差减少，而且使原模型均值参数估计的均方误差也不增加和进一步减少。我们还找到了岭参数仅仅依赖于样本的估计。这样既将岭估计方法推进至方差分量模型，也改进了方差分量模型参数的离差均值对应方法。、方差分量岭估计的构造与性质岭估计和Stein估计是减少均方误差的行之有效的方法。已有文献将岭估计推广到多元线性回归模型和设计阵是列降秩的情况。本节则试图将它们推广到方差分量模型，以期改进离差均值对应方法。设有一般方差分量模型yXU11Umm（）这里X（np,）,Ui（npi）,i1,m为已知设计阵，血i）为固定效应向量，i（pii）为随机效应向量。假定E（i）0,Var（i）iIpi,i1,m;Cov（i,j）0,ij记VUiUi，i1,m,（!，,m）,1Vi则模型（431）可记为E（y）X,Var（y）（）按照两步估计法我们先求方差分量6的估计，为此采取离差均值对应。记QInXX,YQy，于是有模型mE(Y)0,Var(Y)i(QVjQ)i1再记ZYY,D(QViQQVmQ),F(QQ)Var(yy)(QQ),因为ZQyyQ=(QQ)yy，可得派生模型E(Z)D,Var(Z)F(434)此时原模型的方差分量成了派生模型的均值参数。对派生模型(4.3.4)用最小二乘法求出均值参数tr(QV1QVJL(DD)1DZtr(QVmQVJtr(QV1QVm)1yQV1Qy()tr(QVmQVm)yQVmQyer的最小二乘估计(LSE):的最佳线性无偏估计?L是e的无偏估计，但由于派生模型是广义线性模型，它未必是e(BLUE)?b。当F已知时,只是在一些特殊情形研究证明了B(DF1D)1DF1Z(4.3.6)不难举例说明?l可能出现负值分量,也不难看出DD很容易呈现病态。当X,V1,，因此需要对?l加以改进。对?b亦然。Vm中有两个接近相等时，DD的列向量就呈复共线。考虑原模型中心化与派生模型(4.3.4)中心化的关系，有引理431若模型(4.3.1)设计矩阵X,U1，Um均已中心化，则派生模型(4.3.4)也已中心化。证明已知1nX1nUl1nUm这里1n(1,1,1),则QViQ(1n1n)QUiUiQ(1nh)(QUjQU).(InQUJ(InQUi)f；(1nUi1nXXUi)(1nUi1nXXUJI.0,i1,m故1n2D0证毕我们不妨设D已经中心化、标准化。当F未知时，对模型(434)定义厅的关于?l的广义岭估计为1?l(K)(DDPKP)DZ()这里K=diag(匕,km),ki0,i=1,m。P是正交方阵致P(DD)P=diag(入1,，入m)=人。当k1=k2=km时，就成为狭义岭估计。我们将看到狭义岭估计虽然能改进?l的均方误差，但不一定能保证第二步估计时？的均方误差不增加。?l(K)的第i个分量记作？Li(K)。于是?l(K)(?L1(K),?Lm(K)()当F已知时，定义的关于？b的广义岭估计为?b(K)(DF1DPKP)1DF1Z(4.3.9)?b(K)(?B1(K),?Bm(K)(4310)这里P(DF-1D)P=A。自然这里K、P、A与(437)不同。引理若?l、?b存在，则?l(K)?L，?B(K)?B证明仅证后式，由？b(K)定义，注意K0,II?b(K)(DF1dPKP)1(DF1D)?b?b|(K)P(K)1P1P?b引理|p?b|?b|证毕说明？l(K)、?b(K)分别是对?L和？B的压缩。但这种压缩是模长的压缩，有的分量可能还有伸长。能否使各分量都压缩呢?我们有引理433设?l存在，0t0使(4.3.13)成立，同时有MSE(?LK)MSE(?L)()或使(4.3.14)成立，同时有MSE(?BK)MSE(?B)()证明因为E?LE?LK，故MSE(k)trVar(X?L1(K)X)1X?J(K)y)由引理及定理，存在K0，使(4.3.13)成立，且有?l(K)t?L,0t1。于是MSE(?lk)1trVarX1?l1XX1?L1ytttrVar(X?l1X)1X?JyMSE(?l)即(4324)成立。同理可证(4314)、(4325)可以同时成立。证毕如果我们找不到引理433的方法，得不到定理434的结论，那么我们的两步估计法对岭估计而言就将失去意义。如果发生了MSE(?lk)MSE(?l)，那么再进一步对卩采用岭估计我们就无法判断最终B的估计的均方误差到底是减少还是增大了。现在有了定理434,我们就可以对广义线性模型(432)继续采用广义岭估计：?LM(K)(X?l1(K)XPKP)1x?L1(K)y(4326)C1这里Kdiag(ki,kp),kj0,i1,po。P是正交方阵致P(X?lTk)X)Pdiag(i,po)。则又存在I，使均方误差进一步缩小：MSE(?lk(K)mse(?LK)()且存在,01,LK(K)tLK(4328)利用这种形式，可以作出不依赖于未知参数的K的估计。对？bk(K)亦然，不再赘述。第四节方差分量模型参数经验Bayes估计在这一节我们将构造误差正态的方差分量模型参数的经验Bayes估计，定义逆拉直运算,利用作者自己曾经证明的多参数指数族EB估计收敛结果，证明所构造的经验Bayes估计的收敛性。要利用的多参数指数族参数经验Bayes估计收敛速度的结果概述如下：设t=(ti,tp)以及t,t(n)为来自多元密度f(t)的当前样本和历史样本，它们都是独立同分布的，f(t)和它的偏导数f(r)(t)都是局部有界的，r+2是构造密度核估计正交多项式空间的Bayes风险为R(G)，经验Bayes维数。密度f(t|e)的参数e有先验分布G(e)，它的估计的风险为Rn(G)。我们证明了:定理1设1/2WWW1，且E2rEtft(r)(t)/f(t)2Etft(t)/f2(t)Rn(G)q2(r1)/(2rp)R(G)O(n),qn、方差分量模型参数经验Bayes估计的构造设有方差分量模型yxU11Um(441)这里X是kxpo矩阵，Ui,i=1,m是kxpi矩阵。卩是poX1固定效应向量,Ei,i=1,m是PiX1随机效应向量。假定Ei互相独立并且iN(0,2Ipi),i1,mViUiUi,ii2ViiiCov(y)这里Vi,i=1,m以及都是kxk矩阵，令U(U1U2Um)m)这里U是kp(pp1Pm)矩阵,1向量，则模型(541)成为且有E(y)X,Cov(y)这里A是对角阵:diag(二m)令向量m)则模型的基本问题是要估计固定效应卩和方差分量(T2。以下我们要将待估参数卩和c2化为多参数指数族的参数。因为已经假定yNk(X,)，所以y的条件密度为f(y|,2)(2)k/2|Ik/2exp(1/2)(yX)1(yX)(2)k/2|k/2exp(1/2)X1Xy1yy1X由于我们可以定义参数变换为:这里B是pxi向量，p=k(k+1)。这里t是pxi向量，令c()(2则t的条件密度化成了多参数指数族:11tr(yy)tr(yy)1(1/2)121X(1,2)取资料变换为yy2yt(,)1yy)k/2k/2exp(1/2)1xf(t|)c()exp(t)如果我们利用样本变换数据t,t(1),t(n)作出了参数e的估计，则我们可以根据参数变换的逆变换作出卩和的估计，进而作出卩和c2的估计。至于如何利用t,t，t(n)作出B的经验Bayes估计，可以参看第七章第三节与第十章第四节，这里从略。二、方差分量模型参数经验Bayes估计的收敛性对一般的线性模型，当c2是一个数的时候，作者曾得到卩2、c2的联立经验Bayes估计的收敛速度。现在我们对方差分量模型，c2是一个向量，无法得到EB估计的收敛速度，但可以推出收敛性。令?dn(t)为参数e的eb估计,?EdG(t)是0的Bayes估计，由定理1知E?E0,(n于是有E?E?B120(n),?E2?B220(n由于均方收敛必定导致依概率收敛，我们就有?B1?B2?E1?E2所以须标注脚标。如现在我们定义拉直运算的逆运算。由于逆拉直运算结果可能不唯一，果对拉直运算有这里B是将A的列向量顺次排列而成，则我们定义逆拉直运算为BBmnmn11 这里A的列向量是将长向量B截成n等分并顺次排列。如34(1234)13(1234)24当有多重运算复合时，我们规定从里层做起。当逆拉直成方阵时，我们省去脚标。如B为n2X1向量，则(B)1表示先将B逆拉直成nxn方阵，再求逆矩阵，再拉直成n2x1的列向量。显然，拉