Location

I start with the boroughs. I do not know NYC very well but the high prices in Manhattan are not surprising. gm (geometric mean) puts the average price back into dollars.

# --- average prices by borough ------------------------
prep1DF %>%
  group_by( neighbourhood_group ) %>%
  summarise( n = n(),
             m = mean(price),
             s = sd(price),
             gm = exp(m) - 1,
             .groups="drop") 

## # A tibble: 5 x 5
##   neighbourhood_group     n     m     s    gm
##   <chr>               <int> <dbl> <dbl> <dbl>
## 1 Bronx                 760  4.26 0.583  70.1
## 2 Brooklyn            14091  4.58 0.635  96.5
## 3 Manhattan           15146  5.01 0.668 149. 
## 4 Queens               3972  4.37 0.570  78.4
## 5 Staten Island         249  4.41 0.696  81.5

Next the smaller areas

# --- average prices by neighbourhood ------------------
prep1DF %>%
  group_by( neighbourhood ) %>%
  summarise( n = n(),
             m = mean(price),
             s = sd(price),
             gm = exp(m) - 1,
             lat = mean(latitude),
             lng = mean(longitude),
             .groups="drop") %>%
  arrange( desc(m))

## # A tibble: 217 x 7
##    neighbourhood         n     m      s    gm   lat   lng
##    <chr>             <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
##  1 Woodrow               1  6.55 NA      700   40.5 -74.2
##  2 Tribeca             128  5.82  0.813  335.  40.7 -74.0
##  3 Prince's Bay          3  5.60  1.38   270.  40.5 -74.2
##  4 Flatiron District    55  5.59  0.668  268.  40.7 -74.0
##  5 Neponsit              2  5.58  0.394  265.  40.6 -73.9
##  6 NoHo                 53  5.54  0.546  254.  40.7 -74.0
##  7 Sea Gate              7  5.46  1.27   233.  40.6 -74.0
##  8 Midtown            1113  5.44  0.634  230.  40.8 -74.0
##  9 Breezy Point          2  5.40  0.175  221.  40.6 -73.9
## 10 West Village        525  5.38  0.554  217.  40.7 -74.0
## # ... with 207 more rows

Some neighbourhoods have very few properties, so the true average price in those areas will be poorly estimated. From the very small sample that I’ve looked at, it seems that the expensive properties are situated fairly close to one another.

Here is a histogram of the (geometric) average price for each neighbourhood in dollars.

# --- distribution of average prices by neighbourhood ---
prep1DF %>%
  group_by( neighbourhood ) %>%
  summarise( m = mean(price),
             .groups="drop") %>%
  mutate( mp = exp(m) - 1) %>%
  ggplot( aes(x=mp)) +
  geom_histogram( bins=100, fill="steelblue") +
  labs(title="Average price in each neighbourhood",
       x = "Geometric mean price per night (dollars)")

# --- prices by room type -------------------------------
prep1DF %>%
  group_by( room_type ) %>%
  summarise( n = n(),
             m = mean(price),
             s = sd(price),
             gm = exp(m) - 1,
             .groups="drop")

## # A tibble: 3 x 5
##   room_type           n     m     s    gm
##   <chr>           <int> <dbl> <dbl> <dbl>
## 1 Entire home/apt 17806  5.15 0.561 171. 
## 2 Private room    15615  4.31 0.506  73.3
## 3 Shared room       797  3.98 0.647  52.6

# --- price by minimum stay ----------------------------
prep1DF %>%
  ggplot( aes(x=log10(minimum_nights), y=price, colour=room_type)) +
  geom_point() +
  geom_smooth() +
  labs(title = "Price by minimum length of stay",
       y = "log(price+1)") +
  theme( legend.position = c(0.85, 0.85))

# --- days available in a given year ---------------------
prep1DF %>%
  ggplot( aes(x=availability_365, y=price)) +
  geom_point() +
  geom_smooth() +
  labs( title="Price and days available for rent in a year",
       y = "log(price+1)", x="days")

Reviews

The number of reviews tells you something about the number of people who have have rented the accommodation and the date of the last review says something about recent activity, but it is not obvious that these will relate to price.

Number of reviews tells you very little since it will be related to the length of time that the property has been on the market. From the plot it appears that cheaper properties have more reviews.

# --- number of reviews -----------------------------------
prep1DF %>%
  ggplot( aes(x=log10(number_of_reviews+1), y=price)) +
  geom_point() +
  geom_smooth() +
  labs(title="Price and the number of reviews",
       y = "log(price+1)", x="total number of reviews")

