Predicting When Kickers Get Iced with {tidymodels}

Tidymodels

My interpretation of the building blocks are recipes which handle how data should pre-processed, specifications which tells tidymodels which algorithms and parameters to use, and workflows that bring them together. Since I’ve done most of the pre-processing in the data gathering piece these recipes will be pretty vanilla. However, this data is heavily imbalanced with only 4% of attempts being iced. So I will have two recipes. The first sets up the formula and one-hot encodes the categorical variables.

rec_norm <- recipe(is_iced ~ ., data = ice_train) %>%
  step_dummy(all_nominal_predictors(), one_hot =T) 

and the second will add a second step that uses step_smote() to create new examples of the minority class to fix the class imbalance problem.

rec_smote <- recipe(is_iced ~ ., data = ice_train) %>%
  step_dummy(all_nominal_predictors(), one_hot = T) %>%
  step_smote(is_iced) 

Then I’ll define my specification. The hyper-parameters that I want to tune are set to tune() and then I tell {tidymodels} that I want to use XGBoost for a classification problem.

xg_spec <- boost_tree(
  trees = tune(), 
  tree_depth = tune(), 
  min_n = tune(), 
  loss_reduction = tune(),                     
  sample_size = tune(), 
  mtry = tune(),         
  learn_rate = tune(), 
  stop_iter = tune()
) %>% 
  set_engine("xgboost") %>% 
  set_mode("classification")

The recipes and specifications are combined in workflows (if using 1 recipe and 1 specification) or workflow sets if wanting to use different combinations. In the workflow_set() function you can specify a list of recipes as preproc and a list of specifications as models. The cross parameter being set to true creates every possible combination. For this analysis I’ll have 2 preproc/model combinations:

wf_sets <- workflow_set(
  preproc = list(norm = rec_norm, 
                 smote = rec_smote),
  models = list(vanilla = xg_spec),
  cross = T
)

Next step is setting up the grid of parameters that will be tried in the model specifications above. Since there a lot of parameters to be tuned and I don’t want this to run forever I’m using grid_latin_hypercube to set 100 combinations of parameters that try to cover the entire parameter space but without running every combination.

grid <- grid_latin_hypercube(
  trees(),
  tree_depth(),
  min_n(),
  loss_reduction(),
  sample_size = sample_prop(),
  finalize(mtry(), ice_train),
  learn_rate(),
  stop_iter(range = c(10L, 50L)),
  size = 100
)

Finally its time to train the various workflows that have been designed. To do this I’ll pass the workflow_set defined above into the workflow_map() function. The “tune_race_anova” specification tells the training process to abandon certain hyper-parameter values if they’re not showing value. More detail can be found Julia Silge’s post. Also passed into this function are the resamples generated from the 5 folds, the grid of parameters, a control set that will save the predictions and workflows so that I can revisit them later on. Finally, I create a metric set of the performance metrics I want to calculate here choosing accuracy , ROC AUC, Multinomial Log Loss, and F1 Measure.

#Set up Multiple Cores
doParallel::registerDoParallel(cores = 4)


tuned_results <- wf_sets %>% 
  workflow_map(
    "tune_race_anova",
    resamples = train_5fold,
    grid = grid,
    control = control_race(save_pred = TRUE,
                           parallel_over = "everything",
                           save_workflow = TRUE),
    metrics = metric_set(f_meas, accuracy, roc_auc, mn_log_loss, pr_auc, precision, recall),
    seed = 20210109
  )

Tidymodels has an autoplot() function which will plot the best scoring model runs for each metric. However, I want a little more customization then what that function (or at least what I know of that function) provides. Using map_dfr() I’m going to stack the top model for each specification for each of the 5 performance metrics on top of each other using rank_results() to get the top model for each config for each metric.

perf_stats <- map_dfr(c('accuracy', 'roc_auc', 'mn_log_loss', 'pr_auc', 'f_meas',
                        'precision', 'recall'),
                ~rank_results(tuned_results, rank_metric = .x, select_best = T) %>% 
                filter(.metric == .x) 
        )

perf_stats %>% 
  ggplot(aes(x = wflow_id, color = wflow_id, y = mean)) +
    geom_pointrange(aes(y = mean, ymin = mean - 1.96*std_err, ymax = mean + 1.96*std_err)) + 
    facet_wrap(~.metric, scales = "free_y") + 
    scale_color_discrete(guide = 'none') + 
    labs(title = "Performance Metric for Tuned Results",
         x = "Model Spec",
         y = "Metric Value",
         color = "Model Config"
    ) + 
    theme_light()

