--- title: "Comparing some of R's random forest packages" author: "D G Rossiter" date: "`r Sys.Date()`" output: word_document: toc: yes html_document: toc: yes toc_float: yes theme: lumen code_folding: show number_section: yes fig_height: 6 fig_width: 6 fig_align: center --- # Purpose We compare several R packages that build random forests: an older package `randomForest` and a much faster implementation, `ranger`, and the `spmodel` package that also models the spatial structure of the random forest residuals and interpolates these to improve predictions. We hope the results are similar. Note that due to randomization there will always be differences, although not between runs of this script because of the use of `set.seed`. This document also explores Quantile Regression Forests option in `ranger` to assess uncertainty of predictions, and local measures of variable importance including Shapley values. We also have some fun with interpretable machine learning... just because we can! # Dataset We use as an example the small "Meuse (Maas) River soil pollution" dataset. See `help(meuse, package = "sp")` for details and references. This dataset is provided with the `sp` "classes and methods for spatial data: points, lines, polygons and grids" package. Athough `sp` has been superseded by `sf`, it is still available and has this famous test dataset. You do not need to load `sp`, just have it installed in order to locate the `sp::meuse` dataset. ```{r} data(meuse, package = "sp") str(meuse) ``` We try to model one of the heavy metal concentration (zinc, Zn) from all possibly relevant predictors: * `ffreq` flooding frequency (3 classes) * `dist.m` distance from the Meuse river (m) * `elev` elevation above mean sea level (m) * `soil` soil type (3 classes) * `lime` whether agricultural lime was applied to the field (yes/no) We also include locations, because this may account for local effects, by building "boxes" around sets of locations. * `x`, `y` coordinates in the RDH (Dutch) grid (m), EPSG code 28992 Make a log10-transformed copy of the Zn metal concentration to obtain a somewhat balanced distribution of the target variable: ```{r fig.width=8, fig.height=4} meuse$logZn <- log10(meuse$zinc) # 锌 hist(meuse$logZn, main = ""); rug(meuse$logZn) ``` There is a good range of values, in a bimodal distribution: low and high values, with fewere medium values. There is still some right skew towards the highest values. # Random forest with `randomForest` This was (one of?) the first R packages to implement Breiman's random forest algorithm.^[Breiman, L. (2001), Random Forests, Machine Learning 45(1), 5-32]. Packages used in this section: `randomForest` for random forests, `randomForestExplainer` for some diagnostics, and `ggplot2` for graphics. ```{r} library(ggplot2, warn.conflicts = FALSE, verbose = FALSE) library(randomForest, warn.conflicts = FALSE, verbose = FALSE) library(randomForestExplainer, warn.conflicts = FALSE, verbose = FALSE) ``` ## Build the forest First build the forest, using the `randomForest` function. Use `set.seed` for reproducibility, so that your results (if you render this script or run this chunk by itself) will be the same as these. The `mtry` optional argument gives the number of variables randomly sampled as candidates at each split. By default this is $\lfloor{p/3}\rfloor$ where $p$ is the number of predictors, here $7$, so $\lfloor{7/3}\rfloor = 2$. We increase this to 3 to get better trees but still include weak predictors; this is matched for `ranger` (below). ```{r compute.rf} set.seed(314) m.lzn.rf <- randomForest(logZn ~ ffreq + dist.m + elev + soil + lime + x + y, data=meuse, # eight permutations per tree to estimate importance importance = TRUE, nperm = 8, na.action = na.omit, mtry = 3) print(m.lzn.rf) ``` The structure of this object: ```{r} str(m.lzn.rf) ``` How successful was the forest? Compute the RMSE and $R^2$ of the fits, based on the out-of-bag (OOB) cross-validation. These are estimates of the predictive power of the fitted model. ```{r} m.lzn.rf.resid <- (meuse$logZn - m.lzn.rf$predicted) print(sqrt(mse <- mean(m.lzn.rf.resid^2))) # RMSE print(rsq <- (1 - mse/var(meuse$logZn))) # pseudeo-R^2 ``` The model is successful overall, with a low root mean squared error (RMSE) `r round(sqrt(mse),4)` and high $R^2=$ `r round(rsq,4)`, these both from the fit, not cross-validation. Compare the RMSE to the mean of the target variable, `r round(mean(meuse$logZn), 4)` Display the cross-validation error rate against the number of trees, to see how many trees were needed for a stable result: ```{r fig.width=8, fig.height=4} plot(m.lzn.rf) ``` Only about 300 trees were needed in this case. ## Variable importance Examine the variable importance. There are two kinds. From the help text: 1. "The first measure is computed from permuting OOB data: For each tree, the prediction error on the out-of-bag portion of the data is recorded (error rate for classification, MSE for regression). Then the same is done after permuting each predictor variable. The difference between the two are then averaged over all trees, and normalized by the standard deviation of the differences." 2. "The second measure is the total decrease in node impurities from splitting on the variable, averaged over all trees. For classification, the node impurity is measured by the Gini index. For regression, it is measured by residual sum of squares." Here we have a regression, so Type 2 is the change in RSS due to the permutation. We compare both: ```{r fig.width=10, fig.height=5} print(data.frame(randomForest::importance(m.lzn.rf, type=1), randomForest::importance(m.lzn.rf, type=2))) par(mfrow = c(1, 2)) varImpPlot(m.lzn.rf, type=1, main = "Importance: permutation") varImpPlot(m.lzn.rf, type=2, main = "Importance: node impurity") par(mfrow = c(1, 1)) ``` The ranks are similar but not identical. After the first two `dist.m` and `elev` (obviously the most important, since the heavy metal is mostly from flooding of the Meuse/Maas River) the others are not in the same order. There is likely substitution of partially correlated predictors. ## Goodness-of-fit Predict back onto the known points, and evaluate the goodness-of-fit: ```{r fig.width=5, fig.height=5} p.rf <- predict(m.lzn.rf, newdata=meuse) summary(r.rpp <- meuse$logZn - p.rf) print(paste("RMSE:", round(rmse.rf <- sqrt(sum(r.rpp^2)/length(r.rpp)), 4))) plot(meuse$logZn ~ p.rf, asp=1, pch=20, xlab="fitted", ylab="actual", xlim=c(2,3.3), ylim=c(2,3.3), main="log10(Zn), Meuse topsoils, Random Forest") grid(); abline(0,1) ``` Quite a good internal fit, but this is because we use all the training points to make each prediction. ## Out-of bag cross-validation A more realistic view of the predictive power of the model is the out-of-bag cross-validation. This summarizes the predictions at training points *not* used in a set of trees -- recall that in the random forest each tree uses sampling with replacement, and some training points are not included. We get the OOB X-validation with `predict` with no dataset specified, i.e., it's the one used to build the model.. ```{r fig.width=5, fig.height=5} p.rf.oob <- predict(m.lzn.rf) summary(r.rpp.oob <- meuse$logZn - p.rf.oob) (rmse.oob <- sqrt(sum(r.rpp.oob^2)/length(r.rpp.oob))) plot(meuse$logZn ~ p.rf.oob, asp=1, pch=20, xlab="Out-of-bag cross-validation estimates", ylab="actual", xlim=c(2,3.3), ylim=c(2,3.3), main="log10(Zn), Meuse topsoils, Random Forest") grid() abline(0,1) ``` There is quite a bit more spread than the internal fit, as expected. The summary statistics are a bit more than doubled. ## Sensitivity Because of randomization, each run of the random forest will be different. To evaluate this, compute the RF several times and collect statistics. This also allows a speed comparison with other RF packages. ```{r sensitive.rf, fig.width=8, fig.height=4} set.seed(314) n <- 48 # number of simulations # data frame to collect results rf.stats <- data.frame(rep=1:10, rsq=as.numeric(NA), mse=as.numeric(NA)) t.rf <- system.time( for (i in 1:n) { model.rf <- randomForest(logZn ~ ffreq + x + y + dist.m + elev + soil + lime, data=meuse, importance=T, na.action=na.omit, mtry=5) summary(model.rf$rsq) summary(model.rf$mse) rf.stats[i, "rsq"] <- median(summary(model.rf$rsq)) rf.stats[i, "mse"] <- median(summary(model.rf$mse)) } ) summary(rf.stats[,2:3]) hist(rf.stats[,"rsq"], xlab="RandomForest R^2", breaks = 16, main = "Frequency of fits (R^2)") rug(rf.stats[,"rsq"]) hist(rf.stats[,"mse"], xlab="RandomForest RMSE", breaks = 16, main = "Frequency of OOB acuracy (RMSE)") rug(rf.stats[,"mse"]) ``` This is fairly stable. It would be much less so with subsets (e.g., randomly remove 10% of training points). # Random forest with `ranger` Packages used in this section: `ranger` for random forests and `vip` for variable importance. (Note that `vip` can be used for many model types, see https://koalaverse.github.io/vip/articles/vip.html). The `mtry` optional argument gives the number of variables randomly sampled as candidates at each split. By default for `ranger` (unlike `randomForest`) this is $\lfloor \sqrt{p} \rfloor$ where $p$ is the number of predictors, here $5$, so $\lfloor \sqrt{p} \rfloor = 2$. We want to try all the variables at each split for more accurate splitting. We increase this to 3 to get better trees but still include weak predictors; this is matched for `randomForest` (above). ## Build the forest First build the forest, using the `ranger` function. Ask for the `permutation` measure of variable importance, to match `randomForest`. The variable importance measures are of several types: 1. "impurity": variance of the responses for regression (as here) 2. "impurity_corrected": "The 'impurity_corrected' importance measure is unbiased in terms of the number of categories and category frequencies" -- not relevant for regression forests 3. "permutation". For this one can specify `scale.permutation.importance = TRUE`, this should match the `randomForest` concept. Compute two ways, for the two kinds of importance. Use the same random seed, the forests will then be identical, but the importance measures will be different. Permutation: ```{r ranger.perm} library(ranger, warn.conflicts = FALSE, verbose = FALSE) set.seed(314) m.lzn.ra <- ranger(logZn ~ ffreq + x + y + dist.m + elev + soil + lime, data = meuse, importance = 'permutation', scale.permutation.importance = TRUE, mtry = 3) print(m.lzn.ra) ``` Impurity: ```{r ranger.impurity} set.seed(314) m.lzn.ra.i <- ranger(logZn ~ ffreq + x + y + dist.m + elev + soil + lime, data = meuse, importance = 'impurity', mtry = 3) print(m.lzn.ra.i) ``` ## Goodness-of-fit: Predict with fitted model, at all observations. For this we need to specify a `data=` argument. ```{r} p.ra <- predict(m.lzn.ra, data=meuse) str(p.ra) summary(r.rap <- meuse$logZn - p.ra$predictions) (rmse.ra <- sqrt(sum(r.rap^2)/length(r.rap))) ``` Compare the fits RMSE with that computed by `RandomForest`: ```{r} c(rmse.ra, rmse.rf) ``` There is a slight difference. Plot the fits: ```{r fig.width=10, fig.height=5} par(mfrow=c(1,2)) plot(meuse$logZn ~ p.ra$predictions, asp=1, pch=20, xlab="fitted", ylab="log10(Zn)", xlim=c(2,3.3), ylim=c(2,3.3), main="ranger") grid(); abline(0,1) plot(meuse$logZn ~ p.rf, asp=1, pch=20, xlab="fitted", ylab="log10(Zn)", xlim=c(2,3.3), ylim=c(2,3.3), main="randomForest") grid(); abline(0,1) par(mfrow=c(1,1)) ``` Almost identical, including which points are poorly-predicted. ## Out-of-bag cross-validation The default model already has the OOB predictions stored in it. Compare to `RandomForest` results. ```{r oob.compare, fig.width=10, fig.height=5} summary(m.lzn.ra$predictions) summary(p.rf.oob) summary(m.lzn.ra$predictions - p.rf.oob) # difference ``` Overall, `ranger` has a slightly higher median OOB prediction `randomForest`. Plot the OOB cross-validation predictions: ```{r oob, fig.width=10, fig.height=5} par(mfrow=c(1,2)) plot(meuse$logZn ~ m.lzn.ra$predictions, asp=1, pch=20, ylab="log10(Zn)", xlab="OOB X-validation estimates", xlim=c(2,3.3), ylim=c(2,3.3), main="ranger") abline(0,1); grid() plot(meuse$logZn ~ p.rf.oob, asp=1, pch=20, xlab="OOB X-validation estimates", ylab="log10(Zn)", xlim=c(2,3.3), ylim=c(2,3.3), main="randomForest") grid(); abline(0,1) par(mfrow=c(1,1)) ``` These are very similar. As with the fits, the same points are poorly-predicted by both models. ## Variable importance: First, for permutation, the sums and the individual importances: ```{r} sum(ranger::importance(m.lzn.ra)) sum(randomForest::importance(m.lzn.rf)[,1]) cbind(ranger = ranger::importance(m.lzn.ra), rF = randomForest::importance(m.lzn.rf)[,1]) ``` Second, for impurity, the sums and the individual importances:: ```{r} sum(ranger::importance(m.lzn.ra.i)) sum(randomForest::importance(m.lzn.rf)[,2]) cbind(ranger =ranger::importance(m.lzn.ra.i), rF = randomForest::importance(m.lzn.rf)[,2]) ``` There is quite some difference between the packages:`randomForest` gives much more importance to the geographic coordinates and much less to the distance to river and elevation. The sums of variable importance are slightly different. ## Sensitivity Compute the `ranger` RF several times and collect statistics, also the timing: ```{r sensitive.ra, fig.width=8, fig.height=4} n <- 48 ra.stats <- data.frame(rep=1:10, rsq=as.numeric(NA), mse=as.numeric(NA)) t.ra <- system.time( for (i in 1:n) { model.ra <- ranger(logZn ~ ffreq + x + y + dist.m + elev + soil + lime, data=meuse, importance="none", mtry=5, write.forest = FALSE) ra.stats[i, "mse"] <- model.ra$prediction.error ra.stats[i, "rsq"] <- model.ra$r.squared } ) summary(ra.stats[,2:3]) hist(ra.stats[,"rsq"], xlab="ranger R^2", breaks = 16, main = "Frequency of fits (R^2)") rug(ra.stats[,"rsq"]) hist(ra.stats[,"mse"], xlab="ranger RMSE", breaks = 16, main = "Frequency of OOB accuracy (RMSE)") rug(ra.stats[,"mse"]) ``` Again, not a large variation. Compare the timings: ```{r} print(t.rf) print(t.ra) print(round(t.rf/t.ra, 1)) ``` `randomForest` is almost 6x slower than `ranger` in this test. Compare the sensitivity statistics to those from `RandomForest`: ```{r} summary(ra.stats[,2:3]) summary(rf.stats[,2:3]) tmp <- round(data.frame(sd(ra.stats[,2]), sd(rf.stats[,2]), sd(ra.stats[,3]), sd(rf.stats[,3])),6) names(tmp) <- c("s.d. R^2 ranger", "s.d. R^2 randomForest","s.d. MSE ranger","s.d. MSE randomForest") print(tmp) ``` No clear pattern as to which is more sensitive. # Spatial random forests with `splmRF` Thanks to Michael McManus (EPA) for this. Here we try a spatially-explicit random forest plus local model using the experimental (as of mid 2025) `spmodel` package, also from the EPA.^[Dumelle, M., Higham, M., & Hoef, J. M. V. (2023). spmodel: Spatial statistical modeling and prediction in R. PLOS ONE, 18(3), e0282524. https://doi.org/10.1371/journal.pone.0282524] First, make the object spatially-explicit in the `sf` "Simple Features" data structure. Note that we specify `remove = FALSE` to keep `x` and `y` also in the data frame, to use in the random forest model. ```{r} library(sf, warn.conflicts = FALSE, verbose = FALSE) meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), remove = FALSE) summary(meuse_sf) ``` Plot the target variable. ```{r} ggplot(data = meuse_sf) + geom_sf(aes(size = logZn, col = logZn), pch=1) ``` The `spmodel` package includes `splmRF` for random forest spatial residual models, using random forest to fit the mean and a spatial linear model (e.g., kriging) to fit the residuals. "A random forest model is fit to the mean portion of the model specified by formula using `ranger::ranger()`. Residuals are computed and used as the response variable in an intercept-only spatial linear model fit using `splm().`" So in this case we can use the coordinates twice: (1) in the random forest model as with the other methods, to make "boxes", and (2) for spatial interpolation of the residuals from the random forest fit. To compare with the other RF methods, we set `mtry = 3`. The spatial covariance is fit with an exponential function, a usual conservative choice. ```{r} set.seed(314) library(spmodel, warn.conflicts = FALSE, verbose = FALSE) m.lzn.splmRF <- splmRF(logZn ~ ffreq + x + y + dist.m + elev + soil + lime, data = meuse_sf, spcov_type = "exponential", importance = 'permutation', estmethod = "sv-wls", weights = "pairs-invrd", scale.permutation.importance = TRUE, mtry = 3) summary(m.lzn.splmRF) ``` The fitted model object has the call, the results from calling `ranger` to fit the random forest, and the parameters of the fitted residual model, i.e., spatial covariance structure. ```{r} print(m.lzn.splmRF$call) print(m.lzn.splmRF$ranger) print(m.lzn.ra) print(m.lzn.splmRF$splm) ``` The `ranger` portion of the model is identical to calling `ranger` directly, due to the same model form, parameters and seed. There is residual spatial correlation, with an effective range of `r round(m.lzn.splmRF$splm$coefficients$spcov["range"]*3)` meters. This is the `range` coefficient of the exponential model, multiplied by 3 to reach 95% of the total sill, as is conventional. The coefficient `de` is the dependent error variance, and `ie` the independent error variance. These correspond to the partial sill and nugget variance in a variogram model. So the nugget-to-total sill ratio is `ie/(ie+de)`: ```{r} round(m.lzn.splmRF$splm$coefficients$spcov["ie"]/(m.lzn.splmRF$splm$coefficients$spcov["ie"]+m.lzn.splmRF$splm$coefficients$spcov["de"]),2) ``` This is a high nugget, around 40% of the residual variance after the random forest model. Compare the variable importance to `ranger`: ```{r} sort(m.lzn.splmRF$ranger$variable.importance, decreasing = TRUE) sort(m.lzn.ra$variable.importance, decreasing = TRUE) ``` These are identical, because of the same random seed. Out-of-bag predictions from the random forest, comparing to `ranger`: ```{r} summary(m.lzn.splmRF$ranger$predictions) summary(m.lzn.ra$predictions) ``` Again, identical. Compare model performance: ```{r} print(m.lzn.splmRF$ranger$r.squared) print(m.lzn.ra$r.squared) ``` Again, identical. Compare the predictions at the training points: ```{r} summary(p.splmRF <- predict(m.lzn.splmRF, newdata = meuse_sf)) summary(p.ra.p <- predict(m.lzn.ra, data = meuse)$predictions) summary(p.rf.p <- predict(m.lzn.rf, newdata = meuse)) ``` The predictions from `randomForest` and `ranger` are very close but not identical. Including the interpolation of residuals with `splmRF` has widened the predictions and slightly lowered many of them. We can see the extremes by adding a `rug` plot to the histograms: ```{r fig.width=8, fig.height=4} hist(p.splmRF, main = "splmRF", xlab = "log10(Zn)") rug(p.splmRF) hist(p.ra.p, main = "ranger", xlab = "log10(Zn)") rug(p.ra.p) hist(p.rf.p, main = "randomForest", xlab = "log10(Zn)") rug(p.rf.p) ``` # Spatial structure of `ranger` random forest residuals We can try to duplicate `splmRF` by modelling the spatial structure of the `ranger` residuals with a variogram model and then interpolating these to modify the ranger predictions. We use the `gstat` package for this. ```{r gstat.load} library(gstat, warn.conflicts = FALSE, verbose = FALSE) ``` Compute the `ranger` residuals and compute their empirical variogram using the cutoff of 1/2 of the diagonal of the bounding box (to match the default in `splmRF`, personal communication Mike Dunelle), and 15 bins: ```{r gstat.v, fig.width=8, fig.height=6} summary(meuse_sf$r.ra <- (meuse_sf$logZn - p.ra.p)) bbox <- sf::st_bbox(meuse_sf) corners <- st_sfc( st_point(c(bbox["xmin"], bbox["ymin"])), st_point(c(bbox["xmax"], bbox["ymax"])), crs = st_crs(meuse_sf) ) (diag <- st_distance(corners[1], corners[2])) v <- variogram(r.ra ~ 1, data = meuse_sf, cutoff = diag/2) plot(v, plot.numbers = TRUE) ``` Model this with an exponential model (as used by `splmRF` in our example) and fit the model with `gstat::fit.variogram` using its default _ad-hoc_ weighted least-squares optimization; we specified that above for `splmRF`. ```{r gstat.vmf,fig.width=8, fig.height=6} (vm <- fit.variogram(v, vgm(psill = 0.002, "Exp", range = 800/3, nugget = 0.003))) plot(v, plot.numbers = TRUE, model = vm) vm[1,"psill"]/(sum(vm[,"psill"])) ``` This has a nugget/total sill ratio of `r round(vm[1,"psill"]/(sum(vm[,"psill"])), 2)`, somewhat higher than estimated by `splmRF`. The effective range is `r round(vm[2,"range"]*3)` m, quite a bit longer than that estimated by `splmRF`. Compare to the `splmRF`-fitted variogram: ```{r} print(vm) print(m.lzn.splmRF$splm$coefficients$spcov) ``` Both the partial sill and nugget are much higher for `splmRF`, it seems the semivariances are much higher overall. Perhaps not divided by 1/2? Compare the proportional nuggets: ```{r} (vm[1,"psill"]/sum(vm[,"psill"])) (round(m.lzn.splmRF$splm$coefficients$spcov["ie"]/(m.lzn.splmRF$splm$coefficients$spcov["ie"]+m.lzn.splmRF$splm$coefficients$spcov["de"]),2)) ``` Not identical, despite using the same variogram cutoff, covariance function, and fitting method. The `gstat` estimate is longer range and a higher proportional nugget. Maybe because of the binning used by `splmRF`? ## Kriging the residuals Krige the residuals to the observation points by LOOCV (otherwise the predictions would be identical to the actual values, since kriging is an exact interpolator): ```{r} kcv <- krige.cv(r.ra ~ 1, nfold = nrow(meuse_sf), locations = meuse_sf, model = vm) head(kcv) summary(kcv$var1.pred) summary(meuse_sf$r.ra) ``` Kriging has greatly smoothed out the residuals. ## Combined prediction: RF + kriging Add these to the `ranger` predictions, and compare to `splmRF1`: ```{r} summary(p.ra$predictions + kcv$var1.pred) summary(p.splmRF) ``` These are not the same. Try with an exponential variogram model identical to that fit by `splmRF`: ```{r} (coef <- m.lzn.splmRF$splm$coefficients$spcov[1:3]) (vm.splmRF <- vgm(psill = coef["de"], model = "Exp", range = coef["range"], nugget = coef["ie"])) kcv <- krige.cv(r.ra ~ 1, nfold = nrow(meuse_sf), locations = meuse_sf, model = vm.splmRF) head(kcv) summary(kcv$var1.pred) summary(meuse_sf$r.ra) ``` These values of the `ranger` residuals are not identical. Why? Maybe `splmRF` kriges to the points and so gets the exact `ranger` residual? ```{r} summary(p.ra$predictions + meuse_sf$r.ra) summary(p.splmRF) ``` No, that's not the reason. Although the predictions are close: ```{r} summary((p.ra$predictions + meuse_sf$r.ra)-p.splmRF) ``` Maybe there is a difference in the radius used for the kriging step? # Quantile Regression Forest (QRF) This is an extension of `ranger` which produces estimates of each quantile of the predictive distribution for each prediction. The method is explained by Meinshausen.^[Meinshausen, N. (2006). Quantile regression forests. Journal of Machine Learning Research, 7, 983-999, DOI:10.5555/1248547.1248582] To compute these, use `ranger`, and specify `quantreg = TRUE`. Alsokeep the in-bag fits to use in the quantile predictions, with `keep.inbag = TRUE`. ```{r qrf} set.seed(314) m.lzn.qrf <- ranger(logZn ~ ffreq + x + y + dist.m + elev + soil + lime, data = meuse, quantreg = TRUE, importance = 'permutation', keep.inbag=TRUE, # needed for QRF scale.permutation.importance = TRUE, mtry = 3) pred.qrf <- predict(m.lzn.qrf, type = "quantiles", # default is c(0.1, 0.5, 0.9) quantiles = c(0.1, 0.25, 0.5, 0.75, 0.9)) summary(pred.qrf$predictions) ``` This shows the range of predictions from the set of trees. There is quite some variability between quantiles. For example, here is the Q9-Q1 range: ```{r} summary(pred.qrf$predictions[, "quantile= 0.9"] - pred.qrf$predictions[, "quantile= 0.1"]) ``` Most are on the order of 0.3 - 0.5 log10(Zn). The 0.5 (median) quantile can be compared to the single prediction from `ranger`, although the QRF prediction is a median, not a mean: ```{r} summary(pred.qrf$predictions[, "quantile= 0.5"]) summary(p.ra.p) ``` The means and the Q1/Q3 of the prediction means and medians are close but not identical. # Local measures of variable importance Both `randomForest` and `ranger` provide local measures of importance, i.e., how much each variable influences an individual prediction. ## `randomForest` "The 'local' (or casewise) variable importance for a random regression forest is the average increase in squared OOB residuals when the variable is permuted, for each observation separately." This is computed if `localImp = TRUE` in the call to `randomForest`. ```{r rf.local} set.seed(314) m.lzn.rf.l <- randomForest(logZn ~ ffreq + x + y + dist.m + elev + soil + lime, data=meuse, localImp = TRUE, nperm = 3, na.action = na.omit, mtry = 3) dim(m.lzn.rf.l$localImportance) row.names(m.lzn.rf.l$localImportance) ``` There is one importance measure for each variable for each observation. View the local importance for the first few observation points: ```{r table.rf.local} m.lzn.rf.l$localImportance[, 1:6] ``` We see a large difference in importance of `dist.m`, very important for points 1 and 2, and much less so for points 3--6. ```{r} meuse[1:6, "dist.m"] ``` These two points are much closer to the river. ## `ranger` For this package we use the `local.importance` option: "Calculate and return local importance values as in (Breiman 2001)^[Breiman, L. (2001). Random forests. Mach Learn, 45:5-32. doi: 10.1023/A:1010933404324]. Only applicable if importance is set to `permutation`". This is the same permutation measure as for the overall variable importance, but computed for each observation. Refit the model with this option, using the same random seed so that the overall model fit and importance match the first section above. The model fit will be identical. ```{r ranger.perm.l} set.seed(314) m.lzn.ra.l <- ranger(logZn ~ ffreq + x + y + dist.m + elev + soil + lime, data = meuse, importance = 'permutation', local.importance = TRUE, scale.permutation.importance = TRUE, mtry = 3) names(m.lzn.ra.l) dim(m.lzn.ra.l$variable.importance.local) ``` The `variable.importance.local` field in the fitted model results has the importance of each variable for each of the predictions. Show the local importances for the first few observation points, compared the global importance. These both measure the increase in OOB RMSE when the predictor is permuted. ```{r} m.lzn.ra.l$variable.importance.local[1:6,] ranger::importance(m.lzn.ra.l) ``` We can see quite a range of importances for the predictors, depending on the point. As with `randomForest`, the predictor `dist.m` is much more imporant for points 1 and 2 than for 3-6. Show the range of local variable importance for each predictor: ```{r} summary(m.lzn.ra.l$variable.importance.local[,1:7]) ``` There is a large range in the importance, even for `dist.m` and `elev`. Notice that the overall importance of these is between the 3rd quartile and maximum, i.e., they are quite influential in less than 1/4 of cases. What is the relation of the importance with the factor itself? Show a scatterplot of the imporatance for `dist.m`, with a horizontal line at the overall importance. ```{r plot.ra.importance.local.elev, fig.width=8, fig.height=6} plot(meuse$dist.m, m.lzn.ra.l$variable.importance.local[,"dist.m"], xlab = "elevation", ylab = "importance of distance to river", pch = 20) abline(h = ranger::importance(m.lzn.ra.l)["dist.m"], col = "red") ``` Quite important for the points near the river, but not for those further away. What about for `elev`? ```{r plot.ra.importance.local.dist.m, fig.width=8, fig.height=6} plot(meuse$elev, m.lzn.ra.l$variable.importance.local[,"elev"], xlab = "dist.m", ylab = "importance of elevation m.a.s.l", pch = 20) abline(h = ranger::importance(m.lzn.ra.l)["elev"], col = "red") ``` Weak relation to no relation. What about flood frequency class? ```{r plot.ra.importance.local.ffreq, fig.width=8, fig.height=6} plot(meuse$ffreq, m.lzn.ra.l$variable.importance.local[,"ffreq"], xlab = "ffreq", ylab = "importance of `ffreq`", pch = 20) abline(h = ranger::importance(m.lzn.ra.l)["ffreq"], col = "red") ``` Again, important for a few points near the river. # Interpretable Machine Learning We would like to look inside the "black box" of the random forest models to see if it can give insights into its predictions. Christoph Molnar has written an online text^[Molnar, C. (2022). Interpretable Machine Learning. https://christophm.github.io/interpretable-ml-book/] explaining this. One set of interpretations is provided by the `iml` package. ```{r iml.load} library(iml, warn.conflicts = FALSE, verbose = FALSE) ``` We use this to interpret one of the random forests built above. ## Create an `iml::Predictor` object The interpretation methods in the `iml` package need the machine learning model to be wrapped in a `Predictor` object." This "... holds any machine learning model (`mlr`, `caret`, `randomForest`, ...) and the data to be used for analyzing the model. Create a wrapper function for `predict.ranger` that matches the format `iml` expects. This is because `predict.ranger` does not follow the typical `predict` argument naming convention. ```{r} # Define a prediction function compatible with iml predict_iml <- function(model, newdata, type = "response", ...) { preds <- predict(model, data = newdata, type = type)$predictions if (type == "response") { return(preds) } else { return(preds[, 2]) # For probability, not used here } } ``` (The `ranger` package defines a method called `predict.ranger` for objects of class `ranger`, but it does not export a function called `predict`. Instead, `predict.ranger` will be dispatched from generic `predict` automatically if the object is of type `ranger`. So don't use `ranger::predict`, just `predict`.) Now use it for the model fitted by `ranger`, i.e., make a `Predictor` object. Note that the `predict.function` argument is needed if the `model` is not from either the `mlr` or `caret` packages. ```{r iml.predictor} vars <- c("ffreq","x","y","dist.m","elev","soil","lime") X <- meuse[, vars] predictor.ra <- iml::Predictor$new(model = m.lzn.ra, data = X, y = meuse$logZn, predict.function = predict_iml) str(predictor.ra) ``` Also make a `Predictor` object for the `ranger` QRF model: ```{r} predictor.qrf <- iml::Predictor$new(model = m.lzn.qrf, data = X, y = meuse$logZn, predict.function = predict_iml) ``` ## Feature importance Feature importance, based on `mae` (mean absolute error) or `mse` (mean squared error). Plot the median and the distribution (5--95% quantiles over all the cases). ```{r iml.plot.featureImp, fig.width=8, fig.height=6} imp.mae <- iml::FeatureImp$new(predictor.ra, loss = "mae") imp.mse <- iml::FeatureImp$new(predictor.ra, loss = "mse") plot(imp.mae) plot(imp.mse) print(imp.mae) print(imp.mse) ``` This shows the range of local effects for each predictor, for the given indicator. Predictor `dist.m` has quite variable local effects, whereas `x` and `y` have similar (and smaller) effects at all points. Another way to see local effects is with the `FeatureEffect` object: 1. accumulated local effect (ALE); "ALE plots are a faster and unbiased alternative to partial dependence plots (PDPs)."^[https://christophm.github.io/interpretable-ml-book/ale.html] 2. Partial Dependence Plots (PDP): this shows the prediction with other factors kept at their medians (?).^[https://christophm.github.io/interpretable-ml-book/pdp.html] 3. Individual Conditional Expectation (ICE) curves -- one per observation "one line per instance that shows how the instance’s prediction changes when a feature changes".^[see https://christophm.github.io/interpretable-ml-book/ice.html]. This is a PDP plot per observation. See with `plot.FeatureEffect`. Create the object with `$new`, specifying a method, then plot. ```{r iml.plot.featureEffect.elev, fig.width=8, fig.height=6} ale <- FeatureEffect$new(predictor.ra, feature = "elev") plot(.value ~ elev, data = ale$results, type = "b") pdp <- FeatureEffect$new(predictor.ra, feature = "elev", method = "pdp") plot(.value ~ elev, data = pdp$results, type = "b") ice <- FeatureEffect$new(predictor.ra, feature = "elev", method = "ice") plot(.value ~ elev, data = ice$results, pch=20, cex = 0.3) ``` In the ICE graph some interesting behaviour in the 7-8.5 m elevation range. But in general this follows the average PDP plot, we don't have any really unusual observations in terms of their dependence on elevation. See this for distance: ```{r iml.plot.featureEffect.dist.m, fig.width=8, fig.height=6} ice.dist.m <- FeatureEffect$new(predictor.ra, feature = "dist.m", method = "ice") plot(.value ~ dist.m, data = ice.dist.m$results, pch=20, cex = 0.3) ``` There is an interesting contrasting behaviour from 0 to 50 m: some with high values at 0 increase substantially at 50, while most decrease, as expected. ### Interactions "The interaction measure regards how much of the variance of f(x) is explained by the interaction. The measure is between 0 (no interaction) and 1 (= 100% of variance of f(x) due to interactions). For each feature, we measure how much they interact with any other feature."^[https://christophm.github.io/interpretable-ml-book/interaction.html] "When features interact with each other in a prediction model, the prediction cannot be expressed as the sum of the feature effects, because the effect of one feature depends on the value of the other feature." ```{r plot.iml.interact, fig.width=8, fig.height=6} interact <- Interaction$new(predictor.ra) interact$results ``` This shows that the predictors interact a fair amount -- none are independent. The obvious example is distance and elevation, which we can see in a 2-way interaction plot: Look at the 2-way interactions with a predictor, e.g., elevation: ```{r plot.iml.interact.elev} interact.elev <- Interaction$new(predictor.ra, feature = "elev") interact.elev$results ``` This shows strong interactions between elevation and the others, especially `dist.m`, in prediction $\log_{10}\mathrm{Zn}$. ## Explain single predictions with a local model This uses the `glmnet` "Lasso and Elastic-Net Regularized Generalized Linear Models" package, also the `gower` "Gower's Distance" package, based on Gower (1971).^[Gower, John C. "A general coefficient of similarity and some of its properties." Biometrics (1971): 857-871] "The ... model fitted by `LocalModel` is a linear regression model and the data points are weighted by how close they are to the data point for w[h]ich we want to explain the prediction." Here 'close' means in multivariate feature space, hence the use of Gower's distance. From `?LocalModel`: "A weighted glm is fitted with the machine learning model prediction as target. Data points are weighted by their proximity to the instance to be explained, using the Gower proximity measure. L1-regularization is used to make the results sparse." (hence the use of `glmnet`). So this is a linear regression model, therefore we can interpret its fitted coefficients to see local effects. However the model only uses "close-by" points in feature space. Look at this for the first point, which happens to be close to the river and with a high level of Zn. ```{r lime.explain.single, fig.width=8, fig.height=4} meuse[1,] lime.explain <- iml::LocalModel$new(predictor.ra, x.interest = X[1, ]) lime.explain$results plot(lime.explain) ``` The local model at this point predicts a lower value than actual prediction based on the full model. The main reason is the elevation of the surrounding points: these differ from the target point but affect the prediction. Look at this for one of the furthest points from the river: ```{r lime.explain.single.maxdist, fig.width=8, fig.height=4} ix.maxdist <- which.max(meuse$dist.m) meuse[ix.maxdist, ] lime.explain <- LocalModel$new(predictor.ra, x.interest = X[ix.maxdist, ]) lime.explain$results plot(lime.explain) ``` Here both distance and elevation have a strong effect but the predictions are the same. ## Shapley values The local importance metrics explored in the previous sections do not account for the full set of interactions at a data point. The Shapley method corrects this. It is available in several packages, including `iml`. An excellent intuitive explanation is given by Molnar in his "Interpretable Machine Learning" on-line text^[https://christophm.github.io/interpretable-ml-book/shapley.html] "[A] method from coalitional game theory named Shapley value. Assume that for one data point, the feature values play a game together, in which they get the prediction as a payout. The Shapley value tells us how to fairly distribute the payout among the feature values. "The 'game' is the prediction task for a single instance of the dataset. The 'gain' is the actual prediction for this instance minus the average prediction for all instances. The 'players' are the feature values of the instance that collaborate to receive the gain (= predict a certain value). "The Shapley value is the *average marginal contribution of a feature value across all possible coalitions*. "The Shapley value is the only explanation method with a solid theory. The axioms – efficiency, symmetry, dummy, additivity – give the explanation a reasonable foundation."^[https://christophm.github.io/interpretable-ml-book/shapley.html] The Shapley value is a solution for computing feature contributions for single predictions for any machine learning model. It is defined via a value function ${val}$ of players in $S$. The Shapley value of a feature value is its contribution to the payout, weighted and summed over all possible feature value combinations: \[\phi_j(val)=\sum_{S\subseteq\{1,\ldots,p\} \backslash \{j\}}\frac{|S|!\left(p-|S|-1\right)!}{p!}\left(val\left(S\cup\{j\}\right)-val(S)\right)\] where $S$ is a subset of the features used in the model, $x$ is the vector of feature values of the instance to be explained and $p$ the number of features. $val_x(S)$ is the prediction for feature values in set $S$ that are marginalized over features that are not included in set $S$: \[val_{x}(S)=\int\hat{f}(x_{1},\ldots,x_{p})d\mathbb{P}_{x\notin{}S}-E_X(\hat{f}(X))\] Note that the Shapley values are dependent on the reference dataset which is used to compute the marginal effects. This can be the entire dataset, but also some meaningful subset, in which case we get the importance of the predictors within only that subset. Compute and display the Shapley values for the predictive model in `predictor.ra`, i.e. the fitted `ranger` random forest model, for the first-listed observation, and for the maximum-distance observation. ```{r plot.iml.shapley, fig.width=8, fig.height=6} set.seed(314) shapley <- iml::Shapley$new(predictor.ra, x.interest = X[1, ]) shapley$plot() shapley.maxdist <- Shapley$new(predictor.ra, x.interest = X[ix.maxdist, ]) shapley.maxdist$plot() ``` This figure shows the actual values of each predictor at the observation point, and the $\phi$ value, i.e., numerical contribution to the difference between actual and average prediction. These sum to the difference. List the results as a table: ```{r} (results <- shapley$results) sum(shapley$results$phi) (results.maxdist <- shapley.maxdist$results) sum(shapley.maxdist$results.maxdist$phi) ``` For the point near the river the actual is higher than expected; this is the sum of the individual $\phi$ values. The main positive contribution is from distance: it is closer than most points. For the distant point the actual is lower than the prediction; this is mainly because of distance: the point is further, and also lower and these have the most influence on the lower actual value. "Be careful to interpret the Shapley value correctly: The Shapley value is the *average contribution of a feature value to the prediction in different coalitions*. The Shapley value is NOT the difference in prediction when we would remove the feature from the model." # SHAP -- (SHapley Additive exPlanations) Another approach to Shapley values is SHapley Additive exPlanations (SHAP) by Lundberg and Lee (2017)^[Lundberg, Scott M., and Su-In Lee. “A unified approach to interpreting model predictions.” Advances in Neural Information Processing Systems (2017)]. Here the Shapley value explanation is represented as an *additive* feature attribution method, i.e., as a linear model. The original method for SHAP is `KernelSHAP` (see next section), but Lundberg et al.^[Lundberg, Scott M., Gabriel G. Erion, and Su-In Lee. “Consistent individualized feature attribution for tree ensembles.” arXiv preprint arXiv:1802.03888 (2018)] proposed `TreeSHAP`, a variant of SHAP for tree-based machine learning models including random forests. It is a fast, model-specific alternative to the more general `KernelSHAP`". One option to compute `TreeSHAP` is the `fastshap` package by Brandon Greenwll^[https://github.com/bgreenwell/fastshap]. This is based on the work of Štrumbelj &Kononenko.^[Štrumbelj, E., & Kononenko, I. (2014). Explaining prediction models and individual predictions with feature contributions. Knowledge and Information Systems, 41(3), 647-665. https://doi.org/10.1007/s10115-013-0679-x] This is an approximation using Monte-Carlo sampling: \[ \hat{\phi}_{j} = \frac{1}{M}\sum_{m=1}^M\left(\hat{f}(\mathbf{x}^{(m)}_{+j}) - \hat{f}(\mathbf{x}^{(m)}_{-j})\right) \] where $\hat{f}(\mathbf{x}^{(m)}_{+j})$ is the prediction for point $\mathbf{x}$, but with a random number of feature values replaced by feature values from a random data point $\mathbf{x}$, except for the respective value of feature $j$, i.e., the feature for which we want the SHAP value, which has its original value. From this estimate is subtracted the same estimate but with with value for feature $j$ also taken from the random data point. This isolates the effect of the feature within the original data point. ## `fastshap.explain` The key data structure for `fastshap` is an object of class `explain`. This is computed by the `fastshap::explain` function, which requires a prediction function. In this case it is `ranger::predict.ranger`, which is automatically called from the generic `predict` when the object is a `ranger` model. Note the use of `nsim`: "To obtain the most accurate results, `nsim` should be set as large as feasibly possible.". ```{r fastshap} library(fastshap, warn.conflicts = FALSE, verbose = FALSE) pfun <- function(object, newdata) { predict(object, data = newdata)$predictions } set.seed(314) system.time( fshap <- fastshap::explain(object = m.lzn.ra, X = X, pred_wrapper = pfun, adjust = TRUE, nsim = 128) ) class(fshap) ``` Each observation has a set of Shapley values. These are the contribution to the difference between the observed value and the average value of all observations. For the first (closest) point, there is a large positive contribution to the difference from `dist.m` and `elev`. How much do these differ from those computed by `iml`? Examine for the first observation ```{r fastshap.vs.iml} set.seed(314) shapley <- iml::Shapley$new(predictor.ra, x.interest = X[1, ]) print(data.frame(fastshap = fshap[1,], iml = shapley$results$phi)) ``` Similar but not identical. The SHAP method using `fastshap` finds more importance for most predictors. ### Parallelization This is a small dataset, but anyway let's see if the computation can be parallelized. ```{r parallel} library(doParallel) registerDoParallel(cores = 8) set.seed(314) system.time( fshap.p <- fastshap::explain(object = m.lzn.ra, X = X, pred_wrapper = pfun, nsim = 128, adjust = TRUE, parallel = TRUE) ) ``` Speed varies -- sometimes faster, sometimes not, when repeated. Are the values the same? ```{r} fshap[1:5,] fshap.p[1:5, ] fshap[1:5, ] - fshap.p[1:5, ] ``` No! Parallelization must be splitting the dataset and sending a piece to different cores. ## Visualing Shapley values The `shapviz` package can be used to visualize the Shapley values. Load the package and create a `shapviz` object from the `fastshap` explanation result. ```{r} library(shapviz, warn.conflicts = FALSE, verbose = FALSE) shv <- shapviz(fshap, X = X) class(shv) ``` Now we can use the several functions in `shapviz`. First, plot the mean Shapley values over all points, i.e., global variable importance but calculated from all individuals. ```{r shapviz.importance, fig.height=6, fig.width=8} sv_importance(shv, kind = "bar", show_numbers = TRUE) ``` We can see all individual importances with the "beeswarm" version of this plot: ```{r shapviz.importance.bs, fig.height=8, fig.width=10} sv_importance(shv, kind = "beeswarm") ``` This shows clearly that at close distances of `dist.m` the prediction is postively influenced by the predictor, and the reverse is true (but less strikingly) at far distances. We can see the effect of each predictor, compared to its value, with a "dependence" plot. Plot the dependence on elevation for all observations, using the `sv_dependence` function. We can get some idea of the interactions by colouring the points by the value of another predictor. Here we select the flood frequency. ```{r fastshap.auto.elev, fig.height=6, fig.width=10} sv_dependence(shv, v = "elev", color_var = "ffreq") sv_dependence(shv, v = "elev", color_var = "dist.m") ``` A consistent story: influence of elevation on the prediction is significant for both low and high points and the relation is almost linear. The value is positive for low elevations; the reverse is true for high elevations. Flood frequency is not well-related to Shapley values for all but the lowest points. And on distance to river for all observations, with both flood frequency and elevation shown as interactions: ```{r fastshap.auto.dist.m, fig.height=6, fig.width=10} sv_dependence(shv, v = "dist.m", color_var = "ffreq") sv_dependence(shv, v = "dist.m", color_var = "elev") ``` We see that for nearby points there is a strong positive influence of distance, i.e., the distance is very important to that prediction. For points $>250$ m there is also an influence but negative. The contribution of each predictor for the closest and furthest points can be visualized using the `sv_waterfall` function, or for more compact display, the `sv_force` function: ```{r sv.waterfall, fig.height=6, fig.width=10} sv_waterfall(shv, row_id = 1) sv_waterfall(shv, row_id = ix.maxdist) ``` ```{r sv.force, fig.height=2, fig.width=10} sv_force(shv, row_id = 1) sv_force(shv, row_id = ix.maxdist) ``` For this simple model these do not differ much from the Shapley values computed by `iml`. # Kernel SHAP Another approach to compute Shapley values is using the `kernelshap` package, which is an efficient implementation of the Kernal SHAP algorithm. See `?kernelshap` for a detailed explanation. This requires (1) a fitted model, (2) a set of features for which to display the values, with observation values at each point, (3) background data to integrate out "switched off" features. By default this is a random sample of 200 rows from the observations. These are the observtions that will be used to randomize and estimate the effects within a coalition. The `X` argument gives the points for which to compute. In a small data set this can be all the points. The `bg_X` argument is the "background" data, if null this requires a `bg_n` "size of background data" argument. Here again we have a small dataset so we could use all the points. More background points, more computation time, so restrict to about 2/3 of the dataset. ```{r kernelshap} library(kernelshap, warn.conflicts = FALSE, verbose = FALSE) # use existing ranger fits dim(X) system.time(s <- kernelshap(m.lzn.ra.l, X = X, bg_X = NULL, bg_n = 100, verbose = FALSE) ) str(s) ``` Compute the visualization object and show two diagnostic plots. Note that `shapviz` understands `kernelshap` objects, so there is no need to re-format the Shapley values matrix. First, the importances of each predictor at each point: ```{r shapviz.kernel.imp, fig.height=6, fig.width=10} sv <- shapviz(s) sv_importance(sv, kind = "bee") ``` Second, selected points: ```{r shapviz.kernel.force, fig.height=2, fig.width=10} sv_force(sv, row_id = 1) sv_force(sv, row_id = ix.maxdist) ``` Third, some dependence plots: ```{r shapviz.kernel.dep, fig.height=6, fig.width=10} # auto-select most important interacting feature sv_dependence(sv, v = "dist.m", color_var = "elev") sv_dependence(sv, v = "dist.m", color_var = "ffreq") ``` All of these are slightly different than the plots from the `fastshap` computation. This is because `kernelshap` uses a different approach to fast computation of the Shapley values.