王淑影,張亞男,程云飛

(長春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,吉林 長春 130012)

一、引言

在經(jīng)濟(jì)學(xué)、生物醫(yī)學(xué)、臨床試驗(yàn)等研究領(lǐng)域中,由于各種原因經(jīng)常會(huì)出現(xiàn)數(shù)據(jù)刪失的情況。例如,在臨床試驗(yàn)過程中,由于患者的異質(zhì)性或者費(fèi)用的限制等原因?qū)е略谠囼?yàn)結(jié)束時(shí),某些個(gè)體還是沒有發(fā)生感興趣的事件。此時(shí)研究人員不能精確觀測到某個(gè)個(gè)體的確切失效(死亡)時(shí)間,只能觀測到生存時(shí)間大于刪失時(shí)間,即右刪失數(shù)據(jù)。針對(duì)右刪失數(shù)據(jù)構(gòu)建模型,最常用的就是比例風(fēng)險(xiǎn)(Proportional Hazards,PH)模型和加速失效時(shí)間(Accelerated Failure Time,AFT)模型。兩者的不同之處在于PH模型假設(shè)對(duì)數(shù)危險(xiǎn)率函數(shù)與協(xié)變量之間存在線性關(guān)系,而AFT模型直接對(duì)生存時(shí)間構(gòu)建模型,研究協(xié)變量與對(duì)數(shù)生存時(shí)間的線性關(guān)系。但是,在不同實(shí)際問題中,PH模型的比例風(fēng)險(xiǎn)假設(shè)可能很難驗(yàn)證,因此本文直接考慮更具有解釋性的AFT模型來對(duì)右刪失數(shù)據(jù)進(jìn)行建模。

在傳統(tǒng)的AFT模型構(gòu)建過程中,往往都是在假設(shè)協(xié)變量給定的情況下進(jìn)行統(tǒng)計(jì)推斷研究,忽略了此過程中的不確定性,導(dǎo)致估計(jì)量的方差被低估,進(jìn)而得到不準(zhǔn)確的預(yù)測結(jié)果。在這種情況下,有必要將傳統(tǒng)的AFT模型與模型平均方法結(jié)合起來,避免忽略模型選擇過程引入的不確定性。因此,本文考慮在右刪失數(shù)據(jù)下,構(gòu)建基于信息準(zhǔn)則的AFT模型平均方法,并且基于估計(jì)和預(yù)測結(jié)合的雙重視角來分析模型平均后AFT模型的擬合效果,以此來避免選擇單一模型潛在的風(fēng)險(xiǎn),減少估計(jì)的均方誤差,提高覆蓋概率,為以后更好地對(duì)AFT模型進(jìn)行研究提供新的思路。

二、文獻(xiàn)述評(píng)

加速失效時(shí)間(AFT)模型最早是由Pieruschka在1961年提出的,后來國內(nèi)外專家學(xué)者對(duì)加速失效時(shí)間模型進(jìn)行了大量的研究[1]。例如Buckley等將加速失效時(shí)間模型應(yīng)用到了右刪失數(shù)據(jù)下,并提出了經(jīng)典的Buckley-James估計(jì)方法[2]。Tian等在現(xiàn)狀和區(qū)間刪失數(shù)據(jù)下,對(duì)加速失效時(shí)間模型提出了通過Wald-type來構(gòu)造的一種新的參數(shù)估計(jì)方法[3]。Arno?t等針對(duì)任意刪失數(shù)據(jù),提出了一種半?yún)?shù)方法來對(duì)加速失效時(shí)間模型進(jìn)行估計(jì)[4]。Li等考慮在帶治愈組的右刪失數(shù)據(jù)下,提出了一種半?yún)?shù)加速失效時(shí)間治愈模型,并針對(duì)該模型開發(fā)了一種EM算法[5]。Scolas等針對(duì)帶治愈組的區(qū)間刪失數(shù)據(jù),提出了一個(gè)十分靈活的擴(kuò)展廣義Gamma加速失效時(shí)間模型[6]?，F(xiàn)有文獻(xiàn)大多考慮具有固定協(xié)變量的加速失效時(shí)間模型,但是這樣做可能會(huì)丟失候選模型中包含的一些有用信息。克服這種問題的有效途徑就是使用模型平均方法,通過對(duì)候選模型的估計(jì)量進(jìn)行加權(quán),可以有效地減少有用信息的遺失。

