張 璇

分層線性模型的最大后驗(yàn)估計(jì)

張 璇

(中國(guó)人民大學(xué)統(tǒng)計(jì)學(xué)院,北京100084)

最大后驗(yàn)估計(jì)(MAPE)和最大似然估計(jì)(ML E)都是重要的參數(shù)點(diǎn)估計(jì)方法。在介紹一般分層線性模型(HLM)MAPE方法的基礎(chǔ)上,給出這種方法的期望最大化算法(EM)的具體步驟,運(yùn)用對(duì)數(shù)似然函數(shù)的二階導(dǎo)數(shù)推導(dǎo)了MAPE估計(jì)的方差估計(jì)量。同時(shí)運(yùn)用數(shù)據(jù)模擬比較了EM算法下的MAPE和ML E。對(duì)于固定效應(yīng)的估計(jì),兩種方法得到的估計(jì)量是一致的。當(dāng)組數(shù)較少時(shí),EM計(jì)算的MAPE的方差協(xié)方差成分比ML E的更靠近真實(shí)值,而且MAPE的迭代次數(shù)明顯小于ML E。

分層線性模型:最大后驗(yàn)估計(jì):最大似然法:期望最大化算法

一、研究背景

自從Lindley和Smith提出分層線性模型[1]HLM(Hierarchical Linear Model)的概念以來, HLM在生物學(xué)、醫(yī)學(xué)、社會(huì)學(xué)、教育學(xué)和經(jīng)濟(jì)管理方面的應(yīng)用都得到了長(zhǎng)足的發(fā)展。目前比較常用的分層線性模型的估計(jì)方法有最大似然法(ML)和完全貝葉斯法(Full Bayes)。理論研究表明大樣本情況下,ML得到的參數(shù)估計(jì)是一致最優(yōu)估計(jì)量, Bayes估計(jì)量是有偏的,但是Bayes估計(jì)考慮了所有輔助參數(shù)的不確定性,因此在樣本較小的情況下能得到較ML更可靠的方差協(xié)方差的估計(jì)值。在計(jì)算方法上,ML運(yùn)用期望最大化EM(Expectation Maximum)算法,Bayes運(yùn)用MCMC算法,當(dāng)分層線性模型形式較復(fù)雜時(shí)(如待定參數(shù)向量增多、層數(shù)增加等),MCMC運(yùn)算并非能得到一個(gè)收斂的Markov鏈,而EM的算法簡(jiǎn)單快捷。

經(jīng)典的頻率學(xué)派和Bayes學(xué)派在統(tǒng)計(jì)推斷方法上的爭(zhēng)論由來已久,在參數(shù)估計(jì)時(shí),如果對(duì)待估參數(shù)的先驗(yàn)信息一無所知,可以用頻率學(xué)派的思想做參數(shù)估計(jì);但是如果事先知道些許先驗(yàn)信息,Bayes估計(jì)可以提高參數(shù)估計(jì)的精度,減小參數(shù)估計(jì)的方差。最大后驗(yàn)估計(jì)方法運(yùn)用了Bayes的分析理論,嚴(yán)格意義上說它的估計(jì)量并非是Bayes估計(jì)量,但是和極大似然相比,它考慮了待估參數(shù)的先驗(yàn)信息,是一種重要的估計(jì)方法。本文將介紹分層線性模型的最大后驗(yàn)估計(jì)方法以及EM算法的具體步驟。

二、分層線性模型的參數(shù)的最大后驗(yàn)估計(jì)

為了清楚地描述整個(gè)估計(jì)過程,運(yùn)用一個(gè)簡(jiǎn)單的兩層的分層模型。

上述模型的矩陣形式是:

其中Yj是nj×1的因變量向量,Xj是nj×(Q+1)的自變量矩陣,(Q為層1模型自變量的個(gè)數(shù),這里Q =1),βj是(Q+1)×1的待定參數(shù)向量,Ij是nj×nj的單位矩陣,rj是nj×1的隨機(jī)誤差向量;Wj是(Q+ 1)×F的自變量矩陣(這里F=4),γ是F×1的固定效應(yīng)向量,uj是(Q+1)×1的層2誤差向量或隨機(jī)效應(yīng)向量,T是(Q+1)×(Q+1)的方差協(xié)方差矩陣,并且通常假定模型中rj和uj是無關(guān)的,每層的自變量和該層的誤差向量以及其它層的誤差向量是無關(guān)的。

