Linear regression performing better than random forest in Caret

Question

I am dealing with fairly highly dimensional data (spectroscopic data), which I preprocess with PCA, then apply the chosen learning algorithm. When I use linear regression, using the Caret package, I get an $R^2=0.83$ on my test set of 300 PCA components however, when I use random forests on the same set:

modelFit <- train(train$age~.,data=trainPC[,1:300],method = "rf")

I get worse performance, with $R^2$ much lower on the test set. RMSE is also not as good for random forests vs linear regression.

For clarity - I am using Caret's in-built cross validation, and am evaluating the performance on an independent test set.

I get similar results when I use xgboost, neural networks etc. I know that different machine learning algorithms can perform better in different circumstances, but I am nonetheless surprised that linear regression appears to be performing best, and makes me suspect that I am approaching things incorrectly, or making a basic error with using the Caret package.

I have also tried using PLS (the main feature selection method used for my particular example) and get an $R^2=0.83$. The literature suggests that neural nets/non-linear-methods should outperform the linear methods I am using, but at the moment, I cannot reproduce this.

Does anyone have any suggests/tips on this? I also attach the data set I am using (I am trying to predict age) - data link

I also attach my current R script that I am using to do the learning - any help would be much appreciated! - R script

Answer 1

Definitely check on how $R^2$ is being evaluated and whether it is uniform across the different algorithms. Beyond that, a few thoughts occur to me:

• If your features have a smooth, nearly linear dependence on the covariates, then linear regression will model the dependence better than random forests, which will basically approximate a linear curve with an ugly irregular step function. If the dependence is multivariate linear and smooth, with $v$ significant covariates producing the dependence, the fit performance of random forests can be expected to get worse and worse for larger and larger $v$. RF has a much greater ability than a single decision tree to model linearity, since we are adding tree predictions together - but still, it's just not very efficient to approximate a high-dimensional linear relationship with a series of step functions. I think this is the most likely theoretical explanation for RF underperforming linear regression.

• ...that said, the step functions will get nicer as you add more trees to the random forest. So you may want to consider increasing the number of trees a lot - as high as you are willing to wait for. Maybe see if caret will allow you to track the $R^2$ as the number of trees increases? I would also expect boosted trees to potentially work well here, but you may have to fit a ton of trees and fiddle a lot with the parameters. There's no guarantee though - if this is a really high-dimensional smooth relationship, it could take an unfathomably large number of trees to approximate it well.

• Preprocessing with PCA is not ideal for high-dimensional learners like random forests. For best results, they should generally get direct access to the full dataset. If the raw dataset has sparsity of main effects, this nice property could easily be destroyed by PCA. This would make even worse the problem of my first bullet, where RF is trying to approximate a high-dimensional linear relationship by step functions.

• Advanced machine learning algorithms are favored for their generalization performance, not their ability to closely fit training data and produce a high $R^2$. You are performing an interesting diagnostic and that is worth doing, but be aware that this diagnostic is not a good measure of final performance.

• Random Forests are not the only useful algorithm out there for high-dimensional data. You could also try elastic net regression either on your PCA processed data or on the raw dataset. Elastic net gives you the ability to model smooth relationships, as well as the ability to perform high-dimensional feature selection using the $\ell_1$ regression coefficient penalty.

Answer 2

Check out caret's function findLinearCombos and run it on your data. The object returned is a list--the second element is a vector of indices which can safely be removed from your data set since they are linear combinations of other columns (which, now that I think about it makes total sense because you only had 500 observations and, as a result, your model was totally over-specified).

When you pass in your data frame the resulting vector has a thousand elements. Basically, you can go right ahead and delete columns 501 to 1500.

Do that from the get-go and try re-building those models.

---

A couple quick side notes:

• If you include the 'data' argument, then only refer to the columns by their names in the formula (i.e. write age ~. instead of data$age ~.).
• If you plan on doing PCA then you might want to consider setting it up via the preProcess argument in the train function (you can tweak the threshold in trainControl function). (Since you're already using carat anyway you might as well have it handle your pre-processing, too--it gives you one less thing to worry about).

---

By the way, try checking out this code.

prepped <- preProcess(x=df[, -1], method=c('center', 'scale', 'pca'), 
                      thresh=0.9999)
pca <- predict(prepped, df[, -1])
pca$age <- df$age
fit <- ln(age ~ ., data=pca)

num.unique <- length(unique(df$age))

cols <- rainbow(num.unique)[factor(df$age)]

plot(fit, col=cols)