王寧寧徐淑一方積乾

利用得分方程S（θ）在θ0處線性近似，有

?

從經(jīng)典似然到等級似然的理論概述和應(yīng)用*

王寧寧1徐淑一2△方積乾3

1.廣州醫(yī)科大學(xué)公共衛(wèi)生學(xué)院統(tǒng)計學(xué)系（511436）

2.中山大學(xué)嶺南學(xué)院經(jīng)濟(jì)學(xué)系

3.中山大學(xué)公共衛(wèi)生學(xué)院流行病與醫(yī)學(xué)統(tǒng)計系

Fisher（1921）［1］第一次使用了似然（likelihood）這個概念，以后有關(guān)似然理論得到了極大發(fā)展和廣泛應(yīng)用。由Lee和Nelder（1996）［2］提出的帶隨機(jī)效應(yīng)的廣義線性模型的等級似然（hierarchicacl likelihood）估計方法，在隨機(jī)效應(yīng)統(tǒng)計模型估計中得到了較大的發(fā)展和應(yīng)用。如LI.DA.HA等（2002）［3］、LI.DA.HA等（2005）［4］、徐淑一等（2007）［5］、王寧寧等（2011，2014）［6-7］將等級似然估計方法用于生存分析中脆弱模型的估計。而且近年來學(xué)術(shù)界關(guān)于等級似然方法的應(yīng)用研究仍在繼續(xù)，如Lee和Jiang等（2011）［8］使用等級似然方法估計預(yù)測疾病測繪時相對風(fēng)險的區(qū)間，模擬顯示其不弱于貝葉斯方法。Noh和Lee等（2011）［9］使用等級似然建立非線性混合效應(yīng)模型處理縱向數(shù)據(jù)。Noh和Lee等（2012）［10］使用等級似然方法提供了一種減弱對數(shù)據(jù)缺失機(jī)制不正確假定的影響并提供了實例和模擬研究。Lee等（2013）［11］把等級似然函數(shù)方法用于缺失數(shù)據(jù)的建模并應(yīng)用于厚尾分布的縱向數(shù)據(jù)分析。Wu和Bentler（2012）［12］使用等級似然解決因子分析模型中二元響應(yīng)變量的情形，使用簡單而且高效。徐淑一等（2009）［13］、HA等（2014）［14］和Lee等（2014）［15］將等級似然估計方法應(yīng)用于競爭風(fēng)險模型的估計，取得了良好的效果。

相比較于貝葉斯方法，等級似然估計方法在很多情形下相對較為簡潔，而且有較為廣泛的應(yīng)用范圍，因此，本文對這一理論進(jìn)行介紹。

經(jīng)典似然理論到擴(kuò)展似然

似然的定義是：假定一個統(tǒng)計模型f（y），其中θ是未知參數(shù)，L（θ）是觀測到的數(shù)據(jù)y的概率，將它看作是θ的函數(shù)，在這個意義下，L（θ）即為似然。似然函數(shù)的作用在于承載未知參數(shù)的信息，關(guān)于使用似然函數(shù)的推論一直以來都是有爭議的，但是，目前似然理論最廣泛的應(yīng)用仍然是基于似然函數(shù)推導(dǎo)的一些統(tǒng)計量和這些統(tǒng)計量的概率特性。Birnbaum（1962）［16］提出：任何一個通過試驗得到的結(jié)論都應(yīng)該來自于似然。

由（1）式可以得到：

如果數(shù)據(jù)是來自正態(tài)總體，那么二次逼近就是精確的。可以發(fā)現(xiàn)，對數(shù)似然的二階近似恰好就對應(yīng)于＾θ的正態(tài)近似。由（2）式直接可得：

利用得分方程S（θ）在θ0處線性近似，有

（5）式就是Barndorff-Nielsen的p-formula（1983）［17］，已經(jīng)被證實這個近似比基于的正態(tài)近似的密度公式更加精確，而（5）式也正是進(jìn)一步向輪廓似然和等級似然拓展的基礎(chǔ)。

當(dāng)模型的參數(shù)比較多時，我們可能只對一部分參數(shù)感興趣。例如在正態(tài)模型里，我們可能只對均值μ感興趣，此時σ2就是一個討厭參數(shù)（nuisance parameter）。一個方法就是關(guān)于感興趣的參數(shù)集中化似然，這就是輪廓似然（profile likelihood）。未知的討厭參數(shù)會帶來的額外的不確定性，尤其是在小樣本下，考慮這個額外的信息是非常重要的。假定（θ，η）是所有的參數(shù)，其中θ是感興趣的參數(shù)，給定聯(lián)合似然函數(shù)L （θ，η），則θ的輪廓似然函數(shù)定義為：