假設(shè)模型的參數(shù)向量θ=(γ,σ2,vech(T))①vech(T)表示將矩陣T的上三角元素按列拉直向量。,θ的先驗(yàn)聯(lián)合密度函數(shù)為p(θ),則θ的后驗(yàn)密度函數(shù)為f(θ|Y)=p(θ)f(Y|θ)/f(Y),f(Y)是Y的邊際密度函數(shù)和θ無關(guān),所以f(θ|Y)∝p(θ)f(Y|θ)。由式(1)知rj和uj都服從正態(tài)分布,所以θ的對(duì)數(shù)對(duì)于p(θ),在沒有太多信息的前提下,采用Jeffry先驗(yàn)密度函數(shù)[2],令p(θ)=C|T|1/2,其中C為常數(shù),則對(duì)數(shù)后驗(yàn)密度

其中Vj=xjT+σ2In,dj=Yj-xjwjγ。為了簡(jiǎn)單起見,設(shè)定平衡的數(shù)據(jù)結(jié)構(gòu),每組樣本的個(gè)數(shù)相等,均為n。

求lnf(θ|y)對(duì)γ,σ2,T的偏導(dǎo)數(shù),計(jì)算中令φ =(τ00,τ10,τ11,σ2)T,根據(jù)Magnus或周杰對(duì)矩陣微積分的定義[3]89-133;[4]87-103,可以推出

因?yàn)?/p>

所以

三、后驗(yàn)密度最大化的EM算法

在模型式(1)中,將βj看作缺失部分,θ和βj的聯(lián)合后驗(yàn)條件密度函數(shù)為f(βj,θ|Yj)=p(θ)f1(Yj|βj,θ)f2(βj|θ)/f3(Y)∝p(θ)f1(Yj|βj,θ)f2(βj|θ),rj和uj都服從正態(tài)分布,所以θ和β1,β2,…,βN的聯(lián)合后驗(yàn)密度函數(shù)為:

令p(θ)=C|T|1/2,其中C為常數(shù),則對(duì)數(shù)后驗(yàn)密度為:

求lnf對(duì)γ,σ2,T的偏導(dǎo)數(shù),并令其為零:

可以算出

^γ,^σ2,^T的估計(jì)值均依賴βj,而βj是缺失數(shù)據(jù),因此用EM算法實(shí)現(xiàn)計(jì)算。

EM算法能夠處理帶有缺失數(shù)據(jù)的模型估計(jì)[5],對(duì)于多峰分布,EM算法會(huì)收斂到觀測(cè)數(shù)據(jù)似然函數(shù)的局部最大值;但是對(duì)于指數(shù)族分布,EM算法是相當(dāng)好的,每一次迭代都會(huì)使觀測(cè)變量的似然值增加,并達(dá)到全局最優(yōu)[6]。在EM算法里,觀測(cè)變量的似然函數(shù)的最大化可以轉(zhuǎn)換為后驗(yàn)似然函數(shù)的最大化,EM算法可以修改用以計(jì)算最大后驗(yàn)估計(jì),進(jìn)行貝葉斯推斷,在指數(shù)族分布下,能夠和以往一樣有效地收斂[7]64-85。

1.E步

為了找到E(βj|yj)和V ar(βj|yj),利用下面的多元正態(tài)分布。因?yàn)?/p>

其中

式(9)的具體推導(dǎo)過程如下:

首先由式(8),可以推出E(^γ|Yj),E(^σ2|Yj)和E(^T|Yj)。

然后用E(γ^|Yj),E(σ^|Yj)和E(^T|Yj)分別替換γ^, σ^2和^T。

EM的計(jì)算過程為:①定初值γ(0),σ2(0),T(0);②將初值代入E步,計(jì)算;③再將的值代入M步,計(jì)算γ(1),σ2(1),T(1);④重復(fù)②～③,直至收斂。

四、參數(shù)估計(jì)值方差的計(jì)算

EM算法簡(jiǎn)單易行,但是參數(shù)估計(jì)的方差不能夠直接從EM的迭代運(yùn)算中產(chǎn)生,下面運(yùn)用對(duì)數(shù)似然二階導(dǎo)數(shù),計(jì)算方差。將式(4)和式(5)分別再對(duì)r和φ求偏導(dǎo),

其中F=(0 0 0 1)。估計(jì)量的方差為

其中^γk表示γ的第k個(gè)元素,diag(A)k表示方陣A的對(duì)角線上第k個(gè)元素,diag(A)last表示方陣A的對(duì)角線上最后一個(gè)元素。

五、數(shù)據(jù)模擬

