威布尔持续时间模型的模拟
simulation of a weibull duration model
我想创建一个遵循 Weibull 分布的事件变量。重要的是变量应该是其他几个观察到的变量的组合。
例如:死亡是事件变量遵循我想模拟的威布尔分布的时间(这里我的时间尺度是年龄)。我已经有(模拟)变量,例如年龄、性别、BMI 和癌症的 4 个阶段(具有 4 个类别的分类变量),因此我想使用这 4 个变量来模拟事件死亡时间变量。
让我知道是否需要澄清
如果我没记错的话,您对加速失效时间 (AFT) Weibull 模型感兴趣。
生存函数为:
S(t) = exp(- lambda t^p)
其中 lambda 和 p 是尺度和形状参数。 objective 是参数化 lambda。如果求解 t,并假设固定概率 S(t) = q 将得到
t = A*B
其中 A = (− log(q))^ 1/p 和 B = (1/lambda)^(1/p)
对于二元处理指标 TREAT,参数化 lambda:B = exp(beta_0 + beta_1*TREAT)。加速因子是 exp(beta_1)(您可以通过处理变量相对于控制变量的 A*B 表达式的比率来了解这一点)。
您可以按照上面的 AB 表达式模拟您的数据,注意系数、随机分量和固定概率分量。特别是,如果您使用正态分布,极值可能导致负时间,这是没有意义的。时间必须是非负数。
set.seed(123)
library(data.table)
library(survival)
# generate data
# (can use base r or dplyr if not familiar with data.table)
n <- 2000
d <- data.table(id=1:n,
age = runif(n,40,80),
male = rbinom(n,1,0.5),
bmi = runif(n,15,30),
cancer = sample(letters[1:4], n, replace = T), # cancer stages
e = runif(n, 0,2) ) # some error, uniform for instance
# you will need to transform the cancer variable into numeric,
# one category will be the comparison group
d[, cancer_a := ifelse(cancer=="a", 1, 0)]
d[, cancer_b := ifelse(cancer=="b", 1, 0)]
d[, cancer_c := ifelse(cancer=="c", 1, 0)]
# add S(t)
shape <- 1
d[, s_tcomp := (-log(0.01))^(1/shape) ]
# generate the time
d[, time := s_tcomp*exp( -0.001*age - 0.1*male + 0.1*bmi + 0.3*cancer_a + 0.2*cancer_b + 0.1*cancer_c + e)]
#' In case you want to add censoring:
#' we measure time only up to a certain period,
#' if didnt die so far then still alive
censor <- quantile(d[,time], 0.9)
d[, dead := ifelse(time<censor, 1, 0) ]
d[, time := pmin(time, censor) ]
m <- survreg( Surv(time, dead) ~ age + male + bmi + cancer_a + cancer_b + cancer_c,
data=d, dist = "weibull", )
summary(m)
Call:
survreg(formula = Surv(time, dead) ~ age + male + bmi + cancer_a +
cancer_b + cancer_c, data = d, dist = "weibull")
Value Std. Error z p
(Intercept) 2.791584 0.098517 28.34 < 2e-16
age -0.000943 0.001091 -0.86 0.3874
male -0.058586 0.024720 -2.37 0.0178
bmi 0.099430 0.003071 32.37 < 2e-16
cancer_a 0.297261 0.034977 8.50 < 2e-16
cancer_b 0.177142 0.034474 5.14 2.8e-07
cancer_c 0.101467 0.034039 2.98 0.0029
Log(scale) -0.650129 0.018555 -35.04 < 2e-16
Scale= 0.522
Weibull distribution
Loglik(model)= -10272.5 Loglik(intercept only)= -10751
Chisq= 956.93 on 6 degrees of freedom, p= 1.8e-203
Number of Newton-Raphson Iterations: 5
n= 2000
另请参阅:
https://cran.r-project.org/web/packages/coxed/vignettes/simulating_survival_data.html
https://cran.r-project.org/web/packages/simsurv/vignettes/simsurv_usage.html
我想创建一个遵循 Weibull 分布的事件变量。重要的是变量应该是其他几个观察到的变量的组合。
例如:死亡是事件变量遵循我想模拟的威布尔分布的时间(这里我的时间尺度是年龄)。我已经有(模拟)变量,例如年龄、性别、BMI 和癌症的 4 个阶段(具有 4 个类别的分类变量),因此我想使用这 4 个变量来模拟事件死亡时间变量。
让我知道是否需要澄清
如果我没记错的话,您对加速失效时间 (AFT) Weibull 模型感兴趣。
生存函数为:
S(t) = exp(- lambda t^p)
其中 lambda 和 p 是尺度和形状参数。 objective 是参数化 lambda。如果求解 t,并假设固定概率 S(t) = q 将得到
t = A*B
其中 A = (− log(q))^ 1/p 和 B = (1/lambda)^(1/p)
对于二元处理指标 TREAT,参数化 lambda:B = exp(beta_0 + beta_1*TREAT)。加速因子是 exp(beta_1)(您可以通过处理变量相对于控制变量的 A*B 表达式的比率来了解这一点)。
您可以按照上面的 AB 表达式模拟您的数据,注意系数、随机分量和固定概率分量。特别是,如果您使用正态分布,极值可能导致负时间,这是没有意义的。时间必须是非负数。
set.seed(123)
library(data.table)
library(survival)
# generate data
# (can use base r or dplyr if not familiar with data.table)
n <- 2000
d <- data.table(id=1:n,
age = runif(n,40,80),
male = rbinom(n,1,0.5),
bmi = runif(n,15,30),
cancer = sample(letters[1:4], n, replace = T), # cancer stages
e = runif(n, 0,2) ) # some error, uniform for instance
# you will need to transform the cancer variable into numeric,
# one category will be the comparison group
d[, cancer_a := ifelse(cancer=="a", 1, 0)]
d[, cancer_b := ifelse(cancer=="b", 1, 0)]
d[, cancer_c := ifelse(cancer=="c", 1, 0)]
# add S(t)
shape <- 1
d[, s_tcomp := (-log(0.01))^(1/shape) ]
# generate the time
d[, time := s_tcomp*exp( -0.001*age - 0.1*male + 0.1*bmi + 0.3*cancer_a + 0.2*cancer_b + 0.1*cancer_c + e)]
#' In case you want to add censoring:
#' we measure time only up to a certain period,
#' if didnt die so far then still alive
censor <- quantile(d[,time], 0.9)
d[, dead := ifelse(time<censor, 1, 0) ]
d[, time := pmin(time, censor) ]
m <- survreg( Surv(time, dead) ~ age + male + bmi + cancer_a + cancer_b + cancer_c,
data=d, dist = "weibull", )
summary(m)
Call:
survreg(formula = Surv(time, dead) ~ age + male + bmi + cancer_a +
cancer_b + cancer_c, data = d, dist = "weibull")
Value Std. Error z p
(Intercept) 2.791584 0.098517 28.34 < 2e-16
age -0.000943 0.001091 -0.86 0.3874
male -0.058586 0.024720 -2.37 0.0178
bmi 0.099430 0.003071 32.37 < 2e-16
cancer_a 0.297261 0.034977 8.50 < 2e-16
cancer_b 0.177142 0.034474 5.14 2.8e-07
cancer_c 0.101467 0.034039 2.98 0.0029
Log(scale) -0.650129 0.018555 -35.04 < 2e-16
Scale= 0.522
Weibull distribution
Loglik(model)= -10272.5 Loglik(intercept only)= -10751
Chisq= 956.93 on 6 degrees of freedom, p= 1.8e-203
Number of Newton-Raphson Iterations: 5
n= 2000
另请参阅:
https://cran.r-project.org/web/packages/coxed/vignettes/simulating_survival_data.html
https://cran.r-project.org/web/packages/simsurv/vignettes/simsurv_usage.html