7008 properties have never had a review, so I’ve made reviews per month equal to zero when there were no reviews. I’ve also truncated reviews_per_month at 20; some properties had an incredibly high number of reviews per month (more than the number of days in a month).

# --- reviews per month -------------------------------------
prep1DF %>%
  mutate( reviews_per_month = ifelse(
                 number_of_reviews == 0, 0, reviews_per_month),
          reviews_per_month = pmin(20, reviews_per_month) ) %>%
  ggplot( aes(x=reviews_per_month, y=price)) +
  geom_point() +
  geom_smooth() +
  labs(title="Price and reviews per month",
       subtitle="truncated at 20",
       y = "log(price+1)", x="reviews per month")

The date of the last review might say something about recent activity.

# --- date of last review -----------------------------------
prep1DF %>%
  mutate( last_review = as.Date(last_review)) %>%
  ggplot( aes(x=last_review, y=price)) +
  geom_point() +
  geom_smooth() +
  labs(title="Price and the date of the last review",
       subtitle="7008 properties had not been reviewed at all",
       y = "log(price+1)", x="date")

It does not look as though the last review is very predictive of price.

Host

It seems to me that the key host related feature is whether someone is running a business with lots of properties on Airbnb, or they are an individual renting out a single property or room. We can see this on a plot of price against the number of properties a host rents out.

# --- number of properties per host --------------------------
prep1DF %>%
  ggplot( aes(x=calculated_host_listings_count, y=price)) +
  geom_point() +
  geom_smooth() +
  labs(title="Price and number of properties that host rents out",
       y = "log(price+1)", x="number of properties")

The cluster of properties over 300 were all rented out by the same person/company. Let’s look at them

# --- host with over 300 properties -------------------------
prep1DF %>%
  filter( calculated_host_listings_count > 300) %>%
  select( host_id, host_name,  calculated_host_listings_count)

## # A tibble: 229 x 3
##      host_id host_name    calculated_host_listings_count
##        <int> <chr>                                 <int>
##  1 219517861 Sonder (NYC)                            327
##  2 219517861 Sonder (NYC)                            327
##  3 219517861 Sonder (NYC)                            327
##  4 219517861 Sonder (NYC)                            327
##  5 219517861 Sonder (NYC)                            327
##  6 219517861 Sonder (NYC)                            327
##  7 219517861 Sonder (NYC)                            327
##  8 219517861 Sonder (NYC)                            327
##  9 219517861 Sonder (NYC)                            327
## 10 219517861 Sonder (NYC)                            327
## # ... with 219 more rows

So, one company (Sonder) is plotted 327 times, because they have 327 properties.

Let’s plot the hosts rather than the properties and place the number on a log-scale

# --- average price by number of host's properties -----------
prep1DF %>%
  group_by( host_id ) %>%
  summarise( count = mean(calculated_host_listings_count),
             price = mean(price),
             .groups="drop") %>%
  mutate( count = pmin(50, count)) %>%
  ggplot( aes(x=count, y=price)) +
  geom_point() +
  geom_smooth() +
  labs(title="Average price and number of properties",
       subtitle="truncated at 50 properties",
       y = "log(price+1)", x="number of properties")

The dip around 5 properties is interesting and may reflect the division between individuals and companies.

name

name is a misnomer, what we actually have is a brief description of the property.

Here are some examples

# --- example of the name field ----------------------------
prep1DF %>%
  select( name) 

## # A tibble: 34,218 x 1
##    name                                    
##    <chr>                                   
##  1 Cute big one bedroom                    
##  2 Feel like you never leave your home     
##  3 Pristine Lower East Side Sanctuary      
##  4 Luxe, Spacious 2BR 2BA Nr Trains        
##  5 1BD brownstone apt in Fort Greene!      
##  6 LOVELY  LARGE SUNNY ROOM    Sunset Park 
##  7 NE..Comfortable Room All Inclusive      
##  8 Bright, Spacious Studio in Brooklyn     
##  9 Upper west Apt close to Central Pk      
## 10 Artist's Jungle Suite + Private Bathroom
## # ... with 34,208 more rows

There are problems with these descriptions. Some use abbreviations and other don’t, some describe the property and some describe the location and some do neither.

First, I convert everything to lower case, remove the punctuation and the digits and replace a few of the abbreviations

