Predicting NFL stadium attendances with tidymodels

Introduction

In this post I’m going to be trying out the new tidymodels framework in R.

I’ve been reading through the corresponding book and it looks to me like a real game-changer for building robust statistical/ML models quickly in the R language. But reading only gets you so far! Time to test it out.

The data I’m going to be using is a #TidyTuesday dataset from last year related to crowd attendances in the American national football league (NFL).

The tidymodels framework covers all stages of the modeling lifecyle, from data preprocessing and resampling to model building, tuning and evaluation. In this post I’m only going to look at data resampling and model building and evaluation. I’ll look at the other stages in future posts.

I owe credit to Julia Silge, whose blog post served as the inspiration for this piece.

Load and inspect the data

attendance <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-02-04/attendance.csv")
results <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-02-04/standings.csv")

Our first dataset is about attendance per team every week since the year 2000:

head(attendance)

## # A tibble: 6 x 8
##   team    team_name  year  total   home   away  week weekly_attendance
##   <chr>   <chr>     <dbl>  <dbl>  <dbl>  <dbl> <dbl>             <dbl>
## 1 Arizona Cardinals  2000 893926 387475 506451     1             77434
## 2 Arizona Cardinals  2000 893926 387475 506451     2             66009
## 3 Arizona Cardinals  2000 893926 387475 506451     3                NA
## 4 Arizona Cardinals  2000 893926 387475 506451     4             71801
## 5 Arizona Cardinals  2000 893926 387475 506451     5             66985
## 6 Arizona Cardinals  2000 893926 387475 506451     6             44296

And our second dataset gives us results stats for each team about things like final standings, wins/losses, points for/against and whether the team made the playoffs for every year since 2000.

head(results) 

## # A tibble: 6 x 15
##   team    team_name  year  wins  loss points_for points_against points_differen…
##   <chr>   <chr>     <dbl> <dbl> <dbl>      <dbl>          <dbl>            <dbl>
## 1 Miami   Dolphins   2000    11     5        323            226               97
## 2 Indian… Colts      2000    10     6        429            326              103
## 3 New Yo… Jets       2000     9     7        321            321                0
## 4 Buffalo Bills      2000     8     8        315            350              -35
## 5 New En… Patriots   2000     5    11        276            338              -62
## 6 Tennes… Titans     2000    13     3        346            191              155
## # … with 7 more variables: margin_of_victory <dbl>, strength_of_schedule <dbl>,
## #   simple_rating <dbl>, offensive_ranking <dbl>, defensive_ranking <dbl>,
## #   playoffs <chr>, sb_winner <chr>

Let’s join these datasets by team and year to get one large dataset combining data on results and attendances which we can use for modelling.

attendance_joined <- attendance %>%
  left_join(results,
    by = c("year", "team_name", "team")
  )

head(attendance_joined)

## # A tibble: 6 x 20
##   team   team_name  year  total   home   away  week weekly_attendan…  wins  loss
##   <chr>  <chr>     <dbl>  <dbl>  <dbl>  <dbl> <dbl>            <dbl> <dbl> <dbl>
## 1 Arizo… Cardinals  2000 893926 387475 506451     1            77434     3    13
## 2 Arizo… Cardinals  2000 893926 387475 506451     2            66009     3    13
## 3 Arizo… Cardinals  2000 893926 387475 506451     3               NA     3    13
## 4 Arizo… Cardinals  2000 893926 387475 506451     4            71801     3    13
## 5 Arizo… Cardinals  2000 893926 387475 506451     5            66985     3    13
## 6 Arizo… Cardinals  2000 893926 387475 506451     6            44296     3    13
## # … with 10 more variables: points_for <dbl>, points_against <dbl>,
## #   points_differential <dbl>, margin_of_victory <dbl>,
## #   strength_of_schedule <dbl>, simple_rating <dbl>, offensive_ranking <dbl>,
## #   defensive_ranking <dbl>, playoffs <chr>, sb_winner <chr>

