Cox比例風險回歸模型的貝葉斯估計方法研究*

張繼巍 高文龍 李學朝 拉扎提·木拉提 秦天燕 李娟生

蘭州大學公共衛生學院流行病與衛生統計學研究所(730000)

Cox比例風險回歸模型的貝葉斯估計方法研究*

張繼巍?高文龍?李學朝 拉扎提·木拉提 秦天燕 李娟生△

蘭州大學公共衛生學院流行病與衛生統計學研究所(730000)

Cox比例風險回歸模型是目前進行多因素生存分析最常用的半參數模型，由于其兼有參數模型和非參數模型的優點，并可以在數據不完全的情況下分析研究對象生存時間的影響因素，因而得到了廣泛的應用[1]。而貝葉斯統計分析方法可以有效利用先驗信息，在小樣本數據推斷中具有明顯優勢[2]，其在生存分析中也得到了廣泛應用，比如Jennifer和Mike實現了貝葉斯在Weibull模型中的應用[3]，Lee等人實現了貝葉斯方法在生存數據指數分布中的應用[4]。因此將兩者結合對Cox回歸模型進行貝葉斯估計可以有效的處理小樣本、數據不完全和運行環境復雜等問題[5]，在生存分析中具有廣闊的發展前景。

本文對生存分析模型中的Cox比例風險回歸模型基于貝葉斯條件下的無信息先驗分布，對小樣本、隨機截尾和刪失的不完全數據，結合馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法中Gibbs抽樣思想，利用OpenBUGS軟件[6]進行估計，研究該模型在OpenBUGS環境下的建模過程，最后結合生存數據對該模型進行應用研究。

貝葉斯Cox比例風險回歸模型原理

貝葉斯Cox比例風險回歸模型是在原Cox回歸模型[7]的基礎上，利用貝葉斯統計分析方法的基本原理對其回歸參數進行估計。該方法將模型中各估計參數看成隨機變量，并對其指定適當的先驗分布，再結合樣本信息，通過Gibbs抽樣MCMC模擬計算其后驗分布，最終實現模型參數的估計。我們假設有i個觀測對象和p個潛在影響因素，則Cox回歸模型的基本形式為[8]：

h(t,Xi)=h0(t)exp(β1X1+β2X2+…βpXp)

(1)

式中h(t,X)是具有協變量X的個體在時刻t時的風險函數；第i個對象的生存時間為ti；協變量X=(X1，X2，…,Xp)是影響生存時間的p個有關因素，這些變量可以是定量的，也可以是定性的，在整個觀測期間內不受時間的影響；h0(t)是所有協變量取值為0時的風險函數，稱為基線風險函數；βi=(β1,β2,…,βp)為Cox模型的回歸系數，是一組待估計的回歸參數；在貝葉斯分析中，需要先對這些參數指定適當的先驗分布，通常不合適的先驗會對結果產生影響。因此，按照Kalbfleisch and Prentice (1980)，Clayton (1991)，Clayton (1994)等人[9]的建議將Cox回歸模型系數βi設定無信息正態先驗分布。

βi(i=1,2,…,p)～N(μ,τ)

(2)

式(2)中τ=1/σ2，考慮到該先驗分布對參數后驗分布的影響，因此設定μ=0，精度τ=0.000001。生存結局Y為0-1指示性變量(如果觀測對象的生存時間刪失，則Y=0，否則Y=1)。

在OpenBUGS中對該模型進行參數估計時，為了得到更加可靠的后驗參數估計，需要設定不同初始值的多條鏈進行模擬迭代，并且先要確定模型的退火參數burn-in的值，以保證模型在收斂狀態下對后驗參數進行抽樣估計；為了降低各條初始值鏈之間的自相關性和蒙特卡洛誤差，需要設置參數thin值大于1，并需要更多的抽樣樣本用于參數估計；模型的收斂情況可以通過觀測遍歷均值、迭代軌跡圖和BGR圖直觀的判斷，當它們都趨于穩定時，可以認為模型收斂[9]。

實例分析

Cox比例風險回歸模型主要用于影響因素分析、校正協變量后的組間比較和多變量生成預測[8]。本研究選取文獻[7]中介紹的例19-5，利用收集的63例某種惡性腫瘤患者的生存數據，通過構建Cox比例風險模型，結合貝葉斯統計分析方法試進行其生存情況的影響因素分析。變量的賦值和數據見原文獻[7]。利用貝葉斯統計方法進行Cox回歸模型參數估計的步驟如下：

第一步：生存時間的設定

基于貝葉斯條件下的Cox回歸模型建立時，我們通常需要把生存時間ti分成J個等距離的時間區間0=a0

本例中，將生存時間ti∈[2,120]分成了16個等距離的區間，即T=16。

