# Is there any glmmTMB package for Julia?

**URL:** <https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868>\
**Category:** Statistics\
**Created:** [January 21, 2019, 12:51pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868 "2019-01-21T12:51:50Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [January 21, 2019, 12:51pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/1 "2019-01-21T12:51:50Z")

</div>

I’m trying to fit mixed-effect models to large dataset.  
R’s lme4 is very slow and can’t deal with more than 50000 rows on my computer. The amount of memory seems to increase with the cube of the number of rows.  
Julia’s Mixedmodels is approximatelly 5 times faster but has the same limit.

There are other R’s packages that sometimes can deal with much larger datasets. For example mgcv but it only works well if you have few levels of the random variable.

The other package I’ve found is glmmTMB (Generalized linear mixed models with various extensions, including zero-inflation. ), I’m trying but for now I’ve got to fit the model with 400000 rows, though it took 850 seconds. I think it will be able to fit larger problems.

Is there any “Template Model Builder” package for Julia? I guess it would be faster than its R counterpart.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 21, 2019, 12:57pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/2 "2019-01-21T12:57:26Z")

</div>

Can you describe the concrete model (the forum supports LaTeX), and your data dimensions (is 400k rows the target, or was that just an experiment)?

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [January 21, 2019, 3:49pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/3 "2019-01-21T15:49:21Z")

</div>

When I wrote he message I got to fit the model with 400000 rows, I’m trying with increasingly more rows and different packages, now I’ve tried again and got to do it with 1.5 million rows (needed 1 hour computation on a computer with 16GB RAM) with glmmTMB. If I try with lme4 and Mixedmodels.jl they just crash even with much smaller problems.

My whole problem has information (from hospitals) of 4 million people through 15 years, repeated measures.  
When I reshape the data to long format I get a table with around 60 million row (if I just keep one measure per year).

I want to fit several models, some linear, some logistic, such as  
(binaryoutcome ~ age + weight + year + sex + treatment1 + treatment2 + race + (1|ID) + (1|hospital), family=bernoulli)  
or maybe a little bit more complex and with some interactions or even nesting the random effects.

I’ve found the main problem is the (1|ID) term because ID has too many levels.

I know I can try different tricks to simplify things but I’d also like to get the results of the whole model in order to compare it.  
And of course I have already standarized the x variables, otherwise it never converges.

---

<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 21, 2019, 4:55pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/4 "2019-01-21T16:55:46Z")

</div>

Sorry if this is slightly off-topic, but given that you’ve been asking quite a few similar questions over the past weeks - have you considered running your models on random samples?

Looking e.g. at the NHS England resource allocation formula, which uses data on 55 million people over 3-5 years, and fits similar models (albeit often with lots more covariates) they seem to get by perfectly fine on a small random sample - you might want to think about just sampling as much as your hardware can handle in a reasonable amount of time?

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [January 21, 2019, 6:02pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/5 "2019-01-21T18:02:18Z")

</div>

Yes, this is one of the things I want to do, and I want to compare the results to the results obtained by using the whole dataset.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 21, 2019, 6:02pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/6 "2019-01-21T18:02:20Z")

</div>

Written in this form, I think your main problem is the very weak or non-existent identification of `ID` (given that the majority of people die at most once).

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [January 21, 2019, 6:07pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/7 "2019-01-21T18:07:43Z")

</div>

Have you tried OLS? The assymptotic properties of the t-statistics are broadly similar. Plus you get the easy benefit of interpretable coefficients.

You should try `fixedeffectmodels.jl` which claims to have something similar to Stata’s `reghdfe`.

---

<div class="post-metadata">

