Heart Failure Survival Analysis

UTS Survival Model

Email         : putri.angelina@student.matanauniversity.ac.id
RPubs       : https://rpubs.com/putriangelinaw/
GitHub      : https://github.com/putriangelinaw/
Jurusan     : Statistika Bisnis
Address    : ARA Center, Matana University Tower
                      Jl. CBD Barat Kav, RT.1, Curug Sangereng, Kelapa Dua, Tangerang, Banten 15810.

Nama Kelompok

Kefas Ronaldo Immanuel - 20194920004

Putri Angelina Windjaya - 20194920010

Deskripsi Data

Berikut ini beberapa penjelasan terhadap variabel yang ada.

Variabel	Deskripsi
Age	Umur Pasien
Anaemia	Penurunan sel darah merah atau hemoglobin (boolean)
Creatinine_phosphokinase	Tipe chest pain, 1: typical angina, 2: atypical angina, 3: non-anginal pain, 4: asymptomatic
Diabetes	Jika pasien menderita diabetes (boolean)
Ejection_fraction	Persentase darah yang meninggalkan jantung pada setiap kontraksi (persentase)
high_blood_pressure	Jika pasien memiliki hipertensi (boolean)
platelets	Trombosit dalam darah (kiloplatelets / mL)
serum_creatinine	Tingkat kreatinin serum dalam darah (mg / dL)
serum_sodium	Tingkat natrium serum dalam darah (mEq / L)
sex	Wanita atau pria
smoking	Jika pasien merokok atau tidak (boolean)
time	Periode tindak lanjut (hari)
DEATH_EVENT	Jika pasien meninggal selama masa tindak lanjut (boolean)

About the Dataset

Penyakit kardiovaskular (CVD) adalah penyebab kematian nomor 1 secara global, merenggut nyawa sekitar 17,9 juta jiwa setiap tahun, yang menyumbang 31% dari semua kematian secara global. Gagal jantung adalah peristiwa umum yang disebabkan oleh CVD dan dataset ini berisi 12 fitur yang dapat digunakan untuk memprediksi kematian akibat gagal jantung.

Sebagian besar penyakit kardiovaskular dapat dicegah dengan mengatasi faktor risiko perilaku seperti penggunaan tembakau, diet dan obesitas yang tidak sehat, aktivitas fisik dan penggunaan alkohol yang berbahaya menggunakan strategi di seluruh populasi.

Orang dengan penyakit kardiovaskular atau yang berisiko kardiovaskular tinggi (karena adanya satu atau lebih faktor risiko seperti hipertensi, diabetes, hiperlipidaemia atau penyakit yang sudah mapan) memerlukan deteksi dan manajemen dini di mana model pembelajaran machine learning dapat sangat membantu.

Pemasukkan Data

Pertama, kami mengimport data yang kami gunakan pada survival analisis ini. Kami menggunakan data heart_failure

data <- read.csv("heart.csv")
datatable(data)

Pembersihan Data

Selanjutnya, kami melakukan Pembersihan data. Pembersihan data yang pertama adalah kami memeriksa struktur data dari ketiga belas variabel yang ada, dengan function str. Dengan memeriksa struktur data ini, maka kami mendapatkan jenis data dari setiap variabel.

str(data)

## 'data.frame':    299 obs. of  13 variables:
##  $ age                     : num  75 55 65 50 65 90 75 60 65 80 ...
##  $ anaemia                 : int  0 0 0 1 1 1 1 1 0 1 ...
##  $ creatinine_phosphokinase: int  582 7861 146 111 160 47 246 315 157 123 ...
##  $ diabetes                : int  0 0 0 0 1 0 0 1 0 0 ...
##  $ ejection_fraction       : int  20 38 20 20 20 40 15 60 65 35 ...
##  $ high_blood_pressure     : int  1 0 0 0 0 1 0 0 0 1 ...
##  $ platelets               : num  265000 263358 162000 210000 327000 ...
##  $ serum_creatinine        : num  1.9 1.1 1.3 1.9 2.7 2.1 1.2 1.1 1.5 9.4 ...
##  $ serum_sodium            : int  130 136 129 137 116 132 137 131 138 133 ...
##  $ sex                     : int  1 1 1 1 0 1 1 1 0 1 ...
##  $ smoking                 : int  0 0 1 0 0 1 0 1 0 1 ...
##  $ time                    : int  4 6 7 7 8 8 10 10 10 10 ...
##  $ DEATH_EVENT             : int  1 1 1 1 1 1 1 1 1 1 ...

