R Markdown

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
library(correlationfunnel)
## Warning: package 'correlationfunnel' was built under R version 4.4.2
## ══ Using correlationfunnel? ════════════════════════════════════════════════════
## You might also be interested in applied data science training for business.
## </> Learn more at - www.business-science.io </>
data <- readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/main/data/2023/2023-08-15/spam.csv')
## Rows: 4601 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): yesno
## dbl (6): crl.tot, dollar, bang, money, n000, make
## 
## ℹ 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.

Explore data

skimr::skim(data)
Data summary
Name data
Number of rows 4601
Number of columns 7
_______________________
Column type frequency:
character 1
numeric 6
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
yesno 0 1 1 1 0 2 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
crl.tot 0 1 283.29 606.35 1 35 95 266.00 15841.00 ▇▁▁▁▁
dollar 0 1 0.08 0.25 0 0 0 0.05 6.00 ▇▁▁▁▁
bang 0 1 0.27 0.82 0 0 0 0.32 32.48 ▇▁▁▁▁
money 0 1 0.09 0.44 0 0 0 0.00 12.50 ▇▁▁▁▁
n000 0 1 0.10 0.35 0 0 0 0.00 5.45 ▇▁▁▁▁
make 0 1 0.10 0.31 0 0 0 0.00 4.54 ▇▁▁▁▁
data_clean <- data %>%
  mutate(yesno = factor(yesno, levels = c("y", "n")))

Comparing the occurrence of the dollar/money sign and how it relates to emails being spam

data_clean %>%
    filter(dollar > 0, money > 0) %>%
    ggplot(aes(money, dollar, color = yesno)) +
    geom_point(alpha = 0.3, size = 1.0)

data_clean %>%
    ggplot(aes(dollar, y = after_stat(density), fill = yesno)) +
    geom_histogram(position = "identity", alpha = 1.0)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

data_clean %>% 
    ggplot(aes(money, y = after_stat(density), fill = yesno)) +
    geom_histogram(position = "identity", alpha = 1.0)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Step 1: Binarize 
data_binarized <- data_clean %>%
    select(-dollar) %>%
    binarize()

data_binarized %>% glimpse()
## Rows: 4,601
## Columns: 15
## $ `crl.tot__-Inf_35` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, …
## $ crl.tot__35_95     <dbl> 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ crl.tot__95_266    <dbl> 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, …
## $ crl.tot__266_Inf   <dbl> 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, …
## $ `bang__-Inf_0.315` <dbl> 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, …
## $ bang__0.315_Inf    <dbl> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, …
## $ money__0           <dbl> 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, …
## $ `money__-OTHER`    <dbl> 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, …
## $ n000__0            <dbl> 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, …
## $ `n000__-OTHER`     <dbl> 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, …
## $ make__0            <dbl> 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, …
## $ make__0.1          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `make__-OTHER`     <dbl> 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, …
## $ yesno__y           <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ yesno__n           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
# Step 2: Correlation
data_correlation <- data_binarized %>%
    correlate(yesno__y)

data_correlation
## # A tibble: 15 × 3
##    feature bin        correlation
##    <fct>   <chr>            <dbl>
##  1 yesno   y               1     
##  2 yesno   n              -1     
##  3 bang    -Inf_0.315     -0.490 
##  4 bang    0.315_Inf       0.490 
##  5 money   -OTHER          0.475 
##  6 money   0              -0.475 
##  7 n000    0              -0.419 
##  8 n000    -OTHER          0.419 
##  9 crl.tot -Inf_35        -0.360 
## 10 crl.tot 266_Inf         0.299 
## 11 make    0              -0.239 
## 12 make    -OTHER          0.223 
## 13 crl.tot 95_266          0.145 
## 14 crl.tot 35_95          -0.0818
## 15 make    0.1             0.0803
# Step 3: Plot
data_correlation %>%
    correlationfunnel::plot_correlation_funnel()

library(tidymodels)
## Warning: package 'tidymodels' was built under R version 4.4.2
## ── 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.2
## ✔ recipes      1.1.0
## Warning: package 'dials' was built under R version 4.4.2
## Warning: package 'infer' was built under R version 4.4.2
## Warning: package 'modeldata' was built under R version 4.4.2
## Warning: package 'parsnip' was built under R version 4.4.2
## Warning: package 'tune' was built under R version 4.4.2
## Warning: package 'workflows' was built under R version 4.4.2
## Warning: package 'workflowsets' was built under R version 4.4.3
## Warning: package 'yardstick' was built under R version 4.4.2
## ── 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()
## • Learn how to get started at https://www.tidymodels.org/start/

Split Data

set.seed(1234)

# data_clean <- data_clean %>% sample_n(100)

data_split <- initial_split(data_clean, strata = yesno)
data_train <- training(data_split)
data_test <- testing(data_split)

data_cv <- rsample::vfold_cv(data_train, strata = yesno)
data_cv
## #  10-fold cross-validation using stratification 
## # A tibble: 10 × 2
##    splits             id    
##    <list>             <chr> 
##  1 <split [3104/346]> Fold01
##  2 <split [3105/345]> Fold02
##  3 <split [3105/345]> Fold03
##  4 <split [3105/345]> Fold04
##  5 <split [3105/345]> Fold05
##  6 <split [3105/345]> Fold06
##  7 <split [3105/345]> Fold07
##  8 <split [3105/345]> Fold08
##  9 <split [3105/345]> Fold09
## 10 <split [3106/344]> Fold10

