ML system (part 1).Rmd

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

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

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

# Construction of Strategy models.
## Tuning parameters for **PORSM**
```{r}
library(tidymodels)
set.seed(777)
# Preparation of data and specifications
porsm_top_data <- porsm_data %>% select(AMH, AFC, POIorDOR, FSH, Age, P, Weight, DBP, WBC, ALT, RBC, Duration, LH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_top_train <- porsm_train %>% select(AMH, AFC, POIorDOR, FSH, Age, P, Weight, DBP, WBC, ALT, RBC, Duration, LH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_top_test <- porsm_test %>% select(AMH, AFC, POIorDOR, FSH, Age, P, Weight, DBP, WBC, ALT, RBC, Duration, LH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_top_rec <- recipe(POR ~ ., data = porsm_top_train) %>% 
  step_dummy(all_nominal_predictors())
porsm_top_mod <- 
  boost_tree(trees = tune(), 
             min_n = tune(), 
             tree_depth = tune(),
             loss_reduction = tune(),
             learn_rate = tune(), 
             sample_size = tune(), 
             mtry = tune()) %>% 
  set_engine("xgboost") %>% 
  set_mode("classification")
porsm_top_workflow <- workflow() %>% 
  add_model(porsm_top_mod) %>% 
  add_recipe(porsm_top_rec)
porsm_top_bayes_params <- parameters(porsm_top_workflow) %>% 
  update(mtry = finalize(mtry(), porsm_top_train),
         trees = trees(range = c(100L, 1000L)))
# Create cross-validation resamples.
set.seed(777)
porsm_top_valset <- vfold_cv(porsm_top_train, v = 5, strata = POR)
# Start tuning.
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
porsm_top_bayes_regs <- tune_bayes(
  porsm_top_workflow,
  resamples = porsm_top_valset,
  param_info = porsm_top_bayes_params,
  iter = 100, 
  metrics = metric_set(roc_auc),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)
```
## Tuning parameters for **HORSM**
```{r}
# Preparation of data and specifications
set.seed(777)
horsm_top_data <- horsm_data %>% select(AMH, AFC, FSH, Age, LH, POIorDOR, PCOS, PLT, Weight, Duration, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_top_train <- horsm_train %>% select(AMH, AFC, FSH, Age, LH, POIorDOR, PCOS, PLT, Weight, Duration, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_top_test <- horsm_test %>% select(AMH, AFC, FSH, Age, LH, POIorDOR, PCOS, PLT, Weight, Duration, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_top_rec <- recipe(HOR ~ ., data = horsm_top_train) %>% 
  step_dummy(all_nominal_predictors())
horsm_top_mod <- 
  boost_tree(trees = tune(), 
             min_n = tune(), 
             tree_depth = tune(),
             loss_reduction = tune(),
             learn_rate = tune(), 
             sample_size = tune(), 
             mtry = tune()) %>% 
  set_engine("xgboost") %>% 
  set_mode("classification")
horsm_top_workflow <- workflow() %>% 
  add_model(horsm_top_mod) %>% 
  add_recipe(horsm_top_rec)
horsm_top_bayes_params <- parameters(horsm_top_workflow) %>% 
  update(mtry = finalize(mtry(), horsm_top_train),
         trees = trees(range = c(100L, 1000L)))
# Create cross-validation resamples.
set.seed(777)
horsm_top_valset <- vfold_cv(horsm_top_train, v = 5, strata = HOR)
# Start tuning.
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
horsm_top_bayes_regs <- tune_bayes(
  horsm_top_workflow,
  resamples = horsm_top_valset,
  param_info = horsm_top_bayes_params,
  iter = 100, 
  metrics = metric_set(roc_auc),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)
```
## Developing strategy models on train data.
```{r}
set.seed(777)
# Strategy models on train data.
porsm_top <- select_best(porsm_top_bayes_regs, "roc_auc") %>% 
  finalize_workflow(porsm_top_workflow, .) %>% 
  fit(., porsm_top_train)
horsm_top <- select_best(horsm_top_bayes_regs, "roc_auc") %>% 
  finalize_workflow(horsm_top_workflow, .) %>% 
  fit(., horsm_top_train)
```
## Calculate AUC and Brier score.
```{r}
set.seed(777)
# AUC and 95%CI.
roc_porsm_top <- porsm_top_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_top, porsm_top_test, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) 
roc_porsm_top %>% pROC::ci.auc(method = "bootstrap")
set.seed(777)
roc_horsm_top <- horsm_top_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_top, horsm_top_test, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_top %>% pROC::ci.auc(method = "bootstrap")

# Brier score.
brier_score <- function(preds, obs) {
  mean((obs - preds)^2)
}
preds_porsm_top <- predict(porsm_top, porsm_top_test, type = "prob") %>% .[[".pred_Yes"]] 
obs_porsm_top <- porsm_top_test %>% select(POR) %>% sjmisc::rec(., rec = "No = 0; Yes = 1") %>% .[["POR_r"]] %>%
  as.character()%>% as.numeric()
brier_score(obs_porsm_top, preds_porsm_top)

preds_horsm_top <- predict(horsm_top, horsm_top_test, type = "prob") %>% .[[".pred_Yes"]] 
obs_horsm_top <- horsm_top_test %>% select(HOR) %>% sjmisc::rec(., rec = "No = 0; Yes = 1") %>% .[["HOR_r"]] %>%
  as.character()%>% as.numeric()
brier_score(obs_horsm_top, preds_horsm_top)
```
## Developing strategy models on the whole derivation data.
### Tuning parameters.
```{r}
# Tuning for PORSM
porsm_top_final_rec <- recipe(POR ~ ., data = porsm_top_data) %>% 
  step_dummy(all_nominal_predictors())
porsm_top_final_workflow <- workflow() %>% 
  add_model(porsm_top_mod) %>% 
  add_recipe(porsm_top_final_rec)
porsm_top_final_bayes_params <- parameters(porsm_top_final_workflow) %>% 
  update(mtry = finalize(mtry(), porsm_top_data),
         trees = trees(range = c(100L, 1000L)))

set.seed(777)
porsm_top_final_valset <- vfold_cv(porsm_top_data, v = 5, strata = POR)
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
porsm_top_final_bayes_regs <- tune_bayes(
  porsm_top_final_workflow,
  resamples = porsm_top_final_valset,
  param_info = porsm_top_final_bayes_params,
  iter = 100, 
  metrics = metric_set(roc_auc),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)

# Tuning for HORSM
horsm_top_final_rec <- recipe(HOR ~ ., data = horsm_top_data) %>% 
  step_dummy(all_nominal_predictors())
horsm_top_final_workflow <- workflow() %>% 
  add_model(horsm_top_mod) %>% 
  add_recipe(horsm_top_final_rec)
horsm_top_final_bayes_params <- parameters(horsm_top_final_workflow) %>% 
  update(mtry = finalize(mtry(), horsm_top_data),
         trees = trees(range = c(100L, 1000L)))

set.seed(777)
horsm_top_final_valset <- vfold_cv(horsm_top_data, v = 5, strata = HOR)
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
horsm_top_final_bayes_regs <- tune_bayes(
  horsm_top_final_workflow,
  resamples = horsm_top_final_valset,
  param_info = horsm_top_final_bayes_params,
  iter = 100, 
  metrics = metric_set(roc_auc),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)
```

