This webpage is designed to explain the overall steps that readers can follow in order to reproduce the results. Here is the overall process:

Reading the file and defining the variable types based on their defination
Preparing the data
- demonstrating the missing parts
- removing the variables that are not relevant
- regrouping the levels of categorical variable
- imputation
Feature selection
Prediction & Analysis

Here we use Thoracic (Heart) transplantation data which is provided United Organ Sharing Network (UNOS). Per data access agreement, we cannot share data. For obtainig the data please contact the data registry. The data has numerous variables. The variable definition is available through the data dictionary or through the Scientific Registry of Transplant Recepients (SRTR).

The R codes is divided in section that are represented in the navigation bar at the left. In order to access the R codes implemented for this study, the reader can click show. In addition to the functions that are loaded from packages, a list of customized functions are created. Documentation of the functions is available in functions.R

Figure 1 provides an overview of our two-stage framework for obtaining monotonic survival probabilities over time.

Figure 1: The proposed two-stage framework for obtaining monotonically decreasing heart transplantation survival probabilities.

1 Data Exploration

In the snippet below, we load list of R packages and the custom functions that are employed in this study. For the convinience, The pacman package checks the packages (if they are not existing, it installs them) then load them.

rm(list = ls()) # clear global environment
graphics.off() # close all graphics
if(require(pacman)==FALSE) install.packages("pacman") # needs to be installed first
# p_load is equivalent to combining both install.packages() and library()
pacman::p_load(htmltools,Biocomb, ranger,car,AUC,stringr , caret, conflicted, DataExplorer, # important for exploratory data analysis 
               dataPreparation, data.table, DT, effects, emmeans, ggplot2, haven, lmtest, lsmeans, magrittr, mltools, party, readxl, tidyverse, # important analytic packages
               sjmisc, sjPlot, snow, varImp, yhat,checkpoint
)

conflict_prefer("arrange", "dplyr")
conflict_prefer("recode", "dplyr")
conflict_prefer("mutate", "dplyr")
conflict_prefer("makeCluster", "snow")

checkpoint("2020-07-16")
source("https://raw.githubusercontent.com/transplantation/unos-ht/master/data/functions.R") # Has our custom functions

1.1 Loading the UNOS Dataset

In this section we change the working directory to a address of a folder in our PC. The reader needs to edit the address of the folder after attaining the dataset from UNOS. Here, we select the heart transplantation records.

df <- read_sas("thoracic_data.sas7bdat")
df$ID <- row.names(df) # Index to be used in creating/tracking the training and test samples
orgTable <- table(df[, 'WL_ORG'])
df %<>% dplyr::filter(WL_ORG == "HR") # filtering to keep only heart transplants

Our original UNOS data consisted of 159318 observations, which contained: (a) 428 unspecified transplant types, (b) 3148 combined heart-lung transplants, (c) 103570 heart transplants, and (d) 52172 lung transplants. Given our focus on heart transplants, we have only retained the observations where only a heart transplant is performed. The resulting data is resaved into our data.frame titled df, which now consists of 103570 observations and 495 variables (which have been unchanged from our filtering operation).

1.2 Exploring the UNOS Dataset

In the following code chunk, we convert the blank cells to “NA” before we study the distribution of missing values in the data.

# The variable YR_ENTRY_US_TCR was recorded with time zone, we removed the time zone.
df$YR_ENTRY_US_TCR <- as.Date(df$YR_ENTRY_US_TCR)
# We exclude the Date variables when coverting the blank cells to NA as these variables were well coded. 
index_date <- which(unlist(lapply(df, class))=="Date")
df[,-index_date] %<>% mutate_all(na_if, "") %>% mutate_all(na_if, " ") %>% mutate_all(na_if, "U")

Fraction Data Missing

naVec <- is.na(df) %>% colSums() %>% 
  `/`(nrow(df)) # To make it percent missing per variable
names(naVec) <- colnames(df)
naVec.quantiles <- quantile(naVec) %>% round(2)

ggplot() + geom_histogram(aes(naVec), bins = 10) + 
  xlab("Fraction Data Missing Per Variable") + ylab("Number of Variables")+theme_minimal()

Based on the heart transplantation data, 457 variables out of the 103570 have at least one data point missing data. Furthermore, we have computed the percent missing data per variable; the values of the first, second, and third quantiles (in terms of percent of data missing per column) correspond to 37%, 63%, and 97%,respectively. Thus, it is imperative to examine how these variables can be preprocessed prior to any analysis.

Missing Data (Sample Columns)

In the plot below, we sample 30 columns at random from the UNOS dataset to show the actual percentage of the data that is missing for each variable. The colors are used to denote the data quality for that column using a traffic light scheme (where green is good and red is bad).

heart.na.plot <- df[,sample(colnames(df), 30)] %>% plot_missing()

Missing Data (Printout for All Variables)

In the chunk below, we report the percent missing per each variable. The printout is arranged in a descending order of missing values.

data.frame(variableName = names(naVec), percentMissing = round(naVec*100, 2), row.names = NULL) %>%
  arrange(desc(percentMissing)) %>% datatable()

Column Types

R initially divides the columns of different types. We summarize these in the table below.

types <- dataTypes(df) # see functions.R file
numeric.index <- which(types$classes.df=="numeric")
character.index <- which(types$classes.df=="character")
datatable(types) # printing the different types

The reader should note that R reads any variable that purely consists of numeric values to be numeric. However, from a data analysis perspective, a lot of these variables will need to be converted into factors since they represent factor levels rather than numeric values (e.g., donor/recipient ethinicity categories, allocation type, HLA mismatch level, education, etc.). Thus, the numeric variables should be examined more closely to ensure that they are properly coded by the software.

2 Data Preperation

In this section, we will detailedly document the data preperation precedures we have done in the study.

2.1 Technically Correct Variables

In this subsection, we will capitalize on the file titled = “STAR File Documentation.xls” provided by UNOS to automatically check and convert any variable that was assigned an incorrect type by the R software. The reader should note that the outcome of this data transformation step can be referred as technically correct data (De Jonge and Van Der Loo 2013). Note that we convert all factor columns to character to facilitate their manipulation in R.

varInfo <- read_excel("STAR File Documentation.xls", sheet = "THORACIC DATA", skip =1, col_types = c(rep("text",3), rep("date",2), rep("guess",7)) )

# fix an error in the variable information document, there were two start dates: 01-Oct-87, 01-Oct-90 for the variable: CMV_DON. We assign 01-Oct-90 for the start date.

varInfo$`VAR START DATE`[51] <- "1990-10-01"

# pulling names of variables that had type char OR having numeric categories (explained in SAF)
factorVars <- varInfo %>% 
  dplyr::filter(`DATA TYPE`=="CHAR" | (`DATA TYPE`=="NUM" & !is.na(`SAS ANALYSIS FORMAT`)) ) %>%
  pull(`VARIABLE NAME`)

# The following variables are qualitative based on the UNOS documentation but they were not identified. 
Other_categorical_vars <- c("BLOOD_INF_DON", "HBV_CORE_DON", "HBV_SUR_ANTIGEN",   "HCV_SEROSTATUS", "HYPERTENS_DUR_DON", "IABP_TRR", "LAST_INACT_REASON", "OTHER_INF_DON", "STERNOTOMY_TRR", "URINE_INF_DON")

factorVars <- c(factorVars, Other_categorical_vars)

df[, factorVars] %<>%  lapply(as.character) # converting factors to character to facilitate preprocessing

Based on the analysis above, we have identified 322 categorical variables. Note the R has previously only identified character/categorical variables. The names of the 322 variables can be accessed in the factorVars.RDS file. Furthermore, we applied the as.factor() to ensure that R recognizes these variables as categorical.

2.2 Rowwise Operations: Additional Filters

In the code chunk below, we first filter to keep adult patients since our study focuses on only adult transplantations. Then we remove individuals whose weights were lower than .01% of adults’ weights in the data (less than 30lbs) or whose heights were less than .01% adults’ heights in the data (lower than 122 cms).

df <- df %>% subset(AGE>=18) %>% # filtering to keep only Adults
# we excluded too light or too short people in the following code
subset(WGT_KG_DON_CALC >= quantile(WGT_KG_DON_CALC, 0.0001, na.rm = TRUE)) %>% 
subset(WGT_KG_TCR >= quantile(WGT_KG_TCR, 0.0001, na.rm = TRUE)) %>%
subset(HGT_CM_DON_CALC >= quantile(HGT_CM_DON_CALC, 0.0001, na.rm = TRUE)) %>% 
subset(HGT_CM_TCR >= quantile(HGT_CM_TCR, 0.0001, na.rm = TRUE))

