MixedModelsMakie.jl API

Coefficient Plots

MixedModelsMakie.coefplotFunction
coefplot(x::Union{MixedModel,MixedModelBootstrap}; kwargs...)::Figure
coefplot!(fig::Union{GridLayoutBase.GridLayout, GridLayoutBase.GridPosition, GridLayoutBase.GridSubposition, Makie.Figure}, x::Union{MixedModel,MixedModelBootstrap}; 
          kwargs...)
coefplot!(ax::Axis, Union{MixedModel,MixedModelBootstrap};
          conf_level=0.95, vline_at_zero=true, show_intercept=true, attributes...)

Create a coefficient plot of the fixed-effects and associated confidence intervals.

Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) are excluded.

attributes are passed onto scatter! and errorbars!.

The mutating methods return the original object.

Note

Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) are excluded.

source
using CairoMakie
using MixedModels
using MixedModelsMakie
using Random

verbagg = MixedModels.dataset(:verbagg)

gm1 = fit(MixedModel,
          @formula(r2 ~ 1 + anger + gender + btype + situ + (1|subj) + (1|item)),
          verbagg,
          Bernoulli();
          progress=false)

coefplot(gm1)
sleepstudy = MixedModels.dataset(:sleepstudy)

fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy; progress=false)
boot = parametricbootstrap(MersenneTwister(42), 1000, fm1)

coefplot(boot; conf_level=0.999, title="Custom Title")

Ridge Plots

MixedModelsMakie.ridgeplotFunction
ridgeplot(x::Union{MixedModel,MixedModelBootstrap}; kwargs...)::Figure
ridgeplot!(fig::Union{GridLayoutBase.GridLayout, GridLayoutBase.GridPosition, GridLayoutBase.GridSubposition, Makie.Figure}, x::Union{MixedModel,MixedModelBootstrap}; 
          kwargs...)
ridgeplot!(ax::Axis, Union{MixedModel,MixedModelBootstrap};
          conf_level=0.95, vline_at_zero=true, show_intercept=true, attributes...)

Create a ridge plot for the bootstrap samples of the fixed effects.

Densities are normalized so that the maximum density is always 1.

The highest density interval corresponding to conf_level is marked with a bar at the bottom of each density. Setting conf_level=missing removes the markings for the highest density interval.

attributes are passed onto coefplot, band! and lines!.

The mutating methods return the original object.

Note

Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) are excluded.

source
ridgeplot(boot)

Ridge 2D Plots

MixedModelsMakie.ridge2dFunction
ridge2d(bs::MixedModelBootstrap; ptype=:β, kwargs...)
ridge2d!(f::Union{GridLayoutBase.GridLayout, GridLayoutBase.GridPosition, GridLayoutBase.GridSubposition, Makie.Figure}, bs::MixedModelBootstrap; ptype=:β)

Plot pairwise bivariate scatter plots with overlain densities for a bootstrap sample.

ptype specifies the set of parameters to examine, e.g. , , .

The mutating methods return the original object.

source
ridge2d(boot)

Random effects and group-level predictions

Caterpillar Plots

MixedModelsMakie.RanefInfoType
RanefInfo

Information on random effects conditional modes/means, variances, etc.

Used for creating caterpillar plots.

Note

This functionality may be moved upstream into MixedModels.jl in the near future.

source
MixedModelsMakie.ranefinfoFunction
ranefinfo(m::MixedModel)

Return a NamedTuple{fnames(m), NTuple(k, RanefInfo)} from model m

source
ranefinfo(m::MixedModel, gf::Symbol)

Return a RanefInfo corresponding to the grouping variable gf model m.

source
MixedModelsMakie.caterpillarFunction
caterpillar(m::MixedModel, gf::Symbol; kwargs...)

Returns a Figure of a "caterpillar plot" of the random-effects means and prediction intervals

A "caterpillar plot" is a horizontal error-bar plot of conditional means and standard deviations of the random effects.

gf specifies which grouping variable is displayed.

kwargs... are passed on to caterpillar!.

source
MixedModelsMakie.caterpillar!Function
caterpillar!(f::Union{Makie.FigureLike,Makie.GridLayout}, r::RanefInfo;
            orderby=1, cols::Union{Nothing,AbstractVector}=nothing,
            dotcolor=(:red, 0.2), barcolor=:black,
            vline_at_zero::Bool=false)
caterpillar!(f::Union{Makie.FigureLike,Makie.GridLayout}, m::MixedModel,
             gf::Symbol=first(fnames(m)); orderby=1,
             cols::Union{Nothing,AbstractVector}=nothing,
             dotcolor=(:red, 0.2), barcolor=:black,
             vline_at_zero::Bool=false)

Add Axes of a caterpillar plot from r to f.