### Construction final strategy submodels.
```{r}
set.seed(777)
porsm_top_final <- select_best(porsm_top_final_bayes_regs, "roc_auc") %>% 
  finalize_model(porsm_top_mod, .) %>% 
  fit(., POR ~., data = porsm_top_final_rec %>% prep() %>% bake(new_data = NULL))
horsm_top_final <- select_best(horsm_top_final_bayes_regs, "roc_auc") %>% 
  finalize_model(horsm_top_mod, .) %>% 
  fit(., HOR ~., data = horsm_top_final_rec %>% prep() %>% bake(new_data = NULL))

## Do not run! (This is just for fun!)
set.seed(777)
porsm_top_data %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_top_final, porsm_top_final_rec %>% prep() %>% bake(new_data = NULL), type = "prob")) %>%
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) 
set.seed(777)
horsm_top_data %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_top_final, horsm_top_final_rec %>% prep() %>% bake(new_data = NULL), type = "prob")) %>%
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) 
```

# Construction of Diagnostic models.
## Tuning parameters for **PORDM**
```{r}
library(tidymodels)
set.seed(777)
# Preparation of data and specifications
pordm_top_data <- porsm_top_data %>% select(-c(Protocol, Initial.FSH, Recombinant, Use.LH))
pordm_top_train <- porsm_top_train %>% select(-c(Protocol, Initial.FSH, Recombinant, Use.LH))
pordm_top_test <- porsm_top_test %>% select(-c(Protocol, Initial.FSH, Recombinant, Use.LH))
pordm_top_rec <- recipe(POR ~ ., data = pordm_top_train) %>% 
  step_dummy(all_nominal_predictors())
pordm_top_mod <- 
  boost_tree(trees = tune(), 
             min_n = tune(), 
             tree_depth = tune(),
             loss_reduction = tune(),
             learn_rate = tune(), 
             sample_size = tune(), 
             mtry = tune()) %>% 
  set_engine("xgboost") %>% 
  set_mode("classification")
pordm_top_workflow <- workflow() %>% 
  add_model(pordm_top_mod) %>% 
  add_recipe(pordm_top_rec)
pordm_top_bayes_params <- parameters(pordm_top_workflow) %>% 
  update(mtry = finalize(mtry(), pordm_top_train),
         trees = trees(range = c(100L, 1000L)))
# Create cross-validation resamples.
set.seed(777)
pordm_top_valset <- vfold_cv(pordm_top_train, v = 5, strata = POR)
# Start tuning.
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
pordm_top_bayes_regs <- tune_bayes(
  pordm_top_workflow,
  resamples = pordm_top_valset,
  param_info = pordm_top_bayes_params,
  iter = 100, 
  metrics = metric_set(roc_auc),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)
```
## Tuning parameters for **HORDM**
```{r}
# Preparation of data and specifications
set.seed(777)
hordm_top_data <- horsm_top_data %>% select(-c(Protocol, Initial.FSH, Recombinant, Use.LH))
hordm_top_train <- horsm_top_train %>% select(-c(Protocol, Initial.FSH, Recombinant, Use.LH))
hordm_top_test <- horsm_top_test %>% select(-c(Protocol, Initial.FSH, Recombinant, Use.LH))
hordm_top_rec <- recipe(HOR ~ ., data = hordm_top_train) %>% 
  step_dummy(all_nominal_predictors())
hordm_top_mod <- 
  boost_tree(trees = tune(), 
             min_n = tune(), 
             tree_depth = tune(),
             loss_reduction = tune(),
             learn_rate = tune(), 
             sample_size = tune(), 
             mtry = tune()) %>% 
  set_engine("xgboost") %>% 
  set_mode("classification")
hordm_top_workflow <- workflow() %>% 
  add_model(hordm_top_mod) %>% 
  add_recipe(hordm_top_rec)
hordm_top_bayes_params <- parameters(hordm_top_workflow) %>% 
  update(mtry = finalize(mtry(), hordm_top_train),
         trees = trees(range = c(100L, 1000L)))
# Create cross-validation resamples.
set.seed(777)
hordm_top_valset <- vfold_cv(hordm_top_train, v = 5, strata = HOR)
# Start tuning.
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
hordm_top_bayes_regs <- tune_bayes(
  hordm_top_workflow,
  resamples = hordm_top_valset,
  param_info = hordm_top_bayes_params,
  iter = 100, 
  metrics = metric_set(roc_auc),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)
```
## Developing dianosis models on train data.
```{r}
set.seed(777)
# Strategy models on train data.
pordm_top <- select_best(pordm_top_bayes_regs, "roc_auc") %>% 
  finalize_workflow(pordm_top_workflow, .) %>% 
  fit(., pordm_top_train)
hordm_top <- select_best(hordm_top_bayes_regs, "roc_auc") %>% 
  finalize_workflow(hordm_top_workflow, .) %>% 
  fit(., hordm_top_train)
```
## Calculate AUC and Brier score.
```{r}
# AUC and 95%CI.
set.seed(777)
roc_pordm_top <- pordm_top_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_top, pordm_top_test, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE)
roc_pordm_top %>% pROC::ci.auc(method = "bootstrap")
set.seed(777)
roc_hordm_top <- hordm_top_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top, hordm_top_test, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) 
roc_hordm_top %>% pROC::ci.auc(method = "bootstrap")

# Brier score.
brier_score <- function(preds, obs) {
  mean((obs - preds)^2)
}
preds_pordm_top <- predict(pordm_top, pordm_top_test, type = "prob") %>% .[[".pred_Yes"]] 
obs_pordm_top <- pordm_top_test %>% select(POR) %>% sjmisc::rec(., rec = "No = 0; Yes = 1") %>% .[["POR_r"]] %>%
  as.character()%>% as.numeric()
brier_score(obs_pordm_top, preds_pordm_top)

preds_hordm_top <- predict(hordm_top, hordm_top_test, type = "prob") %>% .[[".pred_Yes"]] 
obs_hordm_top <- hordm_top_test %>% select(HOR) %>% sjmisc::rec(., rec = "No = 0; Yes = 1") %>% .[["HOR_r"]] %>%
  as.character()%>% as.numeric()
brier_score(obs_hordm_top, preds_hordm_top)
```

