Obtaining knowledge from a random forest

Question

Random forests are considered to be black boxes, but recently I was thinking what knowledge can be obtained from a random forest?

The most obvious thing is the importance of the variables, in the simplest variant it can be done just by calculating the number of occurrences of a variable.
The second thing I was thinking about are interactions. I think that if the number of trees is sufficiently large then the number of occurrences of pairs of variables can be tested (something like chi square independence). The third thing are nonlinearities of variables. My first idea was just to look at a chart of a variable Vs score, but I'm not sure yet whether it makes any sense.

Added 23.01.2012
Motivation

I want to use this knowledge to improve a logit model. I think (or at least I hope) that it is possible to find interactions and nonlinearities that were overlooked.

Answer 1

Random Forests are hardly a black box. They are based on decision trees, which are very easy to interpret:

#Setup a binary classification problem
require(randomForest)
data(iris)
set.seed(1)
dat <- iris
dat$Species <- factor(ifelse(dat$Species=='virginica','virginica','other'))
trainrows <- runif(nrow(dat)) > 0.3
train <- dat[trainrows,]
test <- dat[!trainrows,]

#Build a decision tree
require(rpart)
model.rpart <- rpart(Species~., train)

This results in a simple decision tree:

> model.rpart
n= 111 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 111 35 other (0.68468468 0.31531532)  
  2) Petal.Length< 4.95 77  3 other (0.96103896 0.03896104) *
  3) Petal.Length>=4.95 34  2 virginica (0.05882353 0.94117647) *

If Petal.Length < 4.95, this tree classifies the observation as "other." If it's greater than 4.95, it classifies the observation as "virginica." A random forest is simple a collection of many such trees, where each one is trained on a random subset of the data. Each tree then "votes" on the final classification of each observation.

model.rf <- randomForest(Species~., train, ntree=25, proximity=TRUE, importance=TRUE, nodesize=5)
> getTree(model.rf, k=1, labelVar=TRUE)
  left daughter right daughter    split var split point status prediction
1             2              3  Petal.Width        1.70      1       <NA>
2             4              5 Petal.Length        4.95      1       <NA>
3             6              7 Petal.Length        4.95      1       <NA>
4             0              0         <NA>        0.00     -1      other
5             0              0         <NA>        0.00     -1  virginica
6             0              0         <NA>        0.00     -1      other
7             0              0         <NA>        0.00     -1  virginica

You can even pull out individual trees from the rf, and look at their structure. The format is slightly different than for rpart models, but you could inspect each tree if you wanted and see how it's modeling the data.

Furthermore, no model is truly a black box, because you can examine predicted responses vs actual responses for each variable in the dataset. This is a good idea regardless of what sort of model you are building:

library(ggplot2)
pSpecies <- predict(model.rf,test,'vote')[,2]
plotData <- lapply(names(test[,1:4]), function(x){
  out <- data.frame(
    var = x,
    type = c(rep('Actual',nrow(test)),rep('Predicted',nrow(test))),
    value = c(test[,x],test[,x]),
    species = c(as.numeric(test$Species)-1,pSpecies)
    )
  out$value <- out$value-min(out$value) #Normalize to [0,1]
  out$value <- out$value/max(out$value)
  out
})
plotData <- do.call(rbind,plotData)
qplot(value, species, data=plotData, facets = type ~ var, geom='smooth', span = 0.5)

I've normalized the variables (sepal and petal length and width) to a 0-1 range. The response is also 0-1, where 0 is other and 1 is virginica. As you can see the random forest is a good model, even on the test set.

Additionally, a random forest will compute various measure of variable importance, which can be very informative:

> importance(model.rf, type=1)
             MeanDecreaseAccuracy
Sepal.Length           0.28567162
Sepal.Width           -0.08584199
Petal.Length           0.64705819
Petal.Width            0.58176828

This table represents how much removing each variable reduces the accuracy of the model. Finally, there are many other plots you can make from a random forest model, to view what's going on in the black box:

plot(model.rf)
plot(margin(model.rf)) 
MDSplot(model.rf, iris$Species, k=5)
plot(outlier(model.rf), type="h", col=c("red", "green", "blue")[as.numeric(dat$Species)])

You can view the help files for each of these functions to get a better idea of what they display.

Answer 2

Some time ago I had to justify a RF model-fit to some chemists in my company. I spent quite time trying different visualization techniques. During the process, I accidentally also came up with some new techniques which I put into an R package (forestFloor) specifically for random forest visualizations.

The classical approach are partial dependence plots supported by: Rminer(data-based sensitivity analysis is reinvented partial dependence), or partialPlot in randomForest package. I find the partial dependence package iceBOX as an elegant way to discover interactions. Have not used edarf package, but seems to have some fine visualizations dedicated for RF. The ggRandomForest package also contain a large set of useful visualizations.

