How to compute R-squared value when doing cross-validation?

Question

I am using multiple linear regression with a data set of 72 variables and using 5-fold cross validation to evaluate the model.

I am unsure what values I need to look at to understand the validation of the model. Is it the averaged R squared value of the 5 models compared to the R squared value of the original data set? In my understanding, the average R squared value of the sampled data needs to be within 2% of the R squared value in the original data set. Is that right? Or are there any other results I should be looking at?

Answer 1

It is neither of them. Calculate mean square error and variance of each group and use formula $R^2 = 1 - \frac{\mathbb{E}(y - \hat{y})^2}{\mathbb{V}({y})}$ to get R^2 for each fold. Report mean and standard error of the out-of-sample R^2.

Please also have a look at this discussion. There are a lots of examples on the web, specifically R codes where $R^2$ is calculated by stacking together results of cross-validation folds and reporting $R^2$ between this chimeric vector and observed outcome variable y. However answers and comments in the discussion above and this paper by Kvålseth, which predates wide adoption of cross-validation technique, strongly recommends to use formula $R^2 = 1 - \frac{\mathbb{E}(y - \hat{y})^2}{\mathbb{V}({y})}$ in general case.

There are several things which might go wrong with the practice of (1) stacking and (2) correlating predictions.

1. Consider observed values of y in the test set: c(1,2,3,4) and prediction: c(8, 6, 4, 2). Clearly prediction is anti-correlated with the observed value, but you will be reporting perfect correlation $R^2 = 1.0$.

2. Consider a predictor that returns a vector which is a replicated mean of the train points of y. Now imagine that you sorted y and before splitting into cross-validation (CV) folds. You split without shuffling, e.g. in 4-fold CV on 16 samples you have following fold ID labels of the sorted y:

foldid = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4)
y = c(0.09, 0.2, 0.22, 0.24, 0.34, 0.42, 0.44, 0.45, 0.45, 0.47, 0.55, 0.63, 0.78, 0.85, 0.92, 1)

When you split you sorted y points, the mean of the train set will anti-correlate with the mean of the test set, so you get a low negative Pearson $R$. Now you calculate a stacked $R^2$ and you get a pretty high value, though your predictors are just noise and the prediction is based on the mean of the seen y. See figure below for 10-fold CV

Answer 2

Update: Revisiting my 'youthful' answer, I agree, this stitching approach is not the right way to compute R-squared metric. Stitching may be useful for visual inspection of residuals. I leave answer as is, as it is mentioned in other answers.

@Is it the averaged R squared value of the 5 models?

-No, it is computed as seen below. You predict k-fold observations, stitch them together to a ordered vector where obs#1 is first and obs#last is last. Calculate then the squared pearson product moment correlation (R²) of this k-fold prediction vector to the response vector(y). CV-correlation to response(y) is lower than a direct MLR-fit. In the example below R²(CV) = .63 and R²(direct fit)=.82. This would suggests simple MLR here is slightly overfitted, and if this bothers you could try to do somewhat better with PLS, ridge-regression or PCR. I have not heard of any 2% rule.

 library(foreach)
 obs=250
 vars=72
 nfolds=5

 #a test data set
 X = data.frame(replicate(vars,rnorm(obs)))
 true.coefs = runif(vars,-1,1)
 y_signal = apply(t(t(X) * true.coefs),1,sum)
 y_noise = rnorm(obs,sd=sd(y_signal)*0.5)
 y = y_signal + y_noise

 #split obs randomly in nfold partitions
 folds = split(sample(obs),1:nfolds)
 #run nfold loops, train, predict..
 #use cbind to stich together predictions of each test set to one
 test.preds = foreach(i = folds,.combine=cbind) %do% {
   Data.train = data.frame(X=X[-i,],y=y[-i])
   Data.test  = data.frame(X=X[i ,],y=y[ i])
   lmf = lm(y~.,Data.train)  
   test.pred = rep(0,obs)
   test.pred[i] = predict(lmf,Data.test)
   return(test.pred)
 }
 CVpreds = apply(test.preds,1,sum)

 cat(nfolds,"-fold CV, pearson R^2=",cor(CVpreds,y)^2,sep="")
 cat("simple MLR fit, pearson R^2=",cor(lm(y~.,data.frame(y,X))$fit,y)^2,sep="")