韓 棟 陳 征△ 陳平雁 Nakamura Tsuyoshi

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

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

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

原理與方法

1．輪廓似然定義

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

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

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

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

實(shí) 例

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

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

因?yàn)椴∷缆使烙?jì)公式只與R有關(guān)，因此上式中R為興趣參數(shù)，其他參數(shù)(h1i，i=1，2，3，…，s)為多余參數(shù)，此時(shí)似然函數(shù)中有(s+1)個(gè)參數(shù)，而且隨著觀察時(shí)間點(diǎn)增多(s增大)，多余參數(shù)個(gè)數(shù)在不斷增加，因此不能直接使用標(biāo)準(zhǔn)最大似然估計(jì)求解參數(shù)?；趯?shí)際數(shù)據(jù)研究〔8〕及Lam〔9〕研究，模型又假設(shè)h1i≡h1為常數(shù)，從而將對(duì)數(shù)似然函數(shù)中的參數(shù)個(gè)數(shù)減至可求解的兩個(gè)(R和h1)。將每個(gè)時(shí)間點(diǎn)的兩個(gè)危險(xiǎn)率均設(shè)為常數(shù)的條件過于苛刻，但無此假設(shè)無法使用MLE估計(jì)參數(shù)。

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

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

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

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

(1)數(shù)值模擬

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

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

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

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

由表1和圖1可以看出，隨著試驗(yàn)次數(shù)增大，Wald CI與PL CI趨于一致，PL CI也逐漸趨于對(duì)稱。試驗(yàn)次數(shù)較小時(shí)(n=3)，PL CI不對(duì)稱性為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，由于兩個(gè)CI上限間的RPL值均大于兩個(gè)CI下限間的RPL值，根據(jù)似然理論以及似然比檢驗(yàn)的原理，Wald CI下限至PL CI下限間包括真實(shí)值的可能性均比Wald CI上限至 PL CI上限間包括真實(shí)值的可能性要低。圖1-B的結(jié)論與此類似，因此PL CI置信區(qū)間更可信。

(2)白鼠毒性實(shí)驗(yàn)

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

表2 白鼠毒性實(shí)驗(yàn)數(shù)據(jù)

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

表3 似然比檢驗(yàn)與Wald檢驗(yàn)

圖2 白鼠毒性實(shí)驗(yàn)中系數(shù)值的相對(duì)輪廓似然值

討 論

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

7．江一濤，胡海蘭，魏巧玲，等．競爭風(fēng)險(xiǎn)模型的發(fā)展與應(yīng)用．中國衛(wèi)生統(tǒng)計(jì)，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．