參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型

胡軍浩，黃熒，楊青

（中南民族大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)學(xué)院，武漢 430074）

競(jìng)爭(zhēng)風(fēng)險(xiǎn)是指一個(gè)研究對(duì)象會(huì)發(fā)生多個(gè)結(jié)局事件，并且一個(gè)結(jié)局事件的發(fā)生會(huì)影響其他結(jié)局事件發(fā)生的概率，這些結(jié)局事件間便存在競(jìng)爭(zhēng)風(fēng)險(xiǎn).競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型作為一種處理具有多種潛在結(jié)局生存數(shù)據(jù)的分析方法，在臨床醫(yī)學(xué)、流行病、經(jīng)濟(jì)學(xué)等領(lǐng)域的研究中越來(lái)越重要.1992年，WEI將線性回歸模型的建模方法引入生存分析領(lǐng)域，提出了相較于經(jīng)典的COX比例風(fēng)險(xiǎn)模型更加簡(jiǎn)單、直觀且易于理解的AFT模型［1］.1999年，F(xiàn)INE和GRAY基于累積發(fā)生函數(shù)和特定失敗類型的邊際風(fēng)險(xiǎn)概率，提出了一種能直接解釋特定終點(diǎn)生存概率的子分布比例風(fēng)險(xiǎn)模型［2］，其準(zhǔn)確性和可靠性引發(fā)了廣泛的研究興趣.2009年，PENG和FINE將分位數(shù)回歸引入到競(jìng)爭(zhēng)風(fēng)險(xiǎn)設(shè)置中，并建立了一個(gè)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型［3］，使得協(xié)變量效應(yīng)能夠得到合理解釋.該模型擴(kuò)展了AFT模型，為子分布比例風(fēng)險(xiǎn)模型提供了一個(gè)有力補(bǔ)充，同時(shí)也為競(jìng)爭(zhēng)風(fēng)險(xiǎn)的研究提供了一個(gè)新的視角.

傳統(tǒng)的競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型采用最小化L1型凸函數(shù)的方法來(lái)估計(jì)其回歸系數(shù)，然而對(duì)不平滑估計(jì)函數(shù)的最小值問(wèn)題求解，可能會(huì)導(dǎo)致解的不唯一性.同時(shí)，一次估計(jì)一個(gè)條件分位數(shù)，不僅估計(jì)效率低，而且估計(jì)結(jié)果難以解釋.文獻(xiàn)［4］通過(guò)計(jì)算標(biāo)準(zhǔn)高斯隨機(jī)擾動(dòng)上的期望值，得到不平滑估計(jì)函數(shù)的平滑近似，文獻(xiàn)［5］利用核光滑方法，對(duì)競(jìng)爭(zhēng)風(fēng)險(xiǎn)下的累積發(fā)生函數(shù)光滑化，進(jìn)行光滑估計(jì)方程的非參數(shù)估計(jì).受文獻(xiàn)［6］的啟發(fā)，本文將系數(shù)函數(shù)參數(shù)化建模的思想，應(yīng)用到競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型中，建立了一個(gè)參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型，并分析了其估計(jì)量的漸近性質(zhì)，給出了模型選擇過(guò)程，最后，通過(guò)仿真模擬和實(shí)際數(shù)據(jù)對(duì)模型進(jìn)行評(píng)價(jià)和實(shí)現(xiàn).

1 競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型

設(shè)T為感興趣的生存時(shí)間，可能受到刪失C的影響，∈={1，…，K}表示失敗的原因，Z是第一維恒為常數(shù)1的p+1維協(xié)變量向量.觀察到n個(gè)獨(dú)立同分 布 的 樣 本{(Xi，δi，Zi)，i=1，…，n}，其 中X=min(T，C)，δ=I(T≤C)∈.存在一個(gè)與T相關(guān)的潛在時(shí)間變量，原因k下的生存時(shí)間T*k=I(∈=k)×T+I(∈≠k)×∞.基于累積發(fā)生函數(shù)，定義Tk*的條件分位數(shù)：

