https://christophm.github.io/interpretable-ml-book/bike-data.html

Data from this company: https://www.capitalbikeshare.com/

instant: record index
dteday : date
season : season (1:winter, 2:spring, 3:summer, 4:autumn)
yr : year (0: 2011, 1:2012)
mnth : month ( 1 to 12)
hr : hour (0 to 23)
holiday : weather day is holiday or not (extracted from [Web Link])
weekday : day of the week
workingday : if day is neither weekend nor holiday is 1, otherwise is 0.
weathersit :
1: Clear, Few clouds, Partly cloudy, Partly cloudy
2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist
3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds
4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog
temp : Normalized temperature in Celsius (0-1). The values are derived via (t-t_min)/(t_max-t_min), t_min=-8, t_max=+39 (only in hourly scale)
atemp: Normalized “feeling” temperature in Celsius (0-1). The values are derived via (t-t_min)/(t_max-t_min), t_min=-16, t_max=+50 (only in hourly scale)
hum: Normalized humidity. The values are divided to 100 (max)
windspeed: Normalized wind speed. The values are divided to 67 (max)
casual: count of casual users
registered: count of registered users
cnt: count of total rental bikes including both casual and registered

library(randomForest)

## randomForest 4.7-1.2

## Type rfNews() to see new features/changes/bug fixes.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine()  masks randomForest::combine()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ ggplot2::margin() masks randomForest::margin()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Mostly we have done classification. We will use regression, and its evaluation metrics.

https://r-charts.com/colors/ R colour picker!

bike_day <- read_csv('./Bike-Sharing-Dataset/day.csv')

## Rows: 731 Columns: 16
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (15): instant, season, yr, mnth, holiday, weekday, workingday, weathers...
## date  (1): dteday
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

bike_hour <- read_csv('./Bike-Sharing-Dataset/hour.csv')

## Rows: 17379 Columns: 17
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (16): instant, season, yr, mnth, hr, holiday, weekday, workingday, weat...
## date  (1): dteday
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# have to remove registered and casual otherwise we get serious leakage (will just do total counts)
# season needs to be a factor

bike_day <- bike_day %>% 
  dplyr::select(-registered, -casual) %>%        
  mutate(season = factor(season, 
                         levels = c(1,2,3,4), 
                         labels = c('Winter', 'Spring', 'Summer', 'Autumn'))) %>% 
  mutate(weathersit = factor(weathersit)) %>%                                                 # convert weather situation to factor
  mutate(days = as.numeric(dteday - as.Date("2011-01-01"))) %>%                               # count number of days
  select(count = cnt, season, holiday, yr, days, workingday, weathersit, temp, hum, windspeed) 

ggplot(bike_day, aes(x = count)) + 
  geom_histogram(binwidth = 400, fill = "cadetblue4", color = "black")                       # plot hist of bike counts/day

Create dummy variables from season and weather situation, then create a matrix, and plot

When you see high negative correlations between dummy variables of the same original categorical variable (e.g., weathersit_1 and weathersit_2), it indicates mutual exclusivity. There are some weathersit_3 in there, but not many

library(corrplot)

## corrplot 0.95 loaded

library(fastDummies)

mcor <- bike_day %>% 
  dummy_cols(select_columns = c('season', 'weathersit')) %>% 
  select_if(is.numeric) %>% 
  cor
corrplot(mcor, method = 'color', diag = F, tl.col = 'grey40',
         order = 'hclust', addrect = 3)

xs_ - data to predict; ys_ - the values we are predicting. Test = 1; train != 1

set.seed(1123)

bike_day$train <- sample(1:10, nrow(bike_day), replace = T)

xs_train <- bike_day %>% 
  filter(train != 1) %>% 
  select(season:windspeed)
ys_train <- bike_day %>% 
  filter(train != 1) %>% 
  pull(count)

xs_test <- bike_day %>% 
  filter(train == 1) %>% 
  select(season:windspeed)
