Como criar dados de sobrevivência de um brinquedo (hora do evento) com a censura correta

12

Desejo criar dados de sobrevivência de brinquedos (tempo até o evento) que sejam corretamente censurados e sigam alguma distribuição com riscos proporcionais e riscos constantes da linha de base.

Criei os dados da seguinte maneira, mas não consigo obter taxas de risco estimadas próximas dos valores reais depois de ajustar um modelo de riscos proporcionais de Cox aos dados simulados.

O que eu fiz errado?

Códigos 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)

Resultados:

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

mas os valores verdadeiros são definidos como

B = c(1.1,1.2,1.3) # hazard ratios (model coefficients)
stats_newb
fonte
1
para a sua tarefa, um início rápido é usar um pacote de simulação existente: cran.r-project.org/web/packages/survsim/index.html
zhanxw

Respostas:

19

Não está claro para mim como você gera seus horários de evento (que, no seu caso, podem ser ) e indicadores de evento:<0

time = rnorm(n,10,2) 
S_prob = S(time)
event = ifelse(runif(1)>S_prob,1,0)

Então, aqui está um método genérico, seguido por algum código R.


Gerando tempos de sobrevivência para simular modelos de riscos proporcionais de Cox

Para gerar tempos de eventos a partir do modelo de riscos proporcionais, podemos usar o método de probabilidade inversa (Bender et al., 2005) : se for uniforme em e se é a função de sobrevivência condicional derivada do modelo de riscos proporcionais, ou seja, então é fato que a variável aleatória possui a função de sobrevivênciaV(0,1)S(|x)

S(t|x)=exp(H0(t)exp(xβ)()
T=S1(V|x)=H01(log(V)exp(xβ))
S(|x). Esse resultado é conhecido como `` a transformação integral de probabilidade inversa ''. Portanto, para gerar um tempo de sobrevivência dado o vetor covariável, basta desenhar de e para fazer a transformação inversa .TS(|x)vVU(0,1)t=S1(v|x)

Exemplo [risco de linha de base Weibull]

Deixe com a forma e a escala . Então e . Seguindo o método de probabilidade inversa, uma realização de é obtida calculando com uma variável uniforme em . Usando resultados em transformações de variáveis ​​aleatórias, pode-se notar que tem uma distribuição Weibull condicional (dadoh0(t)=λρtρ1ρ>0λ>0H0(t)=λtρ TS(H01(t)=(tλ)1ρt = ( - log ( v )TS(|x) v(0,1)Txρλexp(xβ)

t=(log(v)λexp(xβ))1ρ
v(0,1)Tx) com forma e escala .ρλexp(xβ)

Código R

A função R a seguir gera um conjunto de dados com uma única covariável binária (por exemplo, um indicador de tratamento). O risco da linha de base tem uma forma Weibull. Os tempos de censura são sorteados aleatoriamente a partir de uma distribuição exponencial.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)
}

Teste

β=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
fonte
Obrigado pela sua excelente resposta. Percebi que tinha atrapalhado os horários dos eventos, obtendo o status dos eventos depois que aleatorizei os horários dos eventos, o que não fazia sentido ... parvo!
stats_newb 27/01
Posso perguntar se existe algum motivo específico para você extrair o tempo de censura de uma distribuição exponencial?
pthao 16/05
@pthao: não há nenhuma razão particular (este foi apenas uma ilustração que eu usei a distribuição exponencial)
Ocram
1
Existe alguma orientação para escolher a distribuição para os tempos de censura?
Pthao 16/05
Curiosamente, quando executo flexsurvreg(Surv(time, status) ~ x, data=dat, dist = "weibull")os mesmos dados simulados, o coeficiente aparece como 0.6212. Por que é isso?
nem-
3


e(λe(xβ)t)ρ

(1/rho)

então eu modifiquei assim

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

se rho = 1, o resultado será o mesmo.

unko
fonte