模型平均方法分為貝葉斯模型平均(Bayesian Model Averaging,BMA)和頻率模型平均(Frequentist Model Averaging,FMA)兩大類。許多早期工作從貝葉斯的角度討論了模型平均方法,比如高磊等對(duì)BMA方法在一些難點(diǎn)問題上的理論進(jìn)展進(jìn)行了綜述,并介紹了BMA方法在大數(shù)據(jù)背景下的應(yīng)用前景[7]。Draper采用BMA方法進(jìn)行預(yù)測研究并通過油價(jià)預(yù)測等兩個(gè)實(shí)例證明了該方法的有效性[8]。但是,考慮到BMA方法存在模型假設(shè)十分復(fù)雜且難以直觀解釋的問題,FMA方法開始受到越來越多的關(guān)注和研究。例如,Zhang等提出了一個(gè)基于函數(shù)線性回歸模型的交叉驗(yàn)證模型平均估計(jì)方法來預(yù)測響應(yīng)變量[9]。Bai等針對(duì)期望回歸模型進(jìn)行了最優(yōu)模型平均估計(jì),并證明了所得模型平均估計(jì)量是漸近最優(yōu)的[10]。Zhang等在固定參數(shù)設(shè)置的條件下,研究了Mallows模型平均估計(jì)和Jackknife模型平均估計(jì)下模型權(quán)重的漸近性質(zhì)等[11]。這些工作為頻率模型平均研究提供了深刻的理論支持和有效的方法。

文獻(xiàn)梳理表明,國內(nèi)外關(guān)于模型平均方法的研究主要聚焦于完整數(shù)據(jù)類型,對(duì)復(fù)雜刪失數(shù)據(jù)的研究不多,且很少有文章研究加速失效時(shí)間模型的模型平均問題。因此無論在候選模型集合構(gòu)建以及如何評(píng)價(jià)最終模型的好壞等方面都有待進(jìn)一步深入研究。

三、數(shù)據(jù)、模型及似然函數(shù)

分布函數(shù)為:

密度函數(shù)為:

則對(duì)數(shù)似然函數(shù)為:

四、基于信息準(zhǔn)則的AFT模型平均方法

在生物統(tǒng)計(jì)領(lǐng)域中,已有研究大多是先選定一個(gè)模型,然后對(duì)選定模型進(jìn)行統(tǒng)計(jì)推斷,這樣可能會(huì)造成估計(jì)結(jié)果不穩(wěn)健。而模型平均是解決這個(gè)問題的有效方法,模型平均方法通過考慮模型的不確定性和優(yōu)化模型權(quán)重來最小化誤差。

在本文中,我們考慮p+q個(gè)協(xié)變量,其中xi中包括p個(gè)確定被包含在模型中的協(xié)變量,zi中包括q個(gè)可能被包含在模型中的協(xié)變量,那么就會(huì)產(chǎn)生2q個(gè)候選模型。對(duì)第m個(gè)候選模型,假設(shè)為如下AFT模型:

均方根誤差是一個(gè)能夠較好地衡量估計(jì)精度的評(píng)價(jià)值,一般來說均方根誤差值越小則估計(jì)精度越高[12]。

在實(shí)際情況中,受試者可能更關(guān)心的指標(biāo)是生存概率,也就是在單位時(shí)間內(nèi)生存的可能性大小。