ys_test <- bike_day %>% 
  filter(train == 1) %>% 
  pull(count)

Linear Regression

# Simple linear regression 

ggplot(bike_day, aes(x=days, y=count)) +
  geom_point() +
  geom_smooth(method="lm")

## `geom_smooth()` using formula = 'y ~ x'

simple_linear_regression <- lm(count ~ days, data=bike_day)
simple_linear_regression   #inspect model

## 
## Call:
## lm(formula = count ~ days, data = bike_day)
## 
## Coefficients:
## (Intercept)         days  
##    2398.730        5.769

# fit model using all the data

bike_lm <- lm(count ~ ., 
              data = filter(bike_day, train != 1) %>% select(count:windspeed))
summary(bike_lm)

## 
## Call:
## lm(formula = count ~ ., data = filter(bike_day, train != 1) %>% 
##     select(count:windspeed))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3695.7  -402.3    68.4   497.4  3450.1 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1578.9208   251.4734   6.279 6.36e-10 ***
## seasonSpring  1232.9517   128.5306   9.593  < 2e-16 ***
## seasonSummer   986.3988   184.8416   5.336 1.32e-07 ***
## seasonAutumn  1742.8844   172.1555  10.124  < 2e-16 ***
## holiday       -576.2673   196.8245  -2.928  0.00354 ** 
## yr            2207.6143   227.7387   9.694  < 2e-16 ***
## days            -0.5475     0.5898  -0.928  0.35365    
## workingday     150.2442    75.1043   2.000  0.04588 *  
## weathersit2   -424.2398    90.3603  -4.695 3.27e-06 ***
## weathersit3  -1902.4227   218.8778  -8.692  < 2e-16 ***
## temp          5058.8139   337.4811  14.990  < 2e-16 ***
## hum          -1372.4762   331.4684  -4.141 3.94e-05 ***
## windspeed    -2698.0109   468.6402  -5.757 1.34e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 853.3 on 631 degrees of freedom
## Multiple R-squared:  0.8154, Adjusted R-squared:  0.8118 
## F-statistic: 232.2 on 12 and 631 DF,  p-value: < 2.2e-16

What do the outputs mean?

Adjusted R² = 0.812:

Indicates that approximately 81.2% of the variability in daily bike rentals (count) can be explained by the model predictors, which is very strong performance for a real-world dataset. Residual standard error (853.3):

Typical prediction errors are around +/- 853 bikes per day.

Temperature, weather conditions, and seasonality (especially Autumn) significantly influence bike rentals. Year-on-year growth indicates increased popularity of bikes. Negative effects of humidity, windspeed, and holidays. The slight daily trend (days) is not statistically significant.

Fit to test and train sets Calculate R2 metric - to indicate how well the independent variables explain the variance of the dependent variable. It is expressed as a percentage, with higher percentages indicating a better fit of the model to the data.

Each point compares the actual observed bike counts to the model predictions.

High and similar R2 in training and test sets indicates strong and stable performance. A significant drop in test R2 would signal overfitting.

d_train <- tibble(predicted = predict(bike_lm), 
                  observed = ys_train, 
                  set = 'train')

d_test <- tibble(predicted = predict(bike_lm, newdata = xs_test), 
                 observed = ys_test, 
                 set = 'test')

d <- bind_rows(d_train, d_test) %>% 
  mutate(set = factor(set, levels = c('train', 'test')))

# Calculate R² for training data
ssres_train <- sum((d_train$observed - d_train$predicted)^2)   #discrepancy between predicted and observed data
sstot_train <- sum((d_train$observed - mean(d_train$observed))^2)    #total variance in the observed data
r2_train <- 100 * (1 - ssres_train/sstot_train)

# Calculate R² for test data
ssres_test <- sum((d_test$observed - d_test$predicted)^2)
sstot_test <- sum((d_test$observed - mean(d_test$observed))^2)
r2_test <- 100 * (1 - ssres_test / sstot_test)

