Manually calculated $R^2$ doesn't match up with randomForest() $R^2$ for testing new data

Question

I know this is a fairly specific R question, but I may be thinking about proportion variance explained, $R^2$, incorrectly. Here goes.

I'm trying to use the R package randomForest. I have some training data and testing data. When I fit a random forest model, the randomForest function allows you to input new testing data to test. It then tells you the percentage of variance explained in this new data. When I look at this, I get one number.

When I use the predict() function to predict the outcome value of the testing data based on the model fit from the training data, and I take the squared correlation coefficient between these values and the actual outcome values for the testing data, I get a different number. These values don't match up.

Here's some R code to demonstrate the problem.

# use the built in iris data
data(iris)

#load the randomForest library
library(randomForest)

# 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, ]

# fit a model to the training set (column 1, Sepal.Length, will be the outcome)
set.seed(42)
model <- randomForest(x=trainset[ ,-1],y=trainset[ ,1])

# predict values for the testing set (the first column is the outcome, leave it out)
predicted <- predict(model, testset[ ,-1])

# what's the squared correlation coefficient between predicted and actual values?
cor(predicted, testset[, 1])^2

# now, refit the model using built-in x.test and y.test
set.seed(42)
randomForest(x=trainset[ ,-1], y=trainset[ ,1], xtest=testset[ ,-1], ytest=testset[ ,1])

Answer 1

The reason that the $R^2$ values are not matching is because randomForest is reporting variation explained as opposed to variance explained. I think this is a common misunderstanding about $R^2$ that is perpetuated in textbooks. I even mentioned this on another thread the other day. If you want an example, see the (otherwise quite good) textbook Seber and Lee, Linear Regression Analysis, 2nd. ed.

A general definition for $R^2$ is $$ R^2 = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2} . $$

That is, we compute the mean-squared error, divide it by the variance of the original observations and then subtract this from one. (Note that if your predictions are really bad, this value can go negative.)

Now, what happens with linear regression (with an intercept term!) is that the average value of the $\hat{y}_i$'s matches $\bar{y}$. Furthermore, the residual vector $y - \hat{y}$ is orthogonal to the vector of fitted values $\hat{y}$. When you put these two things together, then the definition reduces to the one that is more commonly encountered, i.e., $$ R^2_{\mathrm{LR}} = \mathrm{Corr}(y,\hat{y})^2 . $$ (I've used the subscripts $\mathrm{LR}$ in $R^2_{\mathrm{LR}}$ to indicate linear regression.)

The randomForest call is using the first definition, so if you do

   > y <- testset[,1]
   > 1 - sum((y-predicted)^2)/sum((y-mean(y))^2)

you'll see that the answers match.