分層線性模型的最大后驗估計

張 璇

分層線性模型的最大后驗估計

張 璇

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

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

分層線性模型:最大后驗估計:最大似然法:期望最大化算法

一、研究背景

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

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

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

為了清楚地描述整個估計過程,運用一個簡單的兩層的分層模型。

上述模型的矩陣形式是:

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

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

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

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

因為

所以

三、后驗密度最大化的EM算法

在模型式(1)中,將βj看作缺失部分,θ和βj的聯(liá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)合后驗密度函數(shù)為:

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

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

可以算出

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

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

1.E步

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

六、結(jié) 論

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

[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īng)濟學(xué)。