楊 杰,丁潔麗

(武漢大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,湖北武漢,430072)

1 引言

生存數(shù)據(jù)應(yīng)用在許多學(xué)科領(lǐng)域,如生物醫(yī)學(xué)、流行病學(xué)、公共衛(wèi)生學(xué)、可靠性工程學(xué)、經(jīng)濟(jì)學(xué)和保險(xiǎn)精算學(xué)等領(lǐng)域.生存分析主要研究某特定事件的發(fā)生時(shí)間(例如:死亡時(shí)間、疾病發(fā)生時(shí)間、系統(tǒng)失效時(shí)間、索賠時(shí)間等等)與重要影響因素和協(xié)變量之間的關(guān)聯(lián).對(duì)于生存數(shù)據(jù)的研究,比例風(fēng)險(xiǎn)模型是應(yīng)用最為廣泛的統(tǒng)計(jì)半?yún)?shù)模型之一(Cox,1972[1];Andersen&Gill,1982[2];Lin,1994[3];Chen&Little,1999[4];Kalb fleisch&Prentice,2002[5]等等).作為比例風(fēng)險(xiǎn)模型的一個(gè)重要替代,加性風(fēng)險(xiǎn)模型假設(shè)基準(zhǔn)風(fēng)險(xiǎn)函數(shù)與協(xié)變量效應(yīng)之間具有一個(gè)加性結(jié)構(gòu).在很多實(shí)際應(yīng)用中,加性風(fēng)險(xiǎn)模型往往能更好地?cái)M合數(shù)據(jù)(Lin&Ying,1994[3];Mckeague&Sasieni,1994[6]).而當(dāng)加性風(fēng)險(xiǎn)模型和比例風(fēng)險(xiǎn)模型均能較好地?cái)M合數(shù)據(jù)時(shí),加性風(fēng)險(xiǎn)模型的回歸參數(shù)更容易解釋其實(shí)際意義(Lin&Ying,1994[3];Zeng&Cai,2010[7]).

在許多大型隊(duì)列研究中,往往涉及對(duì)大量研究個(gè)體的長(zhǎng)期追溯和隨訪.當(dāng)重要影響因素或協(xié)變量的采集非常昂貴時(shí),采用傳統(tǒng)的簡(jiǎn)單隨機(jī)抽樣可能會(huì)導(dǎo)致實(shí)驗(yàn)過于昂貴而超過預(yù)算.因此,發(fā)展和研究能節(jié)約成本和提高效率的抽樣機(jī)制具有非常重要的意義.對(duì)于帶有刪失的生存數(shù)據(jù),病例隊(duì)列設(shè)計(jì)(Case-Cohort design)是應(yīng)用最為廣泛的有偏抽樣機(jī)制之一.其主要機(jī)制是:首先,從全隊(duì)列中隨機(jī)地抽取一個(gè)子隊(duì)列(subcohort),全隊(duì)列中所有發(fā)生了感興趣事件的個(gè)體稱為病例(case).然后,子隊(duì)列和子隊(duì)列之外所有的病例組成病例隊(duì)列樣本.最后,僅對(duì)病例隊(duì)列樣本中的個(gè)體進(jìn)行昂貴協(xié)變量的采集和觀測(cè).對(duì)病例隊(duì)列設(shè)計(jì)相關(guān)統(tǒng)計(jì)方法的研究已有大量和廣泛的工作,比例風(fēng)險(xiǎn)模型(Prentice,1986[8];Self&Prentice,1988[9];Chen&Lo,1999[10];Lin&Ying,1993[11];Kulich&Lin,2004[12]),加性風(fēng)險(xiǎn)模型(Kulich&Lin,2000[13];Sun et al,2004[14]),加速失效模型(Kong et al,2004[15];Lu&Tsiatis,2006[16]).近來,基于多元失效時(shí)間的病例隊(duì)列設(shè)計(jì)的研究越來越廣泛(Kong&Cai,2009[17];Kang et al,2013[18];Kim et al,2018[19];Maitra et al,2020[20]).

