Internship ML Prep: Data with spatial coordinates and text features

In this episode of SLICED, contestants are challenged to predict the price of Airbnb listings in NYC.

The evaluation metric for this competition is Root Mean Squared Logarithmic Error.

Split Data

set.seed(2021)
# Load required libraries
library(tidyverse)
library(tidymodels)
library(here)
library(textrecipes)
library(stacks)
library(broom)

# Import data
play_data <- read_csv("train.csv", show_col_types = FALSE) %>% 
  # Account for 0 by adding 1 
  mutate(price = log(price + 1))

holdout_data <- read_csv("test.csv", show_col_types = FALSE)

# Function for predicting on the hold out data
predict_holdout <- function(wf){
  wf %>% augment(holdout_data) %>% 
    mutate(.pred = exp(.pred) - 1) %>% 
    select(id, price = .pred)
}


# Function to make augment() work with stacks
augment.model_stack <- function(stack_object, data){
  data %>% 
    bind_cols(predict(stack_object, .))
}


# Split data
data_split <- initial_split(play_data, prop = 0.75)

# Extract data from split data
train <- training(data_split)
test <- testing(data_split)

EDA

# Brief summary of the variables
train %>% 
  glimpse()

## Rows: 25,669
## Columns: 16
## $ id                             <dbl> 26155827, 2298373, 2954785, 18550867, 6~
## $ name                           <chr> "Greenpoint room", "Luxury private apar~
## $ host_id                        <dbl> 5213686, 10149317, 15081041, 128785279,~
## $ host_name                      <chr> "Charlotte", "Lana", "Jaclyn", "Dusean"~
## $ neighbourhood_group            <chr> "Brooklyn", "Queens", "Manhattan", "Bro~
## $ neighbourhood                  <chr> "Greenpoint", "Rego Park", "Murray Hill~
## $ latitude                       <dbl> 40.73442, 40.72988, 40.74551, 40.85878,~
## $ longitude                      <dbl> -73.95854, -73.85773, -73.97757, -73.83~
## $ room_type                      <chr> "Private room", "Entire home/apt", "Ent~
## $ price                          <dbl> 4.262680, 4.488636, 5.135798, 5.298317,~
## $ minimum_nights                 <dbl> 5, 3, 5, 2, 2, 3, 3, 20, 30, 2, 1, 29, ~
## $ number_of_reviews              <dbl> 2, 33, 0, 81, 203, 2, 15, 0, 1, 6, 165,~
## $ last_review                    <date> 2018-08-09, 2019-04-20, NA, 2019-06-26~
## $ reviews_per_month              <dbl> 0.16, 0.77, NA, 3.17, 1.96, 1.13, 1.18,~
## $ calculated_host_listings_count <dbl> 1, 5, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, ~
## $ availability_365               <dbl> 0, 262, 0, 286, 86, 4, 365, 40, 342, 5,~

Distribution of what we are trying to predict

theme_set(theme_light())
train %>% 
  # Has 0s in it so we add 1
  ggplot(mapping = aes(x = price)) +
  geom_histogram(fill = "midnightblue", alpha = 0.7) 

# Find median and average prices
summary_price <- function(tbl){
  tbl %>% 
    summarize(avg_price = exp(mean(price))-1,
              n = n()) %>% arrange(desc(n))
              
}

train %>% 
  summary_price()

# How average price varies per neighbourhood_group
train %>% 
  mutate(neighbourhood_group = fct_reorder(neighbourhood_group, price)) %>% 
  ggplot(mapping = aes(x = exp(price), y = neighbourhood_group )) +
  geom_boxplot() +
  scale_x_log10()

# How average price varies per neighbourhood
train %>% 
  mutate(neighbourhood = fct_lump(neighbourhood, 40),
    neighbourhood = fct_reorder(neighbourhood, price)) %>% 
  ggplot(mapping = aes(x = exp(price), y = neighbourhood )) +
  geom_boxplot() +
  scale_x_log10()

# How price varies with room type
train %>% 
  mutate(room_type = fct_reorder(room_type, price)) %>% 
  ggplot(mapping = aes(x = exp(price), y = room_type )) +
  geom_boxplot() +
  scale_x_log10()

# How price varies with reviews per month
train %>% 
  ggplot(mapping = aes(x = number_of_reviews, y = price )) +
  geom_point() +
  scale_x_log10() +
  geom_smooth(method = "lm")

# How price varies with reviews per month
train %>% 
  ggplot(mapping = aes(x = reviews_per_month, y = price )) +
  geom_point() +
  scale_x_log10() +
  geom_smooth(method = "lm")