2.3 Columnwise Operations: Dropping Variables

In the code chunk below, we fist filter/remove any observations with missing GSTATUS or GTIME. Our rationale for doing this prior to removing uniformative/irrelevant/unreliable variables was based on the fact that these variables would be used in computing our transplantation OUTCOME and hence, we did not want to lose them. After performing remove any variables that meet any of the following criteria:

Variables with a starting date on or later than Jan 01, 2000, except the variables that are corresponding to prior cardiac surgery (non-transplant) or lung surgery as these variables have been identified as important information in the literature (Vega, Schroder, and Nicoara (2017), Chen et al. (2018), Guven et al. (2018), Medved et al. (2018), Zhang et al. (2019)). These variables are:
PRIOR_CARD_SURG_TCR, PRIOR_CARD_SURG_TRR, PRIOR_CARD_SURG_TYPE_OSTXT_TCR, PRIOR_CARD_SURG_TYPE_OSTXT_TRR, PRIOR_CARD_SURG_TYPE_TCR, PRIOR_CARD_SURG_TYPE_TRR, and PRIOR_LUNG_SURG_TYPE_OSTXT_TRR.
Any variables having an end date for data collection since they would not be collected in any future transplant event.
Any post-transplant variables, with the exception of GTIME and GSTATUS (which we use later to record the transplantation outcome at 1-year post transplant)
Removing the following variables, capturing different types of antigen alleles:
DA1, DA2, RA1, RA2, DB1, DB2, RB1, RB2, RDR1, RDR2, DDR1, and DDR2.
Note that each of those variables had three possible values: 0 indicating no mismatch, 1 indicating one allele matching, and 2 indicating that both alleles mismatched.
Our rationale for removing those variables was based on the observation that the variables HLAMIS, AMIS, BMIS and DRMIS recorded summaries of matches in these antigen alleles. The interested reader is referred to Parham (2014) for a discussion of HLAMIS, and Weisdorf et al. (2008) for a detailed introduction on those four variables.
Removing the two variables DEATH_CIRCUM_DON and DEATH_MECH_DON since a cross-tabulation of their values indicated that they are not consistent, which reduces the interpretation/utility from using those variables.
Removing the variable TX_YEAR since we wanted to develop a predictive and not an explanatory model. For more details on the difference, the reader is referred to Shmueli and others (2010).
Any variables where the number of missing observations exceeded 90%
Any numeric variables where the number of missing observations exceeded 30%

varsVec <- read_rds(gzcon(url("https://raw.githubusercontent.com/My-favor-research/ISR-Supplementary-Materials.github.io/master/dropVars.RDS"))) # a vector of 69 vars to be dropped
alleles = c("DA1","DA2","RA1","RA2","DB1","DB2","RB1","RB2","RDR1","RDR2","DDR1","DDR2")

VarstoBeDropped = varInfo %>% 
dplyr::filter(`VAR START DATE` >= '2000-01-01' | 
!is.na(`VAR END DATE`) | 
str_detect(FORM, "TRF/TRR|TRR/TRF-CALCULATED|TRR/TRF|TRF") |
str_detect(`FORM SECTION`, "POST TRANSPLANT CLINICAL INFORMATION") |
str_detect(DESCRIPTION, "IDENTIFIER| DATE") |
`VARIABLE NAME` %in% varsVec |
`VARIABLE NAME` %in% alleles |
`VARIABLE NAME` %in% c("DEATH_CIRCUM_DON","DEATH_MECH_DON", "TX_YEAR") ) %>% 
pull(`VARIABLE NAME`)

# include variables corresponding to prior cardiac surgery (non-transplant) or lung surgery
VarstoBeDropped <- setdiff(VarstoBeDropped, c("PRIOR_CARD_SURG_TCR", "PRIOR_CARD_SURG_TRR", "PRIOR_CARD_SURG_TYPE_OSTXT_TCR", "PRIOR_CARD_SURG_TYPE_OSTXT_TRR",
"PRIOR_CARD_SURG_TYPE_TCR", "PRIOR_CARD_SURG_TYPE_TRR", "PRIOR_LUNG_SURG_TYPE_OSTXT_TRR"))


df %<>% dplyr::filter(!is.na(GSTATUS) & !is.na(GTIME) ) %>% # remove obs missing graft status/time
select(!all_of(VarstoBeDropped)) %>% # exclude variables to be dropped
discard(~sum(is.na(.x))/length(.x)* 100 >= 90) %>% # remove vars with >= 90 % missing
dropNumericMissing(percent = 30) # drop numeric columns w/ >= 30% missing (see custom_functions.R)

2.4 Creating Variables based on the Literature and Our Exploration of the Data

In the code chunk below, we create the following variables:

ISCHTIME was recoded to be in minutes to make the interpretation consistent with Medved et al. (2018)
PVR, which calculates the recipent’s pulmonary vascular resistance
An ECMO index variable was created to indicate whether the ECMO machine was utilized at either registration or at transplant
CARD_SURG, which captures whether the recipent had prior cardiac surgeries (non-transplant)
BMI_CHNG, which captures the percent change in the recipient’s BMI from transplantation time to registration/listing
WGT_CHNG, which captures the percent change in the recipient’s weight from transplantation time to registration/listing
HGT_CHNG, which captures the percent change in the recipient’s height from transplantation time to registration/listing
AGE_DIFF, which captures the absolute value of the difference between the recipient’s and deceased donor’s age
BMI_DIFF, which captures the absolute value of the difference between the recipient’s and deceased donor’s BMI
ETH_MAT, which captures whether the recipent and deceased donor were of the same ethinicity
GEN_MAT, which captures whether the recipent and deceased donor were of the same gender
ANCEF, which captures whether the deceased donor was prescribed Ancef within 24 hours or procurement
DOPAMINE, which captures whether the deceased donor was prescribed Dopamine within 24 hours or procurement
HEPARIN, which captures whether the deceased donor was prescribed Herapin within 24 hours or procurement
ZOSYN, which captures whether the deceased donor was prescribed Zosyn within 24 hours or procurement
OUTCOME, which captures whether the recipient has survived at one-year post transplant based on combining information from GTIME and GSTATUS. The censoring procedure used herein is identical to that applied in Dag et al. (2016); Dag et al. (2017).

After these variables were created, we dropped both GTIME and GSTATUS as well as the four variables capturing the medications prescribed/applied on the deceased donor (i.e. PT_OTH1_OSTXT_DON – PT_OTH4_OSTXT_DON). Furthemore, we converted all indicator columns into characters.

df %<>% mutate(ISCHTIME = ISCHTIME*60) %>%
mutate(PVR = (HEMO_PA_MN_TRR - HEMO_PCW_TRR)*79.72/HEMO_CO_TRR) %>%
mutate(ECMO = ifelse(ECMO_TCR + ECMO_TRR == 0, 0, 1)) %>% 
mutate(CARD_SURG = if_else(PRIOR_CARD_SURG_TCR=="Y"|PRIOR_CARD_SURG_TRR=="Y", true="Y", false=ifelse(PRIOR_CARD_SURG_TCR=="N"&PRIOR_CARD_SURG_TRR=="N", "N", NA))) %>%
mutate(BMI_CHNG = 100*(BMI_CALC- INIT_BMI_CALC)/INIT_BMI_CALC) %>% 
mutate(WGT_CHNG = 100*(WGT_KG_CALC - INIT_WGT_KG_CALC)/INIT_WGT_KG_CALC) %>% 
mutate(HGT_CHNG = 100*(HGT_CM_CALC - INIT_HGT_CM_CALC)/INIT_HGT_CM_CALC) %>% 
mutate(AGE_DIFF = abs(AGE - AGE_DON)) %>% 
mutate(BMI_DIFF = abs(BMI_CALC - BMI_DON_CALC)) %>% 
mutate(ETH_MAT = if_else(ETHCAT == ETHCAT_DON, true = "Y", false = "N")) %>% 
mutate(GEN_MAT = if_else(GENDER == GENDER_DON, true = "Y", false = "N")) %>%
mutate(ANCEF = if_else(str_detect(PT_OTH1_OSTXT_DON, 'ANCEF') |
str_detect(PT_OTH2_OSTXT_DON, 'ANCEF') |
str_detect(PT_OTH3_OSTXT_DON, 'ANCEF'),
true = "Y", false = "N") ) %>% 
mutate(DOPAMINE = if_else(str_detect(PT_OTH1_OSTXT_DON, 'DOPAMINE') | 
str_detect(PT_OTH2_OSTXT_DON, 'DOPAMINE') |
str_detect(PT_OTH3_OSTXT_DON, 'DOPAMINE'),
true = "Y", false = "N") ) %>% 
mutate(HEPARIN = if_else(str_detect(PT_OTH1_OSTXT_DON, 'HEPARIN') | 
str_detect(PT_OTH2_OSTXT_DON, 'HEPARIN') |
str_detect(PT_OTH3_OSTXT_DON, 'HEPARIN'),
true = "Y", false = "N") ) %>% 
mutate(ZOSYN = if_else(str_detect(PT_OTH1_OSTXT_DON, 'ZOSYN') | 
str_detect(PT_OTH2_OSTXT_DON, 'ZOSYN') |
str_detect(PT_OTH3_OSTXT_DON, 'ZOSYN'),
true = "Y", false = "N") ) %>% 
rowwise() %>% 
mutate(OUTCOME = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 1)) %>% # see custom_functions.R
select(-c(starts_with("PT_OTH")))