# Add R2 annotation coordinates
r2 <- tibble(
  x   = c(quantile(d$observed, 0.05), quantile(d$observed, 0.05)),
  y   = c(max(d$predicted), max(d$predicted) - 500),
  lbl = c(
    paste0("R2 (train) = ", round(r2_train, 2), "%"),
    paste("R2 (test) =", round(r2_test, 2), "%")
  ),
  set = factor(c('train', 'test'), levels = c('train', 'test'))
)

# Plot observed vs predicted
ggplot(d, aes(x = observed, y = predicted)) +
  geom_point(color = 'springgreen4', alpha = 0.5) +
  geom_smooth(method = 'lm', color = 'grey10') +
  geom_text(data = r2, aes(x, y, label = lbl), color = 'black', size = 4, hjust = 0) +
  facet_wrap(~set) +
  labs(x = "Observed Count", y = "Predicted Count",
       title = "Linear regression: Observed vs. Predicted Bike Counts",
       subtitle = "Comparing training and test set predictions") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

This is really easy to interpret.

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

ggcoef(bike_lm)

Random Forest for regression

Need to make sure importance = T and localImp = T are both enabled for calculating variable importance

set.seed(5813)
bike_rf <- randomForest(x = xs_train, y = ys_train, 
                        xtest = xs_test, ytest = ys_test,
                        importance = T, localImp = T, keep.forest = T)
print(bike_rf)

## 
## Call:
##  randomForest(x = xs_train, y = ys_train, xtest = xs_test, ytest = ys_test,      importance = T, localImp = T, keep.forest = T) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 440719.8
##                     % Var explained: 88.59
##                        Test set MSE: 207756.6
##                     % Var explained: 92.87

d_train <- tibble(predicted = predict(bike_rf), observed = ys_train, set = 'train')
d_test <- tibble(predicted = predict(bike_rf, newdata = xs_test), observed = ys_test, set = 'test')
d <- bind_rows(d_train, d_test) %>% 
  mutate(set = factor(set, levels = c('train', 'test')))

# Calculate R² for training data
ssres_train <- sum((d_train$observed - d_train$predicted)^2)   #discrepancy between predicted and observed data
sstot_train <- sum((d_train$observed - mean(d_train$observed))^2)    #total variance in the observed data
r2_train <- 100 * (1 - ssres_train/sstot_train)

# Calculate R² for test data
ssres_test <- sum((d_test$observed - d_test$predicted)^2)
sstot_test <- sum((d_test$observed - mean(d_test$observed))^2)
r2_test <- 100 * (1 - ssres_test / sstot_test)

# Add R2 annotation coordinates
r2 <- tibble(
  x   = c(quantile(d$observed, 0.05), quantile(d$observed, 0.05)),
  y   = c(max(d$predicted), max(d$predicted) - 500),
  lbl = c(
    paste0("R2 (train) = ", round(r2_train, 2), "%"),
    paste("R2 (test) =", round(r2_test, 2), "%")
  ),
  set = factor(c('train', 'test'), levels = c('train', 'test'))
)

ggplot(d, aes(x = observed, y = predicted)) +
  geom_point(color = 'darkseagreen4', alpha = 0.5) +
  geom_smooth(method = 'lm', color = 'grey10') +
  geom_text(data = r2, aes(x, y, label = lbl), color = 'black', size = 4, hjust = 0) +
  facet_wrap(~set) +
  labs(x = "Observed Count", y = "Predicted Count",
       title = "Random Forest: Observed vs. Predicted Bike Counts",
       subtitle = "Comparing training and test set predictions") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

XGBoost for regression

XGBoost (Extreme Gradient Boosting) is a powerful, tree-based ensemble method that iteratively builds multiple weak predictive models to create a stronger predictive model.

library(xgboost)

## 
## Attaching package: 'xgboost'

## The following object is masked from 'package:dplyr':
## 
##     slice

