Survival Model: Competing Risk Analysis

Pemilihan Model dan Menilai Kecukupan dan Kesesuaiannya

Email: imeldasianturi85@gmail.com
RPubs: https://rpubs.com/imelda123

---

1 Analisis Risiko Bersaing

Dalam Part 2 kita telah membahas mengenai “Pengenalan Konsep dan Metode Analisis Survival”. Sekarang, kita akan fokus pada Analisis Risiko Bersaing, dalam bahasa Inggris biasa dikenal sebagai Competing Risk Analysis (CRA). Kita akan memilih model dan menilai kecukupan dan kecocokannya. Package utama yang akan digunakan dalam analisis risiko bersaing (CRA) adalah cmprsk.

library(cmprsk)

1.1 Apa itu Analisis Risiko Bersaing/ Competing Risk Analysis (CRA)

Ketika subjek memiliki beberapa kemungkinan peristiwa yang dapat terjadi, dalam pengaturan waktu menuju peristiwa tertentu.

Contoh:

• Kambuh
• Meninggal karena suatu penyakit
• Meninggal karena penyebab lainnya
• Respon dari pengobatan yang dijalani

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

1.2 Jadi, Apa Masalahnya?

Ketergantungan yang tidak terobservasi (tidak teramati) diantara waktu terjadinya peristiwa adalah masalah mendasar yang membutuhkan pertimbangkan khusus. Sebagai contoh, dapat dibayangkan bahwa pasien yang kambuh lebih mungkin meninggal, oleh karena itu waktu untuk kambuh dan waktu kematian bukanlah kejadian yang independen.

1.3 Latar Belakang dari Analisis Risiko Bersaing/ Competing Risk Analysis (CRA)

Dua pendekatan analisis dengan adanya beberapa hasil yang potensial:

1. Penyebab-spesifik suatu bahaya dari peristiwa tertentu: hal ini mewakili tingkat per unit dari waktu terjadinya suatu peristiwa, diantara mereka yang tidak gagal dari peristiwa lain.

2. Insiden kumulatif dari suatu peristiwa tertentu: hal ini mewakili rasio per unit dari waktu terjadinya suatu peristiwa, serta pengaruh peristiwa yang bersaing.

Masing-masing pendekatan ini hanya dapat menjelaskan satu aspek penting dari suatu data, sementara itu pendekatan ini mungkin tidak bisa menjelaskan aspek yang lain. Pendekatan yang dipilih harus bergantung pada pertanyaan yang menjadi fokus observasi.

1.4 Catatan & Referensi Analisis Risiko Bersaing/ Competing Risk Analysis (CRA)

Saat kita mengerjakan Analisis Risiko Bersaing, berikut beberapa catatan dan referensi tambahan yang dapat menjadi bahan pertimbangan kita:

• Ketika peristiwa bersifat independen (hampir tidak pernah benar), penyebab-spesifik dari bahaya bersifat tidak bias.

• Saat peristiwa bersifat dependen, berbagai hasil dapat diperoleh bergantung pada pengaturan.

• Insiden kumulatif menggunakan Kaplan-Meier selalu >= insiden kumulatif menggunakan metode risiko yang bersaing (Competing Risk Analysis/ CRA), sehingga hanya dapat menyebabkan perkiraan berlebih dari insiden kumulatif. Jumlah overestimasi tergantung pada tingkat kejadian dan ketergantungan antar kejadian

• Untuk menetapkan bahwa kovariat memang bekerja pada peristiwa yang diinginkan, penyebab-spesifik dari bahaya mungkin lebih disukai untuk pengobatan atau pengujian efek penanda prognostik (prognostik: disiplin teknik yang berfokus pada prediksi waktu di mana sistem atau komponen tidak lagi menjalankan fungsi yang di maksudkan)

• Untuk menetapkan manfaat secara keseluruhan, bahaya subdistribusi mungkin lebih disukai untuk membangun nomogram prognostik atau mempertimbangkan efek kesehatan ekonomi, untuk mendapatkan pemahaman yang lebih baik tentang pengaruh pengobatan dan kovariat lainnya pada skala absolut

Beberapa Referensi yang dapat menjadi bahan pembelajaran untuk kita:

1. Dignam JJ, Zhang Q, Kocherginsky M. The use and interpretation of competing risks regression models. Clin Cancer Res. 2012;18(8):2301-8.

2. 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.

3. 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.

4. 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 Analisis Risiko Bersaing (Competing Risk Analysis/ CRA)

• Estimasi non-parametrik dari insiden kumulatif
• Memperkirakan insiden kumulatif peristiwa yang menjadi fokus
• Kapan pun itu, jumlah insiden kumulatif setiap peristiwa yang terjadi akan sama dengan total insiden kumulatif dari peristiwa apapun (tidak benar dalam pengaturan penyebab-spesifik)
• Tes Gray adalah tes Chi-square yang dimodifikasi, yang digunakan untuk membandingkan 2 atau lebih kelompok

2 Contoh: Kumpulan Data Melanoma

Kita menggunakan data Melanoma dari package MASS untuk membuat ilustrasi dari konsep pembelajaran kita. Variabel yang terdapat dalam data Melanoma adalah:

• time: Waktu survival (bertahan hidup) dalam harian, mungkin terjadi sensoring

• status:

  • 1 = Meninggal karena Melanoma
  • 2 = Tetap Hidup
  • 3 = Meninggal karena penyakit lain

• sex:

  • 1 = Laki-Laki
  • 0 = Perempuan

• age: Umur dalam tahun

