# Interpreting a generalised linear mixed model with binomial data

I have a generalised linear mixed model with binomial response data, the model:

model <- glmer(RespYN ~ Treatment + Gender + Length + (1 | Anim_ID),
data = animDat,


I am no statistician (I'm a biologist) so I have no idea how to interpret the data. With a linear mixed model I understand, due to the mean differences, etc. being means of the response variables. With binomial GLMMs I am unsure.

How do I prove that the treatment is causing/not causing the response?

Here is my output (sorry for dumping all of it):

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial  ( logit )
Formula: RespYN ~ Treatment + Gender + Length + (1 | Anim_ID)
Data: animDat

AIC      BIC   logLik deviance df.resid
142.1    158.1    -66.1    132.1      176

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.6913 -0.3021 -0.2059  0.4435  3.8066

Random effects:
Groups       Name        Variance Std.Dev.
Cockroach_ID (Intercept) 0        0
Number of obs: 181, groups:  Cockroach_ID, 10

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  5.03046    3.62723   1.387    0.165
TreatmentPo -4.06399    0.48900  -8.311   <2e-16 ***
GenderM      0.13323    0.49365   0.270    0.787
Length      -0.05896    0.05758  -1.024    0.306
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) TrtmnP GendrM
TreatmentPo -0.175
GenderM     -0.412 -0.039
Length      -0.995  0.139  0.348
convergence code: 0

• Could you explain what Anim_ID and Cockroach_ID represent in your experiment and how they are related? Did you happen to get any warning messages like "boundary (singular) fit: see ?isSingular"? – whuber Jan 14 at 21:20

The interpretation is the same as for a generalised linear model, except that the estimates of the fixed effects are conditional on the random effects.

Since this is a generalized linear mixed model, the coefficient estimates are not interpreted in the same way as for a linear model. In this case you have a binary outcome with a logit link, so the raw estimates are on the log-odds scale. The estimated coefficient for the intercept, 5.03046, is the log odds of RespYN being 1 (or whatever non-reference value it is coded as) when Lengthis equal to zero, and Treatment and Gender take their reference value. A value of zero for dLength might not make sense in your sample, since presumably it will never be negative and is always far above zero, and if so, you might want to consider centering it so that a zero value for the centered variable is more meaningful.

The estimate for Length of -0.05896 means that a 1 unit increase in Length is associated with a 0.05896 decrease in the log-odds of RespYN being 1, compared to RespYN being 0. If we exponentiate this number then we obtain the odds ratio of 0.9427445, which means that for a 1 unit increase in Length we expect to see (approximately) a 6% decrease in the odds of RespYN being 1.

The estimate for TreatmentPo of -4.06399 means that Treatment = Po is associated with 4.06399 lower log-odds than the other treatment group of RespYN being 1, compared to RespYN being 0. This can be exponentiated as above to obtain an odds ratio. The same analysis applies to Gender.

How do I prove that the treatment is causing/not causing the response?

Nothing can be proven with statistics, especially with observational studies. You can say that, while controlling for Gender,Length and the repeated measures within Anim_ID, you have found evidence that the association of Treatment with the the outcome is not zero. You could also say that, if the association of Treatment with the the outcome is actually zero, then the probability of observing the data that you have, or data more extreme, is less than 0.0000000000000002

Lastly, I notice that you have specified random intercepts for Anim_ID in your model formula, yet the model output says that Cockroach_ID is the grouping variable. This is rather odd, normally they would be the same. Moreover, the convergence code is zero which indicates that the model has not converged, and the estimated variance for the random effect is zero. This means that potentially there is no variation within Anim_ID. It would be a good idea to fit a model with glm() (ie without random intercepts but with Anim_ID as a fixed effect) and see how the model estimates compare.

• I think OP understands the "mixed" part of a GLMM. It's the "generalized" part that is the problem – Hong Ooi Jan 15 at 12:25
• @HongOoi yes, that's what the 2nd, 3rd and 4th paragraph are about......perhaps I should make that more clear. Will edit shortly. – Robert Long Jan 15 at 12:35
• @RobertLong for gender does this mean that if there's a 1.1 % chance increase of male giving respYN = 1? Does this mean that females are 98.9 %? Does the logic for treatment work in the same way? From what I get from what you're saying it is a ratio of response 1 vs 0 - not between male/female, etc.? Basically how do I get the results of the opposite treatment and gender? Also, how does this relate to the Pr(>|z|) value given? Thank you! – s33ds yesterday
• So can I say that treatment has an effect on output, with treatmentPo decreasing the likelihood of a positive response by 1-exp(-4.06399) = 98%? Does this just mean that the other treatment is significant? – s33ds yesterday
• @s33ds you can say that treatment = Po is associated with a statistically significant lower odds than the other treatment group, of RespYN being 1. It doesn't mean that the other treatment is significant - it tells you there is a statistically significant difference – Robert Long 23 hours ago

A couple of extra notes on top of what @RobertLong already answered:

• As Robert also noted, the interpretation of the coefficients from generalized linear mixed models are conditional on the random effects. Most often this is not the interpretation you are looking for. For more info on this check here.
• You have fitted the model with the default Laplace approximation. This may be inaccurate, especially for dichotomous data. You should better fit the model with the adaptive Gaussian quadrature by specifying a value much greater than one for the nAGQ argument of glmer().