df %<>% ungroup() %>% 
mutate_at(c('ECMO', 'CARD_SURG', 'ETH_MAT', 'GEN_MAT', 'ANCEF', 'DOPAMINE', 'HEPARIN', 'ZOSYN', 'OUTCOME'), as.character)

2.5 Reducing Number of Levels for Categorical Columns

2.5.1 Number of Distinct Values within Each Categorical Variable

Despite the application of the nearZeroVar() from the caret package, several of the remaining character/factor variables still exhibit a large number of unique values. The purpose of the code chunk below is to quantify the variability in the remaining character/factor variables.

df %>% select(-ID) %>%  summarise_if(is.character, n_distinct, na.rm = TRUE) -> numLevels

If we were to exclude the ID column, the number of distinct values in each character variable (equivallently the number of non-NA levels if these variables were encoded as factors) has the following summary statistics:

Min. : 1.000
1st Qu.: 2.000
Median : 2.000
Mean : 4.975
3rd Qu.: 4.750
Max. :49.000

Accordingly, we will have to examine how to reduce the number of levels within each variable to be able to use these variables in the machine learning stage; we do not want to generate a very large number of dummy variables and we will group some of those levels by capitalizing on the similarities in interpretting some of these levels.

2.5.2 Grouping Operations

While recipes package provides an excellent pipeline for data pre-processing, we have elected to manually go through the process of grouping variables since this would allow us to (a) utilize the information in the factor levels to combine similar categories (whether due to location, education status or medical diagnosis codes); and (b) the integration of the recipes package into the training of machine learning models would require to train the models repeatdly for each recipe per the documentation of Max (2020). In the code chunk below, we provide the details for how we regrouped several values/levels of each categorical/character vector with the hope of improving both the efficiency and prediction performance of the different machine learning models. In addition, we have also used the recode() to rename some of the categories since R can produce errors when one-hot encoding is applied if the name of the category is numeric; therefore, some of the operations performed below merely changed the name of the category to avoid future errors when processing the data. After regrouping, we performed two additional steps of variable elimination: (a) removing variables that have more than 90% of the non-missing observations in one level; (b) using the caret nearZeroVar() function to remove variables that the percentage of unique values in the variable is less than 0.01%.

df$ABO %<>% recode(A = 'A', A1 = 'A', A2 = 'A', B = 'B', O = 'O', AB = 'AB', A1B = 'AB', A2B = 'AB', .default = "OTHER", .missing = NA_character_) 
df$ABO_DON %<>% recode(A = 'A', A1 = 'A', A2 = 'A', B = 'B', O = 'O', AB = 'AB', A1B = 'AB', A2B = 'AB', .default = "OTHER", .missing = NA_character_)

df$ABO_MAT %<>% recode(`1` = "IDENTICAL", `2` = "NOT_IDENTICAL", `3` = "NOT_IDENTICAL", .default = "OTHER", .missing = NA_character_)

df$AMIS %<>% recode(`0` = "NO_MISMATCH", `1` = "ONE_MISMATCHED", `2` = "TWO_MISMATCHED", .default = "OTHER", .missing = NA_character_)

df$BMIS %<>% recode(`0` = "NO_MISMATCH", `1` = "ONE_MISMATCHED", `2` = "TWO_MISMATCHED", .default = "OTHER", .missing = NA_character_)

