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

蘭州大學公共衛(wèi)生學院流行病與衛(wèi)生統(tǒng)計學研究所(730000)

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

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

蘭州大學公共衛(wèi)生學院流行病與衛(wèi)生統(tǒng)計學研究所(730000)

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

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

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

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

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

(1)

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

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

(2)

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

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

實例分析

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

第一步：生存時間的設定

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

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

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

定義指示變量dNij為第i個觀察對象在第j個區(qū)間的“生存”情況，則當dNij=0時，個體i在第j個區(qū)間的數(shù)據(jù)是刪失的；否則，dNij=1，并取dNij服從參數(shù)為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個時間區(qū)間的觀測情況，通過step(e)這個函數(shù)來表示，如果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]在時間區(qū)間[t,t+dt)上的增量，和O[i,j]意義一樣，當dNij=0時，個體i在第j個區(qū)間的數(shù)據(jù)是刪失的；否則，dNij=1。

# Y[i]表示第i個研究對象的生存結局，取值為0或1。如果Y[i]=1，表示出現(xiàn)生存結局；否則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為基線風險函數(shù)在第j個時間區(qū)間上的增量，則有 h0[t]=sum(dh0[1:j])，j=1，2，…，T。

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

}

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

mu[j] <- dh0.e[j] * c # dh0.e[j]是未知風險函數(shù)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表示生存時間分段區(qū)間的尺寸，即Δt=t[j+1] - t[j]

#其他先驗信息

for (k in 1:6)

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

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

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 自相關函數(shù)圖

圖4 迭代診斷圖

圖5 迭代軌跡圖

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

討 論

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

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

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

[2] 韋來生編著.貝葉斯統(tǒng)計.北京:高等教育出版社,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軟件介紹及應用.中國衛(wèi)生統(tǒng)計,2017,34(1):170-172.

[7] 孫振球,徐勇勇主編.醫(yī)學統(tǒng)計學.第4版.北京:人民衛(wèi)生出版社,2014,291-294.

[8] 方積乾主編.衛(wèi)生統(tǒng)計學.第7版.北京:人民衛(wèi)生出版社,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] 張堯庭,陳漢峰.貝葉斯統(tǒng)計推斷.北京:科學出版社,1991.

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

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

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

張 悅)