南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)系(211166)

蔡麗馨# 仲子航# 楊 旻 于全驥 周佳薇 倪森淼 于 浩△ 柏建嶺△

【提 要】 目的 探討生成具有指定刪失比例的模擬數(shù)據(jù)的新方法，并編寫完整的模擬數(shù)據(jù)產(chǎn)生的R代碼，方便使用。方法 基于Cox比例風(fēng)險模型框架分別推導(dǎo)考慮和不考慮協(xié)變量情況下的刪失參數(shù)的閉式表達(dá)式或核密度估計，并基于R 4.0.2軟件編寫函數(shù)模擬產(chǎn)生滿足指定刪失比例的生存數(shù)據(jù)。通過1000次模擬數(shù)據(jù)來驗證生成數(shù)據(jù)的刪失比例是否與指定的刪失比例一致。結(jié)果 函數(shù)CenDatNoCov、CenDatBin、CenDatNorm和CenDatMixed可分別返回不考慮協(xié)變量，協(xié)變量為二項分布、正態(tài)分布以及混合分布情況下的刪失數(shù)據(jù)，其中后三者內(nèi)嵌的函數(shù)CensProbBin，CensProbNorm，CensProbMixed可分別得到相應(yīng)情況下的刪失參數(shù)值。模擬生成數(shù)據(jù)的刪失比例與指定的刪失比例近似相等。結(jié)論 本研究所提供的方法及編寫的R函數(shù)可有效生成指定刪失比例的數(shù)據(jù)。

生存數(shù)據(jù)的相關(guān)模擬近年來越來越受到醫(yī)學(xué)研究者的關(guān)注。一方面，隨著臨床試驗的設(shè)計方法越來越復(fù)雜，《藥物臨床試驗適應(yīng)性設(shè)計指導(dǎo)原則(征求意見稿)》[1]中明確指出了統(tǒng)計模擬試驗的優(yōu)點及重要性，而在許多臨床試驗，尤其是腫瘤藥物研究中，主要終點指標(biāo)通常為生存數(shù)據(jù)，如OS，PFS等[2-5]。另一方面，生存數(shù)據(jù)統(tǒng)計方法的研究也依賴于模擬試驗對相關(guān)性能進(jìn)行評價[6-7]。對于生存數(shù)據(jù)而言，在實際情況下，受人力、財力、資源等限制，通常會有部分受試者的終點事件在研究結(jié)束時仍未被觀察到，從而被標(biāo)記為“刪失”。研究表明刪失會帶來許多問題，同時對不同的分析方法產(chǎn)生一定的影響[6，8]，因此，為反映所研究方法的真實性能，能夠有效控制刪失比例在預(yù)設(shè)的水平顯得十分必要。針對這一需求，我們通過文獻(xiàn)檢索發(fā)現(xiàn)目前國內(nèi)對此研究較少，有研究者提出在產(chǎn)生的完整數(shù)據(jù)中隨機抽取指定比例的觀測作為刪失，并對此部分觀測從新的刪失分布中產(chǎn)生刪失時間作為最終的生存數(shù)據(jù)，這一方法較為理想[9]。更常用的方法是基于特定的分布同時生成事件發(fā)生時間和刪失時間兩個變量，通過邏輯判斷來決定最終狀態(tài)，在此框架下，一些學(xué)者通過選取不同參數(shù)值進(jìn)行重復(fù)手工調(diào)試以確定平均刪失比例滿足要求的參數(shù)值[10]，這一方法對計算機性能的要求很高，同時會消耗較長的時間；Fei Wan等人提出一種方法[11]，基于Cox比例風(fēng)險模型在考慮了刪失分布、協(xié)變量分布及基線風(fēng)險后可以給出分布參數(shù)值的閉式表達(dá)，從而使得模擬的生存數(shù)據(jù)具有指定的刪失比例。本文基于此方法，給出不考慮協(xié)變量和考慮協(xié)變量的不同情形下產(chǎn)生指定刪失比例的生存數(shù)據(jù)的方法，同時編寫了完整的R代碼以供使用。

方法與原理

為了模擬刪失數(shù)據(jù)，我們常用的一種方法即生成事件發(fā)生時間T和刪失時間C兩個變量，并記δ=I(T≥C)為刪失狀態(tài)變量，即當(dāng)T

