# Instability in glmer model at zero random effect variance

I am getting very different results from glmer (using version 0.999375.37 of lme4) when running in different environments (I can't pinpoint the source of discrepancy, but let's assume some floating point precision difference in the input data). Here is the invocation:

glmer(data=data, n ~ 1 + offset(log(mu)) + (1 | id),


Sometimes the output is this:

  0:     532.77427: 0.610172 -0.0811043
1:     525.91180: 0.531918 -0.102716
2:     511.77820: 0.373325 -0.137516
3:     495.51600:  0.00000 -0.127691
4:     478.25662:  0.00000 -0.100781
5:     474.47031:  0.00000 -0.0806405
6:     474.46819:  0.00000 -0.0811088
7:     474.46819:  0.00000 -0.0811044
8:     474.46819:  0.00000 -0.0811043
9:     474.46819:  0.00000 -0.0811043
Generalized linear mixed model fit by the Laplace approximation
Formula: n ~ 1 + offset(log(mu)) + (1 | id)
Data: data
AIC           BIC        logLik    deviance
478.4681939 488.556259695 -237.23409695 474.4681939
Random effects:
Groups Name        Variance Std.Dev.
id     (Intercept)  0        0
Number of obs: 1146, groups: id, 160

Fixed effects:
Estimate       Std. Error  z value   Pr(>|z|)
(Intercept) -0.0811043131501  0.0100762538900 -8.04905 8.3436e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


and sometimes it's this:

  0:     532.77427: 0.610172 -0.0811043
1:     525.91180: 0.531918 -0.102716
2:     511.77820: 0.373325 -0.137516
3:     495.51540:  0.00000 -0.127690
4:     478.25643:  0.00000 -0.100780
5:     474.47031:  0.00000 -0.0806405
6:     474.46819:  0.00000 -0.0811088
7:     474.46819:  0.00000 -0.0811044
8:     474.46819:  0.00000 -0.0811043
9:     474.46819:  0.00000 -0.0811043
10:     474.46819:  0.00000 -0.0811043
11:     474.46819: 1.29719e-05 0.568362
Generalized linear mixed model fit by the Laplace approximation
Formula: n ~ 1 + offset(log(mu)) + (1 | id)
Data: data
AIC  BIC logLik deviance
5700 5710  -2848     5696
Random effects:
Groups      Name        Variance   Std.Dev.
id          (Intercept) 1.6827e-10 1.2972e-05
Number of obs: 1146, groups: id, 160

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.568362   0.007282   78.05   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Note that fitting the model without the random effect gives:

Call:  glm(formula = n ~ 1 + offset(log(mu)), family = poisson(link = log),
data = data, weights = w)

Coefficients:
(Intercept)
-0.0811043131378


First question: Why is the result unstable? IIUC the model is of the form

n ~ Pois(mu * exp(a + b_id))

where the b_id are drawn from a normal distribution N(0, theta). In the verbose output the third column is a, and the second column is related to theta. As theta approaches zero, the b_id's must go to zero, so I would expect a to converge to the maximum likelihood estimate for the model without random effects. This is exactly what we see, except for iteration #11 in the second outcome.

Second question: How can I automatically detect and avoid this? Is it sufficient to fall back to the simple model (without random effect) when theta is below some threshold?

Note that I'm interested in the estimated values for each level, not the variability (the opposite of the situation discussed in this FAQ). I also want some penalty on differences between the values for different id's, so I don't want to model id as a fixed effect.

• Sorry I missed this the first time around. Do you still see this with recent versions of lme4?? – Ben Bolker Oct 19 '14 at 1:10