Email : elvrianae@gmail.com
Instagram : elvriana_e
RPubs : https://rpubs.com/elvriana_e
GitHub : https://github.com/elvrianae

1 Analisis Risiko Bersaing

Pada Bagian 2, kita telah membahas tentang pengenalan konsep dan metode analisis survival. Di sini, kita fokus pada Competing Risk Analysis (CRA) dengan memilih model dan menilai kecukupan dan kecocokannya. Paket utama yang digunakan dalam analisis risiko bersaing adalah cmprsk.

library(cmprsk)

1.1 Apakah CRA itu?

Ketika subjek memiliki beberapa kemungkinan peristiwa dalam pengaturan waktu ke peristiwa.

Contohnya:

Kambuh
Kematian karena penyakit
Kematian karena penyebab lain
Respon pengobatan

Semua atau sebagian dari ini (di antara yang lain) mungkin merupakan peristiwa yang mungkin terjadi dalam pembelajaran tertentu.

1.2 Jadi Apakah Masalahnya?

Ketergantungan yang tidak teramati di antara waktu peristiwa adalah masalah mendasar yang mengarah pada kebutuhan akan pertimbangan khusus. Sebagai contoh, dapat dibayangkan bahwa pasien yang kambuh lebih mungkin meninggal, dan oleh karena itu waktu untuk kambuh dan waktu kematian bukanlah peristiwa yang independen.

1.3 Latar Belakang CRA

Dua pendekatan analisis dengan adanya beberapa hasil potensial:

Bahaya penyebab khusus dari peristiwa tertentu: ini mewakili tingkat per unit waktu peristiwa di antara mereka yang tidak gagal dari peristiwa lain
Insiden kumulatif dari peristiwa tertentu: ini mewakili rasio per unit dari waktu peristiwa serta pengaruh dari peristiwa yang bersaing

Masing-masing pendekatan ini hanya dapat menerangi satu aspek penting dari data sementara mungkin mengaburkan yang lain, dan pendekatan yang dipilih harus bergantung pada pertanyaan yang menarik.

1.4 Catatan dan Referensi CRA

Saat kita mengerjakan Analisis Risiko Bersaing, berikut beberapa catatan dan referensi tambahan yang harus kita pertimbangkan:

Ketika peristiwanya independen (hampir tidak pernah benar), penyebab spesifik risiko tidak bias
Saat peristiwa bergantung, berbagai hasil dapat diperoleh bergantung pada pengaturan
Insiden kumulatif menggunakan Kaplan-Meier selalu \(\geq\) insiden kumulatif menggunakan metode risiko yang bersaing, sehingga hanya dapat menyebabkan perkiraan berlebih dari insiden kumulatif, jumlah overestimasi tergantung pada tingkat peristiwa dan ketergantungan antar peristiwa
Untuk menetapkan bahwa kovariat memang bekerja pada peristiwa yang diinginkan, penyebab spesifik risiko mungkin lebih disukai untuk pengobatan atau pengujian efek penanda pronostik
Untuk menetapkan manfaat secara keseluruhan, bahaya subdistribusi mungkin lebih disukai untuk membangun nomogram prognostik atau mempertimbangkan efek ekonomi kesehatan untuk mendapatkan pemahaman yang lebih baik tentang pengaruh pengobatan dan kovariat lainnya pada skala absolut

Beberapa referensi yang bisa kita dekati untuk mempelajari lebih lanjut:

Dignam JJ, Zhang Q, Kocherginsky M. The use and interpretation of competing risks regression models. Clin Cancer Res. 2012;18(8):2301-8.
Kim HT. Cumulative incidence in competing risks data and competing risks regression analysis. Clin Cancer Res. 2007 Jan 15;13(2 Pt 1):559-65.
Satagopan JM, Ben-Porat L, Berwick M, Robson M, Kutler D, Auerbach AD. A note on competing risks in survival data analysis. Br J Cancer. 2004;91(7):1229-35.
Austin, P., & Fine, J. (2017). Practical recommendations for reporting Fine‐Gray model analyses for competing risk data. Statistics in Medicine, 36(27), 4391-4400.

