I am seeking help getting my latent class linear mixed models to converge using the LCMM package in R.
I am studying common patterns of heroin and cocaine use over the life course of adults who have injected drugs. My data come from an active cohort study which has collected nearly 30 years of data on the drug use behaviors of adults in Baltimore, Maryland. The goal is to identify common trajectories of use over the life course.
To do this, I want to use latent class linear mixed models that model the conditional probability of past six-month use (a binary outcome) as a function of participant age. The outcome is modeled with a threshold (probit) link function. I included linear, quadratic, and cubic age terms that can vary by latent class, as well as a few other covariates (e.g., gender, race, cohort recruitment wave, decade) that are held constant across cohorts. I included a random intercept to account for clustering of observations within participants. I fit models with as few as g=1 to as many as g=6 latent classes.
Here is the code I used to estimate a 1-class model (i.e. a simple multilevlel model):
fit.1class <- lcmm(fixed = heroin ~ age.std + age2.std + age3.std + female + black + cohort + visit.era, random = ~1, link = "thresholds", subject = "newid", data = alive.complete, verbose = TRUE)
Here is the code I used to estimate an 3-class model, using the 1-class model for starting values (note that, in the code below, i = 3).
fit <- lcmm(fixed = heroin ~ age.std + age2.std + age3.std + female + black + cohort + visit.era, random = ~1, mixture = ~ age.std + age2.std + age3.std, ng = i, link = "thresholds", subject = "newid", data = alive.complete, maxiter = 500, B = fit.1class, verbose = TRUE)
This 3-class model did not converge in 500 iterations. Further, my simple understanding is that the convergence algorithm is "getting stuck" somehow, since the second derivative convergence criteria = 1, and this did not change after running for 500 more iterations (not shown).
One possible problem is that one of the latent classes that is being estimated is essentially 100% probability of use at every age. This is actually a reasonable class -- there are some participants who never stop using -- but it seems like this could cause problems for parameter estimation since the predicted model values are at the boundary.
Any suggestions for how to get these models with g>1 classes to converge? Here are some ideas I am considering, in no particular order:
- Fit a g>1-class model without random effects, and then use that model to generate starting values (instead of a one-class model without random effects to generate starting values).
- Allow the random intercept variance to vary across classes.
- Use natural smoothing splines (as generated by ns() in the "splines" package) instead of simple linear, quadratic, and cubic terms, to force the function to be linear outside the range of the data.
- Drop the quadratic and cubic terms entirely.
- Impose a reasonable constraint on some of the model parameters. For example, I could fix the conditional probability of use at age 20 to be 95% in all classes, a reasonable constraint since participants must have injected drugs to enroll in the study.
- Reconsider whether this modeling approach is really appropriate for my data. Maybe latent classes are just not a good way to model drug use over the lifecourse?
Any suggestions you have would be greatly appreciated.