Feature Engineering and Preprocessing

The first feature I would like add is the club size for the participants who signed up as part of a club, as it seems like there is some difference between finshers and DNFs:

#engineer club size
results_preprocessed <- results %>%
  #filter out runners which droped out before KM21
  filter(!(dnf=="DNF" & is.na(hm))) %>% 
  #drop not-used features or features with data leakage (eg. start delay is only available for finishers) 
  select(-overall_place, -birth_year, -gender_place, -age_group_place, -brutto_time, -netto_time, -city, -start_delay, -age_bins,
         -x25km, -x30km, -x35km, -x40km, -split_25km, -split_30km, -split_35km, -split_40km, -split_finish, -age_group) %>%
  #group by year and club to geth club size per year
  group_by(year, club) %>%
  #set club size to 0 for people without a club
  mutate(club_starters = case_when(!is.na(club) ~ as.numeric(n()),
                                   is.na(club) ~ 0)) %>% 
  ungroup() %>% 
  select(-club)

#plot regarding clubsize
club_size_p <- results_preprocessed %>%
  ggplot(aes(x=dnf, y=club_starters, fill=dnf))+
  geom_boxplot()+
  scale_y_continuous(trans='log10')+
  theme_light()+
  labs(y="Number of starters from the same club",
       x=NULL)+
  theme(plot.background = element_rect(fill = "#F4F4F4"),
        panel.background = element_rect(fill = "#F4F4F4",
                                colour = "#F4F4F4",
                                size = 0.5, linetype = "solid"))
#plot
club_size_p

As the next feature the standard variation across the relative split times is calculated for each runner:

results_preprocessed <- results_preprocessed %>%
  mutate(pace_variation = pmap_dbl(.l = select(., starts_with("split_")), .f = ~sd(c(...), na.rm = TRUE)))

results_preprocessed %>% 
  ggplot(aes(x=dnf, y=pace_variation, fill=dnf))+
  geom_boxplot()+
  theme_light()+
  labs(y="Relative split times standard deviation",
       x=NULL)+
  theme(plot.background = element_rect(fill = "#F4F4F4"),
        panel.background = element_rect(fill = "#F4F4F4",
                                colour = "#F4F4F4",
                                size = 0.5, linetype = "solid"))

In the final feature engineering step I focus on the relative and absolute split times. First I compare the pace of the relative split times to the average pace during the first half of the marathon. On these new features rel_pace... I’m conducting a linear regression in order to get the slope across the new feature (–> did a runner keep a steady pace or did the pace increas/decrease during the first half of the marathon).
A similar approach is applied to the perecentile rank on the relative and absolut split times in order to see if a runner got faster or slower compared to the other runners.
As the last preprocessing step the distance to the start from a participant’s hometown is discretized into for roughly equal groups and a “missing” category added.

results_preprocessed <- results_preprocessed %>%
  #drop observations with missing split times in order to be able to do linear regression to get relative pace slope
  drop_na(starts_with("split_")) %>%
  #get the relative pace compared to the average pace of the first half of the marathon
  mutate(rel_pace_5km = (split_5km/5)/(hm/21.0975),
         rel_pace_10km = (split_10km/5)/(hm/21.0975),
         rel_pace_15km = (split_15km/5)/(hm/21.0975),
         rel_pace_hm = (split_hm/6.0975)/(hm/21.0975)) %>%
  #calculate the slope of the relative pace
  mutate(rel_pace_slope = pmap_dbl(.l = select(., starts_with("rel_pace_")), .f = ~lm(rel_pace~split,
                                                                                      tibble(rel_pace = c(...),
                                                                                             split = c(1:4)))$coeff[[2]])) %>%

  #calculate the percentile per year on absolute and relative split times
  group_by(year) %>%
  mutate(rel_percentile_5km = percent_rank(split_5km),
         rel_percentile_10km = percent_rank(split_10km),
         rel_percentile_15km = percent_rank(split_15km),
         rel_percentile_hm = percent_rank(split_hm),
         percentile_5km = percent_rank(split_5km),
         percentile_10km = percent_rank(x10km),
         percentile_15km = percent_rank(x15km),
         percentile_hm = percent_rank(hm)) %>%
  ungroup() %>%
  select(-year) %>%
  #calculate the slope of the percentiles
  mutate(rel_percentile_slope = pmap_dbl(.l = select(., starts_with("rel_percentile_")), .f = ~lm(rel_percentile~split,
                                                                                      tibble(rel_percentile = c(...),
                                                                                             split = c(1:4)))$coeff[[2]]),
         percentile_slope = pmap_dbl(.l = select(., starts_with("percentile_")), .f = ~lm(percentile~split,
                                                                                      tibble(percentile = c(...),
                                                                                             split = c(1:4)))$coeff[[2]]))