Memeriksa data hilang

Lalu, dari data yang kami gunakan, kami akan memeriksa data hilang di dalamnya dengan function sebagai berikut.

sum(is.na(data))

## [1] 0

Didapatkan bahwa pemeriksaan data hilang memiliki hasil nol. Hal ini berarti tidak adanya data yang hilang.

Memeriksa data duplikat

Langkah selanjutnya, kami memeriksa data duplikat. Terkadang di suatu data, terdapat jumlah data yang duplikat. Oleh karena itu, kami akan melakukan pemeriksaan data duplikat dengan function sebagai berikut.

check.duplicate <- data.frame(
                   row_of_data        = data %>% nrow (),
                   row_of_unique.data = data %>% distinct() %>% nrow())
check.duplicate

##   row_of_data row_of_unique.data
## 1         299                299

Hasil dari pemeriksaan data duplikat adalah pada bagian row_of_data dan row_of_unique.data memiliki hasil yang sama, yaitu 299. Hal ini berarti tidak adanya data yang duplikat.

Eksplorasi Data

Setelah melakukan pembersihan data, kami melakukan eksplorasi data menggunakan visualisasi supaya dapat melihat lebih jelas bentuk datanya seperti apa.

pl_age <- ggplot(data, aes(x=age)) + geom_histogram(alpha = 0.5, fill = "orange", colour = "black")

pl_platelets <-  ggplot(data, aes(x=platelets)) + 
                 geom_histogram(alpha = 0.5, fill = "orange", colour = "black")

pl_anaemia <- data %>%
                  mutate(anaemia = factor(ifelse(anaemia == 1, "anaemic", "non-anaemic"))) %>%
                  group_by(anaemia) %>%
                  tally(name = "count") %>%
                  ggplot(aes(x = anaemia, y = count)) + 
                  geom_bar(stat = "Identity", fill = "orange", width = 0.2, colour = "black")

pl_creatinine <- ggplot(data, aes(x=creatinine_phosphokinase)) +
                 geom_histogram(alpha = 0.5, fill = "orange", colour = "black")

pl_serum_creatinine <- ggplot(data, aes(x=serum_creatinine)) +
                       geom_histogram(alpha = 0.5, fill = "orange", colour = "black")

pl_serum_sodium <- ggplot(data, aes(x=serum_sodium)) +
                   geom_histogram(alpha = 0.5, fill = "orange", colour = "black")

pl_ejection_fraction <- ggplot(data, aes(x=ejection_fraction)) +
                 geom_histogram(alpha = 0.5, fill = "orange", colour = "black")

pl_high_blood_pressure <-  data %>%
                  mutate(HBP = factor(ifelse(high_blood_pressure == 1, "high-bp", "non-high-bp"))) %>%
                  group_by(HBP) %>%
                  tally(name = "count") %>%
                  ggplot(aes(x = HBP, y = count)) + 
                  geom_bar(stat = "Identity", fill = "orange", width = 0.2, colour = "black")

pl_sex <- data %>%
                  mutate(Gender = factor(ifelse(sex == 0, "female", "male"))) %>%
                  group_by(Gender) %>%
                  tally(name = "count") %>%
                  ggplot(aes(x = Gender, y = count)) + 
                  geom_bar(stat = "Identity", fill = "orange", width = 0.2, colour = "black")

pl_smoking <- data %>%
                  mutate(Smoking = factor(ifelse(smoking == 0, "non-smoker", "smoker"))) %>%
                  group_by(Smoking) %>%
                  tally(name = "count") %>%
                  ggplot(aes(x = Smoking, y = count)) + 
                  geom_bar(stat = "Identity", fill = "orange", width = 0.2, colour = "black")

