Bagaimana cara menginterpretasikan output dari predict.coxph?

17

Setelah memasang model cox, dimungkinkan untuk membuat prediksi dan mengambil risiko relatif dari data baru. Apa yang saya tidak mengerti adalah bagaimana risiko relatif dihitung untuk seorang individu dan apa itu relatif terhadap (yaitu rata-rata populasi)? Adakah rekomendasi untuk sumber daya untuk membantu memahami (saya tidak terlalu maju dalam analisis kelangsungan hidup sehingga semakin sederhana semakin baik)?

pengguna4673
sumber

Jawaban:

32

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)=dalamh0(t)+β1X1++βhalXhal=dalamh0(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:XsayaXsayat

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 ' biib1,,bpβ1,,βpeXibeXib .

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 referenceopsi yang dijelaskan di halaman bantuan.

caracal
sumber
2
Memberi +1 karena tidak jelas untuk mendapatkan apa yang diprediksi.coxph sebenarnya dari halaman bantuan!
ocram
tadi sangat menyenangkan! Sangat mudah dimengerti!
user4673
meanFin <- mean(as.numeric(Rossi$fin) - 1)tidak masuk akal, karena finkategorikal. Apakah Anda tidak perlu modeFin <- get_Mode(Rossi$fin)dalam hal ini?
Zhubarb
1
@ Zhubarb finadalah 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 yang predict.coxph()dilakukannya.
caracal
Dengan kata lain, bagaimana Anda menginterpretasikan rasio bahaya 3.01 (misalnya relRisk [2])?
RNB