3

I am trying to access AIC, BIC , logLik and deviance data from a model summary of an HLM fitted using maximum likelihood (ML) in lme4::lmer, and combine with essentially the same model fitted with restricted maximum likelihood (REML). The structure of the objects returned from lmer and summary is a mess, and I am unable to find out where/how this data is stored.

[Update:] Based on the responses I've gotten, I have updated the code to reflect the progress made:

Code example:

# Least working example
library(lme4)
library(lmerTest)
df <- lme4::sleepstudy
names(df)
# Example model
model <- lmer(Reaction ~ (1|Subject), df, REML = TRUE)
information_criterion <- data.frame(
            "AIC" = AIC(model),
            "BIC" = BIC(model),
            "logLik" = logLik(model),
            "deviance" = deviance(model, REML=FALSE),
            "df.residual" = df.residual(model)
            )
mod_sum <- list(summary(model), information_criterion)
I essentially want to modify the output to resemble the output of summary if REML = FALSE (not working):
> mod_sum

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Reaction ~ (1 | Subject)
   Data: df

## Information criterion injected here: ##########################

     AIC      BIC   logLik deviance df.resid   # <-- THESE ARE THE LINES I WANT
  1916.5   1926.1   -955.3   1910.5      177   # <-- 

##################################################################

REML criterion at convergence: 1904.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4983 -0.5501 -0.1476  0.5123  3.3446 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 1278     35.75   
 Residual             1959     44.26   
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   298.51       9.05  17.00   32.98   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Robert Long
  • 5,722
  • 5
  • 29
  • 50
Pål Bjartan
  • 793
  • 1
  • 6
  • 18
  • 1
    I think this might be useful? https://stats.stackexchange.com/questions/131272/lme4-why-is-aic-no-longer-displayed-when-using-reml – sjp Jan 09 '22 at 14:59

1 Answers1

2

There are a few points:

  1. You have a typo here:
    m2sum[["information_criterion"]] <- summary(model1)$information_criterion

It should be m2_sum

  1. Instead of summary(model1)$information_criterion you can use:
     AIC(model1)

So, the following should work:

m2_sum[["information_criterion"]] <- AIC(model1)

Update following change to the OP.

This should work, although please see my last comment, because this may not be a sensible thing to do:

> m2_sum$AICtab <- m1_sum$AICtab
> m2_sum

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Reaction ~ (1 | Subject)
   Data: df

     AIC      BIC   logLik deviance df.resid 
  1916.5   1926.1   -955.3   1910.5      177 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.4983 -0.5501 -0.1476  0.5123  3.3446
Robert Long
  • 5,722
  • 5
  • 29
  • 50
  • Thanks for the feedback. I realized I can add all desired information criterion the same way, without actually fitting an ML-model. I'd be interested in how to add in the output as illustrated though. – Pål Bjartan Jan 09 '22 at 15:40
  • I don't know what you mean *"how to add in the output as illustrated though"* ? – Robert Long Jan 09 '22 at 15:43
  • You have an object `m2_sum` and you asked how to assign a value to the `information_criterion` slot of that object, and that's what I have answered. – Robert Long Jan 09 '22 at 15:51
  • Sorry for the ambiguity. I have updated the model based on your input, and I try show where to display the information. Essentially, I want the output similar to the output of `summary` if `REML = FALSE`. – Pål Bjartan Jan 09 '22 at 16:08
  • 2
    OK, I see what you want to do. Please note that this doesn't really make sense because a model fitted with REML does not really have a likelihood to compute AIC with; and this is a good thing because it helps to avoid the big mistake of comparing likelihoods from different models fitted with REML. If you have questions about this I would suggest asking on Cross Validated. I have updated my answer with a possible solution. – Robert Long Jan 09 '22 at 17:21