set.seed(21)
xgb_train <- xs_train %>% 
  dummy_cols(select_columns = c('season', 'weathersit')) %>% 
  select_if(is.numeric) %>% 
  as.matrix %>% 
  xgb.DMatrix(label = ys_train)

xgb_test <- xs_test %>% 
  dummy_cols(select_columns = c('season', 'weathersit')) %>% 
  select_if(is.numeric) %>% 
  as.matrix %>% 
  xgb.DMatrix(label = ys_test)

wl <- list(train = xgb_train, eval = xgb_test)

param <- list(max_depth = 2, eta = 0.2, verbose = 0, nthread = 1, gamma = 0,
              subsample = 0.6, min_child_weight = 2,
              objective = "reg:squarederror", eval_metric = "rmse")

bike_xgb <- xgb.train(param, data = xgb_train, watchlist = wl, nrounds = 200)

## [17:53:06] WARNING: src/learner.cc:767: 
## Parameters: { "verbose" } are not used.
## 
## [1]  train-rmse:4000.517337  eval-rmse:3855.403782 
## [2]  train-rmse:3275.426055  eval-rmse:3136.945451 
## [3]  train-rmse:2700.194489  eval-rmse:2568.293041 
## [4]  train-rmse:2241.666322  eval-rmse:2113.831071 
## [5]  train-rmse:1882.073422  eval-rmse:1754.352200 
## [6]  train-rmse:1591.830922  eval-rmse:1472.676955 
## [7]  train-rmse:1383.392032  eval-rmse:1281.057337 
## [8]  train-rmse:1212.645194  eval-rmse:1104.215222 
## [9]  train-rmse:1091.612446  eval-rmse:978.811837 
## [10] train-rmse:994.534216   eval-rmse:875.424933 
## [11] train-rmse:915.374653   eval-rmse:785.852818 
## [12] train-rmse:863.528413   eval-rmse:741.754129 
## [13] train-rmse:825.409271   eval-rmse:700.960775 
## [14] train-rmse:794.867112   eval-rmse:661.580114 
## [15] train-rmse:770.693687   eval-rmse:636.887450 
## [16] train-rmse:752.036894   eval-rmse:620.310632 
## [17] train-rmse:732.857237   eval-rmse:608.756593 
## [18] train-rmse:718.233709   eval-rmse:595.045954 
## [19] train-rmse:705.054252   eval-rmse:581.373365 
## [20] train-rmse:691.572103   eval-rmse:564.809276 
## [21] train-rmse:679.514139   eval-rmse:549.416468 
## [22] train-rmse:671.848629   eval-rmse:544.576247 
## [23] train-rmse:659.643673   eval-rmse:546.666242 
## [24] train-rmse:649.994604   eval-rmse:545.277975 
## [25] train-rmse:643.525918   eval-rmse:540.190898 
## [26] train-rmse:636.008347   eval-rmse:537.219853 
## [27] train-rmse:631.362480   eval-rmse:532.554819 
## [28] train-rmse:625.727083   eval-rmse:529.723228 
## [29] train-rmse:619.725374   eval-rmse:526.167054 
## [30] train-rmse:611.927234   eval-rmse:523.590609 
## [31] train-rmse:607.111156   eval-rmse:521.007212 
## [32] train-rmse:600.377714   eval-rmse:522.297452 
## [33] train-rmse:595.699787   eval-rmse:527.841983 
## [34] train-rmse:592.727433   eval-rmse:526.442363 
## [35] train-rmse:590.268573   eval-rmse:526.890748 
## [36] train-rmse:585.754605   eval-rmse:517.213895 
## [37] train-rmse:580.413517   eval-rmse:517.524122 
## [38] train-rmse:576.994691   eval-rmse:513.774987 
## [39] train-rmse:573.992829   eval-rmse:515.223932 
## [40] train-rmse:571.539274   eval-rmse:512.950770 
## [41] train-rmse:567.830665   eval-rmse:514.003692 
## [42] train-rmse:565.521320   eval-rmse:511.000659 
## [43] train-rmse:561.208585   eval-rmse:513.380999 
## [44] train-rmse:559.205495   eval-rmse:510.006225 
## [45] train-rmse:556.342344   eval-rmse:511.291711 
## [46] train-rmse:553.197676   eval-rmse:507.164870 
## [47] train-rmse:551.067384   eval-rmse:508.063383 
## [48] train-rmse:548.222717   eval-rmse:510.179073 
## [49] train-rmse:546.597095   eval-rmse:509.340544 
## [50] train-rmse:545.259899   eval-rmse:507.940024 
## [51] train-rmse:543.706238   eval-rmse:506.834845 
## [52] train-rmse:541.897583   eval-rmse:508.151842 
## [53] train-rmse:537.420828   eval-rmse:504.563698 
## [54] train-rmse:534.378459   eval-rmse:502.301156 
## [55] train-rmse:532.632905   eval-rmse:500.286914 
## [56] train-rmse:531.058164   eval-rmse:500.565173 
## [57] train-rmse:530.091695   eval-rmse:499.780702 
## [58] train-rmse:528.922636   eval-rmse:497.613297 
## [59] train-rmse:527.102985   eval-rmse:499.387353 
## [60] train-rmse:526.520744   eval-rmse:498.861298 
## [61] train-rmse:524.251277   eval-rmse:497.018630 
## [62] train-rmse:523.039742   eval-rmse:497.479337 
## [63] train-rmse:521.009042   eval-rmse:494.426885 
## [64] train-rmse:518.809218   eval-rmse:496.457045 
## [65] train-rmse:516.434181   eval-rmse:493.515744 
## [66] train-rmse:515.518709   eval-rmse:492.846757 
## [67] train-rmse:513.011140   eval-rmse:492.350046 
## [68] train-rmse:511.111661   eval-rmse:489.747242 
## [69] train-rmse:509.763299   eval-rmse:489.681018 
## [70] train-rmse:508.005028   eval-rmse:485.004529 
## [71] train-rmse:507.642217   eval-rmse:483.483009 
## [72] train-rmse:505.652062   eval-rmse:483.058101 
## [73] train-rmse:502.838600   eval-rmse:479.116776 
## [74] train-rmse:501.853880   eval-rmse:478.537248 
## [75] train-rmse:499.168388   eval-rmse:480.309072 
## [76] train-rmse:494.950226   eval-rmse:477.010928 
## [77] train-rmse:493.396504   eval-rmse:477.277386 
## [78] train-rmse:491.314297   eval-rmse:480.597697 
## [79] train-rmse:490.518953   eval-rmse:480.352242 
## [80] train-rmse:489.508460   eval-rmse:480.684900 
## [81] train-rmse:486.963415   eval-rmse:481.948369 
## [82] train-rmse:486.260998   eval-rmse:481.192520 
## [83] train-rmse:484.437094   eval-rmse:480.552516 
## [84] train-rmse:482.402816   eval-rmse:475.585829 
## [85] train-rmse:481.592364   eval-rmse:479.397852 
## [86] train-rmse:480.815353   eval-rmse:477.245984 
## [87] train-rmse:479.318902   eval-rmse:476.450456 
## [88] train-rmse:477.118560   eval-rmse:474.645283 
## [89] train-rmse:476.380529   eval-rmse:473.318065 
## [90] train-rmse:475.950145   eval-rmse:470.510122 
## [91] train-rmse:475.681715   eval-rmse:469.593586 
## [92] train-rmse:474.520841   eval-rmse:468.388491 
## [93] train-rmse:472.970862   eval-rmse:471.769860 
## [94] train-rmse:470.483042   eval-rmse:471.175769 
## [95] train-rmse:469.608084   eval-rmse:471.901552 
## [96] train-rmse:467.024317   eval-rmse:471.879536 
## [97] train-rmse:466.031454   eval-rmse:473.447554 
## [98] train-rmse:464.187232   eval-rmse:471.951701 
## [99] train-rmse:463.483610   eval-rmse:471.512792 
## [100]    train-rmse:462.643861   eval-rmse:469.282955 
## [101]    train-rmse:460.602518   eval-rmse:469.316733 
## [102]    train-rmse:459.840642   eval-rmse:464.862199 
## [103]    train-rmse:458.538138   eval-rmse:464.616870 
## [104]    train-rmse:457.903606   eval-rmse:465.075458 
## [105]    train-rmse:455.106541   eval-rmse:462.239356 
## [106]    train-rmse:453.831239   eval-rmse:459.569699 
## [107]    train-rmse:452.067645   eval-rmse:465.243414 
## [108]    train-rmse:451.292178   eval-rmse:463.364892 
## [109]    train-rmse:450.505899   eval-rmse:462.637717 
## [110]    train-rmse:448.639314   eval-rmse:464.771777 
## [111]    train-rmse:447.086876   eval-rmse:464.280816 
## [112]    train-rmse:445.431110   eval-rmse:462.111204 
## [113]    train-rmse:444.960872   eval-rmse:461.106021 
## [114]    train-rmse:444.221621   eval-rmse:460.510170 
## [115]    train-rmse:443.462692   eval-rmse:460.815877 
## [116]    train-rmse:441.855617   eval-rmse:460.875310 
## [117]    train-rmse:440.486479   eval-rmse:469.254014 
## [118]    train-rmse:439.700668   eval-rmse:469.344693 
## [119]    train-rmse:439.133036   eval-rmse:470.006718 
## [120]    train-rmse:437.442051   eval-rmse:466.933628 
## [121]    train-rmse:436.807171   eval-rmse:470.617022 
## [122]    train-rmse:436.197210   eval-rmse:467.624041 
## [123]    train-rmse:435.087233   eval-rmse:465.425807 
## [124]    train-rmse:434.526000   eval-rmse:465.794685 
## [125]    train-rmse:432.910973   eval-rmse:464.177353 
## [126]    train-rmse:431.974343   eval-rmse:461.345122 
## [127]    train-rmse:431.635994   eval-rmse:465.410033 
## [128]    train-rmse:431.160050   eval-rmse:465.780363 
## [129]    train-rmse:429.512656   eval-rmse:466.252120 
## [130]    train-rmse:428.914124   eval-rmse:465.886454 
## [131]    train-rmse:427.638646   eval-rmse:468.488841 
## [132]    train-rmse:427.221285   eval-rmse:471.082272 
## [133]    train-rmse:425.368939   eval-rmse:469.537366 
## [134]    train-rmse:424.122014   eval-rmse:469.629758 
## [135]    train-rmse:422.875803   eval-rmse:465.961674 
## [136]    train-rmse:421.996793   eval-rmse:467.441784 
## [137]    train-rmse:421.099111   eval-rmse:472.290355 
## [138]    train-rmse:420.638818   eval-rmse:473.266655 
## [139]    train-rmse:419.962873   eval-rmse:474.793927 
## [140]    train-rmse:419.084747   eval-rmse:470.050893 
## [141]    train-rmse:419.090617   eval-rmse:471.413328 
## [142]    train-rmse:418.886157   eval-rmse:467.759134 
## [143]    train-rmse:418.486806   eval-rmse:467.170522 
## [144]    train-rmse:417.880029   eval-rmse:466.815110 
## [145]    train-rmse:417.965917   eval-rmse:468.901582 
## [146]    train-rmse:417.530607   eval-rmse:467.655607 
## [147]    train-rmse:416.436857   eval-rmse:464.244836 
## [148]    train-rmse:415.346875   eval-rmse:462.251605 
## [149]    train-rmse:414.345238   eval-rmse:462.946753 
## [150]    train-rmse:413.945478   eval-rmse:461.464911 
## [151]    train-rmse:412.629563   eval-rmse:463.833928 
## [152]    train-rmse:411.182855   eval-rmse:459.966438 
## [153]    train-rmse:410.282951   eval-rmse:458.822014 
## [154]    train-rmse:407.656354   eval-rmse:457.995978 
## [155]    train-rmse:407.072151   eval-rmse:455.995722 
## [156]    train-rmse:405.799114   eval-rmse:456.136619 
## [157]    train-rmse:404.850220   eval-rmse:459.181244 
## [158]    train-rmse:404.490024   eval-rmse:460.491808 
## [159]    train-rmse:403.755405   eval-rmse:458.945488 
## [160]    train-rmse:403.599610   eval-rmse:455.050476 
## [161]    train-rmse:402.326936   eval-rmse:453.023926 
## [162]    train-rmse:402.028937   eval-rmse:451.882093 
## [163]    train-rmse:401.382881   eval-rmse:453.565756 
## [164]    train-rmse:401.439630   eval-rmse:452.692159 
## [165]    train-rmse:400.961409   eval-rmse:452.880255 
## [166]    train-rmse:400.000300   eval-rmse:452.481728 
## [167]    train-rmse:398.666285   eval-rmse:452.773818 
## [168]    train-rmse:397.974994   eval-rmse:456.190994 
## [169]    train-rmse:397.568744   eval-rmse:456.954788 
## [170]    train-rmse:396.539878   eval-rmse:456.653592 
## [171]    train-rmse:395.244207   eval-rmse:452.396998 
## [172]    train-rmse:394.873342   eval-rmse:451.121022 
## [173]    train-rmse:393.949696   eval-rmse:450.767528 
## [174]    train-rmse:393.833471   eval-rmse:450.040574 
## [175]    train-rmse:392.955855   eval-rmse:451.035356 
## [176]    train-rmse:391.034692   eval-rmse:453.015038 
## [177]    train-rmse:390.091064   eval-rmse:452.312737 
## [178]    train-rmse:389.765522   eval-rmse:456.871824 
## [179]    train-rmse:388.573118   eval-rmse:457.080656 
## [180]    train-rmse:387.431505   eval-rmse:455.889257 
## [181]    train-rmse:386.890690   eval-rmse:456.549255 
## [182]    train-rmse:386.753477   eval-rmse:457.898874 
## [183]    train-rmse:386.349542   eval-rmse:458.383267 
## [184]    train-rmse:384.230594   eval-rmse:454.658805 
## [185]    train-rmse:383.628532   eval-rmse:454.811702 
## [186]    train-rmse:383.317044   eval-rmse:454.637728 
## [187]    train-rmse:382.456520   eval-rmse:454.864979 
## [188]    train-rmse:381.441904   eval-rmse:453.803657 
## [189]    train-rmse:381.515919   eval-rmse:454.930186 
## [190]    train-rmse:380.658963   eval-rmse:451.225198 
## [191]    train-rmse:380.166198   eval-rmse:451.344093 
## [192]    train-rmse:379.307168   eval-rmse:451.286994 
## [193]    train-rmse:378.066386   eval-rmse:450.024330 
## [194]    train-rmse:377.757636   eval-rmse:449.537267 
## [195]    train-rmse:377.124710   eval-rmse:449.579836 
## [196]    train-rmse:377.072358   eval-rmse:447.232580 
## [197]    train-rmse:376.772553   eval-rmse:448.592397 
## [198]    train-rmse:375.121269   eval-rmse:447.318155 
## [199]    train-rmse:374.187661   eval-rmse:447.599269 
## [200]    train-rmse:374.155352   eval-rmse:448.221270

