--- title: "Getting Started with randomPlantedForest" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting Started with randomPlantedForest} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5 ) # The classification and recipes sections use suggested packages has_suggests <- requireNamespace("palmerpenguins", quietly = TRUE) && requireNamespace("recipes", quietly = TRUE) has_plot <- requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("patchwork", quietly = TRUE) ``` `randomPlantedForest` implements Random Planted Forest ([Hiabu, Mammen & Meyer](https://arxiv.org/abs/2012.14563)), a directly interpretable tree ensemble. Its defining feature is that the maximum order of interaction between predictors can be bounded via `max_interaction`: the fitted model then decomposes exactly into main effects, two-way interactions, and so on, a functional ANOVA decomposition, with no post-hoc approximation needed. ```{r setup} library(randomPlantedForest) ``` ## Fitting a regression forest To see what the model recovers, we simulate data where the true components are known. The outcome depends on `x1`, `x2` and `x3` through main effects of different shapes and on `x1` and `x2` jointly through a product interaction with `x4` having no effect at all: $$y = 1 + \sin(\pi x_1) + 2\left(x_2^2 - \tfrac{1}{3}\right) + \tfrac{x_3}{2} + 2 x_1 x_2 + \varepsilon$$ With uniform predictors on $[-1, 1]$, each term has mean zero and the interaction has no main-effect part, so the terms are exactly the components the model is meant to estimate. ```{r simulate} set.seed(1) n <- 1000 sim <- data.frame( x1 = runif(n, -1, 1), x2 = runif(n, -1, 1), x3 = runif(n, -1, 1), x4 = runif(n, -1, 1) ) true_effects <- list( x1 = function(x) sin(pi * x), x2 = function(x) 2 * (x^2 - 1 / 3), x3 = function(x) x / 2 ) sim$y <- 1 + true_effects$x1(sim$x1) + true_effects$x2(sim$x2) + true_effects$x3(sim$x3) + 2 * sim$x1 * sim$x2 + rnorm(n, sd = 0.3) ``` ```{r fit-regression} rpfit <- rpf(y ~ ., data = sim, max_interaction = 2, ntrees = 50) rpfit ``` `max_interaction = 2` restricts the forest to main effects and two-way interactions. This is the central tuning choice of the method: lower values give a more interpretable, more constrained model, while `max_interaction = 0` allows interactions of any order. `ntrees` sets the number of tree families in the ensemble. Prediction works as usual: ```{r predict-regression} predict(rpfit, new_data = head(sim)) ``` Note that `rpf()` does not handle missing values: rows containing `NA` must be imputed or removed beforehand (see the `recipes` section below for one way to do this as part of the model specification). ## Inside a tree family The forest averages `ntrees` tree families. Unlike a single CART tree, a family consists of several trees, each splitting on its own set of variables, and it predicts a point by adding up the values of all leaves containing it. New trees grow from existing ones by splitting on an additional variable, and `max_interaction = 2` limits each tree to at most two variables. To look at one family, we fit a forest with a single, small family to the same data: ```{r fit-small, eval = has_plot, fig.width = 10, fig.height = 6, out.width = "100%"} small <- rpf(y ~ ., data = sim, max_interaction = 2, ntrees = 1, splits = 10) plot(small) ``` The trees are arranged in rows by the number of variables they split on, with the leaves coloured by their value. `plot(type = "boxes")` draws the same leaves in variable space instead, one panel per tree: ```{r boxes-small, eval = has_plot, fig.width = 10, fig.height = 6, out.width = "100%"} plot(small, type = "boxes") ``` These trees are not yet the components of the decomposition. The leaves of different trees overlap, so an effect can be spread over several trees: the quadratic effect of `x2`, for instance, can just as well be carried by trees on `x2` together with another variable as by a tree on `x2` alone. Purification moves each effect to the lowest-order tree it belongs to, without changing the family's predictions: ```{r boxes-purified, eval = has_plot, fig.width = 10, fig.height = 6, out.width = "100%"} purify(small) plot(small, type = "boxes") ``` Now each panel is one component: the main effects resemble the true shapes, `x1:x2` shows the sign pattern of a product, and the root value is the intercept, close to the true value of 1. Averaging these components over all tree families gives the decomposition of the whole forest, which `predict_components()` computes. With a single family and few splits the estimates are coarse, which is why the next section uses the forest of 50 families fitted above. ## Decomposing predictions Because the interaction order is bounded, the prediction function is a sum of low-dimensional components: $$\hat{m}(x) = \hat{m}_0 + \sum_k \hat{m}_k(x_k) + \sum_{k < l} \hat{m}_{kl}(x_k, x_l)$$ `predict_components()` returns exactly this decomposition: an intercept and one column per main effect and interaction term, evaluated on the supplied data. The forest is *purified* internally to make the decomposition unique, which is why the print method above reports the purification state. ```{r components} components <- predict_components(rpfit, new_data = sim) components$intercept names(components$m) ``` The components sum exactly to the model prediction: ```{r components-sum} pred_sum <- components$intercept + rowSums(components$m) pred_direct <- predict(rpfit, new_data = sim)$.pred all.equal(pred_sum, pred_direct) ``` Each column of `components$m` is the isolated contribution of one predictor (or one pair) to the prediction, on the scale of the outcome. Their standard deviations give a simple importance measure: the three main effects and the `x1:x2` interaction dominate, while `x4` and the other interactions are close to zero, as they should be. ```{r components-sd} round(sort(sapply(components$m, sd), decreasing = TRUE), 3) ``` The estimated main effects follow the true functions (blue lines). As with any tree-based model, they flatten out towards the edges of the data, where few observations are available to determine them. ```{r main-effects, fig.height = 3} op <- par(mfrow = c(1, 3), mar = c(4, 4, 1, 1)) for (v in names(true_effects)) { ylim <- range(components$m[[v]], true_effects[[v]](sim[[v]])) plot( sim[[v]], components$m[[v]], pch = 20, col = "darkgray", ylim = ylim, xlab = v, ylab = "component" ) curve(true_effects[[v]](x), add = TRUE, lwd = 2, col = "blue") } par(op) ``` The interaction is a function of two variables, so it is shown as a heatmap. We evaluate the estimated `x1:x2` component on a grid. Note `x3` and `x4` must be present in the new data but do not affect this component. ```{r interaction, fig.height = 3.5} grid_points <- seq(-1, 1, length.out = 50) grid <- expand.grid(x1 = grid_points, x2 = grid_points, x3 = 0, x4 = 0) estimated <- predict_components(rpfit, new_data = grid, predictors = c("x1", "x2"))$m[["x1:x2"]] surfaces <- list( true = outer(grid_points, grid_points, function(x1, x2) 2 * x1 * x2), estimated = matrix(estimated, nrow = length(grid_points)) ) zlim <- c(-1, 1) * max(abs(unlist(surfaces))) op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1)) for (s in names(surfaces)) { image( grid_points, grid_points, surfaces[[s]], zlim = zlim, col = hcl.colors(51, "Blue-Red 3"), xlab = "x1", ylab = "x2", main = paste(s, "x1:x2") ) } par(op) ``` The plots above are drawn manually from `predict_components()` for illustration. The [`glex`](https://plantedml.com/glex/) package provides `ggplot2` versions of them: `glex()` computes the decomposition for a fitted `rpf` object, and `plot_main_effect()` and `plot_twoway_effects()` draw main effects and two-way interactions. ## Classification For a factor outcome, `rpf()` automatically fits a classification forest. We illustrate with the `penguins` data from `palmerpenguins`, removing missing values first. ### Binary ```{r classif-binary, eval = has_suggests} penguins <- na.omit(palmerpenguins::penguins) rpfit_sex <- rpf(sex ~ ., data = penguins, ntrees = 50) predict(rpfit_sex, head(penguins), type = "prob") predict(rpfit_sex, head(penguins), type = "class") ``` With the default `loss = "exponential"` (or `"logit"`), `type = "link"` additionally gives raw predictions on the log-odds scale. ### Multiclass Multiclass classification works the same way: probability predictions contain one column per class: ```{r classif-multiclass, eval = has_suggests} rpfit_species <- rpf(species ~ ., data = penguins, ntrees = 50) predict(rpfit_species, head(penguins), type = "prob") predict(rpfit_species, head(penguins), type = "class") ``` The default loss for classification is `"exponential"`, which gives proper probability estimates. `"logit"` gives similar results but is slower, while `"L1"` and `"L2"` fit class indicators directly and do not yield proper probabilities. Note that with `loss = "logit"`, multiclass models are fit on a reference-class (log-odds) scale, so link predictions and components have one column less than there are classes: see `?predict.rpf` for details. ## Using `recipes` `rpf()` also accepts a [`recipes`](https://recipes.tidymodels.org/) recipe in place of a formula, so preprocessing can be bundled with the model specification. This makes it easy to add preprocessing steps, such as imputing missing values, which `rpf()` itself does not accept: ```{r recipes, message = FALSE, eval = has_suggests} library(recipes) # rows with a missing outcome still need to be dropped by hand peng_raw <- subset(palmerpenguins::penguins, !is.na(body_mass_g)) rec <- recipe(body_mass_g ~ ., data = peng_raw) |> step_impute_median(all_numeric_predictors()) |> step_impute_mode(all_nominal_predictors()) rpfit_rec <- rpf(rec, data = peng_raw, ntrees = 50) predict(rpfit_rec, head(peng_raw)) ``` The preprocessing steps are stored with the model and applied automatically in `predict()`. ## Saving and restoring models An `rpf` object holds its forest in an external pointer to a C++ object, which `saveRDS()` cannot preserve, meaning that a naively saved and restored model is not usable: ```{r saveRDS-trap, error = TRUE} tmp <- tempfile(fileext = ".rds") saveRDS(rpfit, tmp) broken <- readRDS(tmp) rpf_is_valid(broken) predict(broken, head(sim)) ``` Use `rpf_marshal()` to convert the model into a plain R list before saving, and `rpf_unmarshal()` after loading: ```{r marshal} saveRDS(rpf_marshal(rpfit), tmp) restored <- rpf_unmarshal(readRDS(tmp)) all.equal( predict(restored, head(sim)), predict(rpfit, head(sim)) ) ``` Alternatively, a [`bundle`](https://rstudio.github.io/bundle/) method wrapping the same mechanism is available via `bundle::bundle()` and `bundle::unbundle()`. ```{r cleanup, include = FALSE} unlink(tmp) ```