事件發(fā)生時間通?？梢酝ㄟ^指數(shù)分布、Weibull分布、Gamma分布、log-normal分布及l(fā)og-logistic分布等來描述，其中Weibull分布Weibull(α，λ)相較于Gamma分布、log-normal分布等形式更為簡潔，相較于指數(shù)分布又可通過額外的尺度參數(shù)(scale parameter)λ和靈活的形狀參數(shù)(shape parameter)α取值來彌補指數(shù)分布的局限性，形狀參數(shù)的不同取值可以描述在研究期間風(fēng)險的變化趨勢而不再是指數(shù)分布假設(shè)下單純的風(fēng)險恒定，因此其應(yīng)用最為廣泛。

刪失時間通常認(rèn)為服從均勻分布uniform(0，θ)，指數(shù)分布exp(θ)或者Weibull分布Weibull(α，θ)，其中θ為所謂的刪失參數(shù)。

本文基于事件發(fā)生時間服從Weibull(α，λ)分布，刪失時間服從于均勻分布uniform(0，θ)來模擬生成指定刪失比例的生存數(shù)據(jù)，考慮協(xié)變量時采用Cox模型。

1.不考慮協(xié)變量的指定刪失比例生存數(shù)據(jù)

(1)事件發(fā)生時間T的參數(shù)分布

Weibull分布的概率密度函數(shù)和風(fēng)險函數(shù)為：

(1)

h(t；α，λ)=αλ-αtα-1

(2)

其中，λ>0，α>0，t≥0

(2)滿足指定刪失比例的刪失時間C的參數(shù)分布

由于所謂刪失參數(shù)θ的取值可直接影響最終的刪失比例p，因此我們在給定一系列參數(shù)及相關(guān)分布的情況下構(gòu)建p與θ的關(guān)系，對θ值進(jìn)行反推以滿足需求。設(shè)定刪失時間C服從均勻分布uniform(0，θ)，即其密度函數(shù)為：

(3)

各個體的刪失概率為：

(4)

2.考慮協(xié)變量的指定刪失比例生存數(shù)據(jù)

(1)事件發(fā)生時間T的參數(shù)分布

基于Cox風(fēng)險比例模型考慮協(xié)變量的影響。對于第i例受試者，其協(xié)變量記為d維向量Xi，回歸系數(shù)向量記為β，基線風(fēng)險函數(shù)記為h0(t)，則風(fēng)險函數(shù)可以表示為：

(5)

上式中，exp(X′iβ)λ-α=[λexp(-X′β/α)]-α，記λi=λ[exp(-X′iβ/α)]=exp[-(-αlogλ+X′iβ)/α]，則考慮協(xié)變量影響所對應(yīng)的事件發(fā)生時間將服從Weibull(α，λi)，密度函數(shù)為：

(6)

(7)

由此可構(gòu)建Xi與λi間的關(guān)系，將復(fù)雜的多維問題簡化為對單變量的處理。不難發(fā)現(xiàn)，當(dāng)Xi中的變量類型相同時，fλ(u)可以準(zhǔn)確給出。

(8)

(9)

當(dāng)Xi為多類型變量的組合時，我們無法準(zhǔn)確地給出該密度函數(shù)表達(dá)式，此時可以通過非參數(shù)核平滑方法得到λi的密度估計函數(shù)，如選用高斯核函數(shù)來估計。

(2)滿足指定刪失比例的刪失時間C的參數(shù)分布

對于刪失參數(shù)的求解步驟同上，但需考慮不同個體的刪失概率推導(dǎo)總體的刪失率。仍假設(shè)刪失時間C服從均勻分布uniform(0，θ)，即其密度函數(shù)為：g(c|θ)=1/θ，0

①描述各個體的刪失概率：

p(δ=1|Xi，θ)=p(C≤T<∞，0≤C<θ)

(10)

②基于個體水平的刪失概率估計總體的平均刪失率：

(11)

其中，fx(·)為d個獨立協(xié)變量X的聯(lián)合密度函數(shù)，D為X的概率空間。前面已介紹可通過估計以λi的概率密度fλi(u)來簡化計算。

