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

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?

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 Kvå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

• Here's an editing tip: use  for code inside a paragraph, but if you have a line (or several) of code you just need to start the like with enough spaces and it will display as code in a nicer way. – Silverfish Apr 15 '16 at 7:16
• I think your formula is missing the square. – Tal Galili Jun 12 '19 at 12:36
• I manually changed the first formula to be (what I suspect is) the correct formula. – Tal Galili Jun 13 '19 at 7:53
• I agree overall. I guess second reason is a little far fetched as it assumes the user performed an incorrect k-fold CV. CV-stitching is not susceptible to sorted targets, as k-fold CV uses random partition/shuffling. Random partition makes any pre-sorting irrelevant. However, CV-stitching should mainly be used for visual inspection. – Soren Havelund Welling Apr 11 at 11:21

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="")

• I have seen this way in many R codes, but no explanation whether and how it actually characterizes the performance of an estimator. sklearn` actually does estimate performance on each subset, and it seems much more logical to me. Is there a reference or it is just R folklore? – Dima Lituiev Apr 15 '16 at 1:38