# R glmer poisson model and random effect singularity

I am running a Poisson model using glmer to look at the effect of a treatment on fat scores (scale from 1-5, hence the Poisson) of an animal. There are multiple timepoints in which fat scores were taken, so since there are repeated measures, I have a random effect of Animal ID added to the model.

The model looks something like this:

glmer(Fat ~ Treatment + Sex + Timepoint + Sex:Treatment + Timepoint:Treatment +
(1 | AnimalID))

However, every single time I try to run the model, it gives me a singularity error. Essentially, it says the random effect's variance is 0. Since I have repeated measures, though, I do have to keep it in the model, correct? I tried finding multiple ways to fix the singularity (taking out individuals that have the same exact fat score throughout the experiment, removing scores that appear to be from one time point, etc.), but nothing is fixing it. The model will still run, but it says the variance is 0, which is puzzling to me. Will anything fix this, or should I just keep the model as is with the singular random effect?

Edit for clarity: I have three different timepoints: one before treatment and two after treatment. There are 2 treatment groups with samples from animals in each treatment group taken during all 3 timepoints (all individuals were sampled 3 times). I have the Timepoint:Treatment interaction in there as I hypothesize that the effect of treatment will be different at each timepoint (ie. there should not be an effect during the first timepoint, but there could be if I happened to have put leaner animals in one treatment vs another). This is the output for this model:

> summary(resultsfat)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: poisson  ( log )
Formula: as.integer(FatScore) ~ Timepoint + Sex + Treatment + Sex:Treatment +
Timepoint:Treatment + (1 | AnimalID)
Data: parentdata

AIC      BIC   logLik deviance df.resid
421.7    447.6   -201.9    403.7      122

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.3137 -0.6158 -0.2114  0.4213  2.0945

Random effects:
Groups   Name        Variance  Std.Dev.
AnimalID (Intercept) 1.738e-16 1.318e-08
Number of obs: 131, groups:  AnimalID, 41

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                 0.5495     0.1831   3.002  0.00268 **
TimepointBtrt1              0.1642     0.2156   0.761  0.44639
TimepointCtrt2              0.5557     0.1993   2.787  0.00531 **
SexM                        0.1287     0.1596   0.807  0.41990
TreatmentS                  0.1545     0.2495   0.619  0.53584
SexM:TreatmentS            -0.1987     0.2370  -0.838  0.40183
TimepointBtrt1:TreatmentS  -0.3985     0.3133  -1.272  0.20334
TimepointCtrt2:TreatmentS  -0.3202     0.2849  -1.124  0.26106
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) TmpnB1 TmpnC2 SexM   TrtmnS SxM:TS TB1:TS
TimpntBtrt1 -0.658
TimpntCtrt2 -0.712  0.597
SexM        -0.485  0.017  0.019
TreatmentS  -0.734  0.483  0.522  0.356
SxM:TrtmntS  0.326 -0.012 -0.013 -0.673 -0.479
TmpntBt1:TS  0.453 -0.688 -0.411 -0.012 -0.613 -0.001
TmpntCt2:TS  0.498 -0.418 -0.700 -0.013 -0.678  0.010  0.536
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
$$`$$
• Have you considered a Bayesian model with a prior on what one believes (possibly a bit more uncertain than that aka "a weakly informative prior") about the size of the between animal SD? If you did that, you would not have to worry about the particular set of data being such as to the maximum likelihood estimate of that SD is 0. Also, this is not really an ideal setting for a Poisson model. A scale from 1-5 misfits that distribution in several ways: 1) limited scale (Poisson is on $[0, \infty)$ and 2) 0 not included. If this is just 1,2,3,4,5, then, why not ordinal logistic regression? Commented Oct 11, 2023 at 7:45
• How many degrees of freedom does your model have? If it is low then the random effects become difficult to compute. Is the Poisson distribution a good model? If it is bad then you get large error terms and the random effects are negligible. Also, the Poisson model will assume a certain variance that may not be realistic (like in this question stats.stackexchange.com/questions/628067). Commented Oct 11, 2023 at 9:26
• Thank you all for your input! As a biologist, I do not understand the detailed math behind the models, but it does seem that Poisson distribution may not be the best way to run this model. I have never used Bayesian models, but from my understanding, they are predictor models and not necessarily for analyzing differences between treatment groups (but again, as someone who is not a math person, I could be wrong). Commented Oct 12, 2023 at 2:58
• To clarify on the timepoints, they were taken based on breeding stage (ie. not a fixed time for every animal). The first timepoint is prior to treatment, and the rest are after treatment (thus I must keep those in the model for consideration). I do expect the effect of treatment at different timepoints to be different. Commented Oct 12, 2023 at 3:00
• I apologize for the confusion; I have updated the question for clarity and included the output. Timepoint is a factor, and the Timepoint:Treatment interaction is ensuring that the effect of treatment is different across timepoints. I do need to ensure that I did not happen to pick leaner animals in one treatment or another, so I believe that I need to keep it that way to test my hypotheses. Thank you all again! Commented Oct 12, 2023 at 14:03

By default, glmer() uses the Laplace approximation to calculate the model's likelihood function, which in some cases may result in such numerical issues. It would be best to switch to using the adaptive Gaussian quadrature in which you can tune the approximation. For more information, check this post.

You can invoke the adaptive Gaussian quadrature in glmer() using the nAGQ argument and setting it to, e.g., 11 or 15. Another package with this capability is GLMMadaptive.

• Thank you for your insight! I am not entirely sure that I understand tuning the approximation, but I will look more into it. Thanks again! Commented Oct 12, 2023 at 3:03

One cause for the singular fit may be the decision to model an ordinal variable with Poisson regression.

As @Björn and @SextusEmpiricus point out in the comments, the outcome — fat scores on a scale from 1 to 5 — is an ordinal variable: the scores are ordered categories, from least fat to most fat. It's appropriate to model these scores with ordinal logistic regression. (We use Poisson regression to model counts — number of occurrences — of an event in a fixed time interval.)

Moreover, there are two post-treatment measurements (+ a baseline measurement taken before the treatment is applied) from each animal, so the observations are clustered. This suggests to use a mixed effects model for ordinal outcomes.

Here are two options to fit such a model.

• Continuation ratio model with the GLMMadaptive package by @DimitrisRizopoulos. The code to specify the model would be something like this:
fit <- mixed_model(
Fat ~ Sex * Treatment + Timepoint * Treatment,
random = ~ 1 | AnimalID,
family = binomial(),
data = your_data
)

• Bayesian proportional odds model with the rmsb package by @FrankHarrell.
fit <- blrm(
Fat ~ Sex * Treatment + Timepoint * Treatment + cluster(AnimalID),
data = your_data
)

See also this section of Regression Modeling Strategies course: Bayesian Proportional Odds Random Effects Model.

Finally, I suggest to consider treating the baseline score as a covariate (on the right side of the model formula) rather than an outcome (on the left side). That is, write the model as:

Fat ~ factor(BaselineFat) + Sex * Treatment + Timepoint * Treatment

after appropriately re-arranging the data. The motivation to do so is that the baseline, which is (likely to be) predictive of the outcome after treatment, is not itself affected by the treatment in a randomized experiment.