ML system (part 2).Rmd

---
title: "A ML system (Part 2)"
author: "Guiquan"
date: "2022/1/30"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Part 2 of XGBoost-top-feature work: Development, validation and stress tests of the ML system which integrated four submodels.

# Superiority analysis: Delong test.
## DeLong test on strategy models.
```{r}
library(pROC)
# PORSMs
roc.test(roc_porsm_top, roc_porsm_age, method = "delong")
roc.test(roc_porsm_top, roc_porsm_fsh, method = "delong")
roc.test(roc_porsm_top, roc_porsm_afc, method = "delong")
roc.test(roc_porsm_top, roc_porsm_amh, method = "delong")
roc.test(roc_porsm_top, roc_porsm_agefsh, method = "delong")
roc.test(roc_porsm_top, roc_porsm_ageafc, method = "delong")
roc.test(roc_porsm_top, roc_porsm_ageamh, method = "delong")
roc.test(roc_porsm_top, roc_porsm_ageamhfsh, method = "delong")
roc.test(roc_porsm_top, roc_porsm_ageamhafc, method = "delong")
roc.test(roc_porsm_top, roc_porsm_ageamhafcfsh, method = "delong")
roc.test(roc_porsm_top, roc_porsm_complete, method = "delong")
# HORSMs
roc.test(roc_horsm_top, roc_horsm_age, method = "delong")
roc.test(roc_horsm_top, roc_horsm_fsh, method = "delong")
roc.test(roc_horsm_top, roc_horsm_afc, method = "delong")
roc.test(roc_horsm_top, roc_horsm_amh, method = "delong")
roc.test(roc_horsm_top, roc_horsm_agefsh, method = "delong")
roc.test(roc_horsm_top, roc_horsm_ageafc, method = "delong")
roc.test(roc_horsm_top, roc_horsm_ageamh, method = "delong")
roc.test(roc_horsm_top, roc_horsm_ageamhfsh, method = "delong")
roc.test(roc_horsm_top, roc_horsm_ageamhafc, method = "delong")
roc.test(roc_horsm_top, roc_horsm_ageamhafcfsh, method = "delong")
roc.test(roc_horsm_top, roc_horsm_complete, method = "delong")
```
## DeLong test on diagnostic models.
```{r}
# PORDMs
roc.test(roc_pordm_top, roc_pordm_age, method = "delong")
roc.test(roc_pordm_top, roc_pordm_fsh, method = "delong")
roc.test(roc_pordm_top, roc_pordm_afc, method = "delong")
roc.test(roc_pordm_top, roc_pordm_amh, method = "delong")
roc.test(roc_pordm_top, roc_pordm_agefsh, method = "delong")
roc.test(roc_pordm_top, roc_pordm_ageafc, method = "delong")
roc.test(roc_pordm_top, roc_pordm_ageamh, method = "delong")
roc.test(roc_pordm_top, roc_pordm_ageamhfsh, method = "delong")
roc.test(roc_pordm_top, roc_pordm_ageamhafc, method = "delong")
roc.test(roc_pordm_top, roc_pordm_ageamhafcfsh, method = "delong")
roc.test(roc_pordm_top, roc_pordm_complete, method = "delong")
# HORDMs
roc.test(roc_hordm_top, roc_hordm_age, method = "delong")
roc.test(roc_hordm_top, roc_hordm_fsh, method = "delong")
roc.test(roc_hordm_top, roc_hordm_afc, method = "delong")
roc.test(roc_hordm_top, roc_hordm_amh, method = "delong")
roc.test(roc_hordm_top, roc_hordm_agefsh, method = "delong")
roc.test(roc_hordm_top, roc_hordm_ageafc, method = "delong")
roc.test(roc_hordm_top, roc_hordm_ageamh, method = "delong")
roc.test(roc_hordm_top, roc_hordm_ageamhfsh, method = "delong")
roc.test(roc_hordm_top, roc_hordm_ageamhafc, method = "delong")
roc.test(roc_hordm_top, roc_hordm_ageamhafcfsh, method = "delong")
roc.test(roc_hordm_top, roc_hordm_complete, method = "delong")
```

