模擬產(chǎn)生具有相關(guān)關(guān)系的終點(diǎn)指標(biāo)數(shù)據(jù)及R軟件實(shí)現(xiàn)*

張玲云 張曉敏 卓仁杰 劉麗亞△

腫瘤藥物療效的金標(biāo)準(zhǔn)評(píng)價(jià)指標(biāo)為總生存期(overall survival,OS),但是應(yīng)用該指標(biāo)作為"最終"臨床研究終點(diǎn)會(huì)使整個(gè)研發(fā)周期延長(zhǎng)、研發(fā)成本提高。實(shí)際工作中，常?？紤]應(yīng)用有效的短期替代終點(diǎn)代替總生存期作為腫瘤藥物療效評(píng)價(jià)的"最終"臨床研究終點(diǎn)。常見(jiàn)腫瘤臨床試驗(yàn)替代終點(diǎn)包括腫瘤大小，無(wú)進(jìn)展生存期(progression-free survival,PFS)，至疾病進(jìn)展生存期(time-to progression,TTP)及無(wú)病生存期(disease-free survival,DFS)等[1-2]。目前，關(guān)于如何驗(yàn)證替代終點(diǎn)是否能夠有效替代真實(shí)終點(diǎn)仍存在許多爭(zhēng)議，但是，毋庸置疑的是替代終點(diǎn)與總生存期具有強(qiáng)相關(guān)關(guān)系是替代終點(diǎn)有效的必要條件[3-4]。一項(xiàng)研究表明，2005年到2007美國(guó)食品藥品監(jiān)督管理局(the U.S.Food and Drug administration,FDA)批準(zhǔn)通過(guò)的53種腫瘤新藥中有9種(17%)腫瘤臨床試驗(yàn)的研究終點(diǎn)指標(biāo)為PFS[5]。由此可知腫瘤藥物臨床試驗(yàn)研究當(dāng)中不同結(jié)局變量間普遍存在相關(guān)關(guān)系[6]。因此，關(guān)于如何模擬產(chǎn)生具有相關(guān)關(guān)系的生存時(shí)間數(shù)據(jù)是腫瘤臨床試驗(yàn)?zāi)M研究亟待解決的問(wèn)題之一。

本文介紹并討論了三種模擬產(chǎn)生具有相關(guān)關(guān)系生存時(shí)間數(shù)據(jù)的方法，并通過(guò)實(shí)例介紹其R軟件實(shí)現(xiàn)[7]。

方法與原理

1.基于Copula函數(shù)法

目前，快速發(fā)展的Copula理論是構(gòu)建具有相關(guān)非正態(tài)變量聯(lián)合分布函數(shù)的有效手段之一，Copula函數(shù)實(shí)際上是一類(lèi)將聯(lián)合分布函數(shù)與它們各自的邊緣分布函數(shù)連接在一起的函數(shù)，因此也有人將它稱為連接函數(shù)。Copula函數(shù)的提出最早可以追溯到1959年，SKlar通過(guò)定理形式將多元分布與Copula函數(shù)聯(lián)系起來(lái)[8]。

Copula函數(shù)種類(lèi)很多，如Gaussian，t，Clayton，F(xiàn)rank，Gumbel，Plackett Copula等，其中由于Gaussian Copula函數(shù)可根據(jù)變量邊緣分布及相關(guān)系數(shù)矩陣唯一確定變量的聯(lián)合分布函數(shù)，故采用Gaussian Copula函數(shù)產(chǎn)生具有相關(guān)關(guān)系的生存時(shí)間數(shù)據(jù)可以根據(jù)實(shí)際情況設(shè)定相關(guān)系數(shù)ρ[9]。

以相關(guān)的二元分布為例介紹如何應(yīng)用Gaussian Copula函數(shù)構(gòu)建具有相關(guān)關(guān)系的二元非正態(tài)變量x1和x2的聯(lián)合分布函數(shù)。

(1)模擬產(chǎn)生樣本量為N相關(guān)系數(shù)為ρ的二元正態(tài)分布(Z1,Z2)，

(1)

(2)分別計(jì)算兩變量Z1和Z2的累積分布用U1和U2表示，

U1=Φ(Z1)

U2=Φ(Z2)

(2)

(U1,U2)表示具有相關(guān)關(guān)系的二元均勻分布，范圍分別為(0,1)。

