# How do I fit generalised linear multilevel models including offsets?

**URL:** <https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286>\
**Category:** Statistics\
**Tags:** glm\
**Created:** [January 9, 2022, 3:03pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286 "2022-01-09T15:03:54Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![timbp](https://avatars.discourse-cdn.com/v4/letter/t/c89c15/32.png) [@timbp](https://discourse.julialang.org/u/timbp)\
**Post date:** [January 9, 2022, 3:03pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/1 "2022-01-09T15:03:54Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [January 9, 2022, 3:09pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/2 "2022-01-09T15:09:49Z")

</div>

> [@Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757):
>
> Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [January 9, 2022, 3:15pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/3 "2022-01-09T15:15:49Z")

</div>

> [@timbp](#):
>
> Turns out `glm` in Julia with `NegativeBinomial` distribution

> [@timbp](#):
>
> So Julia itself is telling me it is not suitable for

> [@timbp](#):
>
> because Julia is not capable of

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

---

<div class="post-metadata">

**Author:** ![timbp](https://avatars.discourse-cdn.com/v4/letter/t/c89c15/32.png) [@timbp](https://discourse.julialang.org/u/timbp)\
**Post date:** [January 9, 2022, 3:23pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/4 "2022-01-09T15:23:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [January 9, 2022, 3:41pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/5 "2022-01-09T15:41:58Z")

</div>

The very first item in [Make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757) 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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 9, 2022, 3:45pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/6 "2022-01-09T15:45:32Z")

</div>

> [@timbp](#):
>
> Perhaps you could explain exactly what was wrong. as yorur post is very generic.

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.

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [January 9, 2022, 3:48pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/7 "2022-01-09T15:48:08Z")

</div>

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

> [@PSA: Julia is not at that stage of development anymore](https://discourse.julialang.org/t/psa-julia-is-not-at-that-stage-of-development-anymore/44872):
>
> You heard about this brand new exciting programming language called Julia and have some suggestions for how it should work. Wonderful! We’re happy to have you here and help you learn the language and listen to your ideas! However, please be aware that “new” is a relative term. Yes, Julia is newer than, say, C, which was first created in 1972 and standardized in 1989. It is not, however, new in the sense that it is still being freely designed and changed in arbitrary ways. Julia development start…

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.

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [January 9, 2022, 5:59pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/8 "2022-01-09T17:59:09Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [January 9, 2022, 6:02pm UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/9 "2022-01-09T18:02:57Z")

</div>

> [@timbp](#):
>
> 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).

Maybe @palday or @dmbates can comment on that warning (introduced by [warning for GLMM with dispersion parameter by palday · Pull Request #373 · JuliaStats/MixedModels.jl · GitHub](https://github.com/JuliaStats/MixedModels.jl/pull/373)). There’s probably a good statistical reason for it.

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [January 10, 2022, 6:00am UTC](https://discourse.julialang.org/t/how-do-i-fit-generalised-linear-multilevel-models-including-offsets/74286/10 "2022-01-10T06:00:16Z")

</div>

There seems to be a lot going on here.

- MixedModels.jl does support GLMMs with offsets (added [in March 2021](https://github.com/JuliaStats/MixedModels.jl/pull/482))
- 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](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](https://github.com/mclements/Splines2.jl) and [BSplines.jl](https://sostock.github.io/BSplines.jl/stable/) 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.