pl_diabetes <- data %>%
                  mutate(Diabetes = factor(ifelse(diabetes == 1, "diabetic", "non-diabetic"))) %>%
                  group_by(Diabetes) %>%
                  tally(name = "count") %>%
                  ggplot(aes(x = Diabetes, y = count)) + 
                  geom_bar(stat = "Identity", fill = "orange", width = 0.2, colour = "black")

pl_deathevent <- data %>%
                  mutate(DeathEvent = factor(ifelse(DEATH_EVENT == 1, "death", "censored"))) %>%
                  group_by(DeathEvent) %>%
                  tally(name = "count") %>%
                  ggplot(aes(x = DeathEvent, y = count)) + 
                  geom_bar(stat = "Identity", fill = "orange", width = 0.2, colour = "black")

figure_plot <- ggarrange(pl_age, pl_creatinine, pl_ejection_fraction, pl_platelets, pl_serum_sodium, pl_serum_creatinine, pl_anaemia, pl_deathevent, pl_diabetes, pl_high_blood_pressure, pl_sex, pl_smoking,
                    ncol = 3, nrow = 2)

figure_plot

## $`1`

## 
## $`2`

## 
## attr(,"class")
## [1] "list"      "ggarrange"

Pada diagram-diagram di atas menunjukkan bahwa variabel yang numerik sudah memenuhi pendekatan normal/ distribusi normal. Lalu, untuk diagram-diagram variabel kategorik hanya menampilkan jumlah dari setiap elemen kategori nya.

Dapat dilihat bahwa ada variabel DeathEvent, yang akan digunakan untuk mengetahui apakah pasien itu sensor atau meninggal. Berikut ini merupakan Time-to-Event chart dari 10 pasien acak.

Time to Event chart

set.seed(456)
plot_censoring <- data %>% group_by(DEATH_EVENT) %>% sample_n(5) %>% ungroup() %>% select(time, DEATH_EVENT)

plot_censoring %>%
    mutate(
        time_start = 0, 
        pasient = factor(c(1:nrow(plot_censoring))),
        death_event = factor(ifelse(DEATH_EVENT == 1, "death", "censored"))
    ) %>%
    pivot_longer(
        cols = c(time, time_start),
        names_to = "source",
        values_to = "time"
    ) %>%
    ggplot(aes(x = time, y = pasient, group = factor(pasient))) + 
    geom_bar(stat = "Identity", aes(fill = death_event), colour = "black", width = 0.3) +
    ggtitle("Kasus sampel dari dataset- Waktu hingga death_event")

Plot di atas menunjukkan 10 pasien acak, dengan 5 pasien meninggal dan 5 pasien sensor (dikeluarkan dari penelitian sebelum kematian dapat diamati). Fakta bahwa beberapa pasien sensor tidak berarti bahwa data mereka tidak ada. Kita masih tahu bahwa mereka tidak meninggal, setidaknya sampai saat mereka sensor.

Tujuan dilakukan analisis ini adalah untuk mengetahui variabel manakah yang signifikan memengaruhi peristiwa pada data ini, yaitu dari waktu hingga meninggalnya pasien gagal jantung.

Kaplan-Meier Model

Pada analisis ini akan digunakan Model Kaplan-Meier. Kaplan-Meier Analysis adalah uji statistik yang memberikan perkiraan fungsi survival sebenarnya dari suatu populasi, perkiraan menjadi lebih baik dengan peningkatan ukuran sampel. Analisis ini dapat menangani sensor dengan kuat, dan dapat berasal dari hazard function menggunakan Estimasi likelihood Maksimum. Analisis ini dapat digunakan untuk perbandingan sederhana tingkat survival antar kelompok (Misalnya, tingkat survival antara perokok dan non-perokok).

Kaplan-Meier curve based on death event

Model Kaplan-Meier menyediakan tabel probabilitas survival kumulatif berdasarkan waktu atau bisa disebut juga Life Table.

set.seed(456)
km_model <- survfit(Surv(time, DEATH_EVENT) ~ 1, data = data)
summary(km_model, times = seq(from = 0, to = 290, by = 30))