df$BRONCHO_LT_DON %<>% recode(`1` = "NO", `2` = "NORMAL", `3` = "ABNORMAL", `4` = "ABNORMAL", `5` = "ABNORMAL",`6` = "ABNORMAL", `7` ="OTHER", `998` = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$BRONCHO_RT_DON %<>% recode(`1` = "NO", `2` = "NORMAL", `3` = "ABNORMAL", `4` = "ABNORMAL", `5` = "ABNORMAL",`6` = "ABNORMAL", `7` ="OTHER", `998` = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$CHEST_XRAY_DON %<>% recode(`0` = "UNKNOWN", `1` = "NO", `2` = "NORMAL", `3` = "ABNORMAL_SINGLE", `4` = "ABNORMAL_SINGLE", `5` = "ABNORMAL_BOTH",
`998` = "UNKNOWN", `999` = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$COD_CAD_DON %<>% recode(`1` = 'ANOXIA', `2` = 'CEREBROVASCULAR_STROKE', `3` = 'HEAD_TRAUMA',`4` = 'OTHER', `999` = 'OTHER', Unknown = 'UNKNOWN', .default = "OTHER", .missing = NA_character_)

df$CORONARY_ANGIO %<>% recode(`1` = "NO", `2` = "YES_NORMAL", `3` = "YES_ABNORMAL", .default = "OTHER", .missing = NA_character_) 

df$CMV_DON %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$DIAB %<>% recode(`1` = 'NO', `2` = 'ONE', `3` = 'TWO', `4` = 'OTHER', `5` = 'OTHER', `998` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_)

df$DIAG %<>% recode(`1000` = 'DILATED_MYOPATHY_IDI', `1001` = 'DILATED_MYOPATHY_OTH', `1002` = 'DILATED_MYOPATHY_OTH', `1003` = 'DILATED_MYOPATHY_OTH',`1004` = 'DILATED_MYOPATHY_OTH', `1005` = 'DILATED_MYOPATHY_OTH',`1006` = 'DILATED_MYOPATHY_OTH', `1007` = 'DILATED_MYOPATHY_ISC', `1049` = 'DILATED_MYOPATHY_OTH',
`1200` = 'CORONARY', `999` = 'OTHER', `1999` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_)

df$TCR_DGN %<>% recode(`1000` = 'DILATED_MYOPATHY_IDI', `1001` = 'DILATED_MYOPATHY_OTH', `1002` = 'DILATED_MYOPATHY_OTH', `1003` = 'DILATED_MYOPATHY_OTH',`1004` = 'DILATED_MYOPATHY_OTH', `1005` = 'DILATED_MYOPATHY_OTH',`1006` = 'DILATED_MYOPATHY_OTH', `1007` = 'DILATED_MYOPATHY_ISC', `1049` = 'DILATED_MYOPATHY_OTH',
`1200` = 'CORONARY', `999` = 'OTHER', `1999` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_)

df$THORACIC_DGN %<>% recode(`1000` = 'DILATED_MYOPATHY_IDI', `1001` = 'DILATED_MYOPATHY_OTH', `1002` = 'DILATED_MYOPATHY_OTH', `1003` = 'DILATED_MYOPATHY_OTH',`1004` = 'DILATED_MYOPATHY_OTH', `1005` = 'DILATED_MYOPATHY_OTH',`1006` = 'DILATED_MYOPATHY_OTH', `1007` = 'DILATED_MYOPATHY_ISC', `1049` = 'DILATED_MYOPATHY_OTH',
`1200` = 'CORONARY', `999` = 'OTHER', `1999` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_)


df$DRMIS %<>% recode(`0` = "NO_MISMATCH", `1` = "ONE_MISMATCHED", `2` = "TWO_MISMATCHED", .default = "OTHER", .missing = NA_character_)

df$EBV_SEROSTATUS %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$EDUCATION %<>% recode(`1` = 'OTHER', `2` = 'GRADE', `3` = 'HIGH', `4` = 'COLLEGE', `5` = 'UNIVERSITY', `6` = 'UNIVERSITY', `996` = 'OTHER', `998` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_)

df$ETHCAT %<>% recode(`1` = 'WHITE', `2` = 'BLACK', `4` = 'HISPANIC', `5` = 'OTHER', `6` = 'OTHER',`7` = 'OTHER', `9` = 'OTHER', `998` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_)

df$ETHCAT_DON %<>% recode(`1` = 'WHITE', `2` = 'BLACK', `4` = 'HISPANIC', `5` = 'OTHER', `6` = 'OTHER',`7` = 'OTHER', `9` = 'OTHER', `998` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_)

df$FUNC_STAT_TCR %<>% recode(`1` = "ABLE", `2` = "ASSISTED", `3` = "DISABLE", `996` = "OTHER", `998` = "UNKNOWN", `2010` = "DISABLE", `2020` = "DISABLE", `2030` = "DISABLE", `2040` = "DISABLE", `2050` = "ASSISTED", `2060` = "ASSISTED" , `2070` = "ASSISTED", `2080` = "ABLE", `2090` = "ABLE", `2100` = "ABLE", .default = "OTHER", .missing = NA_character_)

df$FUNC_STAT_TRR %<>% recode(`1` = "ABLE", `2` = "ASSISTED", `3` = "DISABLE", `996` = "OTHER", `998` = "UNKNOWN", `2010` = "DISABLE", `2020` = "DISABLE", `2030` = "DISABLE", `2040` = "DISABLE",`2050` = "ASSISTED", `2060` = "ASSISTED" , `2070` = "ASSISTED", `2080` = "ABLE", `2090` = "ABLE", `2100` = "ABLE", .default = "OTHER", .missing = NA_character_)

df$HBV_CORE %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$HBV_CORE_DON %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$HBV_SUR_ANTIGEN %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$HCV_SEROSTATUS %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$HEP_C_ANTI_DON %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$HIST_DIABETES_DON %<>% recode(`1` = 'NO', `2` = 'YES', `3` = 'YES', `4` = 'YES' , `5` = 'YES', `998` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_) 

df$HYPERTENS_DUR_DON %<>% recode(`1` = 'NO', `2` = 'YES', `3` = 'YES', `4` = 'YES' , `5` = 'YES', `998` = 'UNKNOWN', .default = "OTHER", .missing = NA_character_) 

df$HIV_SEROSTATUS %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$HLAMIS %<>% recode(`0` = 'LOWEST', `1` = 'LOWEST', `2` = 'LOWEST', `3` = 'LOW', 
`4` = 'MEDIUM', `5` = 'HIGH', `6` = 'HIGHEST', .default = "OTHER", .missing = NA_character_)

df$HTLV2_OLD_DON %<>% recode(C = "UNKNOWN", I = "UNKNOWN", N = "NEGATIVE", ND = "UNKNOWN", P = "POSITIVE", PD = "UNKNOWN", U = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$INIT_STAT %<>%  recode(`1010` = 'ONE', `1020` = 'ONE', `1030` = 'TWO', `1090` ='ONE', `1999` = 'TEMPORARILY_INACTIVE', `2010` = 'ONE', `2020` = 'ONE', `2030` = 'TWO',`2090` = 'ONE', `2999` = 'TEMPORARILY_INACTIVE', `7010` = 'ACTIVE',
`7999` = 'TEMPORARILY_INACTIVE', .default = "OTHER", .missing = NA_character_)

df$END_STAT %<>%  recode(`1010` = 'ONE', `1020` = 'ONE', `1030` = 'TWO', `1090` ='ONE', `1999` = 'TEMPORARILY_INACTIVE', `2010` = 'ONE', `2020` = 'ONE', `2030` = 'TWO',`2090` = 'ONE', `2999` = 'TEMPORARILY_INACTIVE', `7010` = 'ACTIVE',
`7999` = 'TEMPORARILY_INACTIVE', .default = "OTHER", .missing = NA_character_)

df$LAST_INACT_REASON %<>% recode(`7`='HEALTH', `9`='HEALTH', `11`='HEALTH', .missing = NA_character_)

df$STERNOTOMY_TRR %<>% recode(`1`="ONE", `2`="MORE", `3`="MORE", `998`="UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$INOTROPES_TCR %<>% recode(`0` = "NO", `1` = "YES", .default = "OTHER", .missing = NA_character_) 

df$INOTROPES_TRR %<>% recode(`0` = "NO", `1` = "YES",  .default = "OTHER", .missing = NA_character_)

df$MED_COND_TRR %<>% recode(`1` = "ICU_HOSPITALIZED", `2` = "HOSPITALIZED", `3` = "NOT_HOSPITALIZED",  .default = "OTHER", .missing = NA_character_)

df$NUM_PREV_TX %<>% recode(`0` = 'ZERO', .default = 'NON_ZERO', .missing = NA_character_)

df$PRI_PAYMENT_TCR %<>% recode(`1` = "PRIVATE", `2` ="MEDICAID", `3` = "MEDICARE_FREE", `4` = "PUBLIC_OTHER",`5` = "PUBLIC_OTHER", `6` = "PUBLIC_OTHER", `7` = "PUBLIC_OTHER", `8` = "OTHER",`9` = "OTHER", `10` = "OTHER", `11` = "OTHER", `12` = "OTHER", `13` = "PUBLIC_OTHER",`14` = "OTHER", `15` = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$PRI_PAYMENT_TRR %<>% recode(`1` = "PRIVATE", `2` ="MEDICAID", `3` = "MEDICARE_FREE", `4` = "PUBLIC_OTHER",`5` = "PUBLIC_OTHER", `6` = "PUBLIC_OTHER", `7` = "PUBLIC_OTHER", `8` = "OTHER",`9` = "OTHER", `10` = "OTHER", `11` = "OTHER", `12` = "OTHER", `13` = "PUBLIC_OTHER", `14` = "OTHER", `15` = "UNKNOWN", .default = "OTHER", .missing = NA_character_)

df$PROC_TY_HR %<>% recode(`1` = "BICAVAL", `2` = "TRADITIONAL", .default = "OTHER", .missing = NA_character_)

df$PULM_INF_DON %<>% recode(`0` = "NO", `1` = "YES", .default = "OTHER", .missing = NA_character_) 

df$REGION %<>% recode(`1` = "NORTH_EAST", `2` = "NORTH_EAST", `3` = "SOUTH_EAST", `4` = "SOUTH_EAST", `5` = "WEST", `6` = "WEST", `7` = "MIDWEST", `8` = "MIDWEST", `9` = "NORTH_EAST", `10` = "MIDWEST",`11` = "SOUTH_EAST", .default = "OTHER", .missing = NA_character_)

df$SHARE_TY %<>% recode(`3` = "LOCAL", `4` = "REGIONAL", `5` = "NATIONAL", .default = "OTHER", .missing = NA_character_) 

df <- df %>% dropLowVar(percent=90) 
nzvColLoc <- nearZeroVar(df, uniqueCut = 0.01)
df <- df[,-nzvColLoc]

# df_iso is used for isotonic regression
df_iso<-df

df_iso %<>%  select(-c(OUTCOME))
df %<>%  select(-c(GSTATUS, GTIME))

Now, we have 131 variables and 48347 cases in the data.

2.6 Dealing with Missing Values & Encoding Strategy

In this section, we provide some options of dealing with missing values in either categorical or numeric variables and two approaches of encoding categorical variables.

Numerical Imputation: No imputation (Drop), Missing values are replaced with the median

Categorical Imputation: No imputation (Drop), Missing values are replaced with mode category, Missing values are labeled as “missing”, Missing values are labeled with the existing “unknown” category.

We suggest two approaches of ecnoding categorical variables: Label Encoding where each label is mapped to an integer, and One-Hot Encoding where each label is mapped to a binary variable.

In the following code chunk, we create a list (df_all) of 16 sublists corresponding to the imputation strategies and encoding strategies mentioned above. For example, df_all[[“Median_Drop_Labe”l]] stores the data after imputing medain for missing values in the numerical variables, dropping missing values in the categorical variables, and mapping each level in each categorical variable to an integer.

# first, we need to remove cases with missing values in OUTCOME as it is the response variable in our study.
df <- df[!is.na(df$OUTCOME),]

# find the indecies for categorical and numerical variables in the data
index_char <- setdiff(which(unlist(lapply(df, class))=="character"), which(colnames(df)=="ID"|colnames(df)=="OUTCOME"))

index_num <- which(unlist(lapply(df,class))=="numeric")


# define three factors
imputation_num <- c("Drop", "Median")
imputation_cat <- c("Drop", "Missing", "Mode", "Unknown")
encoding <- c("Label", "OneHot")

# Create a list of 16 elements that will store the 16 data sets for all scenarios
df_all <- rep(list(NA), 16)
names(df_all) <- apply(as.data.frame(crossing(imputation_num, imputation_cat, encoding)),1,function(x) paste(x, collapse="_"))

# prepare for imputation cases
char.impute_mode <- df[,index_char] %>% mutate_all(~ifelse(is.na(.x), names(which.max(table(.x))), .x)) 
char.impute_missing <- df[,index_char] %>% mutate_all(~ifelse(is.na(.x), "Missing", .x)) 
char.impute_unknown <- df[,index_char] %>% mutate_all(~ifelse(is.na(.x), "Unknown", .x)) 
num.impute_median <- df[,index_num] %>% mutate_all(~ifelse(is.na(.x), median(.x, na.rm = TRUE), .x)) 

# prepare for encoding with categorical imputations
mode_Label <- encode_category(char.impute_mode, "Label") 
mode_Onehot <- encode_category(char.impute_mode, "Onehot") 
missing_Label <- encode_category(char.impute_missing, "Label") 
missing_Onehot <- encode_category(char.impute_missing, "Onehot") 
unknown_Label <- encode_category(char.impute_unknown, "Label") 
unknown_Onehot <- encode_category(char.impute_unknown, "Onehot")

# prepare for no imputation case
drop_drop <- search_complete(df)  
index_Num_DD <- which(unlist(lapply(drop_drop, class))=="numeric")
index_Char_DD <- setdiff(which(unlist(lapply(drop_drop, class))=="character"), which(colnames(drop_drop)=="ID"|colnames(drop_drop)=="OUTCOME"))
drop_drop_Label <- encode_category(drop_drop[,index_Char_DD], "Label") 
drop_drop_Onehot <- encode_category(drop_drop[,index_Char_DD], "Onehot") 

# prepare for imputing median for numerical variables but dropping categorical missing values
median_notdrop <- bind_cols(df[,c("ID", "OUTCOME")], num.impute_median, df[,index_char])
median_drop <- search_complete(median_notdrop)
index_Num_MD <- which(unlist(lapply(median_drop, class))=="numeric")
index_Char_MD <- setdiff(which(unlist(lapply(median_drop, class))=="character"), which(colnames(median_drop)=="ID"|colnames(median_drop)=="OUTCOME"))
median_drop_Label <- encode_category(median_drop[,index_Char_MD], "Label") 
median_drop_Onehot <- encode_category(median_drop[,index_Char_MD], "Onehot") 

# obtain data for 16 scenarios
df_all$Drop_Drop_Label <- bind_cols(drop_drop[,c("ID", "OUTCOME")], drop_drop[,index_Num_DD], drop_drop_Label)
df_all$Drop_Drop_OneHot <- bind_cols(drop_drop[,c("ID", "OUTCOME")], drop_drop[,index_Num_DD], drop_drop_Onehot)
df_all$Drop_Missing_Label <- bind_cols(df[,c("ID", "OUTCOME")], df[,index_num], missing_Label) %>% drop_na()
df_all$Drop_Missing_OneHot <- bind_cols(df[,c("ID", "OUTCOME")], df[,index_num], missing_Onehot) %>% drop_na()
df_all$Drop_Mode_Label <- bind_cols(df[,c("ID", "OUTCOME")], df[,index_num], mode_Label) %>% drop_na()
df_all$Drop_Mode_OneHot <- bind_cols(df[,c("ID", "OUTCOME")], df[,index_num], mode_Onehot) %>% drop_na()
df_all$Drop_Unknown_Label <- bind_cols(df[,c("ID", "OUTCOME")], df[,index_num], unknown_Label) %>% drop_na()
df_all$Drop_Unknown_OneHot <- bind_cols(df[,c("ID", "OUTCOME")], df[,index_num], unknown_Onehot) %>% drop_na()
df_all$Median_Drop_Label <- bind_cols(median_drop[,c("ID", "OUTCOME")], median_drop[,index_Num_MD], median_drop_Label)
df_all$Median_Drop_OneHot <- bind_cols(median_drop[,c("ID", "OUTCOME")], median_drop[,index_Num_MD], median_drop_Onehot)
df_all$Median_Missing_Label <- bind_cols(df[,c("ID", "OUTCOME")], num.impute_median, missing_Label)
df_all$Median_Missing_OneHot <- bind_cols(df[,c("ID", "OUTCOME")], num.impute_median, missing_Onehot)
df_all$Median_Mode_Label <- bind_cols(df[,c("ID", "OUTCOME")], num.impute_median, mode_Label)
df_all$Median_Mode_OneHot <- bind_cols(df[,c("ID", "OUTCOME")], num.impute_median, mode_Onehot)
df_all$Median_Unknown_Label <- bind_cols(df[,c("ID", "OUTCOME")], num.impute_median, unknown_Label)
df_all$Median_Unknown_OneHot <- bind_cols(df[,c("ID", "OUTCOME")], num.impute_median, unknown_Onehot)

saveRDS(df_all, "all_data.RDS")

We found that if we simply dropped all missing values in the cateorical variables, there were no observations in the data. Thus, we try to maximize remaining cells in the data through a heuristic algorithm: search_complete(), which drops rows and columns based on the percentage of missing values until just non-empty cells are remained. This function can be found in custom_functions.R.

The following table shows the dimensions for the data (including ID) we created for all scenarios.

D <- lapply(df_all, dim)
D <- do.call(rbind.data.frame, D)
colnames(D) <- c("Number_of_Observations", "Number_of_Variables")
D <- data.frame(Scenario=names(df_all), D)
datatable(D)

3 All the scenarios in the first year

Since the overarching question examined is “can we quantify the combined effect of decisions made uring the KDDM process on the predictive performance of the model?”, we will explore the six factors summarized in the following.

Categorical Imputation: Drop, Mode, Missing, Unknown
Numerical Imputation: Drop, Median
Encoding: Label, One-Hot
Subsampling: None, Down, Up, SMOTE, ROSE
Feature Selection: FFS, LASSO, RF
Algorithm: LR, LDA, ANN, CART, SVM, RF and XGB

Readers are encouraged to refer our paper for the descriptions of each factor.

3.1 Feature Selection

In the following code chunk, we show how the feature selection approaches: Fast Feature Selection, LASSO, and Random Forests are utilized. If there is an error occurs related to “JAVA_HOME cannot be determined from the registry”, the error is often resolved by installing a Java version (i.e. 64-bit Java or 32-bit Java) that fits to the type of R version that you are using (i.e. 64-bit R or 32-bit R). And set the JAVA_HOME path using one of the commends below. “jre1.8.0_251” needs to be adjusted depends on the folder’s name.

Sys.setenv(JAVA_HOME=‘C:/Program Files/Java/jre1.8.0_251’) for 64-bit version

Sys.setenv(JAVA_HOME=‘C:/Program Files (x86)/Java/jre1.8.0_251’) for 32-bit version

We need to create training and holdout data since the feature selection procedure should be applied to the training data only. We will use IDs to distinguish the training and holdout data. Since we would like to have 5 disjoint holdout datasets for the cross validation (corresponding to 5 training datasets), we sample cases by using ID without replacement and obtatin each 20% of the sample cases for the holdout data.

set.seed <- 1234
all_data_ID <- lapply(df_all, function(x) sample(x$ID, nrow(x)))
holdout_index <- Get_holdout_index(all_data_ID)

Since there will be \(5\times 16=80\) training datasets, we will use a parallel algorithm to conduct the feature selection procedure for these datasets simultaneously. We apply each feature selection method to each dataset individually and the result is stored in a list with the following structure. The list contains 5 sub-lists which contain the results corresponding to the five fold cross validation training sets. Additionally, each of these 5 sub-lists contains 16 elements correponding to the 16 data scenarios. Then we store 3 lists: features_FFS, features_LASSO, features_RF into a list: all_features.

For example, all_features[[“LASSO”]][[3]][[“Median_Drop_Label”]] stores the features after applying LASSO to the third training datasets from the data: imputing medain for missing values in the numerical variables, dropping missing values in the categorical variables, and mapping each level in each categorical variable to an integer.

numCores <- 4 # this depends on the number of cores and the installed memory (RAM) available in the computer
# Fast Feature selection
features_FFS <- vector(mode="list", 5)
startTime <- proc.time()[3] # set up timer

for (i in 1:5){
  cl <- makeCluster(numCores, type="SOCK")
  features_FFS[[i]] <- parSapply(cl, 1:16, select_variables, df_all, "OUTCOME", "FFS", all_data_ID, holdout_index[[i]], 2090)
  stopCluster(cl)
}

runTime = proc.time()[3] - startTime
print(paste("It took", round(runTime), "seconds to obtain all features used using FFS."))

# Lasso Feature selection for Binomial TARGETS
features_LASSO <- vector(mode="list", 5)
startTime <- proc.time()[3] # set up timer

for (i in 1:5){
  cl <- makeCluster(numCores, type="SOCK")
  features_LASSO[[i]] <- parSapply(cl, 1:16, select_variables, df_all, "OUTCOME", "LASSO", all_data_ID, holdout_index[[i]], 2090)
  stopCluster(cl)
}

runTime = proc.time()[3] - startTime
print(paste("It took", round(runTime), "seconds to obtain all features used using LASSO."))

# Random Forest Feature Selection
features_RF <- vector(mode="list", 5)
startTime <- proc.time()[3] # set up timer

for (i in 1:5){
  cl <- makeCluster(numCores, type="SOCK")
  features_RF[[i]] <- parSapply(cl, 1:16, select_variables, df_all, "OUTCOME", "RF", all_data_ID, holdout_index[[i]], 2090)
  stopCluster(cl)
}

runTime = proc.time()[3] - startTime
print(paste("It took", round(runTime), "seconds to obtain all features used using Random Forest."))

all_features <- list(features_FFS, features_LASSO, features_RF)
names(all_features) <- c("FFS", "LASSO", "RF")

For the reference, the approximate run time (in seconds) for each method to obtain the list using a Windows 10 machine with Intel(R) Core(TM) i7 Processor (4 cores, 8GB RAM) is reported below.

Time <- c(374, 2328, 5012)
names(Time) <- c("FFS", "LASSO", "RF")
data.table(t(Time))

3.2 Obtain Results for All Scenarios

In the following code chunk, a parallel algorithm will be used to obtain the prediction performance of each combination of factors and we save the result in the file: all_scenarios.csv. We report five common performance metrics that are suited for two class classification problems: AUC, accuracy, sensitivity, specificity and G-mean where G-mean \(= \sqrt{\mbox{sensitivity} \times \mbox{specificity}}\).

We use a function DOE_function() to conduct this study. This function can be found in the custom_functions.R. The output of this function could be either a vector or a list depends on if the predicted probability for Survival is returned.

fold <- 1:5
feature_selection <- c("FFS", "LASSO","RF")
subsampling <- c("none", "down", "Up", "smote", "rose")
algorithm <- c("glm" , "lda","nnet" ,"rpart" ,"svm", "ranger" , "xgbDART")
#The Subsampling and algorithm names are the names used in the CARET package.


Design <- crossing(fold, imputation_num, imputation_cat, encoding, feature_selection, subsampling, algorithm)

numCores <- 4 # this depends on the number of cores and the installed memory (RAM) available in the computer

cl <- makeCluster(numCores, type="SOCK")

all_scenarios <- parSapply(cl, 1:nrow(Design), DOE_function, Design, df_all, "OUTCOME", all_data_ID, holdout_index, all_features, return.prediction=FALSE, 2020)

stopCluster(cl)

all_scenarios <- t(all_scenarios)
colnames(all_scenarios) <- c("fold", "imputation_num", "imputation_cat", "encoding", "feature_selection", "subsampling", "algorithm", "auc", "sen", "spec", "accu", "gmean")
all_scenarios <- tbl_df(all_scenarios)
all_scenarios[,c("fold", "auc", "sen", "spec", "accu", "gmean")] %<>% mutate_all( as.numeric)

#write.csv(all_scenarios, "../Results/all_scenarios.csv", row.names = FALSE)

3.3 Selecting the best case

Here the performance of algorithms are compared. Since a majority of the “svm” and “GB” cases has not merged to a result after 24 hours implementation over supercomputer, they are not presented in here. Logistic Regression (LR) outperformed all other algorithms in both AUC and G-mean. We use the combination of factors (data preparation, feature selection, resampling, algorithm) for predicting other time stamps. We refer the readers to our paper for the justifications.

full_fact_res<-read.csv("https://raw.githubusercontent.com/transplantation/unos-ht/master/data/rmarkdown/two_stage/best-gmean-repeat-CV.csv")

DT::datatable(full_fact_res)

4 Section: Create the training and holdout datasets

In this snippet, we create the training and holdout datasets at each time point of interest. In order to both validate models on the holdout datasets and use isotonic regression to calibrate survival probabilities, we make sure all patients selected in the 10th year holdout dataset were included in previous time points. This is not a necessary procedure for either prediction purpose or calibration of survival probabilities. However, we want to evaluate calibration with isotonic regression for the patients that their survival is predicted in all time stamps. Therefore, we need to use the patients IDs as an indicator for checking their data availability in all the time stamps.

We use a list object keep_NA to record the following IDs:

IDs for the dataset at each time point of interest
IDs for the holdout dataset at each time point of interest
IDs for the training dataset at each time point of interest

# find the indecies for categorical and numerical variables in the data
index_char <- setdiff(which(unlist(lapply(df_iso, class))=="character"), which(colnames(df_iso)=="ID"|colnames(df_iso)=="GTIME"|colnames(df_iso)=="GSTATUS"))

index_num <- setdiff(which(unlist(lapply(df_iso,class))=="numeric"),
 which(colnames(df_iso)=="ID"|colnames(df_iso)=="GTIME"|colnames(df_iso)=="GSTATUS"))                    

char.impute_unknown <- df_iso[,index_char] %>% mutate_all(~ifelse(is.na(.x), "Unknown", .x))

iso.char.impute_unknown <- df_iso[,index_char] %>% mutate_all(~ifelse(is.na(.x), "Unknown", .x)) 
iso.num.impute_median <- df_iso[,index_num] %>% mutate_all(~ifelse(is.na(.x), median(.x, na.rm = TRUE), .x)) 
iso.unknown_Onehot <- encode_category(iso.char.impute_unknown, "Onehot")

df_iso<-bind_cols(df_iso[,c("ID","GSTATUS","GTIME")], iso.num.impute_median, iso.unknown_Onehot)


# here we made 11 consecutive dependent variables, month1, year1, year2, year3, ..., year10

df_iso %<>%
mutate(year0 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 1/12)) %>%
mutate(year1 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 1))    %>%
mutate(year2 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 2))    %>%
mutate(year3 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 3))    %>%
mutate(year4 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 4))    %>%
mutate(year5 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 5))    %>%
mutate(year6 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 6))    %>%
mutate(year7 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 7))    %>%
mutate(year8 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 8))    %>%
mutate(year9 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 9))    %>%
mutate(year10 = recodeGSTATUS(gStatus = GSTATUS, gTime = GTIME, targetYear = 10))  %>%
select(-c(GSTATUS, GTIME))
iso_vars=names(df_iso)

