0
$\begingroup$

Hi this is my first post here. I reported an issue on Github and got valuable help from Dr. Bolker. However, I did some other analyses on the same data and got confused about some of the results. I want to know what causes these issues to happen, and I am stuck. I read through some other resources written by Dr. Boller. I found them are very useful, such as lme4 convergence warnings: troubleshooting, GLMM FAQ, Generalized linear mixed models in R: nitty-gritty. Besides that, I also refer to some useful posts such as post 1, post 2, post 3, post 4, post 5, etc.

Here I posted some of the results below, as well as my questions. Here is the data. I will deeply appreciate any comments on one or all of the issues.

The data structure looks like this, and there is no complete separation issue.

> dim(data)
[1] 75  4
> head(data, 10)
   y            x1           x2 group
1  1   0.818448879  4.060243218     1
2  1   7.990440255  1.857443185     1
3  1  -6.595283621  4.715403083     2
4  0   0.675370785 -3.273230423     2
5  0 -16.391619950  9.016722634     3
6  1  -1.124991928  3.728152698     3
7  1   8.931848938  0.784097814     3
8  0   3.445347058 -4.436738943     3
9  1  -3.969142249  4.440858374     4
10 0   2.033157618  0.850871635     4
> glm(y ~ 1 + x1 + x2, data = data, family = binomial, method="detect_separation")
Separation: FALSE 
Existence of maximum likelihood estimates
(Intercept)          x1          x2 
          0           0           0 
0: finite value, Inf: infinity, -Inf: -infinity
Warning message:
'detect_separation' will be removed from 'brglm2' at version 0.8. A new version of 'detect_separation' is now maintained in the 'detectseparation' package.  
> table(table(data$group)) 

 1  2  3  4 
 6 10  7  7

I changed the combination of two optimizers (the default and bobyqa) and three integration methods (Laplace, Penalized quasi-likelihood(PQL), adaptive Gauss-Hermite quadrature(AGQ)) when fitting the data using binomial glmer.

#### Setting 1: Optimization technique: Nelder-Mead and bobyqa (default) & Integration method: Laplace (defualt)
summary(glmer1 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial)) # Here nAGQ =1 (default) which is laplace approximation
# error messsage

#### Setting 2: Optimization technique: Nelder-Mead and bobyqa (default) & Integration method: Penalized quasi-likelihood (PQL)
summary(glmer2 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial, nAGQ=0)) # Let nAGQ =0. 

#### Setting 3: Optimization technique: Nelder-Mead and bobyqa (default) & Integration method: adaptive Gauss-Hermite quadrature (nAGQ=20)
summary(glmer3 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial, nAGQ=20))

#### Setting 4: Optimization technique: bobyqa & Integration method: Laplace (default) 
summary(glmer4 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial, control=glmerControl(optimizer="bobyqa")))
# warning message 

## Setting 5: Optimization technique: bobyqa & Integration method: Penalized quasi-likelihood 
summary(glmer5 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial, nAGQ=0, control=glmerControl(optimizer="bobyqa")))
# same result as glmer2

## Setting 6: Optimization technique: bobyqa & Integration method: adaptive Gauss-Hermite quadrature (nAGQ=20)
summary(glmer6 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial, nAGQ=20, control=glmerControl(optimizer="bobyqa")))
# same result as glmer3 (Note: not exactly same, but same parameter estimates till the fifth decimal places)

The default setting 1 causes an error as below.

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev, : (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

All the other settings work with only setting 4 produces warning messages as below.

Warning message: In checkConv(attr(opt, "derivs"), opt$par, ctrl = > control$checkConv, : Model failed to converge with max|grad| = 0.0262097 (tol = 0.001, component 1)

Question 1: Settings 2 and 5 both use the PQL integration method, and there are no warnings or errors for the analysis results. However, their analysis results are different from Setting 3 and 6, which use the adaptive Gauss-Hermite quadrature technique. How do you know which method is more accurate?

Setting 2 has the same result as setting 5. The output is as below

> glmer2
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data
     AIC      BIC   logLik deviance df.resid 
 66.3493  75.6193 -29.1747  58.3493       71 
Random effects:
 Groups Name        Std.Dev. 
 group  (Intercept) 0.6954595
Number of obs: 75, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
 -0.8870778    0.3750967    0.7150721  
>
> glmer5
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data
     AIC      BIC   logLik deviance df.resid 
 66.3493  75.6193 -29.1747  58.3493       71 
Random effects:
 Groups Name        Std.Dev. 
 group  (Intercept) 0.6954595
Number of obs: 75, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
 -0.8870778    0.3750967    0.7150721

Setting 3 has the same result as setting 6 till the fifth decimal places.

> glmer3
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 20) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data
     AIC      BIC   logLik deviance df.resid 
 65.8476  75.1175 -28.9238  57.8476       71 
Random effects:
 Groups Name        Std.Dev.
 group  (Intercept) 2.163064
Number of obs: 75, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
 -1.1612282    0.6115683    1.1224949  
> glmer6
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 20) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data
     AIC      BIC   logLik deviance df.resid 
 65.8476  75.1175 -28.9238  57.8476       71 
