Fit and predict with tidymodels for #TidyTuesday bird baths in Australia

This is the latest in my series of screencasts demonstrating how to use the tidymodels packages, from just getting started to tuning more complex models. Today’s screencast is good for folks who are newer to modeling or tidymodels; it focuses on how to use feature engineering together with a model algorithm and how to fit and predict, with this week’s #TidyTuesday dataset on bird baths in Australia. 🐦

Here is the code I used in the video, for those who prefer reading instead of or in addition to video.

Explore data

Our modeling goal is to predict whether we’ll see a bird at a bird bath in Australia, given info like what kind of bird we’re looking for and whether the bird bath is in an urban or rural location.

library(tidyverse)

bird_baths <- readr::read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-08-31/bird_baths.csv")

bird_baths %>%
  count(urban_rural)

## # A tibble: 3 × 2
##   urban_rural      n
##   <chr>        <int>
## 1 Rural        49686
## 2 Urban       111202
## 3 <NA>           169

Notice that there are some summary rows in the dataset with NA values for urban_rural, survey_year, etc. We can use that to choose some top bird types to focus on, instead of all the many bird types included in this dataset.

top_birds <-
  bird_baths %>%
  filter(is.na(urban_rural)) %>%
  arrange(-bird_count) %>%
  slice_max(bird_count, n = 15) %>%
  pull(bird_type)

top_birds

##  [1] "Noisy Miner"        "Australian Magpie"  "Rainbow Lorikeet"  
##  [4] "Red Wattlebird"     "Superb Fairy-wren"  "Magpie-lark"       
##  [7] "Pied Currawong"     "Crimson Rosella"    "Eastern Spinebill" 
## [10] "Spotted Dove"       "Lewin's Honeyeater" "Satin Bowerbird"   
## [13] "Crested Pigeon"     "Grey Fantail"       "Red-browed Finch"

How likely were the citizen scientists who collected this data to see birds of different types, in different locations?

bird_parsed <-
  bird_baths %>%
  filter(
    !is.na(urban_rural),
    bird_type %in% top_birds
  ) %>%
  group_by(urban_rural, bird_type) %>%
  summarise(bird_count = mean(bird_count), .groups = "drop")

p1 <-
  bird_parsed %>%
  ggplot(aes(bird_count, bird_type)) +
  geom_segment(
    data = bird_parsed %>%
      pivot_wider(
        names_from = urban_rural,
        values_from = bird_count
      ),
    aes(x = Rural, xend = Urban, y = bird_type, yend = bird_type),
    alpha = 0.7, color = "gray70", size = 1.5
  ) +
  geom_point(aes(color = urban_rural), size = 3) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = "Probability of seeing bird", y = NULL, color = NULL)

p1

Superb fairy-wrens are more rural, while noisy miners are more urban.

Let’s build a model to predict this probability of seeing a bird using just these two predictors.

bird_df <-
  bird_baths %>%
  filter(
    !is.na(urban_rural),
    bird_type %in% top_birds
  ) %>%
  mutate(bird_count = if_else(bird_count > 0, "bird", "no bird")) %>%
  mutate_if(is.character, as.factor)

Build a first model

Let’s start our modeling by setting up our “data budget.” We are going to use a simple logistic regression model that is unlikely to overfit, but let’s still split our data into training and testing, and then create resampling folds.

library(tidymodels)

set.seed(123)
bird_split <- initial_split(bird_df, strata = bird_count)
bird_train <- training(bird_split)
bird_test <- testing(bird_split)

set.seed(234)
bird_folds <- vfold_cv(bird_train, strata = bird_count)
bird_folds

## #  10-fold cross-validation using stratification 
## # A tibble: 10 × 2
##    splits              id    
##    <list>              <chr> 
##  1 <split [9637/1072]> Fold01
##  2 <split [9638/1071]> Fold02
##  3 <split [9638/1071]> Fold03
##  4 <split [9638/1071]> Fold04
##  5 <split [9638/1071]> Fold05
##  6 <split [9638/1071]> Fold06
##  7 <split [9638/1071]> Fold07
##  8 <split [9638/1071]> Fold08
##  9 <split [9639/1070]> Fold09
## 10 <split [9639/1070]> Fold10

We’ll make a couple of attempts at fitting models here, but they will all use straightforward logistic regression.

glm_spec <- logistic_reg()

For this first model, let’s set up our feature engineering recipe with our outcome and two predictors, and begin with only one preprocessing step to transform our nominal (factor or character, like urban_rural and bird_type) predictors to dummy or indicator variables. Then let’s put our preprocessing recipe together with our model specification in a workflow.

rec_basic <-
  recipe(bird_count ~ urban_rural + bird_type, data = bird_train) %>%
  step_dummy(all_nominal_predictors())

wf_basic <- workflow(rec_basic, glm_spec)

We could fit this one time to the training data, but to get better estimates of performance, let’s fit 10 times to our 10 resampling folds.

doParallel::registerDoParallel()
ctrl_preds <- control_resamples(save_pred = TRUE)
rs_basic <- fit_resamples(wf_basic, bird_folds, control = ctrl_preds)

How did this turn out? If we look at some overall metrics, accuracy does not look so bad:

collect_metrics(rs_basic)