## Call: survfit(formula = Surv(time, DEATH_EVENT) ~ 1, data = data)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0    299       0    1.000  0.0000        1.000        1.000
##    30    264      35    0.882  0.0187        0.846        0.920
##    60    239      19    0.817  0.0225        0.774        0.863
##    90    189      15    0.763  0.0250        0.715        0.813
##   120    145       7    0.730  0.0268        0.680        0.785
##   150    118       5    0.703  0.0285        0.649        0.761
##   180    106       8    0.654  0.0313        0.596        0.719
##   210     62       4    0.622  0.0337        0.559        0.692
##   240     34       2    0.594  0.0378        0.524        0.673
##   270      6       1    0.576  0.0407        0.501        0.661

ggsurvplot(km_model, data = data, risk.table = T,
           break.time.by = 50, size = 0.3, tables.height = 0.4, ggtheme = theme_bw())

Pada kurva Kaplan-Meier di atas menunjukkan survival probability secara keseluruhan tanpa memikirkan kategori tertentu.Pada kurva, peristiwa penyensoran dilambangkan dengan tanda plus “+”. Dapat dilihat bahwa pada hari ke-0 survival probabilitynya adalah 1 karena belum terjadi peristiwa meninggal. Lalu survival probability terus turun seiring berjalannya waktu menjadi 0.576 (dilihat dari Life Table).

Apabila ingin melihat kurva Kaplan-Meier dengan memperhatikan beberapa kategori, berikut ada beberapa kurva Kaplan-Meier yang menampilkan hal tersebut.

Kaplan-Meier curve based on high blood pressure

# Kaplan-Meier curve based on high blood pressure
km_model <- data %>%
    mutate(
        high_blood_pressure =  factor(ifelse(high_blood_pressure == 1, "high-bp", "non-high-bp"))
    ) %>%
    survfit(Surv(time, DEATH_EVENT) ~ high_blood_pressure, data = .)
summary(km_model, times = seq(from = 0, to = 290, by = 30))

## Call: survfit(formula = Surv(time, DEATH_EVENT) ~ high_blood_pressure, 
##     data = .)
## 
##                 high_blood_pressure=high-bp 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0    105       0    1.000  0.0000        1.000        1.000
##    30     88      16    0.846  0.0355        0.779        0.918
##    60     75      13    0.718  0.0444        0.636        0.811
##    90     57       2    0.696  0.0458        0.612        0.792
##   120     39       3    0.651  0.0496        0.561        0.756
##   150     32       1    0.633  0.0514        0.540        0.742
##   180     29       2    0.592  0.0557        0.492        0.712
##   210     11       2    0.529  0.0652        0.416        0.674
##   240      3       0    0.529  0.0652        0.416        0.674
##   270      1       0    0.529  0.0652        0.416        0.674
## 
##                 high_blood_pressure=non-high-bp 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0    194       0    1.000  0.0000        1.000        1.000
##    30    176      19    0.902  0.0214        0.861        0.945
##    60    164       6    0.870  0.0242        0.824        0.919
##    90    132      13    0.798  0.0294        0.743        0.858
##   120    106       4    0.771  0.0313        0.712        0.835
##   150     86       4    0.739  0.0338        0.676        0.809
##   180     77       6    0.687  0.0376        0.617        0.765
##   210     51       2    0.665  0.0396        0.592        0.747
##   240     31       2    0.630  0.0444        0.549        0.724
##   270      5       1    0.609  0.0476        0.523        0.710

ggsurvplot(km_model, data = data, risk.table = T, pval = T,
           break.time.by = 50, size = 0.3, tables.height = 0.4, ggtheme = theme_bw())

Kurva Kaplan-Meier berdasarkan tekanan darah tinggi (high blood pressure) menunjukkan bahwa kurva survival pasien dengan non-tekanan darah tinggi berada di atas kurva survival pasien dengan tekanan darah tinggi, atau pasien non-tekanan darah tinggi lebih survive dibanding pasien dengan tekanan darah tinggi. Dengan melihat Life Table, survival probability non tekanan darah tinggi mencapai 0.609, sedangkan survival probability tekanan darah tinggi mencapai 0.529. Kurva di atas juga menunjukkan bahwa p-value = 0.036.

Kaplan-Meier curve based on diabetes

