Email:
RPubs: https://rpubs.com/Ayoung1403

library(dplyr)      #data manipulation
library(ggplot2)    #data visualization
library(DT)         #beautify the data table
library(plotly)     #make a pie chart
library(knitr)

1 Menyelesaikan Analisis Resiko

Di sesi terakhir part 2, kita sudah berdiskusi tentang pengenalan kepada konsep dan metode didalam analisis survival. Disini, kita fokus kepada Analisis Resiko Bersaing (CRA / Competing Risk Analysis) dengan memilih suatu model dan menghitung nilai kecukupan dan nilai kesesuaiannya. Paket utama yang digunakan untuk analisis resiko bersaing adalah cmprsk.

library(cmprsk)

1.1 Apa itu CRA?

Disaat suatu subjek memiliki beberapa kemungkinan peristiwa yang terjadi dalam pegaturan waktu ke peristiwa.

Contoh:

  • Kambuh

  • Kematian karena penyakit

  • Kematian karena penyebab lain

  • Respon dari suatu pengobatan

Semua atau beberapa dari contoh ini (diantara yang lain) mungkin menjadi peristiwa yang terjadi didalam studi tertentu.

1.2 Jadi Apa Masalahnya?

Ketergantungan yang tidak ditemukan diantara waktu terjadinya suatu peristiwa adalah masalah mendasar yang menyebabkan perlunya dilakukan pertimbangan yang khusus. Untuk contoh, dapat kita bayangkan bahwa pasien-pasien yang penyakitnya kambuh lebih memungkinkan akan meninggal, dan oleh sebab itu, waktu kambuhnya suatu penyakit dan waktu kematian bukanlah suatu bentuk peristiwa yang independen.

1.3 Latar belakang CRA

Dua analisis pendekatan dengan adanya beberapa hasil yang potensial:

  1. Penyebab bahaya yang spesifik dari suatu peristiwa tertentu: ini mewakili tingkatan per unit waktu didalam suatu peristiwa diantara mereka yang tidak gagal dari peristiwa lainnya.

  2. Insiden kumulatif dari peristiwa tertentu: ini mewakili peluang per unit waktu didalam suatu peristiwa serta pengaruh dari peristiwa yang saling bersaing.

Masing-masing dari pendekatan ini mungkin hanya menerangi satu aspek penting dari datanya disaat terdapat kemungkinan menghalang aspek yang lain, dan pendekatan yang dipilih harus bergantung kepada pertanyaan yang menarik.

1.4 Catatan dan Referensi CRA

Disaat kita sedang mengerjakan Analisis Resiko Bersaing, aberikut ini beberapa catatan dan referensi tambahan yang perlu kita pertimbangkan:

  • Disaat peristiwanya independen (hampir tidak pernah benar), penyebab bahaya yang spesifik bersifat tidak bias.

  • Disaat peristiwanya dependen, variasi dari hasilnya dapat diperoleh dengan bergantung kepada pengaturannya.

  • Insiden kumulatif yang menggunakan metode Kaplan-Meier selalu lebih besar atau sama dengan insiden kumulatif yang menggunakan metode resiko bersaing, sehingga hanya dapat menyebabkan perkiraan yang berlebih dari insiden kumulatif tersebut, dan jumlah perkiraan yang berlebih tersebut bergantung kepada peluang terjadinya peristiwa dan bergantung antar peristiwa-peristiwa lainnya.

  • Untuk menetapkan bahwa kovariat memang bekerja didalam suatu peristiwa yang diinginkan, penyebab bahaya yang spesifik mungkin lebih disukai untuk pengobatan atau pengujian efek pendanda pronostik.

  • Untuk menetapkan keseluruhan keuntungan, bahaya subdistribusi mungkin lebih disukai untuk membangun nomogram prognostik atau mempertimbangkan efek kesehatan ekonomi untuk pendapatkan pemahaman yang lebih baik tentang pengaruh pengobatan dan kovariat lainnya didalam skala yang absolut.

Beberapa referensi yang bisa kita temukan untuk belajar lebih lanjut adalah:

  • 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 dari suatu peristiwa yang diminati.

  • Berapapun jumlah insiden kumulatif setiap peristiwa, sama dengan total insiden kumulatif peristiwa apapun (tidak benar dalam penyebab pengaturan spesifik).

  • Uji Gray adalah uji Chi-squared yang digunakan untuk membandingkan dua kelompok atau lebih.

2 Contoh: Dataset Melanoma

Kita menggunakan data Melanoma dari paket MASS untuk mengilustrasikan konsep ini. Ini berisikan variabel-variabel:

  • time waktu bertahan hidup dalam hari, kemungkinan tersensor.

  • status 1 meninggal dari melanoma, 2 hidup, 3 meninggal karena penyebab lainnya.

  • sex 1 = laki-laki, 0 = perempuan.

  • age umur dalam tahun.

  • year tahun dari operasinya.

  • thickness ketebalan tumor dalam milimeter.

  • ulcer 1 = presensi, 0 = absensi.

