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:


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,
              family = Gamma(link = "log"), 
              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:

              family = Gamma(link = "log"), 
              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 

WAIP Histogram WAIPSD Histogram WAIPSSD Histogram Does anyone have any suggestion in moving forward?

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

1 Answer 1



This is Emily from the SCC by the way. I think [1] has some helpful information for convergence issues in general with these types of statistical models.

I believe I found the root of your issue: you have a great deal of WAIPSSD scores equal to 0. This is your response variable, which is supposed to follow a Gamma distribution given the data. However, the support of the Gamma distribution is (0, infinity) which does not include 0. The mass number of 0's in your model mean the Gamma regression model is not appropriate for your data.

  • $\begingroup$ 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.... $\endgroup$ May 9 at 7:50
  • $\begingroup$ .... A problem might be when the model starts estimating the mean and variance as zero, then the weights become infinity. This is a bit similar to logistic regression breaking when there is perfect separation. $\endgroup$ May 9 at 7:52

Your Answer

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

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