其中上面的最大化是固定θ進(jìn)行的。固定θ以后，得到的η的估計是一個與θ有關(guān)的函數(shù)，記做L但輪廓似然并不是一個正確的似然函數(shù)，它并不是基于觀測數(shù)據(jù)的概率。對于大部分問題，精確的邊際似然和條件似然雖然理論上存在，實際卻很難得到。設(shè)是（θ，η）的極大似然估計令是固定θ時η的極大似然估計，是相應(yīng)的觀測Fisher信息，那么的逼近密度為

于是，關(guān)于θ的條件似然為：

等級似然理論

對于更一般的包括不可觀測的隨機(jī)效應(yīng)的模型，F(xiàn)isher似然需要拓展。一直以來，包括不可觀測的隨機(jī)變量的模型的似然的定義，統(tǒng)計學(xué)家們沒有達(dá)成一致。Bj?rnstad（1996）［18］建立了擴(kuò)展似然（extended likelihood）理論，表明一般似然的一個特殊定義包含了固定的和隨機(jī)的參數(shù)的所有信息。Lee和Nelder （1996）［2］介紹了等級GLM模型的等級似然與傳統(tǒng)的似然有很大不同。關(guān)于似然一個重要的特性是對變換的不變性，這也是區(qū)分?jǐn)U展的似然與等級似然的重要方面，盲目的最優(yōu)化擴(kuò)展似然可能會使估計缺乏不變性，參數(shù)的不同的尺度（形式）會導(dǎo)致不同的估計，依賴于參數(shù)的尺度使得擴(kuò)展似然面臨一些批評，而這個問題，等級似然則通過定義一個特殊的參數(shù)的尺度而解決。

關(guān)于對似然的拓展，很多學(xué)者都做了嘗試，如Lauritzen（1974）［19］、Butler（1986）［20］以及Bj?rnstad （1996）［18］等。對似然合理的拓展應(yīng)該能夠處理未知參數(shù)θ，不可觀測的隨機(jī)效應(yīng)v和觀測數(shù)據(jù)y，就數(shù)據(jù)生成過程來說，隨機(jī)效應(yīng)v來自于概率函數(shù)fθ（v），當(dāng)v生成之后，固定v，由概率分布函數(shù)fθ（y|v）產(chǎn)生隨機(jī)數(shù)據(jù)y。那么聯(lián)合的隨機(jī)模型就是：

給定數(shù)據(jù)y，可以利用邊際似然L（θ；y）對θ進(jìn)行推理；給定θ之后，再利用條件似然L（θ，v；v|y）=fθ（v |y）對隨機(jī)效應(yīng)v進(jìn)行推理。關(guān)于（θ，v）的擴(kuò)展似然定義為：

兩個過程的聯(lián)系通過L（θ，v；y，v）=fθ（y，v）給出。擴(kuò)展似然的定義式中，左邊是y固定，（θ，v）變動，右邊是θ固定，（y，v）變動。在擴(kuò)展似然框架下，v作為隨機(jī)實現(xiàn)出現(xiàn)在數(shù)據(jù)生成過程中，在參數(shù)估計中又作為未知參數(shù)待估計。而在經(jīng)典似然框架下，僅有一類隨機(jī)數(shù)據(jù)y。在兩類的參數(shù)（固定參數(shù)和隨機(jī)參數(shù)）都未知的情況下，擴(kuò)展似然理論沒有告訴我們應(yīng)該如何推理每一部分的參數(shù)。

記擴(kuò)展的似然為le（θ，v）=logL（θ，v；y，v），經(jīng)典似然記做l（θ）=logL（θ；y）。由擴(kuò)展似然的定義顯然有：le（θ，v）=l（θ）+ logfθ（v|y），對于固定參數(shù)采用l（θ）是經(jīng)典似然方法，對隨機(jī)效應(yīng)參數(shù)采用logfθ（v|y）是經(jīng)驗貝葉斯方法。由數(shù)據(jù)生成的過程邊際似然可以經(jīng)由積分得到，于是：

