# Binomial glmm with a categorical variable with full successes

I am running a glmm with a binomial response variable and a categorical predictor. The random effect is given by the nested design used for the data collection. The data looks like this:

m.gen1$treatment [1] sucrose control protein control no_injection ..... Levels: no_injection control sucrose protein m.gen1$emergence
[1]  1  0  0  1  0  1  1  1  1  1  1  0  0....
> m.gen1$nest [1] 1 1 1 2 2 3 3 3 3 4 4 4 ..... Levels: 1 2 3 4 5 6 8 10 11 13 15 16 17 18 20 22 24  The first model I run looks like this m.glmm.em.<-glmer(emergence~treatment + (1|nest),family=binomial,data=m.gen1) I get two warnings that look like this: Warning messages: 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0240654 (tol = 0.001, component 4) 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control\$checkConv,  :
Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?


The model summary shows that one of the treatments has a unusually large standard error, which you can see here:

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)         2.565      1.038   2.472   0.0134 *
treatmentcontrol   -1.718      1.246  -1.378   0.1681
treatmentsucrose   16.863   2048.000   0.008   0.9934
treatmentprotein   -1.718      1.246  -1.378   0.1681


I tried the different optimizers from glmer control and functions from other packages, and I get a similar output. I have run the model using glm ignoring the random effect, and the problem persist. While exploring the data I realized that the treatment with a high Std. error has only successes in the response variable. Just to check whether that could be causing the problem I added a fake data point with a "failure" for that treatment and the model runs smoothly, and gives reasonable standard error. You can see that here:

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)        3.4090     1.6712   2.040   0.0414 *
treatmentcontrol  -1.8405     1.4290  -1.288   0.1978
treatmentsucrose  -0.2582     1.6263  -0.159   0.8738
treatmentprotein  -2.6530     1.5904  -1.668   0.0953 .


I was wondering if my intuition is right about the lack of failures for that treatment preventing a good estimation, and how can I work around this issue.

Your intuition is exactly right. This phenomenon is called complete separation. You can find quite a lot (now that you know its name) Googling around ... It is fairly thoroughly discussed here in a general context, and here in the context of GLMMs. The standard solution to this problem is to add a small term that pushes the parameters back toward zero -- in frequentist contexts this is called a penalized or bias-corrected method. The standard algorithm is due to Firth (1993, "Bias reduction of maximum likelihood estimates" Biometrika 80, 27-38), and is implemented in the logistf package on CRAN. In Bayesian contexts this is framed as adding a weak prior to the fixed-effect parameters.

To my knowledge Firth's algorithm hasn't been extended to GLMMs, but you can use the Bayesian trick by using the blme package, which puts a thin Bayesian layer over the top of the lme4 package. Here's an example from the above-linked GLMM discussion:

cmod_blme_L2 <- bglmer(predation~ttt+(1|block),data=newdat,
family=binomial,
fixef.prior = normal(cov = diag(9,4)))


The first two lines in this example are exactly the same as we would use in the standard glmer model; the last specifies that the prior for the fixed effects is a multivariate normal distribution with a diagonal variance-covariance matrix. The matrix is 4x4 (because we have 4 fixed-effect parameters in this example), and the prior variance of each parameter is 9 (corresponding to a standard deviation of 3, which is pretty weak -- that means +/- 2SD is (-6,6), which is a very large range on the logit scale).

The very large standard errors of the parameters in your example are an example of a phenomenon closely related to complete separation (it occurs whenever we get extreme parameter values in a logistic model) called the Hauck-Donner effect.

Two more potentially useful references (I haven't dug into them yet myself):

• Gelman A, Jakulin A, Pittau MG and Su TS (2008) A weakly informative default prior distribution for logistic and other regression models. Annals of Applied Statistics, 2, 1360–383.
• José Cortiñas Abrahantes and Marc Aerts (2012) A solution to separation for clustered binary data Statistical Modelling 12(1):3–27 doi: 10.1177/1471082X1001200102

A more recent Google scholar search for "bglmer 'complete separation'" finds:

• Quiñones, A. E., and W. T. Wcislo. “Cryptic Extended Brood Care in the Facultatively Eusocial Sweat Bee Megalopta genalis.” Insectes Sociaux 62.3 (2015): 307–313.
• wow thanks very much!! This makes perfect sense, and the model now runs smoothly with the bglmer. I would just have one more question, can I use the methods as in lme4 to asses the random and fixed effects, in other words to compare different models?
– AtiQP
Jan 8, 2015 at 10:09
• I would say so, but I don't know if there's any formal and/or peer-reviewed support for my opinion ... Jan 8, 2015 at 14:07
• Thanks! This is exactly my problem as well. A quick follow-up: in contrast with your example, which has one factor with 4 levels, I have a 2 x 2 design where each factor has 2 levels (so the total is still 4 levels). Can I also use diag(9,4) for my model? I am not well-versed with matrices so I wanted to double-check. Relatedly, to justify this solution in my paper, should I cite Firth (1993) or is there a more directly relevant paper, which has implemented your solution using bglmer()?
– Sol
Dec 30, 2016 at 14:01
• see updated answer. Dec 30, 2016 at 19:52
• I think so - it should only matter how many fixed effect parameters there are in total. Dec 30, 2016 at 22:46