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!