Survival Analysis of Non-alcoholic Fatty Liver Disease (NAFLD)
Novia Permatasari
05/22/2020
Analisis ini disusun untuk menyelesaikan Workshop Survival Analysis with R Studio Cloud oleh Data Science Indonesia (DSI).
Install and Load package
library(survival)
library(survminer)
library(concordance)
options(scipen=999)
library(dplyr)Dataset
Data yang digunakan adalah dataset nafld1 yang ada pada package survival. Dataset ini merupakan hasil studi penyakit perlemakan hati pada pasien non-alkoholik (Non-alcohol Fatty Liver Disease - NAFLD).
df <- nafld1
# 10 data teratas
df %>% head(10)## id age male weight height bmi case.id futime status
## 3631 1 57 0 60.0 163 22.69094 10630 6261 0
## 8458 2 67 0 70.4 168 24.88403 14817 624 0
## 6298 3 53 1 105.8 186 30.45354 3 1783 0
## 15398 4 56 1 109.3 170 37.83010 6628 3143 0
## 13261 5 68 1 NA NA NA 1871 1836 1
## 14423 6 39 0 63.9 155 26.61559 7144 1581 0
## 9518 7 49 0 66.2 161 25.51934 7507 3109 0
## 15366 8 30 1 NA 180 NA 13417 1339 0
## 8827 9 47 1 110.8 188 31.42360 9 1671 0
## 688 10 79 0 56.6 155 23.71853 13518 2239 1# gambaran data
df %>% glimpse()## Rows: 17,549
## Columns: 9
## $ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ age <dbl> 57, 67, 53, 56, 68, 39, 49, 30, 47, 79, 47, 59, 54, 41, 51, 7…
## $ male <dbl> 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1…
## $ weight <dbl> 60.0, 70.4, 105.8, 109.3, NA, 63.9, 66.2, NA, 110.8, 56.6, NA…
## $ height <dbl> 163, 168, 186, 170, NA, 155, 161, 180, 188, 155, NA, 182, 167…
## $ bmi <dbl> 22.69094, 24.88403, 30.45354, 37.83010, NA, 26.61559, 25.5193…
## $ case.id <int> 10630, 14817, 3, 6628, 1871, 7144, 7507, 13417, 9, 13518, 970…
## $ futime <dbl> 6261, 624, 1783, 3143, 1836, 1581, 3109, 1339, 1671, 2239, 48…
## $ status <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…Data ini terdiri atas 17549 observasi dan 9 variabel sebagai berikut:
- id : subject identifier
- age : age at entry to the study
- male : 0=female, 1=male
- weight : weight in kg
- height : height in cm
- bmi : body mass index
- case.id : the id of the NAFLD case to whom this subject is matched
- futime : time to death or last follow-up
- status : 0= alive at last follow-up, 1=dead
Analisis
Berikut urutan analisis yang dilakukan:
- Melihat Peluang Survival Pasien dengan Kaplan-Meier
- Melihat perbedaan tingkat survival pada pasien NAFDL dengan kategori yang berbeda dengan Uji Log Rank
- Melakukan pemodelan dengan beberapa model Non-Parametric dan model Parametrik
- Menentukan model terbaik berdasarkan nilai Concordance Indeks
Melihat Peluang Survival Pasien dengan Kaplan-Meier
fit <- survfit(Surv(futime, status)~1, data = df)
surv_table <- data.frame(time=fit$time,
n_risk=fit$n.risk,
n_event=fit$n.event,
n_censor=fit$n.censor,
sur_prob=fit$surv)
# Table Kaplan-Meier Estimator
surv_table %>% head()## time n_risk n_event n_censor sur_prob
## 1 7 17549 0 3 1.000000
## 2 9 17546 0 3 1.000000
## 3 10 17543 1 4 0.999943
## 4 11 17538 1 2 0.999886
## 5 12 17535 0 1 0.999886
## 6 13 17534 0 2 0.999886# Kaplan Meier Curve
ggsurvplot(fit,
surv.median.line = "hv",
ggtheme = theme_bw(),
color="blue",
legend="none",
xlab = "Time in Days") +
ggtitle("Kaplan Meier Curve Lung Dataset")
Berdasarkan kurva diaras, dapat dilihat bahwa nilai survival probability cukup tinggi, yang menunjukkan tingkat kematian akibat penyakit nafld rendah.
min(surv_table$sur_prob)## [1] 0.6929272min((surv_table%>%filter(sur_prob<0.6929272))$time)## [1] 6966max((surv_table%>%filter(sur_prob<0.6929272))$time)## [1] 7268Peluang seseorang dapat survive terkecil adalah 69.29%, pada hari ke-6966 hingga ke-7268.
Melihat perbedaan tingkat survival pada pasien NAFLD berdasarkan gender dengan Uji Log Rank
fit2 <- survfit(Surv(futime, status)~male, data = df)
ggsurvplot(fit2,
surv.median.line = "hv",
pval = TRUE,
ggtheme = theme_bw(),
legend="right",
xlab = "Time in Days") +
ggtitle("Kaplan Meier Curve NAFLD Dataset")
Nilai p-value < alfa=0.05 menunjukkan bahwa terdapat perbedaan sigfinikan pada peluang survive penderita NADLF laki-laki dan perempuan. Pada waktu (hari) kurang dari 6700-6800, perempuan cenderung memiliki kemampuan survive yang lebih tinggi daripada laki-laki. Sedangkan pada hari yang lebih dari itu, kemampuan survive laki-laki cenderung lebih tinggi.
Melakukan pemodelan dengan beberapa model Non-Parametric dan model Parametrik
Model Semi-Parametrix: CoX Proportional Hazard
df$male <- as.factor(df$male)
model_coxph <- coxph(Surv(futime, status)~age + male + weight + height + bmi, data = df)
model_coxph %>% summary()## Call:
## coxph(formula = Surv(futime, status) ~ age + male + weight +
## height + bmi, data = df)
##
## n= 12588, number of events= 1018
## (4961 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.098330 1.103327 0.002771 35.483 < 0.0000000000000002 ***
## male1 0.550169 1.733547 0.091513 6.012 0.00000000183 ***
## weight 0.009734 1.009781 0.012206 0.797 0.425
## height -0.022267 0.977979 0.013629 -1.634 0.102
## bmi -0.011134 0.988927 0.033234 -0.335 0.738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.1033 0.9063 1.0974 1.109
## male1 1.7335 0.5769 1.4489 2.074
## weight 1.0098 0.9903 0.9859 1.034
## height 0.9780 1.0225 0.9522 1.004
## bmi 0.9889 1.0112 0.9266 1.055
##
## Concordance= 0.828 (se = 0.007 )
## Likelihood ratio test= 1758 on 5 df, p=<0.0000000000000002
## Wald test = 1489 on 5 df, p=<0.0000000000000002
## Score (logrank) test = 1729 on 5 df, p=<0.0000000000000002Uji Simultan Likelihood Ratio Test menghasilkan nilai pvalue kurang dari alfa =5% yang menunjukkan bahwa minimal terdapat satu variabel yang signifikan dan model telah valid. Terdapat dua variabel yang secara signifikan yang berpengaruh, yaitu age dan male1, sedangkan variabel weight, height, dan bmi tidak signifikan. Seperti kita ketahui bahwa bmi merupakan turunan dari weight dan height, sehingga kedua variabel dapat kita hilangkan.
model_coxph2 <- coxph(Surv(futime, status)~age + male + bmi, data = df)
model_coxph2 %>% summary()## Call:
## coxph(formula = Surv(futime, status) ~ age + male + bmi, data = df)
##
## n= 12588, number of events= 1018
## (4961 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.100601 1.105835 0.002650 37.967 < 0.0000000000000002 ***
## male1 0.366140 1.442158 0.062849 5.826 0.00000000569 ***
## bmi 0.017022 1.017168 0.004945 3.443 0.000576 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.106 0.9043 1.100 1.112
## male1 1.442 0.6934 1.275 1.631
## bmi 1.017 0.9831 1.007 1.027
##
## Concordance= 0.827 (se = 0.007 )
## Likelihood ratio test= 1750 on 3 df, p=<0.0000000000000002
## Wald test = 1472 on 3 df, p=<0.0000000000000002
## Score (logrank) test = 1699 on 3 df, p=<0.0000000000000002Uji Simultan Likelihood Ratio Test menghasilkan nilai pvalue kurang dari alfa =5% yang menunjukkan bahwa minimal terdapat satu variabel yang signifikan dan model telah valid. Uji partial juga menunjukkan bahwa semua variabel independen telah signifikan. Nilai exp(male1) menunjukkan bahwa resiko meninggal pasien laki-laki lebih tinggi 1,442 kali lebih cepat daripada perempuan
Exponential Model
model_exp <- survreg(Surv(futime, status)~age + male + bmi, data = df, dist="exponential")
summary(model_exp)##
## Call:
## survreg(formula = Surv(futime, status) ~ age + male + bmi, data = df,
## dist = "exponential")
## Value Std. Error z p
## (Intercept) 16.77488 0.25084 66.88 < 0.0000000000000002
## age -0.09497 0.00253 -37.50 < 0.0000000000000002
## male1 -0.33232 0.06271 -5.30 0.00000012
## bmi -0.01555 0.00490 -3.18 0.0015
##
## Scale fixed at 1
##
## Exponential distribution
## Loglik(model)= -10674.1 Loglik(intercept only)= -11502.9
## Chisq= 1657.73 on 3 degrees of freedom, p= 0
## Number of Newton-Raphson Iterations: 7
## n=12588 (4961 observations deleted due to missingness)Uji Simultan menghasilkan nilai pvalue kurang dari alfa =5% yang menunjukkan bahwa minimal terdapat satu variabel yang signifikan dan model telah valid. Uji partial juga menunjukkan bahwa semua variabel independen telah signifikan.
exp(model_exp$coefficients)## (Intercept) age male1 bmi
## 19285813.4964394 0.9094037 0.7172556 0.9845663Jenis kelamin laki-laki memiliki peluang survive lebih kecil dibandingkan dengan perempuan.
Menentukan model terbaik berdasarkan nilai Concordance Indeks
data.frame(CoxPh = model_coxph2$concordance[6], Exp = concordance(model_exp)$concordance)## CoxPh Exp
## concordance 0.8269613 0.8268536Berdasarkan nilai corcondance tersebut, dapat dilihat bahwa model CoxPH sedikit lebih baik dari model regresi exponensial karena nilai concordance yang lebih tinggi.
Kesimpulan
Model yang terbaik adalah model semi-parametrik Cox-PH karena memiliki nilai concordance paling tinggi yaitu sebesar 0.8269. Variabel yang mempengaruhi peluang survive pasien adalah jenis kelamin, umur, dan bmi. Semakin tinggi umur dan bmi seorang pasien, maka resiko meninggal pasien tersebut juga semakin tinggi. Selain itu, pasien laki-laki memiliki resiko meninggal lebih tinggi pasien perempuan.