# --- clean the name field ---------------------------------
prep1DF %>%
  mutate( name = tolower(name),
          name = ifelse( is.na(name), "", name),
          name = gsub('[[:punct:] ]+',' ',name),
          name = gsub('[[:digit:] ]+',' ',name),
          name = str_replace(name, "bd ", "bedroom "),
          name = str_replace(name, "br ", "bedroom "),
          name = str_replace(name, "nr ", "near"),
          name = str_replace(name, " rm", " room"),
          name = str_replace(name, "pk", "park"),
          name = str_replace(name, "lux ", "luxury "),
          name = str_replace(name, " ft", " foot"),
          name = str_replace(name, "sq ", "square "),
          name = str_replace(name, " st ", " street "),
          name = str_replace(name, " apt", " apartment")) -> prep2DF

Next, I split the sentences into words (max 12) and make a long column out of all of the words. Of course, some (123) entries have more than 12 words and the extra words will be lost, c’est la vie.

For the first Spiced challenge, s01-e01, I wrote a function that does just what we need, you can find the code in that post. I placed the function in a package called myText and I’ll use it again here.

library(myText)

prep2DF %>%
  select_indicators(col=name, response=price, 
                    nTerms=12, minCount=100, sep=" ") %>%
  filter( str_length(term) > 1 ) %>%
  filter( !(term %in% c("the", "an", "to", "in",
                        "of", "by", "with", "and",
                        "at", "from", "for", "on",
                        "w/")) ) %>%
  filter( count > 100 & p_value < 0.001 ) %>%
  print(n=20) -> wordsDF

## # A tibble: 174 x 5
##    term       count   p_value mWithout mWith
##    <chr>      <int>     <dbl>    <dbl> <dbl>
##  1 room       16343 6.00e-323     4.87  4.59
##  2 private     5100 1.37e-303     4.79  4.44
##  3 cozy        3569 3.95e-157     4.77  4.48
##  4 village     1605 4.11e-156     4.72  5.14
##  5 luxury      1345 8.97e-150     4.72  5.27
##  6 bushwick    1025 2.97e-141     4.75  4.26
##  7 blueground   158 2.01e-129     4.73  5.71
##  8 bedroom     9034 1.32e-115     4.69  4.88
##  9 sonder       299 3.72e-107     4.73  5.43
## 10 bed        10834 1.14e- 97     4.68  4.85
## 11 apartment   8056 2.11e- 88     4.70  4.86
## 12 west        1132 1.94e- 84     4.72  5.13
## 13 art         9981 2.32e- 81     4.69  4.84
## 14 chelsea      501 1.40e- 75     4.73  5.27
## 15 own         2991 2.69e- 75     4.71  4.97
## 16 midtown      914 1.21e- 73     4.73  5.17
## 17 shared       307 3.65e- 68     4.74  4.02
## 18 studio      2917 1.34e- 67     4.72  4.90
## 19 doorman      317 2.79e- 67     4.73  5.37
## 20 gym          316 3.81e- 66     4.73  5.43
## # ... with 154 more rows

So we have 174 potential features in the form of words that occur in the name variable. Some duplicate location information e.g. bushwick, others duplicate host information e.g. sonder or room type information e.g. bedroom. In my judgement it is not worth the effort of trying to edit the words any further and I take the 174 features asis.

Neighbourhoods

I need to sort out the problem of neighbourhoods with very few properties. Such neighbourhoods will produce unreliable estimates in the model.

I will employ a simple clustering algorithm that takes any neighbourhood with fewer that 25 properties and merges it with its closest neighbouring neighbourhood, where distance is calculated from the latitude and longitude.

I checked and found that some of the neighbourhoods in the test data do not appear in the training data and one neighbourhood (Bay Terrace) appears in the test data as Bay Terrace, Staten Island. I reclassify the neighbourhoods based on the combined test and training data in order to ensure that the classifications are consistent.

Here is the code

# --- create DF of unique neighbourhoods ------------
# merge is a flag that I'll use in the subsequent analysis
# it will be set to the row number of the neighbourhood
# that this neighbourhood is merged with
prep2DF %>%
  select( neighbourhood, latitude, longitude) %>%
  bind_rows( testRawDF %>%
               mutate( neighbourhood = str_replace(neighbourhood,
                      "Bay Terrace, Staten Island", "Bay Terrace")) %>%
             select( neighbourhood, latitude, longitude) ) %>%
  group_by( neighbourhood) %>%
  summarise( n = n(),
             lat = mean(latitude),
             lng = mean(longitude),.groups = "drop") %>%
  mutate( merge = 0 ) %>%
  print() -> nbdDF

