Email: widiyantihoen00@gmail.com
RPubs: https://rpubs.com/widyaoen/

1 Analisis Risiko Bersaing

Pada bagian terakhir bab 2 kita telah membahas tentang pengenalan konsep dan metode pada analisis kelangsungan hidup. Sekarang kita akan fokus pada Analisi Risiko Bersaing (CRA) dengan memilih model dan menilai kecukupan dan kecocokannya , package utama yang digunakan dalam analisis risiko bersaing adalah cmprsk.

library(cmprsk)

## Loading required package: survival

1.1 CRA adalah…

Ketika subjek memiliki beberapa kemungkinan dalam pengaturan waktu kejadian

Contoh :

Penyakit berulang
Kematian karena penyakit
Kematian karena penyebab lain
Respon pengobatan

Semua atau sebagian (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 peristiwa adalah masalah mendasar yang mengarah pada kebutuhan akan pertimbangan khusus. Sebagai contoh, dapat dibayangkan bahwa pasie yang kambuh dan kemungkinan meninggal, oleh karena itu waktu untuk kambuh dan waktu kematian tidak termasuk kejadian yang independen.

1.3 Latar belakang CRA

Dua pendekatan analisis dengan adanya beberapa hasil potensial :

1.Bahaya khusus penyebab peristiwa tertentu, ini mewakili tingkat per unit waktu peristiwa di antara mereka yang tidak gagal dari peristiwa lain

2.Insiden kumulatif dari peristiwa tertentu : ini mewakili rasio per unit waktu kejadian serta pengaruh peristiwa yang bersaing

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

1.4 Catatan dan referensi CRA

Saat kami mengerjakan analisis risiko bersaing, berikut beberapa

catatan dan referensi tambahan yang harus kami pertimbangkan :

Ketika kejadian independen (hampir tidak pernah benar), penyebab bahaya spesifik tidak bias.
Saat kejadian bergantung, berbagai hasil dapat diperoleh bergantung pada pengaturan.
Insiden kumulatif selalu menggunakan Kaplan-Meier dengan menggunakan metode risiko yang bersaing, sehingga hanya dapat menyebabkan perkiraan berlebih dari insiden 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 mungkin lebih disukai untuk pengobatan atau pengujian efek penanda pronostik
Untuk menetapkan manfaat secara keseluruhan, bahaya subdistribusi mungkin lebih disukai untuk membangun nomogram pronostik atau mempertimbangkan efek ekonomi kesehatan untuk mendapatkan pemahaman yang lebih baik tentang pengaruh pengobatan dan kovariat lainnya pada skala absolute. Beberapa referensi yang dapat kita dapatkan untuk lebih mempelajari :
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 untuk insiden kumulatif
Estimasi insiden kumulatif dari kejadian bunga
Kapan pun jumlah insiden kumulatif setiap peristiwa sama dengan total insiden kumulatif peristiwa apapun (tidak benar dalam pengaturan penyebab khusus)
Tes Gray adalah tes Chi-Square yang di modifikasi dan digunakan untuk membandingkan dua atau lebih kelompok.

2 Contoh : Dataset Melanoma

Kita menggunakan data Melanom dari package MASS untuk mengilustrasikan konsep ini, terdapat beberapa variable :

Time : waktu kelangsungan hidup dalam hari, kemungkinan sensoring
Status : 1 meninggal dari melanoma, 2 hidup, 3 meninggal dari penyebab lain.
Sex : 1 laki-laki, 0 perempuan
Age : umur dalam tahun
year : Tahun beroperasi
Thickness : ketebalan tumor dalam mm
Ulcer : 1 hadir , 0 tidak hadir

data(Melanoma,package = "MASS")

2.1 Kumulatif insiden

Estimasi kumulatif insiden dalam konteks risiko bersaing menggunakan fungsi cuminc

Note : dalam data Melanoma, pasien yang disensor diberi kode 2 untuk status, jika tidak dapat menggunakan opsi kode 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 dari insiden kumulatif

Hasilkan plot R dengan semua default

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

Dalam legenda

Angka pertama mengindikasikan kelompok, pada kasus ini hanya terdapat keseluruhan perkiraan menjadi 1 untuk keduanya
Angka kedua menunjukan jenis kejadia, 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 membuat plot kumulatif insiden menggunakan fungsi ggscompetingrisks dari package survminer. Dalam hal ini kita mendapatkan panel yang diberi label sesuai dengan kelompok dan legenda berlabel peristiwa, yang menunjukan jenis peristiwa untuk setiap baris.

Note

Kamu dapat menggunakan opsi multiple_panels=FALSE untuk memiliki semua kelompok plot pada panel tunggal
Tidak seperti basis R, sumbu y tidak pergi ke 1 secara default, jadi anda harus mengubahnya secara manual

library(survminer)

## Loading required package: ggplot2

## Loading required package: ggpubr

ggcompetingrisks(ci_fit)

2.3 Membandingkan insiden kumulatif

Dalam cuminc tes Gray menggunakan diantara kelompok tes.

Sebagai contoh, bandingkan hasil Melanoma menurut ulcer, ada tidaknya ulserasi. Hasil tes dapat ditemukan dalam 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 memvisualisasikan atau membuat plot dari insiden kumulatif berdasarkan kelompok 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))