Random effects:
 Groups Name        Std.Dev.
 group  (Intercept) 2.163058
Number of obs: 75, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
 -1.1612252    0.6115669    1.1224928

Question 2: There are different analysis results even after deleting only one other row of data. But deleting one row will not change the data structure. What causes this to happen? Code is as follows.

## Update setting 1
# update setting 1 by delete row 55 in data
glmer(y ~ 1 + x1 + x2 + (1|group), data = data[-55,], family = binomial) 
# error message

# update setting 1 by delete row 60 in data
glmer(y ~ 1 + x1 + x2 + (1|group), data=data[-60,], family=binomial)
# warning message

# update setting 1 by delete row 63 in data
glmer(y ~ 1 + x1 + x2 + (1|group), data=data[-63,], family=binomial)
# boundary (singular) fit

The output is as follows.

> glmer(y ~ 1 + x1 + x2 + (1|group), data = data[-55,], family = binomial) 
Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  : 
  pwrssUpdate did not converge in (maxit) iterations
> 
> glmer(y ~ 1 + x1 + x2 + (1|group), data=data[-60,], family=binomial)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data[-60, ]
     AIC      BIC   logLik deviance df.resid 
 46.5400  55.7562 -19.2700  38.5400       70 
Random effects:
 Groups Name        Std.Dev.
 group  (Intercept) 113.7748
Number of obs: 74, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
  -20.29790     15.95948     29.44503  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
> 
> glmer(y ~ 1 + x1 + x2 + (1|group), data=data[-63,], family=binomial)
boundary (singular) fit: see ?isSingular
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data[-63, ]
     AIC      BIC   logLik deviance df.resid 
 65.0592  74.2754 -28.5296  57.0592       70 
Random effects:
 Groups Name        Std.Dev.         
 group  (Intercept) 0.000000009351224
Number of obs: 74, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
 -0.8793082    0.3655369    0.7093192  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings

As we know that setting 3 (nAGQ=20) works fine with the data, but it also produces different results after deleting one other row. The code is as below.

## Update setting 3
# update setting 3 by delete row 55 in data
summary(update(glmer3, data=data[-55,]))

# update setting 3 by delete row 60 in data
summary(update(glmer3, data=data[-60,]))

# update setting 3 by delete row 63 in data
summary(update(glmer3, data=data[-63,]))
# boundary (singular) fit

And output is as follows

> update(glmer3, data=data[-55,])
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 20) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data[-55, ]
     AIC      BIC   logLik deviance df.resid 
 65.8068  75.0230 -28.9034  57.8068       70 
Random effects:
 Groups Name        Std.Dev.
 group  (Intercept) 2.141286
Number of obs: 74, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
 -1.1577292    0.6048679    1.1109775  
>
> update(glmer3, data=data[-60,])
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 20) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data[-60, ]
     AIC      BIC   logLik deviance df.resid 
 59.3522  68.5684 -25.6761  51.3522       70 
Random effects:
 Groups Name        Std.Dev.
 group  (Intercept) 1.514945
Number of obs: 74, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
 -1.0437079    0.6169493    1.1426508 
> 
> update(glmer3, data=data[-63,])
boundary (singular) fit: see ?isSingular
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 20) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data[-63, ]
     AIC      BIC   logLik deviance df.resid 
 65.0592  74.2754 -28.5296  57.0592       70 
Random effects:
 Groups Name        Std.Dev.         
 group  (Intercept) 0.000000009351224
Number of obs: 74, groups:  group, 30
Fixed Effects:
(Intercept)           x1           x2  
 -0.8793082    0.3655369    0.7093192  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings

Question 3: I restarted the default setting 1 (Laplace) from the estimates produced by other settings. Restarting with the estimates from setting 3 works fine with some warning messages while restarting with the estimates from setting 6 causes an error. However, the parameter estimates between setting 3 and setting 6 are the same till the fifth decimal places. Why does this happen? The code is as follows.

## restart setting 1 with estimates from setting 3
summary(glmer1.3 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial,
                      start=getME(glmer3, c("theta","fixef")))) 
# warning message

## restart setting 1 with estimates from setting 6
summary(glmer1.6 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial,
                      start=getME(glmer6, c("theta","fixef")))) 
# error message

The output is below.

> summary(glmer1.3 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial,
+                       start=getME(glmer3, c("theta","fixef")))) 
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: y ~ 1 + x1 + x2 + (1 | group)
   Data: data

     AIC      BIC   logLik deviance df.resid 
    50.0     59.2    -21.0     42.0       71 

Scaled residuals: 
         Min           1Q       Median           3Q          Max 
-0.127764107 -0.000000015 -0.000000015  0.000000015  0.109007665 

Random effects:
 Groups Name        Variance Std.Dev.
 group  (Intercept) 13832.97 117.6137
Number of obs: 75, groups:  group, 30

Fixed effects:
                  Estimate     Std. Error   z value               Pr(>|z|)    
(Intercept) -19.4439989773   0.0009608065 -20237.17 < 0.000000000000000222 ***
x1           15.4149024304   0.0009644601  15982.93 < 0.000000000000000222 ***
x2           28.7509079765   0.0009637793  29831.42 < 0.000000000000000222 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
   (Intr) x1    
