Emails \(\hspace{0.6cm}:\) 1.
\(\hspace{1.9cm}\) 2.
Instagram \(\hspace{0.03cm}:\) https://www.instagram.com/j.hodawya/
Linked In \(\hspace{0.2cm}:\) https://www.linkedin.com/in/juenzy-hodawya-a310ab1a3/
RPubs \(\hspace{0.67cm}:\) https://rpubs.com/JHodawya
Kaggle \(\hspace{0.6cm}:\) https://www.kaggle.com/juenzyhodawya
GitHub \(\hspace{0.52cm}:\) https://github.com/JuenzyHodawya

1 Analisis Risiko yang Bersaing

Pada bagian terakhir Bagian 2 kita telah membahas tentang pengenalan konsep dan metode analisis survival. Di sini, kita akan fokus pada Analisis Risiko yang Bersaing (Competing Risk Analysis) dengan memilih model dan menilai kecukupan (adequacy) dan kesesuaiannya. Package utama untuk digunakan dalam analisis risiko yang bersaing ini adalah cmprsk.

library(cmprsk)

1.1 Apakah CRA itu?

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

Contoh:

  • kambuh
  • kematian karena penyakit
  • kematian karena sebab lain
  • respon perlakuan

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

1.2 Jadi apa masalahnya?

Ketergantungan yang tidak teramati di antara waktu kejadian adalah masalah mendasar yang mengarah pada kebutuhan akan pertimbangan khusus. Sebagai contoh, dapat dibayangkan bahwa pasien yang kambuh lebih besar kemungkinannya meninggal, dan oleh karena itu waktu untuk kambuh dan waktu kematian merupakan kejadian yang dependen.

1.3 Latar belakang CRA

Dua pendekatan analisis dengan adanya beberapa hasil potensial:

\(\hspace{0.5cm}\) 1. Penyebab spesifik hazard dari peristiwa tertentu: ini mewakili tarif per unit waktu peristiwa di antara mereka yang tidak gagal dari peristiwa lain.
\(\hspace{0.5cm}\) 2. Kejadian kumulatif peristiwa tertentu: ini menunjukkan tarif per unit waktu peristiwa serta pengaruh peristiwa yang bersaing.

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

1.4 Catatan & Referensi CRA

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

  • Ketika kejadiannya independen (hampir tidak pernah benar), penyebab spesifik hazard tidak bias.
  • Jika kejadiannya dependen, berbagai hasil dapat diperoleh bergantung pada pengaturannya.
  • Kejadian kumulatif menggunakan Kaplan-Meier selalu \(\ge\) kejadian kumulatif menggunakan metode risiko yang bersaing, sehingga hanya dapat menyebabkan perkiraan yang berlebihan dari kejadian kumulatif, jumlah perkiraan yang terlalu tinggi tergantung pada tingkat kejadian dan ketergantungan antar kejadian.
  • Untuk menetapkan bahwa kovariat memang bekerja pada peristiwa yang diinginkan, penyebab spesifik hazard mungkin lebih disukai untuk pengobatan atau pengujian efek penanda pronostik.
  • Untuk menetapkan manfaat keseluruhan, subdistribusi hazard mungkin lebih disukai untuk membangun nomogram prognostik atau mempertimbangkan efek ekonomi kesehatan untuk mendapatkan pemahaman yang lebih baik tentang pengaruh perlakuan dan kovariat lainnya dalam skala absolut.

Beberapa referensi yang dapat kita gunakan sebagai pembelajaran 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 kejadian yang diminati.
  • Pada titik waktu mana pun, jumlah kejadian kumulatif setiap peristiwa sama dengan total kejadian kumulatif dari setiap peristiwa (tidak benar dalam pengaturan penyebab khusus).
  • Uji Gray adalah tes Chi-squared yang dimodifikasi yang digunakan untuk membandingkan 2 atau lebih kelompok.

2 Contoh: Dataset Melanoma

Kita menggunakan data Melanoma dari package MASS untuk mengilustrasikan konsep ini. Ini berisi variabel:

  • time waktu bertahan hidup dalam beberapa hari, mungkin akan terjadi sensor.
  • status 1 = meninggal karena melanoma, 2 = hidup, 3 = meninggal karena sebab lain.
  • sex 1 = pria, 0 = wanita.
  • age umur dalam tahun.
  • year tahun beroperasi.
  • thickness ketebalan tumor dalam mm.
  • ulcer 1 = kehadiran, 0 = tidak hadir.