為便于后續(xù)討論，其他事件存在的情況下，僅對(duì)原因1事件的發(fā)生感興趣.對(duì)于τ∈[τL，τU]，0＜τL，τU＜1，廣義線性競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型為：

其中β(τ)為p+1維未知回歸系數(shù)向量，g(·)是一個(gè)已知的單調(diào)連接函數(shù).

當(dāng)競(jìng)爭(zhēng)風(fēng)險(xiǎn)數(shù)據(jù)不存在刪失時(shí)，即X=T，δ=∈，可直接通過(guò)最小絕對(duì)偏差方法來(lái)估計(jì)系數(shù).β(τ)最小化目標(biāo)函數(shù)為：

也可以表示為相應(yīng)梯度函數(shù)的零點(diǎn)：

當(dāng)競(jìng)爭(zhēng)風(fēng)險(xiǎn)數(shù)據(jù)存在刪失C時(shí)，假設(shè)在給定協(xié)變量Z下，刪失時(shí)間C條件獨(dú)立于生存時(shí)間T.文獻(xiàn)［3］采用刪失加權(quán)的逆概率方法，得到調(diào)整后的加權(quán)估計(jì)方程：

2 參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型

考慮為模型（1）中的系數(shù)函數(shù)β(τ)引入一個(gè)有限維參數(shù)θ，定義：

其中b(τ)=[1，b1(τ)，…，br(τ)]T為τ的r+1維已知函數(shù)向量，θ是元素為θjh，j=0，…，p，h=0，…，r的(p+1)×(r+1)維矩陣，指定j=0，h=0分別表示Z和b(τ)的常數(shù)項(xiàng).第j個(gè)協(xié)變量Zj對(duì)應(yīng)的回歸系數(shù)為βj(τ|θ)=θj0+b1(τ)θj1+…+br(τ)θjr，j=1，2，…，p.

由此定義參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型如下：

模型（6）可以通過(guò)控制θ的一些元素為0，來(lái)實(shí)現(xiàn)對(duì)部分系數(shù)函數(shù)的參數(shù)化，或參數(shù)化為b(τ)的子集函數(shù).例如：假定系數(shù)β1與b(τ)的常數(shù)項(xiàng)無(wú)關(guān)，β2不參數(shù)化，指定系數(shù)函數(shù)由二次多項(xiàng)式建模，b(τ)=[1，τ，τ2]T，因此設(shè)置θ3×3，其中θ10=θ20=θ22=0.估計(jì)的分位數(shù)函數(shù)Q1(τ|Z，θ^)關(guān)于τ的一階導(dǎo)數(shù)小于零，則會(huì)發(fā)生分位數(shù)交叉，文獻(xiàn)［7］利用一種強(qiáng)制單調(diào)性的約束優(yōu)化算法解決了分位數(shù)交叉問(wèn)題，并在R軟件qrcm包中得以實(shí)現(xiàn).在確保模型不發(fā)生分位數(shù)交叉的前提下，b(τ)可以靈活選擇，文獻(xiàn)［8］中給出了更多參數(shù)函數(shù)的例子.

由于競(jìng)爭(zhēng)和刪失的存在，大的分位數(shù)可能無(wú)法識(shí)別.本文在分位數(shù)區(qū)間[τL，τU]上，將文獻(xiàn)［6］提出的積分損失最小化（ILM）方法，應(yīng)用于模型（6）的參數(shù)估計(jì)過(guò)程.估計(jì)為如下積分梯度函數(shù)的零點(diǎn)：

其中τi=F1(Xi|Zi，θ)為θ下原因1的累積發(fā)生函數(shù).類似地，結(jié)合式（3）與式（7），得到刪失數(shù)據(jù)下的積分梯度函數(shù)：

