0
$\begingroup$

This is the type of vague question that's likely to get downvoted, but here we go...

I am trying estimate a logistic panel regression to estimate the probability of experiencing a disease before and during the pandemic. I have an unbalanced panel of about 500,000 patients over 5 years, for a total of 1.8m observations. In each year, each patient can take a value of 0 or 1: 1 if they experienced the disease at all during the year, 0 otherwise. I'm afraid this is secure data and I cannot share it, and therefore no reproducible example.

As a baseline I estimated a normal glm, pooling all the observations, and this computed quickly producing the expected values for the regressors. I then managed to estimate a panel regression, with random intercepts for the individual patients and the only time-varying patient-level characteristic, age. It looks like this in lme4 syntax:

glmer(disease_present ~ study_year + dummy2020 + dummy2021 + age_group + (1 | individual_id),
      data = my_panel, family = "binomial")

This specification needed a bit of "help" to converge, based on the troubleshooting advice from the authors and here on StackExchange. In the end doing an additional 20,000 iterations using an alternative optimiser (nloptwrap instead of the default) did the trick, and the model emerged with very similar coefficients to the pooled model.

As a final step, I want to add the time-invariant patient characteristics back to the random effects model: sex, ethnicity, urban/rural and IMD decile (a measure of deprivation in the UK). With the exception of the last one, these all enter as categorical predictors. (I'm not expecting this to affect the coefficients on trend, since I expect them all to be accounted for in the patient random effects, but I want to check the coefficients on these fixed effects are similar to the pooled model.) Here is the full model:

glmer(disease_present ~ study_year + dummy2020 + dummy2021 + age_group + sex + ethnicity + urban_rural + imd_decile + (1 | individual_id),
      data = my_panel, family = "binomial")

At this point, using update to continue from one run to the next, I have put this model through about half a million iterations, using a range of optimisers (with 100,000 iterations taking a full day to run). It shows no sign of converging, with weird coefficients on most of the predictors and the maximum gradient / scaled gradient actually rising on recent runs (how is that possible?)

Depending on the optimiser used, I get different warnings, including:

  • Model failed to converge
  • Model is nearly unidentifiable: very large eigenvalue / large eigenvalue ratio: rescale variables?
  • Unable to evaluate scaled gradient
  • Degenerate Hessian with negative Eigenvalues

I'm running out of ideas for how to get this to work. I don't see how I can rescale categorical predictors, although I see the authors' warning that a GLMM with lots of fixed effects is particularly prone to convergence problems. I don't think the model is mis-specified given the problem, and I should have ample degrees of freedom to estimate these predictors in the random effects framework. I have the feeling that having lots of categorical predictors is less than ideal, but that's the data I have. I do have ideas for alternative fall-back models, but I'll feel like I've failed if I can't get this to run.

Does anyone have any ideas on what could be wrong here, and any ways I could get the model to converge?

$\endgroup$
2
  • $\begingroup$ Have you looked at correlations/collinearity between your predictors? $\endgroup$
    – Roland
    Oct 10, 2022 at 4:46
  • $\begingroup$ Thanks @Roland - there is some association amongst the predictors but nothing approaching perfect collinearity... $\endgroup$ Oct 19, 2022 at 15:09

1 Answer 1

1
$\begingroup$

So the short answer to this was set nAGQ = 0 as described in this SO post. We only discovered this when we successfully ran the same specification in a few minutes in Stata using its melogit command.

I think it's helpful to point out that this is just a simple parameter, it does not need to go inside a glmerControl function, so my new specification is simply:

glmer(disease_present ~ study_year + dummy2020 + dummy2021 + age_group + sex + ethnicity + urban_rural + imd_decile + (1 | individual_id),
  data = my_panel, family = "binomial", nAGQ = 0)

Whether you should simplify the problem in the way reducing the nAGQ parameter entails is a trickier debate, but there's a discussion here on Cross Validated. I'm not sure our data really meets the criteria suggested by Ben Bolker in his reply there, but we really had no choice, it was the only way to get the model to converge.

$\endgroup$

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.