Cómo crear un juguete de supervivencia (tiempo para el evento) con los datos de derecho de censurar

Question

Quiero crear un juguete de supervivencia (tiempo para el evento) de datos que está a la derecha censurado y sigue alguna distribución de riesgos proporcionales y la constante de referencia de peligro.

He creado los datos de la siguiente manera, pero soy incapaz de obtener calcula los cocientes de riesgo que están cerca de los verdaderos valores después de la colocación de riesgos proporcionales de Cox modelo a los datos simulados.

¿Qué hice mal?

R códigos:

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

pero los verdaderos valores se establecen como

B = c(1.1,1.2,1.3) # hazard ratios (model coefficients)

Answer 1

No está claro para mí cómo generar tu evento de veces (que podría ser $<0$) y el evento indicadores:

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

Así que aquí es un método genérico, seguido por algunos R código.

---

La generación de los tiempos de supervivencia para simular modelos de riesgos proporcionales de Cox

Para generar las horas de los eventos desde el modelo de riesgos proporcionales, podemos utilizar la inversa de la probabilidad de método (Bender et al., 2005): si $V$ es uniforme en $(0, 1)$ e si $S(\cdot \,|\, \mathbf{x})$ es el condicional de la supervivencia de la función derivada a partir de el modelo de riesgos proporcionales, es decir, $$ S(t \,|\, \mathbf{x}) = \exp \left( -H_0(t) \exp(\mathbf{x}^\prime \mathbf{\beta}) \vphantom{\Big(} \right) $$ entonces es un hecho de que la variable aleatoria $$ T = S^{-1}(V \,|\, \mathbf{x}) = H_0^{-1} \left( - \frac{\log(V)}{\exp(\mathbf{x}^\prime \mathbf{\beta})} \right) $$ tiene la función de sobrevivencia $S(\cdot \,|\, \mathbf{x})$. Este resultado es conocido como `la inversa de la probabilidad de transformación integral". Por lo tanto, para generar un tiempo de supervivencia $T \sim S(\cdot \,|\, \mathbf{x})$ dado el vector covariante, es suficiente para dibujar $v$$V \sim \mathrm{U}(0, 1)$, y para hacer la transformación inversa de a $t = S^{-1}(v \,|\, \mathbf{x})$.

---

Ejemplo [de Weibull de referencia de peligro]

Deje $h_0(t) = \lambda \rho t^{\rho - 1}$ forma $\rho > 0$ y la escala de la $\lambda > 0$. A continuación,$H_0(t) = \lambda t^\rho$$H^{-1}_0(t) = (\frac{t}{\lambda})^{\frac{1}{\rho}}$. Tras la probabilidad inversa método, una realización de $T \sim S(\cdot \,|\, \mathbf{x})$ se obtiene mediante el cálculo de $$ t = \left( - \frac{\log(v)}{\lambda \exp(\mathbf{x}^\prime \mathbf{\beta})} \right)^{\frac{1}{\rho}} $$ con $v$ un uniforme de la variable aleatoria en $(0, 1)$. El uso de los resultados en las transformaciones de variables aleatorias, uno puede notar que $T$ tiene un condicional de la distribución de Weibull (dado $\mathbf{x}$) con forma de $\rho$ y la escala de la $\lambda \exp(\mathbf{x}^\prime \mathbf{\beta})$.

---

R código de

La siguiente función de R genera un conjunto de datos con un único binario covariable $x$ (por ejemplo, un indicador de tratamiento). La línea de base peligro tiene forma de Weibull. Censurar a veces son seleccionados al azar de una distribución exponencial.

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

---

Prueba

Aquí está una rápida simulación con $\beta = -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