# Kaplan-Meier curve based on diabetes
km_model <- data %>%
    mutate(
        diabetes = factor(ifelse(diabetes == 1, "diabetic", "non-diabetic"))
    ) %>%
    survfit(Surv(time, DEATH_EVENT) ~ diabetes, data = .)
summary(km_model, times = seq(from = 0, to = 290, by = 30))

## Call: survfit(formula = Surv(time, DEATH_EVENT) ~ diabetes, data = .)
## 
##                 diabetes=diabetic 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0    125       0    1.000  0.0000        1.000        1.000
##    30    112      15    0.879  0.0293        0.823        0.938
##    60     97      10    0.796  0.0364        0.728        0.871
##    90     78       4    0.762  0.0386        0.690        0.842
##   120     63       1    0.752  0.0395        0.678        0.833
##   150     53       5    0.690  0.0449        0.607        0.784
##   180     46       3    0.650  0.0478        0.563        0.751
##   210     32       1    0.633  0.0497        0.542        0.738
##   240     18       1    0.602  0.0557        0.503        0.722
##   270      3       0    0.602  0.0557        0.503        0.722
## 
##                 diabetes=non-diabetic 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0    174       0    1.000  0.0000        1.000        1.000
##    30    152      20    0.885  0.0243        0.838        0.933
##    60    142       9    0.832  0.0285        0.778        0.890
##    90    111      11    0.763  0.0329        0.701        0.830
##   120     82       6    0.715  0.0362        0.648        0.790
##   150     65       0    0.715  0.0362        0.648        0.790
##   180     60       5    0.659  0.0411        0.584        0.745
##   210     30       3    0.615  0.0457        0.531        0.711
##   240     16       1    0.591  0.0497        0.501        0.697
##   270      3       1    0.552  0.0600        0.446        0.683

ggsurvplot(km_model, data = data, risk.table = T, pval = T,
           break.time.by = 50, size = 0.3, tables.height = 0.4, ggtheme = theme_bw())

Kurva Kaplan-Meier berdasarkan diabetes menunjukkan bahwa kurva survival pasien dengan diabetes berada di atas kurva survival pasien dengan non-diabetes, atau pasien dengan diabetes lebih survive dibanding pasien non-diabetes. Dengan melihat Life Table, survival probability diabetes mencapai 0.602, sedangkan survival probability non-diabetes mencapai 0.552. Kurva di atas juga menunjukkan bahwa p-value = 0.84.

Kaplan-Meier curve based on gender

# Kaplan-Meier curve based on gender
km_model <- data %>%
    mutate(
        sex = factor(ifelse(sex == 0, "female", "male"))
    ) %>%
    survfit(Surv(time, DEATH_EVENT) ~ sex, data = .)
summary(km_model, times = seq(from = 0, to = 290, by = 30))

## Call: survfit(formula = Surv(time, DEATH_EVENT) ~ sex, data = .)
## 
##                 sex=female 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0    105       0    1.000  0.0000        1.000        1.000
##    30     94      10    0.903  0.0291        0.848        0.962
##    60     85       8    0.825  0.0375        0.754        0.902
##    90     69       6    0.762  0.0426        0.683        0.850
##   120     50       3    0.723  0.0461        0.638        0.819
##   150     43       2    0.693  0.0488        0.603        0.795
##   180     38       2    0.660  0.0516        0.567        0.770
##   210     23       3    0.594  0.0590        0.489        0.722
##   240     13       0    0.594  0.0590        0.489        0.722
##   270      2       0    0.594  0.0590        0.489        0.722
## 
##                 sex=male 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0    194       0    1.000  0.0000        1.000        1.000
##    30    170      25    0.871  0.0241        0.825        0.919
##    60    154      11    0.813  0.0281        0.760        0.870
##    90    120       9    0.763  0.0310        0.705        0.826
##   120     95       4    0.735  0.0329        0.673        0.802
##   150     75       3    0.708  0.0351        0.643        0.780
##   180     68       6    0.650  0.0394        0.577        0.732
##   210     39       1    0.638  0.0406        0.563        0.722
##   240     21       2    0.592  0.0494        0.502        0.697
##   270      4       1    0.562  0.0551        0.464        0.681