對(duì)應(yīng)累積生存率為0.5時(shí)的生存時(shí)間叫做中位生存時(shí)間。這是一個(gè)非常重要的指標(biāo),有時(shí)比平均生存時(shí)間重要[13]。

其中,r為T分布的50%分位點(diǎn)。

如果給定了一組數(shù)據(jù),但是產(chǎn)生數(shù)據(jù)的真實(shí)模型未知時(shí),需要找到所有可能的模型,這時(shí)就要通過比較研究哪個(gè)模型的擬合效果較好來決定,這個(gè)過程就是模型選擇。在模型選擇中,最常用的是基于Akaike信息準(zhǔn)則(AIC)的模型選擇和基于Bayesian信息準(zhǔn)則(BIC)的模型選擇,表達(dá)式分別如下所示:

AICm=-2logm+2lm,BICm=-2logm+lmlogn

通過模型選擇可能得到存在有局部誤差的模型,進(jìn)而導(dǎo)致估計(jì)結(jié)果不準(zhǔn)確,這時(shí)模型平均方法就會(huì)被考慮。模型平均方法最核心的問題就在于權(quán)重的選擇,本文采用基于Smoothed AIC(S-AIC)和Smoothed BIC(S-BIC)的模型平均方法來計(jì)算組合權(quán)重,權(quán)重的表達(dá)式如下:

其中,m代表第m個(gè)模型,ωm是各個(gè)候選模型的權(quán)重,xIC表示AIC或BIC。設(shè)ω=(ω1,ω2,…,ωm)T是m個(gè)模型的權(quán)重向量,則:

對(duì)每個(gè)候選模型進(jìn)行極大似然估計(jì)后會(huì)得到各個(gè)參數(shù)對(duì)應(yīng)的估計(jì)值,則模型平均后的未知參數(shù)估計(jì)的表達(dá)形式如下:

模型平均后的生存概率的具體形式如下:

模型平均后的中位生存時(shí)間表達(dá)形式如下:

其中,r為T分布的50%分位點(diǎn)。

五、模擬研究

其中,i=1,2,…,n,β1、β2、γ1、γ2、γ3、γ4是回歸系數(shù),σ是尺度參數(shù),εi是服從標(biāo)準(zhǔn)正態(tài)分布的誤差項(xiàng),μ是截距項(xiàng)。xi1、xi2是確定包含在模型中的協(xié)變量,zi1、zi2、zi3、zi4是可能包含在模型中的協(xié)變量,那么候選模型為24=16個(gè)。各個(gè)協(xié)變量由以下分布生成:xi1～N(0,1),xi2～N(0,1),zi1～N(0.5,1),zi2～U(0,1),zi3～P(1),zi4～B(1,0.5)。C服從均勻分布U(0,cl),其中cl控制刪失率。本文考慮兩組參數(shù)真值,設(shè)置A為窄模型(即只包含確定協(xié)變量的模型):β1=1,β2=1,γ1=0,γ2=0,γ3=0,γ4=0,μ=0.5,σ=0.4;設(shè)置B為寬模型(即包含確定協(xié)變量和所有可能協(xié)變量的模型):β1=-0.2,β2=0.2,γ1=0.2,γ2=-0.3,γ3=0.2,γ4=-0.2,μ=0.5,σ=0.4。刪失比例為30%,模擬循環(huán)1 000次。

表1 設(shè)置A兩組樣本量下三個(gè)感興趣指標(biāo)的均方根誤差

表2 設(shè)置B兩組樣本量下三個(gè)感興趣指標(biāo)的均方根誤差

