Introduction to supervised machine learning in R: Classification

1 Exploratory Data Analysis (EDA)

Note: save and load Rstudio project (.Rproj) in directory containing the dataset

# setwd("<write/path/to/dataset/>") # if not loading Rproject file

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

data <- read_csv("countdata_large_normed.csv")

## Rows: 11204 Columns: 24
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (7): id, author, born, category, sex, title_monogr, topics
## dbl (17): num_sentences, subst, subst_prop, verb, adj, kolon, punkt, komma, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

inspect the first 6 rows of the data

head(data)

## # A tibble: 6 × 24
##   num_sentences subst subst_prop  verb    adj   kolon  punkt  komma     spm
##           <dbl> <dbl>      <dbl> <dbl>  <dbl>   <dbl>  <dbl>  <dbl>   <dbl>
## 1        0.0817 0.203    0.0949  0.152 0.102  0.0155  0.0596 0.0442 0.00442
## 2        0.0813 0.218    0.119   0.147 0.0923 0       0.0703 0.0462 0      
## 3        0.0622 0.237    0.0705  0.151 0.102  0.00207 0.0560 0.0353 0      
## 4        0.0721 0.232    0.0762  0.140 0.102  0       0.0621 0.0361 0.00200
## 5        0.0764 0.186    0.0306  0.205 0.0939 0.00437 0.0655 0.0393 0.00655
## 6        0.120  0.133    0.00451 0.228 0.0632 0.00226 0.115  0.0497 0      
## # ℹ 15 more variables: utrop <dbl>, overskrift <dbl>, ukjent <dbl>,
## #   mean_sent_length <dbl>, num_words <dbl>, id <chr>, section_number <dbl>,
## #   author <chr>, born <chr>, category <chr>, n_authors <dbl>,
## #   publication_year <dbl>, sex <chr>, title_monogr <chr>, topics <chr>

1.0.1 Description of dataset

data <- data |> 
  rename(noun = subst,
         prop_noun = subst_prop,
         colon = kolon,
         period = punkt,
         comma = komma,
         exclamation = utrop)

names(data)

##  [1] "num_sentences"    "noun"             "prop_noun"        "verb"            
##  [5] "adj"              "colon"            "period"           "comma"           
##  [9] "spm"              "exclamation"      "overskrift"       "ukjent"          
## [13] "mean_sent_length" "num_words"        "id"               "section_number"  
## [17] "author"           "born"             "category"         "n_authors"       
## [21] "publication_year" "sex"              "title_monogr"     "topics"

table(data$category)

## 
##   AV   SA   SK 
## 5205 3000 2999

names(data)

##  [1] "num_sentences"    "noun"             "prop_noun"        "verb"            
##  [5] "adj"              "colon"            "period"           "comma"           
##  [9] "spm"              "exclamation"      "overskrift"       "ukjent"          
## [13] "mean_sent_length" "num_words"        "id"               "section_number"  
## [17] "author"           "born"             "category"         "n_authors"       
## [21] "publication_year" "sex"              "title_monogr"     "topics"

input <- data |> 
  dplyr::select(category, num_sentences:comma, exclamation) 

names(input)

##  [1] "category"      "num_sentences" "noun"          "prop_noun"    
##  [5] "verb"          "adj"           "colon"         "period"       
##  [9] "comma"         "exclamation"

str(input)

## tibble [11,204 × 10] (S3: tbl_df/tbl/data.frame)
##  $ category     : chr [1:11204] "AV" "AV" "AV" "AV" ...
##  $ num_sentences: num [1:11204] 0.0817 0.0813 0.0622 0.0721 0.0764 ...
##  $ noun         : num [1:11204] 0.203 0.218 0.237 0.232 0.186 ...
##  $ prop_noun    : num [1:11204] 0.0949 0.1187 0.0705 0.0762 0.0306 ...
##  $ verb         : num [1:11204] 0.152 0.147 0.151 0.14 0.205 ...
##  $ adj          : num [1:11204] 0.1015 0.0923 0.1017 0.1022 0.0939 ...
##  $ colon        : num [1:11204] 0.01545 0 0.00207 0 0.00437 ...
##  $ period       : num [1:11204] 0.0596 0.0703 0.056 0.0621 0.0655 ...
##  $ comma        : num [1:11204] 0.0442 0.0462 0.0353 0.0361 0.0393 ...
##  $ exclamation  : num [1:11204] 0 0 0 0 0 0 0 0 0 0 ...

input$category <- factor(input$category)

ggplot(data=input, aes(x=category, y=prop_noun, fill=category))+
    geom_boxplot()

ggplot(data=input, aes(x=category, y=prop_noun, fill=category))+
    geom_boxplot()+
    ylim(0, 0.3)

## Warning: Removed 24 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

1.0.1.1 Exercise

ggplot(data=data, aes(x=category, y=noun, fill=category))+
    geom_boxplot()

facet_plot <- input |> 
  tidyr::pivot_longer(!category, names_to = "variable", values_to = "value")

facet_plot

## # A tibble: 100,836 × 3
##    category variable       value
##    <fct>    <chr>          <dbl>
##  1 AV       num_sentences 0.0817
##  2 AV       noun          0.203 
##  3 AV       prop_noun     0.0949
##  4 AV       verb          0.152 
##  5 AV       adj           0.102 
##  6 AV       colon         0.0155
##  7 AV       period        0.0596
##  8 AV       comma         0.0442
##  9 AV       exclamation   0     
## 10 AV       num_sentences 0.0813
## # ℹ 100,826 more rows

ggplot(data=facet_plot, aes(x=variable, y=value, fill=category))+
    geom_boxplot()

ggplot(data=facet_plot, aes(x=variable, y=value, fill=category))+
    geom_boxplot()+
    facet_wrap(~variable, scales="free")+
    theme(axis.text.x = element_blank(), #optional: remove x axis label (redundant) 
          axis.ticks.x = element_blank()) #optional: remove x axis tick

pairs(input[,2:ncol(input)], pch = 19, lower.panel = NULL) #pch is symbol type (maybe leave out)

if ( "GGally" %in% rownames(installed.packages()) ){
  
  library(GGally)
  
  input |>
  dplyr::select(-category) |>
  ggpairs(aes(alpha = 0.4))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
     
}

2 Binary Classification

2.1 Data preparation

input_bin <- input |> 
  filter(category != "SK")

table(input_bin$category)

## 
##   AV   SA   SK 
## 5205 3000    0

input_bin$category <- factor(input_bin$category)
table(input_bin$category)

## 
##   AV   SA 
## 5205 3000

2.2 Training and test sets

library(tidymodels)

## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──

## ✔ broom        1.0.6     ✔ rsample      1.2.1
## ✔ dials        1.3.0     ✔ tune         1.2.1
## ✔ infer        1.0.7     ✔ workflows    1.1.4
## ✔ modeldata    1.4.0     ✔ workflowsets 1.1.0
## ✔ parsnip      1.2.1     ✔ yardstick    1.3.1
## ✔ recipes      1.1.0

## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ scales::discard() masks purrr::discard()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ recipes::fixed()  masks stringr::fixed()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ yardstick::spec() masks readr::spec()
## ✖ recipes::step()   masks stats::step()
## • Search for functions across packages at https://www.tidymodels.org/find/

train_set_prop = 0.8
set.seed(561) #for reproducibility (i.e. ensures the "split" is the same for each rerun)

train_test_split <- rsample::initial_split(input_bin, prop=train_set_prop, strata=category)
train <- training(train_test_split)
test <- testing(train_test_split)

table(train$category)

## 
##   AV   SA 
## 4164 2400

table(test$category)

## 
##   AV   SA 
## 1041  600

2.3 XGBoost: Training and evaluating our machine learning model

library(xgboost)

## 
## Attaching package: 'xgboost'

## The following object is masked from 'package:dplyr':
## 
##     slice

#Note: within tidymodels, you get same result from
#xboost across runs
#set.seed(169) #XGBoost is "stochastic"

xgb_fit <-
  boost_tree(trees=50) |> #specify that you wanted a boosted tree model
  set_engine("xgboost") |> #specify that you want to use xgboost
  set_mode("classification") |> #specify the "objective" is classification
  fit(category ~ ., data=train)


xgb_bin <- xgb_fit #saving for later

xgb_pred = predict(xgb_fit, test)
head(xgb_pred)

## # A tibble: 6 × 1
##   .pred_class
##   <fct>      
## 1 AV         
## 2 AV         
## 3 AV         
## 4 AV         
## 5 AV         
## 6 AV

xgb_pred_prob = predict(xgb_fit, test, type="prob")
head(xgb_pred_prob)

## # A tibble: 6 × 2
##   .pred_AV .pred_SA
##      <dbl>    <dbl>
## 1    0.990  0.00957
## 2    0.919  0.0810 
## 3    0.793  0.207  
## 4    0.997  0.00309
## 5    0.993  0.00742
## 6    0.979  0.0215

xgb_res <- bind_cols(test$category, xgb_pred, xgb_pred_prob)

## New names:
## • `` -> `...1`

names(xgb_res) <- c("obs", "pred", "pred_AV", "pred_SA")

xgb_res$pred <- factor(xgb_res$pred) # Note: this seems to be "necessary" for R v4.0 and upwards

xgb_conf <- xgb_res |> 
  conf_mat(obs, pred)

xgb_conf

##           Truth
## Prediction   AV   SA
##         AV 1020   37
##         SA   21  563

2.3.1 Accuracy

xgb_acc <- sum(diag(xgb_conf$table))/sum(xgb_conf$table)
xgb_acc

## [1] 0.9646557

comp_acc <- function(conf_matrix){
  accuracy <- sum(diag(conf_matrix))/sum(conf_matrix)
  return(accuracy)
}

comp_acc(xgb_conf$table)

## [1] 0.9646557

2.3.2 Precision

comp_prec <- function(conf_matrix){
  precision <- diag(conf_matrix) / rowSums(conf_matrix)
  return(precision)
}

xgb_prec <- comp_prec(xgb_conf$table)
xgb_prec

##        AV        SA 
## 0.9649953 0.9640411

2.3.3 Recall

comp_recall <- function(conf_matrix){
  recall <- diag(conf_matrix) / colSums(conf_matrix) 
  return(recall)
}

xgb_rec <- comp_recall(xgb_conf$table)
xgb_rec

##        AV        SA 
## 0.9798271 0.9383333

xgb_res |> 
  accuracy(obs, pred) 

## # A tibble: 1 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.965

xgb_res |> 
  recall(obs, pred) 

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 recall  binary         0.980

xgb_res |> 
  precision(obs, pred) 

## # A tibble: 1 × 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 precision binary         0.965

levels(xgb_res$obs)

## [1] "AV" "SA"

xgb_res |> 
  recall(obs, pred, event_level = "second")

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 recall  binary         0.938

xgb_res |> 
  precision(obs, pred, event_level = "second") 

## # A tibble: 1 × 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 precision binary         0.964

xgb_res |> 
  metrics(obs, pred)

## # A tibble: 2 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.965
## 2 kap      binary         0.923

xgb_res |> 
  metrics(obs, pred, pred_AV)

## # A tibble: 4 × 3
##   .metric     .estimator .estimate
##   <chr>       <chr>          <dbl>
## 1 accuracy    binary         0.965
## 2 kap         binary         0.923
## 3 mn_log_loss binary         0.102
## 4 roc_auc     binary         0.992

2.3.4 Plots

ggplot(xgb_res, aes(x=obs, y=pred_SA, fill=obs))+
    geom_boxplot(outlier.shape = NA)+
    geom_jitter(width=0.1,alpha=0.1)+
    geom_hline(yintercept = 0.5, color = "blue", linetype = "dashed", alpha=0.6)+
    theme(axis.text.x = element_blank(), #ev. ekstra
          axis.ticks.x = element_blank())+ #ev. ekstra
    labs(fill="actual category",
        x="")    

ggplot(xgb_res, aes(x=obs, y=pred_AV, fill=obs))+
    geom_boxplot(outlier.shape = NA)+
    geom_jitter(width=0.1,alpha=0.2)+
    geom_hline(yintercept = 0.5, color = "blue", linetype = "dashed", alpha=0.6)+
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank())+
    labs(fill="actual category",
        x="")    

knitr::include_graphics("overfitting_runs.svg")

#https://commons.wikimedia.org/wiki/File:2d-epochs-overfitting.svg

2.3.5 Thresholds

threshold_new = 0.8
xgb_pred_new <- ifelse(xgb_pred_prob$.pred_SA >= threshold_new, "SA", "AV")

xgb_res_new <- bind_cols(test$category, xgb_pred_new, xgb_pred_prob)

## New names:
## • `` -> `...1`
## • `` -> `...2`

names(xgb_res_new) <- c("obs", "pred", "pred_AV", "pred_SA")

xgb_res_new$pred <- factor(xgb_res_new$pred) # Note: this seems to be "necessary" for R v4.0 and upwards

xgb_conf_new <- xgb_res_new |> 
  conf_mat(obs, pred)

xgb_conf_new

##           Truth
## Prediction   AV   SA
##         AV 1032   69
##         SA    9  531

comp_acc(xgb_conf_new$table) 

## [1] 0.952468

comp_prec(xgb_conf_new$table)

##        AV        SA 
## 0.9373297 0.9833333

comp_recall(xgb_conf_new$table)

##        AV        SA 
## 0.9913545 0.8850000

2.4 ROC and AUC

2.4.1 ROC

roc_curve(xgb_res, obs, pred_AV) |> 
    autoplot()

roc_curve(xgb_res, obs, pred_AV) |>
  ggplot(aes(x = 1 - specificity, y = sensitivity)) +
  geom_line()+ #in this case, same result with geom_path (order in dataframe, NOT x-axis)
  geom_abline(lty = 3) +
  coord_equal() +
  theme_bw()

