Additional Functionality in Other Packages

Several packages extend the functionality of MixedModels.jl, both in ways specific to mixed models and in ways applicable to more general regression models. In the following, we will use the models from the previous sections to showcase this functionality.

using MixedModels
progress = isinteractive()

false

insteval = MixedModels.dataset("insteval")
ie1 = fit(MixedModel,
          @formula(y ~ 1 + studage + lectage + service + (1|s) + (1|d) + (1|dept)),
          insteval; progress)

ie2 = fit(MixedModel,
          @formula(y ~ 1 + studage + lectage + service +
                      (1 | s) +
                      (1 + service | d) +
                      (1 + service | dept)),
          insteval; progress)

sleepstudy = MixedModels.dataset("sleepstudy")
ss1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), sleepstudy; progress)

ss2 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy; progress)

using DataFrames
contra = DataFrame(MixedModels.dataset("contra"))
contra[!, :anych] .= contra[!, :livch] .!= "0"
contrasts = Dict(:livch => EffectsCoding(; base="0"),
                 :urban => HelmertCoding(),
                 :anych => HelmertCoding())
gm1 = fit(MixedModel,
          @formula(use ~ 1 + urban + anych * age + abs2(age) + (1 | dist & urban)),
          contra,
          Bernoulli();
          contrasts,
          progress)

MixedModelsExtras.jl

https://palday.github.io/MixedModelsExtras.jl/v2

MixedModelsExtras.jl is a collection of odds-and-ends that may be useful when working with mixed effects models, but which we do not want to include in MixedModels.jl at this time. Some functions may one day migrate to MixedModels.jl, when we are happy with their performance and interface (e.g. vif), but some are intentionally omitted from MixedModels.jl (e.g. r2, adjr2).

using MixedModelsExtras

r2(ss2; conditional=true)

0.8263130612096021

r2(ss2; conditional=false)

0.28647138741126765

icc(ie2)

0.2885292866744761

icc(ie2, :dept)

0.016119194950382547

vif(ie1)

9-element Vector{Float64}:
 1.5141903374729053
 1.7354060220617624
 1.7822316986634812
 1.449378975046207
 1.4380891514780725
 1.5948966178188415
 1.4634020911560608
 1.8267103207894357
 1.0161785415795814

DataFrame(; coef=fixefnames(ie1)[2:end], VIF=vif(ie1))

gvif(ie1)

3-element Vector{Float64}:
 1.3110872223511254
 1.325731162491792
 1.0161785415795814

DataFrame(; term=termnames(ie1)[2][2:end], GVIF=gvif(ie1))

RegressionFormulae.jl

https://github.com/kleinschmidt/RegressionFormulae.jl

RegressionFormulae.jl provides a few extensions to the somewhat more restricted variant of the Wilkinson-Roger notation found in Julia. In particular, it adds / for nested designs within the fixed effects and ^ for computing interactions only up to a certain order.

using RegressionFormulae

fit(MixedModel,
          @formula(y ~ 1 + service / (studage + lectage) +
                      (1 | s) +
                      (1 | d) +
                      (1 | dept)),
          insteval; progress)

fit(MixedModel,
          @formula(y ~ 1 + (studage + lectage + service)^2 +
                      (1 | s) +
                      (1 | d) +
                      (1 | dept)),
          insteval; progress)

BoxCox.jl

https://palday.github.io/BoxCox.jl/v0.3/

BoxCox.jl implements a the Box-Cox transformation in an efficient way. Via package extensions, it supports specializations for MixedModels.jl and several plotting functions, but does not incur a dependency penalty for this functionality when MixedModels.jl or Makie.jl are not loaded.

using BoxCox

bc = fit(BoxCoxTransformation, ss2)

Box-Cox transformation

estimated λ: -1.0749
resultant transformation:

 y^-1.1 - 1
------------
    -1.1

using CairoMakie
boxcoxplot(bc; conf_level=0.95)

The estimated λ is very close to -1, i.e. the reciprocal of reaction time, which has a natural interpretation as speed. In other words, the Box-Cox transformation suggests that we should consider modelling the sleepstudy data as speed (reaction per unit time) instead of reaction time:

fit(MixedModel, @formula(1000 / reaction ~ 1 + days + (1 + days|subj)), sleepstudy)

(We multiply by 1000 to get the responses per second instead of the responses per millisecond.)

Effects.jl

https://beacon-biosignals.github.io/Effects.jl/v1.2/

Effects.jl provides a convenient method to compute effects, i.e. predictions and associated prediction intervals computed at points on a reference grid. For models with a nonlinear link function, Effects.jl will also compute appropriate errors on the response scale based on the difference method.

For MixedModels.jl, the predictions are computed based on the fixed effects only.

The functionality of Effects.jl was inspired by the effects and emmeans packages in R and the methods within are based on @fox:effect:2003.

using Effects

design = Dict(:age => -15:1:20,
              :anych => [true, false])

eff_logit = effects(design, gm1; eff_col="use", level=0.95)

first(eff_logit, 10)

eff_prob = effects(design, gm1; eff_col="use", level=0.95, invlink=AutoInvLink())

first(eff_prob, 10)

Effects are particularly nice for visualizing the model fit and its predictions.