在表1和表2中,AIC代表基于AIC準(zhǔn)則的模型選擇方法,BIC代表基于BIC的模型選擇方法,S-AIC代表基于Smoothed AIC準(zhǔn)則的模型平均方法,S-BIC代表基于Smoothed BIC準(zhǔn)則的模型平均方法(下同)。在設(shè)置A中,窄模型是真實(shí)模型。從理論上來說,BIC是一個(gè)一致的模型選擇器,通常適用于具有少量參數(shù)的模型,因此在真實(shí)模型為窄模型的情況下,BIC優(yōu)于AIC,S-BIC優(yōu)于S-AIC。在指標(biāo)(a)和指標(biāo)(c)下,基于模型平均方法得到的結(jié)果與基于模型選擇方法得到的結(jié)果相差較大,基于信息準(zhǔn)則的AFT模型平均方法在參數(shù)數(shù)量較少時(shí)的估計(jì)精度更高。對(duì)于設(shè)置B,如果真實(shí)模型中的非零參數(shù)較多,S-BIC還是最好的方法。這可能是因?yàn)锳IC假設(shè)真實(shí)模型是高維的,需要在有限的情況下盡可能多的參數(shù)來描述它。Sakamoto等指出,AIC不是估計(jì)真實(shí)模型的標(biāo)準(zhǔn),而是最佳模型擬合的標(biāo)準(zhǔn)[15]。也就是說,AIC力求找到最佳的近似模型,如果數(shù)據(jù)較少模型的維數(shù)將會(huì)降低,隨著可用信息的增多,模型的維數(shù)將會(huì)增加。相比之下,BIC是一致的。它假設(shè)真實(shí)模型存在且是低維的,可對(duì)真實(shí)模型進(jìn)行一致性估計(jì)。

總的來說,無論真實(shí)模型為窄模型還是寬模型,基于Smoothed BIC準(zhǔn)則的AFT模型平均方法的表現(xiàn)都是最好的。這主要是因?yàn)槟Ｐ瓦x擇方法在使用過程中是將最終選擇的模型視為提前選擇,忽略了模型選擇過程帶來的額外不確定性,進(jìn)而導(dǎo)致了較大的估計(jì)誤差。隨著樣本量大小的增加,在局部錯(cuò)誤指定的假設(shè)下,參數(shù)將會(huì)更接近真實(shí)值,所有候選模型都更接近真實(shí)模型,因此模型平均方法和模型選擇方法的估計(jì)精度都在逐漸提高。

六、肺癌數(shù)據(jù)分析

使用本文提出的模型平均方法對(duì)一組肺癌數(shù)據(jù)進(jìn)行研究,以此來進(jìn)一步證明所提出方法的有效性。這組完整的數(shù)據(jù)集來自于R軟件包survival中,是為比較接受兩種不同治療方案的肺癌患者生存時(shí)間而設(shè)計(jì)的一個(gè)臨床試驗(yàn)。在這項(xiàng)研究中,137名肺癌患者被隨機(jī)分為兩組,分別接受兩種不同的治療方案,感興趣的是肺癌患者的生存時(shí)間。數(shù)據(jù)集中有6個(gè)協(xié)變量,包括治療方案(x1)、細(xì)胞類型(x2)、卡諾斯基評(píng)分(z1)、從診斷到治療的時(shí)間(z2)、年齡(z3)以及是否有過既往治療(z4)。我們感興趣的問題是基于上述6個(gè)協(xié)變量,對(duì)數(shù)據(jù)進(jìn)行預(yù)測。Kalbfleisch等人的研究結(jié)果表明,不同的治療方案和細(xì)胞類型會(huì)對(duì)肺癌患者的生存時(shí)間產(chǎn)生顯著影響[16]。所以,選定治療方案和細(xì)胞類型為必選協(xié)變量,其他4個(gè)協(xié)變量則為候選協(xié)變量,這樣共有24=16個(gè)不同的候選模型。

為了比較AIC、BIC、S-AIC和S-BIC四種方法的效果,隨機(jī)地將數(shù)據(jù)分為兩部分,分別記為訓(xùn)練集D(0)和測試集D(1)[17]。訓(xùn)練集的樣本量n0依次選取了80,90,100,110和120,而測試集的樣本量大小為n-n0,其中n=137是整體的樣本量。首先對(duì)n0個(gè)樣本進(jìn)行極大似然估計(jì),其次對(duì)n-n0個(gè)樣本進(jìn)行預(yù)測,并計(jì)算樣本的均方預(yù)測誤差(Mean Squared Prediction Error,MSPE),即:

從表3中MSPE的均值和中位數(shù)的角度來看,最優(yōu)的方法是S-AIC,它比其他三種方法的MSPE的均值和中位數(shù)都要小。S-AIC和S-BIC方法的MSPE的均值和中位數(shù)都比相應(yīng)的AIC和BIC方法要小,這說明基于模型平均方法的AFT模型預(yù)測結(jié)果比模型選擇的預(yù)測結(jié)果更準(zhǔn)確。此外S-AIC方法略優(yōu)于S-BIC方法,AIC方法略優(yōu)于BIC方法。盡管表3中四個(gè)方法的MSPE的均值和中位數(shù)有較大的不同,但是可以清楚地分辨出各個(gè)方法的表現(xiàn)。為了檢驗(yàn)所提方法得出的結(jié)論,進(jìn)一步考慮了Diebold-Mariano檢驗(yàn),通過R軟件包forecast中的dm.test函數(shù)計(jì)算出兩種方法構(gòu)造的模型之間的DM統(tǒng)計(jì)量和p值,以此來觀察兩種模型預(yù)測能力的優(yōu)劣性[18]。DM檢驗(yàn)的原假設(shè)為H0:DM統(tǒng)計(jì)量均值為0,即兩種方法構(gòu)造的模型預(yù)測能力相同;備擇假設(shè)為H1:DM統(tǒng)計(jì)量的均值不為0,也就是說兩種方法構(gòu)造的模型具有不同的預(yù)測能力。正的DM值表明分母對(duì)應(yīng)的方法優(yōu)于分子對(duì)應(yīng)的方法,若DM檢驗(yàn)的p值小于0.05,拒絕原假設(shè),那么就說明兩種方法構(gòu)造的模型具備不同的預(yù)測能力。

表3 肺癌患者生存時(shí)間的MSPE

表4 肺癌患者生存時(shí)間MSPE的DM統(tǒng)計(jì)量

在表4的前3列中,我們可以看到DM的值都是負(fù)的,且p值幾乎為0,這表明S-AIC方法顯著優(yōu)于其他三種方法。第4列和第6列中的DM值都是正的,且p值幾乎為0,說明S-BIC方法顯著優(yōu)于AIC和BIC方法。第5列中的DM值都是負(fù)的,且p值幾乎為0,說明AIC方法優(yōu)于BIC方法。這些結(jié)論與表3中的結(jié)論是一致的。

七、總結(jié)及討論

本文主要考慮右刪失數(shù)據(jù)下,提出基于AFT模型的頻率模型平均方法,使用極大似然方法進(jìn)行參數(shù)估計(jì)及預(yù)測,其中誤差項(xiàng)服從標(biāo)準(zhǔn)正態(tài)分布。模擬研究表明,所提方法的表現(xiàn)更好。模型平均方法考慮了所有協(xié)變量對(duì)響應(yīng)變量的影響,而模型選擇是只用選擇出來的部分協(xié)變量來解釋響應(yīng)變量。由于模型平均方法是根據(jù)協(xié)變量的個(gè)數(shù)來構(gòu)建候選模型,所以當(dāng)協(xié)變量比較多的時(shí)候,具體實(shí)現(xiàn)計(jì)算就會(huì)變得十分復(fù)雜[19]。近些年來頻率模型平均方法快速發(fā)展,提出了各種最優(yōu)模型平均的方法,但本文僅考慮了基于信息準(zhǔn)則的模型平均方法。如果使用漸近最優(yōu)的模型平均方法,如Mallows模型平均方法、Jackknife模型平均方法、OPT模型平均方法等,也許會(huì)得到更好的估計(jì)結(jié)果。這些問題值得我們進(jìn)一步研究。