#Preprocess Data

library(themis)
## Warning: package 'themis' was built under R version 4.4.3
data_rec <- recipes::recipe(yesno ~., data = data_train) %>%
    step_YeoJohnson(all_nominal_predictors())

data_rec %>% prep() %>% juice() %>% glimpse()
## Rows: 3,450
## Columns: 7
## $ crl.tot <dbl> 26, 65, 1866, 35, 150, 28, 16, 791, 26, 27, 50, 16, 28, 13, 11…
## $ dollar  <dbl> 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,…
## $ bang    <dbl> 0.149, 0.262, 0.000, 0.000, 0.000, 0.393, 0.729, 0.149, 0.000,…
## $ money   <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
## $ n000    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
## $ make    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
## $ yesno   <fct> n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n,…

Build Model

library(usemodels)
## Warning: package 'usemodels' was built under R version 4.4.2
data <- data%>% 
  mutate(yesno = as.numeric(yesno))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `yesno = as.numeric(yesno)`.
## Caused by warning:
## ! NAs introduced by coercion
usemodels::use_xgboost(yesno ~ ., data = data_train)
## xgboost_recipe <- 
##   recipe(formula = yesno ~ ., data = data_train) %>% 
##   step_zv(all_predictors()) 
## 
## xgboost_spec <- 
##   boost_tree(trees = tune(), min_n = tune(), tree_depth = tune(), learn_rate = tune(), 
##     loss_reduction = tune(), sample_size = tune()) %>% 
##   set_mode("classification") %>% 
##   set_engine("xgboost") 
## 
## xgboost_workflow <- 
##   workflow() %>% 
##   add_recipe(xgboost_recipe) %>% 
##   add_model(xgboost_spec) 
## 
## set.seed(96026)
## xgboost_tune <-
##   tune_grid(xgboost_workflow, resamples = stop("add your rsample object"), grid = stop("add number of candidate points"))
xgboost_spec <- 
  boost_tree(trees = tune()) %>%
  set_mode("classification") %>% 
  set_engine("xgboost") 

xgboost_workflow <- 
  workflow() %>% 
  add_recipe(data_rec) %>% 
  add_model(xgboost_spec) 

Tune Hyperparameters

doParallel::registerDoParallel()

set.seed(50926)
xgboost_tune <-
  tune_grid(xgboost_workflow, 
            resamples = data_cv, 
            grid = 5, 
            control = control_grid(save_pred = TRUE))

Model Evaluation

Identify optimal value

collect_metrics(xgboost_tune) 
## # A tibble: 15 × 7
##    trees .metric     .estimator  mean     n std_err .config             
##    <int> <chr>       <chr>      <dbl> <int>   <dbl> <chr>               
##  1   253 accuracy    binary     0.872    10 0.00421 Preprocessor1_Model1
##  2   253 brier_class binary     0.103    10 0.00333 Preprocessor1_Model1
##  3   253 roc_auc     binary     0.912    10 0.00344 Preprocessor1_Model1
##  4   601 accuracy    binary     0.868    10 0.00368 Preprocessor1_Model2
##  5   601 brier_class binary     0.110    10 0.00314 Preprocessor1_Model2
##  6   601 roc_auc     binary     0.904    10 0.00313 Preprocessor1_Model2
##  7  1190 accuracy    binary     0.864    10 0.00396 Preprocessor1_Model3
##  8  1190 brier_class binary     0.115    10 0.00312 Preprocessor1_Model3
##  9  1190 roc_auc     binary     0.900    10 0.00318 Preprocessor1_Model3
## 10  1512 accuracy    binary     0.863    10 0.00419 Preprocessor1_Model4
## 11  1512 brier_class binary     0.117    10 0.00320 Preprocessor1_Model4
## 12  1512 roc_auc     binary     0.898    10 0.00327 Preprocessor1_Model4
## 13  1976 accuracy    binary     0.862    10 0.00440 Preprocessor1_Model5
## 14  1976 brier_class binary     0.119    10 0.00326 Preprocessor1_Model5
## 15  1976 roc_auc     binary     0.897    10 0.00339 Preprocessor1_Model5
collect_predictions(xgboost_tune) %>%
  group_by(id) %>%
  roc_curve(yesno, .pred_y) %>%  # Ensure you're using the probability for class 1
  autoplot()

xgboost_last <- xgboost_workflow %>%
    finalize_workflow(select_best(xgboost_tune, metric = "accuracy")) %>%
    last_fit(data_split)
## Warning: package 'xgboost' was built under R version 4.4.2
collect_metrics(xgboost_last)
## # A tibble: 3 × 4
##   .metric     .estimator .estimate .config             
##   <chr>       <chr>          <dbl> <chr>               
## 1 accuracy    binary        0.871  Preprocessor1_Model1
## 2 roc_auc     binary        0.914  Preprocessor1_Model1
## 3 brier_class binary        0.0993 Preprocessor1_Model1
collect_predictions(xgboost_last) %>%
    yardstick::conf_mat(yesno, .pred_class) %>%
    autoplot()

library(vip)
## Warning: package 'vip' was built under R version 4.4.3
## 
## Attaching package: 'vip'
## The following object is masked from 'package:utils':
## 
##     vi
xgboost_last %>%
    workflows::extract_fit_engine() %>%
    vip()