指定刪失比例的生存數據模擬及R實現*

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

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

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

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

方法與原理

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

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

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

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

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

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

Weibull分布的概率密度函數和風險函數為：

(1)

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

(2)

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

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

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

(3)

各個體的刪失概率為：

(4)

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

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

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

(5)

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

(6)

(7)

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

(8)

(9)

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

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

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

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

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

(10)

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

(11)

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

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

(12)

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

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

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

(13)

z=(c/λi)α

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

給定回歸系數β=<β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實現與應用

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

我們構建CenDatNoCov()函數來產生指定刪失比例的生存數據，輸入參數包括Weibull(α，λ)的形狀alpha和尺度參數lambda，指定的刪失比例cens.p及樣本量size。由于不考慮協(xié)變量，我們僅需通過uniroot()函數對式(6)進行求根即可得到刪失參數theta，其中下不完全Gamma函數γ(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)

#據方法與原理的描述產生刪失的生存數據

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

2.考慮協(xié)變量

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

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

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

alpha=alpha # weibull分布中的形狀參數

lambda=lambda # weibull分布中的形狀參數

cens.p=cens.p#預設的刪失比例

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

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

n=size#產生數據的樣本量

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

CensProbBin <- function(theta){

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

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

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

CombInd <-list(NULL)

ncomb <- c()

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

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

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

CombInd[[i]]<- ind

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

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

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

}

}

# 根據各協(xié)變量的下標組合進行賦值，得到最終的取值組合

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

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

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

}

# 據公式(9)構建單變量，簡化計算

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

#計算在對應協(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)

}

#求解刪失參數值

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()函數中scale參數值為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)))

#據方法與原理的描述產生刪失的生存數據

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，即為了達到30%的刪失比例，刪失時間的分布uniform(0，θ)的參數應設為11.116?？赏ㄟ^上述方法生成1000次模擬數據來驗證結果，當種子數設定為20200810時，得到的刪失比例平均值為0.30132，滿足預設的要求。

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

我們構建CenDatNorm()來產生均為正態(tài)分布協(xié)變量的刪失數據，輸入參數基本同上。內嵌計算刪失參數的相關函數為CensProbNorm()，通過求根即可獲得theta值：

CensProbNorm<-function(theta){

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

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

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

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

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

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

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

return(cens.Prob)

}

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

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

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

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

set.seed(20200810)

#高斯核密度估計lambda.i的密度函數

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

# 生成構成lambda.i的相關變量值以提供擬合聯(lián)合分布的數據

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)

# 構建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)

}

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

# 得到刪失概率的表達式

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，即為了達到30%的刪失率，刪失時間的分布uniform(0，θ)的參數應設為22.698，同理，模擬1000次得到的平均刪失率為0.29706。

討 論

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

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

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

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

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