Data exploration

As usual, we need to look at our data first to get a sense of it before starting any modelling. I like to do this using the skim function from R’s skimr package, which gives us an overview of the data:

skimr::skim(attendance_joined)

Variable type: character

Variable type: numeric

This looks like a pretty kind dataset, the only obvious problem being the 638 missing values which we will need to deal with. None of the features look particularly skewed (scroll to the right on the table above to see this).

Let’s see how weekly attendance varied over the 20 years of data we have by utilising box and whisker plots. Has going to NFL matches become more popular in recent years?

attendance_joined %>%
  mutate(year = factor(year)) %>%
  ggplot(aes(year, weekly_attendance, fill = year)) +
  geom_boxplot(show.legend = FALSE, outlier.alpha = 0.5) +
  labs(
    x = "Year of NFL season",
    y = "Weekly NFL game attendance"
  ) +
  scale_y_continuous(labels = comma)

Not much variation really. The 2020 figures would be interesting - presumably be much lower due to the pandemic - but we don’t have them unfortunately.

I wonder if teams that got to the playoffs had a higher weekly attendance? Let’s use a histogram to check out the respective distributions.

attendance_joined %>%
  filter(!is.na(weekly_attendance)) %>% 
  group_by(team_name, year, margin_of_victory, playoffs) %>% 
  summarise(avg_weekly_attendance = as.integer(mean(weekly_attendance))) %>% 
  ggplot(aes(avg_weekly_attendance, fill = playoffs)) +
  geom_histogram(position = "identity", alpha = 0.4) +
  labs(
    x = "Average weekly attendance",
    y = "Number of teams",
    fill = NULL
  )

It looks like the Playoffs distribution is slightly to the right of the No Playoffs. Let’s see if we can view this difference more explicitly, using a boxplot.

attendance_joined %>%
  filter(!is.na(weekly_attendance)) %>% 
  group_by(team_name, year, margin_of_victory, playoffs) %>% 
  summarise(avg_weekly_attendance = as.integer(mean(weekly_attendance)))  %>% 
  ggplot(aes(playoffs, avg_weekly_attendance, fill = playoffs)) +
  geom_boxplot() +
  labs(
    x = "Playoff reached",
    y = "Average weekly attendance"
  )

The lower quartile and median are definitely higher for the teams that made it to the Playoffs. We’ll need to take this into account later when we modelling.

In a more comprehensive analysis I would do further EDA and then systematic feature selection, but that’s not the purpose of this exercise - I just want to practise building a model to predict weekly_attendance using several features of my choosing.

Let’s prepare the data for modelling by removing any rows without a weekly attendance and keeping only columns that I think may be predictors for attendance:

attendance_df <- attendance_joined %>%
  filter(!is.na(weekly_attendance)) %>%
  select(
    weekly_attendance, team_name, year, week,
    margin_of_victory, strength_of_schedule, playoffs
  )

head(attendance_df)

## # A tibble: 6 x 7
##   weekly_attendance team_name  year  week margin_of_victory strength_of_schedule
##               <dbl> <chr>     <dbl> <dbl>             <dbl>                <dbl>
## 1             77434 Cardinals  2000     1             -14.6                 -0.7
## 2             66009 Cardinals  2000     2             -14.6                 -0.7
## 3             71801 Cardinals  2000     4             -14.6                 -0.7
## 4             66985 Cardinals  2000     5             -14.6                 -0.7
## 5             44296 Cardinals  2000     6             -14.6                 -0.7
## 6             38293 Cardinals  2000     7             -14.6                 -0.7
## # … with 1 more variable: playoffs <chr>

Building models

library(tidymodels)

set.seed(1234)
attendance_split <- attendance_df %>%
  initial_split(strata = playoffs)

nfl_train <- training(attendance_split)
nfl_test <- testing(attendance_split)

