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

  • $\begingroup$ 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? $\endgroup$ – B_Miner Aug 5 '11 at 17:55
  • $\begingroup$ @B_Miner Precisely. $\endgroup$ – 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:


# 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.