第四章--方差分量線性回歸模型

1、第四章 方差分量線性回歸模型本章考慮的線性模型不僅有固定效應(yīng)、隨機(jī)誤差，而且有隨機(jī)效應(yīng)。我們先從隨機(jī)效應(yīng) 角度理解回歸概念， 導(dǎo)出方差分量模型， 然后研究模型三種主要解法。 最后本章介紹關(guān)于方差 分量模型的兩個(gè)前沿研究成果，是作者近期在應(yīng)用數(shù)學(xué)學(xué)報(bào)與國(guó)際數(shù)學(xué)雜志 Communications in Statistics上發(fā)表的。第一節(jié) 隨機(jī)效應(yīng)與方差分量模型一、隨機(jī)效應(yīng)回歸模型前面所介紹的回歸模型不僅都是線性的，而且自變量看作是固定效應(yīng)。我們從資料對(duì)YX" Xpj；出發(fā)建立回歸模型，過(guò)去一直是把 Y看作隨機(jī)的，X,，X)看作非隨機(jī)的。 但是實(shí)際上， 自變量也經(jīng)常是隨機(jī)的， 而并不是我

2、們可以事先設(shè)計(jì)好的設(shè)計(jì)矩陣。 我們把自變 量也是隨機(jī)變量的回歸模型稱為隨機(jī)效應(yīng)回歸模型。究竟一個(gè)回歸模型的自變量是隨機(jī)的還是非隨機(jī)的, 要視具體情況而定。 比如一般情況下 消費(fèi)函數(shù)可寫為C C0 b(X T)()這里X是居民收入，T是稅收，C0是生存基本消費(fèi)，b是待估系數(shù)。加上隨機(jī)擾動(dòng)項(xiàng)，就是一 元線性回歸模型C C0 b(X T)()那么自變量到底是固定效應(yīng)還是隨機(jī)效應(yīng)?那要看你采樣情況。如果你是按一定收入的家庭去調(diào)查他的消費(fèi), 那是取設(shè)計(jì)矩陣, 固定效應(yīng)。 如果你是隨機(jī)抽取一些家庭, 不管他收入如何都 登記他的收入與消費(fèi),那就是隨機(jī)效應(yīng)。對(duì)于隨機(jī)效應(yīng)的回歸模型, 我們可以從條件期望的角度推

3、導(dǎo)出與最小二乘法則等價(jià)的回歸 函數(shù)。我們希望通過(guò) X預(yù)測(cè)Y,也就是要尋找一個(gè)函數(shù) Y M(X) M (X1, ,Xp)，當(dāng)X的觀 察值為 x 時(shí),這個(gè)預(yù)測(cè)的誤差平均起來(lái)應(yīng)達(dá)到最小,即EY M (X)2 min EY L(X)2()這里min是對(duì)一切X的可測(cè)函數(shù)L(X)取極小。由于當(dāng)M (X) E(Y|X)( )時(shí)，容易證明EY M(X)M(X)L(X)0()故當(dāng) M(X) E(Y | X)時(shí)，EY L(X)2 EY M (X)2 EM (X)L(X)2()要使上式左邊極小，只有取 L(X) M (X) E(Y | X)。這個(gè)結(jié)果告訴我們，預(yù)測(cè)函數(shù)取作條件期望E(YX)時(shí)，可使預(yù)測(cè)誤差最小。我們