然而對于非正態(tài)模型，這種積分很難處理。一種獲得固定參數(shù)θ的邊際極大似然估計的方法是EM算法。EM算法的E步需要E（le）在給定y條件下的解析式，在M步將其最大化。眾所周知EM算法收斂得很慢，對于非正態(tài)模型，E步的期望很難計算，為此一些模擬算法，比如蒙特卡羅EM方法（Vaida和Meng 2004）［21］和Gibbs抽樣（Gelfand和Smith 1990）［22］等都可以用來計算條件期望，但是這些方法的計算都很復(fù)雜。另外的經(jīng)由Gauss-Hermite求積（Crouch和Spiegelman 1990）［23］的數(shù)值積分方法可以直接獲得ML估計，但是當(dāng)隨機(jī)效應(yīng)參數(shù)增加時，計算會相當(dāng)繁重。

在擴(kuò)展似然框架下，對固定參數(shù)的正確推理要使用由le（θ，v）積分掉隨機(jī)效應(yīng)參數(shù)獲得的邊際似然l （θ），由于積分難以處理，邊際似然可以采用拉普拉斯近似：

在估計隨機(jī)效應(yīng)參數(shù)往往叫做隨機(jī)效應(yīng)的預(yù)測，Nelder認(rèn)為應(yīng)該稱之為估計，因為一旦數(shù)據(jù)y已經(jīng)生成，此時隨機(jī)效應(yīng)就是固定的，關(guān)于隨機(jī)效應(yīng)的最優(yōu)無偏預(yù)測也可以認(rèn)為是最優(yōu)無偏估計。令θ1和θ2分別是固定參數(shù)θ的任意兩個值，這兩個參數(shù)的信息都包括在似然比L（θ1；y）/ L（θ2；y）中，設(shè)存在一個隨機(jī)效應(yīng)的尺度v，使θ1和θ2的似然比保持不變：

如果存在上述定義的規(guī)范尺度的參數(shù)v，那么就可以馬上得到等級似然，然而，并不是所有的統(tǒng)計問題都存在規(guī)范參數(shù)，因此需要拓展它的定義。事實上，等級似然的定義正是擴(kuò)展似然的一個特殊情形，也就是說，當(dāng)v是規(guī)范的時候，擴(kuò)展似然L（θ，v；y，v）就是等級似然，下面用H（θ，v）表示等級似然，用h（θ，v）表示對數(shù)等級似然。對數(shù)等級似然可以被看作通常的對數(shù)似然函數(shù)，可以對它關(guān)于固定參數(shù)和隨即參數(shù)（θ，v）一起求導(dǎo)計算相應(yīng)的Fisher信息，在一般的統(tǒng)計問題中，什么樣的隨機(jī)效應(yīng)參數(shù)應(yīng)該是規(guī)范尺度的隨機(jī)變量參數(shù)不是很明顯，然而檢驗一個特定尺度是否是規(guī)范的卻很容易。如果規(guī)范尺度的隨機(jī)效應(yīng)參數(shù)存在，關(guān)于參數(shù)的推理將會大大簡化。

關(guān)于這些性質(zhì)的進(jìn)一步討論可以參考Lee，Nelder 和Pawitan（2006）［24］。假設(shè)v是規(guī)范尺度的，非線性變換u≡v（u）將擴(kuò)展似然變?yōu)椋?/p>

由于雅克比項|J（u）|的存在，u不是規(guī)范尺度的。也就是說，在線性變換意義下，規(guī)范尺度是唯一的。

上述等級似然的定義比較嚴(yán)格，規(guī)范尺度是關(guān)于所有固定參數(shù)定義的，需要進(jìn)一步擴(kuò)展。假定固定參數(shù)分為兩部分（θ，φ），成立：

這種情形下，尺度v僅僅是與θ的信息無關(guān)，而關(guān)于φ則不然。使用等級似然的推理僅適用于（θ，v），而關(guān)于φ的推理則需要邊際似然，關(guān)于φ的邊際似然可以使用調(diào)整的輪廓似然近似：

其中，D（h，v）是φ的函數(shù)，與θ無關(guān)。至此，完成了模型中所有參數(shù)（θ，v，φ）的估計。

