# 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. Any help is appreciated.

-

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 R2.

See section 5 of the caret vignette. (page 39)

-

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. –  mbq Aug 5 '11 at 19:02

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

-