ggsurvplot(km_model, data = data, risk.table = T, pval = T,
           break.time.by = 50, size = 0.3, tables.height = 0.4, ggtheme = theme_bw())

Kurva Kaplan-Meier berdasarkan jenis kelamin menunjukkan bahwa kurva survival pasien wanita berada di atas kurva survival pasien pria, atau pasien wanita lebih survive dibanding pasien pria. Dengan melihat Life Table, survival probability wanita mencapai 0.594, sedangkan survival probability pria mencapai 0.562. Kurva di atas juga menunjukkan bahwa p-value = 0.95.

Kaplan-Meier curve based on presence/absence of smoking

# Kaplan-Meier curve based on presence/absence of smoking
km_model <- data %>%
    mutate(
        smoking = factor(ifelse(smoking == 0, "non-smoker", "smoker"))
    ) %>%
    survfit(Surv(time, DEATH_EVENT) ~ smoking, data = .)
summary(km_model, times = seq(from = 0, to = 290, by = 30))

## Call: survfit(formula = Surv(time, DEATH_EVENT) ~ smoking, data = .)
## 
##                 smoking=non-smoker 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0    203       0    1.000  0.0000        1.000        1.000
##    30    178      25    0.876  0.0233        0.831        0.923
##    60    163      11    0.820  0.0272        0.768        0.875
##    90    131       8    0.777  0.0297        0.721        0.838
##   120     99       6    0.737  0.0325        0.676        0.803
##   150     84       4    0.705  0.0347        0.641        0.777
##   180     75       6    0.654  0.0380        0.583        0.733
##   210     44       4    0.608  0.0416        0.532        0.696
##   240     22       2    0.569  0.0478        0.482        0.670
##   270      2       0    0.569  0.0478        0.482        0.670
## 
##                 smoking=smoker 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     0     96       0    1.000  0.0000        1.000        1.000
##    30     86      10    0.896  0.0312        0.837        0.959
##    60     76       8    0.812  0.0399        0.737        0.894
##    90     58       7    0.731  0.0462        0.646        0.828
##   120     46       1    0.717  0.0475        0.630        0.816
##   150     34       1    0.696  0.0506        0.604        0.802
##   180     31       2    0.654  0.0556        0.553        0.772
##   210     18       0    0.654  0.0556        0.553        0.772
##   240     12       0    0.654  0.0556        0.553        0.772
##   270      4       1    0.594  0.0759        0.463        0.763

ggsurvplot(km_model, data = data, risk.table = T, pval = T,
           break.time.by = 50, size = 0.3, tables.height = 0.4, ggtheme = theme_bw())

Kurva Kaplan-Meier berdasarkan merokok menunjukkan bahwa kurva survival pasien yang merokok berada di atas kurva survival pasien yang tidak merokok, atau pasien yang merokok lebih survive dibanding pasien yang tidak merokok. Dengan melihat Life Table, survival probability merokok mencapai 0.594, sedangkan survival probability non-merokok mencapai 0.569. Kurva di atas juga menunjukkan bahwa p-value = 0.96.

Kesimpulan

Berdasarkan analisis yang telah dilakukan ada beberapa hasil yang cukup mengejutkan yaitu pada bagian kurva Kaplan-Meier diabetes dan merokok. Dikatakan bahwa pasien diabetes lebih survive dibanding pasien non-diabetes. Hal ini terjadi karena sebenarnya pasien non-diabetes memiliki kemungkinan bertahan hidup yang tinggi hanya diawal waktu, tetapi kelangsungan hidupnya cenderung rendah untuk jangka waktu yang lebih lama. Hal ini juga terjadi pada pasien non-merokok. Pasien non-merokok cenderung memiliki kemungkinan bertahan hidup yang lebih tinggi pada awalnya tetapi kelangsungan hidup yang lebih rendah untuk jangka waktu yang lebih lama, dibandingkan pasien merokok.

Dari hasil p-value setiap kurva Kaplan-Meier menunjukkan bahwa semua variabel tersebut signifikan, hanya saja yang paling signifikan memengaruhi pasien gagal jantung adalah pasien yang merokok atau tidak.