# Calculate SMs' and DMs' Brier score on external validation cohort.
```{r}
library(sjmisc)
# Specifications of bs function.
brier_score <- function(preds, obs) {
  mean((obs - preds)^2)
}
# Strategy models.
preds_porsm <- predict(porsm_top_final, porsm_top_val, type = "prob") %>% .[[".pred_Yes"]] 
obs_porsm <- porsm_top_val %>% select(POR) %>% rec(., rec = "No = 0; Yes = 1") %>% .[["POR_r"]] %>%
  as.character() %>% as.numeric()
brier_score(obs_porsm, preds_porsm)

preds_horsm <- predict(horsm_top_final, horsm_top_val, type = "prob") %>% .[[".pred_Yes"]] 
obs_horsm <- horsm_top_val %>% select(HOR) %>% rec(., rec = "No = 0; Yes = 1") %>% .[["HOR_r"]] %>%
  as.character() %>% as.numeric()
brier_score(obs_horsm, preds_horsm)

# Diagnosis models.
preds_pordm <- predict(pordm_top_final, pordm_top_val, type = "prob") %>% .[[".pred_Yes"]] 
obs_pordm <- pordm_top_val %>%  
  select(POR) %>% rec(., rec = "No = 0; Yes = 1") %>% .[["POR_r"]] %>%
  as.character() %>% as.numeric()
brier_score(obs_pordm, preds_pordm)

preds_hordm <- predict(hordm_top_final, hordm_top_val, type = "prob") %>% .[[".pred_Yes"]] 
obs_hordm <- hordm_top_val %>% 
  select(HOR) %>% rec(., rec = "No = 0; Yes = 1") %>% .[["HOR_r"]] %>%
  as.character() %>% as.numeric()
brier_score(obs_hordm, preds_hordm)
```

# Calibration curve.
```{r}
# Preparation data for plotting.
plot_porsm <- predict(porsm_top_final, porsm_top_val, type = "prob") %>% 
  bind_cols(porsm_top_val %>% select(POR)) %>% 
  mutate(POR_r = if_else(POR == "Yes", 1, 0)) %>% select(-POR) %>% 
  arrange(.pred_Yes) %>% 
  mutate(pos = row_number()/n(),
         bin = ceiling(pos*10)) %>% 
  group_by(bin) %>% 
  summarise(estimated_prob = mean(.pred_Yes),
            actual_prob = mean(POR_r)) %>% 
  mutate(model = "PORSM (0.091)")

plot_pordm <- predict(pordm_top_final, pordm_top_val, type = "prob") %>% 
  bind_cols(porsm_top_val %>% select(POR)) %>% 
  mutate(POR_r = if_else(POR == "Yes", 1, 0)) %>% select(-POR) %>% 
  arrange(.pred_Yes) %>% 
  mutate(pos = row_number()/n(),
         bin = ceiling(pos*10)) %>% 
  group_by(bin) %>% 
  summarise(estimated_prob = mean(.pred_Yes),
            actual_prob = mean(POR_r)) %>% 
  mutate(model = "PORRM (0.092)")
plot_pordm <- mutate(plot_pordm, model = "PORRM (0.092)")

plot_horsm <- predict(horsm_top_final, horsm_top_val, type = "prob") %>% 
  bind_cols(horsm_top_val %>% select(HOR)) %>% 
  mutate(HOR_r = if_else(HOR == "Yes", 1, 0)) %>% select(-HOR) %>% 
  arrange(.pred_Yes) %>% 
  mutate(pos = row_number()/n(),
         bin = ceiling(pos*10)) %>% 
  group_by(bin) %>% 
  summarise(estimated_prob = mean(.pred_Yes),
            actual_prob = mean(HOR_r)) %>% 
  mutate(model = "HORSM (0.064)")
plot_hordm <- predict(hordm_top_final, horsm_top_val, type = "prob") %>% 
  bind_cols(horsm_top_val %>% select(HOR)) %>% 
  mutate(HOR_r = if_else(HOR == "Yes", 1, 0)) %>% select(-HOR) %>% 
  arrange(.pred_Yes) %>% 
  mutate(pos = row_number()/n(),
         bin = ceiling(pos*10)) %>% 
  group_by(bin) %>% 
  summarise(estimated_prob = mean(.pred_Yes),
            actual_prob = mean(HOR_r)) %>% 
  mutate(model = "HORRM (0.063)")
plot_hordm <- mutate(plot_hordm, model = "HORRM (0.063)")

plot_calibration <- bind_rows(plot_porsm, plot_pordm, plot_horsm, plot_hordm)
# Plot calibration curve.
theme_set(theme_bw(base_size = 10, base_family = "sans", base_rect_size = .5))
p_calibration <- ggplot(data = plot_calibration, 
                        aes(x = estimated_prob, y = actual_prob, fill = model, 
                            color = model, shape = model)) +
  geom_abline(slope = 1, intercept = 0, color = "darkgrey", linetype = "dashed", size = 1) +
  geom_point(size = 3, aes(fill = model)) +
  geom_smooth(method = "gam", alpha = 0.2) +
  scale_fill_manual(values = c("#003652", "#007dbc", "#8b1d00", "#e39e21")) +
  scale_color_manual(values = c("#003652", "#007dbc", "#8b1d00", "#e39e21")) +
  scale_shape_manual(values = c(21, 22, 23, 24, 25)) +
  guides(fill = guide_legend("Submodel (Brier score)"),
         color = guide_legend("Submodel (Brier score)"),
         shape = guide_legend("Submodel (Brier score)")) +
  labs(x = "Predicted probability",
       y = "Observed probability") + 
  coord_cartesian(xlim = c(0, 0.75), ylim = c(0, 0.75)) +
  theme(legend.background = element_blank(),
        legend.title.align = 0.5,
        legend.title = element_text(face = "bold"),
        legend.position = c(1, 0),
        legend.key.height = unit(5, "mm"),
        legend.justification = c(1, 0),
        legend.box.background = element_rect(fill = "White"),
        plot.margin = margin(c(0, 0, 0, 0), "mm"))
ggsave(p_calibration, filename = "Figure S3.png", device = "png", width = 112, height = 80, units = "mm", dpi = 600)
```