data(Melanoma, package = "MASS")

2.1 Kejadian Kumulatif

Perkirakan kejadian kumulatif dalam konteks risiko yang bersaing menggunakan fungsi cuminc.

Catatan: dalam data Melanoma, pasien yang terjadi sensor diberi kode sebagai 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 Kejadian Kumulatif

Membuat plot R dasar dengan semua defaultnya.

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

Di dalam legend:

  • Angka pertama menunjukkan grup, dalam hal ini hanya ada perkiraan keseluruhan jadi 1 untuk keduanya.
  • Angka kedua menunjukkan jenis kejadian, 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 kejadian kumulatif menggunakan fungsi ggscompetingrisks dari package survminer.

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

Catatan

  • Anda dapat menggunakan opsi multiple_panels = FALSE untuk membuat semua grup diplot pada satu panel.
  • tidak seperti basis R, sumbu y tidak mengarah ke 1 secara default, jadi Anda harus mengubahnya secara manual.
library(survminer)
ggcompetingrisks(ci_fit)

2.3 Bandingkan Kejadian Kumulatif

Dalam cuminc, pengujian Gray digunakan untuk pengujian 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 kejadian kumulatif menurut kelompok dengan menggunakan package 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 pribadi merasa fungsi ggcompetingrisks masih kurang dalam penyesuaian, terutama dibandingkan dengan ggsurvplot. Saya biasanya melakukan plotting sendiri, dengan pertama-tama merapikan dataset dari hasil fit cuminc, 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 kejadian yang menarik, meskipun kita masih ingin memperhitungkan kejadian yang bersaing. Dalam hal ini kejadian yang menarik dapat diplotkan sendiri. Sekali lagi, saya melakukan ini secara manual dengan terlebih dahulu merapikan data set dari hasil fit cuminc, 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 Tambahkan Angka di Tabel Risiko

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

  • Dapatkan plot dari basis 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 nomor pada tabel risiko dari ggsurvplot menggunakan survfit di mana semua peristiwa dihitung sebagai titik akhir gabungan tunggal,
    • Paksa sumbu memiliki batas dan jeda dan judul yang sama,
    • Pastikan warna/linetypes cocok dengan 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 package cowplot untuk ini,
    • Terima kasih kepada beberapa pembaca yang telah mengirimkan saya email berisi tip tentang cara mengubah ukuran teks yang bertuliskan “Jumlah berisiko”! Saya menggunakan yang disarankan oleh Charles Champeaux, diimplementasikan di atas di baris tables.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 Regresi Resiko yang Bersaing

Dua pendekatan:

  • Penyebab spesifik hazard
    • tingkat kejadian seketika dari jenis peristiwa tertentu dalam subjek yang saat ini tidak terjadi suatu peristiwa
    • diperkirakan menggunakan regresi Cox (fungsi coxph)
  • Subdistribusi hazard
    • tingkat kejadian sesaat dari jenis peristiwa tertentu pada subjek yang belum mengalami peristiwa jenis tersebut
    • diperkirakan menggunakan regresi Fine-Grey (fungsi crr)

2.6 Subdistribusi Hazard

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

Catatan:

  • crr membutuhkan spesifikasi kovariat sebagai matriks
  • Jika lebih dari satu peristiwa yang menarik, Anda dapat meminta hasil untuk peristiwa yang berbeda dengan menggunakan opsi failcode, secara default hasil dikembalikan untukfailcode = 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

Pada contoh sebelumnya, baik sex danage dikodekan sebagai variabel numerik. Fungsi crr tidak dapat menangani variabel karakter secara alami, dan Anda akan mendapatkan error, jadi jika ada variabel karakter, kita harus membuat variabel dummy menggunakan model.matrix.

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

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

# Sekarang kita dapat meneruskannya ke argumen cov1, dan ini 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 output dari crr tidak didukung baik oleh broom::tidy() atau gtsummary::tbl_regress() untuk saat ini. Sebagai alternatif, coba Mvcrrres dari paketezfun saya.

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 hazard

Menyensor semua subjek yang tidak memiliki kejadian 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 hazard terhadap risiko yang bersaing.

Hasil dapat diformat dengan broom::tidy() atau gtsummary::tbl_regress().

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 HR1 95% CI1 p-value
sex 1.82 1.08, 3.07 0.025
age 1.02 1.00, 1.03 0.056

1 HR = Hazard Ratio, CI = Confidence Interval