(3)假定待估計(jì)具有相關(guān)關(guān)系的二元非正態(tài)變量X1和X2的邊際分布函數(shù)分別為，F(xiàn)1(·),F2(·)則

(3)

假定生存時(shí)間服從指數(shù)分布，同時(shí)考慮截尾時(shí)間變量，則根據(jù)Gaussian Copula函數(shù)模擬產(chǎn)生具有相關(guān)關(guān)系的生存時(shí)間變量的具體步驟如下，

(1)模擬產(chǎn)生樣本量為N相關(guān)系數(shù)為ρ的二元正態(tài)分布(Z1,Z2)；

(2)依次計(jì)算二元正態(tài)分布(Z1,Z2)每個(gè)變量的累積分布，用(U1,U2)表示，即二元均勻分布，范圍是(0,1)；

(3)生存時(shí)間服從指數(shù)分布，則生存時(shí)間X的分布函數(shù)及其反函數(shù)分別為F(t)=1-e-θt，F(xiàn)-1(t)=-ln(1-F1(t))/θ1。將U1看作是1-F1(t)代入變量分布函數(shù)的反函數(shù)估計(jì)生存時(shí)間X。同理，估計(jì)生存時(shí)間Y；

(4)模擬產(chǎn)生截尾時(shí)間變量C，C～Uniform[0,a]；

(5)根據(jù)上述步驟模擬產(chǎn)生的(X,Y,C)及公式Zx=min(X,Y,C),Zy=min(Y,C),δx=I(X≤min(Y,C))、δy=I(Y≤C)產(chǎn)生TTP(即Zx)與OS(即Zy)變量以及相應(yīng)截尾變量，I表示指示函數(shù)。

例1 生成1000個(gè)相關(guān)關(guān)系為0.5的TTP與OS，要求TTP與OS的中位生存時(shí)間分別為2個(gè)月、4.5個(gè)月且截尾時(shí)間變量服從均勻分布[0,36]。

說(shuō)明：題目設(shè)定的相關(guān)系數(shù)是Copula函數(shù)指定的二元正態(tài)分布相關(guān)系數(shù)，而非TTP與OS相關(guān)系數(shù)。

程序如下：

install.packages("mvtnorm")

library(mvtnorm)

n=1000 #樣本量為1000

corr=0.5 #相關(guān)關(guān)系為0.5

sigma=matrix(c(1,corr,corr,1),ncol=2) #相關(guān)系數(shù)矩陣

x=rmvnorm(n=n,mean=c(0,0),sigma=sigma) #模擬產(chǎn)生二元正態(tài)分布隨機(jī)變量

p1=pnorm(x[,1],mean=0,sd=1,lower.tail=TRUE,log.p=FALSE) #變量Z1的累積分布函數(shù)U1

p2=pnorm(x[,2],mean=0,sd=1,lower.tail=TRUE,log.p=FALSE) #變量Z2的累積分布函數(shù)U2，

TTP0=qgamma(p1,1,rate=log(2)/2,lower.tail=TRUE,log.p=FALSE) #TTP服從中位生存時(shí)間為2的指數(shù)分布

OS0=qgamma(p2,1,rate=log(2)/4.5,lower.tail=TRUE,log.p=FALSE) #OS服從中位生存時(shí)間為4.5的指數(shù)分布

censor=runif(n=n,min=0,max=36) #截尾變量censor服從均勻分布[0,36]

OS=pmin(OS0,censor) #實(shí)際的OS為OS0與censor的最小值

TTP=pmin(TTP0,OS0,censor) #實(shí)際的TTP為T(mén)TP0,OS0和censor的最小值

TTP_cen=(TTP0<=OS) #如TTP0小于等于實(shí)際的OS則認(rèn)為疾病進(jìn)展事件發(fā)生，表示疾病進(jìn)展在死亡前、觀察截止前被觀測(cè)到

OS_cen=(OS0<=censor) #如OS0小于等于截尾時(shí)間則認(rèn)為終點(diǎn)事件發(fā)生，表示死亡事件在觀察截止前被觀測(cè)到

Data_EX1=cbind(TTP,OS,TTP_cen,OS_cen) #最終生成的數(shù)據(jù)集

另上述程序假定生存時(shí)間均服從指數(shù)分布，如生存時(shí)間服從Weibull分布的，上述生成TTP0、OS0的步驟可修改為，

TTP0=qweibull(p1,1,2/log(2),lower.tail=TRUE,log.p=FALSE) #TTP服從中位生存時(shí)間為2的Weibull分布

