How to extract the goodness-of-fit statistics from lmer() model for msummary from modelsummary package

Question

I am using lmerTest::lmer() to perform linear regression with repeated measures data.

My model contains a fixed effect (factor with 5 levels) and a random effect (subject):

library(lmerTest) 
model_lm  <- lmer(likertscore ~ task.f + (1 | subject), data = df_long)  

I would like to include the total number of observations, the number of subjects, total R^2, and the R^2 of the fixed effects in the regression table which I generate with modelsummary().

I tried to extract these and build a gof_map as described by the author of the package but did not succeed. Below my model output from lmerTest::lmer() the performance measures obtained:

Linear mixed model fit by REML ['lmerModLmerTest']
Formula: likertscore ~ factor + (1 | subject)
   Data: df_long
REML criterion at convergence: 6674.915
Random effects:
 Groups   Name        Std.Dev.
 subject  (Intercept) 1.076   
 Residual             1.514   
Number of obs: 1715, groups:  subject, 245
Fixed Effects:
                      (Intercept)                         factor1                         factor2  
                           3.8262                             1.5988                             0.3388  
                      factor3                             factor4                         factor5  
                          -0.7224                            -0.1061                            -1.1102  

library("performance")
performance::model_performance(my_model)

# Indices of model performance

AIC     |     BIC | R2 (cond.) | R2 (marg.) |  ICC | RMSE | Sigma
-----------------------------------------------------------------
6692.91 | 6741.94 |       0.46 |       0.18 | 0.34 | 1.42 |  1.51

Answer 1

The problem is that one of your statistics is not available by default in glance or performance, which means that you will need to do a bit of legwork to customize the output.

First, we load the libraries and estimate the model:

library(modelsummary)
library(lmerTest) 
mod  <- lmer(mpg ~ hp + (1 | cyl), data = mtcars) 

Then, we check what goodness-of-fit statistics are available out-of-the-box using the get_gof function from the modelsummary package:

get_gof(mod)
#>        aic      bic r2.conditional r2.marginal       icc     rmse    sigma nobs
#> 1 181.8949 187.7578      0.6744743   0.1432201 0.6200592 2.957141 3.149127   32

You'll notice that there is no N (subject) statistic there, so we need to add it manually. One way to do this in a replicable way is to leverage the glance_custom mechanism described in the modelsummary documentation. To do this, we need to know what the class of our model is:

class(mod)[1]
#> [1] "lmerModLmerTest"

Then, we need to define a method for this class name. This method should be called glance_custom.CLASSNAME. In lmerModLmerTest models, the number of groups can be retrieved by getting the ngrps object in the summary. So we do this:

glance_custom.lmerModLmerTest <- function(x, ...) {
  s <- summary(x)
  out <- data.frame(ngrps = s$ngrps)
  out
}

Finally, we use the gof_map argument to format the result how you want it:

gm <- list(
  list(raw = "nobs", clean = "N", fmt = 0),
  list(raw = "ngrps", clean = "N (subjects)", fmt = 0),
  list(raw = "r2.conditional", clean = "R2 (conditional)", fmt = 0),
  list(raw = "r2.marginal", clean = "R2 (marginal)", fmt = 0),
  list(raw = "aic", clean = "AIC", fmt = 3)
)

modelsummary(mod, gof_map = gm)

	Model 1
(Intercept)	24.708
	(3.132)
hp	-0.030
	(0.015)
N	32
N (subjects)	3
R2 (conditional)	1
R2 (marginal)	0
AIC	181.895