d_train <- tibble(predicted = predict(bike_lm), 
                  observed = ys_train, 
                  set = 'train')

d_test <- tibble(predicted = predict(bike_lm, newdata = xs_test), 
                 observed = ys_test, 
                 set = 'test')

d <- bind_rows(d_train, d_test) %>% 
  mutate(set = factor(set, levels = c('train', 'test')))

# Calculate R² for training data
ssres_train <- sum((d_train$observed - d_train$predicted)^2)   #discrepancy between predicted and observed data
sstot_train <- sum((d_train$observed - mean(d_train$observed))^2)    #total variance in the observed data
r2_train <- 100 * (1 - ssres_train/sstot_train)

# Calculate R² for test data
ssres_test <- sum((d_test$observed - d_test$predicted)^2)
sstot_test <- sum((d_test$observed - mean(d_test$observed))^2)
r2_test <- 100 * (1 - ssres_test / sstot_test)

# Add R2 annotation coordinates
r2 <- tibble(
  x   = c(quantile(d$observed, 0.05), quantile(d$observed, 0.05)),
  y   = c(max(d$predicted), max(d$predicted) - 500),
  lbl = c(
    paste0("R2 (train) = ", round(r2_train, 2), "%"),
    paste("R2 (test) =", round(r2_test, 2), "%")
  ),
  set = factor(c('train', 'test'), levels = c('train', 'test'))
)

