Alur prediksi berbeda dari survival coxph dan rms cph

9

Saya telah membuat versi termplot saya sendiri yang sedikit ditingkatkan yang saya gunakan dalam contoh ini, Anda dapat menemukannya di sini . Saya sebelumnya telah diposting di SO tetapi semakin saya memikirkannya, saya percaya bahwa ini mungkin lebih terkait dengan interpretasi model bahaya Cox Proportional daripada dengan pengkodean yang sebenarnya.

Masalah

Ketika saya melihat plot Rasio Bahaya saya berharap memiliki titik referensi di mana interval kepercayaan secara alami adalah 0 dan ini adalah kasus ketika saya menggunakan cph () dari rms packagetetapi tidak ketika saya menggunakan coxph () dari survival package. Apakah perilaku yang benar oleh coxph () dan jika demikian apa titik rujukannya? Juga, variabel dummy dalam coxph () memiliki interval dan nilainya selain ?e0

Contoh

Ini kode pengujian saya:

# Load libs
library(survival)
library(rms)

# Regular survival
survobj <- with(lung, Surv(time,status))

# Prepare the variables
lung$sex <- factor(lung$sex, levels=1:2, labels=c("Male", "Female"))
labels(lung$sex) <- "Sex"
labels(lung$age) <- "Age"

# The rms survival
ddist <- datadist(lung)
options(datadist="ddist")
rms_surv_fit <- cph(survobj~rcs(age, 4)+sex, data=lung, x=T, y=T)

Plot cph

Kode ini:

termplot2(rms_surv_fit, se=T, rug.type="density", rug=T, density.proportion=.05,
          se.type="polygon", yscale="exponential", log="y", 
          xlab=c("Age", "Sex"), 
          ylab=rep("Hazard Ratio", times=2),
          main=rep("cph() plot", times=2),
          col.se=rgb(.2,.2,1,.4), col.term="black")

berikan plot ini:

cph () termplot2

Plot coxph

Kode ini:

termplot2(surv_fit, se=T, rug.type="density", rug=T, density.proportion=.05, 
          se.type="polygon", yscale="exponential", log="y", 
          xlab=c("Age", "Sex"), 
          ylab=rep("Hazard Ratio", times=2),
          main=rep("coxph() plot", times=2),
          col.se=rgb(.2,.2,1,.4), col.term="black")

berikan plot ini:

coxph () termplot2

Memperbarui

Sebagai @Frank Harrell menyarankan dan setelah menyesuaikan saran dalam komentarnya baru-baru ini saya mendapat:

p <- Predict(rms_surv_fit, age=seq(50, 70, times=20), 
             sex=c("Male", "Female"), fun=exp)
plot.Predict(p, ~ age | sex,
             col="black",
             col.fill=gray(seq(.8, .75, length=5)))

Ini memberi plot yang sangat bagus ini:

Plot kisi

Saya telah melihat kontras.rms lagi setelah komentar dan mencoba kode ini yang memberikan plot ... walaupun mungkin ada lebih banyak yang dapat dilakukan :-)

w <- contrast.rms(rms_surv_fit, 
                  list(sex=c("Male", "Female"), 
                       age=seq(50, 70, times=20)))

xYplot(Cbind(Contrast, Lower, Upper) ~ age | sex, 
       data=w, method="bands")

Memberi plot ini:

Alur kontrasnya

PEMBARUAN 2

Prof. Thernau berbaik hati mengomentari plot yang kurang percaya diri:

Splines smoothing dalam coxph, seperti yang ada di gam, dinormalisasi sehingga jumlah (prediksi) = 0. Jadi saya tidak punya satu titik tetap yang variansinya ekstra kecil.

Walaupun saya belum terbiasa dengan GAM, ini sepertinya menjawab pertanyaan saya: ini sepertinya masalah interpretasi.

Max Gordon
sumber
3
Beberapa komentar. Pertama baca biostat.mc.vanderbilt.edu/Rrms untuk perbedaan antara rms dan paket Desain. Kedua, gunakan plot () alih-alih plot. Prediksi untuk menyimpan pekerjaan. Ketiga, Anda dapat dengan mudah membuat plot untuk kedua jenis kelamin, misalnya menggunakan Predict (fit, age, sex, fun = exp) # exp = anti-log; kemudian plot (hasil) atau plot (hasil, ~ usia | jenis kelamin). Anda tidak menggunakan "x = NA" di Predict. rms menggunakan grafik kisi sehingga parameter grafik par biasa dan mfrow tidak berlaku. Lihat contoh di handout kursus rms saya di biostat.mc.vanderbilt.edu/rms . Untuk contrast.rms pelajari dokumentasi lebih lanjut.
Frank Harrell
1
Terima kasih banyak atas masukan Anda. Saya telah memperbarui kode dengan contoh yang lebih baik dan menambahkan prof. Tanggapan Thernau. PS Saya sangat senang bahwa Anda merencanakan versi baru buku ini, memperluas bagian bias cutpoint akan sangat berguna sebagai referensi
Max Gordon
1
Anda dapat menggunakan plotdan contrastbukannya plot.Predictdan contrast.rms. Saya akan menggunakan byatau lengthdi dalam seqdaripada timesdan akan memberikan contrastdua daftar sehingga Anda menentukan dengan tepat apa yang sedang kontras. Anda juga dapat menggunakan shading dengan xYplotuntuk band kepercayaan diri.
Frank Harrell
1
Terima kasih. Saya suka menggunakan plot. Prediksi karena saya mendapatkan bantuan yang tepat di RStudio - sesuatu yang dalam kasus saya jauh lebih penting daripada waktu yang diperlukan untuk menulis nama fungsi lengkap (dengan menggunakan pelengkapan otomatis (tab) saya sebenarnya tidak kehilangan banyak waktu).
Max Gordon

Jawaban:

5

Saya pikir pasti ada titik di mana interval kepercayaan nol lebar. Anda juga dapat mencoba cara ketiga yaitu hanya menggunakan fungsi rms. Ada contoh di bawah file bantuan untuk contrast.rms untuk mendapatkan plot rasio bahaya. Dimulai dengan komentar # tunjukkan perkiraan terpisah berdasarkan perawatan dan jenis kelamin. Anda harus anti-log untuk mendapatkan rasio.

Frank Harrell
sumber
1
Terima kasih atas jawaban Anda. Apakah Anda pikir saya harus menyebutkan masalah ini kepada prof. Terry Therneau jika itu dianggap sebagai bug / salah tafsir? Saya juga telah melihat ke solusi grafis dalam paket rms, saya tidak bisa mengerti penggunaan contrast.rms untuk plot. The plot.Predict tampaknya melakukan keluaran yang mirip termplot tapi saya tidak bisa mendapatkannya untuk melakukan apa yang saya inginkan ... lihat pembaruan saya untuk pertanyaan.
Max Gordon
2
Akan baik untuk menulis kepadanya untuk menanyakan, dan katakan padanya terima kasih atas perjalanan ke bandara yang dia berikan kepada saya beberapa menit yang lalu. Saya akan mengomentari pertanyaan-pertanyaan lain di atas.
Frank Harrell