0
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ What's the format of the variables such as county, date, labcaseno, etc? Are they factors? Are you able to fit a glm without random effects? How do you tie an overdose to an individual bust? $\endgroup$ – jsk Jul 16 at 2:49
  • $\begingroup$ County is categorical, with 88 levels. Date is categorical with 60 levels, but glmer is able to fit it as quantitative too (1-60), and LabCaseNo is categorical with many, many levels (probably about 40,000). When we tried glmer without random effects, the number of dummy variables breaks the model. We also can't tie individual overdoses to individual busts. Only deaths (per month and per county) is available. $\endgroup$ – David White Jul 16 at 3:41
  • $\begingroup$ Labcaseno with that many levels will probably break glmer as well. If you only have deaths per month and county, then how can you predict deaths from labcaseno? $\endgroup$ – jsk Jul 16 at 5:24
  • $\begingroup$ @jsk I think you are hitting on the main problem. My method had been to make a column, Deaths, containing the number of deaths in the month when the drug seizure happened. But, because we can easily have 2,000 seizures in a month, that death value is repeated 2,000 times. The answer below says this is what's breaking glmer. $\endgroup$ – David White Jul 16 at 14:45
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks for your answer. We did both of these suggestions earlier in our paper. We began with a "plain, boring old lm()" with total_fentanyl_seizures, and were advised not to use total_grams because one high weight seizure should not be the same as a bunch of low weight seizures (indeed, low weight seizures tend to be deadlier, because of how fentanyl is mixed into the drug supply). In the next section, we fit an ARIMA model with deaths and total_fentanyl_seizures. We wanted the GLMM as an omnibus model to finish the paper. Can you suggest ways to correct the 2 problems and still use a GLMM? $\endgroup$ – David White Jul 16 at 3:39
  • $\begingroup$ I agree with John Davis - you shouldn't be repeating your outcome over and over again within one model. Some further feature engineering of your "plain, boring old lm()" should fix this problem. For example, you could create a new variable called "number of low weight seizures". Use your expert knowledge to decide on appropriate predictors and thresholds when defining them. You just need to make sure you're clear about how each predictor was defined, and give good explanations for how/why you defined them as such, when reporting your results. $\endgroup$ – rw2 Jul 16 at 13:17
  • $\begingroup$ @rw2, thanks for your feedback. Since we've already demonstrated the importance of weight in an earlier section, I might end up just deleting this GLMM section. $\endgroup$ – David White Jul 16 at 14:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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