2.4.2 AUC

auc_xgb_bin <- roc_auc(xgb_res, obs, pred_AV)
auc_xgb_bin

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 roc_auc binary         0.992

2.5 Feature importance (binary classification)

library(vip)

## 
## Attaching package: 'vip'

## The following object is masked from 'package:utils':
## 
##     vi

vip(xgb_fit)

xgb_features <- vip::vi(xgb_bin)
xgb_features

## # A tibble: 9 × 2
##   Variable      Importance
##   <chr>              <dbl>
## 1 colon             0.366 
## 2 exclamation       0.226 
## 3 prop_noun         0.217 
## 4 num_sentences     0.0475
## 5 noun              0.0374
## 6 adj               0.0314
## 7 verb              0.0277
## 8 period            0.0253
## 9 comma             0.0212

xgb_feature_bin <-
ggplot(data=xgb_features, aes(x=reorder(Variable, -Importance),
                              y=Importance))+
    geom_col(aes(fill=Variable))+
    theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
    labs(x="",
        title="Binary classification\n('AV' vs 'SA')")

xgb_feature_bin

3 Multiclass classification

train_set_prop = 0.8
set.seed(231) #for reproducibility

train_test_multi <- rsample::initial_split(input, prop=train_set_prop, strata=category)
train <- training(train_test_multi)
test <- testing(train_test_multi)

#set.seed(1245) #XGBoost is "stochastic"

table(train$category)

## 
##   AV   SA   SK 
## 4164 2400 2399

table(test$category)

## 
##   AV   SA   SK 
## 1041  600  600

xgb_multi <-
  boost_tree(trees=50) |>
  set_engine("xgboost") |>
  set_mode("classification") |>
  fit(category ~ ., data=train)

xgb_pred = predict(xgb_multi, test)
xgb_pred_prob = predict(xgb_multi, test, type="prob")

xgb_pred

## # A tibble: 2,241 × 1
##    .pred_class
##    <fct>      
##  1 AV         
##  2 AV         
##  3 AV         
##  4 AV         
##  5 AV         
##  6 AV         
##  7 AV         
##  8 AV         
##  9 AV         
## 10 AV         
## # ℹ 2,231 more rows

xgb_pred_prob

## # A tibble: 2,241 × 3
##    .pred_AV .pred_SA  .pred_SK
##       <dbl>    <dbl>     <dbl>
##  1    0.964 0.0357   0.000779 
##  2    0.940 0.0490   0.0107   
##  3    0.980 0.0193   0.000341 
##  4    0.993 0.00636  0.000583 
##  5    0.916 0.0260   0.0581   
##  6    0.999 0.000663 0.000166 
##  7    0.999 0.000678 0.0000611
##  8    0.958 0.0202   0.0221   
##  9    0.809 0.152    0.0392   
## 10    0.972 0.0259   0.00175  
## # ℹ 2,231 more rows

xgb_res <- bind_cols(test$category, xgb_pred, xgb_pred_prob)

## New names:
## • `` -> `...1`

names(xgb_res) <- c("obs", "pred", "pred_AV", "pred_SA", "pred_SK")

xgb_res$pred <- factor(xgb_res$pred) # Note: this seems to be "necessary" for R v4.0 and upwards

xgb_conf <- xgb_res |> 
  conf_mat(obs, pred)

xgb_conf

##           Truth
## Prediction   AV   SA   SK
##         AV 1015   59    9
##         SA   15  477   35
##         SK   11   64  556

3.0.1 Accuracy

xgb_res |> 
    accuracy(obs, pred)

## # A tibble: 1 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy multiclass     0.914

3.0.2 Precision

multi_prec <- comp_prec(xgb_conf$table)
multi_prec

##        AV        SA        SK 
## 0.9372114 0.9051233 0.8811410

3.0.3 Recall

multi_recall <- comp_recall(xgb_conf$table)
multi_recall 

##        AV        SA        SK 
## 0.9750240 0.7950000 0.9266667

3.0.4 Macro-averaged metrics

macroPrecision <- mean(multi_prec)
macroPrecision

## [1] 0.9078253

precision(xgb_res, obs, pred)

## # A tibble: 1 × 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 precision macro          0.908

macroRecall <- mean(multi_recall)
macroRecall

## [1] 0.8988969

recall(xgb_res, obs, pred)

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 recall  macro          0.899

3.0.5 Multiclass AUC

xgb_res |> 
    roc_auc(obs, pred_AV:pred_SK) 

## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 roc_auc hand_till      0.977

xgb_res |> 
  metrics(obs, pred, pred_AV:pred_SK)

## # A tibble: 4 × 3
##   .metric     .estimator .estimate
##   <chr>       <chr>          <dbl>
## 1 accuracy    multiclass     0.914
## 2 kap         multiclass     0.865
## 3 mn_log_loss multiclass     0.233
## 4 roc_auc     hand_till      0.977

3.0.6 Optional: Plotting predicted probability with actual category

test_plot <- xgb_res |> 
  dplyr::select(obs,pred_AV, pred_SA, pred_SK)

test_plot_long <- test_plot |> 
  tidyr::pivot_longer(!obs, names_to = "variable", values_to = "value")

head(test_plot_long)

## # A tibble: 6 × 3
##   obs   variable    value
##   <fct> <chr>       <dbl>
## 1 AV    pred_AV  0.964   
## 2 AV    pred_SA  0.0357  
## 3 AV    pred_SK  0.000779
## 4 AV    pred_AV  0.940   
## 5 AV    pred_SA  0.0490  
## 6 AV    pred_SK  0.0107

ggplot(data=test_plot_long, aes(x=obs, y=value, fill=obs))+
    geom_boxplot(outlier.shape = NA)+
    geom_jitter(width=0.1,alpha=0.05)+
    geom_hline(yintercept = 1/3, color = "blue", linetype = "dashed", alpha=0.6)+
    facet_wrap(~variable)+
    labs(y="Predicted probability",
        fill="Actual category")

3.1 Feature importance (multiclass classification)

vip(xgb_multi)

xgb_features <- vip::vi(xgb_multi)

xgb_feature_multi <-
ggplot(data=xgb_features, aes(x=reorder(Variable, -Importance),
                              y=Importance))+
    geom_col(aes(fill=Variable))+
    theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
    labs(x="",
        title="Multiclass")


xgb_feature_multi

library("patchwork", quietly = TRUE)

features <- xgb_feature_bin + xgb_feature_multi & theme(legend.position = "bottom")


features + plot_layout(guides = "collect") +
    plot_annotation(title = "Feature importance from XGB")

3.2 Exercise (optional)

min_n = 2 #choose an integer from range 1-10 (default = 1)
tree_depth = 9  #choose integer from range 1-10 (default = 6)
learn_rate =  0.15 # choose decimal value between 0.01 and 0.5 (default = 0.3)

