華東理工大學(xué)理學(xué)院(200237) 干晨晨 錢夕元

國(guó)內(nèi)對(duì)惡性腫瘤患者的生存分析研究大部分是以醫(yī)院數(shù)據(jù)為基礎(chǔ)的臨床生存分析。為了獲得惡性腫瘤患者出院后的生存情況，傳統(tǒng)的方式是對(duì)病人進(jìn)行隨訪。由于隨訪工作普遍存在花費(fèi)大，時(shí)間久，效率低等缺點(diǎn)，尋求更加簡(jiǎn)便高效的方法成為了一個(gè)亟待解決的難題。Fernando Colchero[1]等人開(kāi)發(fā)了一種專門用于生存分析的R軟件包BaSTA(Bayesian survival trajectory analysis)，即貝葉斯生存軌跡分析。該包采用捕獲-再捕獲方法，利用貝葉斯框架，對(duì)種群的年齡別死亡率和生存率進(jìn)行估計(jì)。本文嘗試著將BaSTA包在惡性腫瘤患者這一特殊種群中進(jìn)行應(yīng)用，借助已有的住院病人醫(yī)療記錄、檔案等對(duì)患者的年齡別死亡率和生存率進(jìn)行估計(jì)。

生存模型及BaSTA包簡(jiǎn)介

1.生存模型簡(jiǎn)介

在研究種群年齡別死亡率和生存率的捕獲-再捕獲實(shí)驗(yàn)中，每個(gè)個(gè)體i(i=1,2,…,n)均涉及到兩組數(shù)據(jù)。第一組是描述該個(gè)體縱向捕獲歷史的二進(jìn)制變量，即yi={yi,t}。若yi,t=1，則表示個(gè)體i在時(shí)間t被捕獲到了；若yi,t=0，則表示沒(méi)有被捕獲到。

第二組數(shù)據(jù)描述了該個(gè)體的出生時(shí)間(bi)與死亡時(shí)間(di)。其中，部分個(gè)體的出生時(shí)間或者死亡時(shí)間是缺失的，這可能是它在實(shí)驗(yàn)之前出生(左側(cè)截?cái)?或者在實(shí)驗(yàn)結(jié)束之后才死亡(右側(cè)刪失)，也可能是它在實(shí)驗(yàn)期間發(fā)生了死亡但未被捕獲到。

通常，對(duì)于捕獲-再捕獲數(shù)據(jù)集的年齡別生存分析，主要是通過(guò)結(jié)合生存模型和捕獲概率模型來(lái)進(jìn)行的。分別用x表示年齡，X表示死亡年齡，θ表示生存參數(shù)。假設(shè)年齡x連續(xù)，則年齡別死亡率為

(1)

因此，生存到年齡x的概率，即生存函數(shù)為

(2)

那么相對(duì)應(yīng)的在年齡x之前發(fā)生死亡的概率為

F(x|θ)=Pr(X

(3)

且死亡年齡的概率密度函數(shù)為

f(x|θ)=Pr(x≤X<(x+Δx))=S(x|θ)μ(x|θ)

(4)

在這里，我們采用貝葉斯方法分別對(duì)生存參數(shù)θ，未知死亡年齡Xu以及捕獲概率π進(jìn)行估計(jì)。

首先，生存參數(shù)的條件后驗(yàn)分布為

p(θ|Xu,Xk)∝p(Xu,Xk|θ)p(θ|θp)∝f(X|θ)p(θ|θp)

其中，p(θ|θp)是生存參數(shù)θ的先驗(yàn)分布。

其次，若個(gè)體i在第一次捕獲之前和最后一次捕獲之后共有ηi年未被捕獲到，則未知死亡年齡Xu的條件后驗(yàn)分布為

p(Xu|Xk,θ,yu,π)∝p(yu|π,Xu)p(Xk|Xu,θ)p(X|θp)∝(1-π)ηf(X|θ)p(X|θp)

最后，由于捕獲概率π的似然函數(shù)為二項(xiàng)分布密度的乘積，同時(shí)取其先驗(yàn)為具有超參數(shù)ρ1和ρ2的Beta先驗(yàn)，則捕獲概率π的共軛Beta條件后驗(yàn)為

在該貝葉斯分層框架下，只需利用Metropolis-within-Gibbs[2]采樣方法就可以獲得上述三者的估計(jì)值。

2.BaSTA包簡(jiǎn)介

BaSTA包采用貝葉斯方法，對(duì)一個(gè)部分?jǐn)?shù)據(jù)缺失的捕獲-再捕獲數(shù)據(jù)集進(jìn)行年齡別生存分析。它不但能夠?qū)ι鎱?shù)，捕獲概率，以及未知的出生和死亡時(shí)間進(jìn)行估計(jì)，還允許用戶測(cè)試一系列不同的生存模型。