第二步：Cox回歸模型在貝葉斯條件下的建立：

定義指示變量dNij為第i個觀察對象在第j個區間的“生存”情況，則當dNij=0時，個體i在第j個區間的數據是刪失的；否則，dNij=1，并取dNij服從參數為Idtij的泊松分布，此時我們建立的Cox模型為：

model{

#指示變量的設定

for(i in 1:N) {

for(j in 1:T) {

O[i,j] <- step(obs.t[i] - t[j] + eps)

# O[i,j]為研究對象i在第j個時間區間的觀測情況，通過step(e)這個函數來表示，如果e>=0，O[i,j]=step(e)=1，表示研究對象能被觀測到；否則O[i,j]=0，表示已刪失。

dN[i,j] <- O[i,j] * step(t[j + 1] - obs.t[i] - eps) * Y[i] } }

# dN[i,j]為O[i,j]在時間區間[t,t+dt)上的增量，和O[i,j]意義一樣，當dNij=0時，個體i在第j個區間的數據是刪失的；否則，dNij=1。

# Y[i]表示第i個研究對象的生存結局，取值為0或1。如果Y[i]=1，表示出現生存結局；否則Y[i]=0，表示刪失。

# Doodle模型的構建

for(j in 1:T) {

for(i in 1:N) {

dN[i,j] ～ dpois(Idt[i,j])

#假設dN[i,j]服從均值為Idt[i,j]的泊松分布

Idt[i,j] <- O[i,j] * exp(β′X) * dh0[j]

#dh0[j]=h0[j]dt為基線風險函數在第j個時間區間上的增量，則有 h0[t]=sum(dh0[1:j])，j=1，2，…，T。

S[i,j] <- pow(exp(-sum(dh0[1:j])),exp(β′X)) #Cox回歸模型的生存率函數，其中函數pow(e1,e2)=e1^e2

}

dh0[j] ～ dgamma(mu[j],c) }

mu[j] <- dh0.e[j] * c # dh0.e[j]是未知風險函數dh0[j]的一個先驗猜測，c為該先驗猜測的精度。

c <- 0.001

r <- 0.1

for (j in 1:T) {

dh0.e[j] <- r * (t[j+1] - t[j]) }

#在該例中，我們設dh0.e[j] <- r *Δt，其中r是對研究對象在每個單位時間內刪失率的猜測，Δt表示生存時間分段區間的尺寸，即Δt=t[j+1] - t[j]

#其他先驗信息

for (k in 1:6)

Beta[k] ～ dnorm(0.0,0.000001) }

第三步：利用OpenBUGS軟件對模型進行貝葉斯估計：通過步驟二構建的模型，結合貝葉斯統計分析工具OpenBUGS軟件[6]，對模型中的未知參數進行貝葉斯估計。本例的數據結構如下：

list( N=63,T=16,eps=1.0E-10,

t=c(2,3,4,5,7,12,15,16,18,19,24,26,29,35,40,66,120),

obs.t=c(52,51,35,103,7,60,58,29,70,67,66,87,85,82,76,74,63,101,100,66,93,24,93,90,15,3,87,120,120,120,120,120,120,40,26,120,120,120,3,120,7,18,120,120,15,4,120,16,24,19,120,24,2,120,12,5,120,120,7,40,108,24,16),

X1=c(54,57,58,43,48,40,44,36,39,42,42,42,51,55,49,52,48,54,38,40,38,19,67,37,43,49,50,53,32,46,43,44,62,40,50,33,57,48,28,54,35,47,49,43,48,44,60,40,32,44,48,72,42,63,55,39,44,42,74,61,45,38,62),

X2=c(0,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1,1,1,0,1,0,0,1,0,1,0,1,1,1,0,1,1,0,1,1,1,1,1,0,1,0,0,1,0,1,0,1,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0),

X3=c(0,1,0,1,1,1,1,0,1,1,1,0,1,1,1,1,1,0,1,1,0,0,0,0,0,0,1,1,1,1,0,0,0,1,0,1,1,0,0,0,1,0,0,1,1,0,1,0,1,0,0,1,0,0,1,0,0,1,0,1,0,1,0),

X4=c(1,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1,1,1,0,1,0,0,1,1,0,0,1,1,1,0,1,1,0,1,0,0,1,0,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,1,1,0,0,1,0,0,1,0,0),

X5=c(1,0,1,1,0,0,0,1,0,0,0,1,0,0,0,0,0,1,0,0,1,1,1,1,1,1,1,0,0,0,1,1,1,0,1,0,0,1,1,1,1,1,1,0,0,1,0,1,0,1,1,1,1,1,0,1,1,0,1,1,1,0,1),

X6=c(0,0,1,0,0,0,0,1,1,1,0,0,0,1,1,1,0,1,0,1,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,0,0,0),

Y=c(0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,1,1,0,0,0,1,1,1,1,0,0,1,1,0,1,1,1,0,1,1,0,1,1,0,0,1,1,0,1,1))

