Je souhaite créer une donnée de survie du jouet (temps avant l'événement) qui est censurée à droite et suit une certaine distribution avec des dangers proportionnels et un danger de base constant.
J'ai créé les données comme suit, mais je ne suis pas en mesure d'obtenir des ratios de risque estimés qui sont proches des valeurs réelles après avoir ajusté un modèle de risques proportionnels de Cox aux données simulées.
Qu'ai-je fait de mal?
Codes 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)
Résultats:
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
mais les vraies valeurs sont définies comme
B = c(1.1,1.2,1.3) # hazard ratios (model coefficients)
survival
cox-model
monte-carlo
stats_newb
la source
la source
Réponses:
Pour moi, la façon dont vous générez vos heures d'événement (qui, dans votre cas, peut être ) et vos indicateurs d'événement ne sont pas clairs :<0
Voici donc une méthode générique, suivie d'un code R.
Génération de temps de survie pour simuler des modèles de risques proportionnels de Cox
Pour générer des temps d'événements à partir du modèle des risques proportionnels, nous pouvons utiliser la méthode de probabilité inverse (Bender et al., 2005) : si est uniforme sur et si est la fonction de survie conditionnelle dérivée du modèle des risques proportionnels, c'est-à-dire alors c'est un fait que la variable aléatoire a la fonction de survieV (0,1) S(⋅|x)
Exemple [danger de base de Weibull]
Soit de forme et d'échelle . Alors et . En suivant la méthode de probabilité inverse, une réalisation de est obtenue en calculant avec une variable uniforme sur . En utilisant des résultats sur les transformations de variables aléatoires, on peut remarquer que a une distribution de Weibull conditionnelle (étant donnéh0(t)=λρtρ−1 ρ>0 λ>0 H0(t)=λtρ H−10(t)=(tλ)1ρ T∼S(⋅|x) v(0,1)Txρλexp(x′β)
Code R
La fonction R suivante génère un ensemble de données avec une seule covariable binaire (par exemple un indicateur de traitement). Le danger de base a une forme de Weibull. Les temps de censure sont tirés au hasard à partir d'une distribution exponentielle.x
Tester
Voici une simulation rapide avec :β=−0.6
la source
flexsurvreg(Surv(time, status) ~ x, data=dat, dist = "weibull")
sur les mêmes données simulées, le coefficient apparaît comme0.6212
. Pourquoi est-ce?Pour la distribution de Weibull,e−(λ∗e(x∗β)∗t)ρ
S (t) =
" " sera uniquement pour le journal (v)(1/rho)
donc j'ai modifié comme ça
si rho = 1, le résultat sera le même.
la source