## Developing diagnosis models on the whole derivation data.
### Tuning parameters.
```{r}
# Tuning for PORDM.
library(tidymodels)
set.seed(777)
pordm_top_final_rec <- recipe(POR ~ ., data = pordm_top_data) %>% 
  step_dummy(all_nominal_predictors())
pordm_top_final_workflow <- workflow() %>% 
  add_model(pordm_top_mod) %>% 
  add_recipe(pordm_top_final_rec)
pordm_top_final_bayes_params <- parameters(pordm_top_final_workflow) %>% 
  update(mtry = finalize(mtry(), pordm_top_data),
         trees = trees(range = c(100L, 1000L)))

set.seed(777)
pordm_top_final_valset <- vfold_cv(pordm_top_data, v = 5, strata = POR)
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
pordm_top_final_bayes_regs <- tune_bayes(
  pordm_top_final_workflow,
  resamples = pordm_top_final_valset,
  param_info = pordm_top_final_bayes_params,
  iter = 100, 
  metrics = metric_set(roc_auc),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)

# Tuning for HORDM.
set.seed(777)
hordm_top_final_rec <- recipe(HOR ~ ., data = hordm_top_data) %>% 
  step_dummy(all_nominal_predictors())
hordm_top_final_workflow <- workflow() %>% 
  add_model(hordm_top_mod) %>% 
  add_recipe(hordm_top_final_rec)
hordm_top_final_bayes_params <- parameters(hordm_top_final_workflow) %>% 
  update(mtry = finalize(mtry(), hordm_top_data),
         trees = trees(range = c(100L, 1000L)))

set.seed(777)
hordm_top_final_valset <- vfold_cv(hordm_top_data, v = 5, strata = HOR)
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
hordm_top_final_bayes_regs <- tune_bayes(
  hordm_top_final_workflow,
  resamples = hordm_top_final_valset,
  param_info = hordm_top_final_bayes_params,
  iter = 100, 
  metrics = metric_set(roc_auc),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)
```

### Construction final diagnostic submodels.
```{r}
set.seed(777)
pordm_top_final <- select_best(pordm_top_final_bayes_regs, "roc_auc") %>% 
  finalize_model(pordm_top_mod, .) %>% 
  fit(., POR ~., data = pordm_top_final_rec %>% prep() %>% bake(new_data = NULL))
hordm_top_final <- select_best(hordm_top_final_bayes_regs, "roc_auc") %>% 
  finalize_model(hordm_top_mod, .) %>% 
  fit(., HOR ~., data = hordm_top_final_rec %>% prep() %>% bake(new_data = NULL))
```

# Superiority analysis: Construction strategy models on conventional combinations of widely used markers and on complete case data.
## Candidate PORSMs with conventional combinations of markers.
```{r}
library(pROC)
library(sjmisc)
set.seed(777)
# Age only.
porsm_age_train <- porsm_train %>% select(Age, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_age_test <- porsm_test %>% select(Age, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_age_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_age_wf <- workflow() %>% add_model(porsm_age_mod) %>% add_formula(POR ~.)
porsm_age <- fit(porsm_age_wf, data = porsm_age_train)
roc_porsm_age <- porsm_age_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_age, porsm_age_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_porsm_age %>% ci.auc(method = "bootstrap")
# FSH only
porsm_fsh_train <- porsm_train %>% select(FSH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_fsh_test <- porsm_test %>% select(FSH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_fsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_fsh_wf <- workflow() %>% add_model(porsm_fsh_mod) %>% add_formula(POR ~.)
porsm_fsh <- fit(porsm_fsh_wf, data = porsm_fsh_train)
roc_porsm_fsh <- porsm_fsh_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_fsh, porsm_fsh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_porsm_fsh %>% ci.auc(method = "bootstrap")
# AFC only.
porsm_afc_train <- porsm_train %>% select(AFC, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_afc_test <- porsm_test %>% select(AFC, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_afc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_afc_wf <- workflow() %>% add_model(porsm_afc_mod) %>% add_formula(POR ~.)
porsm_afc <- fit(porsm_afc_wf, data = porsm_afc_train)
roc_porsm_afc <- porsm_afc_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_afc, porsm_afc_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)
roc_porsm_afc %>% ci.auc(method = "bootstrap")
# AMH only.
porsm_amh_train <- porsm_train %>% select(AMH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_amh_test <- porsm_test %>% select(AMH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_amh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_amh_wf <- workflow() %>% add_model(porsm_amh_mod) %>% add_formula(POR ~.)
porsm_amh <- fit(porsm_amh_wf, data = porsm_amh_train)
roc_porsm_amh <- porsm_amh_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_amh, porsm_amh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_porsm_amh %>% ci.auc(method = "bootstrap")
# Age + FSH
porsm_agefsh_train <- porsm_train %>% select(Age, FSH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_agefsh_test <- porsm_test %>% select(Age, FSH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_agefsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_agefsh_wf <- workflow() %>% add_model(porsm_agefsh_mod) %>% add_formula(POR ~.)
porsm_agefsh <- fit(porsm_agefsh_wf, data = porsm_agefsh_train)
roc_porsm_agefsh <- porsm_agefsh_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_agefsh, porsm_agefsh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)
roc_porsm_agefsh %>% ci.auc(method = "bootstrap")
# Age+ AFC
porsm_ageafc_train <- porsm_train %>% select(Age, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageafc_test <- porsm_test %>% select(Age, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageafc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_ageafc_wf <- workflow() %>% add_model(porsm_ageafc_mod) %>% add_formula(POR ~.)
porsm_ageafc <- fit(porsm_ageafc_wf, data = porsm_ageafc_train)
roc_porsm_ageafc <- porsm_ageafc_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_ageafc, porsm_ageafc_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)  
roc_porsm_ageafc %>% ci.auc(method = "bootstrap")
# Age + AMH
porsm_ageamh_train <- porsm_train %>% select(Age, AMH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageamh_test <- porsm_test %>% select(Age, AMH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageamh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_ageamh_wf <- workflow() %>% add_model(porsm_ageamh_mod) %>% add_formula(POR ~.)
porsm_ageamh <- fit(porsm_ageamh_wf, data = porsm_ageamh_train)
roc_porsm_ageamh <- porsm_ageamh_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_ageamh, porsm_ageamh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)
roc_porsm_ageamh %>% ci.auc(method = "bootstrap")
# Age +AMH +FSH
porsm_ageamhfsh_train <- porsm_train %>% select(Age, AMH, FSH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageamhfsh_test <- porsm_test %>% select(Age, AMH, FSH, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageamhfsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_ageamhfsh_wf <- workflow() %>% add_model(porsm_ageamhfsh_mod) %>% add_formula(POR ~.)
porsm_ageamhfsh <- fit(porsm_ageamhfsh_wf, data = porsm_ageamhfsh_train)
roc_porsm_ageamhfsh <- porsm_ageamhfsh_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_ageamhfsh, porsm_ageamhfsh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_porsm_ageamhfsh %>% ci.auc(method = "bootstrap")
# Age + AMH + AFC
porsm_ageamhafc_train <- porsm_train %>% select(Age, AMH, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageamhafc_test <- porsm_test %>% select(Age, AMH, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageamhafc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_ageamhafc_wf <- workflow() %>% add_model(porsm_ageamhafc_mod) %>% add_formula(POR ~.)
porsm_ageamhafc <- fit(porsm_ageamhafc_wf, data = porsm_ageamhafc_train)
roc_porsm_ageamhafc <- porsm_ageamhafc_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_ageamhafc, porsm_ageamhafc_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_porsm_ageamhafc %>% ci.auc(method = "bootstrap")
# Age + AMH + FSH +AFC
porsm_ageamhafcfsh_train <- porsm_train %>% select(Age, AMH, FSH, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageamhafcfsh_test <- porsm_test %>% select(Age, AMH, FSH, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, POR)
porsm_ageamhafcfsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
porsm_ageamhafcfsh_wf <- workflow() %>% add_model(porsm_ageamhafcfsh_mod) %>% add_formula(POR ~.)
porsm_ageamhafcfsh <- fit(porsm_ageamhafcfsh_wf, data = porsm_ageamhafcfsh_train)
roc_porsm_ageamhafcfsh <- porsm_ageamhafcfsh_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_ageamhafcfsh, porsm_ageamhafcfsh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)  
roc_porsm_ageamhafcfsh %>% ci.auc(method = "bootstrap")
```