4、還可以證明，此時(shí) MX>=YX)與Y具有最大相關(guān)，即(Y,M(X) max (Y, L(X)()這里p表示相關(guān)系數(shù)。這是因?yàn)楫?dāng) M(X) E(Y | X)時(shí)，易證 Cov(Y,L(X) Cov(M (X), L(X),同時(shí)Cov(Y,M (X) Cov(M (X), M (X),于是2 22 (Y L(X) Cov (丫丄(X) Cov (M 丄(X)'D(Y)DL(X)D(Y)DL(X)2Cov (M (X), L(X) DM (X) DM (X)DM (X)DL(X) D(Y) DM (X)2(M(X),L(X)2(Y,M(X)2(Y,M(X)等號(hào)當(dāng)且僅當(dāng)| (M (X),L

5、(X) | 1()時(shí)成立，此時(shí)L(X)是MX>的線性函數(shù)。與表達(dá)了 M(X) E(Y| X)的極好性質(zhì)，我們稱Y M (X) E(Y| X)()為Y關(guān)于X的回歸曲線。上面的L(X)可取一切函數(shù)。如果限定 L(X)是X的線性函數(shù)，即要限定E|Y ( 01X1mXm) |(4.1.17 )mLin(4.1.10 )這里min是對(duì)X的一切線性函數(shù)取極小，則稱滿足上式的線性函數(shù)為Y關(guān)于X的回歸直線。我們可以求出0,1, m的解。記（1,m),則L( 0,) E|Y ( 01X1mXm)|2!b2Rxx2 RxyD(Y)(4.1.11 )這里b E(Y) ( 01EX1mEXm)(4.1.12 )

6、Rxx E( XEX )(XEX)DX1Cov(X1,X2)Cov(X1,Xm)(4.1.13 )COV(Xm,X1)Cov(Xm,X2)D(Xm)Rxy (Cov(Y, X 1), Cov(Y, X m)(4.114 )對(duì)L（卩o,卩）求微分（矩陣微商公式（XAX） 2AX ）得:Xb 0RxxRxy(4.1.15 )解得(4.1.16 )EY 'E(X)Rxx Rxy這里當(dāng)然假定RxX存在，否則使用廣義逆。此時(shí)的預(yù)測(cè)誤差方差是L( J) E|Y ('mXm)|22 ?RxyD(Y)Rxy Rxx RxyXY丄(RXYRXXRXY)2 / Y)為復(fù)相關(guān)系數(shù)。它指出了 Y與多元

7、變量X Xi,Xm之間的線性相關(guān)程度，是一元相關(guān)系數(shù)Cov(X,Y)DX DY的推廣。從條件期望角度我們導(dǎo)出的隨機(jī)效應(yīng)回歸模型的回歸直線表達(dá)式，與從最小二乘角度導(dǎo)出的固定效應(yīng)的回歸方程，表達(dá)式是等價(jià)的，所以從計(jì)算角度，我們不怎么區(qū)分。、方差分量模型概念上段我們建立了隨機(jī)效應(yīng)概念，將自變量也視作隨機(jī)變量，這就可以導(dǎo)出方差分量模型。 方差分量模型研究工作的奠基人是我國(guó)最早的統(tǒng)計(jì)學(xué)家許寶馭馬先生。還是剛才提到的消費(fèi)函數(shù)回歸模型，我們作隨機(jī)抽樣?？紤]居民按職業(yè)的分類，如工人、教師、醫(yī)生、律師、店員等等，記為Xi,i 1, ,m，我們從這些職業(yè)中隨機(jī)抽取了n個(gè)樣本，則模型可寫為Cij Co b(Xi T

8、i) ij, j 1, ,n,i 1, ,m()這里X可看作是第i種職業(yè)對(duì)收入的效應(yīng)。如果我們事先安排好取哪個(gè)職業(yè)的，當(dāng)然X是固定效應(yīng)?？墒俏覀儸F(xiàn)在對(duì)職業(yè)選取是隨機(jī)的，而且我們還想研究職業(yè)效應(yīng)的方差，這就導(dǎo)入了Var(Cj)方差分量模型，因?yàn)楝F(xiàn)在 C的方差由兩部分組成：為了數(shù)學(xué)符號(hào)統(tǒng)一，我們將經(jīng)濟(jì)學(xué)中的符號(hào)改過(guò)來(lái)，剛才建立的模型是YjU1 1iij, i 1,m, j 1, n(4.1.22它有一項(xiàng)固定效應(yīng) 口，一項(xiàng)隨機(jī)效應(yīng)E 影響，還可以加進(jìn)第二個(gè)隨機(jī)效應(yīng) E 2,1，一項(xiàng)隨機(jī)誤 于是可得模型£。如果還要考慮地區(qū)因素對(duì)消費(fèi)的U1 1 U2(4.1.23這次我們省掉了取值的標(biāo)記，Y的

9、方差由三項(xiàng)組成。一般地，我們建立方差分量模型如下：Y X U1 1(4.1.24這里有固定效應(yīng)向量 卩，隨機(jī)效應(yīng)向量m)(4.1.25并且將隨機(jī)誤差項(xiàng)&也并入了隨機(jī)效應(yīng)向量去。設(shè)計(jì)矩陣X以及U (U1,U2, ,U m)(4.1.26 )都是已知的。對(duì)于隨機(jī)效應(yīng)1,m ,合理的假定是E(i)0,Cov( i, j)0,i j2(4.1.27 )D(i)i ,i 1, ,m當(dāng)然以后有時(shí)還可以考慮E i是向量的情況，不過(guò)這里假疋每個(gè)1是一維變量。記ViUiUi,i1,2、/ 2X/,m,1 V1mVm ,(4.1.28 )則方差分量模型可記為E(Y) X , Var(Y)()模型的主要任務(wù)

10、是要估計(jì)固定效應(yīng)向量卩與方差分量；，I, , m。和一般的多元線性回歸模型相比，就是待估的方差多了。通過(guò)這些介紹，我們就可以方便地將各種經(jīng)濟(jì)方面的普通線性回歸模型改造成方差分量模型，當(dāng)然要根據(jù)實(shí)際。(4.2.1)第二節(jié)方差分量模型的解法對(duì)于方差分量模型Y XU11U m mn 1n p p 1 n P1 P1 1n Pm Pm 1mE(Y) X , Var(Y)Q ii 1般都采用二步估計(jì)法,首先估計(jì)方差分量然后再估計(jì)固定效應(yīng)卩。按照廣義最小二乘(422 )* (X ? 1X) X ? 1Y其中2?iUiUii 1所以方差分量模型解法的關(guān)鍵是估計(jì)方差分量。以下介紹的方法，也都是針對(duì)方差分量估計(jì)

11、方法而言的。、方差分析法先從一個(gè)簡(jiǎn)單的模型結(jié)合數(shù)據(jù)結(jié)構(gòu)形象地說(shuō)明方法。考慮模型Yij0 i ij ,i 1, m, j 1, n卩o為總平均，是固定效應(yīng)，E 1,，E m是隨機(jī)效應(yīng)，E i 0, Cov( i, j) 0,i j ,Var( i) A,i 1, ,m。對(duì)于隨機(jī)誤差0, Var( )2。這個(gè)模型如果記作方差分量模型的標(biāo)準(zhǔn)形式是YX 0 U(4.2.5 )其中設(shè)計(jì)陣X=(1,1, 1)'，隨機(jī)效應(yīng)矩陣為10010k01U011(4.2.6 )01mkm我們手中資料只有Y(丫11 , Y1k，丫21 ,，丫2k,，丫mk)我們采用(4.2.4)記法方便一些，將資料 丫排成表1

12、2k組內(nèi)平均1Y11Y12YkY12Y21Y22Y2k丫2mYm1Ym2YmkYm方差分析主要掌握三點(diǎn)， 自由度。先計(jì)算總平均：精品文檔一是計(jì)算組內(nèi)差、 組間差，二是作平方和分解， 三是計(jì)算各自的m kmk i1 j 1總變差（全體資料與總平均的偏差平方和）:mkSTi 1(Yijj 1Y)2各組平均（各組資料橫向相加并平均）-1 kY kjY|j '1i1, ,m組間差（各組平均數(shù)與總平均數(shù)的偏差平方和)kmSaj 1(Yii 1Y)組內(nèi)差（各組數(shù)據(jù)與本組平均數(shù)的偏差平方和)mkSi 1(Yjj 1Yi)2則必有平方和分解STSASYYj(427 )(428 )(429 )(4210

13、 )(4211 )(4212 )將各平方和除以各自的自由度。St有一個(gè)約束丫（427），自由度為n 1 mk 1 ； Sa有m組差，1個(gè)約束，自由度為 m 1； Se有mk組差，m個(gè)約束，自由度為 mkm注意有自由 度分解：fTfAfe,mk 1 (m 1) (mk m)于是算出均方：Qt1STmk1Qa1SAm 1Q1Smkm因?yàn)榧俣殡S機(jī)效應(yīng)，可以算出各均方的均值:(4213 )(4214 )(4215 )(4216 )精品文檔E(QA ) k 2 2(4.2.17 )(4.2.18 )E(Q ) 2以 QA 代者 E(QA) ，Q 代替 E(Q ) ，得方程組：k A22 Q A(4.2

14、.19 )2Q解得?2 Qe,?2A (QA Q ) / k(4.2.20 )這樣就作好了方差分量的估計(jì)，然后可以按(422)作出卩的估計(jì)。因?yàn)檫@里的方差分量是由方差分析法作出的，故稱為方差分析法。推廣到一般的方差分量模型時(shí)，基本原則是類似的。我們不妨考慮方差分量模型Y X U 1 1 U 2 2Cov(Y) 12U1U1 22U 2U22I先對(duì)總平方和Y' Y作平方和分解YY S S1 S 2 S其中Sb是在模型Y=X3 + &中，卩的回歸平方和：S SES( ) YX(X X) XY4.2.21 )4.2.22 )4.2.23 )S 1 是在模型 Y XU 1 1中，消去B

15、影響后E i的平方和類似地， S 2 是在模型 Y XS 1SES( , 1)SES ( )4.2.24 )中消去B和E i影響后，E 2的平方和:S 2EES( , 1, 2 )SES ( , 1)最后的S，為殘差平方和S Y YSES (, 1, 2 )可以驗(yàn)證S Y (ID)YU1 1 U 2 24.2.25 )4.2.26 )4.2.27 )4.2.28 )S1 Y (D D1)YS 2 Y (D1 D12 )Y(4.2.29 )S Y D12Y(4.2.30 )這里D I X(XX) X I PX(4.2.31 )D1 D DU1(U1DU 1) U1D D PDU1(4.2.32

16、)D12D 1D1U 2 (U 2 D1U 2 ) U 2 D1D1PD1U 2(4.2.33 )這里P表示關(guān)于*的投影陣。下面計(jì)算各平方和的均值。E(S 1) X (D D1)X tr(D D1 )U1U1 12 U2U2 22 2I X (D D1) X tr(U1DU 1) 1222tr(U1D1U1) 12 tr(U 2DU 2) 2222tr(U 2 D1U 2 ) 22 tr(D D1) 2 ()因?yàn)?DX 0, D1X 0，所以上式第一項(xiàng)為 0。在第三項(xiàng)中，tr(U1D1U1) trU1DU 1 U1DU 1(U1DU 1) U1DU 1 0 () 在第六項(xiàng)中tr(D D1) t

