Several single-species GLMMs -> A multi-species GLMM (random effects syntax?) I'm looking to describe annual relative abundance for several species over several years by analyzing point count data. I'm using a Poisson GLMM in glmmTMB for each species, it looks like this:
m[[i]] <- glmmTMB(count ~ year_numeric + (1|year_factor) + (1|site), 
                  family = "poisson", data = data[[i]])

As an aside: This has a linear trend for year, plus year random effects, which is conceptually similar to the Sauer & Link models from the Breeding Bird Survey (though they use a Poisson log-linear regression, and it's Bayesian).
Instead of having a separate GLMM for each species, though, I'd like to have one single GLMM for all species (a multi-species model). I know I could do part of this by using year_numeric*species, but: what would be the syntax for the random effects?
Thanks!
 A: As Diogo suggested, you can make use of Joint Species Distribution Models (JSDM). JSDMs are multivariate mixed-effects models that account for correlation between the responses (including, but not limited to, species). To do so they include a random effect per site and species, that is structured by a covariance matrix of correlations between the responses (sometimes referred to or interpreted as a "matrix of species associations").
However, JSDMs suffer from the fact that the the number of parameters in the covariance grows quadratically in the number of responses, so that they are infeasible to fit using standard  mixed-effects formulation for a large number of responses.
A good review on this matter is Warton et al. (2015), but a more comprehensive resource is Ovaskainen and Abrego (2020).
Using a reduced-rank structure for the covariance matrix the number of parameters can be kept in check, so that this type of model can be fitted with the glmmTMB R-package using the reduced-rank covariance structure, or with R-packages such as HMSC, Boral, or gllvm. Alternatively, sjSDM takes an elastic net approach to fitting JSDMs. The ecoCopula R-package uses a Gaussian copula approach to fitting multispecies models, which mostly makes it incredibly fast, but also allows it to be used to determine species associations, and for unconstrained ordination.
The first four approaches are Generalized Linear Latent Variable Models (GLLVMs) that can also be used for ordination purposes (the first three for model-based unconstrained ordination, while the gllvm R-package also supports ordination with predictors and additional random-effects), whereas sjSDM is less suited for that type of application. In GLLVMs the rank of the covariance matrix needs to be a-priori set (i.e. the number of latent variables), in sjSDM that is not needed. However, HMSC including a nifty trick known as "the infinite factor model", which can automatically select the rank for you.
Each software implementation and approach thus has it pros and cons, but arguably glmmTMB is easiest to implement. An example for a multispecies model with the reduced-rank structure in glmmTMB is available in their vignette.
Disclaimer: I am one of the co-authors/developers of the gllvm R-package.
A: First, you create an additional column containing a factor, species, which contains as levels the different species. Next, you have to decide which effects you want to be common for all species and which should be species-specific. E.g., if you want the linear effect of year_numeric to be species-specific, in your formula you have to replace year_numeric with year_numeric*species, as you have already suggested:
glmmTMB(count ~ year_numeric*species + (1|year_factor) + (1|site), 
              family = "poisson", data = data)

Note, that with this formula you will get e.g. the same random site effect for all species, while, when fitting a new model for each species i as before, you would usually get for each species a different random site effect.
Finally, if you want to use a random effect for year_numeric for each species, you replace this term with (0 + year_numeric|species):
glmmTMB(count ~ (0 + year_numeric|species) + (1|year_factor) + (1|site), 
              family = "poisson", data = data)

The added zero 0 + in (0 + year_numeric|species) indicates that no species-dependent random effects offset is to be computed. If you do want that, leave out the 0 +  and just use (year_numeric|species).

To be more precise: the term (year_numeric|species) implicitly means (1 + year_numeric|species), i.e each species would have its own extra offset, a random effect offset, and those offsets would be coupled, both with each other and with the random effect slope of year_numeric. If you don't want this implicit 1 +, you have to explicitly say 0 +.
In the equation:
y = a*year_numeric + b

the parameter a is a slope and the parameter b is an offset. And in this situation here we have a different offset for each species, e.g. we have one offset $b_{hippo}$ for hippos and one offset $b_{lion}$ for lions.
Finally, to get a random offset for each combination of site and species (both I presume to be factors), I would suggest the term
(1 | species:site )

However, keep in mind that the main danger with those mixed effect models is that one is easily tempted to overdo it, creating unnecessarily complex models that are not easily interpretable anymore.
A: You can also consider running a Joint Species Distribution Model (JSDM), which is basically a multivariate GLMM. You can include temporal autocorrelation structures in the model implemented in Hierarchical Model of Speceis Communities (HMSC) and Generalized Linear Latent Variable Models (GLLVM)
