Can a random forest be used for feature selection in multiple linear regression?

Question

Since RF can handle non-linearity but can't provide coefficients, would it be wise to use random forest to gather the most important features and then plug those features into a multiple linear regression model in order to obtain their coefficients?

Answer 1

Since RF can handle non-linearity but can't provide coefficients, would it be wise to use Random Forest to gather the most important Features and then plug those features into a Multiple Linear Regression model in order to explain their signs?

I interpret OP's one-sentence question to mean that OP wishes to understand the desirability of the following analysis pipeline:

1. Fit a random forest to some data
2. By some metric of variable importance from (1), select a subset of high-quality features.
3. Using the variables from (2), estimate a linear regression model. This will give OP access to the coefficients that OP notes RF cannot provide.
4. From the linear model in (3), qualitatively interpret the signs of the coefficient estimates.

I don't think this pipeline will accomplish what you'd like. Variables that are important in random forest don't necessarily have any sort of linearly additive relationship with the outcome. This remark shouldn't be surprising: it's what makes random forest so effective at discovering nonlinear relationships.

Here's an example. I created a classification problem with 10 noise features, two "signal" features, and a circular decision boundary.

set.seed(1)
N  <- 500
x1 <- rnorm(N, sd=1.5)
x2 <- rnorm(N, sd=1.5)

y  <- apply(cbind(x1, x2), 1, function(x) (x%*%x)<1)

plot(x1, x2, col=ifelse(y, "red", "blue"))
lines(cos(seq(0, 2*pi, len=1000)), sin(seq(0, 2*pi, len=1000))) 

And when we apply the RF model, we are not surprised to find that these features are easily picked out as important by the model. (NB: this model isn't tuned at all.)

x_junk   <- matrix(rnorm(N*10, sd=1.5), ncol=10)
x        <- cbind(x1, x2, x_junk)
names(x) <- paste("V", 1:ncol(x), sep="")

rf <- randomForest(as.factor(y)~., data=x, mtry=4)
importance(rf)

    MeanDecreaseGini
x1         49.762104
x2         54.980725
V3          5.715863
V4          5.010281
V5          4.193836
V6          7.147988
V7          5.897283
V8          5.338241
V9          5.338689
V10         5.198862
V11         4.731412
V12         5.221611

But when we down-select to just these two, useful features, the resulting linear model is awful.

summary(badmodel <- glm(y~., data=data.frame(x1,x2), family="binomial"))

The important part of the summary is the comparison of the residual deviance and the null deviance. We can see that the model does basically nothing to "move" the deviance. Moreover, the coefficient estimates are essentially zero.

Call:
glm(formula = as.factor(y) ~ ., family = "binomial", data = data.frame(x1, 
    x2))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.6914  -0.6710  -0.6600  -0.6481   1.8079  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.398378   0.112271 -12.455   <2e-16 ***
x1          -0.020090   0.076518  -0.263    0.793    
x2          -0.004902   0.071711  -0.068    0.946    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 497.62  on 499  degrees of freedom
Residual deviance: 497.54  on 497  degrees of freedom
AIC: 503.54

Number of Fisher Scoring iterations: 4

What accounts for the wild difference between the two models? Well, clearly the decision boundary we're trying to learn is not a linear function of the two "signal" features. Obviously if you knew the functional form of the decision boundary prior to estimating the regression, you could apply some transformation to encode the data in a way that regression could then discover... (But I've never known the form of the boundary ahead of time in any real-world problem.) Since we're only working with two signal features in this case, a synthetic data set without noise in the class labels, that boundary between classes is very obvious in our plot. But it's less obvious when working with real data in a realistic number of dimensions.

Moreover, in general, random forest can fit different models to different subsets of the data. In a more complicated example, it won't be obvious what's going on from a single plot at all, and building a linear model of similar predictive power will be even harder.

Because we're only concerned with two dimensions, we can make a prediction surface. As expected, the random model learns that the neighborhood around the origin is important.

M                 <- 100
x_new             <- seq(-4,4, len=M)
x_new_grid        <- expand.grid(x_new, x_new)
names(x_new_grid) <- c("x1", "x2")
x_pred            <- data.frame(x_new_grid, matrix(nrow(x_new_grid)*10, ncol=10))
names(x_pred)     <- names(x)

y_hat             <- predict(object=rf, newdata=x_pred, "vote")[,2]

library(fields)
y_hat_mat         <- as.matrix(unstack(data.frame(y_hat, x_new_grid), y_hat~x1))

image.plot(z=y_hat_mat, x=x_new, y=x_new, zlim=c(0,1), col=tim.colors(255), 
           main="RF Prediction surface", xlab="x1", ylab="x2")