OS0=qweibull(p2,1,4.5/log(2),lower.tail=TRUE,log.p=FALSE) #OS服從中位生存時(shí)間為4.5的Weibull分布。

2.混合指數(shù)模型法

腫瘤病人在死亡前常經(jīng)歷數(shù)次疾病進(jìn)展，如圖1所示，然而疾病進(jìn)展后的死亡風(fēng)險(xiǎn)一般高于疾病進(jìn)展前，因此構(gòu)建產(chǎn)生生存時(shí)間數(shù)據(jù)模型時(shí)應(yīng)充分考慮疾病發(fā)生發(fā)展機(jī)制。采用混合指數(shù)模型法分別模擬產(chǎn)生疾病進(jìn)展前、后生存時(shí)間，不同階段的死亡風(fēng)險(xiǎn)可依據(jù)實(shí)際情況設(shè)定不同參數(shù)[10]。

圖1 腫瘤患者疾病進(jìn)展示意

(4)

如果k≠i則λi≠λk。

假定受試者只經(jīng)歷一次疾病進(jìn)展后死亡，即n=2時(shí)，t的概率密度函數(shù)可表示為:

(5)

則根據(jù)混合指數(shù)模型模擬產(chǎn)生具有相關(guān)關(guān)系生存時(shí)間數(shù)據(jù)的具體步驟如下:

(1)模擬產(chǎn)生樣本量為N，參數(shù)為λ1的指數(shù)分布，TTP；

(2)模擬產(chǎn)生樣本量為N，參數(shù)為λ2的指數(shù)分布，PPS；

(3)OS=TTP+PPS，則將步驟(1)與步驟(2)產(chǎn)生的模擬數(shù)據(jù)相加得到OS。

基于該法不能考慮受試者在疾病進(jìn)展前發(fā)生死亡事件，F(xiàn)leischer等人于 2009年對(duì)該法進(jìn)行了改進(jìn)[11]。首先，假定TTP與OS存在相關(guān)關(guān)系，TTP服從參數(shù)為λ1的指數(shù)分布。其次，引入X變量，表示未經(jīng)歷疾病進(jìn)展的患者從隨機(jī)化分組至死亡的時(shí)間，服從參數(shù)為λ2的指數(shù)分布，X～Exp(λ2)。由于TTP與X變量相互獨(dú)立，設(shè)定PFS=min(TTP,X)，則PFS仍服從指數(shù)分布，PFS～Exp(λ1+λ2)。最后，假定疾病進(jìn)展后生存期(post-progression survival,PPS)仍服從參數(shù)為λ3指數(shù)分布，PPS～Exp(λ3)。因此，OS可表示為，

(6)

此時(shí)OS風(fēng)險(xiǎn)函數(shù)不再是一個(gè)常量，隨著疾病進(jìn)展而改變。一般來(lái)講，疾病進(jìn)展后的風(fēng)險(xiǎn)率高于進(jìn)展前(λ3>λ2，即疾病進(jìn)展前死亡風(fēng)險(xiǎn)為λ2，疾病進(jìn)展后死亡風(fēng)險(xiǎn)為λ3)。基于上述假定，PFS與OS的相關(guān)關(guān)系為:

(7)

在其他參數(shù)不變前提下，隨著λ1或λ2的增加PFS與OS的相關(guān)關(guān)系下降，隨著λ3的增加PFS與OS的相關(guān)關(guān)系增加，見(jiàn)圖2。

圖2 PFS與OS在不同參數(shù)組合下相關(guān)關(guān)系情況(lamda2=1)

*：lamda1，表示TTP服從參數(shù)為lamda1的指數(shù)分布；lamda2，表示未經(jīng)歷疾病進(jìn)展患者從隨機(jī)化分組至死亡的時(shí)間服從參數(shù)為lamda2的指數(shù)分布；lamda3，表示PPS服從參數(shù)為lamda3的指數(shù)分布；correlation表示PFS與OS相關(guān)系數(shù)。

例2 模擬產(chǎn)生1000個(gè)具有相關(guān)關(guān)系的PFS與OS數(shù)據(jù)，要求TTP、PPS的中位時(shí)間分別為5個(gè)月及15個(gè)月，未經(jīng)歷疾病進(jìn)展患者的中位生存時(shí)間為7.5個(gè)月。

程序如下：

n=1000 #樣本量為1000