# Plot observed vs predicted
ggplot(d, aes(x = observed, y = predicted)) +
  geom_point(color = 'rosybrown', alpha = 0.5) +
  geom_smooth(method = 'lm', color = 'grey10') +
  geom_text(data = r2, aes(x, y, label = lbl), color = 'black', size = 4, hjust = 0) +
  facet_wrap(~set) +
  labs(x = "Observed Count", y = "Predicted Count",
       title = "XGBoost: Observed vs. Predicted Bike Counts",
       subtitle = "Comparing training and test set predictions") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

Explain the results

Two resources on interpretable machine learning (how to understand your models and how they are making their predictions):

http://uc-r.github.io/iml-pkg https://christophm.github.io/interpretable-ml-book/

# The interpretable machine learning package
library(iml)

Need to create a function that outputs a value with the model and new data. Doesn’t matter what model you do, the prediction function needs to work this way.

xgb_pred <- function(model, newdata) {
  xs <- newdata %>% 
    dummy_cols(select_columns = c('season', 'weathersit')) %>% 
    select_if(is.numeric) %>% 
    as.matrix %>% 
    xgb.DMatrix()
  predict(model, newdata = xs)
}

Linear regression is easy to see how predictions are made. xgboost is more opaque

bike_xgb_predictor <- Predictor$new(bike_xgb,
                                    data = xs_train,
                                    y = ys_train,
                                    predict.fun = xgb_pred)

