3
$\begingroup$

I fitted a mixed logit model with crossed random effects in lme4_1.1-7::glmer (R version 3.1.1 / OS X 10.9.4 Mavericks).

Had to simplify the maximal random-effect structure justified by the design due to failed convergence; the final model is estimated without any problems:

fitted_1 <- glmer(DV ~ IV1.d*IV2.d + (IV1.d*IV2.d| SubjN) + (1|Items) +  
                       (0+IV1.d|Items) + (0+IV2.d|Items) + (0+IV1.d:IV2.d|Items), 
                  glmerControl(optimizer='bobyqa', optCtrl = list(maxfun=20000)), 
                  data=myPP, family=binomial) 

DV is the binary response variable

IV1.d and IV2.d are two within-subjects within-items categorical predictors, two levels each, deviation-contrast coded (values: -.5/.5)

I tried to compute confidence intervals for the beta parameters using profile likelihood via confint.merMod() but the computation seems to be failing. For all betas, I got values (-Inf Inf) and warning messages of non-monotonic profiles. Reading on [R-sig-ME], this latter issue should mean there is something wonky with the profile.

I tried to simplify the random structure of the model until profile confidence intervals could be computed. Here is the model:

fitted_4 <- glmer(DV ~ IV1.d + IV2.d + IV1.d:IV2.d + (IV1.d + IV2.d| SubjN) +
                       (1|Items) + (0+IV1.d|Items) + (0 +IV2.d|Items), data=myPP, 
                  glmerControl(optimizer='bobyqa', optCtrl = list(maxfun=20000)),
                  family=binomial)
  1. I'm not understanding what causes the profile likelihood method to fail for the original fitted_1 model but not for fitted_4.

  2. Is there any other way I could obtain profile CI's for fitted_1?


summary(fitted_1)

##       AIC      BIC   logLik deviance df.resid 
##    1074.0   1168.1   -519.0   1038.0     1362 

##  Scaled residuals: 
##      Min      1Q  Median      3Q     Max 
##  -2.2673 -0.3611 -0.2500 -0.1378  4.5826 

##  Random effects:
##   Groups  Name        Variance  Std.Dev.  Corr             
##   SubjN   (Intercept) 2.424e+00 1.557e+00                  
##           IV1.d       1.990e+00 1.411e+00  0.17            
##           IV2.d       6.065e-01 7.788e-01 -0.97 -0.29      
##           IV1.d:IV2.d 2.172e+00 1.474e+00 -0.19 -0.81  0.39
##   Items   (Intercept) 4.615e-03 6.793e-02                  
##   Items.1 IV1.d       3.233e-13 5.686e-07                  
##   Items.2 IV2.d       9.442e-01 9.717e-01                  
##   Items.3 IV1.d:IV2.d 4.801e-01 6.929e-01                  
##  Number of obs: 1380, groups:  SubjN, 88; Items, 12

##  Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
##  (Intercept) -2.40604    0.23196 -10.373   <2e-16 ***
##  IV1.d        0.08249    0.36355   0.227   0.8205    
##  IV2.d        1.11046    0.43579   2.548   0.0108 *  
##  IV1.d:IV2.d  0.16386    0.71246   0.230   0.8181    

##  Correlation of Fixed Effects:
##              (Intr) IV1.d  IV2.d 
##  IV1.d        0.118              
##  IV2.d       -0.434 -0.083       
##  IV1.d:IV2.d -0.090 -0.628  0.064