## Candidate HORSMs with conventional combinations of markers.
```{r}
set.seed(777)
# Age only.
horsm_age_train <- horsm_train %>% select(Age, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_age_test <- horsm_test %>% select(Age, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_age_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_age_wf <- workflow() %>% add_model(horsm_age_mod) %>% add_formula(HOR ~.)
horsm_age <- fit(horsm_age_wf, data = horsm_age_train)
roc_horsm_age <- horsm_age_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_age, horsm_age_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_age %>% ci.auc(method = "bootstrap")
# FSH only.
horsm_fsh_train <- horsm_train %>% select(FSH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_fsh_test <- horsm_test %>% select(FSH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_fsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_fsh_wf <- workflow() %>% add_model(horsm_fsh_mod) %>% add_formula(HOR ~.)
horsm_fsh <- fit(horsm_fsh_wf, data = horsm_fsh_train)
roc_horsm_fsh <- horsm_fsh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_fsh, horsm_fsh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_fsh %>% ci.auc(method = "bootstrap")
# AFC only.
horsm_afc_train <- horsm_train %>% select(AFC, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_afc_test <- horsm_test %>% select(AFC, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_afc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_afc_wf <- workflow() %>% add_model(horsm_afc_mod) %>% add_formula(HOR ~.)
horsm_afc <- fit(horsm_afc_wf, data = horsm_afc_train)
roc_horsm_afc <- horsm_afc_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_afc, horsm_afc_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_afc %>% ci.auc(method = "bootstrap")
# AMH only.
horsm_amh_train <- horsm_train %>% select(AMH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_amh_test <- horsm_test %>% select(AMH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_amh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_amh_wf <- workflow() %>% add_model(horsm_amh_mod) %>% add_formula(HOR ~.)
horsm_amh <- fit(horsm_amh_wf, data = horsm_amh_train)
roc_horsm_amh <- horsm_amh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_amh, horsm_amh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_amh %>% ci.auc(method = "bootstrap")
# Age + FSH
horsm_agefsh_train <- horsm_train %>% select(Age, FSH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_agefsh_test <- horsm_test %>% select(Age, FSH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_agefsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_agefsh_wf <- workflow() %>% add_model(horsm_agefsh_mod) %>% add_formula(HOR ~.)
horsm_agefsh <- fit(horsm_agefsh_wf, data = horsm_agefsh_train)
roc_horsm_agefsh <- horsm_agefsh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_agefsh, horsm_agefsh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_agefsh %>% ci.auc(method = "bootstrap")
# Age+ AFC
horsm_ageafc_train <- horsm_train %>% select(Age, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageafc_test <- horsm_test %>% select(Age, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageafc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_ageafc_wf <- workflow() %>% add_model(horsm_ageafc_mod) %>% add_formula(HOR ~.)
horsm_ageafc <- fit(horsm_ageafc_wf, data = horsm_ageafc_train)
roc_horsm_ageafc <- horsm_ageafc_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_ageafc, horsm_ageafc_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_ageafc %>% ci.auc(method = "bootstrap")
# Age + AMH
horsm_ageamh_train <- horsm_train %>% select(Age, AMH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageamh_test <- horsm_test %>% select(Age, AMH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageamh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_ageamh_wf <- workflow() %>% add_model(horsm_ageamh_mod) %>% add_formula(HOR ~.)
horsm_ageamh <- fit(horsm_ageamh_wf, data = horsm_ageamh_train)
roc_horsm_ageamh <- horsm_ageamh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_ageamh, horsm_ageamh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_ageamh %>% ci.auc(method = "bootstrap")
# Age +AMH +FSH
horsm_ageamhfsh_train <- horsm_train %>% select(Age, AMH, FSH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageamhfsh_test <- horsm_test %>% select(Age, AMH, FSH, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageamhfsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_ageamhfsh_wf <- workflow() %>% add_model(horsm_ageamhfsh_mod) %>% add_formula(HOR ~.)
horsm_ageamhfsh <- fit(horsm_ageamhfsh_wf, data = horsm_ageamhfsh_train)
roc_horsm_ageamhfsh <- horsm_ageamhfsh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_ageamhfsh, horsm_ageamhfsh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_ageamhfsh %>% ci.auc(method = "bootstrap")
# Age + AMH + AFC
horsm_ageamhafc_train <- horsm_train %>% select(Age, AMH, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageamhafc_test <- horsm_test %>% select(Age, AMH, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageamhafc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_ageamhafc_wf <- workflow() %>% add_model(horsm_ageamhafc_mod) %>% add_formula(HOR ~.)
horsm_ageamhafc <- fit(horsm_ageamhafc_wf, data = horsm_ageamhafc_train)
roc_horsm_ageamhafc <- horsm_ageamhafc_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_ageamhafc, horsm_ageamhafc_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_ageamhafc %>% ci.auc(method = "bootstrap")
# Age + AMH + FSH +AFC
horsm_ageamhafcfsh_train <- horsm_train %>% select(Age, AMH, FSH, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageamhafcfsh_test <- horsm_test %>% select(Age, AMH, FSH, AFC, Protocol, Initial.FSH, Recombinant, Use.LH, HOR)
horsm_ageamhafcfsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
horsm_ageamhafcfsh_wf <- workflow() %>% add_model(horsm_ageamhafcfsh_mod) %>% add_formula(HOR ~.)
horsm_ageamhafcfsh <- fit(horsm_ageamhafcfsh_wf, data = horsm_ageamhafcfsh_train)
roc_horsm_ageamhafcfsh <- horsm_ageamhafcfsh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_ageamhafcfsh, horsm_ageamhafcfsh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_ageamhafcfsh %>% ci.auc(method = "bootstrap")
```