③通過對等式γ(θ)=p(δ=1|θ)-p=0進(jìn)行求解得到θ值，其中：

(12)

若p(δ=1|λi，θ)有閉式表達(dá)式，則式(11)在簡單的積分同樣可得到明確的表達(dá)式進(jìn)而求得解析解，但若無法直接寫出該條件刪失概率的表達(dá)式，則需對式(11)進(jìn)行二重積分后求解。

在事件發(fā)生時間滿足Weibull分布的情況下，記其基線風(fēng)險函數(shù)h0(t)=αλ-αtα-1，相互獨立的協(xié)變量為X，則事件發(fā)生時間可認(rèn)為來自Weibull(α，λi)，刪失時間服從uniform(0，θ)，則可求得各個體刪失的概率為：

p(δ=1|Xi，θ)=p(δ=1|λi，θ)=p(C≤T<∞，0≤C<θ)

(13)

z=(c/λi)α

其中γ(1/α，(θ/λi)α)為下不完全Gamma函數(shù)。

給定回歸系數(shù)β=<β1，β2，β3>，記β0=-αlogλ：

①若X均為連續(xù)性變量，如4個均來自N(0，1)，則：

(14)

②若X均為二分類變量，如分別來自B(p1)，B(p2)，B(p3)，B(p4)，則

(15)

③若X為不同類型的變量，如分別來自4個不同分布，包括N(0，1)，B(p2)，Poisson(n3)，Uniform(0，1)，則需要通過核平滑估計來估計λi的概率密度fλi(u)并帶入計算。

R實現(xiàn)與應(yīng)用

1.不考慮協(xié)變量

我們構(gòu)建CenDatNoCov()函數(shù)來產(chǎn)生指定刪失比例的生存數(shù)據(jù)，輸入?yún)?shù)包括Weibull(α，λ)的形狀alpha和尺度參數(shù)lambda，指定的刪失比例cens.p及樣本量size。由于不考慮協(xié)變量，我們僅需通過uniroot()函數(shù)對式(6)進(jìn)行求根即可得到刪失參數(shù)theta，其中下不完全Gamma函數(shù)γ(a，x)在R中可定義為pgamma(x，a)*gamma(a)。

CenDatNoCov<- function(alpha，lambda，cens.p，size){

theta <- round(uniroot(function(x)lambda/(alpha*x)*pgamma((x/lambda)^alpha，1/alpha)*gamma(1/alpha)-cens.p，c(0.00001，1000))$root，3)

#據(jù)方法與原理的描述產(chǎn)生刪失的生存數(shù)據(jù)

data<-data.frame(

T=rweibull(size，alpha，lambda)，

C=runif(size，0，theta))

cens.data<- data.frame(time=ifelse(data$T<=data$C，data$T，data$C)，

status=ifelse(data$T<=data$C，1，0))

return(list(theta=theta，cens.data=cens.data))

}

例1 假設(shè)事件發(fā)生時間T來自Weibull(2，4)，刪失時間C服從uniform(0，θ)，指定的刪失比例為30%，樣本量為200，即CenDatNoCov(alpha=2，lambda=4，cens.p=0.3，size=200)。求得θ=11.816即為了使得最終的刪失比例達(dá)到30%，需設(shè)定刪失參數(shù)θ為11.816。為了檢驗該參數(shù)是否可以滿足需求，我們可以生成1000組T和C并進(jìn)行刪失狀態(tài)的判斷計算其平均刪失率，設(shè)定隨機種子數(shù)為20200810，結(jié)果為0.30107，符合我們的預(yù)設(shè)的刪失比例。

2.考慮協(xié)變量

(1)二分類協(xié)變量

我們構(gòu)建CenDatBin()來產(chǎn)生均為二分類協(xié)變量的刪失數(shù)據(jù)，輸入?yún)?shù)包括Weibull(α，λ)的形狀alpha和尺度參數(shù)lambda，指定的刪失比例cens.p，事件發(fā)生時間的風(fēng)險模型回歸系數(shù)beta，各協(xié)變量的分布參數(shù)p，樣本量size及隨機種子數(shù)seed；返回列表包括刪失參數(shù)theta及所生成的生存數(shù)據(jù)。其中內(nèi)置的CensProbBin()函數(shù)可返回對應(yīng)的刪失參數(shù)theta值。