As implied by our abysmal model output, the prediction surface for the reduced-variable logistic regression model is basically flat.

bad_y_hat     <- predict(object=badmodel, newdata=x_new_grid, type="response")
bad_y_hat_mat <- as.matrix(unstack(data.frame(bad_y_hat, x_new_grid), bad_y_hat~x1))
image.plot(z=bad_y_hat_mat, x=x_new, y=x_new, zlim=c(0,1), col=tim.colors(255), 
           main="Logistic regression prediction surface", xlab="x1", ylab="x2")

HongOoi notes that the class membership isn't a linear function of the features, but that it a linear function is under a transformation. Because the decision boundary is $1=x_1^2+x_2^2,$ if we square these features, we will be able to build a more useful linear model. This is deliberate. While the RF model can find signal in those two features without transformation, the analyst has to be more specific to get similarly helpful results in the GLM. Perhaps that's sufficient for OP: finding a useful set of transformations for 2 features is easier than 12. But my point is that even if a transformation will yield a useful linear model, RF feature importance won't suggest the transformation on its own.

Answer 2

The answer by @Sycorax is fantastic. In addition to those fully described aspects of the problem related to model fit, there is another reason not to pursue a multi-step process such as running random forests, lasso, or elastic net to "learn" which features to feed to traditional regression. Ordinary regression would not know about the penalization that properly went on during the development of the random forest or the other methods, and would fit unpenalized effects that are badly biased to appear too strong in predicting $Y$. This would be no different than running stepwise variable selection and reporting the final model without taking into account how it arrived.

Answer 3

A properly executed random forest applied to a problem that is more "random forest appropriate" can work as a filter to remove noise, and make results that are more useful as inputs to other analysis tools.

Disclaimers:

• Is it a "silver bullet"? No way. Mileage will vary. It works where it works, and not elsewhere.
• Are there ways you can badly wrongly grossly use it and get answers that are in the junk-to-voodoo domain? youbetcha. Like every analytic tool, it has limits.
• If you lick a frog, will your breath smell like frog? likely. I don't have experience there.

I have to give a "shout out" to my "peeps" who made "Spider". (link) Their example problem informed my approach. (link) I also love Theil-Sen estimators, and wish I could give props to Theil and Sen.

My answer isn't about how to get it wrong, but about how it might work if you got it mostly right. While I use "trivial" noise, I want you to think about "non-trivial" or "structured" noise.

One of the strengths of a random forest is how well it applies to high-dimensional problems. I can't show 20k columns (aka a 20k dimensional space) in a clean visual way. It is not an easy task. However, if you have a 20k-dimensional problem, a random forest might be a good tool there when most others fall flat on their "faces".

This is an example of removing noise from signal using a random forest.

#housekeeping
rm(list=ls())

#library
library(randomForest)

#for reproducibility
set.seed(08012015)

#basic
n <- 1:2000
r <- 0.05*n +1 
th <- n*(4*pi)/max(n)

#polar to cartesian
x1=r*cos(th) 
y1=r*sin(th)

#add noise
x2 <- x1+0.1*r*runif(min = -1,max = 1,n=length(n))
y2 <- y1+0.1*r*runif(min = -1,max = 1,n=length(n))

#append salt and pepper
x3 <- runif(min = min(x2),max = max(x2),n=length(n)/2)
y3 <- runif(min = min(y2),max = max(y2),n=length(n)/2)

x4 <- c(x2,x3)
y4 <- c(y2,y3)
z4 <- as.vector(matrix(1,nrow=length(x4)))

#plot class "A" derivation
plot(x1,y1,pch=18,type="l",col="Red", lwd=2)
points(x2,y2)
points(x3,y3,pch=18,col="Blue")
legend(x = 65,y=65,legend = c("true","sampled","false"),
col = c("Red","Black","Blue"),lty = c(1,-1,-1),pch=c(-1,1,18))

Let me describe what is going on here. This image below shows training data for class "1". Class "2" is uniform random over the same domain and range. You can see that the "information" of "1" is mostly a spiral, but has been corrupted with material from "2". Having 33% of your data corrupt can be a problem for many fitting tools. Theil-Sen starts to degrade at about 29%. (link)

Now we separate out the information, only having an idea of what noise is.

#Create "B" class of uniform noise
x5 <- runif(min = min(x4),max = max(x4),n=length(x4))
y5 <- runif(min = min(y4),max = max(y4),n=length(x4))
z5 <- 2*z4 

#assemble data into frame 
data <- data.frame(c(x4,x5),c(y4,y5),as.factor(c(z4,z5)))
names(data) <- c("x","y","z")

#train random forest - I like h2o, but this is textbook Breimann
fit.rf <- randomForest(z~.,data=data,
                       ntree = 1000, replace=TRUE, nodesize = 20)