# We keep IDs of non NA observations in each time stamps
keep_NA <- rep(list(NA), 33)
keep_NA[1:11]<-lapply(1:11,function(x){pull(df_iso[!is.na(df_iso[paste0("year", (x-1))]),"ID"],ID)})
names(keep_NA)[1:11] <- paste0("ID",seq(0,10))

all_sizes <- unname(unlist(lapply(keep_NA[1:11], length)))
test_sizes <- round(all_sizes*0.2)
train_sizes <- all_sizes - test_sizes

# We select the holdout sets that contain the year 10 test data
# Then we exclude those IDs from each year to find IDs for other training sets
set.seed(2019)
keep_NA[[12]] <- sample(keep_NA$ID10,test_sizes[11])
names(keep_NA)[12] <- "ID_holdout10"

keep_NA[13:22]<-lapply(2:11,function(x){
  c(keep_NA$ID_holdout10, sample(setdiff(keep_NA[[(12-x)]], keep_NA$ID_holdout10), (test_sizes[(12-x)]-length(keep_NA$ID_holdout10))))
})
names(keep_NA)[13:22] <- paste0("ID_holdout", seq(9,0))


keep_NA[23:33]<-lapply(1:11,function(x){
  setdiff(keep_NA[[x]], keep_NA[[paste0("ID_holdout",(x-1))]])})
