Why does the rank order of models differ for R squared and RMSE?

Question

I am comparing $R^2$ and RMSE of different models. Interestingly, the rank ordering of the models with respect to $-R^2$ and RMSE is different and I do not understand why.

Here is an example in R:

library(caret) 

set.seed(0)
d<-SLC14_1(n=1000)

folds<-createMultiFolds(d$y,k=10,times=1)
tc<-trainControl(index=folds,returnResamp="all")
t1<-train(y~.,data=d,method="glmnet",trControl=tc) 
order(t1$results$RMSE)==order(-t1$results$Rsquared)

Output:

[1]  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE

Thus, the order if different for $-R^2$ suqared and $RMSE$.

The question is, why.

Let $SS_{res}$ be the sum of squared residuals $\sum(y_i-f_i)^2$.

$RMSE$ is defined as $\sqrt{SS_{res}/n}$.

$R^2$ is defined as $1-SS_{res}/SS_{tot}$ where $SS_{tot}$ is $\sum(y_i-\overline{y})^2$.

Since $SS_{res}=n*(RMSE)^2$, we can write $R^2$ as $1-n*(RMSE)^2/SS_{tot}$. Since $n$ and $SS_{tot}$ are constant and the same for all models, $-R^2$ and $RMSE$ should strictly positively related. However, they are not since the ranking order is in practice not identical (see example code).

What is wrong with my argument?

Answer 1

In caret, the calculation for results$RMSE and results$Rsquared is not as simple as what you've indicated. They are in fact the average of RMSE and $R^2$ over the ten holdout sets.

To confirm this, run the summary:

> t1
glmnet 

1000 samples
  20 predictors

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 900, 900, 900, 900, 900, 900, ... 
Resampling results across tuning parameters:

  alpha  lambda      RMSE      Rsquared 
  0.10   0.01065054  17.93931  0.1655746
  0.10   0.10650539  17.93720  0.1656599
  0.10   1.06505391  17.89291  0.1678166
  0.55   0.01065054  17.93838  0.1657046
  0.55   0.10650539  17.91755  0.1668356
  0.55   1.06505391  17.84962  0.1731936
  1.00   0.01065054  17.93824  0.1657245
  1.00   0.10650539  17.90045  0.1678998
  1.00   1.06505391  17.92535  0.1710923

RMSE was used to select the optimal model using  the smallest value.
The final values used for the model were alpha = 0.55 and lambda = 1.065054.

For the optimal parameter combination alpha = 0.55 and lambda = 1.065054 the performance on each held-out set is seen in the object t1$resample:

> t1$resample
       RMSE   Rsquared Resample
1  18.42848 0.04479504   Fold05
2  21.17820 0.10500276   Fold08
3  18.27933 0.20858027   Fold04
4  17.31308 0.19080079   Fold07
5  16.60865 0.21812706   Fold10
6  20.07291 0.18737052   Fold02
7  16.48082 0.24041654   Fold03
8  17.18363 0.18379930   Fold06
9  17.29819 0.13669866   Fold09
10 15.65289 0.21634546   Fold01

