Correct binomial GLMM for temporal trends in species occurrences

Question

I have a dataset with 80 species which were sampled in about 120 water bodies at two time periods (historical / recent). Only presence/absence of the species in each water body is considered. The data looks like this in long format:

water_no_name              species         success  time    
Water body 1               Species A             0  t1      
Water body 2               Species A             0  t1      
Water body 3               Species A             1  t1      
Water body 4               Species A             0  t1      
Water body 5               Species A             0  t1      
Water body 1               Species A             1  t2      
Water body 2               Species A             0  t2      
Water body 3               Species A             0  t2      
Water body 4               Species A             0  t2      
Water body 5               Species A             0  t2      
Water body 1               Species B             1  t1      
Water body 2               Species B             0  t1      
Water body 3               Species B             0  t1      
Water body 4               Species B             1  t1      
Water body 5               Species B             1  t1      
Water body 1               Species B             1  t2      
Water body 2               Species B             0  t2      
Water body 3               Species B             1  t2      
Water body 4               Species B             1  t2      
Water body 5               Species B             0  t2   

I would like to check which species have significantly increased or decreased in frequency over time.

Therefore I decided to calculate binomial GLMMs using the water body (= site) as random factor to account for the repeated measures design. I am unsure to decide which of these two approaches would be more appropriate:

1) Calculate separate GLMMs for each species:

mod <- glmer(success ~ time + (1 | water_no_name), family = binomial)

With this approach, I would then make a p-value correction (FDR) for the number of species.

2) Calculate combined GLMM for all species using interactions between time and species:

mod <- glmer(success ~ time * species + (1 | water_no_name), family = binomial)

To analyse the interactions, with this approach I would run a post-hoc analysis to see for which species the contrast t1/t2 is significant:

emmeans(mod, pairwise ~ time | species)

When calculating both models, the results are largely similar, with some differences for species with lower counts.

I understand, that one major difference is, that the complex model (2) uses the same random site effect for all species, whereas using separate models (1) I get a different random site effect for each species (see Several single-species GLMMs -> A multi-species GLMM (random effects syntax?))

Still I am not sure which of the approaches would be more appropriate in my case. I find the simpler approach (1) more intuitive to interpret and also sufficient to answer my question. I am also not sure if it would be better to have the random site effect for all species.

Can anyone give me advice on which approach I should prefer here?

Or do you recommend another approach?

Answer 1

It looks like you have a multivariate dataset, with a species composition matrix of 80 species by 120 sites. The question seems to be if there's difference in species composition, and which individual species contributed the most, between two time periods. I highly recommend you to switch to a model-based approach, perhaps using mvabund and modelling the species composition matrix a function of time, specifying binomial distribution.