**Author:** ![jkbest2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jkbest2/32/7350_2.png) [@jkbest2](https://discourse.julialang.org/u/jkbest2)\
**Post date:** [January 21, 2019, 8:09pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/8 "2019-01-21T20:09:41Z")

</div>

To address the original question, there isn’t anything equivalent to TMB in Julia yet. For those who aren’t familiar, TMB ([https://arxiv.org/pdf/1509.00660.pdf](https://arxiv.org/pdf/1509.00660.pdf); [GitHub - kaskr/adcomp: AD computation with Template Model Builder (TMB)](https://github.com/kaskr/adcomp)) is a software package that fits statistical models, typically using maximum likelihood (MCMC is also available). You write your model as a templated C++ function that returns a log-likelihood (log-posterior density if you’re Bayesian). in that sense, it is similar in spirit to @Tamas_Papp’s `DynamicHMC.jl`. TMB uses the CppAD package to take derivatives and can integrate out arbitrary sets of parameters (typically any random effects) using the Laplace approximation. The primary interface is through R. This function is then maximized, typically using `R`’s `optim` or similar. `glmmTMB` defines a few common models.

We have many of the pieces in place to do something similar in Julia. `Distributions.jl` defines most of the probability density functions you could ever want, `Optim.jl` has a robust set of optimization algorithms, including those that can use gradient and Hessian information. I’m not quite sure what the higher-order autodiff story is at this point. TMB uses up to third derivatives to get gradients of the Laplace-approximated expected densities. I think the biggest missing piece at this point is a package for doing Laplace approximations.

---

<div class="post-metadata">

**Author:** ![ohadle](https://avatars.discourse-cdn.com/v4/letter/o/7cd45c/32.png) [@ohadle](https://discourse.julialang.org/u/ohadle)\
**Post date:** [January 21, 2019, 9:49pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/9 "2019-01-21T21:49:17Z")

</div>

Going further in the direction of subsampling, I’d suggest plotting a “learning curve”, or some model convergence metric vs. the number of samples in a dataset you can easily fit. This would address the question “how much data does my model need?”, and whether more data is expected to help.  
As someone who works with large medical data, I’d be surprised if your current numbers aren’t more than enough for the model you described, or even much, much more complex ones.

Though that TMB idea does sound pretty cool.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [January 21, 2019, 11:15pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/10 "2019-01-21T23:15:59Z")

</div>

I was doing that MixedModels.jl and lme4 can’t deal with that amount of data.

Sumbsampling sounds nice but sometimes it doesn’t work: when some combinations of the predictive variables happen very rarely.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [January 21, 2019, 11:19pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/11 "2019-01-21T23:19:12Z")

</div>

All right, I told death but I meant any binary outcome.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 22, 2019, 6:37am UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/12 "2019-01-22T06:37:06Z")

</div>

> [@jkbest2](#):
>
> I’m not quite sure what the higher-order autodiff story is at this point.

From the paper it looks like they use CppAD _and_ optimize for sparsity in the AD graph, which is a neat trick, but I am not aware of any explicit support for this in Julia’s AD packages.

I am curious about the accuracy of Laplace approximations, it is a pity that they only compare to another package using the same methodology.

---

<div class="post-metadata">

**Author:** ![jkbest2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jkbest2/32/7350_2.png) [@jkbest2](https://discourse.julialang.org/u/jkbest2)\
**Post date:** [January 22, 2019, 3:25pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/13 "2019-01-22T15:25:35Z")

</div>

Yeah, I think the sparsity optimization is key for performance.

There is a comparison of random effects estimated via the Laplace approximation vs. MCMC [in this paper](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0197954).

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [January 22, 2019, 3:42pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/14 "2019-01-22T15:42:15Z")

</div>

Will this be integrated on MixedModels.jl?

---

<div class="post-metadata">

**Author:** ![jkbest2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jkbest2/32/7350_2.png) [@jkbest2](https://discourse.julialang.org/u/jkbest2)\
**Post date:** [January 22, 2019, 5:27pm UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/15 "2019-01-22T17:27:40Z")

</div>

That’s probably unlikely; these techniques are general and would fit better in a dedicated package. I also don’t think that `MixedModels.jl` depends on any autodiff packages, so this would be a large change with limited payoff.

I think that you should reconsider subsampling your data. Remember that if you have enough covariates, every covariate combination is rare. You want to think carefully about _how_ you subsample; don’t take random rows, but subsample `ID`, then take the rows that correspond to those `ID`s. That way you keep entire patient histories. It’s unlikely that you’ll lose much.

---

<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:** [April 28, 2022, 4:44am UTC](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/16 "2022-04-28T04:44:31Z")

</div>

It’s a long time since, but dealing with blocking variables with many levels is one of the reasons we introduced the `Grouping()` pseudo-constrast. So you use e.g. `contrasts=Dict(:ID => Grouping())` and this avoids constructing a real contrast matrix (which is what causes the memory problem) that’s never actually used when constructing the random effects matrices.