式（8）和式（9）均為θ的平滑估計(jì)函數(shù)，可通過(guò)Newton Raphson方法或文獻(xiàn)［9］提出的梯度搜索法進(jìn)一步求解，有效避免了不連續(xù)函數(shù)的系數(shù)估計(jì)過(guò)程.與文獻(xiàn)［3］一次估計(jì)一個(gè)分位數(shù)對(duì)應(yīng)的回歸系數(shù)不同，本文通過(guò)引入?yún)?shù)，找到了分位數(shù)與回歸系數(shù)之間連接的橋梁，允許一次估計(jì)整個(gè)分位數(shù)過(guò)程，在進(jìn)行平滑估計(jì)的同時(shí)，極大地提高了估計(jì)效率.同樣的想法也被延伸到刪失和截?cái)鄶?shù)據(jù)的分位數(shù)回歸［10］、縱向數(shù)據(jù)的分位數(shù)回歸［11］.進(jìn)一步地，文獻(xiàn)［12］在分位數(shù)框架下，討論了協(xié)變量和參數(shù)規(guī)格均為非線性的情況，是對(duì)系數(shù)函數(shù)參數(shù)化建模方法的一個(gè)很好的發(fā)展.

3 漸近性分析

引理1［13］在獨(dú)立同分布的數(shù)據(jù)下，Θ為緊集，a(Zi，θ)在每個(gè)θ∈Θ上連續(xù)的概率為1，且存在d(Z)，對(duì)所有θ∈Θ，都有‖ ‖a(Z，θ)≤d(Z)，E[d(Z)]＜∞，則可得到E[a(Z，θ)]是連續(xù)的，且

定理1若滿足如下條件：

（i）Θ為緊集；

（ii）存 在 一 個(gè) 正 數(shù)v，有P(C=v)＞0，P(C＞v)=0；

有界；

（iv）Z，B(τ)，Bˉ(τ)均為一致有界；

證明條件（ii）和條件（iii）為獨(dú)立刪失數(shù)據(jù)的分位數(shù)回歸中的標(biāo)準(zhǔn)假設(shè)，且根據(jù)拉格朗日中值定理易證得為θ的連續(xù)函數(shù)，從而在每個(gè)θ∈Θ上連續(xù)的概率為1，這也意味著b(τ)為定義良好的連續(xù)函數(shù).在條件（iv）下，可知積分梯度函數(shù)一致有界，即：

結(jié)合條件（i）和引理1，則可知Sˉ0(θ)連續(xù)，且：

再根據(jù)式（10）和式（11），可得：

對(duì)任意ε＞0，由式（12）可得：

因此結(jié)合式（13）～（15），對(duì)任意ε＞0，有：

設(shè)N為任意包含θ0的Θ的開子集，則Θ∩Nc為緊集.結(jié)合條件（v），有：

即Q0(θ)在θ0處取得唯一最大值，對(duì)于θ~∈Θ∩Nc，則有：

定理2在滿足定理1的條件下，若滿足如下條件：

（i）θ0為Θ的內(nèi)點(diǎn)；

（v）Z和b(τ)均 為 非 奇 異 的.對(duì) 于G=則有漸近正態(tài)性：

證明類似于文獻(xiàn)［13］的證明思想，在大樣本中估計(jì)量近似等于樣本均值的線性組合，由中心極限定理可得到漸近正態(tài)性.根據(jù)條件（i）和（ii），一階條 件 滿 足的 概 率 接 近1，其 中為Hessian矩陣.在一階條件下，利用拉格朗日中值定理，在θ0周圍擴(kuò)展并乘以得到：

條件（v）意味著G為滿秩矩陣，結(jié)合逆矩陣的連續(xù)性，由式（20）得到：

通過(guò)Slutsky定理，證得估計(jì)量的漸近正態(tài)性：

4 擬合優(yōu)度測(cè)量

文獻(xiàn)［14］和文獻(xiàn)［15］均指定累積發(fā)生函數(shù)為特定的分布，提出了競(jìng)爭(zhēng)風(fēng)險(xiǎn)場(chǎng)景下的擬合優(yōu)度程序，然而關(guān)于參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型的擬合優(yōu)度測(cè)量目前沒(méi)有一個(gè)正式的數(shù)值檢驗(yàn)方法.評(píng)估式（4）的最小值意味著對(duì)數(shù)據(jù)整體擬合優(yōu)度情況的度量，式（4）的值越小，模型與數(shù)據(jù)的擬合情況越好［16］.式（4）的最小值類似于似然函數(shù)的最大似然值，據(jù)此，本文結(jié)合式（4）與信息準(zhǔn)則提出參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型的擬合優(yōu)度評(píng)估過(guò)程.使用文獻(xiàn)［17］提出的AIC信息準(zhǔn)則的調(diào)整版本：