lm_spec <- linear_reg() %>% 
  set_engine(engine = "lm")

lm_spec

## Linear Regression Model Specification (regression)
## 
## Computational engine: lm

lm_fit <- lm_spec %>% 
  fit(weekly_attendance ~ .,
      data = nfl_train)

lm_fit

## parsnip model object
## 
## Fit time:  26ms 
## 
## Call:
## stats::lm(formula = weekly_attendance ~ ., data = data)
## 
## Coefficients:
##          (Intercept)        team_nameBears      team_nameBengals  
##             -81107.9               -2879.8               -4875.5  
##       team_nameBills      team_nameBroncos       team_nameBrowns  
##               -361.1                2805.2                -324.1  
##  team_nameBuccaneers    team_nameCardinals     team_nameChargers  
##              -3063.7               -6139.8               -6489.3  
##      team_nameChiefs        team_nameColts      team_nameCowboys  
##               1974.3               -3392.8                6068.7  
##    team_nameDolphins       team_nameEagles      team_nameFalcons  
##                139.7                1259.2                -204.2  
##      team_nameGiants      team_nameJaguars         team_nameJets  
##               5447.1               -3095.5                4044.2  
##       team_nameLions      team_namePackers     team_namePanthers  
##              -3480.7                1114.1                1227.3  
##    team_namePatriots      team_nameRaiders         team_nameRams  
##               -214.2               -6324.7               -2884.8  
##      team_nameRavens     team_nameRedskins       team_nameSaints  
##               -398.9                6447.1                 381.0  
##    team_nameSeahawks     team_nameSteelers       team_nameTexans  
##              -1405.9               -3567.8                 264.1  
##      team_nameTitans      team_nameVikings                  year  
##              -1118.2               -3183.1                  74.7  
##                 week     margin_of_victory  strength_of_schedule  
##                -72.8                 137.6                 230.7  
##     playoffsPlayoffs  
##               -427.9

rf_spec <- rand_forest(mode = "regression") %>%
  set_engine("ranger")

rf_spec

## Random Forest Model Specification (regression)
## 
## Computational engine: ranger

rf_fit <- rf_spec %>%
  fit(weekly_attendance ~ .,
    data = nfl_train
  )

rf_fit

## parsnip model object
## 
## Fit time:  7.4s 
## Ranger result
## 
## Call:
##  ranger::ranger(x = maybe_data_frame(x), y = y, num.threads = 1,      verbose = FALSE, seed = sample.int(10^5, 1)) 
## 
## Type:                             Regression 
## Number of trees:                  500 
## Sample size:                      7656 
## Number of independent variables:  6 
## Mtry:                             2 
## Target node size:                 5 
## Variable importance mode:         none 
## Splitrule:                        variance 
## OOB prediction error (MSE):       7.5e+07 
## R squared (OOB):                  0.082

Evaluating the models

Now we need to measure how effective our models were and compare them. One way to do this is to predict() the weekly attendance for the test set.

results_test <- lm_fit %>%
  predict(new_data = nfl_test) %>%
  mutate(
    truth = nfl_test$weekly_attendance,
    model = "lm"
  ) %>%
  bind_rows(rf_fit %>%
    predict(new_data = nfl_test) %>%
    mutate(
      truth = nfl_test$weekly_attendance,
      model = "rf"
    ))

This gives us a tibble with the predictions and true values for both models

results_test %>% slice(c(1:5, (n() -5):n()))

## # A tibble: 11 x 3
##     .pred truth model
##     <dbl> <dbl> <chr>
##  1 59162. 65356 lm   
##  2 59017. 50289 lm   
##  3 58871. 37452 lm   
##  4 58798. 65711 lm   
##  5 66880. 73025 lm   
##  6 67178. 58613 rf   
##  7 68914. 66910 rf   
##  8 65276. 60926 rf   
##  9 62456. 65265 rf   
## 10 66187. 71504 rf   
## 11 67666. 90646 rf

