Question about GLMM in R

Question

We have a dataset containing drug seizures every month from 2014 - 2018 in Ohio (60 months). For each drug seizure, we have the following variables in a dataframe I’m calling df, with some 130,000 rows.

• Date = month/year it was seized,
• Weight = the weight of the seizure in grams,
• County = the county in Ohio where the seizure occurred, the chemical composition (only used to create the variable Have_Fentanyl),
• Have_Fentanyl = 1 if fentanyl is present and 0 otherwise, LabCaseNo = lab case identifier

If two drugs were seized in the same bust, they have the same lab case identifier. We also have deaths due to drug overdose every month/year, in a column with a zillion repetitions (because 60 values spread out over 130,000 rows).

We are trying to fit a generalized linear mixed model to predict deaths from the other variables. Because seizures that happened as part of the same bust are not independent, we have one level for seizures, one for busts, one for county, and one for month/year. In R, we fit this model via the following code

library(lme4)
summary(glmer(Death ~ log(Weight) + Have_Fentanyl 
    +  (1 | LabCaseNo) + (1 | County) + Date,
    data = df[complete.cases(df$Death),], 
    control = glmerControl(optimizer = "bobyqa"), 
    family = "poisson"))

(1) As the response variable is count data, we are using Poisson residuals. How would we modify this code to make this quasi-Poisson or negative binomial?

When we run this code we get an error that says Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables? When we start dropping variables like Date and County we get other error messages like Model failed to converge with max|grad| = 2.64655 (tol = 0.001, component 1) Model is nearly unidentifiable: very large eigenvalue - Rescale variables?

We get these errors for every choice of optimizer (“bobyqa” above) that we have tried. The source we learned about the GLMM model from also had “failed to converge” errors, and ignored them.

(2) Do we need to worry about any of these errors?

(3) Why are we getting these errors? Should we use a different GLMM command?

We’ve already taken the log of our weight variable. All the other variables are categorical. When we run the model on only 100 data points we don’t get errors, so it might have to do with the number of rows.

Lastly, we are a bit confused about which variables should be “fixed” effects and which “random” effects. We’ve read other posts on Cross Validated, blogs, course webpages, and the writings of Gelman. I had envisioned a four level model with a level for time, one for county, and one for LabCaseNo, but I don’t actually know how to fit that model with glmer or lmer. If anyone has a tip please share!

Answer 1

I think there are two problems going on here. First, I'll hazard a guess that your model is far too complex—all of these nestings are probably the source of the errors. Second, and more importantly, you are fitting different data to the same outcome variables many times over.

Your outcome (please correct me if I am wrong here) is drug overdose deaths in Ohio during a given month, so you have 60 outcome datapoints. Also, while these deaths per month are technically count data, if they are "big enough" you could reasonably ditch the poisson assumption.

I think a better strategy might be to try a little feature engineering to make some new variables that are more informative.

For example, why not transform your variables weight, county and date from their current form (130,000 rows) into something that summarizes the total amount of drugs seized across the state in a particular month, like total_grams_of_drugs_siezed? (forgive my verbose variable names). And maybe sum up all of your Have_Fentanyl instances during a given month into a new variable total_fentanyl_seizures?

Then you could have something like deaths ~ total_grams + total_fentanyl_seizures in a plain, boring old lm(). That's a good start. If you wanted to get fancy, you could apply some longitudinal techniques (a GAM with a s(time) term, maybe?) but I'd be worried that you're actually regressing out something you care about—yes, drug overdose deaths are changing over time, but don't you think they might be changing because there are more fentanyl-laced drugs?

A final note, another type of variable that you might want to try constructing would be a "time-lagged" variable—drug deaths in October might be predicted by drug deaths in September, for example.