glmer model convergence question

We are working with a longitudinal dataset, with three variables: WAIP, BPSRRI and group. WAIP and BPSRRI are measured repeatedly for 10 times and group refers to the group assignment of our subjects in a treatment study.

The goal is to predict stability of WAIP score from BPSRRI of previous time point. To do this, we first imputed missing data and then created a WAIPSSD (WAIP squared successive difference variable) using the following code:

data<-read.csv("data.csv")%>%
group_by(imputation,ID)%>%
mutate(WAIPint=round(WAIP))%>%
mutate(WAIPLag=lag(WAIPint))%>%
mutate(WAIPSD=WAIPint-WAIPLag)%>%
mutate(WAIPSSD=WAIPSD^2+1)


In the last step, we added 1 to all WAIPSD (WAIP successive difference) to avoid 0s. Our goal is to predict WAIPSSD at Session t from BPSRRI at Session (t-1).

Following Jahng et al., 2008, we modeled the data with generalized multilevel model with a gamma distribution and loglink. Following Wang & Maxwell., 2015 we controlled for linear time trends by adding time variable (Session) in to the model. We modeled random slopes and intercepts for BPSRRI and fixed effects for group and BPSRRI and their interaction.

The final model is detailed below. We ran into some convergence issue:

summary(glmer(WAIPSSD~BPSRRI*gp+(1+BPSRRI|ID)+Session, data=data,
control=glmerControl(optCtrl = list(maxfun=2e9))))

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$$par, ctrl = control$$checkConv,  :
Model failed to converge with max|grad| = 0.0574011 (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?


Following this link, we rescaled and centered parameters and then checked singularity.

Rescaling and Centering:

m1<-glmer(WAIPSSD2~BPSRRI*gp+(1+BPSRRI|ID)+Session,data=data,
control=glmerControl(optCtrl = list(maxfun=2e9)))
numcols <- grep("^c\\.",names(data))
dfs <- data
dfs[,numcols] <- scale(dfs[,numcols])
m1_sc <- update(m1,data=dfs)


This did not change the warning messages. When we checked singularity, the number seems large enough (the link suggests <1.00E-6)

> tt <- getME(m1_sc,"theta")
> ll <- getME(m1_sc,"lower")
> min(tt[ll==0])
[1] 0.1771667


We tried changing the optimizer to nloptwrap but it didn't solve the convergence issue.

Increasing maxfun to 2e10 produced following warning messages:

Warning messages:
1: In (function (par, fn, lower = -Inf, upper = Inf, control = list(),  :
NAs introduced by coercion to integer range
2: In optwrap(optimizer, devfun, start, rho$$lower, control = control, : convergence code 1 from bobyqa: bobyqa -- maximum number of function evaluations exceeded 3: In nM$$setMaxeval(cc$$maxfun) : NAs introduced by coercion to integer range 4: In checkConv(attr(opt, "derivs"), opt$$par, ctrl = control$$checkConv, : unable to evaluate scaled gradient 5: In checkConv(attr(opt, "derivs"), opt$$par, ctrl = control\$checkConv,  :
Model failed to converge: degenerate  Hessian with 1 negative eigenvalues


This is a multiply imputed file and we get convergence error in some but not all imputation. Unfortunately we cannot post the data. Below is the descriptive of the WAIP variables from an imputation that did not converge:

> summary(data[which(data$$imputation==4),]$$WAIP)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
13.20   49.99   58.59   58.57   68.26   84.00
> summary(data[which(data$$imputation==4),]$$WAIPSD)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
-70.000  -3.250   0.000   0.487   4.000  36.000      60
> summary(data[which(data$$imputation==4),]$$WAIPSSD)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
1.00    5.00   17.00   83.03   65.00 4901.00      60


Does anyone have any suggestion in moving forward?

• Side note: the use of tidyverse for this example is hugely inefficient from a programming time standpoint. Dec 28, 2023 at 8:01

• For certain values of the gamma function a support including zero should be ok, like the special case of an exponential distribution. A problem is maybe when the shape parameter is below 1. For a glm, which just models the mean and variance, it should also not matter. If one looks at the glm family object Gamma()[] then there is a stop at the initialisation if any of the values are non-positive. But if you manually remove this then the fitting is not influenced, except only for the case of conputing the AIC which makes a call to dgamma and can lead to infinity.... May 9 at 7:50