How do I fit generalised linear multilevel models including offsets?

Specific Domains Statistics
#1

I first heard of Julia a couple of years ago. It was presented as an alternative to R, SAS, and Stata.
I had a look and I really liked it.The syntax just made sense. I could see where I might need to go to other languages, but it seemed worth pursuing.

I started in a new position, and decided to use Julia for data management. Stopped immediately. My supervisor did not want to learn a new language.

But I like Julia (apart from the need to scatter flyspecks to get a command to run).

So I kept reading the forums. I am really put off by all the posts about trying to make the syntax more like python. I hate python and I worry that by the time I am in a position to choose to work in Julia it will be so like Python that I won’t want to use it.

I have been on holidays the past couple of weeks.

I had run a series of Negative models in R (glmmTBB) just before Christmas. Negative binomial was the best fit.

I decided to tey to replicate the analysis in Julia.

NOT POSSIBLE. (Or if it is, I cannot work out how, or I am not prepared to put in the work).

While I ultimately wanted group as a random effect, I first tried a fixed effect. But the estimates did not match, even after I sorted out how to include an offset in the mode.

Turns out glm in Julia with NegativeBinomial distribution does not actually fit a negative binomial model. It takes an assumed dispersion parameter. To estimate the dispersion parameter, you are supposed to use the negbin function. But, try as I may, I could not work out how to provide the offset to the negbin function.

But I actually wanted the group to be a random effect.
So I tried fitting in MixedModels.jl. And I got a warning that I should not try fitting models with a dispersion parameter in MixedModels. So Julia itself is telling me it is not suitable for my work (that’s why I love this language).

And that was before I tried to work out how to fit a spline to a continuous covariate. (Turns out doing so is too complicated for my simple brain.)

Conclusion

1 I love the current Julia syntax
2 I hate all the posts about trying to turn Julia into Python
3 Julia is useless to me until I can easily fit generalised linear multilevel models including offsets.

I believe my models are possible in Turing, but my supervisor was very much against any Bayesian approach .And I don’t think I should be forced to a Bayesian analysis just because Julia is not capable of the chosen frequentist analysis.

#2
#3

it’s one thing that the ecosystem is incomplete or lacking (I have a feeling what you want is possible just not well out-of-box), it’s another thing if the language itself is bad for the job etc.

not saying you should do anything to solve this, but many packages in Julia were written because users from a different language saw an opportunity to fill the gap in Julia – as you’ve noticed, an awesome tool for doing the job in many aspects.

for the immediate question, I suggest you to produce a minimal working example in R showing what you’re doing and a few domain-specific Julia packages you’ve tried that should but doesn’t work, the experts from those packages should be able to answer it one way or the other

#4

I have been reading this forum since Covid started. I have looked at all the guidelines.

I tried to do everything right when writing my post. Perhaps you could explain exactly what was wrong. as yorur post is very generic.

#5

The very first item in Make it easier to help you the list reads

  1. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of something generic like “help me with plotting” . Make use of tags and categories (you can edit those and the title after posting ). These will make it more likely that your question will be seen by the people who can potentially help.

I think that a title such as:

glm with NegativeBinomial distribution fails to to fit a negative binomial model.

would stand a better chance of finding qualitifed participants for this discussion. The more so that your post is not short and reading through a half of it still does not give a clue.

2 Likes
#6

If you want to post a specific technical question like, “How do I fit generalised linear multilevel models including offsets?”, that is likely to lead to constructive discussions.

But “Julia is useless to me until there is a Julia package X with feature Y” is not likely to lead to useful discussions. No language has every possible software package freely available, and while it’s perfectly reasonable to choose a language based on some package that is personally important to you, it’s not really other people’s obligation to tailor Julia packages for your needs.

6 Likes
#7

Also just to clearly say there is not concerted (or any, really, as far as I’m aware) effort to “turn Julia into Python”.

You might want to look at

Julia syntax has been very stable since 1.0 (as one might guess from the fact that 1.0 was released and Julia follows SemVer), so if you like the current syntax, there is no reason not to use Julia out of a perceived fear that the language might force you to write Python syntax in future.

1 Like
#8

