calculate_clif_c_score <- function(tb_admit, peak_creatinine, he_grade_admission, inr_admit, map, pressor, fio2_admission, spo2_admission, age_admission, wbc_admit) {
clif_score <- 0
# TB Admit Score
clif_score <- clif_score + case_when(
is.na(tb_admit) ~ 0,
tb_admit < 6 ~ 1,
tb_admit >= 6 & tb_admit < 12 ~ 2,
tb_admit >= 12 ~ 3,
TRUE ~ 0
)
# Peak Creatinine Score
clif_score <- clif_score + case_when(
is.na(peak_creatinine) ~ 0,
peak_creatinine < 2 ~ 1,
peak_creatinine >= 2 & peak_creatinine < 3.5 ~ 2,
peak_creatinine >= 3.5 ~ 3,
TRUE ~ 0
)
# HE Grade Admission Score
clif_score <- clif_score + case_when(
is.na(he_grade_admission) ~ 0,
he_grade_admission == 1 ~ 1,
he_grade_admission %in% c(2, 3) ~ 2,
he_grade_admission %in% c(4, 5) ~ 3,
TRUE ~ 0
)
# INR Admit Score
clif_score <- clif_score + case_when(
is.na(inr_admit) ~ 0,
inr_admit < 2 ~ 1,
inr_admit >= 2 & inr_admit < 2.5 ~ 2,
inr_admit >= 2.5 ~ 3,
TRUE ~ 0
)
# MAP and Pressor Score
clif_score <- clif_score + case_when(
is.na(map) ~ 0,
map >= 70 & pressor != 1 ~ 1,
map < 70 & pressor != 1 ~ 2,
!is.na(map) & pressor == 1 ~ 3,
TRUE ~ 0
)
# FIO2 and SPO2 Admission Score
if (!is.na(fio2_admission)) {
ratio <- spo2_admission / fio2_admission
clif_score <- clif_score + case_when(
ratio > 357 ~ 1,
ratio > 214 & ratio <= 357 ~ 2,
ratio <= 214 ~ 3,
TRUE ~ 0
)
}
# Calculate CLIF-C Score
CLIF_C_Score <- 10 * (0.33 * clif_score + 0.04 * age_admission + 0.63 * log(wbc_admit) - 2)
return(CLIF_C_Score)
}
final_master <- final_master %>%
mutate(across(
c(map_terli_day0_time0, map_terli_day0_time1, map_terli_day0_time2, map_terli_day0_time3),
~ na_if(as.numeric(.), 0)
))
map_terli_day0 <- final_master %>%
group_by(subjectid) %>%
summarise(
map_terli_day0 = median(
c(map_terli_day0_time0, map_terli_day0_time1, map_terli_day0_time2, map_terli_day0_time3),
na.rm = TRUE
),
.groups = "drop"
)
final_master <- final_master %>% mutate(across(all_of(c("bili_terli_day0", "scr_terli_day0", "he_terli_day0","alb_terli_day0",
"inr_terli_day0", "map_terli_day0_time0", "norepi_terli_day0",
"fio2_terli_day0", "spo2_terli_day0", "age", "wbc_terli_day0")),
~ ifelse(is.na(.), median(., na.rm = TRUE), .)))
final_master <- final_master %>% merge(map_terli_day0,by="subjectid",all.x = T)
final_master <- final_master %>% mutate(
clif_c_score_day0= mapply(calculate_clif_c_score, bili_terli_day0, scr_terli_day0, he_terli_day0,
inr_terli_day0, map_terli_day0, norepi_terli_day0, fio2_terli_day0, spo2_terli_day0,
age, wbc_terli_day0))
final_master <- final_master %>%
mutate(dc_days_to_death=ifelse(disposition == 5,0,dc_days_to_death)) %>%
mutate(
death_hosp = if_else(disposition == 5, 1L, 0L)
) %>%
mutate(
days_to_dc_init = days_to_discharge - days_to_terli_initiation_day0,
days_to_los_init = days_to_lost_to_followup + days_to_discharge - days_to_terli_initiation_day0,
days_to_lastencounter_init = dc_days_lastencounter + days_to_discharge - days_to_terli_initiation_day0,
days_to_after_dc_death_init= dc_days_to_death + days_to_discharge - days_to_terli_initiation_day0,
death_status_check = case_when(
disposition == 5 | dc_death_90days == 1 ~ 1,
TRUE ~ 0),
time_to_death_check = case_when(
disposition == 5 ~ days_to_dc_init,
dc_death_90days == 2 ~ days_to_los_init,
dc_death_90days == 1 ~ days_to_after_dc_death_init,
dc_death_90days == 0 & is.na(days_to_lastencounter_init) ~ days_to_los_init,
dc_death_90days == 0 & !is.na(days_to_lastencounter_init) ~ days_to_lastencounter_init,
TRUE ~ NA_real_),
death_status_90days_init=ifelse(time_to_death_check<=90 & death_status_check==1,1,0),
time_to_death_90days_init=ifelse(time_to_death_check<=90 ,time_to_death_check, 90)) %>%
mutate(
# ---- Median impute each needed column ----
dc_scr = ifelse(is.na(dc_scr), median(dc_scr, na.rm = TRUE), dc_scr),
dc_tbili = ifelse(is.na(dc_tbili), median(dc_tbili, na.rm = TRUE), dc_tbili),
dc_inr = ifelse(is.na(dc_inr), median(dc_inr, na.rm = TRUE), dc_inr),
dc_sna = ifelse(is.na(dc_sna), median(dc_sna, na.rm = TRUE), dc_sna),
dc_albumin = ifelse(is.na(dc_albumin), median(dc_albumin, na.rm = TRUE), dc_albumin)) %>%
mutate(
dc_meld = round(9.57 * log(dc_scr) + 3.78 * log(dc_tbili) + 11.20 * log(dc_inr) + 6.43),
dc_meld_na = round(dc_meld + 1.32 * (137 - pmin(pmax(dc_sna, 125), 137)) - (0.033 * dc_meld * (137 - pmin(pmax(dc_sna, 125), 137)))),
dc_meld_na = ifelse(dc_meld_na > 11, dc_meld_na, NA), # MELDNa applied only if MELD > 11
sex_male=as.numeric(sex_male),
female_factor = if_else(sex_male == 0, 1.33, 0), # Assuming a column 'sex' is present, adjust accordingly
dc_meld_3 = round(
female_factor +
4.56 * log(dc_tbili) +
0.82 * (137 - dc_sna) -
0.24 * (137 - dc_sna) * log(dc_tbili) +
9.09 * log(dc_inr) +
11.14 * log(dc_scr) +
1.85 * (3.5 - dc_albumin) -
1.83 * (3.5 - dc_albumin) * log(dc_scr) +
6
)
) %>%
mutate(day0_meld = round(9.57 * log(scr_terli_day0) + 3.78 * log(bili_terli_day0) + 11.20 * log(inr_terli_day0) + 6.43),
day0_meld_na = round(day0_meld + 1.32 * (137 - pmin(pmax(na_terli_day0, 125), 137)) - (0.033 * day0_meld * (137 - pmin(pmax(na_terli_day0, 125), 137)))),
day0_meld_na = ifelse(day0_meld_na > 11, day0_meld_na, NA), # MELDNa applied only if MELD > 11
sex_male=as.numeric(sex_male),
female_factor = if_else(sex_male == 0, 1.33, 0), # Assuming a column 'sex' is present, adjust accordingly
day0_meld_3 = round(
female_factor +
4.56 * log(bili_terli_day0) +
0.82 * (137 - na_terli_day0) -
0.24 * (137 - na_terli_day0) * log(bili_terli_day0) +
9.09 * log(inr_terli_day0) +
11.14 * log(scr_terli_day0) +
1.85 * (3.5 - alb_terli_day0) -
1.83 * (3.5 - alb_terli_day0) * log(scr_terli_day0) +
6
))
test <- final_master %>%
select(
scr_terli_day0,
bili_terli_day0,
inr_terli_day0,
na_terli_day0,
alb_terli_day0,
female_factor,
day0_meld,
day0_meld_na,
day0_meld_3
)
cox day0_meld
## Call:
## coxph(formula = Surv(time_to_death_90days_init, death_status_90days_init) ~
## day0_meld, data = final_master, x = TRUE)
##
## n= 243, number of events= 102
##
## coef exp(coef) se(coef) z Pr(>|z|)
## day0_meld 0.05121 1.05254 0.01312 3.903 9.5e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## day0_meld 1.053 0.9501 1.026 1.08
##
## Concordance= 0.626 (se = 0.028 )
## Likelihood ratio test= 14.65 on 1 df, p=1e-04
## Wald test = 15.23 on 1 df, p=1e-04
## Score (logrank) test = 15.2 on 1 df, p=1e-04
cox day0_meld_3
## Call:
## coxph(formula = Surv(time_to_death_90days_init, death_status_90days_init) ~
## day0_meld_3, data = final_master, x = TRUE)
##
## n= 243, number of events= 102
##
## coef exp(coef) se(coef) z Pr(>|z|)
## day0_meld_3 0.05447 1.05598 0.01330 4.096 4.2e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## day0_meld_3 1.056 0.947 1.029 1.084
##
## Concordance= 0.634 (se = 0.027 )
## Likelihood ratio test= 16.88 on 1 df, p=4e-05
## Wald test = 16.78 on 1 df, p=4e-05
## Score (logrank) test = 16.68 on 1 df, p=4e-05
cox day0_meld_na
## Call:
## coxph(formula = Surv(time_to_death_90days_init, death_status_90days_init) ~
## day0_meld_na, data = final_master, x = TRUE)
##
## n= 243, number of events= 102
##
## coef exp(coef) se(coef) z Pr(>|z|)
## day0_meld_na 0.06372 1.06580 0.01505 4.233 2.3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## day0_meld_na 1.066 0.9383 1.035 1.098
##
## Concordance= 0.643 (se = 0.027 )
## Likelihood ratio test= 17.39 on 1 df, p=3e-05
## Wald test = 17.92 on 1 df, p=2e-05
## Score (logrank) test = 17.76 on 1 df, p=3e-05
cox clif_c_score_new
## Call:
## coxph(formula = Surv(time_to_death_90days_init, death_status_90days_init) ~
## clif_c_score_new, data = final_master, x = TRUE)
##
## n= 243, number of events= 102
##
## coef exp(coef) se(coef) z Pr(>|z|)
## clif_c_score_new 0.05691 1.05856 0.01229 4.631 3.64e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## clif_c_score_new 1.059 0.9447 1.033 1.084
##
## Concordance= 0.632 (se = 0.028 )
## Likelihood ratio test= 21.12 on 1 df, p=4e-06
## Wald test = 21.45 on 1 df, p=4e-06
## Score (logrank) test = 21.71 on 1 df, p=3e-06
ROC_day0_meld
## Time-dependent-Roc curve estimated using IPCW (n=243, without competing risks).
## Cases Survivors Censored AUC (%) se
## t=7 11 227 5 57.37 8.82
## t=15 40 193 10 64.40 4.54
## t=30 72 157 14 68.79 3.55
## t=60 93 135 15 66.01 3.56
## t=89 102 119 22 63.35 3.72
##
## Method used for estimating IPCW:marginal
##
## Total computation time : 0.06 secs.
ROC_day0_meld_3
## Time-dependent-Roc curve estimated using IPCW (n=243, without competing risks).
## Cases Survivors Censored AUC (%) se
## t=7 11 227 5 61.30 8.47
## t=15 40 193 10 65.64 4.49
## t=30 72 157 14 69.09 3.59
## t=60 93 135 15 66.84 3.55
## t=89 102 119 22 63.95 3.72
##
## Method used for estimating IPCW:marginal
##
## Total computation time : 0.06 secs.
ROC_day0_meld_na
## Time-dependent-Roc curve estimated using IPCW (n=243, without competing risks).
## Cases Survivors Censored AUC (%) se
## t=7 11 227 5 60.92 8.70
## t=15 40 193 10 67.06 4.51
## t=30 72 157 14 70.14 3.53
## t=60 93 135 15 67.56 3.53
## t=89 102 119 22 64.91 3.71
##
## Method used for estimating IPCW:marginal
##
## Total computation time : 0.06 secs.
ROC_clif_c_score_new
## Time-dependent-Roc curve estimated using IPCW (n=243, without competing risks).
## Cases Survivors Censored AUC (%) se
## t=7 11 227 5 71.39 6.94
## t=15 40 193 10 62.94 5.00
## t=30 72 157 14 67.49 3.90
## t=60 93 135 15 65.78 3.67
## t=89 102 119 22 65.60 3.67
##
## Method used for estimating IPCW:marginal
##
## Total computation time : 0.07 secs.
ROC_day0_meld,ROC_day0_meld_3
## $p_values_AUC
## t=7 t=15 t=30 t=60 t=89
## Non-adjusted 0.3622644 0.5121694 0.8400498 0.5500238 0.6605852
## Adjusted 0.7608672 0.9083621 0.9992873 0.9326755 0.9784685
##
## $Cor
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.6551798 0.4504780 0.3492054 0.3249587
## [2,] 0.6551798 1.0000000 0.6951947 0.5187546 0.4840265
## [3,] 0.4504780 0.6951947 1.0000000 0.7581478 0.7038661
## [4,] 0.3492054 0.5187546 0.7581478 1.0000000 0.9299945
## [5,] 0.3249587 0.4840265 0.7038661 0.9299945 1.0000000
ROC_day0_meld_na,ROC_day0_meld_3
## $p_values_AUC
## t=7 t=15 t=30 t=60 t=89
## Non-adjusted 0.8839211 0.2451665 0.2766896 0.4377083 0.2906298
## Adjusted 0.9998795 0.6034371 0.6562911 0.8591194 0.6782232
##
## $Cor
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.5673338 0.3642811 0.2562463 0.2441758
## [2,] 0.5673338 1.0000000 0.6645594 0.4889851 0.4446916
## [3,] 0.3642811 0.6645594 1.0000000 0.7734885 0.6965963
## [4,] 0.2562463 0.4889851 0.7734885 1.0000000 0.9041193
## [5,] 0.2441758 0.4446916 0.6965963 0.9041193 1.0000000
ROC_clif_c_score_new,ROC_day0_meld_3
## $p_values_AUC
## t=7 t=15 t=30 t=60 t=89
## Non-adjusted 0.2432789 0.6578371 0.7194313 0.7913167 0.6829593
## Adjusted 0.5867714 0.9773633 0.9902902 0.9974349 0.9835655
##
## $Cor
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.0000000 0.3557890 0.2276747 0.2060306 0.1859800
## [2,] 0.3557890 1.0000000 0.7203923 0.6307713 0.5756174
## [3,] 0.2276747 0.7203923 1.0000000 0.8746429 0.8001932
## [4,] 0.2060306 0.6307713 0.8746429 1.0000000 0.9169460
## [5,] 0.1859800 0.5756174 0.8001932 0.9169460 1.0000000
Table 1 function
Table 1
Table death_status_90days_init==1
##
##
## Cell Contents
## |-------------------------|
## | N |
## |-------------------------|
##
##
## Total Observations in Table: 102
##
##
## | final_master_death$meld_na_cat3
## final_master_death$meld_3_cat3 | <25 | >35 | 25 - 35 | Row Total |
## -------------------------------|-----------|-----------|-----------|-----------|
## <25 | 2 | 0 | 1 | 3 |
## -------------------------------|-----------|-----------|-----------|-----------|
## >35 | 0 | 46 | 16 | 62 |
## -------------------------------|-----------|-----------|-----------|-----------|
## 25 - 35 | 1 | 0 | 36 | 37 |
## -------------------------------|-----------|-----------|-----------|-----------|
## Column Total | 3 | 46 | 53 | 102 |
## -------------------------------|-----------|-----------|-----------|-----------|
##
##
logistic regression
roc_meld
roc(final_master$hrs_responders_cat_2, final_master$pred_day0_meld, ci=TRUE,direction = "<")
## Setting levels: control = 0, case = 1
##
## Call:
## roc.default(response = final_master$hrs_responders_cat_2, predictor = final_master$pred_day0_meld, direction = "<", ci = TRUE)
##
## Data: final_master$pred_day0_meld in 109 controls (final_master$hrs_responders_cat_2 0) < 134 cases (final_master$hrs_responders_cat_2 1).
## Area under the curve: 0.5921
## 95% CI: 0.5197-0.6645 (DeLong)
roc.test(roc_meld, roc_meld_3, method = "delong")
##
## DeLong's test for two correlated ROC curves
##
## data: roc_meld and roc_meld_3
## Z = -0.59681, p-value = 0.5506
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
## -0.03475716 0.01853095
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.5920854 0.6001985
roc_meld_3
roc(final_master$hrs_responders_cat_2, final_master$pred_day0_meld_3, ci=TRUE,direction = "<")
## Setting levels: control = 0, case = 1
##
## Call:
## roc.default(response = final_master$hrs_responders_cat_2, predictor = final_master$pred_day0_meld_3, direction = "<", ci = TRUE)
##
## Data: final_master$pred_day0_meld_3 in 109 controls (final_master$hrs_responders_cat_2 0) < 134 cases (final_master$hrs_responders_cat_2 1).
## Area under the curve: 0.6002
## 95% CI: 0.5281-0.6723 (DeLong)
roc_na
roc(final_master$hrs_responders_cat_2, final_master$pred_day0_meld_na, ci=TRUE,direction = "<")
## Setting levels: control = 0, case = 1
##
## Call:
## roc.default(response = final_master$hrs_responders_cat_2, predictor = final_master$pred_day0_meld_na, direction = "<", ci = TRUE)
##
## Data: final_master$pred_day0_meld_na in 109 controls (final_master$hrs_responders_cat_2 0) < 134 cases (final_master$hrs_responders_cat_2 1).
## Area under the curve: 0.5927
## 95% CI: 0.5203-0.6651 (DeLong)
roc.test(roc_na, roc_meld_3, method = "delong")
##
## DeLong's test for two correlated ROC curves
##
## data: roc_na and roc_meld_3
## Z = -0.88085, p-value = 0.3784
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
## -0.024178179 0.009184341
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.5927016 0.6001985
roc_clif
roc(final_master$hrs_responders_cat_2, final_master$pred_clif, ci=TRUE,direction = "<")
## Setting levels: control = 0, case = 1
##
## Call:
## roc.default(response = final_master$hrs_responders_cat_2, predictor = final_master$pred_clif, direction = "<", ci = TRUE)
##
## Data: final_master$pred_clif in 109 controls (final_master$hrs_responders_cat_2 0) < 134 cases (final_master$hrs_responders_cat_2 1).
## Area under the curve: 0.5694
## 95% CI: 0.4959-0.6429 (DeLong)
roc.test(roc_clif, roc_meld_3, method = "delong")
##
## DeLong's test for two correlated ROC curves
##
## data: roc_clif and roc_meld_3
## Z = -0.79561, p-value = 0.4263
## alternative hypothesis: true difference in AUC is not equal to 0
## 95 percent confidence interval:
## -0.10658804 0.04503799
## sample estimates:
## AUC of roc1 AUC of roc2
## 0.5694235 0.6001985