2. Calculate model evaluation metrics (accuracy, macro-precision, macro-recall, auc)

4 K-fold cross-validation

4.1 Cross-validation for binary classification

knitr::include_graphics("CV_image.svg")

#https://upload.wikimedia.org/wikipedia/commons/1/1c/K-fold_cross_validation_EN.jpg

xgb_mod <- boost_tree(trees=50) |>
  set_engine("xgboost") |>
  set_mode("classification")

wflow_xgb <-
  workflow() |> 
  add_model(xgb_mod) |> 
  add_formula(category ~ .)

set.seed(156) 
folds <- vfold_cv(input_bin, v = 10, strata = category) #strata for stratified cv
xgb_cv <- fit_resamples(wflow_xgb, folds)

xgb_cv

## # Resampling results
## # 10-fold cross-validation using stratification 
## # A tibble: 10 × 4
##    splits             id     .metrics         .notes          
##    <list>             <chr>  <list>           <list>          
##  1 <split [7384/821]> Fold01 <tibble [3 × 4]> <tibble [0 × 3]>
##  2 <split [7384/821]> Fold02 <tibble [3 × 4]> <tibble [0 × 3]>
##  3 <split [7384/821]> Fold03 <tibble [3 × 4]> <tibble [0 × 3]>
##  4 <split [7384/821]> Fold04 <tibble [3 × 4]> <tibble [0 × 3]>
##  5 <split [7384/821]> Fold05 <tibble [3 × 4]> <tibble [0 × 3]>
##  6 <split [7385/820]> Fold06 <tibble [3 × 4]> <tibble [0 × 3]>
##  7 <split [7385/820]> Fold07 <tibble [3 × 4]> <tibble [0 × 3]>
##  8 <split [7385/820]> Fold08 <tibble [3 × 4]> <tibble [0 × 3]>
##  9 <split [7385/820]> Fold09 <tibble [3 × 4]> <tibble [0 × 3]>
## 10 <split [7385/820]> Fold10 <tibble [3 × 4]> <tibble [0 × 3]>

4.1.1 Accuracy

collect_metrics(xgb_cv)

## # A tibble: 3 × 6
##   .metric     .estimator   mean     n  std_err .config             
##   <chr>       <chr>       <dbl> <int>    <dbl> <chr>               
## 1 accuracy    binary     0.960     10 0.00199  Preprocessor1_Model1
## 2 brier_class binary     0.0312    10 0.00137  Preprocessor1_Model1
## 3 roc_auc     binary     0.990     10 0.000949 Preprocessor1_Model1

# the same as
# xgb_cv |> 
#   collect_metrics()

collect_metrics(xgb_cv, summarize=FALSE) |> 
  filter(.metric =="accuracy")

## # A tibble: 10 × 5
##    id     .metric  .estimator .estimate .config             
##    <chr>  <chr>    <chr>          <dbl> <chr>               
##  1 Fold01 accuracy binary         0.950 Preprocessor1_Model1
##  2 Fold02 accuracy binary         0.960 Preprocessor1_Model1
##  3 Fold03 accuracy binary         0.954 Preprocessor1_Model1
##  4 Fold04 accuracy binary         0.957 Preprocessor1_Model1
##  5 Fold05 accuracy binary         0.962 Preprocessor1_Model1
##  6 Fold06 accuracy binary         0.972 Preprocessor1_Model1
##  7 Fold07 accuracy binary         0.956 Preprocessor1_Model1
##  8 Fold08 accuracy binary         0.963 Preprocessor1_Model1
##  9 Fold09 accuracy binary         0.961 Preprocessor1_Model1
## 10 Fold10 accuracy binary         0.966 Preprocessor1_Model1

#the same as
# xgb_cv |> 
#   collect_metrics(summarize=FALSE) |> 
#   filter(.metric =="accuracy")

xgb_cv |> 
  collect_metrics(summarize=FALSE) |> 
  filter(.metric =="accuracy") |> 
  summarise(max = max(.estimate), min = min(.estimate))

## # A tibble: 1 × 2
##     max   min
##   <dbl> <dbl>
## 1 0.972 0.950

xgb_cv_res <- 
  xgb_cv |> 
  collect_metrics(summarize=FALSE) |> 
  filter(.metric =="accuracy") |> 
  summarise(max = max(.estimate), min = min(.estimate), 
            mean = mean(.estimate), sd = sd(.estimate))

xgb_cv_res$mean + xgb_cv_res$sd

## [1] 0.9664288

xgb_cv_res$mean - xgb_cv_res$sd

## [1] 0.953868

4.1.2 Other metrics

my_metrics <- metric_set(accuracy,roc_auc, precision,
                         recall)

xgb_cv_rev <- fit_resamples(wflow_xgb, folds,
                            metrics = my_metrics)



xgb_cv_rev |> 
  collect_metrics()

## # A tibble: 4 × 6
##   .metric   .estimator  mean     n  std_err .config             
##   <chr>     <chr>      <dbl> <int>    <dbl> <chr>               
## 1 accuracy  binary     0.960    10 0.00199  Preprocessor1_Model1
## 2 precision binary     0.956    10 0.00257  Preprocessor1_Model1
## 3 recall    binary     0.982    10 0.00178  Preprocessor1_Model1
## 4 roc_auc   binary     0.990    10 0.000949 Preprocessor1_Model1

4.1.3 Out-of-sample predictions

keep_pred <- control_resamples(save_pred = TRUE, save_workflow = TRUE)

set.seed(156) 
xgb_cv_rev <- fit_resamples(wflow_xgb, folds,
                            metrics = my_metrics, control = keep_pred)

collect_predictions(xgb_cv_rev)

## # A tibble: 8,205 × 7
##    .pred_class .pred_AV .pred_SA id      .row category .config             
##    <fct>          <dbl>    <dbl> <chr>  <int> <fct>    <chr>               
##  1 AV             0.991 0.00885  Fold01     8 AV       Preprocessor1_Model1
##  2 AV             0.969 0.0309   Fold01    16 AV       Preprocessor1_Model1
##  3 AV             0.979 0.0213   Fold01    30 AV       Preprocessor1_Model1
##  4 AV             0.999 0.000631 Fold01    31 AV       Preprocessor1_Model1
##  5 AV             1.00  0.000283 Fold01    41 AV       Preprocessor1_Model1
##  6 AV             0.998 0.00226  Fold01    47 AV       Preprocessor1_Model1
##  7 AV             0.998 0.00230  Fold01    68 AV       Preprocessor1_Model1
##  8 AV             0.998 0.00204  Fold01    90 AV       Preprocessor1_Model1
##  9 AV             0.974 0.0263   Fold01    96 AV       Preprocessor1_Model1
## 10 AV             0.967 0.0331   Fold01   101 AV       Preprocessor1_Model1
## # ℹ 8,195 more rows

xgb_cv_conf <- xgb_cv_rev |> 
  collect_predictions() |> 
  conf_mat(category, .pred_class)