## # A tibble: 220 x 5
##    neighbourhood         n   lat   lng merge
##    <chr>             <int> <dbl> <dbl> <dbl>
##  1 Allerton             42  40.9 -73.9     0
##  2 Arden Heights         4  40.6 -74.2     0
##  3 Arrochar             21  40.6 -74.1     0
##  4 Arverne              77  40.6 -73.8     0
##  5 Astoria             900  40.8 -73.9     0
##  6 Bath Beach           17  40.6 -74.0     0
##  7 Battery Park City    70  40.7 -74.0     0
##  8 Bay Ridge           141  40.6 -74.0     0
##  9 Bay Terrace           8  40.7 -73.9     0
## 10 Baychester            7  40.9 -73.8     0
## # ... with 210 more rows

Next I write a function that takes a target neighbourhood and merges it with its nearest neighbour

# --- merge_neighbourhoods --------------------------
# function to merge with the nearest neighbour
#    target ... name of neighbourhood to me merged
merge_neighbourhoods <- function(target) {
  # --- row to be merged ----------------------------
  row <- which(nbdDF$neighbourhood == target & nbdDF$merge == 0)
  # --- distances to other neighbourhoods -----------
  distance <- (nbdDF$lat - nbdDF$lat[row])^2 +
    (nbdDF$lng - nbdDF$lng[row])^2 
  # --- make distance large if already merged -------
  distance <- ifelse( nbdDF$merge > 0, 999, distance )
  # --- make distance to itself large ---------------
  distance[row] <- 999
  # --- find the closest ----------------------------
  closest <- which( distance == min(distance) )[1]
  # --- find new total number -----------------------
  nbdDF$n[closest] <- nbdDF$n[row] + nbdDF$n[closest]
  # --- register the identity of the merge ----------
  nbdDF$merge[row] <- closest
  return(nbdDF)
}

Finally I make the merges by working through nbdDF looking for neighbourhoods with fewer than 25 properties.

# --- How many small areas are left unmerged? ------
nSmall <- sum( nbdDF$n < 25 & nbdDF$merge == 0)
# --- merge while small areas remain ---------------
while( nSmall > 0 ) {
  # --- find a target ------------------------------
  nbdDF %>%
    filter( n < 25 & merge == 0 )  %>%
    dplyr::slice(1) %>%
    pull(neighbourhood) -> target
  # --- merge using the function -------------------
  nbdDF <- merge_neighbourhoods( target)
  
  # --- How many small areas are left --------------
  nSmall <- sum( nbdDF$n < 25 & nbdDF$merge == 0)
}

print(nbdDF)

## # A tibble: 220 x 5
##    neighbourhood         n   lat   lng merge
##    <chr>             <int> <dbl> <dbl> <dbl>
##  1 Allerton             42  40.9 -73.9     0
##  2 Arden Heights         4  40.6 -74.2   101
##  3 Arrochar             21  40.6 -74.1    46
##  4 Arverne              77  40.6 -73.8     0
##  5 Astoria             900  40.8 -73.9     0
##  6 Bath Beach           17  40.6 -74.0    17
##  7 Battery Park City    70  40.7 -74.0     0
##  8 Bay Ridge           153  40.6 -74.0     0
##  9 Bay Terrace           8  40.7 -73.9   125
## 10 Baychester            7  40.9 -73.8   154
## # ... with 210 more rows

So Allerton did not need to be merged but Arden Heights got merged with row 101 (Huguenot). I’m afraid that the names of the neighbourhoods mean little to me, but the merges look sensible on Google maps.

Next I reclassify the properties and create a new variable called neighbourhood2

# --- create a second neighbourhood variable ------
prep2DF$neighbourhood2 <- prep2DF$neighbourhood
# --- fill neighbourhood2 with the larger group ---
for( i in 1:nrow(nbdDF)) {
  if( nbdDF$merge[i] > 0 ) {
    from <- nbdDF$neighbourhood[i]
    to   <- nbdDF$neighbourhood[nbdDF$merge[i]]
    j <- which(prep2DF$neighbourhood2==from)
    prep2DF$neighbourhood2[j] <- to
  }
}

We can check that neighbourhood2 looks sensible

# --- look at the new variable -------------------
prep2DF %>%
  group_by( neighbourhood2) %>%
  summarise( n = n(),
             .groups = "drop") 

