How to motivate the definition of $R^2$ in `sklearn.metrics.r2_score`?

Question

TLDR: What motivates the definition of $R^2$ in the Python function sklearn.metrics.r2_score?

DETAILS

The Python machine learning package sklearn implements an $R^2$ using the following formula.

$$ R^2=1-\left(\dfrac{ \overset{N}{\underset{i=1}{\sum}}\left( y_i-\hat y_i \right)^2 }{ \overset{N}{\underset{i=1}{\sum}}\left( y_i-\bar y \right)^2 }\right) $$

This is fine for in-sample assessments of $R^2$, but when it comes to out-of-sample assessments, the definition seems unmotivated and lacking in meaning.

In simple linear regression (in-sample), there are multiple equivalent definitions of $R^2$.

1. Squared Pearson correlation between the feature and the outcome
2. Squared Pearson correlation between the outcome and the predictions
3. Proportion of variance explained
4. Comparison of the square loss incurred by the model to the square loss incurred by an intercept-only baseline model that always predicts the sample mean, $\bar y$

(There may be more, and perhaps it is me only thinking of these four that prevents me from seeing what this Python function means.)

In regression models with multiple features, definition 1 does not make sense, so it is not a contender for a generalization of $R^2$ to out-of-sample assessments in complicated situations that are likely to have multiple features.

Out-of-sample, squaring the correlation misses prediction bias. For instance, $y=(1,2,3)$ is perfectly correlated with $\hat y=(11,12,13)$, yet $\hat y$ is composed of terrible predictions of $y$. While this kind of bias cannot happen in-sample in ordinary least squares linear regression, out-of-sample, all bets are off (and that’s without even getting into what could go wrong in fancier regression models like support vectors or neural networks). Consequently, definition 2 is out as a generalization of $R^2$ to out-of-sample assessments.

Except in special situations, the usual methods for calculating $R^2$ do not correspond to the proportion of variance explained, and if we twist the definition to force it, we usually lose the connection to the square loss that we aim to minimize. Thus, definition 3 does not seem like a good contender for a generalization of $R^2$.

Finally, definition 4 makes sense. We have some kind of baseline model (naïvely predict $\bar y$ every time, always using the marginal mean as our guess of the conditional mean) and compare our predictions to the predictions made by that baseline model. To draw an analogy to flipping a coin, if someone guesses which side will land up and gets correct predictions less than half the time, that person is a poor predictor. If they are right more than half the time, they are at least improving somewhat upon the naïve “Gee, I don’t know how it’ll land, so I guess I’ll just say heads every time (or tails, or alternate between the two) and get it right about half the time.”

(Put bluntly, why pay a data scientist or statistician a lot of money to predict conditional means when you could do better just by predicting one number every time that you can calculate by AVERAGE(A:A) in Excel?)

In-sample, I am totally on board with the Python function. Out-of-sample, I have a problem. In the $R^2$ formula above, the Python implementation uses the $\bar y$ from the given data. That is, if you input an out-of-sample $y$ and the corresponding predictions $\hat y$, the function uses the out-of-sample $\bar y$ in the denominator.

This corresponds to comparing your model to a model that predicts the out-of-sample mean every time. However, we cannot have access to such a model, since it would require knowledge of the true out-of-sample values.

QUESTION: What motivates this definition that is used in sklearn?

An argument about ease of use does not seem legitimate to me (even if that is why the function is defined this way), because the following is even easier to implement.

def r2_score():
    return(0.5)
    # or return(np.random.uniform(0, 1, 1))

While this is an extreme example, the function does not return useful information about the regression model. Consequently, making a function easy to use at the expense of returning unhelpful information seems like a poor design or a function that we should not use.

EDIT

I have referred to this question so many times yet never noticed until now that I never wrote explicitly what the competing notions of out-of-sample $R^2$ are. That ends now!