Since I do care about the positive class more than the negative class but I don’t have a strong preference to false positive vs. false negative being more costly I’m going to use F1-Score as the performance metric I care most about. As expected the plain vanilla specification had a higher accuracy than the version using SMOTE to correct for imbalance. But it had lower values for F1, PR AUC and ROC AUC. I can also use rank_results() to show the top models for the F1 measure across the different specifications:

rank_results(tuned_results, rank_metric = 'f_meas') %>%
  select(wflow_id, .config, .metric, mean, std_err) %>%
  filter(.metric == 'f_meas') %>% 
  kable()

The top 7 models by F1 are all various configurations of the SMOTE recipe. The best model specification had an average F1 of 0.411 across the five folds. To get a better understanding of what this model’s specification actually was I extract the model configuration that has the best F1-score by using extract_workflow_set_result() with the workflow id and then select_best() with the metric I care about:

##Get Best Model
best_set <- tuned_results %>% 
  extract_workflow_set_result('smote_vanilla') %>% 
  select_best(metric = 'f_meas')

kable(best_set)

The best model in this case had 5 random predictors, 1641 trees, and so on.

Now that I know which model configuration is the best one, the last step is to final the model using the full training data and predict on the test set. The next block of code extracts the workflow, sets the parameters to be those from the best_set defined above using finalize_workflow, and then last_fit() does the final fitting using the full training set and prediction on the testing data when we pass it the workflow and the split object.

final_fit <- tuned_results %>%
  extract_workflow('smote_vanilla') %>%
  finalize_workflow(best_set) %>%
  last_fit(ice_split, metrics=metric_set(accuracy, roc_auc, mn_log_loss, 
                                         pr_auc, f_meas, precision, recall))

Then with collect_metrics() I can see the final results when the model was applied to the test set that had been unused thus far.

collect_metrics(final_fit) %>% 
  kable()

The F1 score is actually higher than in the training at 0.43% with a precision of 34%, a recall of 59%, and a ROC AUC of 0.91%.

Tidymodels also makes it very easy to display ROC curves using collect_predictions to get the predictions from the final model and test set and roc_curve to calculate the sensitivity and specificity.

collect_predictions(final_fit) %>%
  roc_curve(is_iced, .pred_iced) %>%
  ggplot(aes(1 - specificity, sensitivity)) +
  geom_abline(lty = 2, color = "gray80", size = 1.5) +
  geom_path(alpha = 0.8, size = 1) +
  coord_equal() +
  labs(color = NULL)

As well as calculate the confusion matrix with collect_predictions and conf_mat.

collect_predictions(final_fit) %>%
  conf_mat(is_iced, .pred_class) %>%
  autoplot(type = 'heatmap')

Variable Importance

Variable Importance plots are one way of understanding which predictor has the largest effect on the model outcomes. There are many ways to measure variable importance but the one I’m using is the default in the {vip} package for XGBoost which is “gain”. Variable importance using gain measures the fractional contribution of each feature to the model based on the total gain of the feature’s splits where gain is the improvement to accuracy brought by a feature to its branches.

The {vip} package provides variable importance when given a model object as an input. To get that I use extract_fit_parsnip() to get the parsnip version of the model object. Then the vip() function does the rest.

library(vip)

extract_fit_parsnip(final_fit) %>%
  vip(geom = "point", include_type = T) + 
  geom_text(aes(label = scales::percent(Importance, accuracy = 1)),
            nudge_y = 0.023) + 
  theme_light()

Unsurprisingly, the regulation time remaining is the most important feature which makes sense because the amount of time remaining dictates whether using the timeout on a kick is worthwhile. Although whether the kicking team is the home team being the 2nd most important feature is a bit more surprising as I would have thought game situation would apply more than home or away team. I thought score difference would be higher.

Partial Dependency

Variable importance tells us “what variables matters” but it doesn’t tell us “how they matter”. Are the relationships between the predictors and predictions linear or non-linear. Is there some magic number where a step function occurs. Variable Importance cannot answer these questions, but partial dependency plots can!

A partial dependency plot shows the effect of a predictor on the model outcome holding everything else constant. The {pdp} package can be used to generate these plots. The {pdp} package is a little less friendly with {tidymodels} since you need to provide the native model object rather than the parsnip version (which is still easily accessible using extract_fit_engine()). Also, the data passed into the partial() function needs to be the same as the data that actually goes into the model object. So I create fitted_data by prep()ing the recipe and then bake()’ing which applies the recipe to the original data set.

The partial function can also take a while to run, so I’m using {furrr} which allows for {purrr} functions to be run in parallel on the {future} backend. In the future_map_dfr function, I’m running partial on every predictor in the data and stacking the results on top of each other so that I can plot them in the final step. The use of prob=T converts the model output to a probability but since XGBoost probabilities are uncalibrated best not to read too much into the values.

