Import data

library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.3.3
## ── 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.3.3
## ══ correlationfunnel Tip #2 ════════════════════════════════════════════════════
## Clean your NA's prior to using `binarize()`.
## Missing values and cleaning data are critical to getting great correlations. :)
data <- readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-04-27/departures.csv')
## Rows: 9423 Columns: 19
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr   (8): coname, exec_fullname, interim_coceo, still_there, notes, sources...
## dbl  (10): dismissal_dataset_id, gvkey, fyear, co_per_rol, departure_code, c...
## dttm  (1): leftofc
## 
## ℹ 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.

Clean data

skimr::skim(data)
Data summary
Name data
Number of rows 9423
Number of columns 19
_______________________
Column type frequency:
character 8
numeric 10
POSIXct 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
coname 0 1.00 2 30 0 3860 0
exec_fullname 0 1.00 5 790 0 8701 0
interim_coceo 9105 0.03 6 7 0 6 0
still_there 7311 0.22 3 10 0 77 0
notes 1644 0.83 5 3117 0 7755 0
sources 1475 0.84 18 1843 0 7915 0
eight_ks 4499 0.52 69 3884 0 4914 0
_merge 0 1.00 11 11 0 1 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
dismissal_dataset_id 0 1.00 5684.10 25005.46 1 2305.5 4593 6812.5 559044 ▇▁▁▁▁
gvkey 0 1.00 40132.48 53921.34 1004 7337.0 14385 60900.5 328795 ▇▁▁▁▁
fyear 0 1.00 2007.74 8.19 1987 2000.0 2008 2016.0 2020 ▁▆▅▅▇
co_per_rol 0 1.00 25580.22 18202.38 -1 8555.5 22980 39275.5 64602 ▇▆▅▃▃
departure_code 1667 0.82 5.20 1.53 1 5.0 5 7.0 9 ▁▃▇▅▁
ceo_dismissal 1813 0.81 0.20 0.40 0 0.0 0 0.0 1 ▇▁▁▁▂
tenure_no_ceodb 0 1.00 1.03 0.17 0 1.0 1 1.0 3 ▁▇▁▁▁
max_tenure_ceodb 0 1.00 1.05 0.24 1 1.0 1 1.0 4 ▇▁▁▁▁
fyear_gone 1802 0.81 2006.64 13.63 1980 2000.0 2007 2013.0 2997 ▇▁▁▁▁
cik 245 0.97 741469.17 486551.43 1750 106413.0 857323 1050375.8 1808065 ▆▁▇▂▁

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
leftofc 1802 0.81 1981-01-01 2998-04-27 2006-12-31 3627 
factors_vec <- data %>% select(departure_code, leftofc) %>% names()

data_clean <- data %>%

    # Address factors imported as numeric
    mutate(across(all_of(factors_vec), as.factor)) %>%
    
    # Remove NA
    drop_na(ceo_dismissal, tenure_no_ceodb, fyear_gone, departure_code, leftofc) %>%
 
    # Drop zero-variance variables
    select(-c(interim_coceo, eight_ks, gvkey, cik, coname, exec_fullname, sources, notes, "_merge",leftofc)) %>%
    
    # Drop still_there due to high missing values
    select(-still_there) %>%
    
    # Convert dismissal_dataset_id to character 
    mutate(dismissal_dataset_id = as.character(dismissal_dataset_id)) %>%
    distinct(dismissal_dataset_id, .keep_all = TRUE) %>% 
    
    mutate(ceo_dismissal = as.factor(ceo_dismissal))
skimr::skim(data_clean)
Data summary
Name data_clean
Number of rows 7476
Number of columns 8
_______________________
Column type frequency:
character 1
factor 2
numeric 5
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
dismissal_dataset_id 0 1 1 6 0 7476 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
departure_code 0 1 FALSE 7 5: 3577, 7: 2029, 3: 1315, 4: 194
ceo_dismissal 0 1 FALSE 2 0: 5993, 1: 1483

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
fyear 0 1 2005.61 7.45 1987 1999.0 2006.0 2012.00 2020 ▁▇▆▇▆
co_per_rol 0 1 21444.46 16350.23 -1 6975.5 18267.5 33414.25 64601 ▇▅▅▂▁
tenure_no_ceodb 0 1 1.03 0.16 1 1.0 1.0 1.00 3 ▇▁▁▁▁
max_tenure_ceodb 0 1 1.05 0.23 1 1.0 1.0 1.00 4 ▇▁▁▁▁
fyear_gone 0 1 2006.54 13.69 1980 2000.0 2006.0 2013.00 2997 ▇▁▁▁▁

Explore data

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