$$ R^2_{\text{out-of-sample, Dave}}= 1-\left(\dfrac{ \overset{N}{\underset{i=1}{\sum}}\left( y_i-\hat y_i \right)^2 }{ \overset{N}{\underset{i=1}{\sum}}\left( y_i-\bar y_{\text{in-sample}} \right)^2 }\right) $$$$ R^2_{\text{out-of-sample, scikit-learn}}= 1-\left(\dfrac{ \overset{N}{\underset{i=1}{\sum}}\left( y_i-\hat y_i \right)^2 }{ \overset{N}{\underset{i=1}{\sum}}\left( y_i-\bar y_{\text{out-of-sample}} \right)^2 }\right) $$

EDIT 2

My preferred calculation now has a reference in the literature (Hawinkel, Waegeman & Maere (2023)) explaining its superiority over that used by sklearn. One of the authors appears to be a Cross Validated contributor, too!

REFERENCE

Stijn Hawinkel, Willem Waegeman & Steven Maere (2023) Out-of-sample $R^2$: estimation and inference, The American Statistician, DOI: 10.1080/00031305.2023.2216252

Answer 1

1. Squared correlation between the feature and the outcome

That would be the case if you have a single feature and the model is linear regression.

2. Squared correlation between the outcome and the predictions

Same as above, but it will hold also if there are more features.

3. Proportion of variance explained

It tells us the proportion of the variance explained, but only for the linear regression.

4. Comparison of the square loss incurred by the model to the square loss incurred

Again, for linear regression using the formula from Scikit-learn is equivalent to the as we can decompose the squared error to TSS = ESS + RSS and get the equivalent formulation.

So in a sense, all the formulations are the same, just have varying degrees of generality.

As for calculating $R^2$ on the test set, you can check this thread and if you search through CrossValidated.com and Scikit-learn GitHub issues, you can find many discussions considering the Scikit-learn's choice of test set mean in the denominator as controversial. As you can learn from this discussion, one of the problems with using the train set mean for the test set $R^2$ is the API, where the metrics usually are defined as metric(y_true, y_pred) and the same interface can be used regardless if it is training or test set. It would also mean that during evaluation time you would need to have access to the training data, which may not be possible in some setups where there is a hard split of the data between training time and validation time.

Notice also that $R^2$ comes from statistics, where we are usually interested in in-sample metrics, so the derivations would as well regard the train set $R^2$.

Answer 2

4. Comparison of the square loss incurred by the model to the square loss incurred by a baseline model

Comparison of the loss seems a lot like the pseudo-$R^2$ value e.g.

$$R^2_{pseudo} = 1 - \frac{D_{null} - D_{model}}{D_{null}}$$

But with the deviance or loss equal to the sum of squared residuals and the null model the mean, then it becomes the same as the regular $R^2$.

In-sample, I am totally on board with the Python function. Out-of-sample, I have a problem. In the $R^2$ formula above, the Python implementation uses the $\hat{y}$ from the given data.

Possibly the problem stems from the use of $R^2$ as a measure for goodness of fit. But the $R^2$ value is not a goodness of fit measure. The value doesn't tell directly whether your model is a good fit or not; A perfect fit of the conditional distribution mean does not need to coincide with a $R^2=1$. Instead, it is a descriptive statistic that tells how large the variance in the noise/randomness is relative to the variance in the deterministic part. We see this in an alternative way to compute $R^2$

$$R^2 = \frac{SS_{model}}{SS_{model}+SS_{residuals}}$$

where $SS_{model} = \sum (\hat{y} - \bar{\hat{y}})^2$ and $SS_{residuals} = \sum (y - \hat{y})^2$ and $\bar{\hat{y}}$ is the mean of the modelled values.

This way of computing $R^2$ will be equivalent to the 'other $R^2$' if you have a linear model with an intercept. But it will be slightly different in other situations, for instance there is no occurrence of cases with negative values as in the question Why is R^2 negative in my multiple linear regression model in python? .