options(scipen=999)
library(tidyverse)
library(tableone)
library(writexl)
library(gtsummary)
library(survival)
library(ggplot2)
library(readxl)
library(xlsx)
library(mice)
library(survminer)
library(tidycmprsk)
library(pec)
library(casebase)
library(cmprsk)
library(crrstep)
## get site data for extral variables
setwd("C:\\Users\\to909\\Desktop\\Andrew final\\Data\\Final Data")
Baylor0913<- read_excel("Baylor0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
Indiana0913 <- read_excel("Indiana0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
Jacksonville0913<- read_excel("Jacksonville0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
Kentukey0913<- read_excel("Kentukey0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
MCW0913<- read_excel("MCW0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
MGH0913<- read_excel("MGH0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
Michigan0913<- read_excel("Michigan0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
Oschner0913<- read_excel("Oschner0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
Rochester0913<- read_excel("Rochester0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
USC0913<- read_excel("USC0913.xlsx") %>% select(study_id) %>% mutate(days_admitted_rrt= NA,rrt_crrt= NA ) #USC doesn't have days_admitted_rrt
Yale0913<- read_excel("Yale0913.xlsx") %>% select(study_id,days_admitted_rrt,rrt_crrt)
all_extral_var <- rbind(Baylor0913,Indiana0913,Jacksonville0913,Kentukey0913,MCW0913,MGH0913,Michigan0913,Oschner0913,Rochester0913,USC0913,Yale0913)
# set wd
setwd("C:/Users/to909/Desktop/Projects/RRT subgroup analysis/RRT-")
# read master data
master <- read_excel("All_edited_new_12142022.xlsx")
master <- merge(master,all_extral_var,by="study_id",all.x = T)
# subgroup rrt
master_rrt <- master %>% filter(rrt==1) ####### data cleaning ####
model_master_rrt <- master_rrt %>% select(time_90days,status_90days,age_admission,sex,White,hispanic_race,liver_transplant_listed,
final_type_of_aki, site, etiology_cirrhosis, albumin_given_admission,
loop_diuretic,
aldosterone_antagonist ,
lactulose,
rifaximin,
prophylactic_antibiotic,
nsaids,
beta_blockers,
diabetes,
cad,
htn,
ckd,
ascites_admission,
encephalopathy_admission,
gi_bleed_admission,
peritonitis_admission,
hcc_admission,
tips_admission,
lvp_admission,
alcoholic_hepatitis_admission,
MELD_Na_baseline,
creatinine_admission,
log_na_admit,
log_k_admit,
log_cl_admit,
log_co2_admit,
log_bun_admit,
log_alt_admit,
log_alkphos_admit,
log_tb_admit,
log_alb_admit,
log_inr_admit,
log_wbc_admit,
log_hct_admit,
log_plt_admit,
sbp_admission,
dbp_admission
)
all_var_model <- names(model_master_rrt)
cat_var_model <- c("sex","White","hispanic_race","liver_transplant_listed",
"site",
"final_type_of_aki",
"encephalopathy_admission",
"etiology_cirrhosis" ,
"albumin_given_admission",
"loop_diuretic" ,
"aldosterone_antagonist" ,
"lactulose" ,
"rifaximin" ,
"prophylactic_antibiotic",
"nsaids",
"beta_blockers",
"diabetes" ,
"cad" ,
"ckd",
"htn",
"ascites_admission",
"encephalopathy_admission",
"gi_bleed_admission" ,
"peritonitis_admission" ,
"hcc_admission" ,
"tips_admission" ,
"lvp_admission",
"alcoholic_hepatitis_admission")
num_var_model <- setdiff(all_var_model,cat_var_model)
model_master_rrt <- model_master_rrt %>% mutate_at(cat_var_model,factor)
model_master_rrt <- model_master_rrt %>% drop_na()covar_names <- c( "age_admission","sex","White","hispanic_race","liver_transplant_listed","site","MELD_Na_baseline")
covar_matrix <- model.matrix(as.formula(paste0("~ ", paste(covar_names, collapse='+'))), data = model_master_rrt)[,-1]
model1 <- cmprsk::crr(ftime=model_master_rrt$time_90days,fstatus=model_master_rrt$status_90days,cov1=covar_matrix,failcode=1,cencode=0)
summary(model1)## Competing Risks Regression
##
## Call:
## cmprsk::crr(ftime = model_master_rrt$time_90days, fstatus = model_master_rrt$status_90days,
## cov1 = covar_matrix, failcode = 1, cencode = 0)
##
## coef exp(coef) se(coef) z p-value
## age_admission 0.007362 1.007 0.00684 1.07654 0.28000000000
## sex2 0.000937 1.001 0.14160 0.00661 0.99000000000
## White1 0.159939 1.173 0.17883 0.89439 0.37000000000
## hispanic_race1 0.093637 1.098 0.24240 0.38629 0.70000000000
## liver_transplant_listed1 -1.465297 0.231 0.23503 -6.23448 0.00000000045
## siteindiana 1.017996 2.768 0.38416 2.64991 0.00810000000
## sitejacksonville 0.024069 1.024 0.49484 0.04864 0.96000000000
## sitekentukey 1.574550 4.829 0.39385 3.99787 0.00006400000
## siteMCW 0.898470 2.456 0.35009 2.56638 0.01000000000
## sitemgh 0.956257 2.602 0.35418 2.69991 0.00690000000
## sitemichigan 0.817159 2.264 0.44388 1.84093 0.06600000000
## siteoschner 1.296138 3.655 0.42188 3.07232 0.00210000000
## siterochester 0.113181 1.120 0.38122 0.29689 0.77000000000
## siteusc 0.256459 1.292 0.45707 0.56109 0.57000000000
## siteyale 0.359825 1.433 0.38372 0.93773 0.35000000000
## MELD_Na_baseline 0.034892 1.036 0.00817 4.26870 0.00002000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## age_admission 1.007 0.993 0.994 1.021
## sex2 1.001 0.999 0.758 1.321
## White1 1.173 0.852 0.827 1.666
## hispanic_race1 1.098 0.911 0.683 1.766
## liver_transplant_listed1 0.231 4.329 0.146 0.366
## siteindiana 2.768 0.361 1.303 5.876
## sitejacksonville 1.024 0.976 0.388 2.702
## sitekentukey 4.829 0.207 2.231 10.449
## siteMCW 2.456 0.407 1.237 4.878
## sitemgh 2.602 0.384 1.300 5.209
## sitemichigan 2.264 0.442 0.949 5.404
## siteoschner 3.655 0.274 1.599 8.356
## siterochester 1.120 0.893 0.530 2.364
## siteusc 1.292 0.774 0.528 3.165
## siteyale 1.433 0.698 0.676 3.040
## MELD_Na_baseline 1.036 0.966 1.019 1.052
##
## Num. cases = 363
## Pseudo Log-likelihood = -1102
## Pseudo likelihood ratio test = 132 on 16 df,model_selection_variables <- paste(names(model_master_rrt), collapse='+')
# remove intercept
model_selection_data <- as.data.frame(model.matrix(as.formula(paste0("~ ", model_selection_variables)), data = model_master_rrt)[,-1])
model_selection_covariates <- paste(setdiff(names(model_selection_data), c("time_90days","status_90days")), collapse='+')
model_section_func <- as.formula(paste0("Hist(time_90days, status_90days) ~ ",model_selection_covariates))
BIC_selection <- pec::selectFGR(model_section_func, cause = 1, data = model_selection_data,rule = "BIC", direction = c("backward","forward"))
BIC_selection$fit ##
## Right-censored response of a competing.risks model
##
## No.Observations: 363
##
## Pattern:
##
## Cause event right.censored
## 1 214 0
## 2 58 0
## unknown 0 91
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## riskRegression::FGR(formula = Hist(time_90days, status_90days) ~
## White1 + liver_transplant_listed1 + siteindiana + sitekentukey +
## siteMCW + sitemgh + siteoschner + etiology_cirrhosis2 +
## etiology_cirrhosis6 + etiology_cirrhosis7 + loop_diuretic1 +
## lactulose1 + rifaximin1 + prophylactic_antibiotic1 +
## ckd1 + peritonitis_admission1 + alcoholic_hepatitis_admission1 +
## log_co2_admit + log_alt_admit + log_alkphos_admit + log_tb_admit +
## log_alb_admit + log_plt_admit + sbp_admission, data = data,
## cause = cause)
##
## coef exp(coef) se(coef) z
## White1 0.3945 1.484 0.16669 2.37
## liver_transplant_listed1 -1.7786 0.169 0.23707 -7.50
## siteindiana 0.6972 2.008 0.24474 2.85
## sitekentukey 1.0228 2.781 0.28766 3.56
## siteMCW 0.7079 2.030 0.24010 2.95
## sitemgh 0.9004 2.461 0.19356 4.65
## siteoschner 0.8827 2.417 0.32230 2.74
## etiology_cirrhosis2 0.6469 1.910 0.21633 2.99
## etiology_cirrhosis6 0.5740 1.775 0.22070 2.60
## etiology_cirrhosis7 0.5448 1.724 0.24369 2.24
## loop_diuretic1 0.2547 1.290 0.15838 1.61
## lactulose1 0.3266 1.386 0.17554 1.86
## rifaximin1 -0.3156 0.729 0.19338 -1.63
## prophylactic_antibiotic1 -0.4111 0.663 0.25088 -1.64
## ckd1 -0.5016 0.606 0.18656 -2.69
## peritonitis_admission1 0.2713 1.312 0.18739 1.45
## alcoholic_hepatitis_admission1 -0.3505 0.704 0.20485 -1.71
## log_co2_admit -0.5287 0.589 0.23384 -2.26
## log_alt_admit 0.1755 1.192 0.08120 2.16
## log_alkphos_admit -0.2376 0.789 0.10373 -2.29
## log_tb_admit 0.2945 1.342 0.06856 4.30
## log_alb_admit -0.6795 0.507 0.33761 -2.01
## log_plt_admit -0.3289 0.720 0.12314 -2.67
## sbp_admission -0.0153 0.985 0.00471 -3.25
## p-value
## White1 0.018000000000000
## liver_transplant_listed1 0.000000000000063
## siteindiana 0.004400000000000
## sitekentukey 0.000380000000000
## siteMCW 0.003200000000000
## sitemgh 0.000003300000000
## siteoschner 0.006200000000000
## etiology_cirrhosis2 0.002800000000000
## etiology_cirrhosis6 0.009300000000000
## etiology_cirrhosis7 0.025000000000000
## loop_diuretic1 0.110000000000000
## lactulose1 0.063000000000000
## rifaximin1 0.100000000000000
## prophylactic_antibiotic1 0.100000000000000
## ckd1 0.007200000000000
## peritonitis_admission1 0.150000000000000
## alcoholic_hepatitis_admission1 0.087000000000000
## log_co2_admit 0.024000000000000
## log_alt_admit 0.031000000000000
## log_alkphos_admit 0.022000000000000
## log_tb_admit 0.000017000000000
## log_alb_admit 0.044000000000000
## log_plt_admit 0.007600000000000
## sbp_admission 0.001200000000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## White1 1.484 0.674 1.070 2.057
## liver_transplant_listed1 0.169 5.922 0.106 0.269
## siteindiana 2.008 0.498 1.243 3.244
## sitekentukey 2.781 0.360 1.583 4.887
## siteMCW 2.030 0.493 1.268 3.249
## sitemgh 2.461 0.406 1.684 3.596
## siteoschner 2.417 0.414 1.285 4.547
## etiology_cirrhosis2 1.910 0.524 1.250 2.918
## etiology_cirrhosis6 1.775 0.563 1.152 2.736
## etiology_cirrhosis7 1.724 0.580 1.069 2.780
## loop_diuretic1 1.290 0.775 0.946 1.760
## lactulose1 1.386 0.721 0.983 1.956
## rifaximin1 0.729 1.371 0.499 1.065
## prophylactic_antibiotic1 0.663 1.508 0.405 1.084
## ckd1 0.606 1.651 0.420 0.873
## peritonitis_admission1 1.312 0.762 0.908 1.894
## alcoholic_hepatitis_admission1 0.704 1.420 0.471 1.052
## log_co2_admit 0.589 1.697 0.373 0.932
## log_alt_admit 1.192 0.839 1.016 1.397
## log_alkphos_admit 0.789 1.268 0.643 0.966
## log_tb_admit 1.342 0.745 1.174 1.536
## log_alb_admit 0.507 1.973 0.262 0.982
## log_plt_admit 0.720 1.389 0.565 0.916
## sbp_admission 0.985 1.015 0.976 0.994
##
## Num. cases = 363
## Pseudo Log-likelihood = -1055
## Pseudo likelihood ratio test = 225 on 24 df,
##
## Convergence: TRUE#cox ph regression model transplant free — outcome death or transplant
model_master_rrt$status_90days_combined <- ifelse(model_master_rrt$status_90days==1|model_master_rrt$status_90days==2,1,0)
cox_model_combined <- coxph(Surv(time_90days,status_90days_combined) ~ age_admission+sex+White+hispanic_race+site+liver_transplant_listed+MELD_Na_baseline, data =model_master_rrt)
summary(cox_model_combined)## Call:
## coxph(formula = Surv(time_90days, status_90days_combined) ~ age_admission +
## sex + White + hispanic_race + site + liver_transplant_listed +
## MELD_Na_baseline, data = model_master_rrt)
##
## n= 363, number of events= 272
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_admission 0.003987 1.003995 0.005153 0.774 0.439109
## sex2 -0.016757 0.983383 0.129261 -0.130 0.896853
## White1 0.041176 1.042035 0.155846 0.264 0.791620
## hispanic_race1 -0.146831 0.863440 0.221962 -0.662 0.508281
## siteindiana 1.024198 2.784862 0.399980 2.561 0.010448 *
## sitejacksonville 1.096194 2.992755 0.414786 2.643 0.008222 **
## sitekentukey 1.468399 4.342277 0.408509 3.595 0.000325 ***
## siteMCW 0.752270 2.121811 0.401326 1.874 0.060867 .
## sitemgh 1.035486 2.816475 0.388739 2.664 0.007729 **
## sitemichigan 0.522844 1.686818 0.508137 1.029 0.303507
## siteoschner 1.351810 3.864412 0.432328 3.127 0.001767 **
## siterochester 0.160980 1.174662 0.408115 0.394 0.693250
## siteusc 1.129754 3.094894 0.411814 2.743 0.006081 **
## siteyale 0.423606 1.527460 0.427582 0.991 0.321831
## liver_transplant_listed1 0.125548 1.133769 0.141228 0.889 0.374019
## MELD_Na_baseline 0.033407 1.033971 0.008042 4.154 0.0000326 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_admission 1.0040 0.9960 0.9939 1.014
## sex2 0.9834 1.0169 0.7633 1.267
## White1 1.0420 0.9597 0.7678 1.414
## hispanic_race1 0.8634 1.1582 0.5589 1.334
## siteindiana 2.7849 0.3591 1.2716 6.099
## sitejacksonville 2.9928 0.3341 1.3274 6.747
## sitekentukey 4.3423 0.2303 1.9498 9.670
## siteMCW 2.1218 0.4713 0.9663 4.659
## sitemgh 2.8165 0.3551 1.3147 6.034
## sitemichigan 1.6868 0.5928 0.6231 4.567
## siteoschner 3.8644 0.2588 1.6561 9.017
## siterochester 1.1747 0.8513 0.5279 2.614
## siteusc 3.0949 0.3231 1.3807 6.937
## siteyale 1.5275 0.6547 0.6607 3.531
## liver_transplant_listed1 1.1338 0.8820 0.8596 1.495
## MELD_Na_baseline 1.0340 0.9671 1.0178 1.050
##
## Concordance= 0.659 (se = 0.018 )
## Likelihood ratio test= 74.04 on 16 df, p=0.000000002
## Wald test = 66.67 on 16 df, p=0.00000004
## Score (logrank) test = 71.31 on 16 df, p=0.000000006