rate1=log(2)/7.5 #未經(jīng)歷疾病進(jìn)展患者的死亡風(fēng)險(xiǎn)函數(shù)

rate2=log(2)/5 #TTP的風(fēng)險(xiǎn)函數(shù)

rate3=log(2)/15 #PPS的風(fēng)險(xiǎn)函數(shù)

x=rexp(n,rate1) #模擬產(chǎn)生未經(jīng)歷疾病進(jìn)展患者的死亡時(shí)間

ttp=rexp(n,rate2) #模擬產(chǎn)生TTP

pfs=pmin(ttp,x) #模擬產(chǎn)生PFS

pps=rexp(n,rate3) #模擬產(chǎn)生PPS

os=pps+pfs #模擬產(chǎn)生OS

cor(pfs,os) #估計(jì)PFS與OS相關(guān)關(guān)系

3.混合Bayes模型法

Michael等人于2002年提出一種能夠精確模擬產(chǎn)生具有相關(guān)關(guān)系的多元分布新方法，稱為混合Bayes模型法[12]。該法應(yīng)用Bayes思想，易于推廣，其本質(zhì)是根據(jù)后驗(yàn)分布模擬產(chǎn)生與先驗(yàn)具有相同邊際分布的樣本，兩個(gè)邊際分布的相關(guān)關(guān)系由似然函數(shù)決定。

假定隨機(jī)變量X1的先驗(yàn)分布為

g(x1;θ)

(8)

參數(shù)θ表示多維變量；在X1=x條件下，根據(jù)抽樣分布K|(X1=x)獲得樣本，構(gòu)建似然函數(shù)可表示為

h(k|x1;η)

(9)

將公式(8)與公式(9)相乘，得變量X1和K聯(lián)合分布

j(x1,k;θ,η)=g(x1;θ)h(k|x1;η)

(10)

對(duì)X1積分得K的邊際分布

(11)

給定K=k條件下，估計(jì)X2的后驗(yàn)分布

p(x2|k;θ,η)=j(x2,k;θ,η)/m(k;θ,η)

(12)

變量X1與X2有相同的邊際分布，X2分布為不同后驗(yàn)分布的加權(quán)平均，準(zhǔn)確地再現(xiàn)了先驗(yàn)分布。假定生存時(shí)間滿足Gamma分布，則變量X1與X2的邊際分布均為Gamma分布。

應(yīng)用混合Bayes模型法模擬產(chǎn)生具有相關(guān)關(guān)系的生存時(shí)間變量的具體步驟如下：

(1)確定先驗(yàn)分布X1～Γ(a,b)，根據(jù)先驗(yàn)分布模擬產(chǎn)生X1，樣本量為N;

(2)由于Poisson分布的共軛先驗(yàn)分布是Gamma分布，因此條件抽樣分布為K|(X1=x)～Poisson(τx)，據(jù)此分布抽樣獲得樣本量為N的K；

(3)根據(jù)先驗(yàn)以及抽樣獲得的K估計(jì)后驗(yàn)分布，X2～Γ(a+k,b+τ)；

(4)根據(jù)每個(gè)后驗(yàn)分布依次模擬產(chǎn)生X2，據(jù)此獲得的二維隨機(jī)變量(X1,X2)具有相關(guān)關(guān)系且邊際分布均為Gamma分布，相關(guān)系數(shù)為

(13)

(5)模擬產(chǎn)生截尾時(shí)間變量C，C～Uniform[0,a]；

(6)根據(jù)上述步驟模擬產(chǎn)生的(X1,X2,C)及公式Zx1=min(X1,X2,C),Zx2=min(X2,C),δx1=I(X1≤min(X2,C))，δx2=I(X2≤C)產(chǎn)生TTP(即Zx1)與OS(即Zx2)變量以及相應(yīng)截尾變量，I表示指示函數(shù)。

例3 模擬產(chǎn)生1000個(gè)相關(guān)系數(shù)為0.78的TTP與OS，要求TTP服從指數(shù)分布中位時(shí)間為5個(gè)月。

說(shuō)明：依據(jù)實(shí)際情況需要模擬產(chǎn)生生存時(shí)間數(shù)據(jù)需考慮截尾時(shí)間。本文介紹混合Bayes模型法模擬產(chǎn)生截尾時(shí)間變量的步驟同例1一致，故R程序一致，同時(shí)由于考慮截尾時(shí)間將會(huì)影響TTP與OS的相關(guān)關(guān)系，故本例未考慮截尾時(shí)間變量。本例模擬數(shù)據(jù)產(chǎn)生后可評(píng)估TTP與OS的相關(guān)系數(shù)與設(shè)定參數(shù)是否一致，從而檢驗(yàn)?zāi)M過(guò)程是否正確。

