имитация модели длительности Вейбулла - PullRequest
Новая
       114

имитация модели длительности Вейбулла

0 голосов
/ 03 мая 2020

Я хочу создать переменную события, которая следует за распределением Вейбулла. Важно то, что переменная должна быть комбинацией нескольких других наблюдаемых переменных.

Например: Смерть - это время, когда переменная события следует распределению Вейбулла, которое я хочу смоделировать (здесь мой шкала времени - возраст). У меня уже есть (смоделированные) переменные, такие как возраст, пол, ИМТ и 4 стадии рака (категориальная переменная с 4 категориями), поэтому, используя эти 4 переменные, я хочу смоделировать переменную времени до события.

Пусть я знаю, есть ли необходимость в разъяснениях

1 Ответ

1 голос
/ 04 мая 2020

Если я не ошибаюсь, вас интересует модель Weibull с ускоренным временем отказа (AFT).

Функция выживания:

S (t) = exp (- лямбда-т ^ p)

с лямбдой и p как шкала и параметр формы. Цель состоит в том, чтобы параметризовать лямбду. Если вы решите для t, и предположите, что фиксированная вероятность S (t) = q получит

t = A * B

где A = (- log (q)) ^ 1 / p и B = (1 / лямбда) ^ (1 / p)

Для двоичного индикатора обработки TREAT параметризовать лямбду: 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

https://www.ms.uky.edu/~mai/Rsurv.pdf

...