# How price varies with last review
train %>% 
  ggplot(mapping = aes(x = last_review, y = price )) +
  geom_point() +
  geom_smooth(method = "lm")

# How price varies with number of minimum nights
train %>% 
  slice_sample(n = 1000) %>% 
  ggplot(mapping = aes(x = minimum_nights, y = price )) +
  geom_point() +
  scale_x_log10() +
  geom_smooth(method = "lm")

Text Analysis

library(tidytext)

# Word tokenization to see relationship between word and price
train %>% 
  unnest_tokens(input = name, output = word, token = "words") %>% 
  group_by(word) %>% 
  summary_price() %>% 
  head(50) %>% 
  mutate(word = fct_reorder(word, avg_price)) %>% 
  ggplot(aes(x = avg_price, y = word, size = n)) +
  geom_point()

Simple map

library(ggmap)
library(ggthemes)

# Aggregate nearby areas
agg_lat_long <- train %>% 
  group_by(latitude = round(latitude, 2),
           longitude = round(longitude, 2)) %>% 
  summarize(price = mean(price), n = n()) %>% 
  filter(n >= 5)

# Find range of coverage
range(train %>% pull(latitude))

## [1] 40.50641 40.91306

range(train %>% pull(longitude))

## [1] -74.24285 -73.71690

bbox <- c(left = -74.24285,  bottom = 40.50641, right = -73.71690, top = 40.91306 )
nyc_map <- get_stamenmap(bbox, zoom = 11)

# Which areas have high prices
ggmap(nyc_map) +
  geom_point(data = agg_lat_long, aes(longitude, latitude, color = exp(price) - 1, size = n)) +
  scale_color_gradient2(low = "blue", high = "red", midpoint = 2, trans = "log10") +
  theme_map() +
  labs(color = "Price")

• Presence of certain words really does matter in terms of pricing

• Areas in certain long and lat have distinctive high prices

• Seems there is a disparity of prices depending on the neighbourhood group and neighbourhood

• The room type clearly partitions the house prices

• The other features don’t seem to influence the price so much ?

Train xgboost model

Creates a series of decision trees forming an ensemble. Each tree depends on the results of previous trees in an attempt to reduce the error. All trees in the ensemble are combined to produce a final prediction.

# Control aspects of a grid search
eval_met <- metric_set(rmse)
grid_control <- control_grid(
  # Save predictions for each model evaluated
  save_pred = TRUE,
  # Save workflow (model+recipe)
  save_workflow = TRUE,
  # Extract a fitted model
  extract = extract_model)

# Resampling techniques for creating simulated data sets
set.seed(2021)
train_5fold_cv <- train %>% 
  vfold_cv(v = 5, repeats = 1)
train_5fold_cv

• Does xgboost require us to normalize or transform variables to be more symmetric? Not really. More research on this

• While tree-based boosting methods generally do not require the creation of dummy variables, models using the xgboost engine do.

# Function that applies preprocessing steps to data (sanity check)
prep_bake <- function(recipe) juice(prep(recipe))

# Data preprocessing with recipes - MOSTLY NUMERIC
xg_rec <- recipe(price ~ minimum_nights + room_type + number_of_reviews + latitude + longitude + neighbourhood_group + reviews_per_month + calculated_host_listings_count + availability_365 + last_review, data = train) %>% 
  # Make variables symmetrical
  #step_log(all_nominal_predictors(), offset = 1) %>% 
  # if not ever reviewed, 0
  # Narrow down to Manhattan as suggested by the VIP plot below
  step_mutate(manhattan_grp = neighbourhood_group == "Manhattan") %>% 
  step_rm(neighbourhood_group) %>% 
  step_mutate(last_review = coalesce(as.integer(Sys.Date() - last_review), 0)) %>% 
  step_dummy(all_nominal_predictors())

# Sanity check
# xg_rec %>% 
#   prep_bake() %>% 
#   slice_head(n = 5)

# Make a model specification
xg_mod <- boost_tree(
  # Randomly selected predictors at each split
  mtry = tune(),
  # Number of trees
  trees = tune(),
  learn_rate = 0.01
) %>% 
  set_engine("xgboost") %>% 
  set_mode("regression")

# Combine model specifcation and recipe into a workflow
xg_wf <- workflow() %>% 
  add_recipe(xg_rec) %>% 
  add_model(xg_mod)