## # A tibble: 2 × 6
##   .metric  .estimator  mean     n   std_err .config             
##   <chr>    <chr>      <dbl> <int>     <dbl> <chr>               
## 1 accuracy binary     0.822    10 0.0000762 Preprocessor1_Model1
## 2 roc_auc  binary     0.601    10 0.00783   Preprocessor1_Model1

This is because there were not many birds overall, though! The model is just saying “no bird” everywhere and getting good accuracy. The ROC curve, on the other hand, looks not so great.

augment(rs_basic) %>%
  roc_curve(bird_count, .pred_bird) %>%
  autoplot()

Add interactions

We know from the plot we made during EDA that there are interactions between whether a bird bath is urban/rural and what kinds of birds we see there; we could model these interactions either with a model type that can handle it natively (like trees) or with explicit interaction terms like this:

rec_interact <-
  rec_basic %>%
  step_interact(~ starts_with("urban_rural"):starts_with("bird_type"))

wf_interact <- workflow(rec_interact, glm_spec)
rs_interact <- fit_resamples(wf_interact, bird_folds, control = ctrl_preds)

How did this do, our same logistic regression model specification but now with interactions?

collect_metrics(rs_interact)

## # A tibble: 2 × 6
##   .metric  .estimator  mean     n   std_err .config             
##   <chr>    <chr>      <dbl> <int>     <dbl> <chr>               
## 1 accuracy binary     0.822    10 0.0000762 Preprocessor1_Model1
## 2 roc_auc  binary     0.669    10 0.00660   Preprocessor1_Model1

The accuracy is about the same (since the model is always predicting “no bird”) but the probabilities look better.

augment(rs_interact) %>%
  roc_curve(bird_count, .pred_bird) %>%
  autoplot()

Evaluate model on new data

Let’s stick with this model, logistic regression together with interactions between urban/rural and bird type. We can fit the model one time to the entire training set.

bird_fit <- fit(wf_interact, bird_train)

Now this trained model is ready to be applied to new data. For example, we can predict the test set, perhaps to get out probabilities.

predict(bird_fit, bird_test, type = "prob")

## # A tibble: 3,571 × 2
##    .pred_bird `.pred_no bird`
##         <dbl>           <dbl>
##  1     0.213            0.787
##  2     0.123            0.877
##  3     0.141            0.859
##  4     0.283            0.717
##  5     0.119            0.881
##  6     0.252            0.748
##  7     0.0380           0.962
##  8     0.123            0.877
##  9     0.129            0.871
## 10     0.119            0.881
## # … with 3,561 more rows

In fact, we can predict on any kind of new data that has the right input variables. Let’s make some ourselves.

new_bird_data <-
  tibble(bird_type = top_birds) %>%
  crossing(urban_rural = c("Urban", "Rural"))

new_bird_data

## # A tibble: 30 × 2
##    bird_type         urban_rural
##    <chr>             <chr>      
##  1 Australian Magpie Rural      
##  2 Australian Magpie Urban      
##  3 Crested Pigeon    Rural      
##  4 Crested Pigeon    Urban      
##  5 Crimson Rosella   Rural      
##  6 Crimson Rosella   Urban      
##  7 Eastern Spinebill Rural      
##  8 Eastern Spinebill Urban      
##  9 Grey Fantail      Rural      
## 10 Grey Fantail      Urban      
## # … with 20 more rows

We can use a helpful function like augment() to take this new data and “augment” it with predicted probabilities and class predictions, and we can use predict() with specific type arguments to return specialized predictions like confidence intervals. Let’s bind these together.

bird_preds <-
  augment(bird_fit, new_bird_data) %>%
  bind_cols(
    predict(bird_fit, new_bird_data, type = "conf_int")
  )

bird_preds

## # A tibble: 30 × 9
##    bird_type urban_rural .pred_class .pred_bird `.pred_no bird` .pred_lower_bird
##    <chr>     <chr>       <fct>            <dbl>           <dbl>            <dbl>
##  1 Australi… Rural       no bird         0.245            0.755           0.193 
##  2 Australi… Urban       no bird         0.287            0.713           0.249 
##  3 Crested … Rural       no bird         0.0826           0.917           0.0526
##  4 Crested … Urban       no bird         0.141            0.859           0.113 
##  5 Crimson … Rural       no bird         0.215            0.785           0.166 
##  6 Crimson … Urban       no bird         0.123            0.877           0.0969
##  7 Eastern … Rural       no bird         0.283            0.717           0.227 
##  8 Eastern … Urban       no bird         0.0973           0.903           0.0736
##  9 Grey Fan… Rural       no bird         0.254            0.746           0.200 
## 10 Grey Fan… Urban       no bird         0.0614           0.939           0.0435
## # … with 20 more rows, and 3 more variables: .pred_upper_bird <dbl>,
## #   .pred_lower_no bird <dbl>, .pred_upper_no bird <dbl>

Now let’s visualize these predictions.

p2 <-
  bird_preds %>%
  ggplot(aes(.pred_bird, bird_type, color = urban_rural)) +
  geom_errorbar(aes(
    xmin = .pred_lower_bird,
    xmax = .pred_upper_bird
  ),
  width = .2, size = 1.2, alpha = 0.5
  ) +
  geom_point(size = 2.5) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = "Predicted probability of seeing bird", y = NULL, color = NULL)

p2

Actually, let’s put this together with our earlier plot!

library(patchwork)

p1 + p2