I found this paper very useful in general:

https://peerj.com/preprints/27320.pdf

Generalized linear models (GLMs) provide a powerful tool for analyzing count data, usually with Poisson-distributed errors (mean = variance), but not all count data meet the assumptions of the Poisson distribution. Thus, we need to test if the variance > mean (over dispersion) or if the number of zeros is greater than expected, in which case you can choose a different model, such as negative binomial :)

I believe Andrew Paul McKenzie Pegman is on the correct track. Isn't this Panel data as described here: https://www.bing.com/search?q=panel+data+modeling&qs=AS&pq=panel+data+mo&sk=AS1&sc=6-13&cvid=06AB4A74A2914F56A73DF01BD803B836&FORM=QBLH&sp=2&ghc=1

and hence you need to model the data as described here, where Pegman's approach is appropriate to one of the data components, see, for example, https://en.wikipedia.org/wiki/Panel_analysis. For some further approaches, you might consult these notes:

Davidian Notes

https://www4.stat.ncsu.edu/~davidian/st732/notes.html

Multivar. Long. Models:

https://www.researchgate.net/post/What_statistical_test_is_appropriate_for_pre_and_post_trails_with_1_before_treatment_value_and_4_after_treatment_values_assessed_weekly

Hope these are of some value, David Booth

I second Stefano Larsen s answer. The paper gives a very nice introduction to GAMs and how to fit both long-term and seasonal trends as well as random effects (e.g. sites) in your model. They even have a fish data example :)

You could also use mgcv::gam() for fitting random effects, e.g.:

mgcv::gam(y ~ s(var1) + s(site, bs = 're'))

Here's the published version:

https://peerj.com/articles/6876.pdf

Also have a look here:

https://www.fromthebottomoftheheap.net/2014/05/09/modelling-seasonal-data-with-gam/

You haven't mentioned glmmTMB in your pro/con description, though you mention it in the title, so I'll add some pros for this package. The advantage of glmmTMB is that you can easily model complex / nested / cross classified random effects structures and you have different correlation options (like AR1 etc.).

Furthermore, you have various diagnostic options for glmmTMB as well, to quote from another thread:

You could use the "performance" package to calculate indices like r2() or icc() [1]. You can also check your model for overdispersion or zero-inflation with the "performance" package (check_overdispersion() or check_zeroinflation()). Furthermore, it can create diagnostic plots (which requires the "see" package for plotting capabilities [2]). Another, imho important, package in this context is "DHARMa" for residual diagnostic plots [3].

Finally, a comment on your model design: when analysing longitudinal data / repeated measurements, I'd suggest adding the time variable both as fixed effect and random slope. And I would treat "year" as continuous, maybe using a quadratic / cubic or spline to model time variation, so you would have something like:

model

Thank you so much Daniel Lüdecke for the informative answer. Yes I see all the pros of glmmTMB. And thank you for commenting on the design. The reason why you would include year also in the random structure (random slope) is because you would expect different trends in different years?

"The reason why you would include year also in the random structure (random slope) is because you would expect different trends in different years"

Yes, for longitudinal / panel data or repeated measurements, it is plausible to assume that the time trend varies by group / subject. This usually gives you more accurate estimates, and this type of design also takes the uncertainty due to varying time trends into account.

There is an older, but still informative and well written overview on different designs and how to write this in lme4 (or glmmTMB) notation: https://rpsychologist.com/r-guide-longitudinal-lme-lmer

Dear Daniel Lüdecke that is a very usefull link. Thank you so much for your answers.

Best regards