在實(shí)際操作中，首先需要在R軟件中安裝BaSTA包，程序如下：

install.packages(“BaSTA”)；library(BaSTA)

實(shí)際操作BaSTA包主要從以下幾個(gè)步驟考慮:

(1) 數(shù)據(jù)讀取：讀取需要進(jìn)行年齡別生存分析的數(shù)據(jù)集

data <- read.table(“l(fā)ocation of the data set.txt”,sep=“ ”,header=TRUE)

讀取的數(shù)據(jù)集需具備表1的組織形式。

其中，每一行對(duì)應(yīng)于一個(gè)個(gè)體，且用ID作為唯一的身份識(shí)別標(biāo)志。對(duì)于每一個(gè)個(gè)體，若可以明確獲知其出生或者死亡的時(shí)間，就在出生時(shí)間列和死亡時(shí)間列記錄具體的年份，若無(wú)法獲知，就記為0。在捕獲-再捕獲實(shí)驗(yàn)中，假設(shè)捕獲周期的數(shù)量為T，時(shí)間間隔為1年，那么時(shí)間1，時(shí)間2，…，時(shí)間T分別表示實(shí)驗(yàn)開(kāi)始的第1年，第2年，…，第T年。若個(gè)體在該時(shí)間被捕獲到了就記為1，否則就記為0。這組長(zhǎng)度為T的數(shù)據(jù)構(gòu)成了個(gè)體的捕獲歷史。

表1 導(dǎo)入數(shù)據(jù)的格式

(2) 模型選擇：選擇用于數(shù)據(jù)集分析的生存模型

multiout<- multibasta(data,studyStart,studyEnd,

models=c(“EX”,“GO”,“WE”,“LO”),

shapes=c(“simple”,“Makeham”,“bathtub”))

studyStart和studyEnd分別表示捕獲-再捕獲實(shí)驗(yàn)開(kāi)始和結(jié)束的時(shí)間。

BaSTA包提供了四個(gè)基本的models(生存模型)，分別為：(a)Exponential，(b)Gompertz，(c)Weibull和(d)Logistic。每個(gè)模型所描述的年齡別死亡率和生存率具有不同的走向趨勢(shì)。此外，BaSTA包還允許用戶對(duì)這些基本模型進(jìn)行擴(kuò)展，從而具有更加復(fù)雜的模型形式。具體而言，就是定義了三種shapes形式：(a)simple，僅使用上述基本模型的死亡率函數(shù)；(b)Makeham，向函數(shù)中添加一個(gè)常數(shù)；(c)bathtub，向函數(shù)中添加一個(gè)常數(shù)和一個(gè)減小的Gompertz函數(shù)。其中，Exponential模型僅能與simple形式相組合。因此，BaSTA包總共包含了10種組合生存模型，具有相當(dāng)大的靈活性。選擇哪一種生存模型來(lái)進(jìn)行分析，可以利用偏差信息準(zhǔn)則(deviance information criterion，DIC)[3]和相關(guān)的統(tǒng)計(jì)量來(lái)對(duì)模型的擬合程度和復(fù)雜程度進(jìn)行評(píng)價(jià)。

