Improving model fit and model comparison using AIC with logistic regression

Question

So I am running some hypothesis testing on a data set with a binary response (presence/absence data). My response variable is the use or lack of use of a given resource (0 or 1 respectively), over ~6 different resources. I am looking across populations using the glmer function in R with my individual as a random effect. I have my models set up as follows:

glmer(use~var1+var2+var3+var4+var5+(1|ID), family=binomial)
glmer(use~var1+var2+var3+var4+(1|ID), family=binomial) 
glmer(use~var1+var2+var4+var5+(1|ID), family=binomial) 
glmer(use~var1+var3+var4+var5+(1|ID), family=binomial) ## .. and so on. 
#NULL: 
glmer(use~{1|ID))

My end goal is to compare the various models to see what is the best fit using AIC, a common practice within my field.

I am being held up on two main questions:

1. in cases where my null model is the best fit what does this actually tell me outside of the fact that I am probably not catching the variables driving use within that population? Is there something I can do other than looking at different explanatory variables, and could it be an issue with the fact that I have too little use data for that group?

2. When you get warnings about model fit (see example below), what can be done outside of scaling variables to try and improve binomial model fit? Because it is a 0 or 1 I assume that I cannot transform the response data, and I have already scaled my environmental variables that had a much larger range, and still received the same error. I guess I am mainly looking to understand more about how to improve model fit of binomial data outside of scaling the explanatory variables being included in the model.

Example of warning:

Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio

• Rescale variables?

I am not looking for a fix to the warning message, which is why I haven't included any actual data or code, just to try and have a better understanding of ways you may improve model fit in logistic regressions.

Answer 1

A reproducible example would be more useful than you think, because it could be helpful to look more holistically at the data. However:

• when the null model is best, some combination of the usual suspects is at work:

  • your data set is too small
  • you haven't measured the important predictors
  • your data are too noisy (binary outcomes are intrinsically low-quality, it takes a larger number of binary observations to reach the same effective sample size: per Harrell Regression Modeling Strategies, the effective sample size of a binary response is min(# successes, # failures))
  • your model is misspecified in some way (you've assumed linearity when the real pattern isn't linear, you've ignored important interactions ...)

  The last case (model misspecification) is the only problem you may be able to check/fix without getting new data, e.g. use diagnostic plots (e.g. performance::check_model() or DHARMa::plot(simulateResiduals())) to see if there are patterns in the residuals that you could address.

2. Beyond scaling explanatory variables, or changing the model specification (see previous point), there's not very much you can do. Probably the best thing to do is to run allFit() (see ?lme4::allFit, ?lme4::convergence) to see that the results are stable across different optimizers.

Here's a little bit of code that might give some more insight into where the problem lies. (The glmmTMB package has a diagnose() function that tries to automate some of this: it hasn't yet been back-ported to lme4, but you could try running some of the problematic models in glmmTMB and running diagnose() on them ...)

library(lme4)
gm1 <- glmer(as.integer(Reaction) ~ Days + (Days|Subject), sleepstudy,
              family = poisson)
H <- gm1@optinfo$derivs$Hessian
nm <- c(names(getME(gm1, "theta")), names(fixef(gm1)))
dimnames(H) <- list(nm, nm)
ee <- eigen(H, symmetric=TRUE)
max(ee$values)/min(ee$values)

The last value above is (I think) the eigenvalue ratio that lme4 tests. Here's how to see what parameters in your model might be problematic (look at the largest and smallest eigenvectors to see if you get any clues):

ev <- ee$vectors
rownames(ev) <- nm
print(ev,digits = 2)
                             [,1]    [,2]   [,3]    [,4]     [,5]
Subject.(Intercept)      -0.00727 -0.0018 -0.050  0.9984  0.02613
Subject.Days.(Intercept) -0.15352  0.0145 -0.987 -0.0500 -0.00184
Subject.Days             -0.98793 -0.0219  0.153  0.0004 -0.00074
(Intercept)               0.00016  0.0339  0.001  0.0263 -0.99908
Days                     -0.01942  0.9991  0.018  0.0016  0.03395

This model isn't problematic, but if it were I would be thinking about the variance of slope across subjects (Subject.Days), which is heavily loaded on the first eigenvector (-0.99), and the intercept, which is heavily loaded on the last (smallest) eigenvector. (Subject.Days does not exactly correspond to the across-subject variance of slopes, but I'm not going to go down that rabbit hole unless necessary.)

PS while I'm not crazy about the whole "fit all models, then present an AIC-ranked table" approach, if you do want to do it (and I appreciate that it's the standard in some fields), MuMIn::dredge() could streamline your modeling workflow considerably.