x1  0.000       
x2  0.000 -0.017
convergence code: 0
Model failed to converge with max|grad| = 0.0736003 (tol = 0.001, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.0736003 (tol = 0.001, component 1)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?
> summary(glmer1.6 <- glmer(y ~ 1 + x1 + x2 + (1|group), data=data, family=binomial,
+                       start=getME(glmer6, c("theta","fixef")))) 
Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  : 
  (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

Question 4: Restarting setting 2 and setting 5 (PQL) with their fits cause an error below.

Error in glmer(formula = y ~ 1 + x1 + x2 + (1 | group), data = data, family = binomial, : should not specify both start$fixef and nAGQ==0

Why can't we specify the starting values when using PQL integration? The code is below.

## restart setting 2 with its fits
summary(glmer2.2 <- update(glmer2, start=getME(glmer2, c("theta","fixef"))))
# error message

## restart setting 5 with its fits
summary(glmer5.5 <- update(glmer5, start=getME(glmer5, c("theta","fixef"))))
# error message

Question 5: When restarting the other settings with the estimates produced by setting 4, the following settings all make the error message.

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  : 
  (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

Why does this happen? Code as below.

## restart setting 3 with estimates from setting 4
summary(glmer3.4 <- update(glmer3, start=getME(glmer4, c("theta","fixef")))) 
# error message

## restart setting 4 with estimates from setting 4
summary(glmer4.4 <- update(glmer4, start=getME(glmer4, c("theta","fixef"))))
# error message

# restart setting 6 with estimates from setting 4
summary(glmer6.4 <- update(glmer6, start=getME(glmer4,c("theta","fixef"))))
# error message

Other Questions: What causes the inconsistency and instability between different settings? Is it due to the data structure or the relationship between the variables? The likelihood function is of no closed-form and is it possible to write down the gradient or Hessian matrix formula? Any other approaches could be addressed?

Thanks!

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
2
$\begingroup$

An incomplete but hopefully useful set of points:

  • nAGQ=0 does not denote PQL. From Doug Bates: "nAGQ=0 doesn't skip the integral, it just includes the fixed effects in the PIRLS optimization process. Both that and nAGQ=1 use the Laplace approximation to the integral." In other words, it uses the Laplace approximation but uses a faster, approximate computation for the optimization (nAGQ=1 does a full nonlinear search over the combination of fixed-effect and random-effect of parameters.)

  • Worrying about lower values of nAGQ is a waste of time and energy. When you have a problem that requires higher numbers of Gauss-Hermite quadrature points, i.e. higher nAGQ settings, the higher-nAGQ answers are simply better approximations to the true fit of the model you have specified. Lower-nAGQ fits are faster approximations that are sometimes perfectly adequate. When they're not, they're not.

    How do you know when you need higher numbers of quadrature points? (1) when you have little information per cluster, e.g. binary data with small numbers of observation per cluster; (2) whenever increasing the nAGQ setting affects your answer. Both of these criteria are satisfied in your case.

  • It is indeed true that closed-form expressions are not available for the log-likelihood or its gradient. In fact, in most cases there aren't even convenient numerical approaches for computing the gradient (they do exist for some cases of linear mixed models, but rarely for GLMMs), so we generally use derivative-free forms of nonlinear optimization to solve the problem. This complexity makes it very hard to answer questions like "why is the fit sensitive to [some particular setting or perturbation]?"!

  • you can easily confirm that all of the available optimizers give similar answers for the full data set with nAGQ=20:

library(lme4)
dd <- read.csv("CV499269.csv")
glmer3 <- glmer(y ~ 1 + x1 + x2 + (1|group),
                data=dd, family=binomial, nAGQ=20)
ff <- allFit(glmer3)
summary(ff)$fixef
                              (Intercept)        x1       x2
bobyqa                          -1.161224 0.6115668 1.122493
Nelder_Mead                     -1.161225 0.6115668 1.122493
nlminbwrap                      -1.161227 0.6115670 1.122493
optimx.L-BFGS-B                 -1.161239 0.6115761 1.122509
nloptwrap.NLOPT_LN_NELDERMEAD   -1.161271 0.6115621 1.122459
nloptwrap.NLOPT_LN_BOBYQA       -1.161217 0.6114868 1.122337

(The summary(ff) object contains lists of various other information about the fits - log-likelihood, random-effects estimates, etc.)

  • you can use ii <- influence(glmer3); car::infIndexPlot(ii) to explore the effects of case deletion more thoroughly.
  • I don't know why your model fits are sensitive to case deletion (on a different platform [Linux], the results are almost identical when deleting cases 55, 60, 63). The points you've chosen don't appear particularly weird/outlier-ish. (One trick that's always worth a shot is scaling and centering your predictors; it might not make a difference, but it sometimes helps and it's easy, so it's worth a try.)
dd$nm <- rownames(dd)
library(ggplot2); theme_set(theme_bw())
library(ggalt)

ggplot(dd,aes(x1,x2,shape=factor(y))) +
    geom_point(size=3) +
    scale_shape_manual(values=c(1,16)) +
    geom_encircle(aes(group=group),colour="gray")+
    geom_text(data=dd[c(55,60,63),], aes(label=nm), colour="red",
              vjust=-1.5)

enter image description here

Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ Thanks so much for your valuable advice, Dr. Bolker. I will get back to you if I have any questions after I try out your suggestions. $\endgroup$
    – Tiantian Yang
    Dec 4, 2020 at 3:30
  • $\begingroup$ Hi Dr. Bolker, thanks again for your help. Here I have a question. From the above plot that you drew, it seems that almost all the y=1 points are on the upper right while the y=0 points are on the lower left. Is this also some type of separation? I wonder if that is the reason that causes the error message. One thing we found is that when the "x1" values are randomly permuted, the error does not occur, and in fact there is a singular fit. $\endgroup$
    – Tiantian Yang
    Dec 5, 2020 at 16:01
  • $\begingroup$ I don't know. It's true that correlation between parameters will make numerical procedures more unstable, but this pattern doesn't look that extreme. You could try orthogonalizing the model matrix (i.e., center x1 and x2, and compute their principal components; use PC1 and PC2 as your predictor variables; use the info from the PCA to convert the coefficients back to the original variables). $\endgroup$
    – Ben Bolker
    Dec 5, 2020 at 23:13
  • 1
    $\begingroup$ OK. It's just that these are very difficult questions to answer. I can give you general suggestions to improve numerical stability (scale, center, orthogonalize); you could also look at the structure of the (log)likelihood surface (i.e. compute the log-likelihood over a range of {beta0, beta1, beta2, theta} values and visualize the results. But this is all getting beyond the scope of a CrossValidated question. $\endgroup$
    – Ben Bolker
    Dec 6, 2020 at 22:22
  • 1
    $\begingroup$ You can also cross-check with other estimation platforms: both the glmmTMB and the GLMMadaptive packages should be able to handle this model $\endgroup$
    – Ben Bolker
    Dec 6, 2020 at 22:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.