## # A tibble: 134 x 2
##    neighbourhood2         n
##    <chr>              <int>
##  1 Allerton              35
##  2 Arverne               53
##  3 Astoria              640
##  4 Battery Park City     44
##  5 Bay Ridge            116
##  6 Bayside               36
##  7 Bedford-Stuyvesant  2595
##  8 Bensonhurst           76
##  9 Boerum Hill          132
## 10 Borough Park          90
## # ... with 124 more rows

We have reduced the number of neighbourhoods to 134 and all of them have at least 25 properties in the combined test and training data, though not necessarily in the training data alone.

Save results of the pre-processing

# --- save the cleaned data -----------------------------
trainRawDF %>%
  filter( price != 0 ) %>%
  mutate( price = log(price + 1)) %>%
  clean( wordsDF, nbdDF) %>%
  saveRDS( file.path(home, "data/rData/processed_train.rds"))

testRawDF %>%
  mutate( neighbourhood = str_replace(neighbourhood,
                "Bay Terrace, Staten Island", "Bay Terrace")) %>%
  clean( wordsDF, nbdDF) %>%
  saveRDS( file.path(home, "data/rData/processed_test.rds"))

Submission

Although this model was not tuned, let’s get an idea of performance by making a submission.

I refit the model using all of the training data and run 200 iterations to give the model time to minimise the loss.

# --- place predictors is matrix X -------------------------
trainDF %>%
  select( minimum_nights, reviews_per_month,
          calculated_host_listings_count, availability_365,
          starts_with("N", ignore.case=FALSE), 
          starts_with("R", ignore.case=FALSE), 
          starts_with("X", ignore.case=FALSE)) %>%
  as.matrix() -> X

# --- place response in vector Y ---------------------------
trainDF %>% 
  pull(price) -> Y

# --- fit the model ----------------------------------------
xgmod <- xgboost(data=X, label=Y, nrounds=200, verbose=0)

xgmod$evaluation_log %>% 
  dplyr::slice(1,10,25,50,100,150,200)

##    iter train_rmse
## 1:    1   3.027281
## 2:   10   0.467436
## 3:   25   0.412192
## 4:   50   0.386695
## 5:  100   0.359233
## 6:  150   0.342308
## 7:  200   0.327474

The in-sample RMSE looks impressive but we know that this is misleading.

Now I create the test submission

# --- read clean test data ---------------------------------
testDF <- readRDS( file.path(home, "data/rData/processed_test.rds"))

# --- predictors: note change in column numbers (no price) -
testDF %>%
  select( minimum_nights, reviews_per_month,
          calculated_host_listings_count, availability_365,
          starts_with("N", ignore.case=FALSE), 
          starts_with("R", ignore.case=FALSE), 
          starts_with("X", ignore.case=FALSE)) %>%
  as.matrix() -> XT

# --- make predictions and save ----------------------------
testDF %>%
  mutate( yhat     = predict(xgmod, newdata=XT ) ) %>%
  # --- transform price to dollars ---------------
  mutate( price = exp(yhat) - 1 ) %>%
  select(id, price) %>%
  write.csv( file.path(home, "temp/submission1.csv"),
             row.names=FALSE)

I submitted this file and got RMSE=0.43089, good enough to put us in second place on the leaderboard.

Number of Iterations

Let’s look at the number of iterations; without doubt, this is the most critical hyperparameter.

It is in the nature of gradient boosting that we progress towards the minimum loss in a series of small steps. At each step we fit a tree targeted on the observations that are currently poorly predicted. Too few iterations and we will not minimise the loss, too many and we will start to pick up random blips in the training data that will not be seen in the test data, i.e we will overfit.

The optimum will depend on the complexity of the problem that you are analysing, so it cannot be predicted in advance. The only thing that you can do, is to try a large number iterations and look at the loss at each step.

Here is what happens if we run 1000 iterations.

# --- fit 1000 iterations --------------------------
xgb.train(data=dtrain, 
        watchlist=list(train=dtrain, test=dtest),
        objective="reg:squarederror",
        nrounds=1000, verbose=0) -> xgmod

# --- progress of the fit -------------------------
xgmod$evaluation_log %>% 
  dplyr::slice(1, seq(100, 1000, by=100))

##     iter train_rmse test_rmse
##  1:    1   3.027657  3.029342
##  2:  100   0.349448  0.423110
##  3:  200   0.315679  0.422535
##  4:  300   0.285025  0.425237
##  5:  400   0.262911  0.428523
##  6:  500   0.241819  0.432565
##  7:  600   0.225361  0.434975
##  8:  700   0.208922  0.437628
##  9:  800   0.194874  0.440214
## 10:  900   0.182076  0.442238
## 11: 1000   0.170291  0.444098

