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. My screencasts lately have focused on xgboost as I have participated in SLICED, a competitive data science streaming show. This past week were the semifinals, where we competed to predict prices of homes in Austin, TX. ? One of the more interesting available variables for this dataset was the text description of the real estate listings, so let’s walk through one way to incorporate text information with boosted tree modeling.

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 the price (binned) for homes in Austin, TX given features about the real estate listing. This is a multiclass classification challenge, where we needed to submit a probability for each home being in each priceRange bin. The main data set provided is in a CSV file called training.csv.

library(tidyverse)
train_raw <- read_csv("train.csv")

train_raw %>%
  count(priceRange)

## # A tibble: 5 × 2
##   priceRange        n
##   <chr>         <int>
## 1 0-250000       1249
## 2 250000-350000  2356
## 3 350000-450000  2301
## 4 450000-650000  2275
## 5 650000+        1819

You can watch this week’s full episode of SLICED to see lots of exploratory data analysis and visualization of this dataset, but let’s just make a few data visualization for context in this blog post.

How is price distributed across Austin?

price_plot <-
  train_raw %>%
  mutate(priceRange = parse_number(priceRange)) %>%
  ggplot(aes(longitude, latitude, z = priceRange)) +
  stat_summary_hex(alpha = 0.8, bins = 50) +
  scale_fill_viridis_c() +
  labs(
    fill = "mean",
    title = "Price"
  )

price_plot

Let’s look at this distribution and compare it to some other variables available in the dataset. We can create a little plotting function using {{}} to quickly iterate through, and put them together with patchwork.

library(patchwork)

plot_austin <- function(var, title) {
  train_raw %>%
    ggplot(aes(longitude, latitude, z = {{ var }})) +
    stat_summary_hex(alpha = 0.8, bins = 50) +
    scale_fill_viridis_c() +
    labs(
      fill = "mean",
      title = title
    )
}

(price_plot + plot_austin(avgSchoolRating, "School rating")) /
  (plot_austin(yearBuilt, "Year built") + plot_austin(log(lotSizeSqFt), "Lot size (log)"))

Notice the east/west gradients as well as the radial changes. I went to grad school in Austin and this all look very familiar to me!

Finding words related to price

The description variable contains text from each real estate listing. We could try to use the text features directly in modeling, as described in our book, but I’ve found that often isn’t great for boosted tree models (which tend to be what works best overall in an environment like SLICED). Let’s walk through another option which may work better in some situations, which is to use some separate analysis to identify important words and then create dummy variables indicating whether any given listing has those words.

Let’s start by tidying the description text.

library(tidytext)

austin_tidy <-
  train_raw %>%
  mutate(priceRange = parse_number(priceRange) + 100000) %>%
  unnest_tokens(word, description) %>%
  anti_join(get_stopwords())

austin_tidy %>%
  count(word, sort = TRUE)

## # A tibble: 17,944 × 2
##    word         n
##    <chr>    <int>
##  1 home     11620
##  2 kitchen   5721
##  3 room      5494
##  4 austin    4918
##  5 new       4772
##  6 large     4771
##  7 2         4585
##  8 bedrooms  4571
##  9 contains  4413
## 10 3         4386
## # … with 17,934 more rows

Next, let’s compute word frequencies per price range for the top 100 words.

top_words <-
  austin_tidy %>%
  count(word, sort = TRUE) %>%
  filter(!word %in% as.character(1:5)) %>%
  slice_max(n, n = 100) %>%
  pull(word)


word_freqs <-
  austin_tidy %>%
  count(word, priceRange) %>%
  complete(word, priceRange, fill = list(n = 0)) %>%
  group_by(priceRange) %>%
  mutate(
    price_total = sum(n),
    proportion = n / price_total
  ) %>%
  ungroup() %>%
  filter(word %in% top_words)


word_freqs

## # A tibble: 500 × 5
##    word       priceRange     n price_total proportion
##    <chr>           <dbl> <dbl>       <dbl>      <dbl>
##  1 access         100000   180       56290    0.00320
##  2 access         350000   365      114853    0.00318
##  3 access         450000   322      116678    0.00276
##  4 access         550000   294      125585    0.00234
##  5 access         750000   248      112073    0.00221
##  6 appliances     100000   209       56290    0.00371
##  7 appliances     350000   583      114853    0.00508
##  8 appliances     450000   576      116678    0.00494
##  9 appliances     550000   567      125585    0.00451
## 10 appliances     750000   391      112073    0.00349
## # … with 490 more rows