程序如下：

n=1000 #樣本量為1000

shape1=1 #形態(tài)參數(shù)為1的Gamma分布是指數(shù)分布

rate1=log(2)/5 #尺度參數(shù)

prior=rgamma(n,shape1,scale=1/rate1) #先驗(yàn)分布為指數(shù)分布，即TTP

tau=0.5 #Poisson分布的τ設(shè)定為0.5

lamda=data.frame(tau*prior) #估計(jì)條件分布K|(X1=x)～Poisson(τx)參數(shù)

likelihood=rpois(n,lamda[,1]) #據(jù)條件分布抽樣獲得K，樣本量為1000

shape2=data.frame(shape1+likelihood) #根據(jù)先驗(yàn)形態(tài)參數(shù)以及抽樣獲得的K估計(jì)1000個(gè)后驗(yàn)分布的形態(tài)參數(shù)

rate2=1/(1/rate1+tau) #依據(jù)相關(guān)關(guān)系公式估計(jì)后驗(yàn)分布的尺度參數(shù)

posterior=rgamma(n,shape2[,1],scale=1/rate2) #據(jù)1000個(gè)后驗(yàn)分布抽樣獲得1000個(gè)后驗(yàn)樣本，即OS

cor(prior,posterior) #估計(jì)TTP與OS相關(guān)系數(shù)

表1 三種模擬產(chǎn)生具有相關(guān)關(guān)系生存時(shí)間數(shù)據(jù)方法的優(yōu)缺點(diǎn)

*：PFS，progression-free survival，無(wú)進(jìn)展生存期；TTP，time-to progression，至疾病進(jìn)展生存期；PPS，post-progression survival，疾病進(jìn)展后生存期。

討 論

關(guān)于如何模擬產(chǎn)生具有相關(guān)關(guān)系的生存時(shí)間數(shù)據(jù)廣泛應(yīng)用于研究與實(shí)踐領(lǐng)域，主要用于評(píng)估統(tǒng)計(jì)分析方法的統(tǒng)計(jì)學(xué)性能。本文介紹了三種模擬產(chǎn)生具有相關(guān)系數(shù)生存時(shí)間數(shù)據(jù)的方法，其相關(guān)系數(shù)參數(shù)均能設(shè)定，但是由于截尾數(shù)據(jù)的存在，模擬設(shè)定的相關(guān)系數(shù)將會(huì)受到影響，且隨截尾變量均勻分布上限的縮短而降低。

在實(shí)際工作中應(yīng)用三種方法產(chǎn)生模擬數(shù)據(jù)時(shí)需注意：(1)Copula函數(shù)法設(shè)定的相關(guān)系數(shù)是二元正態(tài)分布相關(guān)系數(shù)而非生存時(shí)間相關(guān)；(2)混合指數(shù)模型法可以控制PPS的長(zhǎng)短，如果模擬試驗(yàn)研究考慮該變量作為研究的變化參數(shù)時(shí)，則只能采用該法產(chǎn)生模擬數(shù)據(jù)；(3)混合Bayes模型法不能設(shè)定中位生存時(shí)間；(4)Copula函數(shù)法和混合Bayes模型法除能模擬產(chǎn)生具有相關(guān)關(guān)系的生存時(shí)間變量外，也可推廣應(yīng)用至其他分布類(lèi)型數(shù)據(jù)的模擬產(chǎn)生，但是應(yīng)用混合Bayes模型法要求具有相關(guān)關(guān)系的先驗(yàn)與后驗(yàn)分布類(lèi)型一致。

本研究為建立一套科學(xué)、可行、又可操作的產(chǎn)生具有相關(guān)關(guān)系的模擬數(shù)據(jù)的方式提供依據(jù)，進(jìn)而解決科研工作者在生存資料模擬研究中關(guān)于產(chǎn)生模擬數(shù)據(jù)的困惑。本研究的局限性在于產(chǎn)生截尾時(shí)間數(shù)據(jù)的方式均基于均勻分布且截尾類(lèi)型僅考慮右側(cè)截尾情況，因此在實(shí)際工作中產(chǎn)生模擬數(shù)據(jù)的方式需依據(jù)實(shí)際情況而定。