## Candidate PORDMs with conventional combinations of markers.
```{r}
library(pROC)
library(sjmisc)
set.seed(777)
# Age only.
pordm_age_train <- porsm_train %>% select(Age, POR)
pordm_age_test <- porsm_test %>% select(Age, POR)
pordm_age_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_age_wf <- workflow() %>% add_model(pordm_age_mod) %>% add_formula(POR ~.)
pordm_age <- fit(pordm_age_wf, data = pordm_age_train)
roc_pordm_age <- pordm_age_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_age, pordm_age_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_pordm_age %>% ci.auc(method = "bootstrap")
# FSH only
pordm_fsh_train <- porsm_train %>% select(FSH, POR)
pordm_fsh_test <- porsm_test %>% select(FSH, POR)
pordm_fsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_fsh_wf <- workflow() %>% add_model(pordm_fsh_mod) %>% add_formula(POR ~.)
pordm_fsh <- fit(pordm_fsh_wf, data = pordm_fsh_train)
roc_pordm_fsh <- pordm_fsh_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_fsh, pordm_fsh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_pordm_fsh %>% ci.auc(method = "bootstrap")
# AFC only.
pordm_afc_train <- porsm_train %>% select(AFC, POR)
pordm_afc_test <- porsm_test %>% select(AFC, POR)
pordm_afc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_afc_wf <- workflow() %>% add_model(pordm_afc_mod) %>% add_formula(POR ~.)
pordm_afc <- fit(pordm_afc_wf, data = pordm_afc_train)
roc_pordm_afc <- pordm_afc_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_afc, pordm_afc_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)
roc_pordm_afc %>% ci.auc(method = "bootstrap")
# AMH only.
pordm_amh_train <- porsm_train %>% select(AMH, POR)
pordm_amh_test <- porsm_test %>% select(AMH, POR)
pordm_amh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_amh_wf <- workflow() %>% add_model(pordm_amh_mod) %>% add_formula(POR ~.)
pordm_amh <- fit(pordm_amh_wf, data = pordm_amh_train)
roc_pordm_amh <- pordm_amh_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_amh, pordm_amh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_pordm_amh %>% ci.auc(method = "bootstrap")
# Age + FSH
pordm_agefsh_train <- porsm_train %>% select(Age, FSH, POR)
pordm_agefsh_test <- porsm_test %>% select(Age, FSH, POR)
pordm_agefsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_agefsh_wf <- workflow() %>% add_model(pordm_agefsh_mod) %>% add_formula(POR ~.)
pordm_agefsh <- fit(pordm_agefsh_wf, data = pordm_agefsh_train)
roc_pordm_agefsh <- pordm_agefsh_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_agefsh, pordm_agefsh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)
roc_pordm_agefsh %>% ci.auc(method = "bootstrap")
# Age+ AFC
pordm_ageafc_train <- porsm_train %>% select(Age, AFC, POR)
pordm_ageafc_test <- porsm_test %>% select(Age, AFC, POR)
pordm_ageafc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_ageafc_wf <- workflow() %>% add_model(pordm_ageafc_mod) %>% add_formula(POR ~.)
pordm_ageafc <- fit(pordm_ageafc_wf, data = pordm_ageafc_train)
roc_pordm_ageafc <- pordm_ageafc_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_ageafc, pordm_ageafc_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)  
roc_pordm_ageafc %>% ci.auc(method = "bootstrap")
# Age + AMH
pordm_ageamh_train <- porsm_train %>% select(Age, AMH, POR)
pordm_ageamh_test <- porsm_test %>% select(Age, AMH, POR)
pordm_ageamh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_ageamh_wf <- workflow() %>% add_model(pordm_ageamh_mod) %>% add_formula(POR ~.)
pordm_ageamh <- fit(pordm_ageamh_wf, data = pordm_ageamh_train)
roc_pordm_ageamh <- pordm_ageamh_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_ageamh, pordm_ageamh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)
roc_pordm_ageamh %>% ci.auc(method = "bootstrap")
# Age +AMH +FSH
pordm_ageamhfsh_train <- porsm_train %>% select(Age, AMH, FSH, POR)
pordm_ageamhfsh_test <- porsm_test %>% select(Age, AMH, FSH, POR)
pordm_ageamhfsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_ageamhfsh_wf <- workflow() %>% add_model(pordm_ageamhfsh_mod) %>% add_formula(POR ~.)
pordm_ageamhfsh <- fit(pordm_ageamhfsh_wf, data = pordm_ageamhfsh_train)
roc_pordm_ageamhfsh <- pordm_ageamhfsh_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_ageamhfsh, pordm_ageamhfsh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_pordm_ageamhfsh %>% ci.auc(method = "bootstrap")
# Age + AMH + AFC
pordm_ageamhafc_train <- porsm_train %>% select(Age, AMH, AFC, POR)
pordm_ageamhafc_test <- porsm_test %>% select(Age, AMH, AFC, POR)
pordm_ageamhafc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_ageamhafc_wf <- workflow() %>% add_model(pordm_ageamhafc_mod) %>% add_formula(POR ~.)
pordm_ageamhafc <- fit(pordm_ageamhafc_wf, data = pordm_ageamhafc_train)
roc_pordm_ageamhafc <- pordm_ageamhafc_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_ageamhafc, pordm_ageamhafc_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE) 
roc_pordm_ageamhafc %>% ci.auc(method = "bootstrap")
# Age + AMH + FSH +AFC
pordm_ageamhafcfsh_train <- porsm_train %>% select(Age, AMH, FSH, AFC, POR)
pordm_ageamhafcfsh_test <- porsm_test %>% select(Age, AMH, FSH, AFC, POR)
pordm_ageamhafcfsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
pordm_ageamhafcfsh_wf <- workflow() %>% add_model(pordm_ageamhafcfsh_mod) %>% add_formula(POR ~.)
pordm_ageamhafcfsh <- fit(pordm_ageamhafcfsh_wf, data = pordm_ageamhafcfsh_train)
roc_pordm_ageamhafcfsh <- pordm_ageamhafcfsh_test %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_ageamhafcfsh, pordm_ageamhafcfsh_test, type = "prob")) %>% 
  rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  roc(POR_r, .pred_No, auc = TRUE)  
roc_pordm_ageamhafcfsh %>% ci.auc(method = "bootstrap")
```