Currently forestFloor supports randomForest objects(support for other RF implementions is on its way). Also feature contributions can be computed for gradient boosted trees, as these trees after training are not much different from random forest trees. So forestFloor could support XGBoost in future. Partial dependence plots are completely model invariant.

All packages have in common to visualize the geometrical mapping structure of a model from feature space to target space. A sine curve y = sin(x) would be a mapping from x to y and can be plotted in 2D. To plot a RF mapping directly would often require too many dimensions. Instead the overall mapping structure can be projected, sliced or decomposed, such that the entire mapping structure is boiled down into a sequence of 2D marginal plots. If your RF model only has captured main effects and no interactions between variables, classic visualizations methods will do just fine. Then you can simplify your model structure like this $y = F(X) \approx f_1(x_1) + f_2(x_2) + ... + f_d(x_d)$. Then each partial function by each variable can be visualized just as the sine curve. If your RF model has captured sizable interactions, then it is more problematic. 3D slices of the structure can visualize interactions between two features and the output. The problem is to know which combination of features to visualize, (iceBOX does address this issue). Also it is not easy to tell if other latent interactions still are not accounted for.

In this paper, I used an very early version of forestFloor to explain what actual biochemical relationship a very small RF model had captured. And in this paper we thoroughly describe visualizations of feature contributions, Forest Floor Visualizations of Random Forests.

I have pasted the simulated example from forestFloor package, where I show how to uncover a simulated hidden function $y = {x_1}^2 + sin(x_2\pi) + 2 * x_3 * x_4 + $ noise

#1 - Regression example:
set.seed(1234)
library(forestFloor)
library(randomForest)

#simulate data y = x1^2+sin(x2*pi)+x3*x4 + noise
obs = 5000 #how many observations/samples
vars = 6   #how many variables/features
#create 6 normal distr. uncorr. variables
X = data.frame(replicate(vars,rnorm(obs)))
#create target by hidden function
Y = with(X, X1^2 + sin(X2*pi) + 2 * X3 * X4 + 0.5 * rnorm(obs)) 

#grow a forest
rfo = randomForest(
  X, #features, data.frame or matrix. Recommended to name columns.
  Y, #targets, vector of integers or floats
  keep.inbag = TRUE,  # mandatory,
  importance = TRUE,  # recommended, else ordering by giniImpurity (unstable)
  sampsize = 1500 ,   # optional, reduce tree sizes to compute faster
  ntree = if(interactive()) 500 else 50 #speedup CRAN testing
)

#compute forestFloor object, often only 5-10% time of growing forest
ff = forestFloor(
  rf.fit = rfo,       # mandatory
  X = X,              # mandatory
  calc_np = FALSE,    # TRUE or FALSE both works, makes no difference
  binary_reg = FALSE  # takes no effect here when rfo$type="regression"
)


#plot partial functions of most important variables first
plot(ff,                       # forestFloor object
     plot_seq = 1:6,           # optional sequence of features to plot
     orderByImportance=TRUE    # if TRUE index sequence by importance, else by X column  
)

#Non interacting features are well displayed, whereas X3 and X4 are not
#by applying color gradient, interactions reveal themself 
#also a k-nearest neighbor fit is applied to evaluate goodness-of-fit
Col=fcol(ff,3,orderByImportance=FALSE) #create color gradient see help(fcol)
plot(ff,col=Col,plot_GOF=TRUE) 

#feature contributions of X3 and X4 are well explained in the context of X3 and X4
# as GOF R^2>.8


show3d(ff,3:4,col=Col,plot_GOF=TRUE,orderByImportance=FALSE)

Lastly the code for partial dependence plots coded by A.Liaw described by J.Friedman. Which do fine for main effects.

par(mfrow=c(2,3))
for(i in 1:6) partialPlot(rfo,X,x.var=names(X)[i])

Answer 3

To supplement these fine responses, I would mention use of gradient boosted trees (e.g. the GBM Package in R). In R, I prefer this to random forests because missing values are allowed as compared to randomForest where imputation is required. Variable importance and partial plots are available (as in randomForest) to aid in feature selection and nonlinear transformation exploration in your logit model. Further, variable interaction is addressed with Friedman’s H-statistic (interact.gbm) with reference given as J.H. Friedman and B.E. Popescu (2005). “Predictive Learning via Rule Ensembles.” Section 8.1. A commercial version called TreeNet is available from Salford Systems and this video presentation speaks to their take on variable interaction estimation Video.

Answer 4

Late answer, but I came across a recent R package forestFloor (2015) that helps you doing this "unblackboxing" task in an automated fashion. It looks very promising!

