Linear regression

To provide a baseline, I first run a linear regression

library(broom)

estimateDF %>%
  mutate( room_type = factor(room_type),
          minimum_nights = log10(1+minimum_nights)) %>%
  lm( price ~ room_type + latitude + longitude + 
        minimum_nights + availability_365, data=.) -> linreg

tidy(linreg) %>% print()

## # A tibble: 7 × 5
##   term                     estimate  std.error statistic   p.value
##   <chr>                       <dbl>      <dbl>     <dbl>     <dbl>
## 1 (Intercept)           -326.       13.7          -23.9  2.75e-119
## 2 room_typePrivate room   -0.791     0.0152       -52.1  0        
## 3 room_typeShared room    -1.19      0.0491       -24.3  1.58e-123
## 4 latitude                 1.43      0.136         10.5  1.60e- 25
## 5 longitude               -3.69      0.163        -22.6  5.83e-108
## 6 minimum_nights          -0.191     0.0201        -9.50 3.27e- 21
## 7 availability_365         0.000728  0.0000566     12.9  2.69e- 37

glance(linreg) %>% print()

## # A tibble: 1 × 12
##   r.squ…¹ adj.r…² sigma stati…³ p.value    df logLik   AIC   BIC devia…⁴ df.re…⁵
##     <dbl>   <dbl> <dbl>   <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl>   <int>
## 1   0.459   0.458 0.512    705.       0     6 -3748. 7513. 7565.   1311.    4993
## # … with 1 more variable: nobs <int>, and abbreviated variable names
## #   ¹​r.squared, ²​adj.r.squared, ³​statistic, ⁴​deviance, ⁵​df.residual

The validation RMSE for this model is

validateDF %>%
  mutate( room_type = factor(room_type),
          minimum_nights = log10(1+minimum_nights)) %>%
  mutate( mu = predict(linreg, newdata = .) ) %>%
  summarise(RMSE = sqrt( mean( (price - mu)^2)))

## # A tibble: 1 × 1
##    RMSE
##   <dbl>
## 1 0.492

The validation RMSE is 0.492. For comparison, the leading entry on the sliced leaderboard scored 0.408 with the full test data. My earlier post used XGBoost and achieved a RMSE of 0.431, which I reduced to 0.418 with hyperparameter tuning.

BART

The BART package provides functions for fitting the BART model to different types of response. In the case of a continuous response the appropriate function is called wbart(). To run this function, the data need to be saved in vectors and matrices.

# --- Place the response and predictors in matrices -----

# --- estimation data -----------------------------------
estimateDF %>%
  pull(price) -> Y
  
estimateDF %>%
    mutate( room_type = as.numeric(factor(room_type)),
            minimum_nights = log10(1+minimum_nights)) %>%
    select(room_type, latitude, longitude, minimum_nights, 
           availability_365) %>%
    as.matrix() -> X

# --- validation data -----------------------------------
validateDF %>%
  pull(price) -> YV
  
validateDF %>%
    mutate( room_type = as.numeric(factor(room_type)),
            minimum_nights = log10(1+minimum_nights)) %>%
    select(room_type, latitude, longitude, minimum_nights, 
           availability_365) %>%
    as.matrix() -> XV

I start by running wbart() with default values for all of its parameters. This combines 200 trees over a chain of length 1000 following a burn-in of 100. To help assess convergence, I run three chains with different seeds. Each chain takes about 20s to run on my desktop.

My function reportBART() summarises the results of a chain. The function’s code is given in the accompanying methods post.

# --- Chain 1 -----------------------------------------
set.seed(4592)
bt1 <- wbart(x.train = X, y.train = Y, x.test = XV)
reportBART(bt1, Y, YV)
# --- Chain 2 -----------------------------------------
set.seed(2893)
bt2 <- wbart(x.train = X, y.train = Y, x.test = XV)
reportBART(bt2, Y, YV)
# --- Chain 3 -----------------------------------------
set.seed(7872)
bt3 <- wbart(x.train = X, y.train = Y, x.test = XV)
reportBART(bt3, Y, YV)

## 200  trees: run length  1000  burnin  100  thin by  1 
##  Posterior mean sigma     0.4498403 
##  Average splits per tree  1.193965 
##  training RMSE            0.4418989 
##  test RMSE                0.4539434

## 200  trees: run length  1000  burnin  100  thin by  1 
##  Posterior mean sigma     0.4502711 
##  Average splits per tree  1.18779 
##  training RMSE            0.4418881 
##  test RMSE                0.4533676

## 200  trees: run length  1000  burnin  100  thin by  1 
##  Posterior mean sigma     0.4508597 
##  Average splits per tree  1.17109 
##  training RMSE            0.442924 
##  test RMSE                0.4539443

The validation RMSE (labelled test RMSE by my function) is 0.454, better than linear regression, but worse than my XGBoost model. In fact, 0.454 is embarrassingly good, if that performance were to carry over to the test data, it would rank 6th on the private leaderboard.

Notice that there is little evidence of over-fitting, as the training and test RMSE are similar. The summary tables also show that the trees generated by the algorithm are very stumpy, they average under 1.2 splits per tree.