## Candidate HORDMs with conventional combinations of markers.
```{r}
set.seed(777)
# Age only.
hordm_age_train <- horsm_train %>% select(Age, HOR)
hordm_age_test <- horsm_test %>% select(Age, HOR)
hordm_age_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_age_wf <- workflow() %>% add_model(hordm_age_mod) %>% add_formula(HOR ~.)
hordm_age <- fit(hordm_age_wf, data = hordm_age_train)
roc_hordm_age <- hordm_age_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_age, hordm_age_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_age %>% ci.auc(method = "bootstrap")
# FSH only.
hordm_fsh_train <- horsm_train %>% select(FSH, HOR)
hordm_fsh_test <- horsm_test %>% select(FSH, HOR)
hordm_fsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_fsh_wf <- workflow() %>% add_model(hordm_fsh_mod) %>% add_formula(HOR ~.)
hordm_fsh <- fit(hordm_fsh_wf, data = hordm_fsh_train)
roc_hordm_fsh <- hordm_fsh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_fsh, hordm_fsh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_fsh %>% ci.auc(method = "bootstrap")
# AFC only.
hordm_afc_train <- horsm_train %>% select(AFC, HOR)
hordm_afc_test <- horsm_test %>% select(AFC, HOR)
hordm_afc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_afc_wf <- workflow() %>% add_model(hordm_afc_mod) %>% add_formula(HOR ~.)
hordm_afc <- fit(hordm_afc_wf, data = hordm_afc_train)
roc_hordm_afc <- hordm_afc_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_afc, hordm_afc_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_afc %>% ci.auc(method = "bootstrap")
# AMH only.
hordm_amh_train <- horsm_train %>% select(AMH, HOR)
hordm_amh_test <- horsm_test %>% select(AMH, HOR)
hordm_amh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_amh_wf <- workflow() %>% add_model(hordm_amh_mod) %>% add_formula(HOR ~.)
hordm_amh <- fit(hordm_amh_wf, data = hordm_amh_train)
roc_hordm_amh <- hordm_amh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_amh, hordm_amh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_amh %>% ci.auc(method = "bootstrap")
# Age + FSH
hordm_agefsh_train <- horsm_train %>% select(Age, FSH, HOR)
hordm_agefsh_test <- horsm_test %>% select(Age, FSH, HOR)
hordm_agefsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_agefsh_wf <- workflow() %>% add_model(hordm_agefsh_mod) %>% add_formula(HOR ~.)
hordm_agefsh <- fit(hordm_agefsh_wf, data = hordm_agefsh_train)
roc_hordm_agefsh <- hordm_agefsh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_agefsh, hordm_agefsh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_agefsh %>% ci.auc(method = "bootstrap")
# Age+ AFC
hordm_ageafc_train <- horsm_train %>% select(Age, AFC, HOR)
hordm_ageafc_test <- horsm_test %>% select(Age, AFC, HOR)
hordm_ageafc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_ageafc_wf <- workflow() %>% add_model(hordm_ageafc_mod) %>% add_formula(HOR ~.)
hordm_ageafc <- fit(hordm_ageafc_wf, data = hordm_ageafc_train)
roc_hordm_ageafc <- hordm_ageafc_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_ageafc, hordm_ageafc_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_ageafc %>% ci.auc(method = "bootstrap")
# Age + AMH
hordm_ageamh_train <- horsm_train %>% select(Age, AMH, HOR)
hordm_ageamh_test <- horsm_test %>% select(Age, AMH, HOR)
hordm_ageamh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_ageamh_wf <- workflow() %>% add_model(hordm_ageamh_mod) %>% add_formula(HOR ~.)
hordm_ageamh <- fit(hordm_ageamh_wf, data = hordm_ageamh_train)
roc_hordm_ageamh <- hordm_ageamh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_ageamh, hordm_ageamh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_ageamh %>% ci.auc(method = "bootstrap")
# Age +AMH +FSH
hordm_ageamhfsh_train <- horsm_train %>% select(Age, AMH, FSH, HOR)
hordm_ageamhfsh_test <- horsm_test %>% select(Age, AMH, FSH, HOR)
hordm_ageamhfsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_ageamhfsh_wf <- workflow() %>% add_model(hordm_ageamhfsh_mod) %>% add_formula(HOR ~.)
hordm_ageamhfsh <- fit(hordm_ageamhfsh_wf, data = hordm_ageamhfsh_train)
roc_hordm_ageamhfsh <- hordm_ageamhfsh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_ageamhfsh, hordm_ageamhfsh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_ageamhfsh %>% ci.auc(method = "bootstrap")
# Age + AMH + AFC
hordm_ageamhafc_train <- horsm_train %>% select(Age, AMH, AFC, HOR)
hordm_ageamhafc_test <- horsm_test %>% select(Age, AMH, AFC, HOR)
hordm_ageamhafc_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_ageamhafc_wf <- workflow() %>% add_model(hordm_ageamhafc_mod) %>% add_formula(HOR ~.)
hordm_ageamhafc <- fit(hordm_ageamhafc_wf, data = hordm_ageamhafc_train)
roc_hordm_ageamhafc <- hordm_ageamhafc_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_ageamhafc, hordm_ageamhafc_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_ageamhafc %>% ci.auc(method = "bootstrap")
# Age + AMH + FSH +AFC
hordm_ageamhafcfsh_train <- horsm_train %>% select(Age, AMH, FSH, AFC, HOR)
hordm_ageamhafcfsh_test <- horsm_test %>% select(Age, AMH, FSH, AFC, HOR)
hordm_ageamhafcfsh_mod <- boost_tree(trees = 100, mode = "classification", engine = "xgboost")
hordm_ageamhafcfsh_wf <- workflow() %>% add_model(hordm_ageamhafcfsh_mod) %>% add_formula(HOR ~.)
hordm_ageamhafcfsh <- fit(hordm_ageamhafcfsh_wf, data = hordm_ageamhafcfsh_train)
roc_hordm_ageamhafcfsh <- hordm_ageamhafcfsh_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_ageamhafcfsh, hordm_ageamhafcfsh_test, type = "prob")) %>% 
  rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_ageamhafcfsh %>% ci.auc(method = "bootstrap")
```