CenDatBin<- function(alpha，lambda，cens.p，beta，p，size，seed=20200812){

alpha=alpha # weibull分布中的形狀參數(shù)

lambda=lambda # weibull分布中的形狀參數(shù)

cens.p=cens.p#預(yù)設(shè)的刪失比例

beta=beta# Cox模型中各協(xié)變量的系數(shù)

p=p# 各二分類協(xié)變量的發(fā)生率

n=size#產(chǎn)生數(shù)據(jù)的樣本量

#構(gòu)建刪失率函數(shù)，參數(shù)為均勻分布中的刪失參數(shù)θ，返回值為給定θ后估計的總體刪失率：

CensProbBin <- function(theta){

beta.0 <--alpha*log(lambda) #計算新常數(shù)以代入公式(7)

#得到協(xié)變量取值的所有組合

LenCovar <- length(beta) #獲得協(xié)變量個數(shù)

CombInd <-list(NULL)

ncomb <- c()

finalInd <- c() #所有組合中，取值為1的協(xié)變量下標(biāo)矩陣

for(i in 1：LenCovar){ #通過循環(huán)得到至少1個協(xié)變量取值為1的所有組合下標(biāo)

ind <- combn(1：LenCovar，i) #其中i個協(xié)變量取值為1的變量下標(biāo)組合

CombInd[[i]]<- ind

for(j in 1：ncomb[i]) {

location <- c(CombInd[[i]][，j]，rep(0，(LenCovar-i))) #以0補齊取值為0的協(xié)變量下標(biāo)

finalInd <- cbind(finalInd，location) #得到組合的下標(biāo)矩陣共LenCovar行，2^LenCovar列

}

}

# 根據(jù)各協(xié)變量的下標(biāo)組合進(jìn)行賦值，得到最終的取值組合

comb <- matrix(0，nrow = LenCovar，ncol = 2^LenCovar) #構(gòu)建協(xié)變量取值為0的矩陣

for(i in 1：2^LenCovar-1){

comb[finalInd[，i]，i]=1 #將下標(biāo)矩陣中不為0的協(xié)變量取值賦為1

}

# 據(jù)公式(9)構(gòu)建單變量，簡化計算

lambda.i<-exp(-(beta.0+apply(beta*comb，2，sum))/alpha)

#計算在對應(yīng)協(xié)變量下的個體條件刪失概率

cond.cens.Prob<-(lambda.i/(alpha*theta))*pgamma((theta/lambda.i)^alpha，1/alpha) * gamma(1/alpha)

#計算該變量的分布律

pdf.lambda.i<-apply(matrix(p，nrow=LenCovar，ncol = 2^LenCovar)**comb*((1-matrix(p，nrow=LenCovar，ncol = 2^LenCovar))**(1-comb))，2，prod)

return(t(cond.cens.Prob)%*%pdf.lambda.i)

}

#求解刪失參數(shù)值

theta<- round(uniroot(function(x)censProbBin(x)-cens.p，c(0.0000001，100))$root，3)

#生成模擬協(xié)變量值

set.seed(seed)

a<-c()

cov<-sapply(p，function(x)cbind(a，rbinom(n，1，x)))

colnames(cov)<- paste(“x”，1：length(p)，sep=“”)

#指定rweibull()函數(shù)中scale參數(shù)值為lambda*exp(-1/alpha*(Xβ))

EVENT<-rweibull(n，alpha，lambda*exp(-1/alpha*(cov%*%beta)))

data<-as.data.frame(cbind(id=1：n，cov=cov，EVENT=EVENT，C=runif(n，0，theta)))

#據(jù)方法與原理的描述產(chǎn)生刪失的生存數(shù)據(jù)

data$time=ifelse(data$EVENT<=data$C，data$EVENT，data$C)

data$status=ifelse(data$EVENT<=data$C，1，0)

return(list(theta=theta，cens.data=data))

}

