如何創建具有正確審查的玩俱生存（事件時間）數據

我希望創建一個玩俱生存（事件發生時間）數據，該數據經過正確審查，並遵循一些具有比例風險和恆定基線風險的分佈。

我按如下方式創建了數據，但在將 Cox 比例風險模型擬合到模擬數據後，我無法獲得接近真實值的估計風險比。

我做錯了什麼？

R代碼：

library(survival)

#set parameters
set.seed(1234)

n = 40000 #sample size


#functional relationship

lambda=0.000020 #constant baseline hazard 2 per 100000 per 1 unit time

b_haz <-function(t) #baseline hazard
 {
   lambda #constant hazard wrt time 
 }

x = cbind(hba1c=rnorm(n,2,.5)-2,age=rnorm(n,40,5)-40,duration=rnorm(n,10,2)-10)

B = c(1.1,1.2,1.3) # hazard ratios (model coefficients)

hist(x %*% B) #distribution of scores

haz <-function(t) #hazard function
{
 b_haz(t) * exp(x %*% B)
}

c_hf <-function(t) #cumulative hazards function
{
 exp(x %*% B) * lambda * t 
}

S <- function(t) #survival function
{
 exp(-c_hf(t))
}

S(.005)
S(1)
S(5)

#simulate censoring

time = rnorm(n,10,2)

S_prob = S(time)

#simulate events

event = ifelse(runif(1)>S_prob,1,0)

#model fit

km = survfit(Surv(time,event)~1,data=data.frame(x))

plot(km) #kaplan-meier plot

#Cox PH model

fit = coxph(Surv(time,event)~ hba1c+age+duration, data=data.frame(x))

summary(fit)            

cox.zph(fit)

結果：

Call:
coxph(formula = Surv(time, event) ~ hba1c + age + duration, data = data.frame(x))

 n= 40000, number of events= 3043 

            coef exp(coef) se(coef)     z Pr(>|z|)    
hba1c    0.236479  1.266780 0.035612  6.64 3.13e-11 ***
age      0.351304  1.420919 0.003792 92.63  < 2e-16 ***
duration 0.356629  1.428506 0.008952 39.84  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
hba1c        1.267     0.7894     1.181     1.358
age          1.421     0.7038     1.410     1.432
duration     1.429     0.7000     1.404     1.454

Concordance= 0.964  (se = 0.006 )
Rsquare= 0.239   (max possible= 0.767 )
Likelihood ratio test= 10926  on 3 df,   p=0
Wald test            = 10568  on 3 df,   p=0
Score (logrank) test = 11041  on 3 df,   p=0

但真實值設置為

B = c(1.1,1.2,1.3) # hazard ratios (model coefficients)

我不清楚您如何生成事件時間（在您的情況下，可能是) 和事件指標：

time = rnorm(n,10,2) 
S_prob = S(time)
event = ifelse(runif(1)>S_prob,1,0)

所以這裡是一個通用方法，後面是一些 R 代碼。

---

生成生存時間以模擬 Cox 比例風險模型

為了從比例風險模型生成事件時間，我們可以使用逆概率法（Bender 等人，2005）：如果是統一的而如果是從比例風險模型導出的條件生存函數，即

那麼事實是隨機變量

有生存功能. 這個結果被稱為“逆概率積分變換”。因此，要生成生存時間給定協變量向量，繪製就足夠了從並進行逆變換.

---

示例 [Weibull 基線危害]

讓有形狀和規模. 然後和. 遵循逆概率方法，實現通過計算得到

和一個統一的變量. 使用隨機變量變換的結果，人們可能會注意到具有條件 Weibull 分佈（給定) 與形狀和規模.

---

R代碼

以下 R 函數生成具有單個二進制協變量的數據集（例如治療指標）。基線危險具有 Weibull 形式。審查時間是從指數分佈中隨機抽取的。

# baseline hazard: Weibull

# N = sample size 
# lambda = scale parameter in h0()
# rho = shape parameter in h0()
# beta = fixed effect parameter
# rateC = rate parameter of the exponential distribution of C

simulWeib <- function(N, lambda, rho, beta, rateC)
{
 # covariate --> N Bernoulli trials
 x <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5))

 # Weibull latent event times
 v <- runif(n=N)
 Tlat <- (- log(v) / (lambda * exp(x * beta)))^(1 / rho)

 # censoring times
 C <- rexp(n=N, rate=rateC)

 # follow-up times and event indicators
 time <- pmin(Tlat, C)
 status <- as.numeric(Tlat <= C)

 # data set
 data.frame(id=1:N,
            time=time,
            status=status,
            x=x)
}

---

測試

這是一些快速模擬:

set.seed(1234)
betaHat <- rep(NA, 1e3)
for(k in 1:1e3)
{
 dat <- simulWeib(N=100, lambda=0.01, rho=1, beta=-0.6, rateC=0.001)
 fit <- coxph(Surv(time, status) ~ x, data=dat)
 betaHat[k] <- fit$coef
}

> mean(betaHat)
[1] -0.6085473

引用自：https://stats.stackexchange.com/questions/135124