names(keep_NA)[23:33] <- paste0("ID_train", seq(0,10))

time_stamps<-c("year0","year1","year2","year3","year4","year5", "year6","year7","year8","year9","year10")


# making the data ready for running feature selection over super computer
server_data<-list()

server_data<-lapply(1:11,function(x){df_iso %>% dplyr::filter(ID %in% keep_NA[[paste0("ID_train", x-1)]])%>%
                      select(c(names(.)[!names(.) %in% c("ID",time_stamps)],paste0("year",x-1)))})

5 Section: Multi Time Points Feature Selection

The following code is what is performed over super computer to achive feature selection with LASSO The LASSO feature selection algorithm is loaded

cl <- makeCluster(min(4,parallel::detectCores()), type="SOCK")
time.begin <- proc.time()[3]

LASSO <- parSapply(cl, 1:11, select_variables_iso,data=server_data,method="LASSO",
                          folds=5,alpha=1,seed=110) 

time.window <- proc.time()[3] - time.begin

stopCluster(cl)

6 Section: Multi Time Points training and prediction

In the following code snippet First we identify train and test set for each time stamp. Then we train logistic regressions and test them over all the identified.

train_data<-list()
holdOut_data<-list()
exclud <- c(paste0("year",seq(0,10)), "ID")
exclud_h <- c(paste0("year",seq(0,10)))
for (i in 1:11){
  train_data[[i]]<-df_iso[df_iso$ID %in% keep_NA[[paste0("ID_train", i-1)]],
               !names(df_iso) %in% exclud[exclud!= paste0("year",i-1)] ]

  
  holdOut_data[[i]]<-df_iso[df_iso$ID %in% keep_NA[[paste0("ID_holdout", i-1)]],
                            !names(df_iso) %in% exclud_h[exclud_h!= paste0("year",i-1)]]

}