本文主要探討病例隊(duì)列設(shè)計(jì)中加性風(fēng)險(xiǎn)模型下參數(shù)的統(tǒng)計(jì)推斷方法和應(yīng)用.首先,我們探討如何應(yīng)用加性風(fēng)險(xiǎn)模型擬合由病例隊(duì)列設(shè)計(jì)獲取的生存數(shù)據(jù),考慮參數(shù)的一種加權(quán)估計(jì)方法并綜述其漸近理論.然后,重點(diǎn)研究這種病例隊(duì)列設(shè)計(jì)下的分析方法在實(shí)際中的應(yīng)用問題.一方面,我們編寫了可實(shí)現(xiàn)這種統(tǒng)計(jì)分析方法的R軟件應(yīng)用程序.通過模擬研究展示了這種方法在有限樣本量下的優(yōu)良表現(xiàn),并評(píng)估了病例隊(duì)列設(shè)計(jì)相較于簡(jiǎn)單隨機(jī)抽樣設(shè)計(jì)的有效性.另一方面,我們應(yīng)用該方法分析了一個(gè)實(shí)際的乳腺癌數(shù)據(jù),展示了其成本效益與應(yīng)用價(jià)值.乳腺癌是女性最常見的惡性腫瘤之一,全世界每年約有120萬女性患上乳腺癌,約有50萬女性患者死亡.中國(guó)癌癥年發(fā)病數(shù)為160萬,現(xiàn)癥病人200多萬.乳腺癌占女性惡性腫瘤發(fā)病率第一位,每年約有4萬女性死于乳腺癌[21].本文將基于來源于美國(guó)國(guó)家癌癥研究所(National Cancer Institute)[22]的乳腺癌數(shù)據(jù),探索影響乳腺癌患者生存時(shí)間的影響因素.首先,我們基于全隊(duì)列數(shù)據(jù)應(yīng)用加性風(fēng)險(xiǎn)回歸方法分析數(shù)據(jù),探索出了一些對(duì)乳腺癌患者生存時(shí)間有顯著性影響的因素.進(jìn)一步,應(yīng)用病例隊(duì)列設(shè)計(jì)抽取樣本,并基于病例隊(duì)列樣本進(jìn)行加性風(fēng)險(xiǎn)回歸分析.結(jié)果表明,病例隊(duì)列設(shè)計(jì)僅用了較小的樣本量就達(dá)到了與全隊(duì)列研究幾乎一致的結(jié)果.當(dāng)協(xié)變量的測(cè)量非常昂貴時(shí),病例隊(duì)列設(shè)計(jì)可提高研究效率和節(jié)約成本.

本文的主要結(jié)構(gòu)為:第2節(jié)介紹病例隊(duì)列設(shè)計(jì)下加性風(fēng)險(xiǎn)模型參數(shù)的加權(quán)估計(jì)方法并綜述其漸近理論,第3節(jié)為模擬計(jì)算,第4節(jié)為乳腺癌數(shù)據(jù)的生存分析.

2 加性風(fēng)險(xiǎn)模型與病例隊(duì)列抽樣

假設(shè)一個(gè)大型隊(duì)列研究包含N個(gè)獨(dú)立的研究個(gè)體.我們用表示第i個(gè)個(gè)體的潛在失效時(shí)間,Ci表示第i個(gè)個(gè)體的刪失時(shí)間或隨訪時(shí)間(i=1,···,N).記觀測(cè)時(shí)間為Ti=min(,Ci),右刪失示性變量為?i=I(≤Ci).用Yi(t)=I(Ti≥t)表示風(fēng)險(xiǎn)過程,Ni(t)=?iI(Ti≤t)表示計(jì)數(shù)過程,其中I(·)表示示性函數(shù).記Zi(t)為第i個(gè)個(gè)體在t時(shí)刻處的p維協(xié)變量.記τ為事件終止時(shí)間.

我們考慮如下加性風(fēng)險(xiǎn)模型:給定協(xié)變量Zi(t)時(shí),失效時(shí)間的風(fēng)險(xiǎn)函數(shù)有如下形式:

和

其中a?2=aa0,以及

在研究的病例隊(duì)列設(shè)計(jì)下,首先,通過簡(jiǎn)單隨機(jī)抽樣方式從全隊(duì)列中抽取一個(gè)樣本容量為n0的子隊(duì)列.子隊(duì)列中的個(gè)體和子隊(duì)列之外的所有病例個(gè)體組成病例隊(duì)列樣本,記其樣本量為n.然后,僅對(duì)病例隊(duì)列樣本中的個(gè)體觀測(cè)其協(xié)變量.具體來說,記ξi為子隊(duì)列示性變量,即:ξi=1表示第i個(gè)個(gè)體被選入子隊(duì)列,ξi=0表示第i個(gè)個(gè)體未被選入子隊(duì)列.假設(shè)每個(gè)個(gè)體入選子隊(duì)列的概率P(ξi=1)==n0/N.因此,病例隊(duì)列設(shè)計(jì)下,當(dāng)ξi=1或?i=1時(shí),觀測(cè)數(shù)據(jù)為(Ti,?i,Zi);當(dāng)ξi=0且?i=0時(shí),觀測(cè)數(shù)據(jù)為(Ti,?i).