1.5 Insiden Kumulatif untuk CRA

Estimasi non-parametrik dari insiden kumulatif
Memperkirakan insiden kumulatif peristiwa yang diminati
Kapan pun jumlah insiden kumulatif setiap peristiwa sama dengan total insiden kumulatif peristiwa apa pun (tidak benar dalam pengaturan spesifik penyebab)
Uji Gray adalah uji Chi-squared yang dimodifikasi yang digunakan untuk membandingkan 2 kelompok atau lebih

2 Contoh: Data Melanoma

Kita menggunakan data Melanoma dari paket MASS untuk menggambarkan konsep ini. Ini berisi variabel:

time
Kelangsungan hidup waktu dalam beberapa hari, mungkin disensor.
status
1 : Meninggal karena melanoma,
2 : Hidup,
3 : Meninggal karena sebab lain.
sex
1 : Laki-laki,
0 : Perempuan.
age
Usia dalam tahun.
year
Tahun operasi.
thickness
Ketebalan tumor dalam mm.
ulcer
1 : Ada,
0 : Tidak ada.

data(Melanoma, package = "MASS")

2.1 Insiden Kumulatif

Perkirakan insiden kumulatif dalam konteks risiko bersaing menggunakan fungsi cuminc.

Catatan: Dalam data Melanoma, pasien yang disensor diberi kode \(2\) untuk status, jadi kita tidak dapat menggunakan opsi cencode default \(0\)

cuminc(Melanoma$time, Melanoma$status, cencode = 2)

## Estimates and Variances:
## $est
##           1000       2000       3000      4000      5000
## 1 1 0.12745714 0.23013963 0.30962017 0.3387175 0.3387175
## 1 3 0.03426709 0.05045644 0.05811143 0.1059471 0.1059471
## 
## $var
##             1000         2000         3000        4000        5000
## 1 1 0.0005481186 0.0009001172 0.0013789328 0.001690760 0.001690760
## 1 3 0.0001628354 0.0002451319 0.0002998642 0.001040155 0.001040155

2.2 Plot Insiden Kumulatif

Hasilkan plot R dasar dengan semua default.

ci_fit <- 
  cuminc(
    ftime = Melanoma$time, 
    fstatus = Melanoma$status, 
    cencode = 2
    )

plot(ci_fit)

Dalam legenda:

Angka pertama menunjukkan kelompok, dalam hal ini hanya ada perkiraan keseluruhan jadi \(1\) untuk keduanya
Angka kedua menunjukkan jenis peristiwa, dalam hal ini garis padat adalah \(1\) untuk kematian akibat melanoma dan garis putus-putus adalah \(3\) untuk kematian akibat penyebab lain

Kita juga dapat memplot insiden kumulatif menggunakan fungsi ggscompetingrisks dari paket survminer.

Dalam hal ini kita mendapatkan panel yang berlabel sesuai dengan grup, dan legenda berlabel peristiwa, yang menunjukkan jenis peristiwa untuk setiap baris.

Catatan

Kita dapat menggunakan opsi multiple_panels = FALSE untuk membuat semua grup diplot pada satu panel
tidak seperti R, sumbu \(y\) tidak pergi ke \(1\) secara default, jadi kita harus mengubahnya secara manual

library(survminer)
ggcompetingrisks(ci_fit)

2.3 Membandingkan Insiden Kumulatif

Dalam fungsi cuminc, tes Gray digunakan untuk tes antar-kelompok.

Sebagai contoh, bandingkan hasil Melanoma menurut ulcer, ada atau tidaknya ulserasi. Hasil tes dapat ditemukan di Tests.

ci_ulcer <- 
  cuminc(
    ftime = Melanoma$time, 
    fstatus = Melanoma$status, 
    group = Melanoma$ulcer,
    cencode = 2
    )

ci_ulcer[["Tests"]]

##        stat           pv df
## 1 26.120719 3.207240e-07  1
## 3  0.158662 6.903913e-01  1

Kita juga dapat memvisualisasikan atau memplot insiden kumulatif berdasarkan kelompok dengan menggunakan paket ggcompetingrisks:

ggcompetingrisks(
  fit = ci_ulcer, 
  multiple_panels = FALSE,
  xlab = "Days",
  ylab = "Cumulative incidence of event",
  title = "Death by ulceration",
  ylim = c(0, 1)
)

Catatan Saya menemukan fungsi ggcompetingrisks kurang dalam penyesuaian, terutama dibandingkan dengan ggsurvplot. Saya biasanya melakukan plotting saya, dengan terlebih dahulu membuat kumpulan data yang rapi dari hasil cuminc fit, dan kemudian memplot hasilnya. Lihat kode sumber untuk presentasi ini untuk rincian kode yang mendasari.

library(tidyverse)
ciplotdat <- 
  ci_ulcer %>% 
  list_modify("Tests" = NULL) %>% 
  map_df(`[`, c("time", "est"), .id = "id") %>% 
  mutate(id = recode(
    id, 
    "0 1" = "Not ulcerated:Death melanoma", 
    "0 3" = "Not ulcerated:Death other causes", 
    "1 1" = "Ulcerated:Death melanoma",
    "1 3" = "Ulcerated:Death other causes")
    ) %>% 
  separate(id, c("Ulceration", "Event"), ":") 

ggplot(ciplotdat, aes(x = time, y = est, color = Ulceration)) +
  geom_step(lwd = 1.2, aes(linetype = Event))  +
  ylim(c(0, 1)) +
  theme_classic() +
  theme(plot.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.position = "bottom") +
  labs(x = "Days", 
       y = "Cumulative incidence",
       title = "Death by ulceration status") +
  annotate("text", x = 0, y = 1, hjust = 0,
           label = paste0(
             "Death melanoma p = ", 
             ifelse(ci_ulcer$Tests[1, 2] < .001, 
                    "<.001", 
                    round(ci_ulcer$Tests[1, 2], 3)))) + 
  annotate("text", x = 0, y = 0.92, hjust = 0,
           label = paste0(
             "Death other causes p = ", 
             ifelse(ci_ulcer$Tests[2, 2] < .001, 
                    "<.001", 
                    round(ci_ulcer$Tests[2, 2], 3))))

Seringkali hanya satu dari jenis acara yang menarik, meskipun kita masih ingin memperhitungkan acara yang bersaing. Dalam hal ini acara yang menarik dapat direncanakan sendiri. Sekali lagi, saya melakukan ini secara manual dengan terlebih dahulu membuat kumpulan data yang rapi dari hasil cuminc fit, dan kemudian memplot hasilnya. Lihat kode sumber untuk presentasi ini untuk rincian kode yang mendasari.

ciplotdat1 <- 
  ci_ulcer %>% 
  list_modify("Tests" = NULL) %>% 
  map_df(`[`, c("time", "est"), .id = "id") %>% 
  filter(id %in% c("0 1", "1 1")) %>% 
  mutate(Ulceration = recode(
    id, 
    "0 1" = "Not ulcerated", 
    "1 1" = "Ulcerated")
    )

ggplot(ciplotdat1, aes(x = time, y = est, color = Ulceration)) +
  geom_step(lwd = 1.2)  +
  ylim(c(0, 1)) +
  theme_classic() +
  theme(plot.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.position = "bottom") +
  labs(x = "Days", 
       y = "Cumulative incidence",
       title = "Melanoma death by ulceration status") +
  annotate("text", x = 0, y = 1, hjust = 0,
           label = paste0(
             "p-value = ", 
             ifelse(ci_ulcer$Tests[1, 2] < .001, 
                    "<.001", 
                    round(ci_ulcer$Tests[1, 2], 3))))

2.4 Menambahkan Angka di Tabel Risiko

Kita mungkin ingin menambahkan jumlah tabel risiko ke plot insiden kumulatif, dan tidak ada cara mudah untuk melakukan ini yang saya ketahui. Lihat kode sumber untuk presentasi ini untuk satu contoh (dengan permintaan populer, kode sumber sekarang disertakan langsung di bawah untuk satu contoh spesifik)

Dapatkan plot dari R, ggcompetingrisks, atau ggplot