# Tune model hyperparameters
doParallel::registerDoParallel()
xg_tune <- xg_wf %>% 
  tune_grid(
    resamples = train_5fold_cv,
    metrics = eval_met,
    control = grid_control,
    grid = crossing(mtry = c(4, 5, 6,  7), trees = seq(250, 1000, 25))
  )

xg_tune %>% autoplot()

# See performance metrics
xg_tune %>% 
  collect_metrics() %>% 
  arrange((mean))

• RMSE of 0.535 with 600 trees and mtry = 2: price ~ minimum_nights + room_type + number_of_reviews

• RMSE of 0.468 with 1000 trees and mtry = 6: price ~ minimum_nights + room_type + number_of_reviews + latitude + longitude + neighbourhood_group - do we need 1000 trees? Overfitting??

• RMSE of 0.442 with 800 trees and mtry = 6: price ~ minimum_nights + room_type + number_of_reviews + latitude + longitude + neighbourhood_group + reviews_per_month + calculated_host_listings_count + availability_365

• RMSE of 0.4418 with 800 trees and mtry = 7: price ~ minimum_nights + room_type + number_of_reviews + latitude + longitude + neighbourhood_group + reviews_per_month + calculated_host_listings_count + availability_365 + last_review — Finalize this

• RMSE of 0.4413 with 1000 trees and mtry = 7: price ~ minimum_nights + room_type + number_of_reviews + latitude + longitude + neighbourhood_group + reviews_per_month + calculated_host_listings_count + availability_365 + last_review — with manhattan as a feature

Finalize version without many categorical values

# Finalize workflow and fit model
xg_fit <- xg_wf %>% 
  finalize_workflow(select_best(xg_tune)) %>% 
  fit(train)

# Evaluate model performance on test set
xg_fit %>% 
  augment(test) %>% 
  rmse(truth = price, estimate = .pred)

# Variable importance plot
important_var <- xgboost::xgb.importance(model = xg_fit$fit$fit$fit)

important_var %>% 
  mutate(Feature = fct_reorder(Feature, Gain)) %>% 
  ggplot(mapping = aes(x = Gain, y = Feature)) +
  geom_col()

# Make predictions on holdout data
xg_fit %>% 
  predict_holdout() %>% 
  write_csv("xg_numeric")

The Manhattan neighbourhood group stands out from the rest of its peers.

Linear version on the categorical

• Additionally, step_other() can be used to analyze the frequencies of the factor levels in the training set and convert infrequently occurring values to a catch-all level of “other”, with a specific threshold that can be specified.

  Lumps infrequently occurring observations into a new level called other.

• Normalize to compare them on the same scale later

# Recipe for preprocessing data
lin_rec <- recipe(price ~ name + room_type +
                    latitude + longitude +
                    neighbourhood_group + neighbourhood +
                    host_id, data = train) %>% 
  # Tokenize names
  step_tokenize(name) %>% 
  # Narrow down to top 100 words?
  step_tokenfilter(name, max_tokens = tune()) %>% 
  # Find if a token is found in the column - dumminise 
  step_tf(name) %>% 
  step_mutate(host_id = factor(host_id)) %>% 
  step_other(host_id, neighbourhood, threshold = tune()) %>% 
  step_dummy(all_nominal_predictors())

# Linear model specification
lin_mod <- linear_reg(penalty = tune()) %>% 
  set_engine("glmnet") %>% 
  set_mode("regression")

# Bundle model spec + recipe into workflow
lin_wf <- workflow() %>% 
  add_recipe(lin_rec) %>% 
  add_model(lin_mod)

# Model tuning
# Tune model hyperparameters
doParallel::registerDoParallel()
lin_tune <- lin_wf %>% 
  tune_grid(
    resamples = train_5fold_cv,
    metrics = eval_met,
    control = grid_control,
    grid = crossing(penalty = 10 ^ seq(-7, -1, .1),
                    threshold = c(0.0008,0.001),
                    max_tokens = c(100, 500)))

autoplot(lin_tune)

# Collect metrics
lin_tune %>% 
  collect_metrics() %>% 
  arrange(mean)

• About 0.48 on room_type + latitude + longitude + neighbourhood_group + neighbourhood + host_id . Seems lower threshold(step_other) does better

• Bringing in top 100 words improved performance. Trying to tune the max_tokens. Settle for a threshold of 0.001 (0.1%?) Resulted in 0.45

# Finalize worklow and train model
lin_fit <- lin_wf %>% 
  finalize_workflow(select_best(lin_tune)) %>% 
  fit(train)

# Make predictions on test set
lin_fit %>% 
  augment(test) %>% 
  rmse(truth = price, estimate = .pred)

