I am trying to make a model for the different amount of species caught in different traps in 3 different locations on 3 height levels, along with 3 transects per location (resulting in 9 traps per location). I collected the contents of the traps every 2 weeks from April to July. This is the model I'm trying to run:

glmm1<- glmer(speciesrichness ~ Height*Meander*Date*Transect + 
               family = poisson, data = SpeciesRichness)

I have been told that this model may be too complex, and that I should try model selection to solve this. Based on the AIC (430.7), I think that this is indeed the case, although I have no experience with this. Here is a preview of my data:

> head(SpeciesRichness)
# A tibble: 6 x 6
# Groups:   Sample_Code, Date, Height, Meander [6]
  Sample_Code Date       Height Meander Transect speciesrichness
  <chr>       <chr>      <chr>  <chr>      <dbl>           <int>
1 M1.1a       19/06/2023 A      M1             1              22
2 M1.1b       3/07/2023  B      M1             1               4
3 M1.1b       5/06/2023  B      M1             1              15
4 M1.1b       8/05/2023  B      M1             1               8
5 M1.1c       19/06/2023 C      M1             1               7
6 M1.1c       22/05/2023 C      M1             1              10

I have tried to make a regression tree using the tree function, but it seems this won't work because Height and Meander are non-numeric. I also tried stepwise model selection based on AIC but got the following error:

stepAIC(glmm1,direction = "both")
Error in `$<-`(`*tmp*`, formula, value = Terms) : 
  no method for assigning subsets of this S4 class

So my question is, how can I deal with these problems in order to correctly simplify my model?

  • 2
    $\begingroup$ "Based on the AIC (430.7), I think that this is indeed the case, although I have no experience with this." AIC is a relative measure for comparing models. A single AIC value has no meaning - the only way AIC provides information is if you fit 2 (or more) models and look at the difference in AIC between them. $\endgroup$ Commented Apr 4 at 14:15
  • $\begingroup$ Rather than trying to automate model selection, I would recommend (a) starting with a simple model without interaction terms and adding complexity, rather than starting with a highly complex model and subtracting complexity; (b) examining the model coefficients to see what is useful in the model and what isn't; and (c) making sure you have a theoretical basis for everything in your model.... $\endgroup$ Commented Apr 4 at 14:22
  • 2
    $\begingroup$ Date is immediately suspect as a model term - the snippet of data you show has no repeated dates. Do you really want Date as one of the main effects in your model? Will your audience be satisfied with a conclusion like "the most important variable impacting the number of species caught is if the date is 19/06/2023 or not. Species richness increased by 10 on June 19, 2023, all other factors being equal." $\endgroup$ Commented Apr 4 at 14:25
  • $\begingroup$ On the programming side, I'll comment on the stepAIC() error - that function is made to work with non-mixed models fit with lm() or glm(). Since you're using a mixed model with glmer, it won't work. You could potentially use stepAIC() or another automated model selection technique like the LASSO without random effects to select your main effects, and then re-fit the model with glmer(). But I can't really comment on the validity of such an approach. $\endgroup$ Commented Apr 4 at 14:33
  • $\begingroup$ As noted already you may want to consider how date appears in your model. 1) if multiple samples were collected at different sites on the same date, or the same site revisted on different dates, then you probably want date as a random effect 2) if you are interested in change over time due to eg seasonality you may want to convert date to a numeric value like Julian date. As is date appearing as a fixed and random effect is odd $\endgroup$
    – N Brouwer
    Commented Apr 5 at 4:22

2 Answers 2


One comment which stuck out to me:

I have been told that this model may be too complex, and that I should try model selection to solve this. Based on the AIC (430.7), I think that this is indeed the case, although I have no experience with this.

AIC doesn't really answer this question by itself. As Gregor noted in the comments, this is typically used to compare AIC values for different models rather than models by themselves. I also agree with Gregor that the date seems to be a completely arbitrary variable. I'm guessing what you actually want is some kind of ordered time-based variable that isn't simply set of characters (such as Year 1, Year 2, Year 3, etc.). I would very carefully consider why you need this interaction in the first place (the beginning of this answer explores why an interaction may exist in a hypothetical scenario with frogs).

I will say that I have a bias against high-level interactions like yours. Once we move beyond a two-level interaction, the results become difficult to explain, the model may not have the proper statistical power necessary to observe a relevant effect, and they are often not determined by previous research.

What's more, your random effects term seem exceedingly complex, and I'm sure that the coding is wrong (height cannot be a random effect, only factor variables can be random intercepts). I would refer to the original lme4 article to get a sense of what this is actually doing, along with this article which notes the issues with random effect complexity.

Finally, you should never use stepwise regression, even for mixed models. Beyond the statistical issues it poses, it replaces your thinking about model selection with a sub-optimal algorithm.


As noted by Gregor above, there are a few different ways to compare models, but all will involve running additional models and comparing the AICs (or BICs) to see if the random effects structure is necessary to model the data. Barr et al. (link here: https://www.sciencedirect.com/science/article/pii/S0749596X12001180) recommend starting with the "maximal model" (i.e., the model with the most possible random effects terms specified) and always keeping it barring convergence errors.

Another way is to start with the maximal model and compare less complex models step-by-step (i.e., model 1 is the maximal, model two removes one random slope, model 3 removes one different random slope, etc. all the way down to a random-intercept only model, or the least complex that makes sense for your data), then using the anova() function to compare them via a chi-square test.

This is similar to the process specified by Bates et al. (2015) (link here: https://www.researchgate.net/publication/278734089_Parsimonious_Mixed_Models).

There isn't necessarily one correct way to specify mixed effects models, so I would recommend reading further and taking the approach which makes the most sense for your data- the maximal model may not be ideal if you have many possible random effects, for example.

  • 1
    $\begingroup$ I completely disagree about fitting the maximal model first. It should go the opposite way...fitting basic random effects first and then, if necessary, work one's way up to more model complexity. The second article you link here argues for the opposite of what you prescribe if my memory is correct. $\endgroup$ Commented Apr 5 at 18:55
  • $\begingroup$ Yes, I think you're right about the second article, thanks! And I was taught to start with the maximal model and work your way down- like I said, there probably isn't one right way to do it. $\endgroup$
    – Paul
    Commented Apr 6 at 20:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.