例2 以四個協(xié)變量分別來自B(0.3)，B(0.4)，B(0.5)，B(0.6)為例。CenDatBin(alpha=2，lambda=4，cens.p=0.3，beta=c(-0.1，0.2，-0.3，0.4)，p=c(0.3，0.4，0.5，0.6)，size=200，seed=20200812)，得到θ=11.116，即為了達(dá)到30%的刪失比例，刪失時間的分布uniform(0，θ)的參數(shù)應(yīng)設(shè)為11.116?？赏ㄟ^上述方法生成1000次模擬數(shù)據(jù)來驗證結(jié)果，當(dāng)種子數(shù)設(shè)定為20200810時，得到的刪失比例平均值為0.30132，滿足預(yù)設(shè)的要求。

(2)正態(tài)分布協(xié)變量

我們構(gòu)建CenDatNorm()來產(chǎn)生均為正態(tài)分布協(xié)變量的刪失數(shù)據(jù)，輸入?yún)?shù)基本同上。內(nèi)嵌計算刪失參數(shù)的相關(guān)函數(shù)為CensProbNorm()，通過求根即可獲得theta值：

CensProbNorm<-function(theta){

beta.0 <- -alpha*log(lambda) # alpha、lambda分別為weibull分布中的形狀、尺度參數(shù)

# 據(jù)公式(8)得到單變量λi的概率密度函數(shù)，簡化計算

PdfLambdai<-function(u)dlnorm(u，-beta.0/alpha，beta%*%beta/alpha^2)

# 據(jù)公式(13)得到條件刪失概率

CondCensProb<-function(u)(u/(alpha*theta))*pgamma((theta/u)^alpha，1/alpha)* gamma(1/alpha)

# 據(jù)公式(14)得到刪失概率的表達(dá)式

cens.Prob<-integrate(function(u){PdfLambdai(u)*CondCensProb(u)}，-Inf，Inf)$value

return(cens.Prob)

}

例3 以四個來自N(0，1)的協(xié)變量為例，給定alpha=2，lambda=4，預(yù)定的刪失比例cens.p=0.3，求解對應(yīng)的刪失參數(shù)，theta<-round(uniroot(function(x)CensProbNorm(x)-cens.p，c(0.0001，1000))$root，3)，得θ=11.849，即為了達(dá)到30%的刪失率，刪失時間的分布uniform(0，θ)的參數(shù)應(yīng)設(shè)為11.849，同理，模擬1000次得到的平均刪失率為0.30831。

(3)混和協(xié)變量

我們構(gòu)建CenDatMixed()來產(chǎn)生混合分布協(xié)變量的刪失數(shù)據(jù)，輸入?yún)?shù)除同上的基本參數(shù)外，還包括描述所有混合變量類型的參數(shù)calss，各分布的對應(yīng)參數(shù)para。由于協(xié)變量類型不統(tǒng)一，我們很難推導(dǎo)出λi的概率密度fλi(u)，因此我們采用高斯核密度估計方法，刪失參數(shù)的估計可通過內(nèi)嵌的CensProbMixed()函數(shù)獲得：

CensProbMixed<- function(class，para，theta){

set.seed(20200810)

#高斯核密度估計lambda.i的密度函數(shù)

n<-10000 # 以10000個點來估計lambda.i的密度分布

# 生成構(gòu)成lambda.i的相關(guān)變量值以提供擬合聯(lián)合分布的數(shù)據(jù)

nNorm<- length(class[which(class==“N”)])

nBin<- length(class[which(class==“B”)])

nPos<- length(class[which(class==“P”)])

nUni<- length(class[which(class==“U”)])

if(nNorm> 0){

a<-as.matrix(mapply(rnorm，n，0，rep(1，nNorm)))

}

if(nBin> 0){

b<-as.matrix(mapply(rbinom，n，1，para$B))

}

if(nPos> 0){

c<-as.matrix(mapply(rpois，n，para$P))

}

if(nUni> 0){

d<-as.matrix(mapply(runif，n，para$U[，1]，para$U[，2]))

}

x <-cbind(a，b，c，d)

# 構(gòu)建lambda.i

beta.0<--alpha*log(lambda)

lambda.i<-exp(-(beta.0+x%*%beta)/alpha)

max.lambda.i<-max(lambda.i)

min.lambda.i<-min(lambda.i)

PdfLambdai<-function(u){

# 得到并返回lambda.i的gauss核密度估計

dens<-density(lambda.i，bw=“nrd0”，kernel=“gaussian”，na.rm=TRUE)

y.loess<-loess(dens$y～dens$x，span=0.1)

pred.y<-predict(y.loess，newdata=u)

return(pred.y)

}

# 據(jù)公式(13)得到條件刪失概率CondCensProb<-function(u)(u/(alpha*theta))*pgamma((theta/u)^alpha，1/alpha)*gamma(1/alpha)

# 得到刪失概率的表達(dá)式

cens.Prob<-integrate(function(u){PdfLambdai(u)*CondCensProb(u)}，min.lambda.i，max.lambda.i)$value

return(cens.Prob)

}