Now let’s use modeling to find the words that are increasing with price and those that are decreasing with price.

word_mods <-
  word_freqs %>%
  nest(data = c(priceRange, n, price_total, proportion)) %>%
  mutate(
    model = map(data, ~ glm(cbind(n, price_total) ~ priceRange, ., family = "binomial")),
    model = map(model, tidy)
  ) %>%
  unnest(model) %>%
  filter(term == "priceRange") %>%
  mutate(p.value = p.adjust(p.value)) %>%
  arrange(-estimate)

word_mods

## # A tibble: 100 × 7
##    word     data             term          estimate std.error statistic  p.value
##    <chr>    <list>           <chr>            <dbl>     <dbl>     <dbl>    <dbl>
##  1 outdoor  <tibble [5 × 4]> priceRange 0.00000325    1.85e-7     17.6  4.37e-67
##  2 custom   <tibble [5 × 4]> priceRange 0.00000214    1.47e-7     14.6  3.98e-46
##  3 pool     <tibble [5 × 4]> priceRange 0.00000159    1.22e-7     13.0  6.12e-37
##  4 office   <tibble [5 × 4]> priceRange 0.00000150    1.46e-7     10.3  6.03e-23
##  5 suite    <tibble [5 × 4]> priceRange 0.00000143    1.39e-7     10.3  4.03e-23
##  6 gorgeous <tibble [5 × 4]> priceRange 0.000000975   1.62e-7      6.02 1.19e- 7
##  7 w        <tibble [5 × 4]> priceRange 0.000000920   9.05e-8     10.2  2.33e-22
##  8 windows  <tibble [5 × 4]> priceRange 0.000000890   1.28e-7      6.95 2.81e-10
##  9 private  <tibble [5 × 4]> priceRange 0.000000889   1.15e-7      7.70 1.08e-12
## 10 car      <tibble [5 × 4]> priceRange 0.000000778   1.66e-7      4.69 1.52e- 4
## # … with 90 more rows

Let’s make something like a volcano plot to see the relationship between p-value and effect size for these words.

library(ggrepel)

word_mods %>%
  ggplot(aes(estimate, p.value)) +
  geom_vline(xintercept = 0, lty = 2, alpha = 0.7, color = "gray50") +
  geom_point(color = "midnightblue", alpha = 0.8, size = 2.5) +
  scale_y_log10() +
  geom_text_repel(aes(label = word), family = "IBMPlexSans")

• Words like outdoor, custom, pool, suite, office increase with price.
• Words like new, paint, carpet, great, tile, close, flooring decrease with price.

These are the words that we’d like to try to detect and use in feature engineering for our xgboost model, rather than using all the text tokens as features individually.

higher_words <-
  word_mods %>%
  filter(p.value < 0.05) %>%
  slice_max(estimate, n = 12) %>%
  pull(word)

lower_words <-
  word_mods %>%
  filter(p.value < 0.05) %>%
  slice_max(-estimate, n = 12) %>%
  pull(word)

We can look at these changes with price directly. For example, these are the words most associated with price decrease.

word_freqs %>%
  filter(word %in% lower_words) %>%
  ggplot(aes(priceRange, proportion, color = word)) +
  geom_line(size = 2.5, alpha = 0.7, show.legend = FALSE) +
  facet_wrap(vars(word), scales = "free_y") +
  scale_x_continuous(labels = scales::dollar) +
  scale_y_continuous(labels = scales::percent, limits = c(0, NA)) +
  labs(x = NULL, y = "proportion of total words used for homes at that price") +
  theme_light(base_family = "IBMPlexSans")

Cheaper houses are “great” but not expensive houses, and apparently you don’t need to mention the location (“close,” “minutes,” “location”) of more expensive houses.

Build a model

Let’s start our modeling by setting up our “data budget,” as well as the metrics (this challenge was evaluate on multiclass log loss).

library(tidymodels)

set.seed(123)
austin_split <- train_raw %>%
  select(-city) %>%
  mutate(description = str_to_lower(description)) %>%
  initial_split(strata = priceRange)