和文獻(xiàn)［18］提出的BIC信息準(zhǔn)則的調(diào)整版本：

5 數(shù)值模擬

5.1 仿真數(shù)據(jù)生成

設(shè)置協(xié)變量向量Z=(1，Z1，Z2)T，其中Z1由標(biāo)準(zhǔn)均勻分布隨機(jī)生成，Z2由概率為0.5的伯努利分布隨機(jī)生成.兩個(gè)相互競(jìng)爭(zhēng)的失敗原因?yàn)椤?1和∈=2，原因1發(fā)生的概率滿足P(∈=1|Z)=p0I(Z2=0)+p1I(Z2=1).生存時(shí)間服從條件分布P(T≤t|∈=1，Z)=Φ(logt-γTZ)，P(T≤t|∈=2，Z)=Φ(logtαTZ).刪失時(shí)間C服從(0，cu)上的均勻分布，通過(guò)點(diǎn)cu來(lái)控制刪失率.在上述設(shè)置下，對(duì)于τ∈[0.1，0.4]，潛在的競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型可以描述為：

分別在兩個(gè)樣本容量n=100，n=200下，設(shè)置兩個(gè)不同的模擬場(chǎng)景.場(chǎng)景1：cu=5，p0=0.8，p1=0.6，a0=(1，0，-0.5)T，γ0=(0，0，-0.5)T，該 場(chǎng) 景下大約有20%的個(gè)體發(fā)生刪失，55%的個(gè)體失敗于原因1，25%的個(gè)體失敗于原因2.場(chǎng)景2：cu=7.6，p0=0.8，p1=0.6，a0=(0，0.5，-0.5)T，γ0=(0，-1，-0.5)T，該場(chǎng)景下大約有10%的個(gè)體發(fā)生刪失，60%的個(gè)體失敗于原因1，30%的個(gè)體失敗于原因2.

5.2 參數(shù)估計(jì)量的有限樣本性質(zhì)

為評(píng)估估計(jì)量的有限樣本性能，通過(guò)800次蒙特卡羅模擬，比較參數(shù)估計(jì)量β^n(τ)=β(τ|θ^n)與競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型中估計(jì)量β^n(τ)的偏差和標(biāo)準(zhǔn)誤，兩個(gè)估計(jì)量的偏差均很小，不作報(bào)告，標(biāo)準(zhǔn)誤的測(cè)量結(jié)果如表1所示.

表1 模擬場(chǎng)景下競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型（CRQR）和參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型（PCRQR）的標(biāo)準(zhǔn)誤對(duì)比Tab.1 Comparison of standard errors of competing risks quantile regression model（CRQR）and parametric competing risks quantile regression model（PCRQR）in simulated scenarios

仿真實(shí)驗(yàn)結(jié)果表明，參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型在不同場(chǎng)景、不同樣本量下，估計(jì)的標(biāo)準(zhǔn)誤都要低于傳統(tǒng)的競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型，在較大的分位數(shù)下優(yōu)勢(shì)更明顯.與期望的結(jié)果相同，參數(shù)結(jié)構(gòu)下的競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型有更好的性能.

6 真實(shí)數(shù)據(jù)分析

將參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型應(yīng)用到文獻(xiàn)［19］有關(guān)濾泡細(xì)胞淋巴瘤疾病的研究數(shù)據(jù)中.該數(shù)據(jù)集由541名疾病早期的濾泡淋巴瘤（I或II）患者組成，接受單純放療（化療=0）或放療和化療的聯(lián)合治療（化療=1）.疾病復(fù)發(fā)或?qū)χ委煙o(wú)反應(yīng)和緩解期死亡為兩個(gè)結(jié)局事件，實(shí)驗(yàn)對(duì)疾病復(fù)發(fā)或?qū)χ委煙o(wú)反應(yīng)（RNT）這個(gè)結(jié)局事件感興趣，緩解期死亡事件視為它的競(jìng)爭(zhēng)事件.患者的年齡（age：平均=57），血紅蛋白水平（hgb：平均=138）也被記錄.其中有272個(gè)感興趣的事件、76個(gè)競(jìng)爭(zhēng)風(fēng)險(xiǎn)事件和193個(gè)刪失事件.首先對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行預(yù)處理，將協(xié)變量分別表示為Z1：年齡/10，Z2：血紅蛋白/100，Z3：臨床階段，Z4：治療方法.用如下參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型描述RNT事件的生存時(shí)間：