17、r D1U1(U1D1U1) U1D1tr(U1D1U1 ) U1D1U1rk(U1D1U1 ) rk(U1 D1)rk(U1X ) rk(X )rk(U1 X) rk (X) ()所以最后有E(S 1 ) c11 (c2 c3 ) 2 r22(4.3.37 )其中c1tr(U1DU1)(4.2.38 )c2tr (U 2 DU 2)(4.2.39 )c3 tr (U 2 D1U 2)r1 rk(X), r1 r2 rk(U1| X)4.2.40 )4.2.41 )于是我們得到方程組解此方程組，就可以得到固定效應(yīng)的估計(jì)。E(S2)E(S ) (n r1r3 rk(X2c2 2r3)U1U 2)

18、 r1 r2S1c121(c222c3 ) 2 r2S2c222r32S(nr1r2r3 ) 224.2.42 )4.2.43 )4.2.44 )4.2.45 )21,22,的估計(jì)。然后進(jìn)入二步估計(jì)的第二步，就可以得到關(guān)于類似還可以求得算例 市場(chǎng)收益率與股利和換手率的關(guān)系 考慮一個(gè)隨機(jī)效應(yīng)的多元線性模型YX n 1 n p p 1nmm1U 的形式如同 。問(wèn)題的實(shí)際背景是，觀測(cè)對(duì)象被分成了 m組，可能存在一個(gè)隨機(jī)效應(yīng)向量對(duì)各組資料有不同的作用。模型也可以寫作YijXij i ij , i 1, ,m, j 1, ,k數(shù)據(jù)結(jié)構(gòu)及具體數(shù)值如下表所示，m=6,k=6。這些資料采自'96上海股