xgb_cv_conf

##           Truth
## Prediction   AV   SA
##         AV 5113  235
##         SA   92 2765

comp_prec(xgb_cv_conf$table)

##        AV        SA 
## 0.9560583 0.9677984

comp_recall(xgb_cv_conf$table)

##        AV        SA 
## 0.9823247 0.9216667

xgb_cv_rev |> 
  collect_metrics(summarize=FALSE) |> 
  filter(.metric=="accuracy") |> 
  summarise(mean = mean(.estimate)) |> 
  pull()

## [1] 0.9601484

comp_acc(xgb_cv_conf$table)

## [1] 0.9601463

4.2 Repeated cross-validation

rep_folds <- vfold_cv(input_bin, v = 10, repeats = 50,
                      strata = category) #strata for stratified cv

xgb_mod_cv <- boost_tree(trees=5) |>
  set_engine("xgboost") |>
  set_mode("classification")

wflow_xgb_cv <-
  workflow() |> 
  add_model(xgb_mod_cv) |> 
  add_formula(category ~ .)

# If doParallel is installed, and your computer has
# several cores, this will make use of parallelization
# (splitting the task into bits to be run simultaneously)

if ("doParallel" %in% rownames(installed.packages())){
  library("doParallel")
  all_cores <- parallel::detectCores(logical = FALSE)
  cl <- makePSOCKcluster(all_cores) #number of cores?
  registerDoParallel(cl) #only version > 1.0.16..
}

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

## Loading required package: iterators

## Loading required package: parallel

set.seed(777)
xgb_cv_rep <- fit_resamples(wflow_xgb_cv, rep_folds, 
                            metrics = my_metrics, control = keep_pred)

if (exists("cl")){
  stopCluster(cl) 
}

xgb_cv_rep

## # Resampling results
## # 10-fold cross-validation repeated 50 times using stratification 
## # A tibble: 500 × 6
##    splits             id       id2    .metrics         .notes   .predictions
##    <list>             <chr>    <chr>  <list>           <list>   <list>      
##  1 <split [7384/821]> Repeat01 Fold01 <tibble [4 × 4]> <tibble> <tibble>    
##  2 <split [7384/821]> Repeat01 Fold02 <tibble [4 × 4]> <tibble> <tibble>    
##  3 <split [7384/821]> Repeat01 Fold03 <tibble [4 × 4]> <tibble> <tibble>    
##  4 <split [7384/821]> Repeat01 Fold04 <tibble [4 × 4]> <tibble> <tibble>    
##  5 <split [7384/821]> Repeat01 Fold05 <tibble [4 × 4]> <tibble> <tibble>    
##  6 <split [7385/820]> Repeat01 Fold06 <tibble [4 × 4]> <tibble> <tibble>    
##  7 <split [7385/820]> Repeat01 Fold07 <tibble [4 × 4]> <tibble> <tibble>    
##  8 <split [7385/820]> Repeat01 Fold08 <tibble [4 × 4]> <tibble> <tibble>    
##  9 <split [7385/820]> Repeat01 Fold09 <tibble [4 × 4]> <tibble> <tibble>    
## 10 <split [7385/820]> Repeat01 Fold10 <tibble [4 × 4]> <tibble> <tibble>    
## # ℹ 490 more rows

xgb_cv_rep |>
  collect_metrics()

## # A tibble: 4 × 6
##   .metric   .estimator  mean     n  std_err .config             
##   <chr>     <chr>      <dbl> <int>    <dbl> <chr>               
## 1 accuracy  binary     0.923   500 0.000418 Preprocessor1_Model1
## 2 precision binary     0.917   500 0.000479 Preprocessor1_Model1
## 3 recall    binary     0.967   500 0.000398 Preprocessor1_Model1
## 4 roc_auc   binary     0.969   500 0.000286 Preprocessor1_Model1

# xgb_cv_rep |>
#   collect_metrics(summarize=FALSE) |>
#   group_by(.metric) |>
#   summarise(max = max(.estimate), min = min(.estimate),
#             mean = mean(.estimate), sd = sd(.estimate))

nrow(xgb_cv_rep)

## [1] 500

xgb_cv_rep |>
  collect_metrics(summarize = FALSE) |> 
  filter(.metric == "accuracy") |> 
  ggplot(aes(x=.estimate))+
  geom_histogram()+
  labs(x="accuracy")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

xgb_cv_rep |> 
    collect_predictions() |> 
    group_by(id) |>  #id stores nth repeat
    accuracy(category, .pred_class) |> 
    ggplot(aes(x=.estimate))+
      geom_histogram()+
      labs(x="accuracy")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

xgb_cv_rep |> 
    collect_predictions() |> 
    group_by(id) |> 
    accuracy(category, .pred_class) |> 
    summarise(min =min(.estimate), max = max(.estimate))

## # A tibble: 1 × 2
##     min   max
##   <dbl> <dbl>
## 1 0.919 0.926

xgb_cv_rep |> 
    collect_predictions() |> 
    group_by(id) |> 
    accuracy(category, .pred_class) |> 
    summarise(sd = sd(.estimate))

## # A tibble: 1 × 1
##        sd
##     <dbl>
## 1 0.00177

4.3 Multiclass classification with cross-validation

rep_folds <- vfold_cv(input, v = 10, repeats = 25,
                      strata = category) #strata for stratified cv


if ("doParallel" %in% rownames(installed.packages())){
  all_cores <- parallel::detectCores(logical = FALSE)
  cl <- makePSOCKcluster(all_cores) #number of cores?
  registerDoParallel(cl) #only version > 1.0.16..
}

set.seed(456)
xgb_cv_rep <- fit_resamples(wflow_xgb_cv, rep_folds, 
                            metrics = my_metrics,
                            control = keep_pred)


if (exists("cl")){
  stopCluster(cl) 
}

Exercise:
1) Get mean accuracy, auc, precision, recall
2)Plot a histogram for all the accuracy estimates

xgb_cv_rep |> 
  collect_metrics()

## # A tibble: 4 × 6
##   .metric   .estimator  mean     n  std_err .config             
##   <chr>     <chr>      <dbl> <int>    <dbl> <chr>               
## 1 accuracy  multiclass 0.874   250 0.000596 Preprocessor1_Model1
## 2 precision macro      0.867   250 0.000663 Preprocessor1_Model1
## 3 recall    macro      0.856   250 0.000684 Preprocessor1_Model1
## 4 roc_auc   hand_till  0.962   250 0.000284 Preprocessor1_Model1

xgb_cv_rep |> 
  collect_predictions() |> 
  group_by(id) |> 
  accuracy(category, .pred_class) |> 
  ggplot(aes(x=.estimate))+
  geom_histogram()+
  labs(x="accuracy")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

4.3.1 Impact of sample size

size = nrow(input) * 0.05
set.seed(126)
small <- dplyr::sample_n(input, size)
nrow(small)

## [1] 560