為了使模型更好地?cái)M合數(shù)據(jù)，根據(jù)第4節(jié)提出的擬合優(yōu)度測(cè)量過(guò)程，首先對(duì)b(τ)的規(guī)格進(jìn)行選擇，不同規(guī)格下模型的擬合優(yōu)度評(píng)估結(jié)果如表2所示.

通過(guò)系數(shù)函數(shù)圖像進(jìn)行模型的初步篩選后，表2列出了7個(gè)不同規(guī)格指定下，模型所含的參數(shù)個(gè)數(shù)及AIC*值、BIC*值.模型1和2指定b(τ)為τ的冪函數(shù)形式.模型3指定為τ的指數(shù)函數(shù)形式，該模型涉及參數(shù)較少，但從AIC*，BIC*的值來(lái)看，模型與數(shù)據(jù)的擬合效果不理想.模型4和模型5分別使用了非對(duì)稱邏輯分布和三角分布的分位數(shù)函數(shù)，相比前面的模型，它們擁有較少的參數(shù)數(shù)量，且擬合效果也是比較令人滿意的.基于上述模型的表現(xiàn)，冪律分布、邏輯分布和三角分布應(yīng)該被考慮，據(jù)此指定了模型6.模型7是對(duì)模型6的進(jìn)一步調(diào)整，雖然涉及參數(shù)個(gè)數(shù)較多，但其AIC*和BIC*的值相較于其他模型來(lái)說(shuō)更小，擬合情況是讓人滿意的.選擇模型7為最終模型，其參數(shù)的估計(jì)結(jié)果及估計(jì)標(biāo)準(zhǔn)誤如表3所示.

表2 擬合優(yōu)度測(cè)量Tab.2 Goodness-of-fit measures

表3 模型7的參數(shù)估計(jì)結(jié)果Tab.3 Parameter estimation results for model 7

表3的最后一列為協(xié)變量零假設(shè)檢驗(yàn)的顯著性，類似的，最后一行為b(τ)分量零假設(shè)的P值，*號(hào)表示顯著性的程度（括號(hào)中為標(biāo)準(zhǔn)誤）.由表3的結(jié)果得到系數(shù)函數(shù)的表達(dá)式，例如截距項(xiàng)對(duì)應(yīng)的估計(jì)系數(shù)函數(shù)為：

協(xié)變量血紅蛋白水平對(duì)應(yīng)的估計(jì)系數(shù)為：

圖1 模型7中截距項(xiàng)和協(xié)變量Z2對(duì)應(yīng)的回歸系數(shù)函數(shù)Fig.1 Regression coefficients functions corresponding to the intercept and covariate Z2 in model 7

7 結(jié)論

參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型是對(duì)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型［3］的發(fā)展與改進(jìn)，不僅解決了非平滑系數(shù)估計(jì)函數(shù)產(chǎn)生的多重根問(wèn)題，而且在參數(shù)結(jié)構(gòu)下有著簡(jiǎn)單、高效和易于解釋的優(yōu)點(diǎn).本文構(gòu)建了參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型，給出了積分損失最小化方法下的參數(shù)求解過(guò)程，其估計(jì)量滿足一致性和漸近正態(tài)性.通過(guò)實(shí)驗(yàn)?zāi)M數(shù)據(jù)集評(píng)估估計(jì)量的有限樣本性能，結(jié)果表明：參數(shù)競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)模型有優(yōu)于競(jìng)爭(zhēng)風(fēng)險(xiǎn)分位數(shù)回歸模型的估計(jì)性能.