Explain away

Importance

We can look at feature importance of XGBoost model. Why are days so important here, but insignificant for the linear regression model?

imp <- FeatureImp$new(bike_xgb_predictor, loss = 'rmse')
imp$plot()

Partial dependence - what has the modelled learned about temperature?

# The influence of varying temperature
eff <- FeatureEffect$new(bike_xgb_predictor, method = 'pdp', 
                         feature = 'temp', grid.size = 100)
eff$plot()

eff$plot() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_line(color = 'grey50') +
  geom_smooth(method = 'loess', se = F, color = 'slateblue') +
  labs(x = 'Temperature\n(0 = -8C, 1 = 39C)') +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

Explain a random day using Shapley values

When applied to machine learning models, it helps in understanding how each feature contributes to the prediction for a specific instance.

A Shapley value for a feature indicates how much that feature pushes the prediction higher or lower relative to the model’s baseline prediction.

Positive SHAP value: Feature increased the prediction. Negative SHAP value: Feature decreased the prediction. Magnitude: Indicates strength of influence on the prediction.

i <- sample(1:nrow(xs_train), 1)

x.interest <- xs_train[i,]
print(x.interest)

## # A tibble: 1 × 9
##   season holiday    yr  days workingday weathersit  temp   hum windspeed
##   <fct>    <dbl> <dbl> <dbl>      <dbl> <fct>      <dbl> <dbl>     <dbl>
## 1 Summer       0     0   223          1 1          0.708 0.415     0.126

shapley <- Shapley$new(bike_xgb_predictor, x.interest = x.interest)
shapley$plot()

Week 09: Biking

Alison

2025-03-13