0

For a prediction problem I am working on I want to calculate the variance in the data explained by the linear effects of my linear mixed model. To evaluate my predictive performance I plan on using five-fold cross-validation.

The common approach to this problem I see on the internet is to a marginal R. This can be done with the r.squaredGLMM function from the MuMln package (https://cran.r-project.org/web/packages/MuMIn/MuMIn.pdf) based on work by Nakagawa et al. (2017). As I understand it, his formula divides the variances explained by the fixed effect by the total model variance:

enter image description here

However, I don't know how to make this formula work in a cross-validation setting as I don't know how to calculate the specific variances for fixed, random, and error in the test dataset. My current solution is to apply the OLS formula (1 - var(fitted)/var(observed)), but this does not match the maximum likelihood approach used in LMMs. It also leads to a negative outcome in some situations.

My question in two stepped:

  1. Is there a way to calculate the required variances for the Nakagawa formula?
  2. If not, can the 'classic' R2 be applied in this setting?
jrcalabrese
  • 2,184
  • 3
  • 10
  • 30
Pieter van der Veere
  • 3
  • 4
  • If you want to fit a predictive model can you just use raw mean-squared-error on the test set as your criterion? i.e., why do you need to mess with R^2 measures? (Can you provide links to the internet examples you see using marginal R^2 measures?) – Ben Bolker Mar 21 '23 at 13:37
  • Using a raw RMSE is my current approach for evaluating my own models, but it makes a comparison with the results from for example Dansson 2021 (doi: 10.1186/s13195-021-00886-5) difficult as express their predictive performance using an R^2. Calculating an marginal R^2 would be the closest analogue as I understand it. – Pieter van der Veere Jun 15 '23 at 09:00

1 Answers1

0

If you are after Nakagawa R2 in mixed models then, library(performance) is very helpful. Most people use library(lme4) for linear mixed models, lmer(), and for generalised ones, glmer(). On lmer() models, performance() defaults to using Nakagawa methods. Then you can use performance_cv() for various cross validation methods.

Example:

library(performance)
library(lme4)
data("iris")
test1 <- lmer(Sepal.Width ~ Petal.Width + (1|Species), data = iris, REML = FALSE)
test1
performance(test1)
performance_cv(test1, method = "k_fold", k = 5, stack = FALSE)

Should give you results:

# Indices of model performance

AIC    |   AICc |     BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
---------------------------------------------------------------------------
91.941 | 92.217 | 103.983 |      0.910 |      0.318 | 0.869 | 0.297 | 0.300
> performance_cv(test1, method = "k_fold", k = 5, stack = FALSE)
# Cross-validation performance (5-fold method)

MSE   | MSE_SE | RMSE | RMSE_SE |   R2 | R2_SE
----------------------------------------------
0.095 |  0.036 |  0.3 |   0.059 | 0.46 |  0.24

That should help you get started, however I am having a few issues dissecting some of the nuances of these results from the CV, and have posted on StackExchange about them here, Nakagawa R2 & Cross Validation issues, if anyone can help with that.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
MrSwaggins
  • 87
  • 8
  • My apologies for the delay in responding. This does solve the problem I had a the time of posting. Thank you! I have now extended the analysis by adding an external validation dataset for my prediction model. If I read the package manual correctly, it cannot calculate an R2 using a different dataset as reference. Do you have a suggestion for calculating an R2 for my prediction model on an external dataset – Pieter van der Veere Jun 15 '23 at 08:56
  • I'm not exactly sure what you mean by "external", but it might be this. In `performance_cv` you can add an argument `data = yourExternalData`, "data Optional. A data frame containing the same variables as model that will be used as the cross-validation sample.", so something like this: `performance_cv(test1, data = yourExternalData)` which I think would compare them. I haven't done this before so your mileage may vary. – MrSwaggins Jun 26 '23 at 02:41
  • This does provide an out of sample R2, but it does not calculate an R2 using the Nakagawa formula (if that is even possible). A conventional R2 leads to very different results for me. Depending on the setting this might work for someone else. – Pieter van der Veere Jul 07 '23 at 10:35
  • The cross-validation `performance_cv()` will give you the CV (Nakagawa) marginal R2, which is the R2 of the fixed factors only. In the example I gave, R2(marg.) 0.318, is compared to cross-validated R2 0.46, which is 0.142 difference. Depending on your data, "big" differences may point to problems with over-fitting, as I understand it. – MrSwaggins Jul 09 '23 at 23:55