multibasta()函數(shù)允許用戶在同一數(shù)據(jù)集上運(yùn)行具有不同函數(shù)形式的模型，進(jìn)行模型比較并可視化輸出運(yùn)行的結(jié)果。如果模型的所有參數(shù)都已經(jīng)收斂，那么BaSTA包會(huì)計(jì)算出DIC值。一般情況下，選擇DIC值較小的模型。

(3) 模型擬合：用選擇的生存模型對(duì)數(shù)據(jù)集進(jìn)行擬合，并且估計(jì)參數(shù)值

fit<- basta(data,studyStart,studyEnd,minAge,model,shape,niter,burnin,thinning,nsim,…)

其中，minAge表示年齡別生存分析的起始年齡，niter表示MCMC算法迭代的次數(shù)，burnin表示刪除的初始迭代的次數(shù)，thinning表示為了減小自相關(guān)所跳過(guò)的MCMC步驟數(shù)，以及nsim表示運(yùn)行的模擬鏈數(shù)目。

(4) 結(jié)果輸出：輸出參數(shù)估計(jì)值，并繪制相關(guān)曲線

summary(fit)

plot(fitt,plot.trace)

summary()函數(shù)可查看模型擬合后的結(jié)果，包括參數(shù)的估計(jì)值，以及相對(duì)應(yīng)的標(biāo)準(zhǔn)差等。plot()函數(shù)可以繪制相關(guān)的曲線圖。若plot.trace=FALSE，則繪制種群的年齡別死亡率和生存率曲線圖；若plot.trace=TRUE，則繪制出生存參數(shù)的迭代過(guò)程。

實(shí)例分析

1.數(shù)據(jù)來(lái)源

本文數(shù)據(jù)來(lái)源于2007年1月至2015年12月上海曙光醫(yī)院住院部腸癌患者的就診記錄。一共搜集到原始數(shù)據(jù)33118條。將搜集到的數(shù)據(jù)用病人的患者主索引號(hào)(EMPI)作為唯一的身份識(shí)別標(biāo)志，且以12個(gè)月為一個(gè)捕獲周期，借助Python軟件對(duì)原始數(shù)據(jù)進(jìn)行處理。整理得到2050例腸癌患者的捕獲-再捕獲數(shù)據(jù)，其中僅有151例可以明確地獲知死亡年齡。一共有女性患者924例，男性患者1126例；直腸癌患者為932例，結(jié)腸癌患者為939例，其他腸癌患者為179例。

2.數(shù)據(jù)分析

由于所獲得的數(shù)據(jù)中，年齡小于30歲的住院患者樣本數(shù)目較少，故本次年齡別生存從最小年齡30歲開(kāi)始進(jìn)行。

(1) 腸癌患者整體年齡別生存分析

根據(jù)multibasta()函數(shù)運(yùn)行的結(jié)果，選擇Gompertz函數(shù)[4]來(lái)對(duì)腸癌患者的整體數(shù)據(jù)進(jìn)行擬合，則年齡別死亡率函數(shù)可以表示為μb(x|b)=eb0+b1x，相對(duì)應(yīng)的生存函數(shù)可以表示為Sb(x|b)=exp[eb0(1-eb1x)/b1]。其中，eb0表示初始死亡率，b1表示死亡率隨著年齡呈指數(shù)增長(zhǎng)的速度。假定捕獲概率隨時(shí)間而變化，構(gòu)造兩條馬爾科夫鏈，且各自迭代40000次。運(yùn)行程序如下：

fit <-basta(data,studyStart=2007,studyEnd=2015,minAge=30,

model=“GO”,shape=“simple”,

niter=40000,burnin=20000,thinning=50,

nsim=2)

程序運(yùn)行之后，得到的生存模型的參數(shù)迭代過(guò)程如圖1所示。從圖中可以觀察到，兩條馬爾科夫鏈最后收斂于同一值。去掉最開(kāi)始的20000次迭代后，得到的參數(shù)估計(jì)值分別為b0=-5.30105，b1=0.04982，對(duì)應(yīng)的標(biāo)準(zhǔn)差分別為0.141918，0.003668，可信水平95%下的可信區(qū)間分別為(-5.59110，-5.05633)，(0.04340，0.05695)。由于參數(shù)的標(biāo)準(zhǔn)差都很小，且可信區(qū)間都很短，即波動(dòng)性很小，因此參數(shù)的估計(jì)是合理的。腸癌患者的年齡別死亡率和生存率曲線如圖2和圖3所示。