rep_folds <- vfold_cv(small, v = 10, repeats = 50,
                      strata = category) #strata for stratified cv



if ("doParallel" %in% rownames(installed.packages())){
  all_cores <- parallel::detectCores(logical = FALSE)
  cl <- makePSOCKcluster(all_cores) #number of cores?
  registerDoParallel(cl) #only version > 1.0.16..
}

set.seed(123)
cv_result_small <- fit_resamples(wflow_xgb_cv, rep_folds, 
                            metrics = my_metrics,
                            control = keep_pred)

if (exists("cl")){
  stopCluster(cl) 
}

cv_result_small |> 
  collect_metrics()

## # A tibble: 4 × 6
##   .metric   .estimator  mean     n std_err .config             
##   <chr>     <chr>      <dbl> <int>   <dbl> <chr>               
## 1 accuracy  multiclass 0.814   500 0.00228 Preprocessor1_Model1
## 2 precision macro      0.805   500 0.00251 Preprocessor1_Model1
## 3 recall    macro      0.802   500 0.00246 Preprocessor1_Model1
## 4 roc_auc   hand_till  0.928   500 0.00134 Preprocessor1_Model1

# cv_result_small |> 
#   collect_metrics(summarize=FALSE) |> 
#   group_by(.metric) |> 
#   summarise(mean = mean(.estimate))

Here you see a large variability in accuracy

cv_result_small |>
  collect_metrics(summarize = FALSE) |> 
  filter(.metric == "accuracy") |> 
  ggplot(aes(x=.estimate))+
  geom_histogram()+
  labs(x="accuracy")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

cv_result_small |> 
  collect_predictions() |> 
  group_by(id) |> 
  accuracy(category, .pred_class) |> 
  ggplot(aes(x=.estimate))+
  geom_histogram()+
    labs(x="accuracy")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

5 Preprocessing data with “recipe”

# if need to split into training and test set again
# train_set_prop = 0.8
# set.seed(231) #for reproducibility
# 
# train_test_multi <- rsample::initial_split(input, prop=train_set_prop, strata=category)
# train <- training(train_test_multi)
# test <- testing(train_test_multi)

norm_recipe <- recipe(category ~ ., data = train) |>
  step_normalize(all_predictors())

library("discrim")

## 
## Attaching package: 'discrim'

## The following object is masked from 'package:dials':
## 
##     smoothness

lda_mod <-
  discrim_linear() |>
  set_engine("MASS")

lda_wfl <- workflow() |> 
  add_recipe(norm_recipe) |> 
  add_model(lda_mod)

lda_fit_norm <- fit(lda_wfl, data = train)

lda_norm_pred = predict(lda_fit_norm, test)
# head(lda_norm_pred)

train_norm <- norm_recipe |> 
  prep(training = train) |> 
  bake(train)

test_norm  <- norm_recipe |> 
  prep(training = train) |> 
  bake(test)

head(train_norm)

## # A tibble: 6 × 10
##   num_sentences    noun prop_noun   verb     adj  colon  period   comma
##           <dbl>   <dbl>     <dbl>  <dbl>   <dbl>  <dbl>   <dbl>   <dbl>
## 1         0.766 -0.0745     0.757 -0.418  0.350   2.37   0.137  -0.196 
## 2         0.746  0.255      1.30  -0.606 -0.0879 -0.743  0.868  -0.0858
## 3         0.236  0.594      0.327 -0.866  0.381  -0.743  0.309  -0.638 
## 4         2.88  -1.67      -1.31   2.40  -1.47   -0.288  3.92    0.106 
## 5         0.171  0.387     -0.651  1.84  -0.702  -0.743 -0.0838  0.0163
## 6         0.997  0.0895     0.304 -1.46   1.03    0.122  1.05   -0.968 
## # ℹ 2 more variables: exclamation <dbl>, category <fct>

head(test_norm)

## # A tibble: 6 × 10
##   num_sentences   noun prop_noun    verb     adj  colon period  comma
##           <dbl>  <dbl>     <dbl>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl>
## 1        -0.315  0.686    0.199  -0.450   0.355  -0.325 -0.108 -0.682
## 2         0.474 -0.473   -0.716   1.55   -0.0131  0.137  0.539 -0.461
## 3         0.313  1.18    -0.0212  0.874  -0.878   0.104  0.227 -1.69 
## 4         1.63   2.04    -0.576   0.184  -0.470  -0.307  2.11  -0.129
## 5        -0.374 -0.231    0.468   0.0384 -1.36   -0.318 -0.618  0.738
## 6        -1.31  -0.667    1.26    0.354   0.781  -0.743 -1.32  -0.740
## # ℹ 2 more variables: exclamation <dbl>, category <fct>

# other_recipe <- recipe(category ~ ., data = input) |>
#   step_normalize(all_predictors())

other_recipe <- recipe(category ~ ., data = input) |>
  step_normalize(all_predictors())

train_other <- other_recipe |> 
  prep(training = input) |> 
  bake(train)

test_other  <- other_recipe |> 
  prep(training = input) |> 
  bake(test)

head(train_other)

## # A tibble: 6 × 10
##   num_sentences    noun prop_noun   verb     adj  colon  period   comma
##           <dbl>   <dbl>     <dbl>  <dbl>   <dbl>  <dbl>   <dbl>   <dbl>
## 1         0.750 -0.0761     0.747 -0.412  0.358   2.33   0.129  -0.194 
## 2         0.731  0.254      1.29  -0.601 -0.0810 -0.736  0.849  -0.0837
## 3         0.227  0.594      0.319 -0.860  0.389  -0.736  0.298  -0.639 
## 4         2.83  -1.67      -1.31   2.40  -1.46   -0.288  3.86    0.110 
## 5         0.163  0.387     -0.654  1.84  -0.697  -0.736 -0.0889  0.0189
## 6         0.978  0.0883     0.295 -1.45   1.04    0.116  1.03   -0.971 
## # ℹ 2 more variables: exclamation <dbl>, category <fct>

head(test_other)

## # A tibble: 6 × 10
##   num_sentences   noun prop_noun    verb      adj   colon period  comma
##           <dbl>  <dbl>     <dbl>   <dbl>    <dbl>   <dbl>  <dbl>  <dbl>
## 1        -0.316  0.686    0.191  -0.445   0.363   -0.324  -0.112 -0.683
## 2         0.462 -0.475   -0.719   1.56   -0.00594  0.131   0.525 -0.461
## 3         0.303  1.19    -0.0275  0.877  -0.873    0.0985  0.217 -1.70 
## 4         1.60   2.04    -0.579   0.189  -0.465   -0.307   2.07  -0.128
## 5        -0.374 -0.233    0.459   0.0431 -1.36    -0.317  -0.616  0.744
## 6        -1.30  -0.670    1.24    0.359   0.791   -0.736  -1.31  -0.742
## # ℹ 2 more variables: exclamation <dbl>, category <fct>

6 Build and evaluate multiple models simultaneously

