predict.coxph()
menghitung rasio bahaya relatif terhadap rata-rata sampel untuk semua variabel prediktor . Faktor dikonversi menjadi prediktor dummy seperti biasa yang rata-rata dapat dihitung. Ingat bahwa model PH Cox adalah model linear untuk log-bahaya ln h ( t ) :plnh(t)
dalamh ( t ) = lnh0( t ) + β1X1+⋯ + βhalXhal= lnh0( t ) + X β
Di mana adalah bahaya garis dasar yang tidak ditentukan. Sama dengan itu, bahaya h ( t ) dimodelkan sebagai h ( t ) = h 0 ( t ) ⋅ e β 1 X 1 + ⋯ + β p X p = h 0 ( t ) ⋅ e X β . Rasio bahaya antara dua orang i dan i ′ dengan nilai prediktorh0( T )h ( t )h ( t ) = h0( t ) ⋅ eβ1X1+ ⋯ + βhalXhal= h0( t ) ⋅ eX βsayasaya′ dan X i ' demikian independen dari bahaya dasar dan independen dari waktut:XsayaXsaya′t
hi(t)hi′(t)=h0(t)⋅eXiβh0(t)⋅eXi′β=eXiβeXi′β
Untuk rasio bahaya diperkirakan antara orang dan i ' , kita hanya pasang di estimasi koefisien b 1 , ... , b p untuk β 1 , ... , β p , memberikan e X i b dan e X i ' bii′b1,…,bpβ1,…,βpeXibeXi′b .
Sebagai contoh dalam R, saya menggunakan data dari lampiran John Fox pada model Cox-PH yang menyediakan teks pengantar yang sangat bagus. Pertama, kami mengambil data dan membangun model Cox-PH sederhana untuk penangkapan tahanan yang dibebaskan ( fin
: faktor - menerima bantuan keuangan dengan pengkodean dummy "no"
-> 0, "yes"
-> 1 age
,: usia pada saat pembebasan, prio
: jumlah hukuman sebelumnya):
> URL <- "http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt"
> Rossi <- read.table(URL, header=TRUE) # our data
> Rossi[1:3, c("week", "arrest", "fin", "age", "prio")] # looks like this
week arrest fin age prio
1 20 1 no 27 3
2 17 1 no 18 8
3 25 1 no 19 13
> library(survival) # for coxph()
> fitCPH <- coxph(Surv(week, arrest) ~ fin + age + prio, data=Rossi) # Cox-PH model
> (coefCPH <- coef(fitCPH)) # estimated coefficients
finyes age prio
-0.34695446 -0.06710533 0.09689320
eXb
meanFin <- mean(as.numeric(Rossi$fin) - 1) # average of financial aid dummy
meanAge <- mean(Rossi$age) # average age
meanPrio <- mean(Rossi$prio) # average number of prior convictions
rMean <- exp(coefCPH["finyes"]*meanFin # e^Xb
+ coefCPH["age"] *meanAge
+ coefCPH["prio"] *meanPrio)
eXb
r1234 <- exp(coefCPH["finyes"]*(as.numeric(Rossi[1:4, "fin"])-1)
+ coefCPH["age"] *Rossi[1:4, "age"]
+ coefCPH["prio"] *Rossi[1:4, "prio"])
Sekarang hitung risiko relatif untuk 4 orang pertama terhadap rata-rata sampel dan bandingkan dengan hasil dari predict.coxph()
.
> r1234 / rMean
[1] 1.0139038 3.0108488 4.5703176 0.7722002
> relRisk <- predict(fitCPH, Rossi, type="risk") # relative risk
> relRisk[1:4]
1 2 3 4
1.0139038 3.0108488 4.5703176 0.7722002
Jika Anda memiliki model bertingkat, perbandingan dalam predict.coxph()
bertentangan dengan rata-rata strata, ini dapat dikontrol melalui reference
opsi yang dijelaskan di halaman bantuan.
meanFin <- mean(as.numeric(Rossi$fin) - 1)
tidak masuk akal, karenafin
kategorikal. Apakah Anda tidak perlumodeFin <- get_Mode(Rossi$fin)
dalam hal ini?fin
adalah biner, jadi representasi numerik dari faktor hanya memiliki nilai 1 dan 2. Mengurangkan 1 memberi kita variabel dummy-kode dengan nilai 0 dan 1 yang juga muncul dalam matriks desain. Perhatikan bahwa ini tidak akan berfungsi untuk faktor dengan lebih dari 2 level. Tentu saja masih bisa diperdebatkan apakah rata-rata variabel dummy masuk akal, tetapi itulah yangpredict.coxph()
dilakukannya.