(Needless to say, the RMSE and Rsquared seen above are evaluated on different CV folds, so they don't rank order the same.) If you average these columns, you'll get:

> mean(t1$resample$RMSE)
[1] 17.84962
> mean(t1$resample$Rsquared)
[1] 0.1731936

...which are the same as the RMSE and Rsquared numbers seen in row 6 of the summary.

---

EDIT: Why does averaging over folds disrupt the rank ordering? Suppose we have split the data into $F$ folds, and we are considering $C$ tuning combinations. For each combo $c$ and held-out fold $f$, the relationship between the $R^2$ and MSE calculated on fold $f$ is: $$\operatorname{Rsquared}(c,f)=1-\frac{\operatorname{MSE}(c,f)}{\operatorname{Var}(f)},\tag1 $$ where $\operatorname{Var}(f)$ is shorthand for the variance of the observed responses in fold $f$. It is certainly true that for a given $f$, if we average over all $c$ then the monotonic relationship between $R^2$ and MSE is preserved, since by linearity: $$\frac1C\sum_c\operatorname{Rsquared}(c,f)=1-\frac{\frac1C\sum_c\operatorname{MSE}(c,f)}{\operatorname{Var}(f)}.\tag2 $$ However, if we average (1) over all $f$ we cannot assert a similar statement, since the denominator $\operatorname{Var}(f)$, which varies with the fold being held out, gets in the way: $$\frac1F\sum_f\operatorname{Rsquared}(c,f)=1-\frac1F\sum_f\left(\frac{\operatorname{MSE}(c,f)}{\operatorname{Var}(f)}\right).\tag3 $$ The RHS of (3) cannot be simplified further to reveal a monotonic relationship between the average $R^2$ over all folds and the average MSE over all folds.

Since MSE is the square of RMSE, the relationship between fold-averaged $R^2$ and fold-averaged RMSE is even less direct. Indeed, for any given fold, there is not even an analog for (2) between combo-averaged $R^2$ and combo-averaged RMSE.

Answer 2

It's because caret calculates R-squared differently than you are. See the answer to this question: How caret calculates R Squared.

To see it in your code,

library(caret) 

set.seed(0)
d<-SLC14_1(n=1000)

folds<-createMultiFolds(d$y,k=10,times=1)
tc<-trainControl(index=folds,returnResamp="all",
             savePredictions = TRUE) # New option 
t1<-train(y~.,data=d,method="glmnet",trControl=tc) 
order(t1$results$RMSE)==order(-t1$results$Rsquared)

library(data.table)
preds <- data.table(t1$pred)
preds[, overall_mean := mean(obs), by = .(lambda, alpha, Resample)]

sum_sq <- preds[, .(SS_res = sum((obs - pred)^2),
                SS_tot = sum((obs - overall_mean)^2),
                n = .N,
                var = var(obs),
                Rsquared_corr = cor(obs, pred)^2),
            by = .(lambda, alpha, Resample)]
sum_sq <- sum_sq[, ':=' (RMSE_Julian = sqrt(SS_res / n),
                         Rsquared_Julian = 1 - (SS_res/SS_tot),
                         Rsquared_traditional = 1 - (SS_res/ ((n-1)*var) ))]
sum_sq <- merge(sum_sq, t1$resample, by = c("lambda", "alpha", "Resample"))
head(sum_sq)

Note the savePredictions = TRUE in the call to traincontrol(). In the final dataset, sum_sq, you can see your result, Rsquared_Julian matches Rsquared_traditional, but these don't match Rsquared_corr which does match the R-squared from caret, Rsquared.

Also in your question, you assume n and SS_tot are constant, but that only holds true for a fold, not across all the cross-validations.

Answer 3

@grand_chat has the correct maths, I'm just growing in a comparative example to help illustrate what the issue is in different terms that will hopefully help understanding.

We're working with fractional terms here, similar to say miles per gallon. If we average mpg over set units of time we get very different results compared to over set units of fuel or distance.

If we travel 10 minutes at 50 mph achieving 50mpg then 10 minutes at 60 achieving 30 mpg and we then want to calculate the average fuel efficiency for the journey.

Time based average (with one minute representing a unit of time) is $(50*10+30*10) /20 = 40 mpg$

But the distance we travel is $50/6 + 60/6 = 18.33 miles $ given that ten minutes is 1/6 th of an hour

The fuel we use is $(50/6)/50 +(60/6)/30=1/6+2/6 = 1/2 gallon$

This means our average mpg is in fact $18.33/(1/2)=36.66$

Because the total variance is different in every fold you would need to account for this in the averaging to maintain the monotonic relationship. Since it is present in the R2 calculation but not the RMSE then you can get rank switching by not accounting for the total variance in each fold