# norm_recipe <- recipe(category ~ ., data = train) |>
#   step_normalize(all_predictors())

# lda_mod <-
#   discrim_linear() |>
#   set_engine("MASS")

xgb_mod <- boost_tree(trees=50) |>
  set_engine("xgboost") |>
  set_mode("classification")


library("kernlab")

## 
## Attaching package: 'kernlab'

## The following object is masked from 'package:scales':
## 
##     alpha

## The following object is masked from 'package:purrr':
## 
##     cross

## The following object is masked from 'package:ggplot2':
## 
##     alpha

svm_mod <- svm_linear(mode="classification") |> 
  set_engine("kernlab")

library(workflowsets)

all_models_set <-
  workflow_set(
    preproc = list(norm = norm_recipe),
    models = list(lda = lda_mod, svm = svm_mod,
                  xgb = xgb_mod)
  )

if ("doParallel" %in% rownames(installed.packages())){
  
  unregister <- function() {
    env <- foreach:::.foreachGlobals
    rm(list=ls(name=env), pos=env)
  }
  
  unregister()
}


folds <- vfold_cv(input, v = 10, strata = category)



all_models <-
  all_models_set |>
  workflow_map(
    fn = "fit_resamples",
    resamples = folds,
    verbose = TRUE,
    seed = 1234
  )

## i 1 of 3 resampling: norm_lda

## ✔ 1 of 3 resampling: norm_lda (850ms)

## i 2 of 3 resampling: norm_svm

## ✔ 2 of 3 resampling: norm_svm (1m 2.5s)

## i 3 of 3 resampling: norm_xgb

## ✔ 3 of 3 resampling: norm_xgb (21.6s)

all_models |>
  collect_metrics()

## # A tibble: 9 × 9
##   wflow_id .config         preproc model .metric .estimator   mean     n std_err
##   <chr>    <chr>           <chr>   <chr> <chr>   <chr>       <dbl> <int>   <dbl>
## 1 norm_lda Preprocessor1_… recipe  disc… accura… multiclass 0.759     10 0.00430
## 2 norm_lda Preprocessor1_… recipe  disc… brier_… multiclass 0.169     10 0.00164
## 3 norm_lda Preprocessor1_… recipe  disc… roc_auc hand_till  0.901     10 0.00190
## 4 norm_svm Preprocessor1_… recipe  svm_… accura… multiclass 0.762     10 0.00342
## 5 norm_svm Preprocessor1_… recipe  svm_… brier_… multiclass 0.165     10 0.00175
## 6 norm_svm Preprocessor1_… recipe  svm_… roc_auc hand_till  0.903     10 0.00204
## 7 norm_xgb Preprocessor1_… recipe  boos… accura… multiclass 0.916     10 0.00262
## 8 norm_xgb Preprocessor1_… recipe  boos… brier_… multiclass 0.0607    10 0.00179
## 9 norm_xgb Preprocessor1_… recipe  boos… roc_auc hand_till  0.982     10 0.00117

autoplot(all_models)

6.1 Statistically compare models

all_auc_wide <- all_models |>
  collect_metrics(summarize = FALSE) |>
  filter(.metric == "roc_auc") |>
  dplyr::select(wflow_id, .estimate, id) |>
  pivot_wider(id_cols = "id", names_from = "wflow_id", values_from = ".estimate")

all_auc_wide

## # A tibble: 10 × 4
##    id     norm_lda norm_svm norm_xgb
##    <chr>     <dbl>    <dbl>    <dbl>
##  1 Fold01    0.903    0.901    0.985
##  2 Fold02    0.888    0.890    0.975
##  3 Fold03    0.902    0.905    0.988
##  4 Fold04    0.906    0.909    0.984
##  5 Fold05    0.903    0.907    0.979
##  6 Fold06    0.898    0.899    0.983
##  7 Fold07    0.898    0.904    0.982
##  8 Fold08    0.899    0.900    0.980
##  9 Fold09    0.911    0.915    0.986
## 10 Fold10    0.899    0.902    0.979

all_auc_wide |>
  with( t.test(norm_lda, norm_svm, paired = TRUE))

## 
##  Paired t-test
## 
## data:  norm_lda and norm_svm
## t = -3.8411, df = 9, p-value = 0.00396
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0038345089 -0.0009920063
## sample estimates:
## mean of the differences 
##            -0.002413258

all_auc_wide |>
  with( t.test(norm_lda, norm_xgb, paired = TRUE))

## 
##  Paired t-test
## 
## data:  norm_lda and norm_xgb
## t = -58.45, df = 9, p-value = 6.326e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.08430949 -0.07802672
## sample estimates:
## mean of the differences 
##              -0.0811681

all_auc_wide |>
  with( t.test(norm_svm, norm_xgb, paired = TRUE))

## 
##  Paired t-test
## 
## data:  norm_svm and norm_xgb
## t = -49.625, df = 9, p-value = 2.749e-12
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.08234492 -0.07516477
## sample estimates:
## mean of the differences 
##             -0.07875485

7 Hyperparameter tuning

7.1 Automatic grid search

xgb_spec <- boost_tree(
  trees = 10, 
  tree_depth = tune(),
  loss_reduction = tune()                     
) |> 
  set_engine("xgboost") |>  
  set_mode("classification")

xgb_wf <- workflow() |>
  add_formula(category ~ .) |>
  add_model(xgb_spec)


rep_folds <- vfold_cv(input_bin, v = 10, strata = category, repeats = 5)


tree_grid <- grid_regular(tree_depth(), loss_reduction(), levels=2)

tree_grid

## # A tibble: 4 × 2
##   tree_depth loss_reduction
##        <int>          <dbl>
## 1          1   0.0000000001
## 2         15   0.0000000001
## 3          1  31.6         
## 4         15  31.6

if ("doParallel" %in% rownames(installed.packages())){
  all_cores <- parallel::detectCores(logical = FALSE)
  cl <- makePSOCKcluster(all_cores) #number of cores?
  registerDoParallel(cl) #only version > 1.0.16..
}

set.seed(234)
xgb_tune <- tune_grid(
  xgb_wf,
  resamples = rep_folds,
  grid = tree_grid
)


if (exists("cl")){
  stopCluster(cl) 
}

collect_metrics(xgb_tune)

## # A tibble: 12 × 8
##    tree_depth loss_reduction .metric     .estimator   mean     n std_err .config
##         <int>          <dbl> <chr>       <chr>       <dbl> <int>   <dbl> <chr>  
##  1          1   0.0000000001 accuracy    binary     0.813     50 2.10e-3 Prepro…
##  2          1   0.0000000001 brier_class binary     0.151     50 6.13e-4 Prepro…
##  3          1   0.0000000001 roc_auc     binary     0.884     50 1.96e-3 Prepro…
##  4         15   0.0000000001 accuracy    binary     0.947     50 1.13e-3 Prepro…
##  5         15   0.0000000001 brier_class binary     0.0408    50 7.12e-4 Prepro…
##  6         15   0.0000000001 roc_auc     binary     0.983     50 6.25e-4 Prepro…
##  7          1  31.6          accuracy    binary     0.813     50 2.10e-3 Prepro…
##  8          1  31.6          brier_class binary     0.151     50 6.13e-4 Prepro…
##  9          1  31.6          roc_auc     binary     0.884     50 1.96e-3 Prepro…
## 10         15  31.6          accuracy    binary     0.938     50 1.29e-3 Prepro…
## 11         15  31.6          brier_class binary     0.0494    50 7.12e-4 Prepro…
## 12         15  31.6          roc_auc     binary     0.979     50 7.17e-4 Prepro…