StephenW.Raudenbush和AnthonyS. Bryk[8]428-429給出了兩層HLM的MLE的具體步驟,雖然REML和Bayes估計(jì)在小樣本時(shí)改進(jìn)了模型的估計(jì),但MLE的估計(jì)量仍可以作為一個(gè)參照的標(biāo)準(zhǔn)估計(jì)。下面將從樣本量的變化和初值設(shè)定的變化來對(duì)最大后驗(yàn)估計(jì)方法進(jìn)行數(shù)據(jù)模擬,并將結(jié)果和MLE的估計(jì)進(jìn)行比較。

(一)模擬樣本的產(chǎn)生

假定模型式(1)的一個(gè)較簡(jiǎn)單形式。從N(0,1)中隨機(jī)生成x和w的樣本,樣本量分別為n3N(n =30,50,100為組內(nèi)個(gè)體數(shù),N=30,50為組數(shù));σ設(shè)為0.7,從N(0,0.72)中隨機(jī)生成r的樣本;為了計(jì)算方便,令τ10=0,τ11=τ00=0.3,隨機(jī)生成u0和u1的樣本;令γ00=1,γ10=γ01=γ11=0.3,由模型(1)得到y(tǒng)的樣本值。以下表格數(shù)據(jù)均為500次模擬計(jì)算的平均值。

(二)MAPE估計(jì)結(jié)果以及和MLE的比較

考慮到樣本量的大小會(huì)影響估計(jì)的結(jié)果,本研究用不同的組數(shù)30,50,100和個(gè)體數(shù)30,50的不同組合進(jìn)行數(shù)據(jù)模擬,每個(gè)樣本量下模擬500次。將MAPE和MLE估計(jì)中的迭代的初值都設(shè)定為1),結(jié)果參見表1～3。

表1 固定效應(yīng)和方差σ的MAPE和ML E及其方差估計(jì)值的比較表

表2 方差協(xié)方差成分的MAPE和ML E及其方差估計(jì)值的比較表

表1列出了MAPE和MLE兩種估計(jì)方法下,固定效應(yīng)參數(shù)γ00至γ11以及σ在不同樣本量下的估計(jì)值和標(biāo)準(zhǔn)差(表中用括號(hào)表示)?？梢郧宄乜吹絻煞N估計(jì)方法對(duì)固定效應(yīng)的系數(shù)和方差σ2的估計(jì)值在不同的樣本量下是一致的,并且500次模擬的平均值和真實(shí)值的絕對(duì)差都在0.006以內(nèi);MAPE估計(jì)量的標(biāo)準(zhǔn)差稍大一點(diǎn),兩者的絕對(duì)差在0.004以內(nèi),差距很小。

表2描述了在不同樣本量下,兩種方法對(duì)方差協(xié)方差成分τ00,τ11的估計(jì)情況。在(N,n)=(30, 30)時(shí),τ00,τ11的MLE估計(jì)值與真實(shí)值絕對(duì)差為01015和0.013,MAPE為0.008和0.005,在(N,n) =(30,50)時(shí),τ00,τ11的ML E估計(jì)值與真實(shí)值絕對(duì)差為0.011和0.014,MAPE為0.003和0.006,當(dāng)組數(shù)和總的樣本數(shù)增加時(shí)ML E的絕對(duì)差也能減小到0.01以下。這說明組數(shù)少、總樣本量小的時(shí)候,τ00, τ11的MAPE的估計(jì)比MLE的估計(jì)更接近真實(shí)值,當(dāng)組數(shù)增加,總樣本量增加時(shí)τ00,τ11的估計(jì)是一致的,并且都靠近真實(shí)值。另外,與固定效應(yīng)的估計(jì)一樣,估計(jì)量方差的MAPE稍大一點(diǎn),但也隨著樣本量的增加而趨于相同。

表3 MAPE和ML E的均方誤差和迭代次數(shù)的比較表

表3比較了兩種方法估計(jì)的均方誤(MSE)和計(jì)算的迭代次數(shù)。在相同的樣本的情況下,MAPE和ML E的MSE是一致的,并且模擬顯示它們都隨樣本的增加而減少。對(duì)迭代次數(shù)而言,ML E的迭代次數(shù)明顯高于MAPE,而且MAPE的迭代次數(shù)不受樣本量的影響,收斂得很快。

