Cara membuat data kelangsungan hidup mainan (waktu ke acara) dengan sensor yang benar

Saya ingin membuat data survival mainan (waktu untuk acara) yang disensor dengan benar dan mengikuti distribusi dengan bahaya proporsional dan bahaya baseline konstan.

Saya membuat data sebagai berikut, tetapi saya tidak dapat memperoleh estimasi rasio bahaya yang mendekati nilai sebenarnya setelah menyesuaikan model bahaya proporsional Cox dengan data yang disimulasikan.

Apa kesalahan yang telah aku perbuat?

Kode 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)

Hasil:

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

tetapi nilai true ditetapkan sebagai

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

Tidak jelas bagi saya bagaimana Anda menghasilkan waktu acara Anda (yang, dalam kasus Anda, mungkin ) dan indikator acara:$<0$

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

Jadi di sini adalah metode generik, diikuti oleh beberapa kode R.

---

Menghasilkan waktu bertahan hidup untuk mensimulasikan model bahaya proporsional Cox

Untuk menghasilkan waktu kejadian dari model bahaya proporsional, kita dapat menggunakan metode probabilitas terbalik (Bender et al., 2005) : jika seragam pada dan jika adalah fungsi survival bersyarat yang berasal dari model bahaya proporsional, yaitu maka itu adalah fakta bahwa variabel acak memiliki fungsi survival$V$$(0, 1)$$S(\cdot \,|\, \mathbf{x})$

$$S(t \,|\, \mathbf{x}) = \exp \left( -H_0(t) \exp(\mathbf{x}^\prime \mathbf{\beta}) \vphantom{\Big(} \right)$$ $$T = S^{-1}(V \,|\, \mathbf{x}) = H_0^{-1} \left( - \frac{\log(V)}{\exp(\mathbf{x}^\prime \mathbf{\beta})} \right)$$ $S(\cdot \,|\, \mathbf{x})$. Hasil ini dikenal sebagai `` transformasi integral probabilitas terbalik ''. Oleh karena itu, untuk menghasilkan waktu survival dengan vektor kovariat, cukup untuk menggambar dari dan untuk membuat transformasi terbalik .$T \sim S(\cdot \,|\, \mathbf{x})$$v$$V \sim \mathrm{U}(0, 1)$$t = S^{-1}(v \,|\, \mathbf{x})$

---

Contoh [Bahaya garis dasar Weibull]

Biarkan dengan bentuk dan skala . Kemudian dan . Mengikuti metode probabilitas terbalik, realisasi diperoleh dengan menghitung dengan variasi yang seragam pada . Dengan menggunakan hasil pada transformasi variabel acak, orang mungkin memperhatikan bahwa memiliki distribusi Weibull bersyarat (diberikan$h_0(t) = \lambda \rho t^{\rho - 1}$$\rho > 0$$\lambda > 0$$H_0(t) = \lambda t^\rho$ T∼S($H^{-1}_0(t) = (\frac{t}{\lambda})^{\frac{1}{\rho}}$t = ( - log ( v $T \sim S(\cdot \,|\, \mathbf{x})$ v(0,1)Txρλexp(x′β

$$t = \left( - \frac{\log(v)}{\lambda \exp(\mathbf{x}^\prime \mathbf{\beta})} \right)^{\frac{1}{\rho}}$$ $v$$(0, 1)$$T$$\mathbf{x}$) dengan bentuk dan skala .$\rho$$\lambda \exp(\mathbf{x}^\prime \mathbf{\beta})$

---

Kode r

Fungsi R berikut menghasilkan set data dengan kovariat biner tunggal (misalnya indikator perawatan). Bahaya garis dasar memiliki bentuk Weibull. Waktu sensor diambil secara acak dari distribusi eksponensial.$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)
}

---

Uji

Berikut ini beberapa simulasi cepat dengan :$\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

Untuk distribusi Weibull,
S (t) =$e^{-(\lambda * e^(x * \beta)*t)^\rho}$

" " hanya untuk log (v)$^{(1/rho)}$

jadi, saya modifikasi seperti ini

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

jika rho = 1, hasilnya akan sama.