Cox比例風(fēng)險(xiǎn)回歸模型的貝葉斯估計(jì)方法研究*

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

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

Cox比例風(fēng)險(xiǎn)回歸模型的貝葉斯估計(jì)方法研究*

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

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

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

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

貝葉斯Cox比例風(fēng)險(xiǎn)回歸模型原理

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

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

(1)

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

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

(2)

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

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

實(shí)例分析

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

第一步：生存時(shí)間的設(shè)定

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

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

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

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

model{

#指示變量的設(shè)定

for(i in 1:N) {

for(j in 1:T) {

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

# O[i,j]為研究對(duì)象i在第j個(gè)時(shí)間區(qū)間的觀測(cè)情況，通過(guò)step(e)這個(gè)函數(shù)來(lái)表示，如果e>=0，O[i,j]=step(e)=1，表示研究對(duì)象能被觀測(cè)到；否則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]在時(shí)間區(qū)間[t,t+dt)上的增量，和O[i,j]意義一樣，當(dāng)dNij=0時(shí)，個(gè)體i在第j個(gè)區(qū)間的數(shù)據(jù)是刪失的；否則，dNij=1。

# Y[i]表示第i個(gè)研究對(duì)象的生存結(jié)局，取值為0或1。如果Y[i]=1，表示出現(xiàn)生存結(jié)局；否則Y[i]=0，表示刪失。

# Doodle模型的構(gòu)建

for(j in 1:T) {

for(i in 1:N) {

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

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

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

#dh0[j]=h0[j]dt為基線風(fēng)險(xiǎn)函數(shù)在第j個(gè)時(shí)間區(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]是未知風(fēng)險(xiǎn)函數(shù)dh0[j]的一個(gè)先驗(yàn)猜測(cè)，c為該先驗(yàn)猜測(cè)的精度。

c <- 0.001

r <- 0.1

for (j in 1:T) {

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

#在該例中，我們?cè)O(shè)dh0.e[j] <- r *Δt，其中r是對(duì)研究對(duì)象在每個(gè)單位時(shí)間內(nèi)刪失率的猜測(cè)，Δt表示生存時(shí)間分段區(qū)間的尺寸，即Δt=t[j+1] - t[j]

#其他先驗(yàn)信息

for (k in 1:6)

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

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

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))

通過(guò)Gibbs抽樣和馬爾科夫迭代，最終得到的結(jié)果見(jiàn)表1和圖1～5。

表1 基于OpenBUGS的貝葉斯估計(jì)結(jié)果

圖1 迭代歷史圖

圖2 核密度圖

圖3 自相關(guān)函數(shù)圖

圖4 迭代診斷圖

圖5 迭代軌跡圖

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

討 論

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

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

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

[2] 韋來(lái)生編著.貝葉斯統(tǒng)計(jì).北京:高等教育出版社,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的貝葉斯生存分析理論及其在可靠性評(píng)估中的應(yīng)用.南京:南京理工大學(xué),2008.

[6] 張繼巍,高文龍,李娟生，等.OpenBUGS軟件介紹及應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2017,34(1):170-172.

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

[8] 方積乾主編.衛(wèi)生統(tǒng)計(jì)學(xué).第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)計(jì)推斷.北京:科學(xué)出版社,1991.

教育部人文社科項(xiàng)目(15XJC910001)；蘭州大學(xué)高校基本科研業(yè)務(wù)費(fèi)(LZUjbky-2016-025)

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

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

張 悅)