How do we extract the contrast/hypothesis matrix for a predictor from a fitted model?

[yjunechoe]

To see what built-in coding schemes like SeqDiffCoding() got translated into for our data/model

[palday]

Let's take a look at a simple model. For now, I won't specify contrasts explicitly, so they should be dummy coding.

julia> fm1 = fit(MixedModel, @formula(rt_trunc ~ 1 + spkr*prec*load + (1|subj) + (1|item)), MixedModels.dataset(:kb07));

I can extract the formula associated with a fitted model via the formula function. This formula differs slightly from the formula you put into the model because it now carries detailed information about contrasts for categorical values and various other meta data related to applying the data schema to the formula.

julia> f = formula(fm1)
FormulaTerm
Response:
  rt_trunc(continuous)
Predictors:
  1
  spkr(DummyCoding:2→1)
  prec(DummyCoding:2→1)
  load(DummyCoding:2→1)
  spkr(DummyCoding:2→1) & prec(DummyCoding:2→1)
  spkr(DummyCoding:2→1) & load(DummyCoding:2→1)
  prec(DummyCoding:2→1) & load(DummyCoding:2→1)
  spkr(DummyCoding:2→1) & prec(DummyCoding:2→1) & load(DummyCoding:2→1)
  (1 | subj)
  (1 | item)

Let's start taking this apart. First, we look at the right hand side (rhs), i.e. the predictors.

julia> f.rhs
1spkr(DummyCoding:2→1)prec(DummyCoding:2→1)load(DummyCoding:2→1)spkr(DummyCoding:2→1) & prec(DummyCoding:2→1)spkr(DummyCoding:2→1) & load(DummyCoding:2→1)prec(DummyCoding:2→1) & load(DummyCoding:2→1)spkr(DummyCoding:2→1) & prec(DummyCoding:2→1) & load(DummyCoding:2→1)
(1 | subj)
(1 | item)

The first entry in the RHS are the fixed effects

julia> f.rhs[1]
1spkr(DummyCoding:2→1)prec(DummyCoding:2→1)load(DummyCoding:2→1)spkr(DummyCoding:2→1) & prec(DummyCoding:2→1)spkr(DummyCoding:2→1) & load(DummyCoding:2→1)prec(DummyCoding:2→1) & load(DummyCoding:2→1)spkr(DummyCoding:2→1) & prec(DummyCoding:2→1) & load(DummyCoding:2→1)

We can now look at the set of terms that make up the fixed effects.

julia> f.rhs[1].terms
1
spkr(DummyCoding:2→1)
prec(DummyCoding:2→1)
load(DummyCoding:2→1)
spkr(DummyCoding:2→1) & prec(DummyCoding:2→1)
spkr(DummyCoding:2→1) & load(DummyCoding:2→1)
prec(DummyCoding:2→1) & load(DummyCoding:2→1)
spkr(DummyCoding:2→1) & prec(DummyCoding:2→1) & load(DummyCoding:2→1)

Now let's zoom in on just one of those terms:

julia> spkr = f.rhs[1].terms[2]
spkr(DummyCoding:2→1)

Now let's look at the contrasts property of that term:

julia> spkr.contrasts
StatsModels.ContrastsMatrix{DummyCoding, Matrix{Float64}, String, String}([0.0; 1.0;;], ["old"], ["new", "old"], DummyCoding(nothing, nothing), Dict("new" => 1, "old" => 2))

Okay, that's a lot of stuff, but let's zoom in again. This is the coefficient matrix associated with that term. Since there are two levels to spkr, we have two rows and one column in the matrix -- i.e. one contrast.

julia> spkr.contrasts.matrix
2×1 Matrix{Float64}:
 0.0
 1.0

This single contrast corresponds to a single coefficient named old because it's the contrast between old and the reference level (new because dummy coding):

julia> spkr.contrasts.coefnames
1-element Vector{String}:
 "old"

The coefficient names are the column names of the contrast matrix; the levels are the rows:

julia> spkr.contrasts.levels
2-element Vector{String}:
 "new"
 "old"

We can also repeat this whole things for effects coding, which I'll present here with less commentary:

julia> contrasts = Dict(:spkr => EffectsCoding(), :load => EffectsCoding(), :prec => EffectsCoding())
Dict{Symbol, EffectsCoding} with 3 entries:
  :spkr => EffectsCoding(nothing, nothing)
  :load => EffectsCoding(nothing, nothing)
  :prec => EffectsCoding(nothing, nothing)

julia> fm2 = fit(MixedModel, @formula(rt_trunc ~ 1 + spkr*prec*load + (1|subj) + (1|item)), MixedModels.dataset(:kb07); contrasts);

julia> f2 = formula(fm2);

julia> spkr2 = f2.rhs[1].terms[2]
spkr(EffectsCoding:2→1)

We can see that the levels are the same (as expected):

julia> spkr2.contrasts.levels
2-element Vector{String}:
 "new"
 "old"

The coefnames are also the same, but the contrast encoded is different because the reference level differs:

julia> spkr2.contrasts.coefnames
1-element Vector{String}:
 "old"

julia> spkr2.contrasts.matrix
2×1 Matrix{Float64}:
 -1.0
  1.0

(Recall that the reference level here -- the "zero point" of these contrasts doesn't correspond to any individual level but rather the grand mean of conditions.)