圖1 腸癌患者整體生存模型參數(shù)的迭代過(guò)程

圖2 腸癌患者整體年齡別死亡率曲線圖

圖3 腸癌患者整體年齡別生存率曲線圖

(2) 不同性別年齡別生存分析

選擇Gompertz函數(shù)來(lái)對(duì)女性腸癌患者的數(shù)據(jù)進(jìn)行擬合，得到的參數(shù)估計(jì)值為b0=-5.39197，b1=0.05283，對(duì)應(yīng)的標(biāo)準(zhǔn)差分別為0.178862，0.004184，可信水平95%下的可信區(qū)間分別為(-5.72823，-5.01510)，(0.04465，0.06090)。經(jīng)檢驗(yàn)，參數(shù)的估計(jì)是合理的。

不同性別的年齡別死亡率和生存率曲線分別如圖4和圖5所示。其中實(shí)線代表為女性患者，虛線代表男性患者。

圖4 不同性別腸癌患者的年齡別死亡率曲線圖

圖5 不同性別腸癌患者的年齡別生存率曲線圖

(3) 直腸和結(jié)腸癌患者年齡別生存分析

選擇logistic函數(shù)加上Makeham形狀的組合模型來(lái)對(duì)結(jié)腸癌患者的數(shù)據(jù)進(jìn)行擬合，得到的參數(shù)估計(jì)值為c=0.009167，b0=-6.037297，b1=0.060407，b2=0.159034，對(duì)應(yīng)的標(biāo)準(zhǔn)差分別為0.004261，0.621924，0.011753，0.136671，可信水平下的95%可信區(qū)間分別為(0.00122，0.01725)，(-7.10853，-4.50925)，(0.033934，0.08072)，(0.00588，0.50675)。 經(jīng)檢驗(yàn)，參數(shù)的估計(jì)是合理的。

直腸和結(jié)腸癌患者的年齡別死亡率和生存率曲線如圖6和圖7所示。其中實(shí)線代表為直腸癌患者，虛線代表結(jié)腸癌患者。

圖6 不同患癌部位腸癌患者的年齡別死亡率曲線圖

圖7 不同患癌部位腸癌患者的年齡別生存率曲線圖

3.結(jié)論分析

腸癌患者整體死亡率隨著年齡的增長(zhǎng)而增長(zhǎng)，其死亡中位年齡為71.55歲。不同性別腸癌患者的死亡率均隨著年齡的增長(zhǎng)而增長(zhǎng)，其中男性患者的年齡別死亡率高于女性患者，這與李德錄[7]，王寧[8]等人得到的結(jié)論相同。死亡中位年齡分別為女性71.65歲，男性69.25歲。

不同患癌部位的腸癌患者死亡率均隨著年齡的增長(zhǎng)而增長(zhǎng)，其中大約在60歲之前，直腸癌患者的年齡別死亡率低于結(jié)腸癌患者，在60歲之后，直腸癌患者的年齡別死亡率高于結(jié)腸癌患者，得到的結(jié)果與李德錄[7]大致相同。死亡中位年齡分別為直腸癌71.25歲，直腸癌68.05歲。

本文通過(guò)實(shí)例介紹了R軟件BaSTA包在貝葉斯框架下，對(duì)捕獲-再捕獲數(shù)據(jù)集進(jìn)行年齡別生存分析的方法。借助該包對(duì)患者進(jìn)行生存分析，可有效分析現(xiàn)有的臨床醫(yī)學(xué)數(shù)據(jù)，給出患者的生存參數(shù)估計(jì)值，并進(jìn)行可視化展示，具有廣泛的應(yīng)用價(jià)值。