# Compare R-squared from two different Random Forest models

I'm using the randomForest package in R to develop a random forest model to try to explain a continuous outcome in a "wide" dataset with more predictors than samples.

Specifically, I'm fitting one RF model allowing the procedure to select from a set of ~75 predictor variables that I think are important.

I'm testing how well that model predicts the actual outcome for a reserved testing set, using the approach posted here previously, namely, ... or in R:

1 - sum((y-predicted)^2)/sum((y-mean(y))^2)


But now I have an additional ~25 predictor variables that I can add. When using the set of ~100 predictors, the R² is higher. I want to test this statistically, in other words, when using the set of ~100 predictors, does the model test significantly better in testing data than the model fit using ~75 predictors. I.e., is the R² from testing the RF model fit on the full dataset significantly higher than the R² from testing the RF model on the reduced dataset.

This is important for me to test, because this is pilot data, and getting those extra 25 predictors was expensive, and I need to know whether I should pay to measure those predictors in a larger follow-up study.

I'm trying to think of some kind of resampling/permutation approach but nothing comes to mind.

Cross-validate! Use the train function in caret to fit your 2 models. Use one value of mtry (the same for both models). Caret will return a re-sampled estimate of RMSE and $R^2$.

See page 3 of the caret vignette (also in the full reference manual)

• I agree except for your omission that OP should use R-squared adjusted, as this is what it is designed for. The models have same Y but different predictor sets, so R-squared adjusted is needed to penalize for the difference in spent degrees of freedom/model complexity. The other problem I see is potentially a sample size issue; OP would require a MASSIVE sample size to have any sort of real hope this isn't just overfitting, even with CV. – LSC Oct 27 '19 at 8:59
• @LSC I would say that if of holds out an additional dataset (not used for selection), they could use R2 without adjustment, as its just normalized RMSE. I agree that they'd need a pretty large dataset to do this. – Zach Jan 3 at 19:27
• new data or old, all else constant, putting 100 terms in a model vs 25 will have a higher unadjusted R-squared when compared to the R-squared from the 25 predictor model. This is easily seen when computing the sum of squared errors between the two (more terms is lower SSE, all else the same). I think people often forget that more terms will never decrease R-squared, but if they stink relative to their value, they can decreased adjusted R-squared which is a better measure to use to look at bang for the buck. – LSC Jan 3 at 20:44
• "more terms will never decrease R-squared" <- This statement is false, when computing R2 out of sample. – Zach Jan 3 at 21:05
• As I mentioned, I am pointing it out because many people asking these questions don't remember this idea in the first place. I think understanding why this is true in sample and may not hold on new data is important to understand what and how something is optimized in the original equation. – LSC Jan 4 at 10:24

I agree with Zach that the best idea is to cross-validate both models and then compare the $R^2$s, for instance by collecting values from each fold and comparing the resulting vectors with Wilcoxon test (paired for k-fold, unpaired for random CV).

The side option is to use all relevant feature selection, what would told you which attributes have a chance to be significantly useful for classification -- thus weather those expensive attributes are worth their price. It can be done for instance with a RF wrapper, Boruta.

• for a paired test, I assume that the 2 models are fit on the same folds? So that the k rows of the matrx are the folds and the two columns are model 1 and model 2? – B_Miner Aug 5 '11 at 17:55
• @B_Miner Precisely. – user88 Aug 5 '11 at 19:02

You may want to think in terms of practical significance rather than statistical significance (or both). With enough data you can find things signifcant statistically that will have no real impact on your usage. I remember analyzing a model one time where the 5-way interactions were statistically significant, but when the predictions from the model including everything up to the 5-way interactions were compared to the predictions from a model including only 2-way interactions and main effects, the biggest difference was less than 1 person (the response was number of people and all interesting values were away from 0). So the added complexity was not worth it. So look at the differences in your predictions to see if the differences are enough to justify the extra cost, if not then why bother even looking for the statistical significance? If the differences are big enough to justify the cost if they are real, then I second the other sugestions of using cross validation.

One option would be to create a confidence interval for the mean squared error. I would use the mean squared error instead of $R^2$ since the denominator is the same for both models. The paper by Dudoit and van der Laan (article and working paper) provides a general theorem for the construction of a confidence interval for any risk estimator. Using the example from the iris data, here is some R code creating a 95% confidence interval using the method:

library(randomForest)
data(iris)
set.seed(42)

# split the data into training and testing sets
index <- 1:nrow(iris)
trainindex <- sample(index, trunc(length(index)/2))
trainset <- iris[trainindex, ]
testset <- iris[-trainindex, ]

# with species
model1 <- randomForest(Sepal.Length ~ Sepal.Width + Petal.Length +
Petal.Width + Species, data = trainset)
# without species
model2 <- randomForest(Sepal.Length ~ Sepal.Width + Petal.Length +
Petal.Width, data = trainset)

pred1 <- predict(model1, testset[, -1])
pred2 <- predict(model2, testset[, -1])

y <- testset[, 1]
n <- length(y)

# psi is the mean squared prediction error (MSPE) estimate
# sigma2 is the estimate of the variance of the MSPE
psi1 <- mean((y - pred1)^2)
sigma21 <- 1/n * var((y - pred1)^2)
# 95% CI:
c(psi1 - 1.96 * sqrt(sigma21), psi1, psi1 + 1.96 * sqrt(sigma21))

psi2 <- mean((y - pred2)^2)
sigma22 <- 1/n * var((y - pred2)^2)
# 95% CI:
c(psi2 - 1.96 * sqrt(sigma22), psi2, psi2 + 1.96 * sqrt(sigma22))


The method can also be extended to work within cross-validation (not just sample-split as shown above).

Since you're already using randomForest after you cross-validate you might emit the chosen fit's computation of the predictor importance values.

> require(randomForest)
> rf.fit = randomForest(Species~.,data=iris,importance=TRUE)
> rf.fit\$importance
setosa   versicolor   virginica MeanDecreaseAccuracy MeanDecreaseGini
Sepal.Length 0.036340893  0.021013369 0.032345037          0.030708732         9.444598
Sepal.Width  0.005399468 -0.002131412 0.007499143          0.003577089         2.046650
Petal.Length 0.319872296  0.297426025 0.290278930          0.299795555        42.494972
Petal.Width  0.343995456  0.309455331 0.277644128          0.307843300        45.286720


I see this question has been asked long time ago; however, no answer points out to the significant shortcomings and misunderstandings in the question yet.