# Laplace approximation for mixed models in Julia (TMB-like functionality)

**URL:** <https://discourse.julialang.org/t/laplace-approximation-for-mixed-models-in-julia-tmb-like-functionality/61796>\
**Category:** Statistics\
**Tags:** question\
**Created:** [May 25, 2021, 5:08pm UTC](https://discourse.julialang.org/t/laplace-approximation-for-mixed-models-in-julia-tmb-like-functionality/61796 "2021-05-25T17:08:27Z")\
**Posts on this page:** 1\
**Showing post:** 7

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [May 25, 2021, 7:12pm UTC](https://discourse.julialang.org/t/laplace-approximation-for-mixed-models-in-julia-tmb-like-functionality/61796/7 "2021-05-25T19:12:23Z")

</div>

Welcome Cole! I gave that [talk](https://discourse.julialang.org/t/an-introduction-to-julia-for-quantitative-ecologists/61112), and can answer a few of your questions.

All the pieces mentioned in that previous TMB [thread](https://discourse.julialang.org/t/is-there-any-glmmtmb-package-for-julia/19868/8) are still in place, and we have a few more as well. I’m not sure about the status of higher-order AD, but as to the other pieces:

- Hessian ~~sparsity detection~~ sparse differentiation: [sort of there, but not quite](https://discourse.julialang.org/t/best-simplest-way-to-calculate-sparse-hessian/58629/3). My impression (maybe @ChrisRackauckas can elaborate) is that it wouldn’t take too much effort to implement, it just needs attention from someone with the time and domain knowledge.

- Laplace approximation: given an optimizer and an AD system, this is in principle just three lines of code (assuming `f` is your log-density function and `N` is the number of variables):

```julia
using Optim, ForwardDiff, LinearAlgebra
opt = optimize(f, zeros(N), LBFGS(), autodiff=:forward)
H = ForwardDiff.hessian(f, opt.minimizer)
integral = -opt.minimum + 0.5 * (log((2π)^N) - logdet(H))

```

@jkbest and I already have working implementations [here](https://github.com/jkbest2/LaplaceApproximations.jl) and [here](https://github.com/ElOceanografo/MarginalLogDensities.jl). The trickier parts are making sure the optimization is fast, and the Hessian is automatically sparse if it can be.

- Sparse Cholesky decomposition in Julia is done via SuiteSparse/CHOLMOD, which I think is the same as TMB. This [PR](https://github.com/JuliaLang/julia/pull/40716) should let those decompositions be used more easily for AD and Laplace approximations after the next Julia release (see also [this issue](https://github.com/TuringLang/DistributionsAD.jl/issues/89)).

Basically, I agree with @jkbest that a TMB equivalent in Julia is very much within reach. My work-in-progress [MarginalLogDensities](https://github.com/ElOceanografo/MarginalLogDensities.jl) package can already take an arbitrary log-density function and integrate out an arbitrary subset of its arguments using the Laplace approximation with automatic sparsity detection (though that part is not optimal). See this [simple demo](https://github.com/ElOceanografo/MarginalLogDensities.jl/blob/master/test/MarginalLogDensities.jl%20Demo.ipynb) for an example. It will just take some sustained effort from a couple/few knowledgeable people to fill in the missing bits and plug all the pieces together.

---

_[View the full topic](https://discourse.julialang.org/t/laplace-approximation-for-mixed-models-in-julia-tmb-like-functionality/61796)._