Lee等（2006）［24］指出，在廣義混合線型模型中，如果隨機(jī)效應(yīng)和協(xié)變量在線性預(yù)測部分中可加，那么這種隨機(jī)效應(yīng)尺度就是弱規(guī)范尺度，弱規(guī)范尺度的隨機(jī)效應(yīng)總可以找到；而且，此時協(xié)變量系數(shù)β的估計和隨機(jī)效應(yīng)v的實現(xiàn)可以通過最大化等級似然獲得，等級似然函數(shù)二階導(dǎo)數(shù)的相反數(shù)之逆仍然給出了協(xié)變量系數(shù)和隨機(jī)效應(yīng)實現(xiàn)估計值的方差。上面的這些性質(zhì)為等級似然的應(yīng)用提供了可行的理論基礎(chǔ)。

應(yīng)用舉例

我們以競爭風(fēng)險下比例危險模型為例說明等級似然估計的應(yīng)用。我們假定兩個競爭風(fēng)險情形。設(shè)觀測數(shù)據(jù)為Tij（i =1，2，…，q，j =1，2，…，ni）表示第i個體的第j個重復(fù)觀測。用Δij=（δij1，δij2）表示第ij個觀測的時間屬性，δij1或δij2為1表示事件發(fā)生，為0表示刪失。設(shè)二元隨機(jī)效應(yīng)變量V =（V1，V2），假設(shè)在給定協(xié)變量X以及V的條件下，T1，T2相互獨立；第ij個個體在不同風(fēng)險下的危險函數(shù)設(shè)為半?yún)?shù)Cox比例危險模型為：

對第ij個個體，可以定義競爭風(fēng)險下比例危險模型的等級似然函數(shù)：它是可觀測的生存時間和隨機(jī)效應(yīng)因子的聯(lián)合對數(shù)密度函數(shù)。記y（k）（k =1，2，…，K）為第k個觀測的持續(xù)期，y（1）＜y（2）＜…＜y（K）?？梢缘玫剑?/p>

其中，R（y（k））表示y（k）時刻的危險集，是y（k）時刻分別因為兩種風(fēng)險事件發(fā)生的數(shù)目，ηij1=xijβ1+ vi1，ηij2=xijβ2+ vi2。如果第ij個觀測在y（k）時刻因第一種風(fēng)險事件發(fā)生，則令否則令0；同樣定義記表示K個時刻第一種風(fēng)險的基本危險率向量，其中l(wèi)og（λ1（y（k）））；同樣定義記exp（ηij2），則

其中（i，j）∈R（y（k）），且parameters）］，并不影響上述結(jié)論。

通過上述討論，發(fā)現(xiàn)競爭風(fēng)險下的Cox脆弱性比例危險模型可以納入混合廣義線性模型的框架，也就是說，根據(jù)Lee等（2006）［24］，競爭風(fēng)險下Cox比例危險模型中隨機(jī)效應(yīng)尺度是弱規(guī)范尺度，從而可以利用第三部分討論的等級似然估計方法和程序?qū)﹄S機(jī)效應(yīng)競爭風(fēng)險模型進(jìn)行估計。顯然，經(jīng)典脆弱性Cox比例危險模型也同樣可以納入混合廣義線性模型框架，利用等級似然理論，可以擴(kuò)展Cox脆弱性模型中隨機(jī)效應(yīng)的分布。

小 結(jié)

本文對近年來模型估計理論中的一個新的熱點問題進(jìn)行了介紹和討論：即等級似然理論及其應(yīng)用。本文從經(jīng)典的似然理論開始，引出其擴(kuò)展，接著介紹了當(dāng)模型的參數(shù)比較多時，估計感興趣參數(shù)的輪廓似然方法，由此過渡到隨機(jī)效應(yīng)模型的等級似然估計方法。在上述基礎(chǔ)上，本文介紹了等級似然估計對于模型中三類參數(shù)的估計方法和程序，即模型中關(guān)注變量的系數(shù)、隨機(jī)效應(yīng)的實現(xiàn)值、隨機(jī)效應(yīng)分布參數(shù)。最后，本文對等級似然理論的應(yīng)用進(jìn)行舉例說明。等級似然理論相對于貝葉斯理論的應(yīng)用較為簡潔，對于隨機(jī)效應(yīng)廣義線性的估計而言，是一種有效的估計方法。對于等級似然估計方法，期待未來更多地應(yīng)用領(lǐng)域。

參考文獻(xiàn)

［1］Fisher RA.On the“Probable Error”of a Coefficient of Correlation Deduced from a Small Sample.Metron，1921，1：3-32.

［2］Lee Y，Nelder JA.Hierarchical generalized linear models.Journal of the Royal Statistical Society.Series B（Methodological），1996，619-678.