austin_train <- training(austin_split)
austin_test <- testing(austin_split)
austin_metrics <- metric_set(accuracy, roc_auc, mn_log_loss)

set.seed(234)
austin_folds <- vfold_cv(austin_train, v = 5, strata = priceRange)
austin_folds

## #  5-fold cross-validation using stratification 
## # A tibble: 5 × 2
##   splits              id   
##   <list>              <chr>
## 1 <split [5996/1502]> Fold1
## 2 <split [5998/1500]> Fold2
## 3 <split [5999/1499]> Fold3
## 4 <split [5999/1499]> Fold4
## 5 <split [6000/1498]> Fold5

For feature engineering, let’s use basically everything in the dataset (aside from city, which was not a very useful variable) and create dummy or indicator variables using step_regex(). The idea here is that we will detect whether these words associated with low/high price are there and create a yes/no variable indicating their presence or absence.

higher_pat <- glue::glue_collapse(higher_words, sep = "|")
lower_pat <- glue::glue_collapse(lower_words, sep = "|")

austin_rec <-
  recipe(priceRange ~ ., data = austin_train) %>%
  update_role(uid, new_role = "uid") %>%
  step_regex(description, pattern = higher_pat, result = "high_price_words") %>%
  step_regex(description, pattern = lower_pat, result = "low_price_words") %>%
  step_rm(description) %>%
  step_novel(homeType) %>%
  step_unknown(homeType) %>%
  step_other(homeType, threshold = 0.02) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_nzv(all_predictors())

austin_rec

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         13
##        uid          1
## 
## Operations:
## 
## Regular expression dummy variable using `outdoor|custom|pool|office|suite|gorgeous|w|windows|private|car|high|full`
## Regular expression dummy variable using `carpet|paint|close|flooring|shopping|new|easy|minutes|tile|great|community|location`
## Delete terms description
## Novel factor level assignment for homeType
## Unknown factor level assignment for homeType
## Collapsing factor levels for homeType
## Dummy variables from all_nominal_predictors()
## Sparse, unbalanced variable filter on all_predictors()

Now let’s create a tunable xgboost model specification, tuning a lot of the important model hyperparameters, and combine it with our feature engineering recipe in a workflow(). We can also create a custom xgb_grid to specify what parameters I want to try out, like not-too-small learning rate, avoiding tree stubs, etc. I chose this parameter grid to get reasonable performance in a reasonable amount of tuning time.

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

xgb_word_wf <- workflow(austin_rec, xgb_spec)

set.seed(123)
xgb_grid <-
  grid_max_entropy(
    tree_depth(c(5L, 10L)),
    min_n(c(10L, 40L)),
    mtry(c(5L, 10L)),
    sample_prop(c(0.5, 1.0)),
    learn_rate(c(-2, -1)),
    size = 20
  )

xgb_grid

## # A tibble: 20 × 5
##    tree_depth min_n  mtry sample_size learn_rate
##         <int> <int> <int>       <dbl>      <dbl>
##  1          7    33     8       0.768     0.0845
##  2         10    33     7       0.928     0.0784
##  3          5    21     6       0.626     0.0868
##  4          9    31     8       0.728     0.0162
##  5          8    35     5       0.666     0.0937
##  6          6    21     5       0.907     0.0105
##  7          6    27     6       0.982     0.0729
##  8          7    33     8       0.936     0.0102
##  9          7    15     5       0.559     0.0182
## 10          6    35     9       0.784     0.0347
## 11          9    39     9       0.737     0.0582
## 12          8    17     8       0.596     0.0818
## 13          9    21     7       0.601     0.0136
## 14          7    15     7       0.763     0.0197
## 15          6    12    10       0.800     0.0569
## 16          9    19     9       0.589     0.0138
## 17         10    14     5       0.829     0.0140
## 18          8    37    10       0.664     0.0202
## 19          5    11     5       0.514     0.0136
## 20         10    38     9       0.962     0.0150

Now we can tune across the grid of parameters and our resamples. Since we are trying quite a lot of hyperparameter combinations, let’s use racing to quit early on clearly bad hyperparameter combinations.

library(finetune)
doParallel::registerDoParallel()