Wow! Almost as good as xg_boost – 0.437

Now, let’s interpret the model and see which variables contributed to prediction

# Interpret model
library(broom)

lin_fit$fit$fit$fit %>% 
  tidy() %>% 
  filter(lambda >= select_best(lin_tune)$penalty) %>% 
  filter(lambda == min(lambda), term!= "(Intercept)") %>% 
  slice_max(order_by = abs(estimate), n = 50) %>% 
  mutate(term = fct_reorder(term, estimate)) %>% 
  ggplot(aes(estimate, term, fill = estimate > 0)) +
  geom_col(show.legend = F)+
  xlab("Coefficient in regularized linear model")

• Name turned out to be important!

• As longitude decreases price decreases? Hidden places fetch more price? 4 bedroom houses are pricier. Makes sense.

Stack a boosted tree model and a linear model

• Model stacks can be thought of as a group of fitted member models and a set of instructions on how to combine their predictions.

• Library stacks requires you to save predictions and model as we did in the control_grid

# Keep only the results from the numerically best combination
lin_best <- lin_tune %>% 
  filter_parameters(parameters = select_best(lin_tune))
xg_best <- xg_tune %>% 
  filter_parameters(parameters = select_best(xg_tune))

# Create a stack of models
lin_xg_stack <- stacks() %>%
  add_candidates(lin_best) %>% 
  add_candidates(xg_best)
  
lin_xg_stack

Under the hood, a data_stack object is really just a tibble with some extra attributes. Checking out the actual data:

lin_xg_stack %>% 
  as_tibble() %>% 
  slice_head(n=5)

The first column gives the first response value, and the remaining columns give the assessment set predictions for each ensemble member.

That’s it! We’re now ready to evaluate how it is that we need to combine predictions from each candidate ensemble member.

Fit the stack

The outputs from each of these candidate ensemble members are highly correlated, so the blend_predictions() function performs regularization to figure out how we can combine the outputs from the stack members to come up with a final prediction.

The blend_predictions function determines how member model output will ultimately be combined in the final prediction by fitting a LASSO model on the data stack, predicting the true assessment set outcome using the predictions from each of the candidate members. Candidates with nonzero stacking coefficients become members.

# Fit the stack
blended_lin_xg <- lin_xg_stack %>% 
  blend_predictions()

blended_lin_xg

## # A tibble: 2 x 3
##   member        type       weight
##   <chr>         <chr>       <dbl>
## 1 xg_best_1_124 boost_tree  0.573
## 2 lin_best_2_42 linear_reg  0.492

Now that we know how to combine our model output, we can fit the candidates with non-zero stacking coefficients on the full training set.

blended_lin_xg_fit <- blended_lin_xg %>% 
  fit_members()

blended_lin_xg_fit

## # A tibble: 2 x 3
##   member        type       weight
##   <chr>         <chr>       <dbl>
## 1 xg_best_1_124 boost_tree  0.573
## 2 lin_best_2_42 linear_reg  0.492

class(blended_lin_xg_fit)

## [1] "linear_stack" "model_stack"  "list"

The object is now ready to predict with new data

# Make predictions on test set
test_stacks <- test %>% 
  bind_cols(predict(blended_lin_xg_fit, .))

# Evaluate stack
test_stacks %>% 
  rmse(truth = price, estimate = .pred)

Yes! We have lowered rmse! Good job!

Little trick by David to get more power out of stacked models. Retrain the stacks on the whole data set

# Make a copy of blended stack
blended_lin_xg_fulldata <- blended_lin_xg

# Substitute the train component with the whole initial data
blended_lin_xg_fulldata$train <- play_data

# Train the stacks object on the whole data
blended_lin_xg_fulldata_fit <- blended_lin_xg_fulldata %>% 
  fit_members()

# Function to make augment() work with stacks
augment.model_stack <- function(stack_object, data){
  data %>% 
    bind_cols(predict(stack_object, .))
}

# Predict on holdout
blended_lin_xg_fulldata_fit %>% 
  predict_holdout() %>% 
  write_csv("lin_xg_stack_pred.csv")

Further reading

https://www.khanacademy.org/math/algebra-home/alg-exp-and-log/alg-logarithmic-scale/v/logarithmic-scale

https://stats.stackexchange.com/questions/298/in-linear-regression-when-is-it-appropriate-to-use-the-log-of-an-independent-va

Inverse Document Frequency: https://www.tidytextmining.com/tfidf.html

Getting started with stacks: https://stacks.tidymodels.org/articles/basics.html