## Construction of candidate strategy and diagnostic models on complete data.
### Prepare data.
```{r}
porsm_complete <- openxlsx::read.xlsx("original_derivation.xlsx") %>% 
  as_tibble() %>% 
  select(AMH, AFC, POIorDOR, FSH, Age, P, Weight, DBP, WBC, ALT, RBC, Duration, LH, Protocol, Initial.FSH, Recombinant, Use.LH, POR) %>% 
  mutate_if(is.character, as.factor) %>% 
  filter(if_all(everything(), ~ !is.na(.x)))
horsm_complete <- openxlsx::read.xlsx("original_derivation.xlsx") %>% 
  as_tibble() %>% 
  select(AMH, AFC, FSH, Age, LH, POIorDOR, PCOS, PLT, Weight, Duration, Protocol, Initial.FSH, Recombinant, Use.LH, HOR) %>% 
  mutate_if(is.character, as.factor) %>% 
  filter(if_all(everything(), ~ !is.na(.x)))
pordm_complete <- porsm_complete %>% select(-c(Protocol, Recombinant, Initial.FSH, Use.LH))
hordm_complete <- horsm_complete %>% select(-c(Protocol, Recombinant, Initial.FSH, Use.LH))
```
### Candidate Strategy models on complete data.
```{r}
# Prepare data and some specifications.
set.seed(777)
porsm_complete_split <- initial_split(porsm_complete, strata = POR, prop = 0.8)
porsm_complete_train <- training(porsm_complete_split)
porsm_complete_test <- testing(porsm_complete_split)
horsm_complete_split <- initial_split(horsm_complete, strata = HOR, prop = 0.8)
horsm_complete_train <- training(horsm_complete_split)
horsm_complete_test <- testing(horsm_complete_split)
porsm_top_complete_rec <- recipe(POR ~., data = porsm_complete_train) %>% 
  step_dummy(all_nominal_predictors())
porsm_top_complete_wf <- workflow() %>% 
  add_model(porsm_top_mod) %>% 
  add_recipe(porsm_top_complete_rec) 
horsm_top_complete_rec <- recipe(HOR ~., data = horsm_complete_train) %>% 
  step_dummy(all_nominal_predictors())
horsm_top_complete_wf <- workflow() %>% 
  add_model(horsm_top_mod) %>% 
  add_recipe(horsm_top_complete_rec)
# Construction of PORM.
porsm_top_complete <- finalize_workflow(porsm_top_complete_wf, select_best(porsm_top_bayes_regs, "roc_auc")) %>% 
  fit(., data = porsm_complete_train)
# Construction of HORM.
horsm_top_complete <- finalize_workflow(horsm_top_complete_wf, select_best(horsm_top_bayes_regs, "roc_auc")) %>% 
  fit(., data = horsm_complete_train) 

# Calculate AUC and CI.
set.seed(777)
roc_porsm_complete <- porsm_complete_test%>% 
  select(POR) %>% 
  bind_cols(predict(porsm_top_complete, porsm_complete_test, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE)
roc_porsm_complete %>% pROC::ci.auc(method = "bootstrap")
roc_horsm_complete <- horsm_complete_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_top_complete, horsm_complete_test, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE)
roc_horsm_complete %>% pROC::ci.auc(method = "bootstrap")
```
### Candidate diagnostic models on complete data.
```{r}
set.seed(777)
pordm_complete_split <- initial_split(pordm_complete, strata = POR, prop = 0.8)
pordm_complete_train <- training(pordm_complete_split)
pordm_complete_test <- testing(pordm_complete_split)
hordm_complete_split <- initial_split(hordm_complete, strata = HOR, prop = 0.8)
hordm_complete_train <- training(hordm_complete_split)
hordm_complete_test <- testing(hordm_complete_split)
pordm_top_complete_rec <- recipe(POR ~., data = pordm_complete_train) %>% 
  step_dummy(all_nominal_predictors())
pordm_top_complete_wf <- workflow() %>% 
  add_model(pordm_top_mod) %>% 
  add_recipe(pordm_top_complete_rec) 
hordm_top_complete_rec <- recipe(HOR ~., data = hordm_complete_train) %>% 
  step_dummy(all_nominal_predictors())
hordm_top_complete_wf <- workflow() %>% 
  add_model(hordm_top_mod) %>% 
  add_recipe(hordm_top_complete_rec)
# Construction of PORM.
pordm_top_complete <- finalize_workflow(pordm_top_complete_wf, select_best(pordm_top_bayes_regs, "roc_auc")) %>% 
  fit(., data = pordm_complete_train)
# Construction of HORM.
hordm_top_complete <- finalize_workflow(hordm_top_complete_wf, select_best(hordm_top_bayes_regs, "roc_auc")) %>% 
  fit(., data = hordm_complete_train) 

# Calculate AUC and CI.
set.seed(777)
roc_pordm_complete <- pordm_complete_test%>% 
  select(POR) %>% 
  bind_cols(predict(pordm_top_complete, pordm_complete_test, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE)
roc_pordm_complete %>% pROC::ci.auc(method = "bootstrap")
roc_hordm_complete <- hordm_complete_test %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top_complete, hordm_complete_test, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE)
roc_hordm_complete %>% pROC::ci.auc(method = "bootstrap")
```