# SHAP related plots.
## Prepare datas.
```{r}
porsm_top_data_forshap <- recipe(POR ~., data = porsm_top_data) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  prep() %>% 
  bake(new_data = NULL) %>% 
  select(-POR) %>% 
  as.matrix()
horsm_top_data_forshap <- recipe(HOR ~., data = horsm_top_data) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  prep() %>% 
  bake(new_data = NULL) %>% 
  select(-HOR) %>% 
  as.matrix()
pordm_top_data_forshap <- recipe(POR ~., data = pordm_top_data) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  prep() %>% 
  bake(new_data = NULL) %>% 
  select(-POR) %>% 
  as.matrix()
hordm_top_data_forshap <- recipe(HOR ~., data = hordm_top_data) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  prep() %>% 
  bake(new_data = NULL) %>% 
  select(-HOR) %>% 
  as.matrix()
or_complete <- openxlsx::read.xlsx("E:/BaiduNetdiskWorkspace/or_Rwork/raw_to_amp/original_final_or.xlsx") %>% 
  as_tibble() %>% 
  select(AMH, AFC, FSH, LH, POIorDOR, Age, Weight, Duration, PCOS, PLT, P, DBP, WBC, ALT, RBC, POR, HOR) %>% 
  mutate_if(is.character, as.factor) %>% 
  filter(if_all(everything(), ~ !is.na(.x)))
# Workflow version of PORDM and HORDM.
set.seed(777)
pordm_top_wf_v <- select_best(pordm_top_final_bayes_regs, "roc_auc") %>% 
  finalize_workflow(pordm_top_final_workflow, .) %>% 
  fit(., pordm_top_data)
hordm_top_wf_v <- select_best(hordm_top_final_bayes_regs, "roc_auc") %>% 
  finalize_workflow(hordm_top_final_workflow, .) %>% 
  fit(., hordm_top_data)

pred_or_complete <- or_complete %>% mutate(id = rep(1: 7808)) %>% 
  select(id, everything()) %>% 
  bind_cols(predict(pordm_top_wf_v, or_complete)) %>% 
  rename(pred_POR = .pred_class) %>% 
  bind_cols(predict(pordm_top_wf_v, or_complete, type = "prob")) %>% 
  select(-.pred_No) %>% rename(prob_POR = .pred_Yes) %>% 
  bind_cols(predict(hordm_top_wf_v, or_complete)) %>% 
  rename(pred_HOR = .pred_class) %>% 
  bind_cols(predict(hordm_top_wf_v, or_complete, type = "prob")) %>% 
  select(-.pred_No) %>% rename(prob_HOR = .pred_Yes) 

instance_por <- pred_or_complete %>% 
  filter(id == 2583)
instance_hor <- pred_or_complete %>% 
  filter(id == 1324)
instance_normal <- pred_or_complete %>% 
  filter(id == 1321)
```

