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,
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:
m1<-glmer(WAIPSSD2~BPSRRI*gp+(1+BPSRRI|ID)+Session,data=data,
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
Does anyone have any suggestion in moving forward?