Let’s visualise the predictions versus the truth for both models:

results_test %>%
  ggplot(aes(x = truth, y = .pred, color = model)) +
  geom_abline(lty = 2) +
  geom_point(alpha = 0.5) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  facet_wrap(~model) +
  labs(y = "Predicted weekly attendance", x = "Actual weekly attendance") +
  coord_obs_pred() +  # scale and size the x- and y-axis uniformly
  theme(panel.spacing = unit(2, "lines"))  # widen the gap between plots for aesthetic purposes

It looks like neither model does a particularly good job of predicting the attendance, though it’s hard to see which is better. Let’s try to establish this numerically rather than visually. The yardstick package contains a metric_set function, from which we can compute any metrics of interest for both models.

As discussed here in the tidymodels textbook, the choice of which metric to examine can be critical. Accuracy, as measured by the RMSE, and correlation, as measured by the R², are not the same thing, and in proper ML projects the choice of which to optimise should be made on a case-by-case basis. For now, we will simply evaluate both.

nfl_metrics <- metric_set(rmse, rsq)

results_test %>% 
  group_by(model) %>% 
  nfl_metrics(truth = truth, estimate = .pred) %>% 
  select(model, .metric, .estimate)

## # A tibble: 4 x 3
##   model .metric .estimate
##   <chr> <chr>       <dbl>
## 1 lm    rmse     8351.   
## 2 rf    rmse     8582.   
## 3 lm    rsq         0.142
## 4 rf    rsq         0.117

The linear model produced a lower error and higher r-squared score, showing that for this particular dataset it’s the better model. Its average error in predicting weekly attendance is ~8350 people, which isn’t too bad given how basic it is!

Presumably the random forest model, being much more complex than the linear one, overfit to the training data, which is why it did a worse job on the test set (further analysis showed this to indeed be the case - code not shown).

Resampling the data

As touched on earlier, in a real-lfie ML project, we would not simply train models on the training set and then immediately apply them to the test set. The test set should not be touched until validation for model selection has been completed. Instead, we need to take a resampling approach using cross-validation.

Below, we use the vfold_cv() function from the rsample package to achieve this and then re-compute our metrics of interest for both models.

Linear model:

set.seed(1234)
nfl_folds <- vfold_cv(nfl_train, strata = playoffs)

lm_res <- fit_resamples(
  lm_spec,
  weekly_attendance ~ .,
  nfl_folds,
  control = control_resamples(save_pred = TRUE)
)

lm_res %>%
  collect_metrics() %>% 
  select(.metric, mean)

## # A tibble: 2 x 2
##   .metric     mean
##   <chr>      <dbl>
## 1 rmse    8335.   
## 2 rsq        0.148

Random forest model:

set.seed(1234)
nfl_folds <- vfold_cv(nfl_train, strata = playoffs)

rf_res <- fit_resamples(
  rf_spec,
  weekly_attendance ~ .,
  nfl_folds,
  control = control_resamples(save_pred = TRUE)
)

rf_res %>%
  collect_metrics() %>% 
  select(.metric, mean)

## # A tibble: 2 x 2
##   .metric     mean
##   <chr>      <dbl>
## 1 rmse    8617.   
## 2 rsq        0.117

So now using only the training set we have got metric estimates very close to the ones we computed from the test set earlier. Using cross-validation in this way allows us to be more economical with our data spending, as well as making it possible to tune the models without the risk of any data leakage from the test set.

Conclusions

In this post we:

1. Loaded and explored the data, which allowed us to make more informed choices about how to later model it.
2. Gained familiarity with the tidymodels framework of packages, such as:

• Using rsample for data splitting
• Using parsnip for model specifying and fitting.
• Using yardstick for model evaluation
• Using rsample again, this time for doing cross-validation.

3. Established that, for this particular dataset, the linear model did a better job of predicting on new data than the random forest, which overfit to the training set.