## PORDM related plots.
```{r}
# SHAP summary plot.
library(SHAPforxgboost)
pordm_xgb <- 
pordm_shap_prep <- shap.prep(xgb_model = extract_fit_engine(pordm_top_final),
                             X_train = pordm_top_data_forshap)
pordm_shap_trans <- pordm_shap_prep %>% select(variable, rfvalue) %>% 
  mutate(ID = rep(1:12012, 13)) %>% 
  group_by(variable) %>% 
  pivot_wider(names_from = variable, values_from = rfvalue) %>% 
  recipe(~.) %>% 
  step_YeoJohnson(-c(ID, POIorDOR_Yes)) %>% 
  step_filter(AFC<9, FSH>0.9, FSH<2.4, ALT>1.15, LH>0.5, LH<2.5, WBC>1.16, RBC>0.975, RBC<1.175, P<1.4, Duration<2.6) %>% 
  step_range(-c(ID, POIorDOR_Yes), min = 0, max = 1) %>% 
  prep() %>% bake(new_data = NULL) %>%
  pivot_longer(-ID, names_to = "variable", values_to = "rfvalue") %>% 
  left_join(., pordm_shap_prep %>% select(variable, value, stdfvalue, mean_value) %>% 
              mutate(ID = rep(1:12012, 13)), by = c("ID", "variable")) %>% 
  arrange(desc(mean_value)) %>% 
  mutate(ID = rep(1:11762, 13)) %>% select(1:2, value, everything()) %>% 
  mutate(variable = forcats::fct_recode(variable, "POI or DOR" = "POIorDOR_Yes", Progesterone = "P"),
         variable = forcats::fct_reorder(as.factor(variable), .x = mean_value))

theme_set(theme_bw(base_size = 10, base_family = "sans", base_rect_size = 0.3))
p_shapsum_pordm <- pordm_shap_trans %>% 
  group_by(variable) %>% 
  slice_sample(prop = .2) %>% 
  ggplot(aes(x = variable, y = value, color = rfvalue)) +
  ggforce::geom_sina(alpha = 0.5, 
                     size = 0.05,
                     method = "density",
                     maxwidth = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  coord_flip() +
  scale_color_gradientn(colors = c("#007dbc", "#eadbeb", "#8b1d00")) +
  labs(x = "Top features", y = "SHAP value (impact on model output)") +
  guides(color = guide_colorbar(title = "Standardized feature value",
                                title.position = "top", barwidth = 6, barheight = 0.7)) +
  theme(legend.title.align = 0,
        plot.margin = margin(c(0, 6, 0, 0), "mm")) 

# Local explanation
library(DALEXtra)
explainer_pordm <- explain_tidymodels(pordm_top_wf_v,
                                      data = pordm_top_data %>% select(-POR),
                                      y = pordm_top_data %>% 
                                        mutate(POR = if_else(POR == "Yes", 1, 0)) %>% 
                                        .[["POR"]],
                                      label = "PORDM")
set.seed(777)
shap_por <- predict_parts(explainer = explainer_pordm,
                          new_observation = instance_por,
                          type = "shap",
                          B = 100)
shap_normal_pordm <- predict_parts(explainer = explainer_pordm,
                          new_observation = instance_normal,
                          type = "shap",
                          B = 100)

theme_set(theme_bw(base_size = 10, base_family = "sans", base_rect_size = 0.5))
p_shaplocal_pordm1 <- shap_por %>% 
  group_by(variable) %>% 
  mutate(mean_val = mean(contribution)) %>% 
  ungroup() %>% 
  mutate(variable = case_when(variable == "POIorDOR = Yes" ~ "POI/DOR = Yes",
                              TRUE ~ as.character(variable))) %>%
  mutate(variable = forcats::fct_reorder(variable, abs(mean_val))) %>% 
  ggplot(aes(contribution, variable, fill = mean_val >0)) +
  geom_col(data = ~distinct(., variable, mean_val),
           aes(mean_val, variable),
           alpha = .5) +
  geom_boxplot(width = .5, size = 0.2,
               outlier.shape = NA) +
  scale_fill_manual(values = c("#8b1d00", "#007dbc"),
                    labels = c("Drive prediction towards POR/HOR", "Drive prediction towards non-POR/non-HOR")) +
  labs(y = NULL, x = NULL) +
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        legend.direction = "horizontal",
        panel.grid = element_blank(),
        axis.text.y = element_text(angle = 45, hjust = 0.5, vjust = 0.5),
        plot.margin = margin(c(4, 4, 0, 0), "mm"))

theme_set(theme_bw(base_size = 10, base_family = "sans", base_rect_size = 0.5))
p_shaplocal_pordm2 <- shap_normal_pordm %>% 
  group_by(variable) %>% 
  mutate(mean_val = mean(contribution)) %>% 
  ungroup() %>% 
  mutate(variable = case_when(variable == "POIorDOR = No" ~ "POI/DOR = No",
                              TRUE ~ as.character(variable))) %>%
  mutate(variable = forcats::fct_reorder(variable, abs(mean_val))) %>% 
  ggplot(aes(contribution, variable, fill = mean_val >0)) +
  geom_col(data = ~distinct(., variable, mean_val),
           aes(mean_val, variable),
           alpha = .5) +
  geom_boxplot(width = .5, size = 0.2,
               outlier.shape = NA) +
  scale_fill_manual(values = c("#8b1d00", "#007dbc")) +
  labs(y = NULL, x = "SHAP value (impact on POR prediction)") +
  theme(legend.position = "none",
        panel.grid = element_blank(),
        axis.text.y = element_text(angle = 45, hjust = 0.5, vjust = 0.5),
        plot.margin = margin(c(4, 4, 0, 0), "mm"))
```