show_best(xgb_tune, metric = "roc_auc")

## # A tibble: 4 × 8
##   tree_depth loss_reduction .metric .estimator  mean     n  std_err .config     
##        <int>          <dbl> <chr>   <chr>      <dbl> <int>    <dbl> <chr>       
## 1         15   0.0000000001 roc_auc binary     0.983    50 0.000625 Preprocesso…
## 2         15  31.6          roc_auc binary     0.979    50 0.000717 Preprocesso…
## 3          1   0.0000000001 roc_auc binary     0.884    50 0.00196  Preprocesso…
## 4          1  31.6          roc_auc binary     0.884    50 0.00196  Preprocesso…

best_model <- select_best(xgb_tune, metric = "roc_auc")

best_model

## # A tibble: 1 × 3
##   tree_depth loss_reduction .config             
##        <int>          <dbl> <chr>               
## 1         15   0.0000000001 Preprocessor1_Model2

final_xgb <- finalize_workflow(
  xgb_wf,
  best_model
)

final_xgb

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: boost_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## category ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Boosted Tree Model Specification (classification)
## 
## Main Arguments:
##   trees = 10
##   tree_depth = 15
##   loss_reduction = 1e-10
## 
## Computational engine: xgboost

7.2 Manual grid search

# alternatively, can use base R expand.grid function (same paramater setup)

manual_grid <- tidyr::crossing(tree_depth = c(5, 10),
                               loss_reduction = c(2, 4))

manual_grid

## # A tibble: 4 × 2
##   tree_depth loss_reduction
##        <dbl>          <dbl>
## 1          5              2
## 2          5              4
## 3         10              2
## 4         10              4

if ("doParallel" %in% rownames(installed.packages())){
  all_cores <- parallel::detectCores(logical = FALSE)
  cl <- makePSOCKcluster(all_cores) #number of cores?
  registerDoParallel(cl) #only version > 1.0.16..
}

set.seed(155)
xgb_tune_manual <- tune_grid(
  xgb_wf,
  resamples = rep_folds,
  grid = manual_grid
)

if (exists("cl")){
  stopCluster(cl) 
}

show_best(xgb_tune_manual, metric = "roc_auc")

## # A tibble: 4 × 8
##   tree_depth loss_reduction .metric .estimator  mean     n  std_err .config     
##        <dbl>          <dbl> <chr>   <chr>      <dbl> <int>    <dbl> <chr>       
## 1         10              2 roc_auc binary     0.984    50 0.000631 Preprocesso…
## 2         10              4 roc_auc binary     0.984    50 0.000577 Preprocesso…
## 3          5              2 roc_auc binary     0.980    50 0.000697 Preprocesso…
## 4          5              4 roc_auc binary     0.980    50 0.000685 Preprocesso…

xgb_tune_manual |>
  collect_metrics() |>
  mutate(tree_depth = factor(tree_depth)) |>
  ggplot(aes(loss_reduction, mean, color = tree_depth)) +
    geom_line(size = 1.5, alpha = 0.6) +
    geom_point(size = 2) +
    facet_wrap(~ .metric, scales = "free", nrow = 2)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

8 Brief example of predicting a continuous variable

mod_cont_recipe <- 
  recipe(num_sentences ~ ., data=input) |> 
  step_dummy(all_nominal())

xgb_mod_cont <- boost_tree(trees=50) |>
  set_engine("xgboost") |>
  set_mode("regression")

wflow_xgb_cont <-
  workflow() |> 
  add_model(xgb_mod_cont) |> 
  #add_formula(num_sentences ~ .)
  add_recipe(mod_cont_recipe)

# input_cont <- input |> 
#   dplyr::select(-c(category))

if ("doParallel" %in% rownames(installed.packages())){
  all_cores <- parallel::detectCores(logical = FALSE)
  cl <- makePSOCKcluster(all_cores) 
  registerDoParallel(cl) 
}

set.seed(156) 
folds_cont <- vfold_cv(input, v = 10) 
cont_metrics <- metric_set(rmse, rsq, mae)


keep_pred <- control_resamples(save_pred = TRUE)

xgb_cv_cont <- fit_resamples(wflow_xgb_cont, folds_cont,
                             metrics = cont_metrics,
                             control = keep_pred)

if (exists("cl")){
  stopCluster(cl) 
}

collect_predictions(xgb_cv_cont)

## # A tibble: 11,204 × 5
##     .pred id      .row num_sentences .config             
##     <dbl> <chr>  <int>         <dbl> <chr>               
##  1 0.0756 Fold01     5        0.0764 Preprocessor1_Model1
##  2 0.0739 Fold01    15        0.0768 Preprocessor1_Model1
##  3 0.0892 Fold01    23        0.0841 Preprocessor1_Model1
##  4 0.0982 Fold01    25        0.0960 Preprocessor1_Model1
##  5 0.105  Fold01    41        0.0905 Preprocessor1_Model1
##  6 0.0715 Fold01    42        0.0720 Preprocessor1_Model1
##  7 0.0689 Fold01    50        0.0602 Preprocessor1_Model1
##  8 0.0704 Fold01    81        0.0725 Preprocessor1_Model1
##  9 0.0911 Fold01    82        0.0864 Preprocessor1_Model1
## 10 0.0689 Fold01    84        0.0669 Preprocessor1_Model1
## # ℹ 11,194 more rows

collect_predictions(xgb_cv_cont) |> 
  ggplot(aes(x=num_sentences, y=.pred))+
    geom_point(alpha=0.5)+
    geom_smooth(method="lm", se = TRUE)+
    coord_fixed()

## `geom_smooth()` using formula = 'y ~ x'

xgb_cv_cont |> 
  collect_metrics()

## # A tibble: 3 × 6
##   .metric .estimator    mean     n   std_err .config             
##   <chr>   <chr>        <dbl> <int>     <dbl> <chr>               
## 1 mae     standard   0.00472    10 0.0000394 Preprocessor1_Model1
## 2 rmse    standard   0.00708    10 0.000108  Preprocessor1_Model1
## 3 rsq     standard   0.849      10 0.00280   Preprocessor1_Model1

cor(collect_predictions(xgb_cv_cont)$num_sentences, 
    collect_predictions(xgb_cv_cont)$.pred)

## [1] 0.9214295