When passing a MixedModel, gf specifies which grouping variable is displayed. Alternatively, ranefinfo may be used to construct the RanefInfo object directly. Constructing RanefInfo directly can be used to avoid re-computing the conditional variances.

The order of the levels on the vertical axes is increasing orderby column of r.ranef, usually the (Intercept) random effects. Setting orderby=nothing will disable sorting, i.e. return the levels in the order they are stored in.

The display can be restricted to a subset of random effects associated with a grouping variable by specifying cols, either by indices or term names.

Note

Even when not sorting the levels, they might have already been sorted during model matrix construction. If you want impose a particular ordering on the levels, then you must sort the relevant fields in the RanefInfo object before calling caterpillar!.

Note

orderby is the $n$th column of the columns specified by cols.

source
using CairoMakie
CairoMakie.activate!(type = "svg")
using MixedModels
using MixedModelsMakie
sleepstudy = MixedModels.dataset(:sleepstudy)
verbagg = MixedModels.dataset(:verbagg)

fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy; progress=false)
gm0 = fit(MixedModel,
          @formula(r2 ~ 1 + anger + gender + btype + situ + (1|subj) + (1|item)),
          verbagg,
          Bernoulli();
          progress=false)

subjre = ranefinfo(fm1)[:subj]

caterpillar!(Figure(; size=(800,600)), subjre)
caterpillar!(Figure(; size=(800,600)), subjre; orderby=2)
caterpillar!(Figure(; size=(800,600)), subjre; orderby=nothing)
caterpillar(gm0, :item)
MixedModelsMakie.qqcaterpillarFunction
qqcaterpillar(m::MixedModel, gf::Symbol=first(fnames(m));
              kwargs...)

Returns a Figure of a "qq-caterpillar plot" of the random-effects means and prediction intervals.

gf specifies which grouping variable is displayed.

kwargs... are passed on to qqcaterpillar!.

source
MixedModelsMakie.qqcaterpillar!Function
qqcaterpillar!(f::Union{Makie.FigureLike,Makie.GridLayout}, r::RanefInfo;
               cols::Union{Nothing,AbstractVector}=nothing,
               dotcolor=(:red, 0.2), barcolor=:black)
qqcaterpillar!(f::Union{Makie.FigureLike,Makie.GridLayout}, m::MixedModel,
               gf::Symbol=first(fnames(m));
               cols::Union{Nothing,AbstractVector}=nothing,
               dotcolor=(:red, 0.2), barcolor=:black,
               vline_at_zero::Bool=false)

Update the figure with a caterpillar plot with the vertical axis on the Normal() quantile scale.

When passing a MixedModel, gf specifies which grouping variable is displayed. Alternatively, ranefinfo may be used to construct the RanefInfo object directly. Constructing RanefInfo directly can be used to avoid re-computing the conditional variances.

The display can be restricted to a subset of random effects associated with a grouping variable by specifying cols, either by indices or term names.

The order of the levels on the vertical axes is increasing orderby column of r.ranef, usually the (Intercept) random effects. Setting orderby=nothing will disable sorting, i.e. return the levels in the order they are stored in.

source
qqcaterpillar(fm1)
qqcaterpillar(gm0, :item)
qqcaterpillar!(Figure(; size=(400,300)), subjre; cols=[1])
qqcaterpillar!(Figure(; size=(400,300)), subjre; cols=[:days])

Shrinkage Plots

MixedModelsMakie.shrinkageplotFunction
shrinkageplot(m::MixedModel, gf::Symbol=first(fnames(m)), θref, args...; kwargs...)

Return a scatter-plot matrix of the conditional means, b, of the random effects for grouping factor gf.

Two sets of conditional means are plotted: those at the estimated parameter values and those at θref. The default θref results in Λ being a very large multiple of the identity. The corresponding conditional means can be regarded as unpenalized.

args... and kwargs... are passed on to shrinkageplot!

source
MixedModelsMakie.shrinkageplot!Function
shrinkageplot!(f::Union{Makie.FigureLike,Makie.GridLayout}, m::MixedModel,
               gf::Symbol=first(fnames(m)), θref;
               ellipse=false, ellipse_scale=1, n_ellipse=5,
               cols::Union{Nothing,AbstractVector}=nothing,
               shrunk_dotcolor=(:blue, 0.25), ref_dotcolor=(:red, 0.25),
               ellipse_color=:green, ellipse_linestyle=:dash)

Return a scatter-plot matrix of the conditional means, b, of the random effects for grouping factor gf.

Two sets of conditional means are plotted: those at the estimated parameter values and those at θref. The default θref results in Λ being a very large multiple of the identity. The corresponding conditional means can be regarded as unpenalized.

The display can be restricted to a subset of random effects associated with a grouping variable by specifying cols, either by indices or term names.