data2 <- predict(fit.rf,newdata=data[data$z==1,c(1,2)],type="response")

#separate class "1" from training data
idx1a <- which(data[,3]==1)

#separate class "1" from the predicted data
idx1b <- which(data2==1)

#show the difference in classes before and after RF based filter
plot(data[idx1a,1],data[idx1a,2])
points(data[idx1b,1],data[idx1b,2],col="Red")

Here is the fitting result:

I really like this because it can show both strengths and weaknesses of a decent method to a hard problem at the same time. If you look near the center you can see how there is less filtering. The geometric scale of information is small and the random forest is missing that. It says something about number of nodes, number of trees, and sample density for class 2. There is also a "gap" near (-50,-50), and "jets" in several locations. In general, however, the filtering is decent.

Compare vs. SVM

Here is the code to allow a comparison with SVM:

#now to fit to svm
fit.svm <-  svm(z~., data=data, kernel="radial",gamma=10,type = "C")

x5 <- seq(from=min(x2),to=max(x2),by=1)
y5 <- seq(from=min(y2),to=max(y2),by=1)

count <- 1
x6 <- numeric()
y6 <- numeric()
for (i in 1:length(x5)){

     for (j in 1:length(y5)){
          x6[count]<-x5[i]
          y6[count]<-y5[j]
          count <- count+1
     }
}

data4 <- data.frame(x6,y6)
names(data4) <- c("x","y")

data4$z <- predict(fit.svm,newdata=data4)

idx4 <- which(data4$z==1,arr.ind=TRUE)


plot(data4[idx4,1],data4[idx4,2],col="Gray",pch=20)
points(data[idx1b,1],data[idx1b,2],col="Blue",pch=20)
lines(x1,y1,pch=18,col="Green", lwd=2)
grid()
legend(x = 65,y=65,
       legend = c("true","from RF","From SVM"),
       col = c("Green","Blue","Gray"),lty = c(1,-1,-1),pch=c(-1,20,15),pt.cex=c(1,1,2.25))

It results in the following image.

This is a decent SVM. The gray is the domain associated with class "1" by the SVM. The blue dots are the samples associated with class "1" by the RF. The RF based filter performs comparably to SVM without an explicitly imposed basis. It can be seen that the "tight data" near the center of the spiral is much more "tightly" resolved by the RF. There are also "islands" toward the "tail" where the RF finds association that the SVM does not.

I am entertained. Without having the background, I did one of the early things also done by a very good contributor in the field. The original author used "reference distribution" (link, link).

EDIT:

Apply random FOREST to this model:
While user777 has a nice thought about a CART being the element of a random forest, the premise of the random forest is "ensemble aggregation of weak learners". The CART is a known weak learner but it is nothing remotely near an "ensemble". The "ensemble" though in a random forest is intended "in the limit of a large number of samples". The answer of user777, in the scatterplot, uses at least 500 samples and that says something about human readability and sample sizes in this case. The human visual system (itself an ensemble of learners) is an amazing sensor and data processor and it finds that value to be sufficient for ease of processing.

If we take even the default settings on a random-forest tool, we can observe the behavior of the classification error increases for the first several trees, and does not reach the one-tree level until there are around 10 trees. Initially error grows reduction of error becomes stable around 60 trees. By stable I mean

x        <- cbind(x1, x2)
plot(rf,type="b",ylim=c(0,0.06))
grid()

Which yields:

If instead of looking at the "minimum weak learner" we look at the "minimum weak ensemble" suggested by a very brief heuristic for default setting of the tool the results are somewhat different.

Note, I used "lines" to draw the circle indicating the edge over the approximation. You can see that it is imperfect, but much better than the quality of a single learner.

The original sampling has 88 "interior" samples. If the sample sizes are increased (allowing ensemble to apply) then the quality of the approximation also improves. The same number of learners with 20,000 samples makes a stunningly better fit.

The much higher quality input information also allows evaluation of appropriate number of trees. Inspection of the convergence suggests that 20 trees is the minimum sufficient number in this particular case, to represent the data well.

Answer 4

Despite the legitimate warnings that this approach might fail in some cases, this should not discourage you from trying it out! Breimann reports an example (Breimann 2001) that selecting features by variable importance from a random forest and plugging them into logistic regresission outperformed variable selections specifically tailored for logistic regression, and others report similar observations, e.g., with using Boruta as a preprocessing variable selection step for logistic regression.

As both Random Forest variable importance computation and Boruta feature selection are readily available in R or other software and thus can be tested without much effort, this is something that should be given a try. If you have enough data, you can even validate the approach by doing both steps on different fractions of the data.