# 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. Commented Apr 15, 2016 at 7:16
• I think your formula is missing the square. Commented Jun 12, 2019 at 12:36
• I manually changed the first formula to be (what I suspect is) the correct formula. Commented Jun 13, 2019 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. Commented Apr 11, 2020 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? Commented Apr 15, 2016 at 1:38