mel_plot <- 
  ggplot(ciplotdat1, aes(x = time, y = est, color = Ulceration)) +
  geom_step(lwd = 1.2)  +
  ylim(c(0, 1)) +
  coord_cartesian(xlim = c(0, 5000)) +
  scale_x_continuous(breaks = seq(0, 5000, 1000)) +
  theme_classic() +
  theme(plot.title = element_text(size = 14),
        legend.title = element_blank(),
        legend.position = "bottom") +
  labs(x = "Days", 
       y = "Cumulative incidence",
       title = "Melanoma death by ulceration status") +
  annotate("text", x = 0, y = 1, hjust = 0,
           label = paste0(
             "p-value = ", 
             ifelse(ci_ulcer$Tests[1, 2] < .001, 
                    "<.001", 
                    round(ci_ulcer$Tests[1, 2], 3))))

Dapatkan tabel nomor berisiko dari ggsurvplot menggunakan survfit di mana semua peristiwa dihitung sebagai titik akhir komposit tunggal
- Paksa sumbu untuk memiliki batas dan jeda dan judul yang sama
- Pastikan warna / jenis garis cocok untuk label grup
- Cobalah untuk mendapatkan ukuran font yang sama

mel_fit <- survfit(
  Surv(time, ifelse(status != 2, 1, 0)) ~ ulcer, 
  data = Melanoma
)

num <- ggsurvplot(
  fit = mel_fit, 
  risk.table = TRUE, 
  risk.table.y.text = FALSE,
  ylab = "Days",
  risk.table.fontsize = 3.2,
  tables.theme = theme_survminer(font.main = 10),
  title = "Test"
  )

Kemudian gabungkan plot dan risktable. Saya menggunakan fungsi plot_grid dari paket cowplot untuk ini
- Terima kasih kepada dosen saya yang telah mengirimi saya tentang cara mengubah ukuran teks yang bertuliskan “Number at risk”! Saya menggunakan yang disarankan oleh Charles Champeaux, diimplementasikan di atas dalam baris tabel.theme = theme_survminer(font.main = 10)!

cowplot::plot_grid(
  mel_plot, 
  num$table + theme_cleantable(), 
  nrow = 2, 
  rel_heights = c(4, 1), 
  align = "v", 
  axis = "b"
  )

2.5 Risiko yang Bersaing

Dua pendekatan:

Penyebab spesifik risiko
- tingkat peristiwa seketika dari jenis peristiwa tertentu dalam subjek yang saat ini bebas peristiwa
- diperkirakan menggunakan regresi Cox (fungsi coxph)
Bahaya subdistribusi
- tingkat kemunculan sesaat dari jenis peristiwa tertentu pada subjek yang belum mengalami peristiwa jenis itu
- diperkirakan menggunakan regresi Fine-Grey (fungsi crr)

2.6 Bahaya Subdistribusi

Katakanlah kita tertarik untuk melihat pengaruh usia dan jenis kelamin pada kematian akibat melanoma, dengan kematian akibat penyebab lain sebagai acara yang bersaing.

Catatan:

crr membutuhkan spesifikasi kovariat sebagai matriks
Jika lebih dari satu peristiwa yang menarik, kita dapat meminta hasil untuk peristiwa yang berbeda dengan menggunakan opsi failcode, secara default hasil dikembalikan untuk failcode = 1

shr_fit <- 
  crr(
    ftime = Melanoma$time,
    fstatus = Melanoma$status,
    cov1 = Melanoma[, c("sex", "age")],
    cencode = 2
    )

shr_fit

## convergence:  TRUE 
## coefficients:
##     sex     age 
## 0.58840 0.01259 
## standard errors:
## [1] 0.271800 0.009301
## two-sided p-values:
##  sex  age 
## 0.03 0.18

Dalam contoh sebelumnya, baik sex dan age diberi kode sebagai variabel numerik. Fungsi crr tidak dapat menangani variabel karakter secara alami, dan kita akan mendapatkan error, jadi jika variabel karakter ada, kita harus membuat variabel dummy menggunakan model.matrix

# Membuat contoh variabel karakter
chardat <- 
  Melanoma %>% 
  mutate(
    sex_char = ifelse(sex == 0, "Male", "Female")
  )