［3］Do Ha I，Lee Y，Song JK.Hierarchical-likelihood approach for mixed linear models with censored data.Lifetime data analysis，2002，8（2）：163-176.

［4］Do Ha I，Lee Y.Multilevel mixed linear models for survival data.Lifetime data analysis，2005，11（1）：131-142.

［5］徐淑一，王寧寧.競爭風(fēng)險下縱列數(shù)據(jù)的隨機(jī)效應(yīng)建模和估計.中山大學(xué)學(xué)報（自然科學(xué)版），2007，46（1）：7-10.

［6］Wang N，Xu S，F(xiàn)ang J.Hierarchical likelihood approach for the Weibull frailty model.Journal of Statistical Computation and Simulation，2011，81（3）：343-356.

［7］王寧寧，徐淑一，方積乾.重復(fù)觀測生存數(shù)據(jù)的AR（1）隨機(jī)效應(yīng)建模和估計.中國衛(wèi)生統(tǒng)計，2014，31（6）：1-4.

［8］Lee Y，Jang M，Lee W.Prediction interval for disease mapping using hierarchical likelihood.Computational Statistics，2011，26（1）：159-179.

［9］Noh M，Lee Y，Kenward MG.Robust estimation of dropout models using hierarchical likelihood.Journal of Statistical Computation and Simulation，2011，81（6）：693-706.

［10］Noh M，Wu L，Lee Y.Hierarchical likelihood methods for nonlinear and generalized linear mixed models with missing data and measurement errors in covariates.Journal of Multivariate Analysis，2012，109：42-51.

［11］Lee D，Lee Y，Paik MC，et al.Robust inference using hierarchical likelihood approach for heavy-tailed longitudinal outcomes with missing data：An alternative to inverse probability weighted generalized estimating equations.Computational Statistics &Data Analysis，2013，59（0）：171-179.

［12］Wu J，Bentler PM.Application of H-likelihood to factor analysis models with binary response data.Journal of Multivariate Analysis，2012，106：72-79.

［13］徐淑一，王寧寧，王美今.競爭風(fēng)險下縱列持續(xù)數(shù)據(jù)隨機(jī)效應(yīng)模型的估計與模擬研究.數(shù)理統(tǒng)計與管理，2009，28（6）：1013-1023.

［14］Ha ID，Lee M，Oh S，et al.Variable selection in subdistribution hazard frailty models with competing risks data.Statistics in medicine，2014，33（26）：4590-4604.

［15］Lee M，Do Ha I，Lee Y.Frailty modeling for clustered competing risks data with missing cause of failure.Statistical methods in medical research，2014，0962280214545639.

［16］Birnbaum A.On the foundations of statistical inference.Journal of the American Statistical Association，1962，57（298）：269-306.

［17］Barndorff-Nielsen O.On a formula for the distribution of the maximum likelihood estimator.Biometrika，1983，70（2）：343-365.

［18］Bj?rnstad JF.On the generalization of the likelihood function and the likelihood principle.Journal of the American Statistical Association，1996，91（434）：791-806.

［19］Lauritzen SL.Sufficiency，prediction and extreme models.Scandinavian Journal of Statistics，1974：128-134.

［20］Butler RW.Predictive likelihood inference with applications.Journal of the Royal Statistical Society.Series B（Methodological），1986：1-38.

［21］Vaida F，Meng XL，Xu R.Mixed effects models and the EM algorithm.Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives：An Essential Journey with Donald Rubin′s Statistical Family，2004：253-264.

［22］Gelfand AE，Smith AF.Sampling-based approaches to calculating marginal densities.Journal of the American Statistical Association，1990，85（410）：398-409.

［23］Crouch EA，Spiegelman D.The Evaluation of Integrals of the form∫+∞-∞f（t）exp（-t2）dt：Application to logistic-Normal Models.Journal of the American Statistical Association，1990，85（410）：464-469.

［24］Lee Y，Nelder JA，Pawitan Y.Generalized linear models with random effects：unified analysis via H-likelihood.CRC Press，2006.

［25］Lee Y，Nelder JA.Likelihood for random-effect models.SORT，2005，29（2）：141-164.

（責(zé)任編輯：郭海強(qiáng)）

·短篇報道·

*基金項目：廣州醫(yī)科大學(xué)項目（L135021）；廣東省軟科學(xué)項目（2013B070206027）

通信作者：△徐淑一