Gradient and Hessian computation

Experimental support for computing the gradient and the Hessian of the objective function (i.e. negative twice the profiled log likelihood) via ForwardDiff.jl and FiniteDiff.jl are provided as package extensions.

via ForwardDiff.jl

The core functionality is provided by defining appropriate methods for ForwardDiff.gradient and ForwardDiff.hessian:

ForwardDiff.gradient — Method
ForwardDiff.gradient(model::LinearMixedModel)

Evaluate the gradient of the objective function at the currently fitted parameter values.

Large allocations

Most of MixedModels.jl relies strongly on in-place methods in order to minimize the amount of memory allocated. In addition to reducing the memory burden (especially for large models), this practice generally speeds up evaluation of the objective. In-place methods, however, generally do not play well with automatic differentiation. For the automatic differentiation support provided here, the developers instead implemented alternative, out-of-place methods. These will generally be slower and much more memory intensive, so use of this functionality is not recommended for large models.

ForwardDiff.jl support is experimental.

Compatibility with ForwardDiff.jl is experimental. The precise structure, including function names and method definitions, is subject to change without being considered a breaking change. In particular, the exact set of parameters included is subject to change. The θ parameter is always included, but whether σ and/or the fixed effects should be included is currently still being decided.

ForwardDiff.hessian — Method
ForwardDiff.hessian(model::LinearMixedModel)

Evaluate the Hessian of the objective function at the currently fitted parameter values.

Large allocations

Most of MixedModels.jl relies strongly on in-place methods in order to minimize the amount of memory allocated. In addition to reducing the memory burden (especially for large models), this practice generally speeds up evaluation of the objective. In-place methods, however, generally do not play well with automatic differentiation. For the automatic differentiation support provided here, the developers instead implemented alternative, out-of-place methods. These will generally be slower and much more memory intensive, so use of this functionality is not recommended for large models.

ForwardDiff.jl support is experimental.

Compatibility with ForwardDiff.jl is experimental. The precise structure, including function names and method definitions, is subject to change without being considered a breaking change. In particular, the exact set of parameters included is subject to change. The θ parameter is always included, but whether σ and/or the fixed effects should be included is currently still being decided.

Exact zero at optimum for trivial models

using MixedModels, MixedModelsDatasets, ForwardDiff
fm1 = lmm(@formula(yield ~ 1 + (1|batch)), MixedModelsDatasets.dataset(:dyestuff2))
Linear mixed model fit by maximum likelihood
 yield ~ 1 + (1 | batch)
   logLik   -2 logLik     AIC       AICc        BIC    
   -81.4365   162.8730   168.8730   169.7961   173.0766

Variance components:
            Column   Variance Std.Dev.
batch    (Intercept)   0.00000 0.00000
Residual              13.34610 3.65323
 Number of obs: 30; levels of grouping factors: 6

  Fixed-effects parameters:
───────────────────────────────────────────────
              Coef.  Std. Error     z  Pr(>|z|)
───────────────────────────────────────────────
(Intercept)  5.6656    0.666986  8.49    <1e-16
───────────────────────────────────────────────
ForwardDiff.gradient(fm1)
1-element Vector{Float64}:
 0.0
ForwardDiff.hessian(fm1)
1×1 Matrix{Float64}:
 28.768680157921427

Approximate zero at optimum for non trivial models

fm2 = lmm(@formula(reaction ~ 1 + days + (1+days|subj)), MixedModelsDatasets.dataset(:sleepstudy))
Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 + days | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -875.9697  1751.9393  1763.9393  1764.4249  1783.0971

Variance components:
            Column    Variance Std.Dev.   Corr.
subj     (Intercept)  565.51080 23.78047
         days          32.68252  5.71686 +0.08
Residual              654.94042 25.59180
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405      6.63226  37.91    <1e-99
days          10.4673     1.50224   6.97    <1e-11
──────────────────────────────────────────────────
ForwardDiff.gradient(fm2)
3-element Vector{Float64}:
 4.4116120307080564e-5
 0.00284014455504078
 0.0017907903212801557
ForwardDiff.hessian(fm2)
3×3 Matrix{Float64}:
 45.4121    35.9342    6.35331
 35.9342   465.731   203.975
  6.35331  203.975   963.918

via FiniteDiff.jl

The core functionality is provided by defining appropriate methods for FiniteDiff.finite_difference_gradient and FiniteDiff.finite_difference_hessian:

FiniteDiff.finite_difference_gradient — Method
FiniteDiff.finite_difference_gradient(model::LinearMixedModel, args...; kwargs...)

Evaluate the gradient of the objective function at the currently fitted parameter values.

FiniteDiff.jl support is experimental.

Compatibility with FiniteDiff.jl is experimental. The precise structure, including function names and method definitions, is subject to change without being considered a breaking change. In particular, the exact set of parameters included is subject to change. The θ parameter is always included, but whether σ and/or the fixed effects should be included is currently still being decided.

FiniteDiff.finite_difference_hessian — Method
FiniteDiff.finite_difference_hessian(model::LinearMixedModel, args...; kwargs...)

Evaluate the Hessian of the objective function at the currently fitted parameter values.

FiniteDiff.jl support is experimental.

Compatibility with FiniteDiff.jl is experimental. The precise structure, including function names and method definitions, is subject to change without being considered a breaking change. In particular, the exact set of parameters included is subject to change. The θ parameter is always included, but whether σ and/or the fixed effects should be included is currently still being decided.

using FiniteDiff
FiniteDiff.finite_difference_gradient(fm2)
3-element Vector{Float64}:
 4.466402795463656e-5
 0.002841769899686239
 0.0017902784269433515
FiniteDiff.finite_difference_hessian(fm2)
3×3 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
  40.8873   31.5147  -14.5312
  31.5147  461.416   183.578
 -14.5312  183.578   867.528