## HORDM related plots.
```{r}
# SHAP summary plot.
library(SHAPforxgboost)
hordm_shap_prep <- shap.prep(xgb_model = extract_fit_engine(hordm_top_final),
                             X_train = hordm_top_data_forshap)
hordm_shap_trans <- hordm_shap_prep %>% select(variable, rfvalue) %>% 
  mutate(ID = rep(1:12012, 10)) %>% 
  group_by(variable) %>% 
  pivot_wider(names_from = variable, values_from = rfvalue) %>% 
  recipe(~.) %>% 
  step_YeoJohnson(-c(ID, POIorDOR_Yes, PCOS_Yes)) %>% 
  step_filter(AFC<9, FSH>0.9, FSH<2.4, LH>0.5, LH<2.5, Duration<2.6, PLT < 23, PLT > 12.5) %>% 
  step_range(-c(ID, POIorDOR_Yes, PCOS_Yes), min = 0, max = 1) %>% 
  prep() %>% bake(new_data = NULL) %>%  
  pivot_longer(-ID, names_to = "variable", values_to = "rfvalue") %>% zz
  left_join(., hordm_shap_prep %>% select(variable, value, stdfvalue, mean_value) %>% 
              mutate(ID = rep(1:12012, 10)), by = c("ID", "variable")) %>% 
  arrange(desc(mean_value)) %>% 
  mutate(ID = rep(1:11743, 10)) %>% select(1:2, value, everything()) %>% 
  mutate(variable = forcats::fct_recode(variable, "POI or DOR" = "POIorDOR_Yes",
                                        "PCOS" = "PCOS_Yes", Platelet = "PLT"),
         variable = forcats::fct_reorder(as.factor(variable), .x = mean_value)) 

theme_set(theme_bw(base_size = 10, base_family = "sans", base_rect_size = 0.3))
p_shapsum_hordm <- hordm_shap_trans %>% 
  group_by(variable) %>% 
  slice_sample(prop = .2) %>% 
  ggplot(aes(x = variable, y = value, color = rfvalue)) +
  ggforce::geom_sina(alpha = 0.5, 
                     size = 0.05, 
                     method = "counts",
                     maxwidth = 0.75,
                     show.legend = FALSE) + 
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  coord_flip() +
  scale_color_gradientn(colors = c("#007dbc", "#eadbeb", "#8b1d00")) +
  labs(x = element_blank(), y = "SHAP value (impact on model output)") +
  guides(color = guide_colorbar(title = "Standardized feature value",
                                title.position = "top", barwidth = 6, barheight = 0.7)) +
  theme(plot.margin = margin(c(0, 0, 0, 4), "mm")) 

# Local explanation
library(DALEXtra)
explainer_hordm <- explain_tidymodels(hordm_top_wf_v,
                                      data = hordm_top_data %>% select(-HOR),
                                      y = hordm_top_data %>% 
                                        mutate(HOR = if_else(HOR == "Yes", 1, 0)) %>% 
                                        .[["HOR"]],
                                      label = "HORDM")
set.seed(777)
shap_hor <- predict_parts(explainer = explainer_hordm,
                          new_observation = instance_hor,
                          type = "shap",
                          B = 100)
shap_normal_hordm <- predict_parts(explainer = explainer_hordm,
                          new_observation = instance_normal,
                          type = "shap",
                          B = 100)

theme_set(theme_bw(base_size = 10, base_family = "sans", base_rect_size = 0.5))
p_shaplocal_hordm1 <- shap_hor %>% 
  group_by(variable) %>% 
  mutate(mean_val = mean(contribution)) %>% 
  ungroup() %>% 
  mutate(variable = case_when(variable == "POIorDOR = No" ~ "POI/DOR = No",
                              TRUE ~ as.character(variable))) %>% 
  mutate(variable = forcats::fct_reorder(variable, abs(mean_val))) %>%
  ggplot(aes(contribution, variable, fill = mean_val >0)) +
  geom_col(data = ~distinct(., variable, mean_val),
           aes(mean_val, variable),
           alpha = .5) +
  geom_boxplot(width = .5, size = 0.2,
               outlier.shape = NA) +
  scale_fill_manual(values = c("#8b1d00", "#007dbc")) +
  labs(y = NULL, x = NULL) +
  theme(legend.position = "none",
        panel.grid = element_blank(),
        axis.text.y = element_text(angle = 45, hjust = 0.5, vjust = 0.5),
        plot.margin = margin(c(4, 0, 0, 4), "mm"))

theme_set(theme_bw(base_size = 10, base_family = "sans", base_rect_size = 0.5))
p_shaplocal_hordm2 <- shap_normal_hordm %>% 
  group_by(variable) %>% 
  mutate(mean_val = mean(contribution)) %>% 
  ungroup() %>% 
  mutate(variable = case_when(variable == "POIorDOR = No" ~ "POI/DOR = No",
                              TRUE ~ as.character(variable))) %>% 
  mutate(variable = forcats::fct_reorder(variable, abs(mean_val))) %>% 
  ggplot(aes(contribution, variable, fill = mean_val >0)) +
  geom_col(data = ~distinct(., variable, mean_val),
           aes(mean_val, variable),
           alpha = .5) +
  geom_boxplot(width = .5, size = 0.2,
               outlier.shape = NA) +
  scale_fill_manual(values = c("#8b1d00", "#007dbc")) +
  labs(y = NULL, x = "SHAP value (impact on HOR prediction)") +
  theme(legend.position = "none",
        panel.grid = element_blank(),
        axis.text.y = element_text(angle = 45, hjust = 0.5, vjust = 0.5),
        plot.margin = margin(c(4, 0, 0, 4), "mm"))
```

## Combinations of plots
```{r}
# Combination of summary plots.
library(ggpubr)
ggarrange(p_shapsum_pordm, p_shapsum_hordm, nrow = 1, widths = c(1.05, 1),
          labels = c("B", "C"),
          font.label = list(size = 12, face = "bold", family = "sans"),
          common.legend = TRUE, legend = "bottom",
          hjust = c(0, 0), vjust = c(1, 1)) %>% 
ggsave(filename = "Figure 4BC.pdf", device = cairo_pdf, width = 172, height = 100, units = "mm")

# Combination of local explanation plots
p_local <- ggarrange(p_shaplocal_pordm1, p_shaplocal_hordm1, p_shaplocal_pordm2,p_shaplocal_hordm2, 
                     ncol = 2, nrow = 2, labels = c("A", "C", "B", "D"),
                     font.label = list(size = 12, face = "bold", family = "sans"),
                     hjust = c(0, -0.5, 0, -0.5), vjust = c(1, 1, 1, 1),
                     legend = "bottom", common.legend = TRUE)
ggsave(p_local, filename = "Figure S4.png", device = "png", width = 172, height = 215, units = "mm", dpi = 600)
```