由于病例隊(duì)列設(shè)計(jì)中,協(xié)變量不是對(duì)每一個(gè)個(gè)體都進(jìn)行了觀測(cè),因此(2.2)中的估計(jì)方程方法不再適用,需要提出新的推斷方法.受Horvitz&Thompson(1951)[23]逆概率加權(quán)思想和Liang&Ziger(1986)[24]廣義估計(jì)方程思想的啟發(fā).在病例隊(duì)列的設(shè)計(jì)下,對(duì)加性風(fēng)險(xiǎn)模型(2.1)中的β和Λ0(t)的推斷可建立如下加權(quán)估計(jì)方程:

其中

上述加權(quán)估計(jì)方程有如下顯式解:

和

其中

這種逆概率加權(quán)的思想由Kalb fleisch&Lawless(1988)[25]首次應(yīng)用于生存分析數(shù)據(jù).基于多類型疾病的病例隊(duì)列設(shè)計(jì),Kang et al(2013)[18]在加性風(fēng)險(xiǎn)模型下為多元失效時(shí)間發(fā)展了逆概率加權(quán)推斷方法.以下我們討論和綜述了上述的漸近性質(zhì).設(shè)β0為β的真值.假設(shè)如下正則條件成立:

(C1)β的參數(shù)空間是緊集,Z的取值空間是緊集.

(C2)給定Zi時(shí),與Ci相互獨(dú)立,ξi與(Ti,?i,Zi)相互獨(dú)立.

(C3)P(Ti≥τ)>0 且 Λ0(t)<∞.

(C4)存在某個(gè)α∈(0,1),使得當(dāng)N→∞時(shí),=n0/N→α.

(C5)矩陣

是非奇異矩陣,其中

進(jìn)一步,有

其中

這里

以及

其中

易得

是A?1{Σ1+Σ2}A?1的一個(gè)相合估計(jì).

3 模擬研究

本節(jié)我們通過一系列模擬研究來展示上節(jié)中討論的加權(quán)估計(jì)方法在有限樣本量下的優(yōu)良表現(xiàn),展示病例隊(duì)列設(shè)計(jì)下加性風(fēng)險(xiǎn)回歸方法的應(yīng)用價(jià)值.考慮如下加性風(fēng)險(xiǎn)模型:

設(shè)定參數(shù)真值β1=0或0.5,β2=0或0.5.協(xié)變量Z1分別從均勻分布U(0,1)和正態(tài)分布N(0,1)中生成,Z2從成功率為0.5的Bernoulli分布中生成.設(shè)定基準(zhǔn)風(fēng)險(xiǎn)函數(shù)λ0(t)=1,則失效時(shí)間可以從風(fēng)險(xiǎn)率為λ0(t)+β1Z1+β2Z2的指數(shù)分布中生成.刪失時(shí)間C從均勻分布U[0,c]中生成,通過挑選c的不同取值從而產(chǎn)生相應(yīng)的刪失率,分別為ρ=80%,85%和90%.對(duì)于病例隊(duì)列設(shè)計(jì),設(shè)定全隊(duì)列樣本總量N=1000,子隊(duì)列樣本量為n0=200.

為了闡明問題,我們比較以下幾種方法:

在每種參數(shù)設(shè)定下,比較上述四種方法的參數(shù)估計(jì)值的均值(Mean),估計(jì)值的樣本標(biāo)準(zhǔn)差(SD),標(biāo)準(zhǔn)差估計(jì)值的均值(SE),95%的正態(tài)區(qū)間覆蓋率(CP)以及估計(jì)的相對(duì)效率(RE),結(jié)果均基于1000次的模擬結(jié)果計(jì)算獲得.模擬結(jié)果請(qǐng)見表1和表2.

表1:參數(shù)β1和β2的模擬結(jié)果,其中Z1～U(0,1).

表2:參數(shù)β1和β2的模擬結(jié)果,其中Z1～N(0,1).

4 乳腺癌數(shù)據(jù)的生存分析

本節(jié)我們研究乳腺癌相關(guān)數(shù)據(jù).數(shù)據(jù)來源于美國(guó)國(guó)家癌癥研究所(National Cancer Institute)[22].我們選取了1975-2017年期間40歲以上的女性乳腺癌患者的共118763條數(shù)據(jù).我們基于加性風(fēng)險(xiǎn)模型對(duì)此數(shù)據(jù)集進(jìn)行生存回歸分析,探索影響乳腺癌患者生存時(shí)間的主要影響因素.進(jìn)一步,我們應(yīng)用病例隊(duì)列設(shè)計(jì)分析數(shù)據(jù)集,展示此種有偏抽樣機(jī)制的實(shí)際應(yīng)用價(jià)值.

