病例隊列設計下加性風險回歸分析及其在乳腺癌數據中的應用

楊 杰,丁潔麗

(武漢大學數學與統計學院,湖北武漢,430072)

1 引言

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

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

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

本文的主要結構為:第2節介紹病例隊列設計下加性風險模型參數的加權估計方法并綜述其漸近理論,第3節為模擬計算,第4節為乳腺癌數據的生存分析.

2 加性風險模型與病例隊列抽樣

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

我們考慮如下加性風險模型:給定協變量Zi(t)時,失效時間的風險函數有如下形式:

和

其中a?2=aa0,以及

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

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

其中

上述加權估計方程有如下顯式解:

和

其中

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

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

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

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

(C4)存在某個α∈(0,1),使得當N→∞時,=n0/N→α.

(C5)矩陣

是非奇異矩陣,其中

進一步,有

其中

這里

以及

其中

易得

是A?1{Σ1+Σ2}A?1的一個相合估計.

3 模擬研究

本節我們通過一系列模擬研究來展示上節中討論的加權估計方法在有限樣本量下的優良表現,展示病例隊列設計下加性風險回歸方法的應用價值.考慮如下加性風險模型:

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

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

在每種參數設定下,比較上述四種方法的參數估計值的均值(Mean),估計值的樣本標準差(SD),標準差估計值的均值(SE),95%的正態區間覆蓋率(CP)以及估計的相對效率(RE),結果均基于1000次的模擬結果計算獲得.模擬結果請見表1和表2.

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

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

4 乳腺癌數據的生存分析

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

我們感興趣的因變量是乳腺癌患者的生存時間,而觀測到的生存時間存在刪失,其刪失率為84.4%.我們考慮如下6個潛在影響因素.患者年齡(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).癌細胞分化程度(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),腫瘤無法評估(T=5).是否并發淋巴癌的N分期分為5種類型:同側腋窩無腫大淋巴結(N=1),同側腋窩有尚可推動的腫大淋巴結(N=2),同側腋窩腫大淋巴結融合或粘連(N=3),有同側胸骨旁淋巴結轉移(N=4),區域淋巴結無法評估(N=5).腫瘤是否轉移的M分期(M)分為2種類型:腫瘤未轉移(M=1),腫瘤轉移(M=2).我們對數據中考慮的上述影響因素進行了描述性統計分析,畫出了其條形圖(圖1)及生存曲線圖(圖2).

圖1 影響因素條形圖

圖2 生存曲線條形圖

為了探究乳腺癌患者生存時間的影響因素.我們建立如下加性風險模型:

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

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

表3:乳腺癌細胞瘤研究數據分析結果