Hitung tingkat insiden menggunakan model poisson: hubungan dengan rasio bahaya dari model Cox PH

Saya ingin menghitung tingkat kejadian untuk disajikan di sepanjang rasio bahaya untuk menyajikan ukuran risiko relatif dan absolut. Saya melihat dalam penelitian lain bahwa tingkat kejadian seperti itu dapat dihitung menggunakan model poisson dengan waktu tindak lanjut dalam model sebagai offset. Jadi saya mencobanya di R sebagai berikut:

library(survival)

# Get example data
data(colon)
colon$status <- ifelse(colon$etype==1,0,1) # set to 0/1 (needed for poisson later on)

# Fit cox model for rx (age + sex adjusted)
coxph(Surv(time,status)~rx+sex+age, data=colon)
# HR (rxLev): 0.92  
# HR (rxLev+5FU): 0.74

# Get incidence rates using poisson models with same terms and log(time) as offset
mod <- glm(status~offset(log(time))+rx+sex+age, data=colon, family=poisson)

# Get rates using predict-function
Obs <- predict(mod, data.frame(time=1, rx="Obs", age=mean(colon$age),
                                   sex=mean(colon$sex)),  type="response")
Lev <- predict(mod, data.frame(time=1, rx="Lev", age=mean(colon$age), 
                                   sex=mean(colon$sex)),  type="response")
Lev5FU <- predict(mod, data.frame(time=1, rx="Lev+5FU", age=mean(colon$age), 
                                      sex=mean(colon$sex)),  type="response")

# Calculate incidence rate ratio's:
Lev/Obs # 0.98
Lev5FU/Obs # 0.84

Saya berharap bahwa rasio tingkat kejadian mirip dengan rasio bahaya dari model Cox PH dengan istilah yang sama, tetapi entah bagaimana mereka berbeda. Apakah saya menggunakan pendekatan yang benar untuk menghitung tingkat kejadian?

Bantuan apa pun akan sangat dihargai!

Sejauh yang saya bisa lihat, tidak ada yang salah dengan kode atau perhitungan Anda. Anda dapat melewati beberapa baris kode, dengan mendapatkan rasio tingkat kejadian dengan . Kedua model membuat asumsi yang berbeda, dan ini berpotensi mengarah pada hasil yang berbeda.${\tt exp(coef(mod))}$

Regresi Poisson mengasumsikan bahaya konstan. Model Cox hanya mengasumsikan bahayanya proporsional. Jika asumsi bahaya konstan terpenuhi pertanyaan ini

Apakah Cox Regression memiliki distribusi Poisson yang mendasarinya?

menjelaskan hubungan antara regresi Cox dan Poisson.

Kita dapat menggunakan simulasi untuk mempelajari dua situasi: bahaya konstan dan bahaya tidak konstan (tetapi proporsional). Pertama mari kita mensimulasikan data dari populasi dengan bahaya konstan. Rasio bahaya memiliki formulir