下面將減小樣本量繼續(xù)進(jìn)行模擬比較,樣本量組數(shù)N和組內(nèi)個(gè)體數(shù)n分別為:N為5,10和20;n為10,20。比較結(jié)果見表4～5。表4描述的模擬結(jié)果和表2類似?？梢郧宄乜吹疆?dāng)組數(shù)減小到5,組內(nèi)個(gè)體數(shù)減小到10時(shí),τ00的MLE的估計(jì)值為0.18,而MAPE為0.293,和真值0.3相比,MAPE的估計(jì)值比ML E更靠近真實(shí)值。隨著組數(shù)的增加,當(dāng)組數(shù)超過50時(shí),兩種估計(jì)值都趨于真實(shí)值。

隨著組數(shù)的減少,MLE的迭代次數(shù)明顯增加。表5顯示出小樣本下的兩種方法的模擬的MSE值和迭代次數(shù)。兩種估計(jì)方法的MSE仍是一致的,但是MAPE估計(jì)的迭代次數(shù)明顯低于MLE的迭代次數(shù)。

表4 更小樣本下方差協(xié)方差成分的MAPE和ML E及其方差估計(jì)值的比較表

表5 更小樣本下MAPE和ML E的均方誤差和迭代次數(shù)的比較表

六、結(jié) 論

本文推導(dǎo)出一個(gè)適合分層線性模型進(jìn)行最大后驗(yàn)估計(jì)的EM算法,并利用似然函數(shù)的二階導(dǎo)數(shù)推導(dǎo)出參數(shù)估計(jì)的方差估計(jì)值。在數(shù)據(jù)模擬中,將新的EM算法的結(jié)果同極大似然估計(jì)的EM算法的結(jié)果進(jìn)行比較,得出的結(jié)論是:1.對(duì)于固定效應(yīng)r和方差σ估計(jì)值,兩種方法是完全一致的。2.組數(shù)較少,樣本量較小時(shí),EM計(jì)算的MAPE的方差協(xié)方差成分的估計(jì)比用EM計(jì)算的ML E的估計(jì)更靠近真實(shí)值。3.EM計(jì)算MAPE迭代次數(shù)小于ML E;尤其當(dāng)樣本量較小時(shí),兩者之間的差距很明顯。

[1] Lindley D V,Smith A F M.Bayes Estimates for Linear Model[J].Journal of the Royal Statistical Society,Series B,1972,34.

[2] Jeffreys H.An Invariant Form for the Prior Probability in Estimation Problems[J].Proceedings of the Royal Society of London,Series A,Mathematical and Physical Sciences,1946,186.

[3] Magnus J R,Neudecker H.Matrix Differential Calculus with Applications in Statistics and Econometrics[M].New York:John Wiley,1988.

[4] 周杰.矩陣分析及應(yīng)用[M].成都:四川大學(xué)出版社,2008.

[5] Dempster A P,Laird N M,Rubin D B.Maximum Likelihood from Incomplete Data Via the EM Algorithm[J].Journal of the Royal Statistical Society,Series B,1977,39.

[6] Wu C F,Jeff.On the Convergence Properties of the EM Algorithm[J].Annals of Statistics,1983,11.

[7] Tanner M A.Tools for Statistical Inference[M].Heidelberg:Springer-Verlag,1996.

[8] Bryk A S,Raudenbush S W.Hierarchical Linear Models[M].Newbury Park,CA:Sage,1992.

Maximum a Posteriori Estimate of Hierarchical Linear Model

ZHANG Xuan
(School of Statistics,Renmin University of China,Beijing 100084,China)

Maximum a posteriori(MAP)estimation is an important method in statistical estimation. This research develops MAP estimation for Hierarchical Linear Model(HLM),and gives a detailed process of estimation and computational program.The expectation maximum(EM)algorithm was adapted to the concrete calculation,and the second derivatives of log likelihood function help us to obtain the esti2 mators of variance and covariance components in this model according to asymptotic distribution theory of maximum likelihood estimator.At last,we compare the results of MAP estimate to the general maximum likelihood(ML)estimate,noting that EM algorithm is exploited into the two different estimates.Simula2 tion shows that the estimators of fixed effects using the two methods are almost consistent,but the MAP estimators of variance and covariance components are closer to real values than the ML estimate counter2 part when the group size is small,at the same time,the iterations needed in MAP estimate are obviously less than those in ML estimate.

hierarchical linear model;maximum a posteriori estimate;maximum likelihood estimate; expectation maximum algorithm

F224.7

A

1007-3116(2011)01-0010-06

(責(zé)任編輯:馬 慧)

2010-08-29

張 璇,女,湖南湘潭人,講師,博士生,研究方向:統(tǒng)計(jì)模型及其計(jì)算、計(jì)量經(jīng)濟(jì)學(xué)。