19、票市場(chǎng)資料總匯。我們研究目的一是看過(guò)去一年的股利收入與當(dāng)年換手率對(duì)當(dāng)年市場(chǎng)收益率有何影響，二是想知道是否存在一個(gè)潛在的尚未觀測(cè)到的隨機(jī)效應(yīng)， 對(duì)行業(yè)有明顯影響。 當(dāng)然這種情況采用方 差分量模型比較合適。要注意本例是兩個(gè)方差量， 上一章第二節(jié)模型 也是兩個(gè)待估的方差量。 它們的隨 機(jī)效應(yīng)作用范圍不一樣，不是一回事。表 年股市資料類別股號(hào)股名1996年收益率%1995年股利%1996年日換手率商 業(yè) 類628新世界64.769203.12631中百一店46.84511.81.68632華聯(lián)商廈41.95811.31.81655豫園商城16.19511.21.10682新百公司79.9115.23.

20、36694大連商場(chǎng)91.3885.84.26電 子 類602真空電子33.112103.52651飛樂(lè)音響8.10801.95800天津磁卡271.76353.74839四川長(zhǎng)虹381.686604.41850華東計(jì)算 機(jī)14.63813.27870廈華電子68.5793.94.20化工類617聯(lián)合化纖-21.8710.51.53618鹵堿化工22.37022.63672廣東化纖11.86004.65688上海石化179.8173.54.38886湖北興化236.32852.34.45889南京化纖-33.1222.44.644.11醫(yī) 藥 類664哈醫(yī)藥191.66654.32671天目藥業(yè)

21、111.135164.11812華北制藥152.0158.44.11842中西藥業(yè)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機(jī) 械 類604二紡機(jī)14.4102.31605輕工機(jī)械6.12203.50610中紡機(jī)0.70102.27806昆明機(jī)床41.8521.14.22841上柴股份66.9812

22、02.20862南通股份41.093201.36首先我們作普通最小二乘回歸，得到(XX) 1XY , 然后計(jì)算Yj*YijXjj ?。此時(shí)的Yj已消除固定效應(yīng)影響，我們將它排成平面表，以作方差分析，計(jì)算A與2。計(jì)算過(guò)程從 至。從下面計(jì)算過(guò)程可以看到， 平方和分解式是滿足的：SST SSA SSE即（St=S+SJ。對(duì)于本例資料，隨機(jī)誤差 4055.6 遠(yuǎn)大于隨機(jī)效應(yīng)方差，組內(nèi)差遠(yuǎn)大于組間差，可以認(rèn)為隨機(jī)效 應(yīng)不明顯，即行業(yè)差別不明顯。對(duì)于選定的方差分量模型，回歸結(jié)果是Y 0.1373 4.7053X1 3.3258 X 2它的標(biāo)準(zhǔn)差很小， 為 1.0084 ，這正是采用方差分量模型廣義最小二乘

23、意義所在。 擬合效果圖 （圖 令人滿意。方差分量模型方差分析法計(jì)算程序，例421.第一列為 Y, 以后各列為 X例 421.D 數(shù)據(jù)文件中 ， m=6， k=6， p=2 要顯示原始資料嗎 ? 0= 不顯示 ， 1= 顯示（ 0）先作普通最小二乘 ， 并打印結(jié)果 :現(xiàn)在作線性回歸顯著性檢驗(yàn) ， 計(jì)算 t，F(xiàn)，R 統(tǒng)計(jì)量請(qǐng)輸入顯著性水平 a， 通常取 a=0.01， 0.05， 0.10， a=?（0.05） 線性回歸分析計(jì)算結(jié)果樣本總數(shù) 36自變量個(gè)數(shù) 2回歸方程 Y = b0+b1*X1+.+b2*X2Y =5.3188 +4.5656 X1 +6.0128 X2回歸系數(shù) b0， b1， b