Note

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

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --

## v tibble  3.0.4     v dplyr   1.0.4
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.0
## v purrr   0.3.4

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

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 kami masih ingin memperhitungkan peristiwa yang bersaing. Dalam hal ini peristiwa yang menarik dapat direncanakan sendiri. Sekali lagi, saya melakukan ini secara manual dengan terlebih dahulu membuat set data yang rapi dari hasil ‘cuminc’ , 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",
       tittle="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 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)

 Peroleh plot dari dasar R, 0ggcompetingrisks dari 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",
       tittle="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))))

 Peroleh nomor table risiko dari ‘ggsurvplot’ menggunakan ‘survfit’ dimana semua kejadian dihitung sebagai titik akhir komposit tunggal

yang memiliki batas dan jeda serta judul yang sama
Pastikan warna / jenis garis cocok untuk label grup
Coba 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")

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

 Kemudian gabungkan plot dan risktable. Saya menggunakan fungsi ‘plot_grid’ dari paket ‘cowplot’ untuk ini

Terima kasih kepada beberapa pembaca yang telah mengirimi saya email berisi tips tentang cara mengubah ukuran teks yang bertuliskan “Jumlah berisiko”! Saya menggunakan yang disarankan oleh Charles Champeaux, diterapkan di atas dalam baris tables.themes=themes_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

Regresi dua pendekatan

 Penyebab bahaya spesifik

tingkat kejadian seketika dari jenis peristiwa tertentu dalam subjek peristiwa bebas saat ini
memperkirakan menggunakan regresi Cox (fungsi ‘coxph’)

 Bahaya subdistribusi

tingkat kemunculan sesaat dari jenis peristiwa tertentu pada subjek yang belum mengalami jenis peristiwa itu
memperkirakan menggunakan regresi Fine-Gray (fungsi ‘crr’)

2.6 Bahaya Subdistribusi

Katakanlah kita tertarik untuk melihat pengaruh usia dan jenis kelamion pada kematian akibat Melanoma, dengan kematian akibat penyebab lain sebagai peristiwa yang bersaing

Notes :  crr membutuhkan spesifikasi kovariat sebagai matriks

 Jika lebih dari satu peristiwa yang menarik, anda dapat meminta hasil untuk peristiwa yang berbeda dengan menggunakan opsi kode failcode, secara default hasil dikembalikan untuk kode ‘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

Pada contoh sebelumnya, antara jenis kelamin dan umur dikodekan sebagai variable numeric. Funsi ‘crr’ tidak bias secara alami mengatur variable karakter, dan kamu akan memperoleh error, sehingga jika terdapar variable karakter kita harus membuat variable dummy menggunakan mode.matrix

#membuat contoh variabel karakter
chardat<-Melanoma %>% mutate(sex_char=ifelse(sex==0,"Male","Female"))
#membuat dummy variables dengan model.matrix
# [,-1] menghilangkan intersept
covs1<-model.matrix(~sex_char+age,data=chardat)[,-1]
# sekarang kita dapat melalui argument cov1, dan itu akan bekerja
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 tidak didukung antara broom::tidy() atau gtsummary::tbl_regression() pada waktu ini. Sebagai alternative , coba (tidak fleksibel, tetapi lebih baik daripada tidak?) mvcrres dari package saya ezfun

library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

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 spesifik bahaya

Sensor semua subjek yang tidak memiliki kejadian yang tidak diinginkan, dalam hal ini kematian akibat Melanoma, dan menggunakan coxph seperti sebelumnya. Jadi, pasien yang meninggal karena penyebab lain sekarang di sensor untuk pendekatan penyebab spesifik bahaya terhadap risiko yang bersaing.

hasil dapat di format dengan broom::tidy() atau gtsummary::tbl_refression()

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

Widi Yantih

Februari 24, 2021