# Can complete separation between a continuous predictor and a random effect cause failure to converge in a logit GLMM?

I’m running a logit mixed-effects model on binary data with a 2x2 within-subjects design, with subjects and items as crossed random effects, and the two independent variables deviation-contrast coded.

Here are model specification and summary:

mod1 <- glmer(DV ~ devX1*devX2 + (devX1*devX2|Subject) + (devX1*devX2|Item),
data=mydata, family=binomial, glmerControl(optimizer='bobyqa',
optCtrl=list(maxfun=400000)))

AIC      BIC   logLik deviance df.resid
628.9    734.3   -290.4    580.9      573

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.0527 -0.5025 -0.2217  0.5654  4.0493

Random effects:
Groups  Name        Variance Std.Dev. Corr
Subject (Intercept) 0.1184   0.3440
devX1       3.5387   1.8812   -0.74
devX2       0.2461   0.4961   -0.54  0.06
devX1:devX2 4.5912   2.1427    0.32 -0.84  0.07
Item    (Intercept) 0.5568   0.7462
devX1       0.2693   0.5190    0.48
devX2       0.3862   0.6215   -0.31 -0.51
devX1:devX2 2.2109   1.4869   -0.57  0.42 -0.31
Number of obs: 597, groups:  Subject, 30; Item, 20

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.47781    0.27602  -5.354 8.60e-08 ***
devX1        2.70622    0.55692   4.859 1.18e-06 ***
devX2        0.08229    0.45801   0.180    0.857
devX1:devX2 -0.41055    0.99645  -0.412    0.680
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) devX1  devX2
devX1       -0.498
devX2       -0.179  0.046
devX1:devX2 -0.021 -0.266 -0.657


The model does converge with full random structure without any problems. (It may be worth mentioning that the binned Pearson residual plot reveals that the model has some issues accounting for y = 0 original data points.)

I'm encountering big convergence issues as soon as I include a centered continuous covariate (Age) as fixed effect. It does not matter how much I simplify the random structure, the model will not converge.

mod1.age <- glmer(DV ~ devX1*devX2*cAge + (devX1*devX2|Subject) + (devX1*devX2|Item),
data=mydata, family=binomial, glmerControl(optimizer='bobyqa',
optCtrl=list(maxfun=400000)))

AIC      BIC   logLik deviance df.resid
624.4    747.4   -284.2    568.4      569

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.9512 -0.5140 -0.2234  0.5361  5.3189

Random effects:
Groups  Name        Variance  Std.Dev.  Corr
Subject (Intercept) 1.037e-11 3.220e-06
devX1       2.692e+00 1.641e+00  0.28
devX2       3.864e-02 1.966e-01  0.08 -0.94
devX1:devX2 4.489e+00 2.119e+00 -0.50 -0.97  0.82
Item    (Intercept) 5.280e-01 7.267e-01
devX1       2.662e-01 5.159e-01  0.79
devX2       3.948e-01 6.284e-01 -0.36 -0.48
devX1:devX2 2.906e+00 1.705e+00 -0.59  0.02 -0.22
Number of obs: 597, groups:  Subject, 30; Item, 20

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)      -1.3832677  0.0019843  -697.1  < 2e-16 ***
devX1             2.4397103  0.0020486  1190.9  < 2e-16 ***
devX2             0.1386076  0.0019838    69.9  < 2e-16 ***
cAge             -0.0091753  0.0016630    -5.5 3.44e-08 ***
devX1:devX2      -0.3524321  0.0028066  -125.6  < 2e-16 ***
devX1:cAge        0.0150530  0.0019310     7.8 6.41e-15 ***
devX2:cAge        0.0121991  0.0018876     6.5 1.03e-10 ***
devX1:devX2:cAge  0.0005894  0.0019504     0.3    0.763
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) devX1  devX2  cAge   dvX1:X2 dvX1:A dvX2:A
devX1       -0.001
devX2       -0.001  0.001
cAge         0.002 -0.002 -0.002
devX1:devX2 -0.002  0.001  0.001 -0.001
devX1:cAge  -0.002  0.001  0.001 -0.043  0.002
devX2:cAge  -0.001  0.001  0.001 -0.040  0.001  -0.017
dvX1:dvX2:A  0.001 -0.001 -0.001 -0.019 -0.001  -0.020 -0.009


There is no problem of complete or quasi complete separation between Age and the binary DV. However, there is an almost perfect 1:1 match between Age and Subject (with Subject specified as a random effect in the models). In other words, for most values of Age, there is only one subject corresponding to that value, which makes Age a sort of another version of Subject.

Could this be what is causing severe convergence problems?

If so, would transforming Age into a categorical variable (e.g., with 3 levels) be a suitable solution? I would like to avoid largely arbitrary choices about model specification.

What makes me doubt about this explanation though is that if I remove Subject as random effect, the resulting model still fails to converge.

mod1.age4 <- glmer(DV ~ devX1*devX2*cAge + (devX1*devX2|Item), data=mydata,
family=binomial, glmerControl(optimizer='bobyqa',
optCtrl=list(maxfun=400000)))

• A few comments and questions: (1) I can't quite tell how many observations you have per subject (the total number of obs. is approx (#subject)*(#items), suggesting each subject receives all items, but then I'm not sure how the treatments are allocated within subjects -- a little more info about your design? (2) How big are your scaled gradient estimates? Maybe false positive -- see rpubs.com/bbolker/lme4trouble1 (3) Binning age does seem reasonable. (4) note that your subject RE fit is singular (intercept var is ~ 0) ... – Ben Bolker Mar 28 '15 at 23:14
• ... (2/2) Fitting 10 var-cov parameters to grouping variables with 20 and 30 levels is very ambitious -- equivalent to trying to fit a 4x4 variance-covariance matrix to that many independent data points. I imagine you're following the Barr et al. "keep it maximal" prescription, but still -- difficult. – Ben Bolker Mar 28 '15 at 23:15
• @BenBolker (1) 20 obs per subject, so 5 obs per cell of the 2x2 within-subjects design per subject. Subjects and items are crossed, each subject sees each item only once. 8 obs per design cell per item. (2) Based on pre-computed information (mod1.age1$optinfo$derivs), max(abs(sc_grad1)) = .0723 and max(pmin(abs(sc_grad1),abs(derivs1\$gradient))) = 0.0723, so rather large. Also, using numDeriv, I get very different results max(pmin(abs(sc_grad2),abs(grad2))) = 5.1735. What does it all mean exactly? – exfalso Mar 29 '15 at 1:14
• ...(2/2) (4) Is a RE var ~ 0 always sign of singularity (even if the model does converged)? And following the singularity of subject RE fit, how would you recommend to simplify the random structure? Forcing independence between by-subject random effects (i.e., ... + (devX1*devX2||Subject) ) or removing higher-order random term (i.e., ... + (devX1 + devX2 | Subject) ) ? (5) Extra info, on top of the warning failed convergence message regarding the max|grad|, I also get the one on degenerate Hessian with 4 negative eigenvalues. Many thanks! – exfalso Mar 29 '15 at 1:18