{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"### Import dependencies\nlibrary(tidyverse)\nlibrary(reshape2)","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2022-05-08T18:02:24.396337Z","iopub.execute_input":"2022-05-08T18:02:24.398708Z","iopub.status.idle":"2022-05-08T18:02:25.782358Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Read the data into memory\nraw_data <- read.table(\"../input/aacc-2022-predicting-pthrp-results//AACC_Kaggle_2022_PTHrP_input_data.txt\", header = T, sep = \"\\t\")\ntrain_target <- read.table(\"../input/aacc-2022-predicting-pthrp-results//AACC_Kaggle_2022_PTHrP_train_results.txt\", header = T, sep = \"\\t\")\ntest_target <- read.table(\"../input/aacc-2022-predicting-pthrp-results//AACC_Kaggle_2022_PTHrP_test_results_masked.txt\", header = T, sep = \"\\t\")","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:02:25.784821Z","iopub.execute_input":"2022-05-08T18:02:25.816755Z","iopub.status.idle":"2022-05-08T18:02:35.971807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Summarize the data\n\nggplot(as.data.frame(sort(table(raw_data$TASK_ASSAY), decreasing = T))[1:10,], aes(x = Var1, y = Freq)) + geom_bar(stat = \"identity\") + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))\nggplot(as.data.frame(sort(table(raw_data$PATIENT_ID), decreasing = T))[1:10,], aes(x = Var1, y = Freq)) + geom_bar(stat = \"identity\") + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) \nggplot(as.data.frame(sort(table(raw_data$ORDERABLE), decreasing = T))[1:10,], aes(x = Var1, y = Freq)) + geom_bar(stat = \"identity\") + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))\nggplot(as.data.frame(sort(table(raw_data$ORDER_PROVIDER), decreasing = T))[1:10,], aes(x = Var1, y = Freq)) + geom_bar(stat = \"identity\") + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:02:35.975707Z","iopub.execute_input":"2022-05-08T18:02:35.977832Z","iopub.status.idle":"2022-05-08T18:02:38.197738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### PTHrP results should be correlated to other labs involved in calcium homeostasis, so we will restrict our dataset to only include some common calcium-related labs\n\nlabs_of_interest <- c(\"Calcium\", \"AST\", \"ALT\", \"Albumin\", \"PTH Intact\", \"Magnesium\", \"Vit D 1,25\", \"Vit D 25 OH\", \"Phos Plas\", \"Alk Phos\", \"Bone Alk Phos\", \"Calcium UConc\", \"Creatinine\")\n\ndata_filtered <- raw_data %>% filter(TASK_ASSAY %in% labs_of_interest)","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:02:38.200141Z","iopub.execute_input":"2022-05-08T18:02:38.201554Z","iopub.status.idle":"2022-05-08T18:02:38.357819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Our result values include some values above and below the limit of quantitation, so we will just reassign the inequalities to be at those limits to make our data completely numeric\n\ndata_filtered$RESULT_VAL <- gsub(data_filtered$RESULT_VAL, pattern = \"<\", replacement = \"\")\ndata_filtered$RESULT_VAL <- gsub(data_filtered$RESULT_VAL, pattern = \">\", replacement = \"\")\n\ndata_filtered$RESULT_VAL <- as.numeric(data_filtered$RESULT_VAL)\ndata_filtered <- data_filtered[which(!is.na(data_filtered$RESULT_VAL)),]","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:02:38.360346Z","iopub.execute_input":"2022-05-08T18:02:38.361730Z","iopub.status.idle":"2022-05-08T18:02:38.522725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Now we make our table. This will include the most recent value for each of our labs of interest, as well as min, max and mean.\n\nIDs <- unique(raw_data$PATIENT_ID)\n\ndata_filtered <- data_filtered %>% group_by(PATIENT_ID) %>% arrange(PERFORMED_DT_TM)\n\n###### Create Cast Data.Frame ######\nPTHrp_melt <- as.data.frame(matrix(NA, nrow = 1, ncol = 3))\ncolnames(PTHrp_melt) <- c(\"PATIENT_ID\", \"TASK_ASSAY\", \"RESULT_VAL\")  \n\nfor(i in IDs){\n  tmp_data_filtered <- data_filtered %>% filter(PATIENT_ID == i) %>% group_by(TASK_ASSAY)\n  PATIENT_ID <- i\n    \n  maxes <- tmp_data_filtered %>% summarise(RESULT_VAL = max(RESULT_VAL))\n  maxes$TASK_ASSAY <- paste(maxes$TASK_ASSAY, \"max\", sep = \"_\")\n  maxes$PATIENT_ID <- PATIENT_ID\n   \n  mins <-  tmp_data_filtered %>% summarise(RESULT_VAL = min(RESULT_VAL))\n  mins$TASK_ASSAY <- paste(mins$TASK_ASSAY, \"min\", sep = \"_\")\n  mins$PATIENT_ID <- PATIENT_ID\n    \n  means <-  tmp_data_filtered %>% drop_na(RESULT_VAL) %>% summarise(RESULT_VAL = mean(RESULT_VAL))\n  means$TASK_ASSAY <- paste(means$TASK_ASSAY, \"mean\", sep = \"_\")\n  means$PATIENT_ID <- PATIENT_ID\n    \n  tmp_melt <- tmp_data_filtered %>% group_by(TASK_ASSAY) %>% select(PATIENT_ID, TASK_ASSAY, RESULT_VAL) %>% slice_tail(n = 1)\n    \n  tmp_melt <- rbind(tmp_melt, mins, means, maxes)\n    \n  PTHrp_melt <- rbind(PTHrp_melt, tmp_melt)\n}\n\nPTHrp_cast <- dcast(PTHrp_melt, formula = PATIENT_ID ~ TASK_ASSAY, value.var = \"RESULT_VAL\")\nPTHrp_cast <- as.data.frame(lapply(PTHrp_cast, as.numeric))","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:02:38.525241Z","iopub.execute_input":"2022-05-08T18:02:38.526625Z","iopub.status.idle":"2022-05-08T18:05:13.624343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Now that we have our tabular data, we will write two versions of input files to use for our machine learning models. One will be the full table, while the other will filter out any column that is so sparsely populated that it won't be imputable in any meaningful way. Arbitrarily, I've chosen a cut-off of 20%.\n\nNA_col_percent <- as.numeric(lapply(PTHrp_cast, function(x) sum(is.na(x))))/nrow(PTHrp_cast)\n\nPTHrp_cast <- PTHrp_cast[,which(colnames(PTHrp_cast) != \"NA.\")]\n\nPTHrp_cast_dense <- PTHrp_cast[,which(NA_col_percent < 0.2)]\n\nwrite.table(PTHrp_cast_dense, \"PTHrp_ML_input_dense.txt\", sep = \"\\t\", row.names = F, quote = F)\nwrite.table(PTHrp_cast, \"PTHrp_ML_input_sparse.txt\", sep = \"\\t\", row.names = F, quote = F)","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:13.628615Z","iopub.execute_input":"2022-05-08T18:05:13.631076Z","iopub.status.idle":"2022-05-08T18:05:13.840664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Now we can build and train our XGBoost model on this tabular data. \n\nlibrary(tidymodels)\nlibrary(themis)\n\ndata <- PTHrp_cast\n\nset.seed(12345)","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:13.843679Z","iopub.execute_input":"2022-05-08T18:05:13.845063Z","iopub.status.idle":"2022-05-08T18:05:17.055893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Add a calcium correction to the most recent calcium results\n\ndata$corrCalcium <- (0.8 * (4.3 - data$Albumin)) + data$Calcium","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:17.059807Z","iopub.execute_input":"2022-05-08T18:05:17.061392Z","iopub.status.idle":"2022-05-08T18:05:17.074799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Summarize the data\n\nfeature_summaries <- t(as_tibble(lapply(data, summary)))\ncolnames(feature_summaries) <- c(\"Min\", \"1st Quartile\", \"Median\", \"Mean\", \"3rd Quartile\", \"Max\", \"NA\")\nfeature_summaries[,'NA'] <- feature_summaries[,'NA']/nrow(data)\nfeature_summaries","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:17.077412Z","iopub.execute_input":"2022-05-08T18:05:17.078901Z","iopub.status.idle":"2022-05-08T18:05:17.214001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Make train and test splits, including some data from the PTHrP orders\n\ntrain <- data %>% filter(PATIENT_ID %in% train_target$PATIENT_ID)\ntrain <- left_join(train, train_target %>% select(PATIENT_ID, MEDICAL_SERVICE, ENCOUNTER_TYPE, NORMALCY))\nrownames(train) <- train$PATIENT_ID\ntrain <- train %>% select(-PATIENT_ID)\ntrain <- train %>% rename(target = NORMALCY)\n\ntest <- data %>% filter(PATIENT_ID %in% test_target$PATIENT_ID)\ntest <- left_join(test, test_target %>% select(PATIENT_ID, MEDICAL_SERVICE, ENCOUNTER_TYPE, NORMALCY))\nrownames(test) <- test$PATIENT_ID\ntest <- test %>% select(-PATIENT_ID)\ntest <- test %>% rename(target = NORMALCY)\n","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:17.216520Z","iopub.execute_input":"2022-05-08T18:05:17.217918Z","iopub.status.idle":"2022-05-08T18:05:17.358775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Make our recipe. This will include converting the categorical data to dummy variables, removing low-variance features, removing any column with more than 20% missing data, a center and scaling step for normalization, an imputation technique (bagged trees), removing any correlated or linear combinations of variables, then a synthetic upsampling technique (AdaSyn) to balance the classes. \n\nrecipe <- recipe(target ~ ., data = train) %>% \n  step_dummy(all_nominal_predictors()) %>%\n  step_nzv(all_numeric_predictors()) %>% \n  step_filter_missing(all_predictors(), threshold = 0.2) %>%\n  step_normalize(all_numeric_predictors()) %>%\n  step_impute_bag(all_numeric_predictors()) %>%\n  step_corr(all_numeric_predictors()) %>% \n  step_lincomb(all_numeric_predictors()) %>% \n  step_adasyn(target, over_ratio = 1) %>%\n  prep()","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:17.362555Z","iopub.execute_input":"2022-05-08T18:05:17.364091Z","iopub.status.idle":"2022-05-08T18:05:47.432060Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Apply preprocessing steps, set up 5-fold cross-validation (repeated 5 times) and update our recipe\n\ntrain_preprocessed <- bake(recipe, NULL)\ntest_preprocessed <- bake(recipe, test)\ncv <- vfold_cv(train_preprocessed, v = 5, repeats = 5, strata = target)\n\nrecipe <- recipe(target ~ ., data = train_preprocessed) ","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:47.434664Z","iopub.execute_input":"2022-05-08T18:05:47.466501Z","iopub.status.idle":"2022-05-08T18:05:47.925093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Define models  \n\nxg <- boost_tree(trees = tune(), tree_depth = tune(), mtry = tune(), learn_rate = tune()) %>% set_engine(\"xgboost\") %>% set_mode(\"classification\") \n\nxg_grid <- grid_latin_hypercube(trees(), tree_depth(), finalize(mtry(), cv), learn_rate(), size = 30)\n\nmodels <- c(\"xg\")\n\n## Define workflows\n\nxg_workflow <- workflow() %>% add_recipe(recipe) %>% add_model(xg)","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:47.928824Z","iopub.execute_input":"2022-05-08T18:05:47.930332Z","iopub.status.idle":"2022-05-08T18:05:47.980357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Tuning our hyperparameters based on the tune() assignments in our model definition, then finalize the workflow with the optimized values\n\nxg_tune_results <- xg_workflow %>% tune_grid(resamples = cv, grid = xg_grid) #This may take a very long time. Reduce the size in xg_grid, the number of features in your training data, or the number of tuning parameters to speed this up.\nxg_best <- xg_tune_results %>% select_best(metric = 'roc_auc')\nxg_workflow <- xg_workflow %>% finalize_workflow(xg_best)\n","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:05:47.983022Z","iopub.execute_input":"2022-05-08T18:05:47.984567Z","iopub.status.idle":"2022-05-08T18:29:36.767312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Fit our model to the training set, then apply it to the test set -- this code is written to allow for multiple models to be run, but we only have our XGBoost, feel free to add more models and workflows, then append the new models to the 'models' list.\n\nfor(model in models){\nassign(paste(model, \"fit\", sep = \"_\"), eval(parse(text = paste(model, \"workflow\", sep = \"_\"))) %>% fit(train_preprocessed))\n}\n\nfor(model in models){\nassign(paste(model, \"test_predictions\", sep = \"_\"), eval(parse(text = paste(model, \"fit\", sep = \"_\"))) %>% predict(test_preprocessed))\n}","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:29:36.769668Z","iopub.execute_input":"2022-05-08T18:29:36.771030Z","iopub.status.idle":"2022-05-08T18:29:40.968627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Convert our predictions to the correct submission format\n\nsubmission <- cbind(rownames(test), xg_test_predictions)\ncolnames(submission) <- c(\"PATIENT_ID\", \"PREDICTION\")\nsubmission <- submission %>% mutate(PREDICTION = case_when(PREDICTION == \"Normal\" ~ '0', PREDICTION == \"Abnormal\" ~ '1'))\nwrite.csv(submission, \"submission_file.csv\", quote = F, row.names = F)","metadata":{"execution":{"iopub.status.busy":"2022-05-08T18:29:40.971055Z","iopub.execute_input":"2022-05-08T18:29:40.972375Z","iopub.status.idle":"2022-05-08T18:29:40.992545Z"},"trusted":true},"execution_count":null,"outputs":[]}]}