# AUC and CI: External validation of strategy models.
## Preparation of data.
```{r}
ningbo <- openxlsx::read.xlsx("validation for OR.xlsx") %>% 
  mutate_if(is.character, as.factor) %>% 
  select(AMH, AFC, POIorDOR, FSH, Age, P, Weight, DBP, WBC, ALT, RBC, Duration, LH, PCOS, PLT, PCOS, Protocol, Recombinant, Initial.FSH, Use.LH, POR, HOR) %>% 
  filter(if_all(everything(), ~ !is.na(.x)))

porsm_top_val <- ningbo %>% select(AMH, AFC, POIorDOR, FSH, Age, P, Weight, DBP, WBC, ALT, RBC, Duration, LH, Protocol, Initial.FSH, Recombinant, Use.LH, POR) %>% 
  recipe(POR ~.) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  step_mutate(Protocol_Short = 0,
              Recombinant_Yes = 0) %>% 
  prep() %>% bake(new_data = NULL)
horsm_top_val <- ningbo %>% select(AMH, AFC, FSH, Age, LH, POIorDOR, PCOS, PLT, Weight, Duration, Protocol, Initial.FSH, Recombinant, Use.LH, HOR) %>% 
  recipe(HOR ~.) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  step_mutate(Protocol_Short = 0,
              Recombinant_Yes = 0) %>% 
  prep() %>% bake(new_data = NULL) 
pordm_top_val <- ningbo %>% select(AMH, AFC, POIorDOR, FSH, Age, P, Weight, DBP, WBC, ALT, RBC, Duration, LH, POR) %>% 
  recipe(POR ~.) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  prep() %>% bake(new_data = NULL) 
hordm_top_val <- ningbo %>% select(AMH, AFC, FSH, Age, LH, POIorDOR, PCOS, PLT, Weight, Duration, HOR) %>% 
  recipe(HOR ~.) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  prep() %>% bake(new_data = NULL)
porsm_amh_val <- porsm_top_val %>% 
  mutate(AMH = NA) %>% mutate_if(is.logical, as.numeric)
horsm_amh_val <- horsm_top_val %>% 
  mutate(AMH = NA) %>% mutate_if(is.logical, as.numeric)
porsm_afc_val <- porsm_top_val %>% 
  mutate(AFC = NA) %>% mutate_if(is.logical, as.numeric)
horsm_afc_val <- horsm_top_val%>% 
  mutate(AFC = NA) %>% mutate_if(is.logical, as.numeric)
pordm_amh_val <- pordm_top_val %>% 
  mutate(AMH = NA) %>% mutate_if(is.logical, as.numeric)
hordm_amh_val <- hordm_top_val %>% 
  mutate(AMH = NA) %>% mutate_if(is.logical, as.numeric)
pordm_afc_val <- pordm_top_val %>% 
  mutate(AFC = NA) %>% mutate_if(is.logical, as.numeric)
hordm_afc_val <- hordm_top_val%>% 
  mutate(AFC = NA) %>% mutate_if(is.logical, as.numeric)
```
## External validation of strategy models.
```{r}
set.seed(777)
# PORSM
roc_porsm_top_external <- porsm_top_val %>% select(POR) %>% 
  bind_cols(predict(porsm_top_final, porsm_top_val, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE)  
roc_porsm_top_external %>% pROC::ci.auc(method = "bootstrap")
# HORSM
roc_horsm_top_external <- horsm_top_val %>% select(HOR) %>% 
  bind_cols(predict(horsm_top_final, horsm_top_val, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) 
roc_horsm_top_external %>% pROC::ci.auc(method = "bootstrap")
```
## External validation of diagnostic models.
```{r}
set.seed(777)
# PORDM
roc_pordm_top_external <- pordm_top_val %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_top_final, pordm_top_val, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) 
roc_pordm_top_external %>% pROC::ci.auc(method = "bootstrap")
# HORDM
roc_hordm_top_external <- hordm_top_val %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top_final, hordm_top_val, type = "prob")) %>%
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) 
roc_hordm_top_external %>% pROC::ci.auc(method = "bootstrap")
```

## Export roc objects for plotting ROC.
```{r}
roc_top_porsm <- roc_porsm_top
roc_top_horsm <- roc_horsm_top
roc_top_pordm <- roc_pordm_top
roc_top_hordm <- roc_hordm_top
roc_top_porsm_external <- roc_porsm_top_external
roc_top_horsm_external <- roc_horsm_top_external
roc_top_pordm_external <- roc_pordm_top_external
roc_top_hordm_external <- roc_hordm_top_external
save(roc_top_porsm, roc_top_horsm, roc_top_pordm, roc_top_hordm, roc_top_porsm_external, roc_top_horsm_external, roc_top_pordm_external, roc_top_hordm_external, 
     file = "roc_top.RData")
```
## AUC and CI: Stress tests.
### **PORSM**
```{r}
set.seed(777)
# Missing AFC
porsm_afc_val %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_top_final, porsm_afc_val, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Missing AMH
porsm_amh_val %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_top_final, porsm_amh_val, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Age >= 35
porsm_top_val %>% filter(Age >= 35) %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_top_final, porsm_top_val %>% filter(Age >= 35), type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Age < 35
porsm_top_val %>% filter(Age < 35) %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_top_final, porsm_top_val %>% filter(Age < 35), type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# AMH >1.1ng/mL or AFC > 5
porsm_top_val %>% filter(AMH > 1.1 | AFC > 5) %>% 
  select(POR) %>% 
  bind_cols(predict(porsm_top_final, porsm_top_val %>% filter(AMH > 1.1 | AFC > 5), type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
```
### **PORDM**
```{r}
set.seed(777)
# Missing AFC
pordm_afc_val %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_top_final, pordm_afc_val, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Missing AMH
pordm_amh_val %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_top_final, pordm_amh_val, type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Age >= 35
pordm_top_val %>% filter(Age >= 35) %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_top_final, pordm_top_val %>% filter(Age >= 35), type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Age < 35
pordm_top_val %>% filter(Age < 35) %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_top_final, pordm_top_val %>% filter(Age < 35), type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# AMH >1.1ng/mL or AFC > 5
pordm_top_val %>% filter(AMH > 1.1 | AFC > 5) %>% 
  select(POR) %>% 
  bind_cols(predict(pordm_top_final, pordm_top_val %>% filter(AMH > 1.1 | AFC > 5), type = "prob")) %>% 
  sjmisc::rec(POR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, POR)) %>% 
  pROC::roc(POR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
```

### **HORSM**
```{r}
set.seed(777)
# Missing AFC
horsm_afc_val %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_top_final, horsm_afc_val, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Missing AMH
horsm_amh_val %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_top_final, horsm_amh_val, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Age >= 35
horsm_top_val %>% filter(Age >= 35) %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_top_final, horsm_top_val %>% filter(Age >= 35), type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Age < 35
horsm_top_val %>% filter(Age < 35) %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_top_final, horsm_top_val %>% filter(Age < 35), type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# AMH <= 3.4ng/mL or AFC <= 24
horsm_top_val %>% filter(AMH <= 3.4 | AFC <= 24) %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm_top_final, horsm_top_val %>% filter(AMH <= 3.4 | AFC <= 24), type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
```
### **HORDM**
```{r}
set.seed(777)
# Missing AFC
hordm_afc_val %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top_final, hordm_afc_val, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Missing AMH
hordm_amh_val %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top_final, hordm_amh_val, type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Age >= 35
hordm_top_val %>% filter(Age >= 35) %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top_final, hordm_top_val %>% filter(Age >= 35), type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# Age < 35
hordm_top_val %>% filter(Age < 35) %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top_final, hordm_top_val %>% filter(Age < 35), type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# AMH <= 3.4ng/mL or AFC <= 24
hordm_top_val %>% filter(AMH <= 3.4 | AFC <= 24) %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top_final, hordm_top_val %>% filter(AMH <= 3.4 | AFC <= 24), type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
# AMH >3.4ng/mL or AFC > 24
hordm_top_val %>% filter(AMH > 3.4 | AFC > 24) %>% 
  select(HOR) %>% 
  bind_cols(predict(hordm_top_final, hordm_top_val %>% filter(AMH > 3.4 | AFC > 24), type = "prob")) %>% 
  sjmisc::rec(HOR, rec = "No = 0; Yes = 1") %>% 
  select(-c(.pred_Yes, HOR)) %>% 
  pROC::roc(HOR_r, .pred_No, auc = TRUE) %>% 
  pROC::ci.auc(method = "bootstrap")
```