using AlgebraOfGraphics # like ggplot2, but an algebra instead of a grammar
using CairoMakie

plt1 = data(eff_logit) * mapping(:age; color=:anych) *
      (mapping(:use) * visual(Lines) +
       mapping(:lower, :upper) * visual(Band; alpha=0.3))
draw(plt1)

plt2 = data(eff_prob) * mapping(:age; color=:anych) *
      (mapping(:use) * visual(Lines) +
       mapping(:lower, :upper) * visual(Band; alpha=0.3))
draw(plt2)

using Statistics: mean
contra_by_age = transform(contra,
                          :age => ByRow(x -> round(Int, x)),
                          :use => ByRow(==("Y"));
                          renamecols=false)
contra_by_age = combine(groupby(contra_by_age, [:age, :anych]),
                        :use => mean => :use)
plt3 = plt2 +
       data(contra_by_age) *
       mapping(:age, :use;
               color=:anych => "children") * visual(Scatter)

draw(plt3;
     axis=(; title="Estimated contraceptive use by age and children",
            limits=(nothing, (0, 1)) # ylim=0,1, xlim=auto
            ))

Effects and estimated marginal (least squares) means are closely related and partially concepts. Effects.jl provides convenience function emmeans and empairs for computing EM means and pairwise differences of EM means.

emmeans(gm1)

empairs(gm1; dof=Inf)

Margins.jl

StandardizedPredictors.jl

https://beacon-biosignals.github.io/StandardizedPredictors.jl/v1/

StandardizedPredictors.jl provides a convenient way to express centering, scaling, and z-standardization as a "contrast" via the pseudo-contrasts Center, Scale, ZScore. Because these use the usual contrast machinery, they work well with any packages that use that machinery correctly (e.g. Effects.jl). The default behavior is to empirically compute the center and scale, but these can also be explicitly provided, either as a number or as a function (e.g. median to use the median for centering.)

using StandardizedPredictors

contrasts = Dict(:days => Center())
fit(MixedModel,
    @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy;
    contrasts)

MixedModelsSmallSample.jl

https://arnostrouwen.github.io/MixedModelsSmallSample.jl/stable/

MixedModelsSmallSample.jl provides the Satterthwaite and Kenward-Roger approximations for denominator degrees of freedom.

using MixedModelsSmallSample

sw = small_sample_adjust(ss1, Satterthwaite())
coeftable(ss1)

──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405     9.50618   26.45    <1e-99
days          10.4673    0.801735  13.06    <1e-38
──────────────────────────────────────────────────

coeftable(sw)

────────────────────────────────────────────────────────────
                Coef.  Std. Error     DenDF      t  Pr(>|t|)
────────────────────────────────────────────────────────────
(Intercept)  251.405     9.50618    24.4905  26.45    <1e-18
days          10.4673    0.801735  162.0     13.06    <1e-26
────────────────────────────────────────────────────────────

# Kenward-Roger degrees of freedom requires REML fit

ss1_reml = refit!(ss1; REML=true)
kr = small_sample_adjust(ss1_reml, KenwardRoger())

coeftable(ss1_reml)

──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405     9.74672   25.79    <1e-99
days          10.4673    0.804221  13.02    <1e-38
──────────────────────────────────────────────────

coeftable(kr)

────────────────────────────────────────────────────────────
                Coef.  Std. Error     DenDF      t  Pr(>|t|)
────────────────────────────────────────────────────────────
(Intercept)  251.405     9.74672    22.8102  25.79    <1e-17
days          10.4673    0.804221  161.0     13.02    <1e-26
────────────────────────────────────────────────────────────

RCall.jl and JellyMe4.jl

https://juliainterop.github.io/RCall.jl/stable/

https://github.com/palday/JellyMe4.jl/

RCall.jl provides a convenient interface for interoperability with R from Julia. JellyMe4.jl extends the functionality of RCall so that MixedModels.jl-fitted models and lme4-fitted models can be translated to each other. In practical terms, this means that you can enjoy the speed of Julia for model fitting, but use all the extra packages you love from R's larger ecosystem.

MixedModelsSerialization.jl

https://juliamixedmodels.github.io/MixedModelsSerialization.jl/stable/api/

MixedModelsSerialization.jl provides a reduced-memory "summary" representation of a fitted model, LinearMixedModelSummary, that discards the model matrices and other data-sized fields while retaining the fixed- and random-effects estimates, θ, the log-likelihood, and other quantities needed to support many StatsAPI and MixedModels methods (e.g. coef, vcov, VarCorr, coeftable). This makes it practical to save and later reload models fit to very large datasets without also serializing the original data. The package is a proving ground for these ideas so that its API can evolve and have breaking releases independently of MixedModels.jl.

MixedModelsSim.jl

https://repsychling.github.io/MixedModelsSim.jl/stable/

MixedModelsSim.jl provides utilities for generating experimental designs, especially designs with crossed grouping factors (e.g. "Subject" and "Item") and both within- and between-unit experimental factors. Combined with simulate/parametricbootstrap from MixedModels.jl, it is commonly used for power analysis: create a design, specify hypothesized effect sizes and variance components, simulate many datasets, and examine the distribution of the resulting test statistics.