Multiple vs Single Predictor Variables for GLMM Pairwise Comparisons

Question

I am running binomial GLMMs in R to determine whether species presence (binary) on a hydrophone is different between seasons (i.e. spring, summer, fall, and winter) and photoperiods (i.e. day, night, dawn, and dusk). My models include a temporal autocorrelation structure with Group given as a single value since I am using a single hydrophone. I understand that I can use the multcomp or emmeans packages to conduct pairwise comparisons on my model(s), but am unsure whether I should run separate GLMMs with a single predictor:

M1 <- glmmTMB(Presence ~ Photoperiod + ou(Time - 1|Group), data = df, family = binomial(link="logit"))
M2 <- glmmTMB(Presence ~ Season + ou(Time - 1|Group), data = df, family = binomial(link="logit"))

or a single model with both predictors:

M3 <- glmmTMB(Presence ~ Photoperiod + Season + ou(Time - 1|Group), data = df, family = binomial(link="logit"))

Using car::Anova(), both photoperiod and season have significant effects on presence in all three models. However, my pairs(emmeans()) results are different enough to effect significance depending on whether I model the predictors together (M3) or separately (M1 and M2). M3 has a slightly lower AIC value than M1 or M2.

Are there any ways to justify using two models with a single predictor vs a single model with multiple predictors if my goal is determining whether species presence differs between photoperiods and seasons? I also have a few environmental covariates (i.e. sea-surface temperature, chlorophyll concentration, and sea level) I was planning on putting in a separate model, but am now wondering if I should model them alongside season and photoperiod?

As you may be able to tell from this post, I am a bit of a modelling novice and so am partial to simpler methods so long as they don't lead to incorrect/misleading results.

Answer 1

I suggest that you use one model with all the predictors. Here are two good reasons:

1. If you have a model with one predictor, you only explain the effects of that one predictor, and any unexplained effects, including the effects due to other predictors, go into the residual variation. If all the predictors that have effects are in one model, then the residual errors will be much smaller, the SEs of your estimates will be smaller, and your statistical tests will be more powerful.

2. Some of those predictors may interact with one another. If so, your one-predictor models will produce misleading estimates. You should try models with interactions and simplify from there if some interactions have high p values in an anova test (say, using car::Anova()).