胡軍浩，黃熒，楊青

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

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

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

1 競爭風(fēng)險分位數(shù)回歸模型

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

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

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

當競爭風(fēng)險數(shù)據(jù)不存在刪失時，即X=T，δ=∈，可直接通過最小絕對偏差方法來估計系數(shù).β(τ)最小化目標函數(shù)為：

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

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

2 參數(shù)競爭風(fēng)險分位數(shù)模型

考慮為模型（1）中的系數(shù)函數(shù)β(τ)引入一個有限維參數(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ù)項.第j個協(xié)變量Zj對應(yīng)的回歸系數(shù)為βj(τ|θ)=θj0+b1(τ)θj1+…+br(τ)θjr，j=1，2，…，p.

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

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

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

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

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

3 漸近性分析

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

定理1若滿足如下條件：

（i）Θ為緊集；

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

有界；

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

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

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

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

對任意ε＞0，由式（12）可得：

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

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

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

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

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

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

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

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

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

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

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

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

5 數(shù)值模擬

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

7 結(jié)論

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