0

I am relatively new to both R and Stack overflow so please bear with me. I am currently using GLMs to model ecological count data under a negative binomial distribution in brms. Here is my general model structure, which I have chosen based on fit, convergence, low LOOIC when compared to other models, etc:

enter image description here

My goal is to characterize population trends of study organisms over the study period. I have created marginal effects plots by using the model to predict on a new dataset where all covariates are constant except year (shaded areas are 80% and 95% credible intervals for posterior predicted means):

enter image description here

I am now hoping to extract trend magnitudes that I can report and compare across species (i.e. say a certain species declined or increased by x% (+/- y%) per year). Because I use poly() in the model, my understanding is that R uses orthogonal polynomials, and the resulting polynomial coefficients are not easily interpretable. I have tried generating raw polynomials (setting raw=TRUE in poly()), which I thought would produce the same fit and have directly interpretable coefficients. However, the resulting models don't really run (after 5 hours neither chain gets through even a single iteration, whereas the same model with raw=FALSE only takes a few minutes to run). Very simplified versions of the model (e.g. count ~ poly(year, 2, raw=TRUE)) do run, but take several orders of magnitude longer than setting raw=FALSE, and the resulting model also predicts different counts than the model with orthogonal polynomials. My questions are (1) what is going on here? and (2) more broadly, how can I feasibly extract the linear term of the quartic polynomial describing response to year, or otherwise get at a value corresponding to population trend?

I feel like this should be relatively simple and I apologize if I'm overlooking something obvious. Please let me know if there is further code that I should share for more clarity–I didn't want to make the initial post crazy long, but happy to show specific predictions from different models or anything else. Thank you for any help.

reco
  • 11
  • 1
  • 1
    One thought is that the 4th order terms may create variables of drastically different magnitude. This can create difficulties for convergence. You may want to center and possibly even scale your numeric date variable. Another thought is that having polynomial terms for both year and month and Time.dec may be overkill. A a 3rd order polynomial for one numeric date variable might be enough to describe a variety of shapes over time. – Arthur Apr 14 '22 at 18:03
  • Another problem that would be created by using `raw=TRUE` is that the resulting polynomials are not orthogonal. In fact, they will be very highly correlated, which may pose problems for the sampler. – DaveArmstrong Apr 14 '22 at 20:15

0 Answers0