confint(fitted_1, method="profile", which='beta_')`

##               2.5 % 97.5 %
##   (Intercept)  -Inf    Inf
##   IV1.d        -Inf    Inf
##   IV2.d        -Inf    Inf
##   IV1.d:IV2.d  -Inf    Inf

##  Warning messages:
##  1: In profile.merMod(object, signames = oldNames, ...) :
##    non-monotonic profile
##  2: In profile.merMod(object, signames = oldNames, ...) :
##    non-monotonic profile
##  3: In profile.merMod(object, signames = oldNames, ...) :
##    non-monotonic profile
##  4: In profile.merMod(object, signames = oldNames, ...) :
##    non-monotonic profile

summary(fitted_4)

##       AIC      BIC   logLik deviance df.resid 
##      1068     1136     -521     1042     1367 

##  Scaled residuals: 
##      Min      1Q  Median      3Q     Max 
##  -2.3575 -0.3555 -0.2522 -0.1613  4.6391 

##  Random effects:
##   Groups  Name        Variance Std.Dev. Corr       
##   SubjN   (Intercept) 2.23144  1.4938              
##           IV1.d       1.53606  1.2394    0.09      
##           IV2.d       0.31120  0.5579   -1.00 -0.18
##   Items   (Intercept) 0.01344  0.1159              
##   Items.1 IV1.d       0.00000  0.0000              
##   Items.2 IV2.d       0.92942  0.9641              
##  Number of obs: 1380, groups:  SubjN, 88; Items, 12

##  Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
##  (Intercept) -2.30252    0.21029 -10.949   <2e-16 ***
##  IV1.d        0.17448    0.27729   0.629   0.5292    
##  IV2.d        0.80072    0.36862   2.172   0.0298 *  
##  IV1.d:IV2.d  0.01351    0.41660   0.032   0.9741    

##  Correlation of Fixed Effects:
##              (Intr) IV1.d  IV2.d 
##  IV1.d        0.012              
##  IV2.d       -0.274 -0.010       
##  IV1.d:IV2.d  0.006 -0.255 -0.038

confint(fitted_4, which='beta_', method='profile')

##                    2.5 %     97.5 %
##  (Intercept) -2.74571641 -1.9052351
##  IV1.d       -0.37989551  0.7320931
##  IV2.d        0.03993436  1.5903197
##  IV1.d:IV2.d -0.80790153  0.8346440

UPDATE

## re-compute profiles for both random and fixed effects

pp <- profile(fitted_1)

## 24 warnings with profile(fitted_1) of the types:
## In profile.merMod(fitted_1) : non-monotonic profile
## In optwrap(optimizer, par = start, fn = function(x) dd(mkpar(npar1,  ... :
   # convergence code 1 from bobyqa: bobyqa -- maximum number of function evaluations exceeded

(c_ci <- confint(pp))

## 2.5 % 97.5 %
## .sig01          0    Inf
## .sig02         -1      1
## .sig03         -1      1
## .sig04         -1      1
## .sig05          0    Inf
## .sig06         -1      1
## .sig07         -1      1
## .sig08          0    Inf
## .sig09         -1      1
## .sig10          0    Inf
## .sig11          0    Inf
## .sig12          0    Inf
## .sig13          0    Inf
## .sig14          0    Inf
## (Intercept)  -Inf    Inf
## IV1.d        -Inf    Inf
## IV2.d        -Inf    Inf
## IV1.d:IV2.d  -Inf    Inf


## plot the profiles (all weird)
    ggplot(as.data.frame(pp),aes(.focal,.zeta))+
    geom_point()+geom_line()+
    facet_wrap(~.par,scale="free_x")+
    geom_hline(yintercept=0,colour="gray")+
    geom_hline(yintercept=c(-1.96,1.96),linetype=2,
               colour="gray")

Plot_profile_fitted1.jpeg

### setting delta to a smaller value to make the profile stepsize smaller
system.time(pp2 <- profile(fitted_1, delta = 0.1))

## user    system   elapsed 
## 64292.282   135.451 75676.403

## Warning messages:
## 1: In profile.merMod(orig.pp, delta = 0.1) : non-monotonic profile
## 2: display list redraw incomplete
## [...]

c_ci2 <- confint(pp2)
c_ci2

## 2.5 % 97.5 %
## .sig01          0    Inf
## .sig02         -1      1
## .sig03         -1      1
## .sig04         -1      1
## .sig05          0    Inf
## .sig06         -1      1
## .sig07         -1      1
## .sig08          0    Inf
## .sig09         -1      1
## .sig10          0    Inf
## .sig11          0    Inf
## .sig12          0    Inf
## .sig13          0    Inf
## .sig14          0    Inf
## (Intercept)  -Inf    Inf
## IV1.d        -Inf    Inf
## IV2.d        -Inf    Inf
## IV1.d:IV2.d  -Inf    Inf


## plot of profiles (delta = 0.1)

ggplot(as.data.frame(pp2),aes(.focal,.zeta))+
    geom_point()+geom_line()+
    facet_wrap(~.par,scale="free_x")+
    geom_hline(yintercept=0,colour="gray")+
    geom_hline(yintercept=c(-1.96,1.96),linetype=2,
               colour="gray")

Plot of profiles delta = 0.1

$\endgroup$
  • $\begingroup$ (1) in what precise way is the computation failing (error message?); (2) have you tried plotting the profile? (3) singularity shouldn't make the profile impossible, but it would be good to check for (th <- getME(fitted,"theta"); lwr <- getME(fitted,"lower"); any(th[lwr==0]<1e-7). (4) Depending on what the profile looks like, you might try setting delta to a small value (0.1?) to make the profile stepsize smaller ... $\endgroup$ – Ben Bolker Oct 7 '14 at 18:51
  • $\begingroup$ @Ben Bolker - Thanks very much for your reply. <br/> 1. No error message, only warnings. Yet profile confidence intervals for all beta parameters got values (-Inf +Inf) which seems a bit wonky, plus it's something I cannot report. 2. Yes, see UPDATE above. They do seem a bit suspicious... 3. FALSE, so singularity shouldn't be at stake. 4. Tried, see results for pp2 under UPDATE above. Not much different from result from the original profiles computation (pp), yet the profile plots look even more weird... Any idea what's going on / how I can deal with it? Tnx! $\endgroup$ – exfalso Oct 9 '14 at 7:49

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.