輪廓似然函數(shù)及其應用*

韓 棟 陳 征△ 陳平雁 Nakamura Tsuyoshi

在回歸分析中，似然函數(shù)通常會含有多個參數(shù)，但有時只有其中一個或幾個是欲研究的參數(shù)，稱為興趣參數(shù)(parameter of interest)，其他參數(shù)就被稱作多余參數(shù)(nuisance parameter)，這些多余參數(shù)對模型的求解有時會有阻礙作用。當存在多個多余參數(shù)時，標準的似然方法無法消除或減少它們，所以變得不可靠或完全無效，而輪廓似然(profile likelihood，PL)作為一種處理多余參數(shù)的方法能夠解決多余參數(shù)過多的問題。1970年，Kalbfleisch和 Sprott等〔1〕首次將輪廓似然方法應用于帶有多余參數(shù)的參數(shù)推斷，并稱最大輪廓似然函數(shù)為最大相對似然函數(shù)(maximum relative likelihood function)。Barndorff-Nielsen〔2〕首先使用“輪廓似然”命名該方法，之后該名字被廣泛接受〔3〕。2000年，Murphy等〔4〕證明了在一般情況下最大輪廓似然點估計等價于最大似然估計。

另外，在興趣參數(shù)呈非正態(tài)分布時，如果計算基于正態(tài)分布的Wald型置信區(qū)間(Wald CI)將會產(chǎn)生偏差〔5〕，尤其在無法計算興趣參數(shù)的標準誤時，Wald CI也無法計算。而輪廓似然置信區(qū)間(PL CI)是基于χ2分布且無需計算標準誤，因此，PL CI能夠解決參數(shù)不服從正態(tài)分布和標準誤無法計算時置信區(qū)間的計算問題。Venzon等〔5〕于1988年提出了簡化輪廓似然置信區(qū)間計算的一種新方法。

本文將描述輪廓似然的定義及其兩個應用，模擬比較PL CI與Wald CI的優(yōu)劣并運用PL方法解決多余參數(shù)過多和參數(shù)呈非正態(tài)分布時的問題。

原理與方法

1．輪廓似然定義

輪廓似然函數(shù)是固定興趣參數(shù)時，對多余參數(shù)求最大化后的函數(shù)，因而不是真正的似然函數(shù)。設θ表示興趣參數(shù)或興趣參數(shù)向量，γ表示多余參數(shù)或多余參數(shù)向量，假設X1，…，Xn為獨立同分布且密度函數(shù)為，然后輪廓似然函數(shù)被定義為pl(θ)=l［θ，^γ(θ)］，其中，^γ(θ)為固定θ時，γ的最大似然估計值(MLE)，即:pl(θ)=maxγl(θ，γ)。

2．輪廓似然置信區(qū)間

Wald CI是根據(jù)一個預先給定的置信水平和參考分布(在線性回歸分析中選用t分布，其他為標準正態(tài)分布)選定分位數(shù)，采用“估計值±分位數(shù)×估計值的標準誤”來計算模型中某個參數(shù)的置信區(qū)間。如果興趣參數(shù)的分布呈偏態(tài)分布或無法計算其標準誤時，Wald CI的結果不可靠，而PL CI對以上特殊情況并不敏感，是一種更加穩(wěn)健的方法。PL方法可應用于所有基于似然理論的統(tǒng)計分析。

興趣參數(shù)θ的95%PL CI是由檢驗水準為0.05時似然比檢驗無統(tǒng)計學意義的所有θ構成，即所有使似然比統(tǒng)計量小于等于3.84)的 θ值。用公式表示為滿足ln［pl(θ)］≥ln［pl(θ^)］－3.84/2=ln［pl(θ^)］－1.92的所有 θ值構成了95%PL CI，其中 θ^是θ的最大輪廓似然估計值。用代替 3.84 可以計算其他置信水平為100(1－α)%的置信區(qū)間。

實 例

1．多個多余參數(shù)出現(xiàn)的問題

在對2003年SARS病死率估計的研究中，陳征等〔6〕基于競爭風險理論〔7〕建立模型:令 ni、di、ci和 ai分別指代在第i點的新增患者、死亡人數(shù)、治愈康復人數(shù)和觀察人數(shù)，h1i、h2i分別表示死亡與治愈的危險率，其中i=1，…，s，表示不同時間點。根據(jù)實際數(shù)據(jù)觀察可假設治愈-死亡危險率比Ri=h2i/h1i≡R是一個常數(shù)，則病死率估計值為(1+R)－1。關于R和h1i的對數(shù)似然函數(shù)為:

因為病死率估計公式只與R有關，因此上式中R為興趣參數(shù)，其他參數(shù)(h1i，i=1，2，3，…，s)為多余參數(shù)，此時似然函數(shù)中有(s+1)個參數(shù)，而且隨著觀察時間點增多(s增大)，多余參數(shù)個數(shù)在不斷增加，因此不能直接使用標準最大似然估計求解參數(shù)。基于實際數(shù)據(jù)研究〔8〕及Lam〔9〕研究，模型又假設h1i≡h1為常數(shù)，從而將對數(shù)似然函數(shù)中的參數(shù)個數(shù)減至可求解的兩個(R和h1)。將每個時間點的兩個危險率均設為常數(shù)的條件過于苛刻，但無此假設無法使用MLE估計參數(shù)。

使用輪廓似然方法解決上述問題:

此處僅假設Ri為常數(shù)，即Ri=h2i/h1i≡R，基于似然函數(shù)公式(1)，解方程組 ?l/?h1i=0，得出 ^h1i=(di+ci)/［ai(1+R)］，然后將 ^h1i代替 h1i代入公式得出對數(shù)輪廓似然函數(shù):則R的近似方差估計是:

本例也驗證了Murphy的結論，即最大輪廓似然點估計(式(2)和(3))與MLE結果〔6〕一致。由于輪廓似然方法的假設相比MLE方法〔6〕的假設弱化了很多，因此當存在多余參數(shù)時，使用輪廓似然方法可以提高方法的適用性。

2．偏態(tài)分布的輪廓似然置信區(qū)間

(1)數(shù)值模擬