set.seed(234)
xgb_word_rs <-
  tune_race_anova(
    xgb_word_wf,
    austin_folds,
    grid = xgb_grid,
    metrics = metric_set(mn_log_loss),
    control = control_race(verbose_elim = TRUE)
  )

xgb_word_rs

## # Tuning results
## # 5-fold cross-validation using stratification 
## # A tibble: 5 × 5
##   splits              id    .order .metrics          .notes          
##   <list>              <chr>  <int> <list>            <list>          
## 1 <split [5996/1502]> Fold1      3 <tibble [20 × 9]> <tibble [0 × 1]>
## 2 <split [5999/1499]> Fold3      1 <tibble [20 × 9]> <tibble [0 × 1]>
## 3 <split [6000/1498]> Fold5      2 <tibble [20 × 9]> <tibble [0 × 1]>
## 4 <split [5999/1499]> Fold4      4 <tibble [10 × 9]> <tibble [0 × 1]>
## 5 <split [5998/1500]> Fold2      5 <tibble [4 × 9]>  <tibble [0 × 1]>

That takes a little while but we did it!

Evaluate results

First off, how did the “race” go?

plot_race(xgb_word_rs)

We can look at the top results manually as well.

show_best(xgb_word_rs)

## # A tibble: 4 × 11
##    mtry min_n tree_depth learn_rate sample_size .metric   .estimator  mean     n
##   <int> <int>      <int>      <dbl>       <dbl> <chr>     <chr>      <dbl> <int>
## 1     5    14         10     0.0140       0.829 mn_log_l… multiclass 0.920     5
## 2     7    15          7     0.0197       0.763 mn_log_l… multiclass 0.921     5
## 3     5    15          7     0.0182       0.559 mn_log_l… multiclass 0.921     5
## 4     9    19          9     0.0138       0.589 mn_log_l… multiclass 0.923     5
## # … with 2 more variables: std_err <dbl>, .config <chr>

Let’s use last_fit() to fit one final time to the training data and evaluate one final time on the testing data, with the numerically optimal result from xgb_word_rs.

xgb_last <-
  xgb_word_wf %>%
  finalize_workflow(select_best(xgb_word_rs, "mn_log_loss")) %>%
  last_fit(austin_split)

xgb_last

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits              id               .metrics  .notes   .predictions .workflow
##   <list>              <chr>            <list>    <list>   <list>       <list>   
## 1 <split [7498/2502]> train/test split <tibble … <tibble… <tibble [2,… <workflo…

How did this model perform on the testing data, that was not used in tuning/training?

collect_predictions(xgb_last) %>%
  mn_log_loss(priceRange, `.pred_0-250000`:`.pred_650000+`)

## # A tibble: 1 × 3
##   .metric     .estimator .estimate
##   <chr>       <chr>          <dbl>
## 1 mn_log_loss multiclass     0.910

This result is pretty good for a single (not ensembled) model and is a wee bit better than what I did during the SLICED competition. I had an R bomb right as I was finishing up tuning a model just like the one I am demonstrating here!

How does this model perform across the different classes?

collect_predictions(xgb_last) %>%
  conf_mat(priceRange, .pred_class) %>%
  autoplot()

We can also visualize this with an ROC curve.

collect_predictions(xgb_last) %>%
  roc_curve(priceRange, `.pred_0-250000`:`.pred_650000+`) %>%
  ggplot(aes(1 - specificity, sensitivity, color = .level)) +
  geom_abline(lty = 2, color = "gray80", size = 1.5) +
  geom_path(alpha = 0.8, size = 1.2) +
  coord_equal() +
  labs(color = NULL)

Notice that it is easier to identify the most expensive homes but more difficult to correctly classify the less expensive homes.

What features are most important for this xgboost model?

library(vip)
extract_workflow(xgb_last) %>%
  extract_fit_parsnip() %>%
  vip(geom = "point", num_features = 15)

The spatial information in latitude/longitude are by far the most important. Notice that the model uses low_price_words more than it uses, for example, whether there is a spa or whether it is a single family home (as opposed to a townhome or condo). It looks like the model is trying to distinguish some of those lower priced categories. The model does not really use the high_price_words variable, perhaps because it is already easy to find the expensive houses.

The two finalists from SLICED go on to compete next Tuesday, which should be fun and interesting to watch! I have enjoyed the opportunity to participate this season.