data(Melanoma, package = "MASS")

2.1 Insiden Kumulatif

Perkirakan insiden kumulatif dalam konteks resiko bersaing dengan menggunakan fungsi kumulatif.

Catatan: didalam data Melanoma, pasien yang disensor diberikan kode sebagai 2 untuk status, jadi kita tidak dapat menggunakan opsi kegagalan cencode dari 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

Buatlah plot dasar R dengan semua kegagalan

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

Dalam penjelasan grafik diatas:

  • Nomor pertama mengindikasi kelompok, dikasus ini, hanya ada keseluruhan perkiraan, jadi 1 untuk keduanya.

  • Nomor kedua mengindikasi tipe peristiwanya, dikasus ini, garis yang tebal adalah 1 untuk kematian dari melanoma dan garis terputus-putus adalah 3 untuk kematian dari penyebab lainnya.

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

Di kasus ini, kita mendapatkan panel yang diberi label yang sesuai dengan kelompoknya, dan judul yang dilabeli peristiwanya, menunjukkan tipe peristiwa untuk setiap baris.

Catatan:

  • Anda dapat menggunakan opsi multiple_panels = FALSE untuk memiliki semua kelompok yang sudah di plotkan disetiap panel.

  • Tidak seperti basis R, sumbu y tidak pergi ke 1 dengan kegagalan, jadi Anda harus menggantinya secara manual.

library(survminer)
ggcompetingrisks(ci_fit)

2.3 Perbandingan Insiden Kumulatif

Sebagai contoh, bandingkan data Melanoma yang dihasilkan yang didasari dari ulcer, ada atau tidaknya Ulserasi. Hasil dari ujinya dapat ditemukan didalam 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 membuat plit insiden kumulatifnya yang berdasarkan kelompoknya 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: Menurut saya pribadi, menemukan fungsi ggcompetingrisks kurang bagus dalam penyesuaian, terutama dibandingkan dengan ggsurvplot. Saya biasanya melakukan plotting saya sendiri, dengan terlebih dahulu membuat dataset yang rapi dari hasil cuminc, dan kemudian memplot hasilnya. Lihat kode sumber untuk presentasi ini untuk rincian kode yang mendasari.

library(ggplot2)
library(tidyr)
library(readr)
library(tibble)
library(stringr)
library(forcats)
library(purrr)
library(dplyr)
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 peristiwa yang menarik, meskipun kita masih ingin memperhitungkan persitiwa yang saling 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, dan kemudian memplot hasilnya. Lihat sumber kode 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 = "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 pada Tabel Resiko

Anda mungkin ingin menambahkan jumlah tabel resiko ke plot insiden kumulatif, dan yang saya ketahui, tidak ada cara yang mudah untuk melakukan ini. Lihat kode sumber presentasi ini untuk satu contoh (dengan permintaan populer, kode sumber sekarang disertakan langsung di bawah untuk satu contoh 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.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 beresiko dari ggsurvplot menggunakan survfit di mana semua peristiwa dihitung sebagai titik akhir komposit tunggal,

  • Paksa sumbu memiliki batas dan jeda serta judul yang sama

  • Pastikan warna / jenis garis cocok dengan label grup

  • Cobalah untuk mendapatkan ukuran font yang sama

library(ggplot2)
library(survminer)
library(ggfortify)
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"
)

Lalu gabungkan plot dan tabel resiko. Saya menggunakan fungsi plot_grid dari paket cowplot untuk ini

Terima kasih kepada dosen saya yang telah mengirimi saya email berisi tip tentang cara mengubah ukuran teks yang bertuliskan “Number at risk”! Saya menggunakan yang disarankan oleh Charles Champeaux, diimplementasikan di atas dalam tabels.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 bahaya

    • Tingkat kejadian 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 tersebut.

    • 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, Anda 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

Di contoh sebelumnya, kedua sex dan age diberi kode sebagai variabel numerik. Fungsi crr tidak dapat menangani variabel karakter secara alami, dan Anda 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] menghilangkan pemotongan
covs1 <- model.matrix(~ sex_char + age, data = chardat)[, -1]

# Sekarang kita dapat melewatinya 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

Keluaran dari crr saat ini tidak didukung oleh broom::tidy() atau gtsummary::tbl_regress(). Sebagai alternatif, coba mvcrrres (tidak fleksibel, tapi lebih baik dari pada tidak sama sekali?) Dari paket ezfun 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 Penyebab munculnya bahaya yang spesifik

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 bahaya spesifik penyebab terhadap resiko yang bersaing.

Hasil dapat diformat dengan groom::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