24、2， .， b2 5.31884.56566.0128殘差平方和 : 149400.60 回歸平方和 : 131759.30誤差方差的估計(jì) : 4150.0170 標(biāo)準(zhǔn)差 = 64.4206線 性 回 歸 顯 著 性 檢 驗(yàn) 顯著性水平 : .050回歸方程整體顯著性 F檢驗(yàn),H0:b0=b仁=b2=0F 統(tǒng)計(jì)量 :14.5517 F 臨界值 F（2， 33） 3.285全相關(guān)系數(shù) R :.6846回歸系數(shù)逐一顯著性 t 檢驗(yàn) , H0:bi=0, i=1,.,2 t 臨界值 t( 33) 1.6924回歸系數(shù) b1-b 2 的 t 值 : 5.6318 1.1073打印方差分析資料 Y(I,

25、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計(jì)算各種平均: 總平均YYba:.0000各組平均

26、 Yba: -16.216 -2.34139.04534.244-18.463 -36.270計(jì)算各種變差: 總變差SST, 組間差 SSA,組內(nèi)差 SSESST= 149400.6000 SSA= 27731.9200 SSE= 121668.7000 打印方差分量的估計(jì)SIGMAA= 248.4600 SIGMAE= 4055.6240 下面計(jì)算協(xié)方差陣及其逆的一塊 , 并作分解計(jì)算廣義最小二乘估計(jì),模型轉(zhuǎn)換Y=PY, X=PX,下面打印矩陣P的一塊,P ' P近的逆SIG1=.1239-.0007-.0007-.0007-.0007-.0007SIG1=.0000.1238-.0

27、008-.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下面打印的是廣義最小二乘的統(tǒng)計(jì)結(jié)果 : 現(xiàn)在作線性回歸顯著性檢驗(yàn) , 計(jì)算 t,F,R 統(tǒng)計(jì)量 請(qǐng)輸入顯著性水平 a, 通常取 a=0.01, 0.05, 0.10, a=? (0.05)線性回歸分析計(jì)算結(jié)果樣本總數(shù) 36 自變量個(gè)數(shù) 2回歸方程 Y = b0+b1

28、*X1+.+b2*X2Y =.7309 +4.5682 X1 +5.7726 X2回歸系數(shù) b0, b1, b2, ., b2 .7309 4.5682 5.7726殘差平方和 : 2290.69 回歸平方和 : 2014.61誤差方差的估計(jì) :63.6303 標(biāo)準(zhǔn)差 = 7.9769線 性 回 歸 顯 著 性 檢 驗(yàn) 顯著性水平 : .050回歸方程整體顯著性F檢驗(yàn),H0:b0=b 仁=b2=0F 統(tǒng)計(jì)量：14.5113F臨界值 F(2, 33) 3.285全相關(guān)系數(shù)R :.6841回歸系數(shù)逐一顯著性t檢驗(yàn),H0:bi=0, i=1,.,2t臨界值 t( 33) 1.6924回歸系數(shù)b1-b

29、 2的t 值:5.6411 1.0683比較殘差平方和：普通最小二乘的：149400.6000:廣義最小二乘的：2290.6910下面打印的是利用廣義最小二乘的回歸系數(shù)去計(jì)算原始資料的回歸擬合結(jié)果 要作回歸預(yù)測(cè)嗎？鍵入0=不預(yù)測(cè)，仁要預(yù)測(cè)（0）回歸系數(shù)：.73094.56825.7726要打印擬合數(shù)據(jù)嗎？ 0=不打印，仁 打印（0）計(jì)算結(jié)束。圖 421.1原始數(shù)據(jù)擬合數(shù)據(jù)計(jì)算中需要一些分析與技巧。因?yàn)锳uu2i需要用它及計(jì)算卩(X?1X) 1X?1Y注意UU是一個(gè)分塊對(duì)角陣，共有 m個(gè)對(duì)角塊，每塊是一 k x k方陣，元素皆為1。這樣工的計(jì)算就容易了。恰好精品文檔X也是分成m塊，每塊為kx p

30、陣。于是XiX ?1X (Xi, ,Xm)diag(?i,?m)XmXi ?iXim Aii 2I m，這就大大簡(jiǎn)化了計(jì)算。如果真的要計(jì)算36X 36的矩陣工的逆，一般微機(jī)是不可能的。這個(gè)程序開始時(shí)調(diào)用了普通多元線性回歸程序?qū)υ假Y料回歸。得到的回歸方程為Y 5.31884.5656Xi 6.0128X264.4206，遠(yuǎn)大于方擬合效果也很好(見圖421.2)。但是，這個(gè)模型的殘差向量的標(biāo)準(zhǔn)差為 差分量模型的1.0084。這可以對(duì)比說(shuō)明方差分量模型的作用。圖 421.2t一原始數(shù)據(jù) -擬合數(shù)據(jù)、最小范數(shù)二次無(wú)偏估計(jì)方法方差分量模型中方差分量的最小范數(shù)二次無(wú)偏估計(jì)(Mi nimum Norm Q

31、uadratic Un biasedEstimator,MINQUE是 提出的，思路類似于他處理奇異廣義線性模型的分塊逆矩陣法，先提出估計(jì)應(yīng)滿足的性質(zhì)，根據(jù)這些性質(zhì)去求解(一般是一個(gè)極小值問(wèn)題)，若能解出，當(dāng)然就有那些設(shè)定的性質(zhì)?？紤]一般的方差分量模型，我們要估計(jì)方差分量；，,m及其線性函數(shù)c , C (c1 , Cm ),( 1 , m )( 4246 )首先考慮 $的估計(jì)應(yīng)具有的形式，因?yàn)槭枪烙?jì)方差，可考慮采用二次型Y' AY的形式，即? c ?2 Y A?Y（4.2.47 ）A為待估對(duì)稱矩陣。再來(lái)考慮 Y'AY應(yīng)具有的性質(zhì)：(1)關(guān)于參數(shù)卩的平移不變性。若參數(shù)卩有平移：0

