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

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?

• @user777 you mean, you're reading "explain" as "obtain" or something like that? Jul 31, 2015 at 3:05
• Since there appears to be a substantial concern about what this question might be trying to ask, could you please edit it to clear up the confusion?
– whuber
Jul 31, 2015 at 15:19

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")
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.

• I think that, while nice, it misses the whole point. Try a slightly more complex model - one I would call sufficiently non-trivial enough to be informative. Consider this (people.kyb.tuebingen.mpg.de/spider/demo_dec_vs_svm_2.jpg) with some salt-and-pepper noise. Now fit it with RF. Now use RF as "filter" to remove "spurious" phenomena, and train your alternate model on the outputs of the RF, not its coefficients. The OP specifies "features" not coefficients. The OP does not say "reduced random forest" but "random forest". The implication is for the output of the complete forest. Jul 31, 2015 at 12:54
• Interestingly, your example is linear in the predictors -- if you transform them first. Your equation is $y = x_1^2 + x_2^2$, so taking the square of $x_1$ and $x_2$ turns it into a linear model. In fact, it's one of the examples that Friedman, Hastie & Tibshirani used to show that boosting fits a low-dimensional model (additive, in the case of stumps). Jul 31, 2015 at 15:18
• @user777, I have used random forests for dimensionality reduction against complex problems for years. My former employer, Intel Semiconductor, has production processes that have 20k columns and 20k steps, and for under 10 defects uses this as part of a suite of analytic tools to go from 20k columns, to 30 columns. They teach courses internally, and lecture on the content externally. (web.stanford.edu/class/ee392m/Lecture3Tuv.pdf) The "forest" is an amazingly powerful tool and you took the one bad road. There are two solidly good ones you missed. Aug 1, 2015 at 14:44
• @EngrStudent "Intel uses RF to go from many to fewer columns." That's exactly what my demonstration does. It also answers OP's particular question about whether those features may be profitably used in a linear model, and my answer is "in general they will not be directly useful, but in this example, you may find them useful under a transformation."
– Sycorax
Aug 1, 2015 at 16:36

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.

• Thanks, Dr. Harrell! In my mind, if OP were to use this analysis pipeline, OP would do the RF down-selection and linear model fitting inside of their (cross-)validation scheme. Is that sufficient to mitigate the bias you describe, or is there another problem lurking here that I'm missing?
– Sycorax
Aug 1, 2015 at 16:08
• I don't know how to do a proper linear model "inside" the scheme. I don't know how to mitigate the overfitting/estimation bias that would result. Perhaps a better approach is model approximation, sometimes called pre-conditioning. Here you use traditional models or single trees (which would require a ton of nodes) to approximate the output of a black box, inheriting the shrinkage of the black box. Aug 1, 2015 at 16:51
• @FrankHarrell Suppose one cross-validates LASSO solutions to arrive at the set of penalized predictors that across many trials best reduces prediction error. This seems to me very different from your example of stepwise variable selection, which is so much more opportunistic and would seem to produce much greater bias. It would also mitigate bias to select predictors with only substantial coefficients in LASSO. Of course choosing those penalized to the verge of zero and using them in regression with no limit on coefficient magnitude there would be counter-productive. Apr 5, 2020 at 17:51
• I don't see how that does anything but cover up the problem with model uncertainty. Lasso simply cannot find the right variables with a high probability. Apr 17, 2020 at 15:42

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)

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.

• How does this prove/disprove that RF can be used to select high-quality features for a linear model? You don't discuss feature selection or linear models in your answer.
– Sycorax
Aug 1, 2015 at 16:17
• the premise of the asker is about using RF as a pre-filter then looking at parameters of the linear model. He was not looking at parameters of the forest. When I read it, because of its lack of depth or specificity, the "multiple linear model" seemed secondary. To me this looked like a "can a RF preprocess/clean data before sticking it into analytic tool x". Aug 1, 2015 at 20:29
• I see horizontal and vertical artifacts. I am tempted to make a "conjugate data" on coordinates rotated 45 degrees, and append the columns to the RF. I bet that a transform like that would reduce the number of artifacts, though I'm not sure how it would change input parameters. Nov 11, 2015 at 14:19
• As an aside, the "rotation forest" is a recent-ish variation on this idea, where PCA is used to determine the new basis. ncbi.nlm.nih.gov/pubmed/16986543
– Sycorax
Nov 16, 2015 at 18:30
• Neat: Liu, Fei Tony, Ting, Kai Ming and Zhou, Zhi-Hua. “Isolation forest.” Data Mining, 2008. ICDM‘08. Eighth IEEE International Conference on. Dec 11, 2017 at 19:16

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.