2
$\begingroup$

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)))
$\endgroup$
  • $\begingroup$ 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) ... $\endgroup$ – Ben Bolker Mar 28 '15 at 23:14
  • $\begingroup$ ... (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. $\endgroup$ – Ben Bolker Mar 28 '15 at 23:15
  • $\begingroup$ @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? $\endgroup$ – exfalso Mar 29 '15 at 1:14
  • $\begingroup$ ...(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! $\endgroup$ – exfalso Mar 29 '15 at 1:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.