• year: Tahun operasi

• thickness: Ketebalan tumor dalam satuan mm

• ulcer:

  • 1 = Terkena Ulcer (Kelainan pada mukosa mulut)
  • 0 = Tidak terdapat Ulcer

data(Melanoma, package = "MASS")

2.1 Insiden Kumulatif

Perkirakan insiden kumulatif dalam konteks risiko bersaing menggunakan fungsi cuminc

Catat!: Dalam data Melanoma, pasien yang tersensor ditandai dengan angka 2 untuk status nya, jadi kita tidak bisa menggunakan opsi default 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 dari Insiden Kumulatif

Hasilkan plot dasar dalam R dengan semua default.

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

plot(ci_fit)

Dalam Penjelasannya:

• 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 (1): Kematian akibat Melanoma.
  • Garis putus-putus (3): Kematian akibat penyebab lain.

Kita juga dapat membuat plot dari insiden kumulatif menggunakan fungsi ggcompetingrisks dari package survminer

Dalam kasus ini, kita mendapatkan panel yang diberi label sesuai grup, dan legenda (penjelasan) berlabel peristiwa, yang menunjukkan jenis peristiwa untuk setiap baris dalam data (setiap pasien yang tercantum dalam data)

Catat! :

• Kita dapat menggunakan opsi multiple_panels = FALSE untuk membuat semua grup diplot dalam satu panel (dalam satu bingkai grafik)
• Tidak seperti basis R, sumbu y tidak pergi ke 1 secara default, jadi kita harus mengubahnya secara manual.

library(survminer)
ggcompetingrisks(ci_fit)

2.3 Membandingkan Indsiden Kumulatif

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

Sebagai contohnya, membandingkan hasil Melanoma berdasarkan Ulcer (Ulser), yaitu 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 insiden kumulatif menurut kelompok dengan menggunakan package ggcompetingrisks

ggcompetingrisks(
  fit = ci_ulcer,
  multiple_panels = FALSE,
  xlab = "Days",
  ylab = "Insiden Kumulatif dari Peristiwa",
  title = "Kematian karena Ulserasi",
  ylim = c(0,1)
)

Catat! : Saya sendiri menemukan fungsi ggcompetingrisks kurang dalam penyesuaian, terutama dibandingkan dengan ggsurvplot. Saya biasanya melakukan plotting sendiri, dengan terlebih dauhlu membuat dataset yang tapi dari hasil fit cuminc, kemudian memplot hasilnya. Perhatikan coding di bawah, untuk presentasi ini dan untuk rincian codingnya.

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 = "Hari",
       y = "Insiden Kumulatif",
       title = "Meninggal berdasarkan status Ulser") +
  annotate("text", x = 0, y = 1, hjust = 0,
           label = paste0(
             "Meninggal karena 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(
             "Meninggal karena penyebab lain p = ",
             ifelse(ci_ulcer$Tests[2,2] < .001,
                    "<.001",
                    round(ci_ulcer$Tests[2,2], 3))))

Seringkali hanya satu dari jenis peristiwa yang menjadi fokus, meskipun kita masih ingin memperhitungkan peristiwa yang bersaing. Dalam hal ini peristiwa yang menjadi fokus dapat direncanakan sendiri. Sekali lagi, kita melakukan ini secara manual dengan terlebih dahulu membuat kumpulan data yang rapi dari hasil cuminc fit, kemudian memplot hasilnya. Lihat coding di bawah ini untuk presentasi ini untuk rincian coding yang mendasar.

ciplotdata1 <-
  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(ciplotdata1, 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 = "Hari",
       y = "Insiden Kumulatif",
       title = "Kematian Melanoma berdasarkan status Ulserasi") +
  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 Risiko

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

• Dapatkan plot dari basis R, ggcompetingrisks, atau ggplot

mel_plot <-
  ggplot(ciplotdata1, aes(x = time, y = est, color = Ulceration)) +
  geom_step(lwd = 2.1) +
  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 = "Hari",
       y = "Insiden Kumulatif",
       title = "Kematian Melanoma berdasarkan status Ulserasi") +
  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 tabel berisiko 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(survival)
mel_fit <- survfit(
  Surv(time, ifelse(status != 2, 1, 0)) ~ ulcer,
  data = Melanoma
)

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

• Kemudian gabungkan plot dan risktable. Untuk hal ini saya menggunakan fungsi plot_grid dari package cowplot.

  • Terima kasih kepada dosen saya yang telah membimbing saya tentang cara mengubah ukuran teks yang bertuliskan “Number at risk”! Dosen saya menggunakan cara yang disarankan oleh Charles Champeaux, diimplementasikan di atas dalam 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 Risiko 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 tertentu.
  • Diperkirakan menggunakan regresi Fine-Grey (fungsi crr)

2.6 Bahaya Subdistribusi

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

Catatan:

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

library(cmprsk)
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, variabel 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
# Coding [,-1] menghapus intercept
covs1 <- model.matrix(~ sex_char + age, data = chardat)[,-1]

# Sekarang, kita bisa meneruskannya ke argumen cov1, dan 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 Hasil Format dari ccr

Hasil dari ccr saat ini tidak di dukung oleh broom::tidy() atau gtsummary:tbl_regression(). Sebagai alternatif, kita bisa mencoba fungsi mvcrres dari package ezfun milik saya. (tidak fleksibel, namun 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 suatu Bahaya

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

Hasilnya dapat diformat menggunakan broom::tidy() atau juga dengan gtsummary::tbl_regression()

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