例4 以四個協(xié)變量分別來自N(0，1)，B(0.5)，Poisson(5)，Uniform(0，1)為例，給定alpha=2，lambda=4，指定的刪失比例cens.p=0.3，theta.mixed<- round(uniroot(function(y)CensProbMixed(class=c(“N”，“B”，“P”，“U”)，para=list(B=c(0.5)，P=c(5)，U=matrix(c(0，1)，ncol=2，byrow=TRUE))，y)-cens.p，c(0.001，100))$root，3)，得θ=22.698，即為了達(dá)到30%的刪失率，刪失時間的分布uniform(0，θ)的參數(shù)應(yīng)設(shè)為22.698，同理，模擬1000次得到的平均刪失率為0.29706。

討 論

為了通過模擬試驗對統(tǒng)計分析方法的表現(xiàn)性能進(jìn)行評價[6，13]，指定刪失比例的生存數(shù)據(jù)在研究中通常被廣泛應(yīng)用。

肖媛媛等人提出的方法基于完全隨機刪失的假設(shè)，其使得整體的刪失狀態(tài)服從二項分布，而每個個體的刪失時間則服從不同的均勻分布。該方法中刪失時間的分布依賴于事件發(fā)生時間，從而使得其產(chǎn)生的模擬數(shù)據(jù)變異性更大，本研究的方法相比之下顯得更為穩(wěn)定。

錢俊等人提出的大量取值通過模擬調(diào)試確定最終參數(shù)的方法雖然簡單，但如取值范圍選擇不當(dāng)將消耗大量的計算機資源和時間[10]?；贔ei Wan提出的框架[11]，本文介紹的方法可通過閉式表達(dá)式或者核密度估計方法得到刪失時間的分布參數(shù)從而生成滿足需求的生存數(shù)據(jù)，極大地提高了計算效率。同時，我們基于最基本的組合情形編寫了完整可用的生成刪失數(shù)據(jù)的函數(shù)：在事件發(fā)生時間和刪失時間分別為Weibull分布和均勻分布的條件下，當(dāng)無需考慮協(xié)變量時，可通過CenDatNocov()快速生成相關(guān)數(shù)據(jù)；當(dāng)考慮協(xié)變量對時間的影響時，針對單純二分類變量的CenDatBin()，針對正態(tài)分布數(shù)據(jù)的CenDatNorm()以及針對混合型協(xié)變量CenDatMixed()均可返回刪失參數(shù)及具有指定刪失比例的包含協(xié)變量的生存數(shù)據(jù)。相關(guān)的參考代碼可通過https：//github.com/zihang1012/simulated-survival-data-with-predefiend-censoring-rate-.git獲得。但是，由于本方法是在尋找表達(dá)式的解，因此會消耗一定的時間，此外，如果參數(shù)設(shè)置較為極端，可能會超過預(yù)先設(shè)定的尋根范圍，需要使用者額外調(diào)整函數(shù)中預(yù)先設(shè)定的界值。

結(jié)合各自需求，讀者可基于本研究進(jìn)一步拓展，如可通過指定基線中位生存時間，反推事件發(fā)生時間的分布參數(shù)得到滿足中位生存時間要求的基線生存數(shù)據(jù)；也可通過引入分組變量并指定回歸系數(shù)取值以模擬具有不同療效的RCT數(shù)據(jù)等。

本研究建立了一套科學(xué)可操作的模擬方法，給出了生存數(shù)據(jù)模擬研究中關(guān)于指定刪失比例的解決方案，具有實際應(yīng)用價值。本研究的局限性在于指定刪失時間分布為均勻分布，其他分布類型的刪失時間值得進(jìn)一步研究。