韩 栋 陈 征△ 陈平雁 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．