data_binarized %>% glimpse()
## Rows: 7,476
## Columns: 29
## $ dismissal_dataset_id__1        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `dismissal_dataset_id__-OTHER` <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ `fyear__-Inf_1999`             <dbl> 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, …
## $ fyear__1999_2006               <dbl> 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, …
## $ fyear__2006_2012               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, …
## $ fyear__2012_Inf                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `co_per_rol__-Inf_6975.5`      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ co_per_rol__6975.5_18267.5     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ co_per_rol__18267.5_33414.25   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ co_per_rol__33414.25_Inf       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ departure_code__1              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ departure_code__2              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ departure_code__3              <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, …
## $ departure_code__4              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ departure_code__5              <dbl> 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, …
## $ departure_code__6              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ departure_code__7              <dbl> 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, …
## $ ceo_dismissal__0               <dbl> 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, …
## $ ceo_dismissal__1               <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, …
## $ tenure_no_ceodb__1             <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ tenure_no_ceodb__2             <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `tenure_no_ceodb__-OTHER`      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ max_tenure_ceodb__1            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ max_tenure_ceodb__2            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `max_tenure_ceodb__-OTHER`     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `fyear_gone__-Inf_2000`        <dbl> 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, …
## $ fyear_gone__2000_2006          <dbl> 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, …
## $ fyear_gone__2006_2013          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, …
## $ fyear_gone__2013_Inf           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
# Step 2: Correlation
data_correlation <- data_binarized %>%
    correlate(ceo_dismissal__0)

data_correlation
## # A tibble: 29 × 3
##    feature        bin         correlation
##    <fct>          <chr>             <dbl>
##  1 ceo_dismissal  0                1     
##  2 ceo_dismissal  1               -1     
##  3 departure_code 3               -0.929 
##  4 departure_code 5                0.476 
##  5 departure_code 7                0.304 
##  6 departure_code 4               -0.273 
##  7 departure_code 6                0.0786
##  8 fyear          -Inf_1999        0.0775
##  9 co_per_rol     -Inf_6975.5      0.0595
## 10 fyear_gone     -Inf_2000        0.0585
## # ℹ 19 more rows
# Step 3: Plot
data_correlation %>%
    correlationfunnel::plot_correlation_funnel()

Model Building

Split Data

library(tidymodels)
## Warning: package 'tidymodels' was built under R version 4.3.3
## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
## ✔ broom        1.0.5     ✔ 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
## Warning: package 'dials' was built under R version 4.3.3
## Warning: package 'infer' was built under R version 4.3.3
## Warning: package 'modeldata' was built under R version 4.3.3
## Warning: package 'parsnip' was built under R version 4.3.3
## Warning: package 'recipes' was built under R version 4.3.3
## Warning: package 'rsample' was built under R version 4.3.3
## Warning: package 'tune' was built under R version 4.3.3
## Warning: package 'workflows' was built under R version 4.3.3
## Warning: package 'workflowsets' was built under R version 4.3.3
## Warning: package 'yardstick' was built under R version 4.3.3
## ── 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/
set.seed(1234)

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

data_cv <- rsample::vfold_cv(data_train, strata = ceo_dismissal)
data_cv
## #  10-fold cross-validation using stratification 
## # A tibble: 10 × 2
##    splits             id    
##    <list>             <chr> 
##  1 <split [5044/562]> Fold01
##  2 <split [5044/562]> Fold02
##  3 <split [5045/561]> Fold03
##  4 <split [5045/561]> Fold04
##  5 <split [5046/560]> Fold05
##  6 <split [5046/560]> Fold06
##  7 <split [5046/560]> Fold07
##  8 <split [5046/560]> Fold08
##  9 <split [5046/560]> Fold09
## 10 <split [5046/560]> Fold10

Preprocess Data

library(themis)
## Warning: package 'themis' was built under R version 4.3.3
# data_train <- data_train %>% mutate(ceo_dismissal = as.factor(ceo_dismissal))

xgboost_rec <- recipes::recipe(ceo_dismissal ~ ., data = data_train) %>% 
    update_role(dismissal_dataset_id, new_role = "ID") %>%
    step_dummy(all_nominal_predictors()) %>%
    step_smote(ceo_dismissal) 

xgboost_rec %>%
    prep() %>%
    juice() %>%
    glimpse()
## Rows: 8,988
## Columns: 15
## $ dismissal_dataset_id <fct> 12, 31, 43, 51, 63, 75, 76, 80, 99, 109, 110, 112…
## $ fyear                <dbl> 1997, 1998, 2001, 1997, 1997, 1993, 2007, 1993, 1…
## $ co_per_rol           <dbl> 1, 6, 11, 16, 22, 33, 34, 43, 60, 66, 68, 71, 73,…
## $ tenure_no_ceodb      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ max_tenure_ceodb     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ fyear_gone           <dbl> 1998, 1998, 2002, 1997, 1998, 1995, 2007, 1993, 2…
## $ ceo_dismissal        <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ departure_code_X2    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ departure_code_X3    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ departure_code_X4    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ departure_code_X5    <dbl> 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0…
## $ departure_code_X6    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
## $ departure_code_X7    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1…
## $ departure_code_X8    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ departure_code_X9    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…

Specify Model

xgboost_spec <- 
  boost_tree(trees = tune()) %>% 
  set_mode("classification") %>% 
  set_engine("xgboost") 

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

Tune Hyperparameters

doParallel::registerDoParallel()

set.seed(65743)
xgboost_tune <-
  tune_grid(xgboost_workflow,
            resamples = data_cv, 
            grid = 5)