library(pdp)

##Get Processed Training Data
model_object <- extract_fit_engine(final_fit)

fitted_data <- rec_smote %>%
  prep() %>%
  bake(new_data = model_data) %>%
  select(-is_iced)

library(furrr)
plan(multisession, workers = 4)

all_partial <- future_map_dfr(
  names(fitted_data), ~as_tibble(partial(
    model_object,
    train = fitted_data,
    pred.var = .x,
    type = 'classification',
    plot = F,
    prob = T, #Converts model output to probability scale
    trim.outliers = T
  )) %>% 
    mutate(var = .x) %>%
    rename(value = all_of(.x)),
  .progress = T,
  .options = furrr_options(seed = 20220109)
)

all_partial %>% 
  #Remove Prior Miss since its one-hot encoded
  filter(!str_detect(var, 'prior_miss|overtime')) %>% 
  ggplot(aes(x = value, y = yhat, color = var)) + 
    geom_line() + 
    geom_smooth(se = F, lty = 2, span = .5) + 
    facet_wrap(~var, scales = "free") + 
    #scale_y_continuous(labels = percent_format(accuracy = .1)) + 
    scale_color_discrete(guide = 'none') +
    labs(title = "Partial Dependency Plots for Whether A Kick Gets Iced?",
         subtitle = "Looking at 19,072 NCAA Field Goal Attempts (2013-2021)",
         x = "Variable Value",
         y = "Prob. of Attempt Getting Iced") + 
    theme_light()

From these plots we can tell that the likelihood of getting iced increases when:

• The Attempted Distance is between 30-50 yards
• When two teams are expected to be somewhat evenly matched (based on pre-game win probabilities)
• When nearly the end of the game or the end of the half (that middle spike in regulation time remaining is halftime since timeouts reset at the beginning of each half)
• When the kicking team is losing by a very small margin (or when the game is within +/- 10 points)

While variable importance told us that Regulation Time Remaining was the most important variable, the partial dependency plot shows us how it affects the model in a non-linear way.

SHAP Values

The next measure of interpretability combines pieces of both variable importance and partial dependency plots. SHAP values are claimed to be the most advanced method to interpret results from tree-based models. They are based on Shaply values from game theory and measure feature importance based on the marginal contribution of each predictor for each observation to the model output.

The {SHAPforxgboost} package provides an interface to getting SHAP values. The plot that will give us overall variable importance is the SHAP summary plot which we’ll get using shap.plot.summary. However, first the data structure needs to be prepped using the model object and the training data in a matrix.

library(SHAPforxgboost)

shap_long <- shap.prep(xgb_model = extract_fit_engine(final_fit), 
                        X_train = fitted_data %>% as.matrix())
                       
shap.plot.summary(shap_long)

In the summary plot, the most important variables are ordered from top to bottom. Within any given variable each point represents an observation. The shading of the point represents which that observation has a higher or lower value for that features. For example, in regulation time remaining lower amounts of remaining time will be orange while higher amounts will be purple. The position on the left or right side of zero represents whether they decrease or increase the likelihood of getting iced. For regulation time remaining notice that the very purple is strongly negative (on the left side) and the very orange is strongly positive (on the right side).

Similar to the variable importance plot, regulation time remaining was the most important feature.

We can also get dependency plots similar to the partial dependency plots with SHAP values using shap.plot.dependence. We’ll look at the regulation time remaining on the x-axis and the SHAP values for regulation time remaining on the y-axis. Since this returns a ggplot object, I’ll add in vertical lines to represent the end of each quarter.

SHAPforxgboost::shap.plot.dependence(data_long = shap_long, x = 'regulation_time_remaining', 
                                     y = 'regulation_time_remaining', 
                                     color_feature = 'regulation_time_remaining') + 
  ggtitle("Shap Values vs. Regulation Time Remaining") + 
  geom_vline(xintercept = 0, lty = 2) + 
    geom_vline(xintercept = 900, lty = 2) + 
    geom_vline(xintercept = 1800, lty = 2) + 
    geom_vline(xintercept = 2700, lty = 2) 

Similar to the summary plot, the less time remaining in the game the more orange the point and the more time remaining the more purple. Again, like in the partial dependency plot, we see a non-linear relationship with increases towards the end of each quarter and heavy spikes in the last 3 minutes of the 2nd and 4th quarters.

This is just an example of the things that can be done with SHAP values but hopefully its usefulness for understanding both what’s important and how its important has been illustrated.