cl <- makeCluster(min(4,parallel::detectCores()), type="SOCK")
time.begin <- proc.time()[3]

pred_LR <- parSapply(cl, 1:11,train_iso ,t_data=train_data,h_data=holdOut_data,features=LASSO,
                          folds=5,resampl_meth="up",alg_used="glm",seed0=2019) 

time.window <- proc.time()[3] - time.begin

stopCluster(cl)

7 Section: Important features & prediction performance before isotonic regression

The following codes demonstrates a report about importance of the features used for prediction & prediction performance in each time stamp. “LR_models” contains information of all the logsitic regression models that have been trained over the all the time stamps.

LR_models<-readRDS(url("https://github.com/transplantation/unos-ht/raw/master/data/pred_LR.rds"))
# here we load name of all the variables in the train set, they are loaded here to create a report them in each time stamps
iso_vars<-readRDS(url("https://github.com/transplantation/unos-ht/raw/master/data/iso_vars.rds"))
 
exclud <- c(paste0("year",seq(0,10)), "ID")
vars<-iso_vars[!iso_vars %in% exclud]

varimp_years<-as.data.frame(matrix(0, ncol = 12, nrow = length(vars)))
colnames(varimp_years)<-c("vars",paste0("year",seq(0,10)))
varimp_years$vars<-vars

perf_years<-as.data.frame(matrix(0, ncol = 13, nrow = 5))
colnames(perf_years)<-c("measure",paste0("year",seq(0,10)),"all_mean")
perf_years$measure<-c("AUC","Specificity","Sensitivity","Accuracy","Geometic_Mean")
 


lapply(0:10, function(x) {
   varimp_years[varimp_years$vars %in% LR_models[4,x+1]$var_imp$vars,paste0("year",x)]<<-1
   perf_years[paste0("year",x)]<<-LR_models[1,x+1]$Performance
   }
 )


# prediction performance
perf_years$all_mean<-rowMeans(perf_years[paste0("year",seq(0,10))])
perf_years[,2:13]<-round(perf_years[,2:13],4)

# feature inportance
varimp_years$no_years<-rowSums(varimp_years[paste0("year",seq(0,10))])
varimp_years<-varimp_years[which(varimp_years$no_years!=0),]

The following table shows if a variable selected by the internal mechanism of the training algorithms in each time stamp. ‘1’ means if a variable is selected and ‘0’ means otherwise.In addition, the total apperance of a variable in all the time stamps is counted.

DT::datatable(varimp_years)

Performance of prediction for the patients that their data is available in all the time stamps is presented below. The performance of the first year is slightly different from the earlier prediction since the holdout set is different.

DT::datatable(perf_years)

8 Section: Calibrating the predictios with isotonic regression

The following code demonstrates the application of isotonic regression. Then performance of predictions before and after isotonic regression is compared.

predictions<-list()
# holdout ID of patients that their survival is known in year 10, should be existed in all previous time stamps but if their data is not known (e.g. year5) their data survival status is not known in the following years
#all_TARGETS: actual survival status of patients in each time stamp
#all_TARGETS: predicted survival status of the patients
#all_preds: predicted probability of survival in each time stamp

in_all_h<-LR_models[2,11]$Predicted$ID
all_preds<-as.data.frame(matrix(0, ncol = 11, nrow = length(in_all_h)))
all_TARGETS<-as.data.frame(matrix(0, ncol = 11, nrow = length(in_all_h)))
colnames(all_preds)<-c(paste0("year",seq(0,10)))
colnames(all_TARGETS)<-c(paste0("year",seq(0,10)))
all_TARGETS_pred<-all_TARGETS


lapply(0:10, function(x) {
  predictions[[x+1]]<<-LR_models[2,x+1]$Predicted[which(LR_models[2,x+1]$Predicted$ID %in%
                                                          in_all_h),]
  if(sum(predictions[[x+1]]$ID==in_all_h)==length(in_all_h)){
    all_preds[paste0("year",x)]<<-predictions[[x+1]]$Probability
    all_TARGETS[paste0("year",x)]<<-predictions[[x+1]]$TARGET_raw[paste0("year",x)]
    all_TARGETS_pred[paste0("year",x)]<<-predictions[[x+1]]$glm
    return()
  }
})


# here we apply isotonic regression 
# all_preds2: calibrated survival probability
# all_TARGETS_pred2: survival status after calibrating survival probability

# here we calibrate the survival probability

all_preds2<- all_preds %>% rowwise() %>%
  mutate(output = list(isoreg(rev(as.numeric(c_across(year0:year10))))$yf)) %>%
  tidyr::unnest_wider(output)%>%
  select(., -starts_with("year")) %>% as_tibble(.)%>% rev(.) %>% 
  rename_(.dots = setNames(names(.), paste0("year", 0:10)))%>% 
  mutate(monotonic=all(.==cummin(.))) 



# here we identify survival status based on the calibrated probabilities
all_TARGETS_pred2= all_preds2 %>%
  select(paste0("year",0:10)) %>%
  mutate_all(funs(replace(., . >= 0.5, "Survival")))%>%
  mutate_all(funs(replace(., . < 0.5, "Death"))) %>% mutate_if(is.character,as.factor)
all_TARGETS_pred2<-as.data.frame(all_TARGETS_pred2)

# here we calculate performance of prediction for the patients that their data is available in all time stamps before/after applying isotonic regression

perf_before_iso<-data.frame(matrix(NA, nrow=5, ncol=0))
perf_before_iso$Measures<-c("Accuracy","Specificity", "Sensitivity","AUC","G-Mean" )
perf_after_iso<-perf_before_iso

