輪廓似然函數及其應用*

韓 棟 陳 征△ 陳平雁 Nakamura Tsuyoshi

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

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

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

原理與方法

1．輪廓似然定義

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

2．輪廓似然置信區間

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

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

實 例

1．多個多余參數出現的問題

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

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

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

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

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

2．偏態分布的輪廓似然置信區間

(1)數值模擬

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

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

表1 輪廓似然置信區間與Wald置信區間的置信水平

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

由表1和圖1可以看出，隨著試驗次數增大，Wald CI與PL CI趨于一致，PL CI也逐漸趨于對稱。試驗次數較小時(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值的區間為0.15～0.36，由于兩個CI上限間的RPL值均大于兩個CI下限間的RPL值，根據似然理論以及似然比檢驗的原理，Wald CI下限至PL CI下限間包括真實值的可能性均比Wald CI上限至 PL CI上限間包括真實值的可能性要低。圖1-B的結論與此類似，因此PL CI置信區間更可信。

(2)白鼠毒性實驗

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

表2 白鼠毒性實驗數據

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

表3 似然比檢驗與Wald檢驗

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

討 論

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

7．江一濤，胡海蘭，魏巧玲，等．競爭風險模型的發展與應用．中國衛生統計，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．