32、（4.2.48 ）則方差分量的估計(jì)應(yīng)該不變。此時(shí)原模型變?yōu)閅 X 0 X U（4.2.49 ）其二次型估計(jì)變?yōu)?(Y X 0) A(Y X 0) ，應(yīng)該有（Y X 0 ） A（YX 0）Y AY（4.2.50 ）這等價(jià)于2Y AX 00X AX00（4.2.51 ）即平移不變性 (4.2.50) 應(yīng)滿足的充要條件是AX0(2) 估計(jì)量的無(wú)偏性。 因?yàn)镋（Y AY） tr（A ）X AXmi1i2 tr （ AUiUi）（4.2.53 ）今要求對(duì)一切 刖成立：E（Y AY）2 cm2 ci i i1（4.2.54 ）則其充要條件應(yīng)為tr（AUiUi）ci,i 1, , m（4.2.55 ）(3)

33、 最小范數(shù)準(zhǔn)則。經(jīng)過(guò)研究，滿足平移不變性與無(wú)偏性，尚不能唯一確定待估對(duì)稱矩陣A。于是可以再加一個(gè)優(yōu)良性質(zhì)。如果隨機(jī)效應(yīng)向量 j,i=1,，m是已知的，則$ =c'異的估計(jì)應(yīng)該為G(丄)C2(2)mm、Cm()P1P2pm這里C1diag(ICm 1 )P1Pm /P1P m現(xiàn)在用Y'AY去估計(jì)0 =c，CT2,在滿足不變性條件AX=0 時(shí)，(4256 )(4257 )YAY (X U )A(X U ) U AU(4258 )我們自然希望與之間相差很小，這只要求矩陣與UAU之間相差很小。我們選用矩陣范數(shù)” U'AU A H來(lái)度量A與U'AU之間的差異，即應(yīng)選 A使

34、U AU | min范數(shù)可選歐氏范數(shù)(4259 )tr(BB)(4260 )若記V=UU，則因A是分塊對(duì)角陣及無(wú)偏性要求:U AUtrU AUU AU 2U AU22tr(AUU AUU ) 2tr(AU U ) tr2tr(AV )m ci 1 PLtr(AUiUi) tritr(AV2)n Ci 1 P2citri2tr(AV )2tr 2tr 22tr(AV)-tr 22(4261 )A是已知的，則極小化范數(shù)|UAU等價(jià)于極小化tr(AV)2?？偨Y(jié)上述三項(xiàng)優(yōu)良性要求，求 0=c'c2的最小范數(shù)無(wú)偏估計(jì)的問(wèn)題，歸結(jié)為求下述極值問(wèn)題：(L1)mi n tr(AV)2AX 0s ttr

35、(AVi) Ci, i 1, ,m因?yàn)槟繕?biāo)函數(shù)是矩陣的跡，所以稱為最小跡問(wèn)題。則模型(L1)可變?yōu)? 1B V2AV2,Z1v 2x,w v(4262)(L2)min2tr(B )BZ 0tr(BW) C,i 1, ,m定理421極值問(wèn)題L2的解為mPz (iWi)FZi 1(4263)其中m為方程組mtr(PzWiPzWj) i Cj,i 11,m(4264 )的解，PZ為L(zhǎng)(Z)的正交補(bǔ)空間L丄(Z)上的投影陣:PzI Z(ZZ) Z(4265 )先證方程組 設(shè)B0滿足極值問(wèn)題證明(4.2.64)(L2)的約束條件，即兼容。BcZ0,tr(B°Wi) Ci, i1,m(4266

36、)則 L(B。) L (Z)。因是Pz往L (Z)上的投影陣B0B0 PzPz B°Pz，記WjPzWjPz , j 1,m，則Cjtr(B°Wj)tr(Pz B°PzWj)tr(B°FZWjPz )、 . 1Vi ,ViU iU i , V正定而V 存在。要解決極值問(wèn)題(L1)的解，可以先對(duì)其簡(jiǎn)化。因?yàn)閂tr(B°W*)(4267 )22引進(jìn)拉直算符，W j表示將 W按列拉成一個(gè)n x 1的向量，B0表示將Bo按列拉直成一個(gè) n x 1G(W1* W2*Wm* )(4.2.68 )L( G) 表示按G的各列展成的子空間，是n2維。定義口 1為

