Comment créer des données de survie de jouet (temps avant l'événement) avec une censure à droite

12

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)
stats_newb
la source
1
pour votre tâche, un démarrage rapide consiste à utiliser un package de simulation existant: cran.r-project.org/web/packages/survsim/index.html
zhanxw

Réponses:

19

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

time = rnorm(n,10,2) 
S_prob = S(time)
event = ifelse(runif(1)>S_prob,1,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)

S(t|x)=exp(H0(t)exp(xβ)()
T=S1(V|x)=H01(log(V)exp(xβ))
S(|x). Ce résultat est connu comme `` la transformation intégrale de probabilité inverse ''. Par conséquent, pour générer un temps de survie étant donné le vecteur covariable, il suffit de tirer de et pour faire la transformation inverse .TS(|x)vVU(0,1)t=S1(v|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λ>0H0(t)=λtρH01(t)=(tλ)1ρTS(|x) v(0,1)Txρλexp(xβ)

t=(log(v)λexp(xβ))1ρ
v(0,1)Tx) avec la forme et l'échelle .ρλ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

# 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)
}

Tester

Voici une simulation rapide avec :β=0.6

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
ocram
la source
Merci pour votre excellente réponse. J'ai réalisé que j'avais gâché les heures des événements en obtenant le statut des événements après avoir aléatoire les heures des événements, ce qui n'avait pas de sens .. idiot!
stats_newb
Puis-je vous demander s'il existe une raison particulière pour laquelle vous tirez le temps de censure d'une distribution exponentielle?
pthao
@pthao: il n'y a pas de raison particulière (c'était juste une illustration où j'ai utilisé la distribution exponentielle)
ocram
1
Existe-t-il des directives pour choisir la distribution des temps de censure?
pthao
@ocram Fait intéressant, lorsque je cours flexsurvreg(Surv(time, status) ~ x, data=dat, dist = "weibull")sur les mêmes données simulées, le coefficient apparaît comme 0.6212. Pourquoi est-ce?
ni-ni
3

Pour la distribution de Weibull,
S (t) =e(λe(xβ)t)ρ

" " sera uniquement pour le journal (v)(1/rho)

donc j'ai modifié comme ça

Tlat <- (- log(v))^(1 / rho) / (lambda * exp(x * beta))

si rho = 1, le résultat sera le même.

unko
la source