Before using this model, I need to check that the algorithm has converged. I discuss convergence checking in the methods post, so here I will just present results.

Here are the trace plots for sigma, the inherent variation about the trend.

Mixing is quite good but the 3 chains are still drifting down, which suggests that I need a longer burn-in.

Here is the trace plot for 3 chains with the burn-in increased for 100 to 5000. These chains each take about 100 seconds to run.

Agreement is better, but not ideal. The validation (test) RMSE has hardly changed; it is still around 0.454.

Using code developed in the methods post, I look at the in-sample predictive performance. With 5000 observations the plot is very crowded, so I only show every 20th observation.

Here is the equivalent plot for the validation sample.

In both cases, the performance is much as expected.

wbart() returns an object that includes the number of times that each predictor is used to split a tree. This provides a measure of variable importance. The plot below shows the average usage per tree based on the first of the chains with a 5000 burn-in.

Room type is used more than the other predictors but not by as much as I would have thought. I guess that this is because there are only 3 levels of room type and once a couple of the trees have split on room type, there is nothing to be gained by reusing that predictor. This is a weakness of using the number of splits as a measure of variable importance.

Sparsity-inducing priors

So far, I have not changed any of the default priors. In my limited experience, the defaults are well-chosen and making small changes to them does not noticeably alter the fitted model.

Even changing the number of trees, only has a small impact on performance. The original papers on the BART algorithm suggested summing over 50 trees, rather than the 200 that BART defaults to. When the algorithm is limited to 50 trees, it tends to make those trees deeper, but ends with a similar predictive performance. I think that 200 trees is a marginally better choice.

One interesting choice made when setting the priors is the probability with which the predictors are chosen when the model seeks to extend a tree. The default is to make each predictor equally likely to be tried, though of course, poorly performing predictors will get rejected and will not make it into the final model.

In 2018, Antonio Linero published a paper that advocated the use of sparsity-inducing priors that encourage the model to use a smaller set of predictors.

Linero, A. R. (2018).
Bayesian regression trees for high-dimensional prediction and variable selection.
Journal of the American Statistical Association, 113(522), 626-636.

Suppose that there are J potential predictors (J=315 for my model). When the default analysis picks a predictor to add to the tree, each one has a probability 1/J of being chosen. Instead, Linero suggests making the selection probabilities follow a Dirichlet distribution with equal parameters \(\theta/J\). The model then places a beta prior on \(\theta\) and treats the selection probabilities of the predictors as parameters to be learnt.

The form of the prior on \(\theta\) is, \[ \frac{\theta}{\theta+\rho} \sim \text{Beta}(a, b) \] where the user chooses a, b and \(\rho\). The defaults are a=0.5, b=1, \(\rho\)=J.

Although, the Dirichlet prior starts by treating all predictors equally, the model is free to learn which predictors are most useful and to use that information when selecting a new predictor to use in the BART model. By changing, a, b and \(\rho\), the user changes the sparsity of the selection. Once BART has identified the important predictors, it tends to stick with them.

I tried changing a and b but found that they have little impact on the airbnb model. However, reducing \(\rho\) does make the predictor selection more sparse. Below, I show the results for \(\rho\)=100.

# --- single sparse chain -----------------------------
set.seed(4592)
bt1 <- wbart(x.train = X, y.train = Y, x.test = XV,
             nskip=10000, sparse=TRUE, rho=100)
reportBART(bt1, Y, YV)

## 200  trees: run length  1000  burnin  10000  thin by  1 
##  Posterior mean sigma     0.408865 
##  Average splits per tree  1.28454 
##  training RMSE            0.3983437 
##  test RMSE                0.4229048

The chain took about 3.5 minutes to run.

The tree depths have increased very slightly and the test RMSE is slightly better. Remember that a default XGBoost model had a RMSE of 0.431 (though that was with the sliced test data).

Here are the frequencies for the top predictors

## # A tibble: 39 × 2
##    var                             freq
##    <chr>                          <dbl>
##  1 calculated_host_listings_count 37.6 
##  2 room_type                      36.2 
##  3 longitude                      21.4 
##  4 availability_365               19.4 
##  5 N78                            16.9 
##  6 latitude                        9.49
##  7 X11                             6.30
##  8 X18                             5.95
##  9 X3                              5.71
## 10 minimum_nights                  5.36
## # … with 29 more rows

Notice that the frequencies have increased dramatically, so that fewer predictors are used.

Finally, here are the selection probabilities that the model learns from its Dirichlet prior. The order is more or less the same as that for the frequencies.

## # A tibble: 21 × 2
##    var                              prob
##    <chr>                           <dbl>
##  1 calculated_host_listings_count 0.135 
##  2 room_type                      0.130 
##  3 longitude                      0.0772
##  4 availability_365               0.0697
##  5 N78                            0.0606
##  6 latitude                       0.0349
##  7 X11                            0.0228
##  8 X18                            0.0212
##  9 X3                             0.0207
## 10 minimum_nights                 0.0197
## # … with 11 more rows