# Different estimates for mixed effects Logistic regression and pwrssUpdate Error message with binomial glmer

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)
 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!

## 1 Answer

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) • 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. – Tiantian Yang Dec 4 '20 at 3:30
• 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. – Tiantian Yang Dec 5 '20 at 16:01
• 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). – Ben Bolker Dec 5 '20 at 23:13
• Are you still getting unsatisfactory results from the fit to the full data set? – Ben Bolker Dec 5 '20 at 23:13
• Hello Dr. Bolker, thanks for your quick reply. It is weird that you got different results on a different platform [Linux], could you please post your result here so that I can have a look? Will the R results sometimes be sensitive to a different platform? – Tiantian Yang Dec 6 '20 at 14:47