I will wrap my plotting code in a function, so that I can reuse it later.

# --- plot RMSE by iteration ----------------------
xgmod$evaluation_log %>%
    ggplot( aes(x=iter, y=train_rmse)) +
    geom_line(colour="blue") +
    geom_line( aes(y=test_rmse), colour="red") +
    scale_y_continuous( limits = c(0, 1), breaks=seq(0, 1, by=0.1)) +
    labs(x="Iteration", y="RMSE", 
         title="In-sample and out-of-sample RMSE")

The in-sample RMSE keeps improving. It will eventually plateau at the minimum possible loss for this tree-based model, but the complexity of the trees is such that by the time that that happens we will be picking up every meaningless blip in the data.

For the Airbnb data, the out-of-sample RMSE is minimised at about 200 iterations and after that it increases slowly. The optimum happens to coincide with the choice that I made for my submission. You will just have to believe me when I say that this was a lucky guess. However, the rise in the out-of-sample RMSE is so slow that we lose very little by having a few extra iterations. Generally, a few too many iterations is better than not enough.

For these data, the actual minimum RMSE occurs at iteration 170

# --- iteration with the minumum RMSE -----------------------
xgmod$evaluation_log %>%
  summarise( minIter = which(test_rmse == min(test_rmse))[1] )

##   minIter
## 1     170

This, of course, is only the result for this particular validation sample. Let’s try some different random estimation/validation splits

# --- optimum performance for different data splits --------
set.seed(4256)
# --- variables for saving the results
minIteration <- minRMSE <- rep(0,5)
for( i in 1:5 ) {
  # --- random split --------------------------
  split <- sample(1:nrow(trainDF), size=10000, replace=FALSE)
  # --- training X and Y matrices -------------
  trainDF[-split, ] %>%
    select( minimum_nights, reviews_per_month,
            calculated_host_listings_count, availability_365,
            starts_with("N", ignore.case=FALSE), 
            starts_with("R", ignore.case=FALSE), 
            starts_with("X", ignore.case=FALSE)) %>%
    as.matrix() -> Xi
  trainDF[-split, ] %>% 
    pull(price) -> Yi
  dtraini <- xgb.DMatrix(data = Xi, label = Yi)
  # --- validation X and Y matrices ----------
  trainDF[ split, ] %>%
    select( minimum_nights, reviews_per_month,
            calculated_host_listings_count, availability_365,
            starts_with("N", ignore.case=FALSE), 
            starts_with("R", ignore.case=FALSE), 
            starts_with("X", ignore.case=FALSE)) %>%
    as.matrix() -> XVi
  trainDF[ split, ] %>% 
    pull(price) -> YVi
  dtesti <- xgb.DMatrix(data = XVi, label=YVi)
  # --- fit the model ------------------------
  xgb.train(data=dtraini, 
            watchlist=list(train=dtraini, test=dtesti),
            objective="reg:squarederror",
            nrounds=300, verbose=0) -> xgmod
  # --- find optimum number of iterations ----
  xgmod$evaluation_log %>%
    summarise( minIter = which(test_rmse == min(test_rmse))[1] ) %>%
    pull(minIter) -> minIteration[i]
  # --- RMSE at the optimum ------------------
  minRMSE[i] <- xgmod$evaluation_log$test_rmse[minIteration[i]]
}
print(minIteration)

## [1] 186 168 115 183 215

print(minRMSE)

## [1] 0.417998 0.431867 0.427539 0.419533 0.426129

A single validation sample only gives an approximation to the optimum number of iterations. In this analysis, the validation sample has size 10,000, which it is quite large. Many problems will have a much smaller validation sample and therefore be far less stable.

The random splits suggest that the required number of iterations for this problem is something between 115 and 215.

My submission was based on 200 iterations of the algorithm, but looking at these results I might have chosen 150 iterations instead. If I had, perhaps my model would have topped the leaderboard or come 5th instead of 2nd.

The five random splits give minimum RMSEs between 0.418 and 0.432. This tells us that there is no point in getting excited about selections of hyperparameters that improve a single measure of the RMSE by say, 0.010.

A cross-validated RMSE would be equivalent to averaging over several validation samples so would be more stable, but it would still have an error associated with it.

Other hyperparameters

There are too many hyperparameters to run through them all here and this post is already very long. I’ll look at some other hyperparameters when I next use xgboost.