Я хочу создать данные о выживаемости игрушек (время до события), которые подвергаются цензуре и следуют некоторому распределению с пропорциональными опасностями и постоянной базовой опасностью.
Я создал данные следующим образом, но я не могу получить расчетные коэффициенты опасности, близкие к истинным значениям, после подбора модели пропорциональных рисков Кокса для смоделированных данных.
Что я сделал не так?
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)
survival
cox-model
monte-carlo
stats_newb
источник
источник
Ответы:
Мне не ясно, как вы генерируете время вашего события (которое в вашем случае может быть ) и индикаторы события:< 0
Итак, вот общий метод, за которым следует R-код.
Генерация времени выживания для имитации моделей пропорциональных рисков Кокса
Чтобы сгенерировать времена событий из модели пропорциональных опасностей, мы можем использовать метод обратной вероятности (Bender et al., 2005) : если равномерно по и если - функция условного выживания, полученная из модели пропорциональных рисков, т.е. тогда это факт, что случайная величина имеет функцию выживания( 0 , 1 ) S ( ⋅В ( 0 , 1 ) S ( тS( ⋅|х )
Пример [базовая опасность Вейбулла]
Пусть с формой и масштабом . Тогда и . Следуя методу обратной вероятности, реализация получается вычислением с равномерной переменной на . Используя результаты о преобразованиях случайных величин, можно заметить, что имеет условное распределение Вейбулла (учитываяh0(t)=λρtρ−1 ρ>0 λ>0 H0(t)=λtρ H−10(t)=(tλ)1ρ t = ( - log ( v )T∼S(⋅|x) v(0,1)Txρλexp(x′β)
Код R
Следующая функция R генерирует набор данных с одним двоичным ковариатом (например, индикатором лечения). Базовая опасность имеет форму Вейбулла. Время цензуры выбирается случайным образом из экспоненциального распределения.x
Тест
Вот небольшая симуляция с :β=−0.6
источник
flexsurvreg(Surv(time, status) ~ x, data=dat, dist = "weibull")
те же данные, смоделированные, коэффициент отображается как0.6212
. Почему это?Для распределения Вейбуллаe−(λ∗e(x∗β)∗t)ρ
S (t) =
« » будет только для log (v)(1/rho)
Итак, я изменил, как это
если rho = 1, результат будет таким же.
источник