#discretize distance into 4 equal sized groups and add "missing" as an explicit group
results_preprocessed <- results_preprocessed %>% 
  mutate(dist_cat = cut_number(dist, n=4) %>% fct_explicit_na()) %>%
  #drop no longer used variables
  select(-dist, -x5km, -x10km, -x15km, -hm, -starts_with("split_"))

Modelling!!

Splitting the data, creating recipe (with SMOTE via the {themis} package due to heavy class imbalance) and create folds for resampling:

library(tidymodels)
library(themis) #for SMOTE

#create initial split with strata
set.seed(1234)
marathon_split <- initial_split(results_preprocessed, prop = .8, strata = dnf)

marathon_train <- training(marathon_split)
marathon_test <- testing(marathon_split)

#create folds for resampling
set.seed(1234)
marathon_folds <- vfold_cv(marathon_train, strata = dnf, v = 5)

#create recipe
marathon_rec <- recipe(dnf ~ ., data = marathon_train) %>%
  step_rm(runner_id) %>% 
  step_dummy(all_nominal(), -all_outcomes()) %>%
  step_smote(dnf) %>%
  step_zv(all_numeric()) %>%
  step_normalize(all_numeric()) %>% 
  prep()

We will create two classification models: One using logistic regression and the other using eXtreme Gradient Boosting (xgboost).
Both models will be tuned on various tuning parameters via a grid search using the {tune} package.

library(workflows)
library(tune)

#set model specification
glm_spec <- logistic_reg(mode = "classification",
                         penalty = tune(),
                         mixture = tune()) %>%
  set_engine("glmnet")

#workflow
glm_wflow <-
  workflow() %>%
  add_model(glm_spec) %>%
  add_recipe(marathon_rec)


#train resamples
set.seed(1291)
glm_res <- 
  tune_grid(
    glm_wflow, 
    resamples = marathon_folds,
    grid = 5,
    metrics = metric_set(roc_auc, sens, spec),
    control = control_resamples(save_pred = TRUE, verbose = TRUE, allow_par = TRUE)
  )

glm_res %>%
  show_best(metric = "roc_auc") %>%
  kable() %>%
  kable_styling(font_size = 8)

Across the resample the best hyperparameter combination for the logistig regression model reached an AUC of 0.73.
Let’s see if the xgboost performs better?

#set model specification
xgb_spec <- boost_tree(mode = "classification",
                       mtry = tune(),
                       trees = tune(),
                       min_n = tune(),
                       tree_depth = tune(),
                       learn_rate = tune(),
                       loss_reduction = tune(),
                       sample_size = tune()) %>%
  set_engine("xgboost")

#workflow
xgb_wflow <-
  workflow() %>%
  add_model(xgb_spec) %>%
  add_recipe(marathon_rec)

#train resamples
#doParallel for paralell processing:
library(doParallel)
cl <- makePSOCKcluster(parallel::detectCores(logical = T)-1)
set.seed(1291)
xgb_res <- 
  tune_grid(
    xgb_wflow, 
    resamples = marathon_folds,
    grid = 20,
    metrics = metric_set(roc_auc, sens, spec),
    control = control_resamples(save_pred = TRUE, verbose = TRUE, allow_par = TRUE)
  )

xgb_res %>%
  show_best(metric = "roc_auc") %>%
  kable() %>%
  kable_styling(font_size = 8)

It seems like the AUC is a bit better than the one of the logistic regression. The ROC plot shows the difference between the two models across all hyperparameter combinations:

Across all hyperparameter combinations on all folds the logistic regression seems to perform better than the xgboost. Let’s fit the best combination of each model on the complete training data and get the predictions for the test data:

Both models perform quite similar on the test data.
Below you can find the confusion matrix for the logistic regression model:

##           Truth
## Prediction  DNF Finisher
##   DNF        47      525
##   Finisher   39     1306

Among the DNFs the model made most of false predicitions for runners within the 0.25 to 0.75 percentiles, while among the Finishers ir didn’t perform so good for the very low and very high ranked runners at the half-way point.

Overall the logistic regression model does quite a good job at identifying DNFs but also falsley a big portion of runners who actually finished.
Such a model could be further enhanced with booking data (eg. how many days before the marathon did a participant sign up?) and historical data (eg. when was the last time a participant ran the Köln Marathon?), but it already gives a great opportunity of identifing struggling participants and eventually give them a bit of extra support in form of energy gels, drinks or cheering!

Let me know how you think this model could be improved! 😄

/Lukas