Correlation ellipses can be added with ellipse=true, with the number of ellipses controlled by n_ellipse. The ellipses are equally spaced between the outer ellipse and the origin (center). The scaling of the ellipses can be adjusted with the multiplicative ellipse_scale. If you are unable to see the ellipses, try increasing ellipse_scale.

Note

For degenerate (singular) models, the correlation ellipse will also be degenerate, i.e., collapse to a point or line.

source
using CairoMakie
using MixedModels
using MixedModelsMakie
sleepstudy = MixedModels.dataset(:sleepstudy)
verbagg = MixedModels.dataset(:verbagg)

fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy; progress=false)
shrinkageplot(fm1)
shrinkageplot(fm1; ellipse=true)
shrinkageplot!(Figure(; size=(400,400)), fm1)
gm1 = fit(MixedModel,
          @formula(r2 ~ 1 + anger + gender + btype + situ + (1|subj) + (1+gender|item)),
          verbagg,
          Bernoulli();
          progress=false)
shrinkageplot(gm1, :item)

Diagnostics

We have also provided a few useful plot recipes for common plot types applied to mixed models. These are especially useful for diagnostics and model checking.

QQ Plots

The methods for qqnorm and qqplot are implemented using Makie recipes. In other words, these are convenience wrappers for calling the relevant plotting methods on residuals(model).

Specify the type of line on the QQ plots with the qqline keyword-argument. The default for qqnorm is :fitrobust, which delivers an R-style line connecting the first and third quartiles. The default for qqplot is :identity, which plots the line with slope = 1 and intercept = 0. The final possibility is :fit, which plots the line of best fit (i.e. regressing the quantiles of the residuals onto the quantiles of the reference distribution).

The reference distribution for qqnorm is the standard normal, which differs from the behavior in previous versions of Makie.

Compat

The options and associated names for the qqline keyword argument changed in Makie 0.16.3 (and were broken in Makie 0.16.0-0.16.2). The equivalent to qqline=:R is qqline=:fitrobust. qqline=:R will be supported for backwards compatibility only until the next breaking release.

using CairoMakie
using MixedModels
using MixedModelsMakie

sleepstudy = MixedModels.dataset(:sleepstudy)

fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy; progress=false)
qqnorm(fm1; qqline=:fitrobust)
# the residuals should have mean 0
# and standard deviation equal to the residual standard deviation
qqplot(Normal(0, fm1.σ), fm1)

Profile Plots

Requires MixedModels 4.14 or above.

using CairoMakie
using MixedModels
using MixedModelsMakie

sleepstudy = MixedModels.dataset(:sleepstudy)
fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy; progress=false)
pr1 = profile(fm1)
zetaplot(pr1)
# show zeta on the absolute value scale with coverage intervals
zetaplot(pr1; absv=true)
# show zeta for the variance components
zetaplot(pr1; absv=true, ptyp='θ')
profiledensity(pr1)
profiledensity(pr1; share_y_scale=false)
profiledensity(pr1; ptyp='σ')

General plots

We also provide a splom or scatter-plot matrix plot for data frames with numeric columns (i.e. a matrix of all pairwise plots). These plots can be used to visualize the joint distribution of, say, the parameter estimates from a simulation.

MixedModelsMakie.splom!Function
splom!(f::Union{Makie.FigureLike,Makie.GridLayout}, df::DataFrame)

Create a scatter-plot matrix in f from the columns of df.

Non-numeric columns are ignored.

source
using CairoMakie
using DataFrames
using LinearAlgebra
using MixedModelsMakie

data = rmul!(randn(100, 3), LowerTriangular([+1 +0 +0;
                                             +1 +1 +0;
                                             -1 -1 +1]))
df = DataFrame(data, [:x, :y, :z])

splom!(Figure(; size=(800, 800)), df)

Meanwhile, splomaxes! provides a lower-level backend for splom!

MixedModelsMakie.splomaxes!Function
splomaxes!(f::Union{Makie.FigureLike,Makie.GridLayout}, labels::AbstractVector{<:AbstractString},
           panel!::Function, args...;
           extraticks::Bool=false, kwargs...)

Populate f with a set of (k*(k-1))/2 axes in a lower triangle for all pairs of labels, where k is the length of labels. The panel! function should have the signature panel!(ax::Axis, i::Integer, j::Integer, args...; kwargs...) and should draw the [i,j] panel in ax.

source
using Statistics

mat = Array(df)
function pfunc(ax, i, j)
    # note that this references mat from the outer scope!
    # [j, i] because:
    # - i is the row in the figure, which corresponds to the y var in each panel
    # - j is the col in the figure, which corresponds to the x var in each panel
    v = view(mat, :, [j, i])
    scatter!(ax, v; color=(:blue, 0.2))
    cc = cor(eachcol(v)...)
    cc = round(cc; digits=2)
    text!(ax, "r=$(cc)")
    return ax
end
splomaxes!(Figure(; size=(800, 800)),
           names(df), pfunc)