此節(jié)對不同偏態(tài)分布情況下PL CI和Wald CI的置信水平進行檢測。為了模擬非正態(tài)分布參數(shù)，選取logistic模型 log(pi/1－pi)=β1+β2xi，并設定 xi分別為(60，65，75，90)，β1= － 6.5，β2=0.1。采用二項分布，每個x下的試驗次數(shù)分別設定為3、8、20，以每一個pi為發(fā)生率，模擬出每個試驗次數(shù)下的事件發(fā)生次數(shù)與失敗次數(shù)，擬合logistic回歸模型并計算PL CI和Wald CI界值在χ2(1)分布下的置信水平。相對輪廓似然值(relative PL，RPL)定義為:輪廓似然值/最大輪廓似然值。根據(jù)似然理論，RPL表示數(shù)據(jù)對兩個參數(shù)估計值支持程度的比值，取值為(0，1］，因此可采用RPL比較不同數(shù)據(jù)情況下的置信限處的似然。輪廓似然不對稱性指標的計算公式〔11〕為:

表示置信限到估計值距離之差占置信區(qū)間長度的百分比，不對稱性越趨近于0，表示PL CI越趨于對稱。模擬結果反映在表1和圖1上。

表1 輪廓似然置信區(qū)間與Wald置信區(qū)間的置信水平

圖1 不同試驗次數(shù)下的相對輪廓似然值(左1－A，n=3，右1－B，n=20)

由表1和圖1可以看出，隨著試驗次數(shù)增大，Wald CI與PL CI趨于一致，PL CI也逐漸趨于對稱。試驗次數(shù)較小時(n=3)，PL CI不對稱性為28.9%，95%Wald CI的置信水平僅為93.0%，由于采用PL方法，95%PL CI的置信水平被控制在95.0%。

圖1-A中，Wald CI下限至PL CI下限間的RPL值在0.03～0.15之間，而Wald CI上限至PL CI上限間RPL值的區(qū)間為0.15～0.36，由于兩個CI上限間的RPL值均大于兩個CI下限間的RPL值，根據(jù)似然理論以及似然比檢驗的原理，Wald CI下限至PL CI下限間包括真實值的可能性均比Wald CI上限至 PL CI上限間包括真實值的可能性要低。圖1-B的結論與此類似，因此PL CI置信區(qū)間更可信。

(2)白鼠毒性實驗

利用PL來分析白鼠毒性實驗〔12〕，ni表示總的白鼠數(shù)，ri表示死亡鼠數(shù)，xi表示毒藥劑量，數(shù)據(jù)如下表:

表2 白鼠毒性實驗數(shù)據(jù)

對以上數(shù)據(jù)擬合logistic回歸模型:log(pi/(1－pi))=β1+β2log xi(i=1，…，4)。結果見表 3，經(jīng) Wald檢驗，毒藥劑量的對數(shù)值對白鼠的死亡率沒有影響(P=0.119)，但由圖2可以看出，β2的輪廓似然函數(shù)值呈正偏態(tài)，不對稱性達到41.2%，因此采用Wald法不可靠。如果采用似然比檢驗，由表3的結果顯示，毒藥劑量的對數(shù)值對白鼠的死亡率的影響有統(tǒng)計學意義(P＜0.001)，毒藥劑量對數(shù)值系數(shù)的 PL CI為(2.283，21.491)。

表3 似然比檢驗與Wald檢驗

圖2 白鼠毒性實驗中系數(shù)值的相對輪廓似然值

討 論

本文就輪廓似然方法及其應用進行了闡述，并用模擬與實例說明輪廓似然在估計參數(shù)值和計算置信區(qū)間等方面都有較強的實用性。除了文中所述的一些性質外，在參數(shù)模型中，對數(shù)輪廓似然函數(shù)的二階導函數(shù)是觀察信息量的估計值，甚至是在輪廓似然函數(shù)不能寫成外顯函數(shù)的情況下，數(shù)值計算方法也可以計算出信息矩陣的估計值。輪廓似然方法還有其他特殊的性質，如利用輪廓似然方法消去普通似然函數(shù)中的基準危險率，從而推導出擬合Cox回歸時使用的偏似然函數(shù)〔4〕;也可以利用輪廓似然方法消去基準危險率后，構造全輪廓似然函數(shù)〔13〕，在中小樣本情況下，最大全輪廓似然估計值比最大偏似然估計值更有用;與標準的似然方法相比，利用輪廓似然方法處理有刪失的生存時間數(shù)據(jù)時，無需對刪失類型進行假設〔14〕。除了輪廓似然方法外，處理多余參數(shù)的方法還有邊際似然、條件似然、聯(lián)合似然等。由于以上三種似然方法的使用都需要依賴一定的特殊結構，而本文所述的輪廓似然沒有這種限制，甚至在輪廓似然函數(shù)不能被寫成顯性函數(shù)的形式時，輪廓似然方法依然適用。因此輪廓似然作為一種處理多余參數(shù)的方法更可行〔15〕。

1．Kalbfleisch JD，Sprott DA．Application of likelihood methods to models involving large numbers of parameters．Journal of the Royal Statistical Society．Series B(Methodological)，1970:175-208．

2．Barndorff-Nielsen O．On a formula for the distribution of the maximum likelihood estimator．Biometrika，1983，70(2):343

3．Bjφrnstad JF．Predictive Likelihood．Encyclopedia of Statistical Sciences，2006，9:6369-6375．

4．Murphy SA，Van der Vaart AW．On profile likelihood．Journal of the A-merican Statistical Association，2000，95(450):449-465．

5．Venzon DJ，Moolgavkar SH．A method for computing profile-likelihoodbased confidence intervals．Applied Statistics，1988，37(1):87-94．

6．陳征，Nakamura T．基于競爭風險理論和概要型數(shù)據(jù)的病死率估計模型．中國衛(wèi)生統(tǒng)計，2010，27(3):249-252．

7．江一濤，胡海蘭，魏巧玲，等．競爭風險模型的發(fā)展與應用．中國衛(wèi)生統(tǒng)計，2009，26(4):445-447．

8．Chen Z，Nakamura T．Statistical evidence for the usefulness of Chinese medicine in the treatment of SARS．Phytotherapy Research，2004，18(7):592-594．

9．Lam KF，Deshpande JV，Lau E，et al．A test for constant fatality rate of an emerging epidemic:with applications to severe acute respiratory syndrome in Hong Kong and Beijing．Biometrics，2008，64(3):869-876．

10．Tsodikov A，Garibotti G．Profile information matrix for nonlinear transformation models．Lifetime data analysis，2007，13(1):139-159．

11．Royston P．Profile likelihood for estimation and confidence intervals．Stata Journal，2007，7(3):376-387．

12．Aitkin M．Statistical modelling:the likelihood approach．The Statistician，1986，35(2):103-113．

13．Ren J，Zhou M．Full likelihood inferences in the Cox model:an empirical likelihood approach．Annals of the Institute of Statistical Mathematics，2011，63(5):1005-1018．

14．Zhang Z．Profile likelihood and incomplete data．International Statistical Review，2010，78(1):102-116．

15．Montoya J，Díaz-Francés E，Sprott D．On a criticism of the profile likelihood function．Statistical Papers，2009，50(1):195-202．