# calculating prediction performance before applying isotonic regression
perf_before_iso<-as.data.frame(lapply(0:10, function(x) {
  perf_before_iso %<>% mutate(!!paste0("year",x) := c(
    (as.data.frame(confusionMatrix(all_TARGETS_pred[,paste0("year",x)],
                                   as.factor(all_TARGETS[,paste0("year",x)]))$overall))[1,],
    caret::specificity(all_TARGETS_pred[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)])), 
    caret::sensitivity(all_TARGETS_pred[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)])),
    AUC::auc(roc(all_preds[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)]))),
    sqrt(caret::specificity(all_TARGETS_pred[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)]))*
           caret::sensitivity(all_TARGETS_pred[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)])))
  ))
  return(perf_before_iso[paste0("year",x) ])
}))

perf_before_iso %<>% mutate(Mean = rowMeans(.[,colnames(perf_before_iso)])) %>% 
  add_rownames("Performance Measure")%>% mutate(`Performance Measure` = c("Accuracy","Specificity", "Sensitivity","AUC","G-Mean"))

# calculating prediction performance after applying isotonic regression

perf_after_iso<-as.data.frame(lapply(0:10, function(x) {
  perf_after_iso %<>% mutate(!!paste0("year",x) := c(
    (as.data.frame(confusionMatrix(all_TARGETS_pred2[,paste0("year",x)],
                                 all_TARGETS[,paste0("year",x)])$overall))[1,],
    caret::specificity(all_TARGETS_pred2[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)])), 
    caret::sensitivity(all_TARGETS_pred2[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)])),
    AUC::auc(roc(pull(all_preds2, paste0("year",x)),as.factor(all_TARGETS[,paste0("year",x)]))),   
    sqrt(caret::specificity(all_TARGETS_pred2[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)]))*
         caret::sensitivity(all_TARGETS_pred2[,paste0("year",x)],as.factor(all_TARGETS[,paste0("year",x)])))
    
  ))
  return(perf_after_iso[paste0("year",x) ])
}))

perf_after_iso %<>% mutate(Mean = rowMeans(.[,colnames(perf_after_iso)])) %>% 
  add_rownames("Performance Measure")%>% mutate(`Performance Measure` = c("Accuracy","Specificity", "Sensitivity","AUC","G-Mean"))

9 Results

In this section, we report the results obtained from the calibration of the predicted survivals.

Here you can find performance of predictions before/after isotonic regression for the patients that their data is available in all the timestamps.

Here the results before isotonic regression is presented:

perf_before_iso[,2:13]<-round(perf_before_iso[,2:13],3)
DT::datatable(perf_before_iso)

Here the results after isotonic regression is presented:

perf_after_iso[,2:13]<-round(perf_after_iso[,2:13],3)
DT::datatable(perf_after_iso)

Here the change after isotonic regression is presented:

perf_delta=perf_after_iso[,1:12]
perf_delta[,2:12]<-round(perf_after_iso[,2:12]-perf_before_iso[,2:12],3)
DT::datatable(perf_delta)

9.1 Survival Probabilites for Random Sample Patients before and after Isotonic Regression

9.2 ROC Curves

Month 1

Year 1

Year 2

Year 3

Year 4

Year 5

Year 6

Year 7

Year 8

Year 9

Year 10

9.3 Observed vs forcasted survival probabities before and afer Applying Isotonic Regression

In this section, we exam how good our models are by checking the survival probability before and after applying Isotonic Regression.

Month 1

Year 1

Year 2

Year 3

Year 4

Year 5

Year 6

Year 7

Year 8

Year 9

Year 10

10 Interactive App

We have created an interactive app that presents two modules for performing the analysis:

Manual Entry, where users can insert the values of predictor variables using several text boxes.
CSV Entry, where users can upload the values of predictor variables using a comma seperated variable (CSV) file.

Visit the app: here.

Figure 2 presents the workflow of the app. The user does not need to have a background in programming nor machine learning to operate the app.

Figure 2: The workflow of the app.

References

Chen, Ching Kit, Cedric Manlhiot, Seema Mital, Steven M Schwartz, Glen S Van Arsdell, Christopher Caldarone, Brian W McCrindle, and Anne I Dipchand. 2018. “Prelisting Predictions of Early Postoperative Survival in Infant Heart Transplantation Using Classification and Regression Tree Analysis.” Pediatric Transplantation 22 (2): e13105.

Dag, Ali, Asil Oztekin, Ahmet Yucel, Serkan Bulur, and Fadel M Megahed. 2017. “Predicting Heart Transplantation Outcomes Through Data Analytics.” Decision Support Systems 94: 42–52.

Dag, Ali, Kazim Topuz, Asil Oztekin, Serkan Bulur, and Fadel M Megahed. 2016. “A Probabilistic Data-Driven Framework for Scoring the Preoperative Recipient-Donor Heart Transplant Survival.” Decision Support Systems 86: 1–12.

De Jonge, Edwin, and Mark Van Der Loo. 2013. An Introduction to Data Cleaning with R. Statistics Netherlands Heerlen.

Guven, Goksel, Milos Brankovic, Alina A Constantinescu, Jasper J Brugts, Dennis A Hesselink, Sakir Akin, Ard Struijs, et al. 2018. “Preoperative Right Heart Hemodynamics Predict Postoperative Acute Kidney Injury After Heart Transplantation.” Intensive Care Medicine 44 (5): 588–97.

Max, Kuhn. 2020. Using Recipes with Train. https://topepo.github.io/caret/using-recipes-with-train.html.

Medved, Dennis, Mattias Ohlsson, Peter Höglund, Bodil Andersson, Pierre Nugues, and Johan Nilsson. 2018. “Improving Prediction of Heart Transplantation Outcome Using Deep Learning Techniques.” Scientific Reports 8 (1): 1–9.

Parham, Peter. 2014. The Immune System. Garland Science.

Shmueli, Galit, and others. 2010. “To Explain or to Predict?” Statistical Science 25 (3): 289–310.

Vega, Eleanor, Jacob Schroder, and Alina Nicoara. 2017. “Postoperative Management of Heart Transplantation Patients.” Best Practice & Research Clinical Anaesthesiology 31 (2): 201–13.

Weisdorf, Daniel, Stephen Spellman, Michael Haagenson, Mary Horowitz, Stephanie Lee, Claudio Anasetti, Michelle Setterholm, et al. 2008. “Classification of HLA-Matching for Retrospective Analysis of Unrelated Donor Transplantation: Revised Definitions to Predict Survival.” Biology of Blood and Marrow Transplantation 14 (7): 748–58.

Zhang, Kan, Richard Sheu, Nicole M Zimmerman, Andrej Alfirevic, Shiva Sale, A Marc Gillinov, and Andra E Duncan. 2019. “A Comparison of Global Longitudinal, Circumferential, and Radial Strain to Predict Outcomes After Cardiac Surgery.” Journal of Cardiothoracic and Vascular Anesthesia 33 (5): 1315–22.

A Two-Stage Machine Learning Framework to Predict Heart Transplantation Survival Probabilities over Time with a Monotonic Probability Constraint

(Authors’ names blinded for peer review)

19 July, 2020

1 Data Exploration

1.1 Loading the UNOS Dataset

1.2 Exploring the UNOS Dataset

Fraction Data Missing

Missing Data (Sample Columns)

Missing Data (Printout for All Variables)

Column Types

2 Data Preperation

2.1 Technically Correct Variables

2.2 Rowwise Operations: Additional Filters

2.3 Columnwise Operations: Dropping Variables

2.4 Creating Variables based on the Literature and Our Exploration of the Data

2.5 Reducing Number of Levels for Categorical Columns

2.5.1 Number of Distinct Values within Each Categorical Variable

2.5.2 Grouping Operations

2.6 Dealing with Missing Values & Encoding Strategy

3 All the scenarios in the first year

3.1 Feature Selection

3.2 Obtain Results for All Scenarios

3.3 Selecting the best case

4 Section: Create the training and holdout datasets

5 Section: Multi Time Points Feature Selection

6 Section: Multi Time Points training and prediction

7 Section: Important features & prediction performance before isotonic regression

8 Section: Calibrating the predictios with isotonic regression

9 Results

9.1 Survival Probabilites for Random Sample Patients before and after Isotonic Regression

9.2 ROC Curves

Month 1

Year 1

Year 2

Year 3

Year 4

Year 5

Year 6

Year 7

Year 8

Year 9

Year 10

9.3 Observed vs forcasted survival probabities before and afer Applying Isotonic Regression

Month 1

Year 1

Year 2

Year 3

Year 4

Year 5

Year 6

Year 7

Year 8

Year 9

Year 10

10 Interactive App

References