Model diagnosis in GLMM model of binary outcome variable

Question

I am trying to do a research article on morbidity, which is categorical, i.e., yes or no. I also consider fixed and random effects, given that all of the fixed effect covariates are categorical. Here is the model.

twoRIM <- glmer(morbidity ~ co_fuel + wealth + san_class + residencey + education + marital + age + parity + Electricity_avaliablity + education * age + (1|region) + (1|participant_id), family=binomial, data=M)

Hence, I try to compare these two random intercept models to a model that has a single random intercept mode. Based on AIC, the model "twoRIM" is the best model because it has the lowest AIC. Then I tried to check the model diagnostics of this model, model "twoRIM", by running the following R codes. I got the following results or plots.

scatter.smooth(fitted(twoRIM), sqrt(abs(resid(twoRIM))), col=6)

qqline(resid(twoRIM))

plot(twoRIM)

qqnorm(resid(twoRIM),main="Residual normal plot",col=4,adj=0.1)

qqnorm(ranef(twoRIM)$"region"
[[1]],main="Regional level random effects",col=2,adj=0.1)

qqnorm(ranef(twoRIM)$"participant_id"[[1]],main="Cluster level random effects",col=6,adj=0.1)

plot(fitted(twoRIM),resid(twoRIM),col=4)

qqnorm(resid(twoRIM))

The plots I got — the residual and other related plots — are quite different from other forms of model diagnosis I knew before, and I faced a bit more difficulty interpreting and understanding the nuance of these plots. Maybe I took the codes and arguments wrongly?

Answer 1

You have a lot of plots here, but let me just save you some time with these. The typical residual plots are here not going to be helpful because the outcome of a logistic regression can only span the range of $[0,1]$. As such, this is why you get horizontal asymptotes and dual lines in your residuals. What's more, the residuals from GLMMs, as opposed to Gaussian mixed models, behave in fairly bizarre ways.

A much better solution is to get residuals based on simulated data from the model. A package that allows for this is the DHARMa package, which has a very useful vignette on this topic here. I recommend running your residual analysis on these types of residuals instead, as they will be a lot more useful to you.