$\lambda(t) = \lambda_0\exp(\beta'x)$ ,

di mana adalah vektor parameter, adalah vektor kovariat dan adalah bilangan positif tetap. Dengan demikian, asumsi bahaya konstan dari regresi Poisson terpenuhi. Sekarang kita mensimulasikan dari model ini dengan menggunakan fakta (ditemukan dalam banyak buku, misalnya Modeling Survival Data oleh Therneau, hal.13) bahwa fungsi distribusi, , dari waktu bertahan hidup dengan sebagai bahaya dapat ditemukan sebagai$\beta$$x$$\lambda_0$$F$$\lambda$

$F(t) = 1 - \exp\left(\int_0^t \lambda(s)\text{ d}s\right)$ .

Dengan ini kita juga dapat menemukan kebalikan dari , . Dengan fungsi ini kami mensimulasikan waktu bertahan hidup dengan bahaya yang benar dengan menggambar variabel yang seragam pada dan mengubahnya menggunakan . Ayo lakukan.$F$$F^{-1}$$(0,1)$$F^{-1}$

library(survival)
data(colon)
data <- with(colon, data.frame(sex = sex, rx = rx, age = age))
n <- dim(data)[1]
# defining linP, the linear predictor, beta*x in the above notation
linP <- with(colon, log(0.05) + c(0.05, 0.01)[as.factor(sex)] + c(0.01,0.05,0.1)[rx] + 0.1*age)

h <- exp(linP)

simFuncC <- function() {
  cens <- runif(n) # simulating censoring times
  toe <- -log(runif(n))/h # simulating times of events
  event <- ifelse(toe <= cens, 1, 0) # deciding if time of event or censoring is the smallest
  data$time <- pmin(toe, cens)
  data$event <- event
  mCox <- coxph(Surv(time, event) ~ sex + rx + age, data = data)
  mPois <- glm(event ~ sex + rx + age, data = data, offset = log(time))
  c(coef(mCox), coef(mPois))
}

sim <- t(replicate(1000, simFuncC()))
colMeans(sim)

Untuk model Cox rata-rata estimasi parameter adalah

        sex       rxLev   rxLev+5FU         age 
-0.03826301  0.04167353  0.09069553  0.10025534 

dan untuk model Poisson

(Intercept)         sex       rxLev   rxLev+5FU         age 
-1.23651275 -0.03822161  0.03678366  0.08606452  0.09812454 

Untuk kedua model, kita melihat bahwa ini dekat dengan nilai sebenarnya, mengingat bahwa perbedaan antara pria dan wanita adalah -0,04, misalnya, dan diperkirakan -0,038 untuk kedua model. Kita sekarang dapat melakukan hal yang sama dengan fungsi bahaya yang tidak konstan

$\lambda(t) = \lambda_0 t \exp(\beta'x)$ .

Kami sekarang mensimulasikan seperti sebelumnya.

simFuncN <- function() {
  cens <- runif(n)
  toe <- sqrt(-log(runif(n))/h)
  event <- ifelse(toe <= cens, 1, 0)
  data$time <- pmin(toe, cens)
  data$event <- event
  mCox <- coxph(Surv(time, event) ~ sex + rx + age, data = data)
  mPois <- glm(event ~ sex + rx + age, data = data, offset = log(time))
  c(coef(mCox), coef(mPois))
}

sim <- t(replicate(1000, simFuncN()))
colMeans(sim)

Untuk model Cox sekarang kita dapatkan

        sex       rxLev   rxLev+5FU         age 
-0.04220381  0.04497241  0.09163522  0.10029121  

dan untuk model Poisson

(Intercept)         sex       rxLev   rxLev+5FU         age 
-0.12001361 -0.01937333  0.02028097  0.04318946  0.04908300

Dalam simulasi ini, rata-rata model Poisson jelas lebih jauh dari nilai sebenarnya daripada model Cox. Ini tidak mengherankan karena kita telah melanggar asumsi bahaya konstan.

Ketika bahaya konstan, fungsi survivor, , berbentuk$S$

$S(t) = \exp(-\alpha * t)$ ,

untuk beberapa positif tergantung pada subjek tertentu, sehingga adalah cembung. Jika kami menggunakan estimator Kaplan-Meier untuk mendapatkan estimasi untuk data asli, kami melihat yang berikut ini.$\alpha$$S$$S$

Fungsi ini terlihat cekung. Ini tidak membuktikan apa-apa, tetapi bisa menjadi petunjuk bahwa asumsi bahaya konstan tidak terpenuhi untuk kumpulan data ini, yang pada gilirannya dapat menjelaskan perbedaan antara kedua model.

Komentar terakhir tentang data, sejauh yang saya tahu menyimpan data tentang waktu untuk kambuhnya kanker dan waktu sampai mati (ada dua pengamatan untuk setiap nilai ). Di atas, kami telah memodelkannya seperti itu adalah hal yang sama. Itu mungkin bukan ide yang bagus.${\tt colon}$${\tt id}$