Could you show what code you tried to pass offsets to negbin? From a quick look at the code, it seems like passing offset like to glm should work, though it’s undocumented (which should be fixed if that’s the case).

1 Like
#9

Maybe @palday or @dmbates can comment on that warning (introduced by https://github.com/JuliaStats/MixedModels.jl/pull/373). There’s probably a good statistical reason for it.

#10

There seems to be a lot going on here.

  • MixedModels.jl does support GLMMs with offsets (added in March 2021)
  • MixedModels.jl does not support negative binomial models and, as of right now, we don’t have any plans to add in that support. Negative binomial is only in the exponential family under certain restrictions and it’s challenging to get that to work well with mixed models. (See also the rather brittle support for negative binomial in R’s lme4, which is the spiritual predecessor to MixedModels.jl. R’s glmmTMB is a little bit better on this front because it uses a more general but less optimized approach; I’ve thought about implementing something similar in Julia as an additional package, but I have a lot of other stuff to do and only so many hours to work on FOSS.) More specifically, the negative binomial distribution is only in the exponential family if you fix the first parameter (the number of failures), which is why GLM.jl requires you to estimate the dispersion separately – you can determine the number of failures from the the location parameter (which is what the GLM framework – not just in GLM.jl but the the general sense used in McCullagh and Nelder’s Generalized Linear Models book – estimates) and dispersion.
  • MixedModels.jl will return a result for models with a dispersion parameter, but it’s likely wrong. There are a couple of thorny issues here that I haven’t had the time to address because it’s a mixture of numerical issues (it seems that we’re failing to find the optimum) and theoretical issues (the approximation we’re using for the log-likelihood doesn’t obviously take the dispersion parameter into account, although several of the quantities evaluated covary with the value of the dispersion parameter), so for a given parameter value, we may be not returning the correct log-likelihood. I’m pretty sure it’s a mixture of these issues and not just one or the other because I’ve also used a very slow but I’m-pretty-sure-correct procedure for evaluating the log likelihood and optimization still failed and this value differed from the optimized procedure. Realistically, I need to spend a few dedicated days working through how the profiled log-likelihood we use for linear mixed model combines with the optimized Laplace approximation we use for GLMM via the local weighted linear approximation to see where the missing term is. Then I need to take the work I’ve been doing to get GLMM opimization to have better starting values (https://github.com/JuliaStats/MixedModels.jl/pull/588) and combine all these things together. Both the theory and practice sections here are things that I really need days blocked off to focus on, but MixedModels.jl is something I do after I get off of work. All that said, for a given set of parameter values, MixedModels.jl will compute the correct predictions, residuals, etc.; but the log-likelihood is probably wrong (and so the optimizer doesn’t find the correct set of parameter values).
  • I imagine you can use Turing.jl and Soss.jl to fit a frequentist model, although neither seem to make it easy out of the box. The trick is that sampling methods like HMC require evaluating the gradient, so if you use the gradient to find the optimum with flat priors, then you have found the maximum likelihood. In other words, you can use a Bayesian PPL and take advantage of the fact that the maximum a posteriori (MAP) estimates under flat priors are equivalent to the maximum likelihood estimates (MLE).
  • Alternatively, you can often reformulate (the inferential question associated with) a negative binomial model as (an inferential question associated with) a binomial or Poisson model, both of which are supported by MixedModels.jl.
  • Splines2.jl and BSplines.jl will help you splines. Creating a package that integrates one of these and StatsModels.jl (which provides the formula interface in Julia) is on my long todo list, because that would help me implement GAMMs, but that’s not happening in the foreseeable future.

And going back to other comments about the OP’s post – the concerns about Julia turning into Python and other general complaints about transitioning to Julia are certainly things that are interesting, but it helps to have a very focused question, especially ones with sample data and minimum working examples (“here’s where I got stuck”). There seem to be two distinct questions in this post:

  1. How do I fit negative binomial models with an offset using GLM.jl?
  2. Why doesn’t MixedModels.jl support models with a dispersion parameter and is this support planned for the future?

Both of these are great questions, but they weren’t immediately obvious in the OP text.

6 Likes