GLM_ridge.Rmd

---
title: "GLM-ridge"
author: "Guiquan"
date: "2021/7/12"
output:
  pdf_document: default
  html_document: default
editor_options:
  chunk_output_type: inline
---

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

Development of GLM-ridge-based candidate strategy models.

# Bayesian Optimization for Hyper Parameters.
## **PORSM**
```{r}
library(tidymodels)
# Data preparation and specifications.
porsm_mod_spec <- 
  logistic_reg(mixture = 0, penalty = tune()) %>% 
  set_engine("glmnet") 
porsm_workflow_spec <- workflow() %>% 
  add_model(porsm_mod_spec) %>% 
  add_formula(POR ~.)
porsm_bayes_params <- parameters(porsm_workflow_spec) %>% 
  update(penalty = penalty(c(-4, 1), trans = log10_trans()))
# Create cross-validation resamples.
set.seed(777)
porsm_valset <- vfold_cv(porsm_train, v = 5, strata = POR)
# Start tuning.
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
porsm_bayes_regs <- tune_bayes(
  porsm_workflow_spec,
  resamples = porsm_valset,
  param_info = porsm_bayes_params,
  iter = 100,
  metrics = metric_set(roc_auc, mn_log_loss),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)
```

## **HORSM**
```{r}
# Data preparation and specifications.
horsm_mod_spec <- 
  logistic_reg(mixture = 0, penalty = tune()) %>% 
  set_engine("glmnet") 
horsm_workflow_spec <- workflow() %>% 
  add_model(horsm_mod_spec) %>% 
  add_formula(HOR ~.)
horsm_bayes_params <- parameters(horsm_workflow_spec) %>% 
  update(penalty = penalty(c(-4, 1), trans = log10_trans()))
# Create cross-validation resamples.
set.seed(777)
horsm_valset <- vfold_cv(horsm_train, v = 5, strata = HOR)
# Start tuning.
cl <- parallel::makeCluster(8)
doParallel::registerDoParallel(cl)
horsm_bayes_regs <- tune_bayes(
  horsm_workflow_spec,
  resamples = horsm_valset,
  param_info = horsm_bayes_params,
  iter = 100,
  metrics = metric_set(roc_auc, mn_log_loss),
  control = control_bayes(no_improve = 50, verbose = TRUE)
)
parallel::stopCluster(cl)
```

# Construction of Strategy models and calculate AUC/Brier score.
## Developing models.
```{r}
set.seed(777)
# Strategy models on train data.
porsm <- select_best(porsm_bayes_regs, "roc_auc") %>% 
  finalize_workflow(porsm_workflow_spec, .) %>% 
  fit(., porsm_train)
horsm <- select_best(horsm_bayes_regs, "roc_auc") %>% 
  finalize_workflow(horsm_workflow_spec, .) %>% 
  fit(., horsm_train)
```

## Calculate AUC and Brier score.
```{r}
# AUC and 95%CI.
set.seed(777)
roc_porsm <- porsm_test %>% 
  select(POR) %>% 
  bind_cols(predict(porsm, porsm_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 %>% pROC::ci.auc(method = "bootstrap")
roc_horsm <- horsm_test %>% 
  select(HOR) %>% 
  bind_cols(predict(horsm, horsm_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 %>% pROC::ci.auc(method = "bootstrap")

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

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

# Export roc objects for plotting ROC.
```{r}
roc_ridge_porsm <- roc_porsm
roc_ridge_horsm <- roc_horsm
save(roc_ridge_porsm, roc_ridge_horsm, 
     file = "roc_ridge.RData")
```