Katakanlah saya memiliki kumpulan data "kateter ginjal". Saya mencoba memodelkan kurva bertahan hidup menggunakan model Cox. Jika saya mempertimbangkan model Cox: saya perlu perkiraan bahaya awal. Dengan menggunakan fungsi paket R bawaan , saya dapat dengan mudah melakukannya seperti ini:
survival
basehaz()
library(survival)
data(kidney)
fit <- coxph(Surv(time, status) ~ age , kidney)
basehaz(fit)
Tetapi jika saya ingin menulis fungsi selangkah demi selangkah dari bahaya baseline untuk perkiraan parameter tertentu, b
bagaimana saya bisa melanjutkan? Saya mencoba:
bhaz <- function(beta, time, status, x) {
data <- data.frame(time,status,x)
data <- data[order(data$time), ]
dt <- data$time
k <- length(dt)
risk <- exp(data.matrix(data[,-c(1:2)]) %*% beta)
h <- rep(0,k)
for(i in 1:k) {
h[i] <- data$status[data$time==dt[i]] / sum(risk[data$time>=dt[i]])
}
return(data.frame(h, dt))
}
h0 <- bhaz(fit$coef, kidney$time, kidney$status, kidney$age)
Tetapi ini tidak memberikan hasil yang sama dengan basehaz(fit)
. Apa masalahnya?
Jawaban:
Rupanya,
basehaz()
sebenarnya menghitung tingkat bahaya kumulatif, bukan tingkat bahaya itu sendiri. Rumusnya adalah sebagai berikut: dengan mana menunjukkan waktu acara yang berbeda, adalah jumlah acara di , dan adalah risiko yang ditetapkan pada mengandung semua individu yang masih rentan terhadap kejadian di .Ayo coba ini. (Kode berikut hanya untuk ilustrasi dan tidak dimaksudkan untuk ditulis dengan sangat baik.)
output parsial:
Saya menduga bahwa sedikit perbedaan mungkin karena perkiraan kemungkinan parsial
coxph()
karena ikatan dalam data ...sumber
kidney$time >= y[l]
status=0
status=1
status=0
coxph
panggilan denganfit<-coxph(Surv(time, status)~age, data=kidney, method="breslow")
akan memperbaiki perbedaan dalam metode.