# Membuat variabel dummy dengan model.matrix
# [, -1] untuk menghapus intersep
covs1 <- model.matrix(~ sex_char + age, data = chardat)[, -1]

# Sekarang kita bisa meneruskannya ke argumen cov1, dan itu akan berhasil
crr(
  ftime = chardat$time,
  fstatus = chardat$status,
  cov1 = covs1,
  cencode = 2
  )

## convergence:  TRUE 
## coefficients:
## sex_charMale          age 
##     -0.58840      0.01259 
## standard errors:
## [1] 0.271800 0.009301
## two-sided p-values:
## sex_charMale          age 
##         0.03         0.18

2.7 Memformat Hasil dari `crr`

Hasil dari crr saat ini tidak didukung oleh broom::tidy() atau gtsummary::tbl_regress(). Sebagai alternatif, coba mvcrrres dari paket ezfun (tidak fleksibel, tetapi lebih baik)

library(kableExtra)
ezfun::mvcrrres(shr_fit) %>% 
  kable()

	HR (95% CI)	p-value
sex	1.8 (1.06, 3.07)	0.03
age	1.01 (0.99, 1.03)	0.18

2.8 Pendekatan Penyebab Spesifik Risiko

Menyensor semua subjek yang tidak memiliki peristiwa yang diinginkan, dalam hal ini kematian akibat melanoma, dan menggunakan coxph seperti sebelumnya. Jadi, pasien yang meninggal karena penyebab lain sekarang disensor untuk pendekatan penyebab spesifik risiko terhadap risiko yang bersaing.

Hasil dapat diformat dengan broom::tidy() atau gtsummary::tbl_regresi()

chr_fit <- 
  coxph(
    Surv(time, ifelse(status == 1, 1, 0)) ~ sex + age, 
    data = Melanoma
    )

broom::tidy(chr_fit, exp = TRUE) %>% 
  kable()

term	estimate	std.error	statistic	p.value
sex	1.818949	0.2676386	2.235323	0.0253961
age	1.016679	0.0086628	1.909514	0.0561958

gtsummary::tbl_regression(chr_fit, exp = TRUE)

Characteristic	HR¹	95% CI¹	p-value
sex	1.82	1.08, 3.07	0.025
age	1.02	1.00, 1.03	0.056
¹ HR = Hazard Ratio, CI = Confidence Interval

Analisis Risiko Bersaing

Pemilihan Model dan Menilai Kecukupan dan Kesesuaiannya

Elvriana Elvani

February 23, 2021

1 Analisis Risiko Bersaing

1.1 Apakah CRA itu?

1.2 Jadi Apakah Masalahnya?

1.3 Latar Belakang CRA

1.4 Catatan dan Referensi CRA

1.5 Insiden Kumulatif untuk CRA

2 Contoh: Data Melanoma

2.1 Insiden Kumulatif

2.2 Plot Insiden Kumulatif

2.3 Membandingkan Insiden Kumulatif

2.4 Menambahkan Angka di Tabel Risiko

2.5 Risiko yang Bersaing

2.6 Bahaya Subdistribusi

2.7 Memformat Hasil dari `crr`

2.8 Pendekatan Penyebab Spesifik Risiko

Analisis Risiko Bersaing

Pemilihan Model dan Menilai Kecukupan dan Kesesuaiannya

Elvriana Elvani

February 23, 2021

1 Analisis Risiko Bersaing

1.1 Apakah CRA itu?

1.2 Jadi Apakah Masalahnya?

1.3 Latar Belakang CRA

1.4 Catatan dan Referensi CRA

1.5 Insiden Kumulatif untuk CRA

2 Contoh: Data Melanoma

2.1 Insiden Kumulatif

2.2 Plot Insiden Kumulatif

2.3 Membandingkan Insiden Kumulatif

2.4 Menambahkan Angka di Tabel Risiko

2.5 Risiko yang Bersaing

2.6 Bahaya Subdistribusi

2.7 Memformat Hasil dari crr

2.8 Pendekatan Penyebab Spesifik Risiko

2.7 Memformat Hasil dari `crr`