library(forestFloor)
library(randomForest)
#simulate data
obs=1000
vars = 18
X = data.frame(replicate(vars,rnorm(obs)))
Y = with(X, X1^2 + sin(X2*pi) + 2 * X3 * X4 + 1 * rnorm(obs))
#grow a forest, remeber to include inbag
rfo=randomForest(X,Y,keep.inbag = TRUE,sampsize=250,ntree=50)
#compute topology
ff = forestFloor(rfo,X)
#ggPlotForestFloor(ff,1:9)
plot(ff,1:9,col=fcol(ff))

Produces the following plots:

It also provides three-dimensional visualization if you are looking for interactions.

Answer 5

As mentioned by Zach, one way of understanding a model is to plot the response as the predictors vary. You can do this easily for "any" model with the plotmo R package. For example

library(randomForest)
data <- iris
data$Species <- factor(ifelse(data$Species=='virginica','virginica','other'))
mod <- randomForest(Species~Sepal.Length+Sepal.Width, data=data)
library(plotmo)
plotmo(mod, type="prob")

which gives

This changes one variable while holding the others at their median values. For interaction plots, it changes two variables. (Note added Nov 2016: plotmo now also supports partial dependence plots.)

The example above uses only two variables; more complicated models can be visualized in a piecemeal fashion by looking at one or two variables at a time. Since the "other" variables are held at their median values, this shows only a slice of the data, but can still be useful. Some examples are in the vignette for the plotmo package. Other examples are in Chapter 10 of Plotting rpart trees with the rpart.plot package.

Answer 6

Late in the game but there are some new developments in this front, for example LIME and SHAP. Also a package worth checking is DALEX (in particular if using R but in any case contains nice cheatsheets etc.), though doesn't seem to cover interactions at the moment. And these are all model-agnostic so will work for random forests, GBMs, neural networks, etc.

Answer 7

I'm very interested in these types of questions myself. I do think there is a lot of information we can get out of a random forest.

About interactions, it seems like Breiman and Cutler have already tried to look at it, especially for classification RFs.

To my knowledge, this has not been implemented in the randomForest R package. Maybe because it might not be as simple and because the meaning of "variable interactions" is very dependent of your problem.

About the nonlinearity, I'm not sure what you are looking for, regression forest are used for nonlinear multiple regression problems without any priors on what type of nonlinear function to use.

Answer 8

A slight modification of random forests that provide more information about the data are the recently-developed causal forest methods. See the GRF R-package and the motivating paper here. The idea is to use the random forest baseline methods to find heterogeneity in causal effects.

An earlier paper (here) gives a detailed approach to a simple causal forest. Page 9 of the paper gives a step-by-step procedure for growing a causal tree, which can then be expanded to a forest in the usual ways.

Equation 4:

Equation 5:

Answer 9

Late answer related to my question here (Can we make Random Forest 100% interpretable by fixing the seed?):

Let $z_1$ be the seed in the creation of boostrapped training set, and $z_2 $ be the seed in the selection of feature's subset (for simplification, I only list 2 kinds of seeds here).

1. From $z_1$, $m$ boostrapped training sets are created: $D_1(z_1)$, $D_2(z_1)$, $D_3(z_1)$, ..., $D_m(z_1)$.
2. From those traning sets, $m$ corresponding decision trees are created, and tuned via cross-validation: $T_1(z_1,z_2)$, $T_2(z_1,z_2)$, $T_3(z_1,z_2)$,..., $T_m(z_1,z_2)$.
3. Let's denote predictions from the ${j^\text{th}}_{(j=1,2,...,m)}$ tree for an individual $x_i$ (from training or testing set, whatever) as $\hat{f}^j(x_i)_{(i \le n, j \le m)}$. Hence the final predictions by the ensemble trees are: $$\hat{F}(x_i) = > \frac{1}{m}\sum\limits_{j=1}^m \hat{f}^j(x_i)$$
4. Once the model is validated, and is stable (meaning $\hat{F}(x_i)$ doesn't depend strongly on the pair $(z_1,z_2)$). I start to create every possible combinations of my features, which give me a very big set ($x'_i$).
5. Applying my forest on each $x'_i$ gives me the corresponding predictions: $$x'_1 \rightarrow \hat{F}(x'_1) \text{ - which is fixed > thanks to $(z_1, z_2)$}$$ $$x'_2 \rightarrow \hat{F}(x'_2) \text{ - > which is fixed thanks to $(z_1, z_2)$}$$ $$x'_3 \rightarrow > \hat{F}(x'_3) \text{ - which is fixed thanks to $(z_1, z_2)$}$$ $$x'_4 > \rightarrow \hat{F}(x'_4) \text{ - which is fixed thanks to $(z_1, > z_2)$}$$ $$....$$
6. The latter can be easily represented in form of a single (huge) tree. For example: $x'_1$: (Age = 18, sex = M, ...), $x'_2$ = (Age = 18, sex = F, ...), ... could be regrouped to create a leaf.

This works also for every ensemble methods based on aggregation of trees.