R軟件BaSTA包在惡性腫瘤患者生存分析中的應(yīng)用*

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

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

生存模型及BaSTA包簡介

1.生存模型簡介

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

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

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

(1)

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

(2)

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

F(x|θ)=Pr(X

(3)

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

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

(4)

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

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

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

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

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

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

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

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

2.BaSTA包簡介

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

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

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

實際操作BaSTA包主要從以下幾個步驟考慮:

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

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

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

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

表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分別表示捕獲-再捕獲實驗開始和結(jié)束的時間。

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

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

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

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

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

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

summary(fit)

plot(fitt,plot.trace)

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

實例分析

1.數(shù)據(jù)來源

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

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

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

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

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

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

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

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

nsim=2)

程序運行之后，得到的生存模型的參數(shù)迭代過程如圖1所示。從圖中可以觀察到，兩條馬爾科夫鏈最后收斂于同一值。去掉最開始的20000次迭代后，得到的參數(shù)估計值分別為b0=-5.30105，b1=0.04982，對應(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ū)間都很短，即波動性很小，因此參數(shù)的估計是合理的。腸癌患者的年齡別死亡率和生存率曲線如圖2和圖3所示。

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

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

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

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

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

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

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

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

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

選擇logistic函數(shù)加上Makeham形狀的組合模型來對結(jié)腸癌患者的數(shù)據(jù)進行擬合，得到的參數(shù)估計值為c=0.009167，b0=-6.037297，b1=0.060407，b2=0.159034，對應(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)檢驗，參數(shù)的估計是合理的。

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

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

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

3.結(jié)論分析

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

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

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