我們感興趣的因變量是乳腺癌患者的生存時(shí)間,而觀測(cè)到的生存時(shí)間存在刪失,其刪失率為84.4%.我們考慮如下6個(gè)潛在影響因素.患者年齡(Age)分為5組:40-49歲(Age=1),50-59歲(Age=2),60-69歲(Age=3),70-79歲(Age=4),80歲以上(Age=5).種族(Race)分為3種:其他種族(Race=1),白種人(Race=2),黑種人(Race=3).癌細(xì)胞分化程度(Grade)分為4種類型:腫瘤高度分化(Grade=1),腫瘤中度分化(Grade=2),腫瘤低分化(Grade=3),腫瘤未分化(Grade=4).腫瘤直徑大小的T分期(T)分為5種類型:腫瘤直徑≤2cm(T=1),2cm<腫瘤直徑≤5cm(T=2),腫瘤直徑>5cm(T=3),腫瘤直接侵犯胸壁或皮膚(T=4),腫瘤無法評(píng)估(T=5).是否并發(fā)淋巴癌的N分期分為5種類型:同側(cè)腋窩無腫大淋巴結(jié)(N=1),同側(cè)腋窩有尚可推動(dòng)的腫大淋巴結(jié)(N=2),同側(cè)腋窩腫大淋巴結(jié)融合或粘連(N=3),有同側(cè)胸骨旁淋巴結(jié)轉(zhuǎn)移(N=4),區(qū)域淋巴結(jié)無法評(píng)估(N=5).腫瘤是否轉(zhuǎn)移的M分期(M)分為2種類型:腫瘤未轉(zhuǎn)移(M=1),腫瘤轉(zhuǎn)移(M=2).我們對(duì)數(shù)據(jù)中考慮的上述影響因素進(jìn)行了描述性統(tǒng)計(jì)分析,畫出了其條形圖(圖1)及生存曲線圖(圖2).

圖1 影響因素條形圖

圖2 生存曲線條形圖

為了探究乳腺癌患者生存時(shí)間的影響因素.我們建立如下加性風(fēng)險(xiǎn)模型:

首先,基于全隊(duì)列118763條數(shù)據(jù)進(jìn)行加性風(fēng)險(xiǎn)回歸分析,結(jié)果總結(jié)在表3中(見Full-Cohort).結(jié)果表明,考慮的6個(gè)協(xié)變量對(duì)乳腺癌患者的生存時(shí)間均有顯著的影響.具體來說,年齡越大的患者死亡率越高.黑種人死亡風(fēng)險(xiǎn)率最高,白種人次之,其他人種最低.癌細(xì)胞分化程度越高,腫瘤大小的T分期等級(jí)越高,區(qū)域淋巴癌的N分期等級(jí)越高,患者的生存率越低.癌細(xì)胞有遠(yuǎn)處轉(zhuǎn)移的患者比無遠(yuǎn)處轉(zhuǎn)移的患者的死亡風(fēng)險(xiǎn)率高出1.476%.

為了評(píng)估病例隊(duì)列設(shè)計(jì)并展現(xiàn)其在實(shí)際應(yīng)用中的可行性與有效性,我們首先從全隊(duì)列中隨機(jī)抽取了一個(gè)容量為n0=30000的子隊(duì)列,然后將子隊(duì)列和子隊(duì)列之外的病例(nc=13903)組成病例隊(duì)列樣本,其樣本容量為43903.我們應(yīng)用本文研究的病例隊(duì)列設(shè)計(jì)下加性風(fēng)險(xiǎn)模型下的加權(quán)估計(jì)方法對(duì)數(shù)據(jù)進(jìn)行分析,其結(jié)果總結(jié)在表3中(見Case-Cohort).結(jié)果表明,病例隊(duì)列設(shè)計(jì)采用了較小的樣本量,但估計(jì)結(jié)果與基于全隊(duì)列數(shù)據(jù)分析得到的估計(jì)結(jié)果十分接近.病例隊(duì)列設(shè)計(jì)僅用了全隊(duì)列約37%的樣本量,對(duì)于Grade,T,N的估計(jì)分別達(dá)到了約52.6%,41.2%和45.8%的效率.這說明在實(shí)際應(yīng)用中,當(dāng)影響因素Grade,T和N的測(cè)量非常昂貴或耗時(shí)時(shí),采用病例隊(duì)列設(shè)計(jì)會(huì)提高研究效率和節(jié)約成本.

表3:乳腺癌細(xì)胞瘤研究數(shù)據(jù)分析結(jié)果