37、向量Bo 在 L( G) 上的投影1PGBo(4.2.69 )口 2為向量B0在L丄(G上的投影2PGBo,*2Wi*o, i1, ,m(4.2.7o )具體形式 PGG(GG) G ,PG IPG , 則12Bo(4.2.71 )這個(gè)情形與最小二乘的基本原理圖標(biāo)(圖一樣，那里是將向量 Y向子空間L(X)投影，于是存在常數(shù)卩1,，卩p,使Y=X3=Xi 3 1+XpBp,即將Y表為X列向量的線性組合?，F(xiàn)在是必存在常數(shù) 0,i 1, ,m，可將Bo表示為G的列向量的線性組合，當(dāng)然也可以將Bo在L(G)的投影口 1表成線性組合：mo*1io Wi*()i1于是mo*Bo12ioWi*2()i1現(xiàn)在

38、看(4267)，由公式tr(AB) A'B得Cjtr(BoWj* ) Bo'Wj*mo *(ioWj*2) Wj*i1*2Wj*mo *ioWj* 'Wj*i1 mmi1o *iotr(Wj*Wj*)iotr(PZ Wi PZ PZWjPZ )i1(4.2.74 )Otr(PzWjPzWj), j 1, ,mi 1這就證明了 1, m是方程組(的一組解。再證明的B*是問(wèn)題(L2)的解。由于PZZ 0故B滿足的約束條件B*Z0。i ,i 1, ,m是從方程組tr( BW)=C,i=1,，m中解出來(lái)的，當(dāng)然滿足這個(gè)約束條件。余下看B是否是tr E2的極小值解。設(shè) B為任一滿

39、足(L2)約束條件的解，記 D=B- B*，貝U D對(duì)稱，DZ=0,tr( DW)=0, PZ D D ,于是mtr(B*D)jtr(PzWjPz D)j 1mmjtr(PzWjD)j 1j 1jtr(WjD) 0(4.2.75 )因此2 *tr(B ) tr(B* * 2D)(B D) tr(B )2*2tr(D ) tr(B )(4.2.76 )這就證明了B是極小值解。現(xiàn)在回到原問(wèn)題(L1)。由于所定義的三個(gè)變換都是可逆變換，故(L1)與(L2)等價(jià)。于是(L1)的解存在，1 1* *A V 2B V 2()對(duì)于最小范數(shù)二次無(wú)偏估計(jì)，我們介紹了它的原理、算法、解的存在性。深入的討論還可 以

40、引出一些問(wèn)題，如解的非負(fù)性，計(jì)算的復(fù)雜性等，我們這里不再討論了。、極大似然法在假定方差分量模型隨機(jī)效應(yīng)服從正態(tài)的情況下，可以使用極大似然法求出參數(shù)估計(jì)。設(shè)模型為Y X U1 1Um m()i N(0, i Ipi),i1, ,m()其它記號(hào)仍同以前。貝UmY N(X , ),2UiUi()Y的條件密度為f(Y|2) (2 嚴(yán) |nI 2 exp1(Y1(Y X)(4281取對(duì)數(shù)，略去常數(shù),得似然函數(shù)方程2)21 n|L(,冷(Y1(Y(4282用矩陣微商公式| A| t XAXIX| A|tr(A(AA)X,A)aXX(4283(4284得似然方程組1YX 1xL2i1Vi)2(Y X2A1t

41、(42851 1Vi(Y0,i1,m(但是在數(shù)值計(jì)算技術(shù)發(fā)達(dá)由于似然方程組沒(méi)有顯式解，統(tǒng)計(jì)學(xué)家提出了一些迭代算法。的今天似乎沒(méi)有必要再介紹這些迭代算法，甚至連求導(dǎo)的似然方程組都沒(méi)有必要。對(duì)于樣本Y(Y1,Yn),X(X1i, , X pi), i 1, ,n，以及已知設(shè)計(jì)陣 V UU，代入到似然函數(shù)中，調(diào)用本軟件所附的計(jì)算極值的程序(由Sargent改進(jìn)的Powell算法)，就可以計(jì)算出2的估計(jì)。第三節(jié) 方差分量模型參數(shù)的廣義嶺估計(jì)在這一節(jié)我們先將方差分量模型的方差分量化為派生模型的均值參數(shù)，分別作出其相對(duì)于LSE和BLUE的廣義嶺估計(jì)，再根據(jù)二步估計(jì)法作出原模型均值參數(shù)的廣義二乘估計(jì)及其進(jìn)一

42、 步的嶺估計(jì)，證明了這樣不僅使方差分量估計(jì)的均方誤差減少，而且使原模型均值參數(shù)估計(jì)的均方誤差也不增加和進(jìn)一步減少。我們還找到了嶺參數(shù)僅僅依賴于樣本的估計(jì)。這樣既將嶺估計(jì)方法推進(jìn)至方差分量模型，也改進(jìn)了方差分量模型參數(shù)的離差均值對(duì)應(yīng)方法。、方差分量嶺估計(jì)的構(gòu)造與性質(zhì)嶺估計(jì)和Stein估計(jì)是減少均方誤差的行之有效的方法。已有文獻(xiàn)將嶺估計(jì)推廣到多元線性回歸模型和設(shè)計(jì)陣是列降秩的情況。本節(jié)則試圖將它們推廣到方差分量模型，以期改進(jìn)離差-均值對(duì)應(yīng)方法。設(shè)有一般方差分量模型(431 )(p0 i)為固定效應(yīng)向量，i(pi i)為隨機(jī)效應(yīng)Jpi,i 1, ,m；Cov( i, j)0,i j 記imi iVi

43、則模型(4.3.1)可記為Var(y)(4.3.2 )(T的估計(jì)，為此采取離差均值對(duì)應(yīng)。記y X U 1 1U m m這里 X (n P0) ,U i(n pi) ,i1, ,m為已知設(shè)計(jì)陣,向 量 。假疋E( i) 0,Var( JVi UiUi,i1,m,(1 , m),E(y) X ,按照兩步估計(jì)法我們先求方差分量Q In XX ,Y Qy，于是有模型mE(Y)0,Var(Y)i(QViQ)i 1(4.3.3 )再記 Z YY ,D (QV1QQVmQ), F (QQ)Var(yy )(QQ),因?yàn)?Z Qyy Q =(Q Q)yy，可得派生模型E(Z) D ,Var(Z) F(4.3.

44、4 )此時(shí)原模型的方差分量成了派生模型的均值參數(shù)。對(duì)派生模型(4.3.4)用最小二一乘法求出均值參數(shù)T的最小一乘估計(jì)(LSE):tr(QV1QV1)1tr(QV1QVm)y QV1 Qy?L(DD)1DZtr(QVmQy)tr(QVmQVm)yQVmQy(4.3.5 )?L是T的無(wú)偏估計(jì)，但由于派生模型是廣義線性模型，它未必是T的最佳線性無(wú)偏估計(jì)(BLUE) ?B。當(dāng)F已知時(shí),1 1 1B (D F D) D F Z(4.3.6)只是在一些特殊情形研究證明了不難舉例說(shuō)明?l可能出現(xiàn)負(fù)值分量，也不難看出D'D很容易呈現(xiàn)病態(tài)。當(dāng) X, V1,，V中有兩個(gè)接近相等時(shí)，D'D的列向量就

45、呈復(fù)共線。因此需要對(duì)?L加以改進(jìn)。對(duì)?B亦然?？紤]原模型中心化與派生模型 引理已中心化。證明中心化的關(guān)系，有已知1nX 1nUinUm這里1n(1,1,1),則1n2 QViQ (1n1n)QUiUiQ (1n 1n)(QUiQU)n若模型設(shè)計(jì)矩陣X, U,，匕均已中心化，則派生模型也(1nQUi)(1nQUi)In(1nUi 1nXX Ui)(1nUi 1 n XXUi)In 0,i1,m證畢1?l(K) (D D PKP ) 1DZ(437 )這里 K=diag( ki,km) , ki>0, i=1,m P是正交方陣致 P'(D'D) F=diag(入 i,入 m)

46、= A。當(dāng)?L的均方誤差，但不kl=k2=km時(shí)，就成為狹義嶺估計(jì)。我們將看到狹義嶺估計(jì)雖然能改進(jìn)定能保證第二步估計(jì)時(shí)？的均方誤差不增加。?L (K)的第i個(gè)分量記作?Li(K)。于是?l(K)(?L1(K), ?Lm(K)(438 )當(dāng)F已知時(shí)，定義 廳的關(guān)于?B的廣義嶺估計(jì)為?b(K)(D F 1Dpkp)1d f 1z(439)?b(K)(?b1(K),?Bm(K)(4.3.10)這里P'D'F1D)P= A。自然這里 K、P、A與(4.3.7)不同。引理?B存在，則L、若?l(K)?l,?b(K)?b證明僅證后式,由？b(K)定義，注意K>0,II ?b(K)(d

47、f 1d pkp) 1(d F 1D)?bP( k) 1 P ?b(K) 1 P ?b引理P ?b?b證畢4.3.2說(shuō)明？l(K)、?b(K)分別是對(duì)?l和?b的壓縮。但這種壓縮是模長(zhǎng)的壓縮,有的分量可能還有伸長(zhǎng)。能否使各分量都?jí)嚎s呢？我們有引理4.3.3 設(shè)?l存在，0<t<1，則存在函數(shù)關(guān)系kiki(吧0,i1, ,m，使?l(K)t?L(4311)類似地有?b (K) t ?b(4312)證明 令ki1,i 1,m。則kidkiki%i 1, m。?l(K)P( k) 1P ?lP diag11 k1Pdiag(t, ,t)P ?lt?L證畢對(duì)？b（K）亦然。此時(shí)廣義嶺估計(jì)實(shí)

48、際成了均勻壓縮估計(jì)-Stein估計(jì)。我們記 ？L（t）t?L,?B（t） t ?B定理431右？L存在，則存在區(qū)域R (k1, ,km)|0 k1k；,0 km km當(dāng)（k1 ,，km) R有MSE( ?l(K)MSE( ?l)(4.3.13 )這里均方誤差MSE( ?)2E ?。類似有MSE( ?b(K)MSE(?b)(4.3.14 )證明注意雖然?L（0）不再是定義中的嶺估計(jì)，但仍是有意義的，？L（0）_?L令f （K）E ?l(K)"(K)tr Var(?L(K),f2(K)2E?L（k）。它們都是 k1,km的函數(shù)。記 B=DD+PKP，貝Uf (K) f1(K) f2(K) tr(D FDB 2) (D DB 2D D