How to compare linear mixed models that do not satisfy the assumption of Gaussian errors?

Question

I have fitted two linear mixed models, A and B, to a data set using lme4::lme: model B has one more predictor variable than model A.

Using check_heteroscedasticity and check_normality from performance package, I have confirmed that the residuals have homogenous variance in both models, however, both models have failed to satisfy the assumption of Gaussian residuals even after (probably) appropriate transformation (as suggested by MASS::boxcox).

How do I go on about deciding which model I should trust if they don't satisfy the normality assumption? If model A is considered a null model, it consists of all the relevant predictor variables (and appropriately transformed) except one based on the experiment through which the data were collected.

Answer 1

You should try to model the regression off the error term that it belongs with. If the regression residuals are not normally distributed, this could affect the accuracy of it's predictions and make it a less useful model. Instead, you should check which distribution the errors look like and try to fit the regression accordingly. I show below with a simulated example. First, I load the requisite libraries in R and set a seed for reproducibility:

#### Set Seed and Load Libraries ####
library(lmerTest)
library(performance)
set.seed(123)

Then I simulate some normally distributed IVs with a binary response:

#### Simulate Data ####
y <- rbinom(n=1000,
            size=1,
            prob=.5)

x1 <- y*.5 + rnorm(n=1000)
x2 <- y*.5 + rnorm(n=1000)

re <- c("school1",
        "school2",
        "school3",
        "school4",
        "school5")

school <- sample(re, 
       size = 1000,
       replace = T)

df <- data.frame(x1,x2,y,school)

After I fit the model and check what it's residuals look like with check_distribution like you did:

#### Fit Model ####
fit.1 <- lmer(y ~ x1 + x2 + (1|school),
            data = df)

#### Check Distribution ####
plot(check_distribution(fit.1))

One can clearly see that the residuals are bimodal and thus probably require a logistic GLMM instead of a Gaussian version. This can be easily achieved with glmer, which only differs in the family argument you specify.

#### Fit to Binomial Family ####
fit.2 <- glmer(y ~ x1 + x2 + (1|school),
             data = df,
             family = binomial)

We can check the performance of the two models and see how they compare at predicting the outcome using the code compare_performance(fit.1,fit.2), which gives this output:

# Comparison of Model Performance Indices

Name  |           Model |  AIC (weights) | AICc (weights) |  BIC (weights) | R2 (cond.) | R2 (marg.) |       ICC |  RMSE | Sigma | Log_loss | Score_log | Score_spherical
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
fit.1 | lmerModLmerTest | 1359.7 (<.001) | 1359.8 (<.001) | 1384.3 (<.001) |      0.098 |      0.097 |     0.002 | 0.475 | 0.476 |          |           |                
fit.2 |        glmerMod | 1292.3 (>.999) | 1292.3 (>.999) | 1311.9 (>.999) |      0.127 |      0.127 | 5.122e-04 | 0.475 | 1.000 |    0.642 |      -Inf |           0.001

The AIC, BIC, and pseudo R2 values are all unsurprisingly more attractive in the logistic model. I'm assuming your two models share similarly distributed outcomes, so you could compare them against teach other with the same code. You could also use an LRT test to see which model is preferable:

#### Drop Predictors ####
fit.3 <- glmer(y ~ x1 + (1|school),
               data = df,
               family = binomial)

#### LRT Test ####
anova(fit.2,fit.3)

If for some reason the the two models do not have the same distribution of errors (probably unlikely), you can still compare their metrics in a similar way, though if the distributions varied dramatically from each other that would be something you would need to explain in your results.

The residuals in GLMMs often act fairly bizarre, so checking them at this point isn't recommended. For this approach, it would be better to use simulated residuals from the DHARMa package, which has a useful vignette here.