通過Gibbs抽樣和馬爾科夫迭代，最終得到的結果見表1和圖1～5。

表1 基于OpenBUGS的貝葉斯估計結果

圖1 迭代歷史圖

圖2 核密度圖

圖3 自相關函數圖

圖4 迭代診斷圖

圖5 迭代軌跡圖

本例采用初始值不同的兩條鏈模擬迭代，拋去前10000次抽樣以保證樣本是從后驗分布中抽取，額外的20000次抽樣用于參數估計；由圖2～5可以看出各鏈混合較好，模型收斂；因此我們可以得到該模型中的參數估計值(表1)，這與我們從SAS和SPSS[7]得到的結果基本一致。

討 論

本文利用Cox回歸模型在生存分析中的優勢，結合貝葉斯統計分析方法，構建基于貝葉斯條件下的Cox回歸模型，通過其統計分析工具OpenBUGS軟件[6]進行Gibbs抽樣和馬爾科夫迭代，實現對模型中參數的貝葉斯估計。一方面，克服了小樣本資料在Cox回歸模型中的限制[7]；另一方面，提高了在數據刪失、不完全的狀態下對研究對象生存影響因素估計的精度[5]。

貝葉斯統計分析方法是近幾十年來發展起來的一種數理統計方法。隨著其統計分析工具OpenBUGS軟件[6]的不斷更新和完善，貝葉斯統計分析方法日益受到重視并在各個領域得到廣泛的應用[11-13]。基于貝葉斯統計框架下的Cox比例風險回歸模型將待估計的參數作為隨機變量而不是常數，并對其指定適當的先驗信息，再結合該參數的樣本信息和總體信息，得到參數的后驗信息進而實現參數的貝葉斯估計[2]。但是，值得注意的是在利用貝葉斯統計分析方法進行推斷時最重要的，也是最受經典統計學派批判的是先驗信息的選擇以及如何利用先驗信息確定先驗分布[2]。由于先驗信息常常會受到主觀因素的影響，盡管許多統計學家提出了多種方法，但至今理論上仍然沒有一個統一的確定先驗信息的方法。而本文對Cox回歸模型中待估計參數指定的無信息正態先驗分布，已經得到了多次驗證[9,14]。除此之外，本文所采用的實例滿足比例風險假定的要求[7]，因此可以利用其進行實證分析。

[1] Cox DR,Hinkly DV.Theoretical Statistics.New York:John Wiley &Sons.1974.

[2] 韋來生編著.貝葉斯統計.北京:高等教育出版社,2016，3.

[3] Jennifer Clarke,Mike West.Bayesian Weibull tree models for survival analysis of clinico-genomic data.Statistical Methodology,2008,5(3):238-262.

[4] Jaeyong Lee,Jinseog Kim,Sin-Ho Jung.Bayesian analysis of paired survival data using a bivariate exponential distribution.Lifetime Data Anal,2007,13(1):119-137.

[5] 林靜.基于MCMC的貝葉斯生存分析理論及其在可靠性評估中的應用.南京:南京理工大學,2008.

[6] 張繼巍,高文龍,李娟生，等.OpenBUGS軟件介紹及應用.中國衛生統計,2017,34(1):170-172.

[7] 孫振球,徐勇勇主編.醫學統計學.第4版.北京:人民衛生出版社,2014,291-294.

[8] 方積乾主編.衛生統計學.第7版.北京:人民衛生出版社,2012,420-422.

[9] Spiegelhalter D,Thomas A,Best N,et al.Open BUGS user manual (version 3.2.3).Cambridge:MRC Biostatistics Unit,2014.

[10] Congdon P.Applied Bayesian Modelling.England:John Wiley and Sons,2003.

[11] Lyle W,Konigsberg,Frankenberg.Bayes in Biological Anthropology.American Journal of Physical Anthropology,2013,152(57):153-184.

[12] Du Changde,Luo Ali,Yang Haifeng.Adaptive stellar spectral subclass classification based on Bayesian SVMs.New Astronomy,2017,51:51-58.

[13] Spence GT,Steinsaltz D,Fanshawe TR.A Bayesian approach to sequential meta-analysis.Statistic in Medicine,2016,35(29):5356-5375.

[14] 張堯庭,陳漢峰.貝葉斯統計推斷.北京:科學出版社,1991.

教育部人文社科項目(15XJC910001)；蘭州大學高?；究蒲袠I務費(LZUjbky-2016-025)

?張繼巍，高文龍為共同第一作者

△通信作者：李娟生，E-mail：lijsh@lzu.edu.cn

張 悅)