Mixed Models Tutorial: Contrast Coding

Author

Reinhold Kliegl

Published

2026-08-23

After working through this page you will be able to:

apply sequential-difference, Helmert, and hypothesis contrast coding to a multi-level factor;
interpret the fixed-effect parameters produced by each coding scheme;
derive a set of contrasts from a PCA of the data.

Before you start

Prerequisites: Transformations of the predictors and the response.

Datasets used: fggk21 (Emotikon; see the dataset catalog).

This script uses a subset of data reported in Fühner et al. (2021).

All children were between 6.0 and 6.99 years at legal keydate (30 September) of school enrollment, that is in their ninth year of life in the third grade. To avoid delays associated with model fitting we work with a reduced data set and less complex models than those in the reference publication. The script requires only a few changes to specify the more complex models in the paper.

The script is structured in three main sections:

Setup with reading and examining the data
Contrasts coding

Effect and sequential difference contrasts
Helmert contrast
Hypothesis contrast
PCA-based contrast

Other topics

LMM goodness of fit does not depend on contrast (i.e., reparameterization)
VCs and CPs depend on contrast
VCs and CPs depend on random factor

1 Setup

1.1 Packages and functions

Code

using CategoricalArrays
using DataFrames
using MixedModels
using Statistics
using MixedModels: likelihoodratiotest
using SMLP2026: fggk21_teaching_sample

const progress = isinteractive()

1.2 Preprocessing

Variable descriptions for dataset(:fggk21) — including the five Test levels and their units — are in the dataset catalog.

We work with the canonical teaching subsample of the Emotikon data: a stratified sample of 1000 boys and 1000 girls (whole children, keeping all of a child’s test scores), with age centered at 8.5 years (a1) and z-scores of score computed within each Test on the full data (zScore). The preprocessing and sampling steps are walked through in the Emotikon capstone, which analyzes the same subsample, so estimates are directly comparable between the two pages.

dat = fggk21_teaching_sample()

9704×8 DataFrame

9679 rows omitted

Row	Test	Cohort	School	Child	Sex	age	a1	zScore
	Cat…	String	String	String	Cat…	Float64	Float64	Float64
1	Speed	2013	S100432	C002372	Boys	7.99452	-0.505476	2.49861
2	PowerUP	2013	S100432	C002372	Boys	7.99452	-0.505476	1.61952
3	PowerLOW	2013	S100432	C002372	Boys	7.99452	-0.505476	-1.06736
4	Coordination	2013	S100432	C002372	Boys	7.99452	-0.505476	0.572288
5	Endurance	2013	S100432	C002372	Boys	7.99452	-0.505476	0.628055
6	Speed	2013	S100924	C002412	Girls	7.99452	-0.505476	-0.861152
7	PowerUP	2013	S100924	C002412	Girls	7.99452	-0.505476	-0.90311
8	PowerLOW	2013	S100924	C002412	Girls	7.99452	-0.505476	-0.188712
9	Coordination	2013	S100924	C002412	Girls	7.99452	-0.505476	0.190742
10	Endurance	2013	S100924	C002412	Girls	7.99452	-0.505476	-1.18768
11	Speed	2013	S102295	C002539	Girls	7.99452	-0.505476	-1.06685
12	PowerUP	2013	S102295	C002539	Girls	7.99452	-0.505476	-1.32355
13	PowerLOW	2013	S102295	C002539	Girls	7.99452	-0.505476	-2.87634
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
9693	Coordination	2018	S106781	C117692	Boys	9.10609	0.606092	0.538978
9694	Endurance	2018	S106781	C117692	Boys	9.10609	0.606092	0.938553
9695	Speed	2018	S111302	C117786	Girls	9.10609	0.606092	0.31086
9696	PowerUP	2018	S111302	C117786	Girls	9.10609	0.606092	-0.622817
9697	PowerLOW	2018	S111302	C117786	Girls	9.10609	0.606092	1.1551
9698	Coordination	2018	S111302	C117786	Girls	9.10609	0.606092	-1.39492
9699	Endurance	2018	S111302	C117786	Girls	9.10609	0.606092	-0.154938
9700	Speed	2018	S130370	C117936	Boys	9.10609	0.606092	0.31086
9701	PowerUP	2018	S130370	C117936	Boys	9.10609	0.606092	-0.482671
9702	PowerLOW	2018	S130370	C117936	Boys	9.10609	0.606092	0.224769
9703	Coordination	2018	S130370	C117936	Boys	9.10609	0.606092	0.572288
9704	Endurance	2018	S130370	C117936	Boys	9.10609	0.606092	0.452557

The subsample carries descriptive labels for the levels of Test; we recode them back to the shorter actual test names to reduce clutter in the model outputs.

recode!(
  dat.Test,
  "Endurance" => "Run",
  "Coordination" => "Star_r",
  "Speed" => "S20_r",
  "PowerLOW" => "SLJ",
  "PowerUP" => "BPT",
)
combine(
  groupby(dat, [:Test, :Sex]),
  :zScore => mean,
  :zScore => std,
)

10×4 DataFrame

Row	Test	Sex	zScore_mean	zScore_std
	Cat…	Cat…	Float64	Float64
1	Run	Boys	0.28737	1.0491
2	Run	Girls	-0.248045	0.894015
3	Star_r	Boys	0.115898	1.04888
4	Star_r	Girls	-0.113522	0.96379
5	S20_r	Boys	0.137753	1.02252
6	S20_r	Girls	-0.170976	0.950767
7	SLJ	Boys	0.183029	0.957072
8	SLJ	Girls	-0.213059	0.915784
9	BPT	Boys	0.333946	0.942954
10	BPT	Girls	-0.354512	0.89236

The main results of relevance here are shown in Figure 2 of Scientific Reports 11:17566; the figure is constructed in the Emotikon capstone and, as a plotting exercise, in Creating multi-panel plots.

2 Contrast coding

Contrast coding is part of StatsModels.jl. Here is the primary author’s (i.e., Dave Kleinschmidt’s) documentation of Modeling Categorical Data.

A couple of general remarks about the following examples. First, all contrasts defined in this tutorial return an estimate of the Grand Mean (GM) in the intercept, that is they are so-called sum-to-zero contrasts. In both Julia and R the default contrast is dummy coding which is not a sum-to-zero contrast, but returns the mean of the reference (control) group - unfortunately for (quasi-)experimentally minded scientists.

Second, the factor Sex has only two levels. We use effects coding (also known as sum coding in R) to estimate the difference of the levels from the Grand Mean. Unlike in R, the default sign of the effect is for the second level (the base level is the first, not the last level), but this can be changed with the base kwarg in the command. Effects coding is a sum-to-zero contrast, but when applied to factors with more than two levels does not yield orthogonal contrasts.

Finally, contrasts for the five levels of the fixed factor Test represent the hypotheses about differences between them. In this tutorial, we use this factor to illustrate various options.

We (initially) include only Test as fixed factor and Child as random factor. More complex LMMs can be specified by simply adding other fixed or random factors to the formula.

2.1 Sequential difference coding: `contr1`

Sequential difference coding was used in the publication. This specification tests pairwise differences between the five neighboring levels of Test, that is:

SDC1: 2-1
SDC2: 3-2
SDC3: 4-3
SDC4: 5-4

The levels were sorted such that these contrasts map onto four a priori hypotheses; in other words, they are theoretically motivated pairwise comparisons. The motivation also encompasses theoretically motivated interactions with Sex. The order of levels can also be explicitly specified during contrast construction. This is very useful if levels are in a different order in the dataframe. We recommend the explicit specification to increase transparency of the code.

The statistical disadvantage of sequential difference coding is that the contrasts are not orthogonal, that is the contrasts are correlated. This is obvious from the fact that levels 2, 3, and 4 are all used in two contrasts. One consequence of this is that correlation parameters estimated between neighboring contrasts (e.g., 2-1 and 3-2) are difficult to interpret. Usually, they will be negative because assuming some practical limitation on the overall range (e.g., between levels 1 and 3), a small “2-1” effect “correlates” negatively with a larger “3-2” effect for mathematical reasons.

Obviously, the tradeoff between theoretical motivation and statistical purity is something that must be considered carefully when planning the analysis.

contr1 = Dict(
    :Sex => EffectsCoding(; levels=["Girls", "Boys"]),
    :Test => SeqDiffCoding(;
      levels=["Run", "Star_r", "S20_r", "SLJ", "BPT"]
    ),
  )

Dict{Symbol, StatsModels.AbstractContrasts} with 2 entries:
  :Test => SeqDiffCoding(["Run", "Star_r", "S20_r", "SLJ", "BPT"])
  :Sex  => EffectsCoding(nothing, ["Girls", "Boys"])

m_ovi_SeqDiff_1 = lmm(@formula(zScore ~ 1 + Test + (1 | Child)),
                      dat; contrasts=contr1, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0069	0.0159	-0.43	0.6637	0.6176
Test: Star_r	-0.0197	0.0252	-0.78	0.4335
Test: S20_r	-0.0161	0.0251	-0.64	0.5204
Test: SLJ	-0.0021	0.0249	-0.08	0.9333
Test: BPT	0.0121	0.0248	0.49	0.6252
Residual	0.7774

In this case, any differences between tests identified by the contrasts would be spurious because each test was standardized (i.e., M=0, \(SD\)=1). The differences could also be due to an imbalance in the number of boys and girls or in the number of missing observations for each test.

The primary interest in this study related to interactions of the test contrasts with and age and Sex. We start with age (linear) and its interaction with the four test contrasts.

m_ovi_SeqDiff_2 = lmm(@formula(zScore ~ 1 + Test * a1 + (1 | Child)),
                      dat; contrasts=contr1, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0232	0.0162	-1.44	0.1511	0.6129
Test: Star_r	-0.0247	0.0257	-0.96	0.3369
Test: S20_r	-0.0204	0.0256	-0.80	0.4253
Test: SLJ	0.0032	0.0254	0.13	0.9000
Test: BPT	-0.0042	0.0253	-0.17	0.8679
a1	0.2683	0.0548	4.90	<1e-06
Test: Star_r & a1	0.0836	0.0869	0.96	0.3362
Test: S20_r & a1	0.0713	0.0866	0.82	0.4105
Test: SLJ & a1	-0.0868	0.0862	-1.01	0.3143
Test: BPT & a1	0.2736	0.0859	3.19	0.0014
Residual	0.7764

With the present subsample, only one pairwise difference in developmental gain is significant: the gain is larger for BPT than for SLJ (0.2736). The differences between the other neighboring tests do not reach significance.

Please note that standard errors of this LMM are anti-conservative because the LMM is missing a lot of information in the RES (e..g., contrast-related VCs snd CPs for Child, School, and Cohort.

Next we add the main effect of Sex and its interaction with the four test contrasts.

m_ovi_SeqDiff_3 = lmm(@formula(zScore ~ 1 + Test * (a1 + Sex) + (1 | Child)),
                      dat; contrasts=contr1, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0226	0.0154	-1.47	0.1421	0.5752
Test: Star_r	-0.0254	0.0255	-1.00	0.3196
Test: S20_r	-0.0203	0.0254	-0.80	0.4236
Test: SLJ	0.0036	0.0253	0.14	0.8880
Test: BPT	-0.0043	0.0252	-0.17	0.8635
a1	0.2611	0.0522	5.01	<1e-06
Sex: Boys	0.2153	0.0151	14.28	<1e-45
Test: Star_r & a1	0.0922	0.0863	1.07	0.2850
Test: S20_r & a1	0.0704	0.0860	0.82	0.4134
Test: SLJ & a1	-0.0877	0.0856	-1.02	0.3056
Test: BPT & a1	0.2677	0.0853	3.14	0.0017
Test: Star_r & Sex: Boys	-0.1583	0.0250	-6.34	<1e-09
Test: S20_r & Sex: Boys	0.0414	0.0249	1.66	0.0966
Test: SLJ & Sex: Boys	0.0414	0.0247	1.67	0.0940
Test: BPT & Sex: Boys	0.1473	0.0246	5.98	<1e-08
Residual	0.7709

The significant interactions with Sex reflect mostly differences related to muscle power, where the physiological constitution gives boys an advantage. The sex difference is smaller when coordination and cognition play a role – as in the Star_r test. (Caveat: SEs are estimated with an underspecified RES.)

The final step in this first series is to add the interactions between the three covariates. A significant interaction between any of the four Test contrasts and age (linear) x Sex was hypothesized to reflect a prepubertal signal (i.e., hormones start to rise in girls’ ninth year of life). However, this hypothesis is linked to a specific shape of the interaction: Girls would need to gain more than boys in tests of muscular power.

# we'll be re-using this formula with other contrast specifications,
# so we assign it to its own variable
f_ovi = @formula(zScore ~ 1 + Test * a1 * Sex + (1 | Child))

FormulaTerm
Response:
  zScore(unknown)
Predictors:
  1
  Test(unknown)
  a1(unknown)
  Sex(unknown)
  Test(unknown) & a1(unknown)
  Test(unknown) & Sex(unknown)
  a1(unknown) & Sex(unknown)
  Test(unknown) & a1(unknown) & Sex(unknown)
  (Child)->1 | Child

m_ovi_SeqDiff = lmm(f_ovi, dat; contrasts=contr1, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0227	0.0154	-1.47	0.1412	0.5752
Test: Star_r	-0.0254	0.0255	-1.00	0.3182
Test: S20_r	-0.0204	0.0254	-0.80	0.4218
Test: SLJ	0.0034	0.0253	0.13	0.8940
Test: BPT	-0.0041	0.0251	-0.16	0.8692
a1	0.2610	0.0522	5.00	<1e-06
Sex: Boys	0.2134	0.0154	13.86	<1e-42
Test: Star_r & a1	0.0920	0.0862	1.07	0.2861
Test: S20_r & a1	0.0704	0.0860	0.82	0.4132
Test: SLJ & a1	-0.0878	0.0856	-1.03	0.3049
Test: BPT & a1	0.2678	0.0853	3.14	0.0017
Test: Star_r & Sex: Boys	-0.1576	0.0255	-6.18	<1e-09
Test: S20_r & Sex: Boys	0.0392	0.0254	1.54	0.1234
Test: SLJ & Sex: Boys	0.0353	0.0253	1.40	0.1624
Test: BPT & Sex: Boys	0.1531	0.0251	6.09	<1e-08
a1 & Sex: Boys	0.0314	0.0522	0.60	0.5467
Test: Star_r & a1 & Sex: Boys	-0.0101	0.0862	-0.12	0.9071
Test: S20_r & a1 & Sex: Boys	0.0358	0.0860	0.42	0.6775
Test: SLJ & a1 & Sex: Boys	0.1031	0.0856	1.21	0.2282
Test: BPT & a1 & Sex: Boys	-0.0987	0.0853	-1.16	0.2468
Residual	0.7707

The results are very clear: Despite an abundance of statistical power there is no evidence for the differences between boys and girls in how much they gain in the ninth year of life in these five tests. The authors argue that, in this case, absence of evidence looks very much like evidence of absence of a hypothesized interaction.

In the next two sections we use different contrasts. Does this have a bearing on this result? We still ignore for now that we are looking at anti-conservative test statistics.

2.2 _Helmert coding: `contr2`

The second set of contrasts uses Helmert coding. Helmert coding encodes each level as the difference from the average of the lower levels. With the default order of Test levels we get the following test statistics which we describe in reverse order of appearance in model output

HeC4: 5 - mean(1,2,3,4)
HeC3: 4 - mean(1,2,3)
HeC2: 3 - mean(1,2)
HeC1: 2 - 1

In the model output, HeC1 will be reported first and HeC4 last.

There is some justification for the HeC4 specification in a post-hoc manner because the fifth test (BPT) turned out to be different from the other four tests in that high performance is most likely not only related to physical fitness, but also to overweight/obesity, that is for a subset of children high scores on this test might be indicative of physical unfitness. A priori the SDC4 contrast 5-4 between BPT (5) and SLJ (4) was motivated because conceptually both are tests of the physical fitness component Muscular Power, BPT for upper limbs and SLJ for lower limbs, respectively.

One could argue that there is justification for HeC3 because Run (1), Star_r (2), and S20 (3) involve running but SLJ (4) does not. Sports scientists, however, recoil. For them it does not make much sense to average the different running tests, because they draw on completely different physiological resources; it is a variant of the old apples-and-oranges problem.

The justification for HeC3 is thatRun (1) and Star_r (2) draw more strongly on cardiorespiratory Endurance than S20 (3) due to the longer duration of the runs compared to sprinting for 20 m which is a pure measure of the physical-fitness component Speed. Again, sports scientists are not very happy with this proposal.

Finally, HeC1 contrasts the fitness components Endurance, indicated best by Run (1), and Coordination, indicated by Star_r (2). Endurance (i.e., running for 6 minutes) is considered to be the best indicator of health-related status among the five tests because it is a rather pure measure of cardiorespiratory fitness. The Star_r test requires execution of a pre-instructed sequence of forward, sideways, and backward runs. This coordination of body movements implies a demand on working memory (i.e., remembering the order of these subruns) and executive control processes, but performance also depends on endurance. HeC1 yields a measure of Coordination “corrected” for the contribution of Endurance.

The statistical advantage of Helmert coding is that the resulting contrasts are orthogonal (uncorrelated). This allows for optimal partitioning of variance and statistical power. It is also more efficient to estimate “orthogonal” than “non-orthogonal” random-effect structures.

contr2 = Dict(
  :Sex => EffectsCoding(; levels=["Girls", "Boys"]),
  :Test => HelmertCoding(;
    levels=["Run", "Star_r", "S20_r", "SLJ", "BPT"],
  ),
);

m_ovi_Helmert = lmm(f_ovi, dat; contrasts=contr2, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0227	0.0154	-1.47	0.1412	0.5752
Test: Star_r	-0.0127	0.0127	-1.00	0.3182
Test: S20_r	-0.0110	0.0073	-1.51	0.1318
Test: SLJ	-0.0047	0.0052	-0.91	0.3635
Test: BPT	-0.0036	0.0040	-0.91	0.3611
a1	0.2610	0.0522	5.00	<1e-06
Sex: Boys	0.2134	0.0154	13.86	<1e-42
Test: Star_r & a1	0.0460	0.0431	1.07	0.2861
Test: S20_r & a1	0.0388	0.0248	1.56	0.1181
Test: SLJ & a1	-0.0026	0.0175	-0.15	0.8837
Test: BPT & a1	0.0520	0.0135	3.85	0.0001
Test: Star_r & Sex: Boys	-0.0788	0.0127	-6.18	<1e-09
Test: S20_r & Sex: Boys	-0.0132	0.0073	-1.80	0.0715
Test: SLJ & Sex: Boys	0.0022	0.0052	0.43	0.6672
Test: BPT & Sex: Boys	0.0320	0.0040	8.02	<1e-14
a1 & Sex: Boys	0.0314	0.0522	0.60	0.5467
Test: Star_r & a1 & Sex: Boys	-0.0050	0.0431	-0.12	0.9071
Test: S20_r & a1 & Sex: Boys	0.0102	0.0248	0.41	0.6798
Test: SLJ & a1 & Sex: Boys	0.0309	0.0175	1.77	0.0766
Test: BPT & a1 & Sex: Boys	-0.0012	0.0135	-0.09	0.9289
Residual	0.7707

We forego a detailed discussion of the effects, but note that again none of the interactions between age x Sex with the four test contrasts was significant.

The default labeling of Helmert contrasts may lead to confusions with other contrasts. Therefore, we could provide our own labels:

labels=["c2.1", "c3.12", "c4.123", "c5.1234"]

Once the order of levels is memorized the proposed labelling is very transparent.

2.3 Hypothesis coding: `contr3`

The third set of contrasts uses hypothesis coding. Hypothesis coding allows the user to specify their own a priori contrast matrix, subject to the mathematical constraint that the matrix has full rank. For example, sport scientists agree that the first four tests can be contrasted with BPT, because the difference is akin to a correction of overall physical fitness. However, they want to keep the pairwise comparisons for the first four tests.

HyC1: BPT - mean(1,2,3,4)
HyC2: Star_r - Run_r
HyC3: Run_r - S20_r
HyC4: S20_r - SLJ

contr3 = Dict(
  :Sex => EffectsCoding(; levels=["Girls", "Boys"]),
  :Test => HypothesisCoding(
    [
      -1 -1 -1 -1 +4
      -1 +1  0  0  0
       0 -1 +1  0  0
       0  0 -1 +1  0
    ];
    levels=["Run", "Star_r", "S20_r", "SLJ", "BPT"],
    labels=["BPT-other", "Star-End", "S20-Star", "SLJ-S20"],
  ),
);

m_ovi_Hypo = lmm(f_ovi, dat; contrasts=contr3, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0227	0.0154	-1.47	0.1412	0.5752
Test: BPT-other	-0.0728	0.0797	-0.91	0.3611
Test: Star-End	-0.0254	0.0255	-1.00	0.3182
Test: S20-Star	-0.0204	0.0254	-0.80	0.4218
Test: SLJ-S20	0.0034	0.0253	0.13	0.8940
a1	0.2610	0.0522	5.00	<1e-06
Sex: Boys	0.2134	0.0154	13.86	<1e-42
Test: BPT-other & a1	1.0407	0.2701	3.85	0.0001
Test: Star-End & a1	0.0920	0.0862	1.07	0.2861
Test: S20-Star & a1	0.0704	0.0860	0.82	0.4132
Test: SLJ-S20 & a1	-0.0878	0.0856	-1.03	0.3049
Test: BPT-other & Sex: Boys	0.6392	0.0797	8.02	<1e-14
Test: Star-End & Sex: Boys	-0.1576	0.0255	-6.18	<1e-09
Test: S20-Star & Sex: Boys	0.0392	0.0254	1.54	0.1234
Test: SLJ-S20 & Sex: Boys	0.0353	0.0253	1.40	0.1624
a1 & Sex: Boys	0.0314	0.0522	0.60	0.5467
Test: BPT-other & a1 & Sex: Boys	-0.0241	0.2701	-0.09	0.9289
Test: Star-End & a1 & Sex: Boys	-0.0101	0.0862	-0.12	0.9071
Test: S20-Star & a1 & Sex: Boys	0.0358	0.0860	0.42	0.6775
Test: SLJ-S20 & a1 & Sex: Boys	0.1031	0.0856	1.21	0.2282
Residual	0.7707

With hypothesis coding we must generate our own labels for the contrasts. The default labeling of contrasts is usually not interpretable. Therefore, we provide our own.

Anyway, none of the interactions between age x Sex with the four Test contrasts was significant for these contrasts.

contr1b = Dict(
  :Sex => EffectsCoding(; levels=["Girls", "Boys"]),
  :Test => HypothesisCoding(
    [
      -1 +1  0  0  0
       0 -1 +1  0  0
       0  0 -1 +1  0
       0  0  0 -1 +1
    ];
    levels=["Run", "Star_r", "S20_r", "SLJ", "BPT"],
    labels=["Star-Run", "S20-Star", "SLJ-S20", "BPT-SLJ"],
  ),
);

m_ovi_SeqDiff_v2 = lmm(f_ovi, dat; contrasts=contr1b, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0227	0.0154	-1.47	0.1412	0.5752
Test: Star-Run	-0.0254	0.0255	-1.00	0.3182
Test: S20-Star	-0.0204	0.0254	-0.80	0.4218
Test: SLJ-S20	0.0034	0.0253	0.13	0.8940
Test: BPT-SLJ	-0.0041	0.0251	-0.16	0.8692
a1	0.2610	0.0522	5.00	<1e-06
Sex: Boys	0.2134	0.0154	13.86	<1e-42
Test: Star-Run & a1	0.0920	0.0862	1.07	0.2861
Test: S20-Star & a1	0.0704	0.0860	0.82	0.4132
Test: SLJ-S20 & a1	-0.0878	0.0856	-1.03	0.3049
Test: BPT-SLJ & a1	0.2678	0.0853	3.14	0.0017
Test: Star-Run & Sex: Boys	-0.1576	0.0255	-6.18	<1e-09
Test: S20-Star & Sex: Boys	0.0392	0.0254	1.54	0.1234
Test: SLJ-S20 & Sex: Boys	0.0353	0.0253	1.40	0.1624
Test: BPT-SLJ & Sex: Boys	0.1531	0.0251	6.09	<1e-08
a1 & Sex: Boys	0.0314	0.0522	0.60	0.5467
Test: Star-Run & a1 & Sex: Boys	-0.0101	0.0862	-0.12	0.9071
Test: S20-Star & a1 & Sex: Boys	0.0358	0.0860	0.42	0.6775
Test: SLJ-S20 & a1 & Sex: Boys	0.1031	0.0856	1.21	0.2282
Test: BPT-SLJ & a1 & Sex: Boys	-0.0987	0.0853	-1.16	0.2468
Residual	0.7707

m_zcp_SeqD = lmm(@formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Child)),
                 dat; contrasts=contr1b, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0230	0.0154	-1.49	0.1363	0.5931
Test: Star-Run	-0.0259	0.0250	-1.03	0.3013	0.4040
Test: S20-Star	-0.0196	0.0236	-0.83	0.4070	0.2097
Test: SLJ-S20	0.0026	0.0230	0.11	0.9089	0.0000
Test: BPT-SLJ	-0.0037	0.0260	-0.14	0.8874	0.5319
a1	0.2607	0.0522	5.00	<1e-06
Sex: Boys	0.2135	0.0154	13.86	<1e-42
Test: Star-Run & a1	0.0937	0.0847	1.11	0.2687
Test: S20-Star & a1	0.0702	0.0800	0.88	0.3801
Test: SLJ-S20 & a1	-0.0874	0.0779	-1.12	0.2618
Test: BPT-SLJ & a1	0.2677	0.0880	3.04	0.0023
Test: Star-Run & Sex: Boys	-0.1573	0.0250	-6.28	<1e-09
Test: S20-Star & Sex: Boys	0.0384	0.0236	1.62	0.1044
Test: SLJ-S20 & Sex: Boys	0.0358	0.0230	1.56	0.1196
Test: BPT-SLJ & Sex: Boys	0.1530	0.0260	5.90	<1e-08
a1 & Sex: Boys	0.0332	0.0522	0.64	0.5243
Test: Star-Run & a1 & Sex: Boys	-0.0126	0.0847	-0.15	0.8822
Test: S20-Star & a1 & Sex: Boys	0.0399	0.0800	0.50	0.6176
Test: SLJ-S20 & a1 & Sex: Boys	0.0996	0.0779	1.28	0.2012
Test: BPT-SLJ & a1 & Sex: Boys	-0.0991	0.0880	-1.13	0.2599
Residual	0.7010

m_zcp_SeqD_2 = lmm(@formula(zScore ~ 1 + Test * a1 * Sex + (0 + Test | Child)),
                   dat; contrasts=contr1b, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0233	0.0154	-1.51	0.1309
Test: Star-Run	-0.0263	0.0262	-1.00	0.3165
Test: S20-Star	-0.0195	0.0244	-0.80	0.4244
Test: SLJ-S20	0.0015	0.0220	0.07	0.9438
Test: BPT-SLJ	-0.0027	0.0250	-0.11	0.9131
a1	0.2601	0.0522	4.98	<1e-06
Sex: Boys	0.2133	0.0154	13.83	<1e-42
Test: Star-Run & a1	0.0929	0.0888	1.05	0.2952
Test: S20-Star & a1	0.0691	0.0825	0.84	0.4022
Test: SLJ-S20 & a1	-0.0871	0.0744	-1.17	0.2422
Test: BPT-SLJ & a1	0.2672	0.0846	3.16	0.0016
Test: Star-Run & Sex: Boys	-0.1565	0.0262	-5.96	<1e-08
Test: S20-Star & Sex: Boys	0.0378	0.0244	1.55	0.1215
Test: SLJ-S20 & Sex: Boys	0.0358	0.0220	1.63	0.1032
Test: BPT-SLJ & Sex: Boys	0.1534	0.0250	6.14	<1e-09
a1 & Sex: Boys	0.0342	0.0522	0.66	0.5119
Test: Star-Run & a1 & Sex: Boys	-0.0152	0.0888	-0.17	0.8637
Test: S20-Star & a1 & Sex: Boys	0.0445	0.0825	0.54	0.5892
Test: SLJ-S20 & a1 & Sex: Boys	0.0932	0.0744	1.25	0.2105
Test: BPT-SLJ & a1 & Sex: Boys	-0.0971	0.0846	-1.15	0.2515
Test: Run					0.7710
Test: Star_r					0.8121
Test: S20_r					0.7850
Test: SLJ					0.7233
Test: BPT					0.6865
Residual	0.5941

m_cpx_0_SeqDiff = lmm(@formula(zScore ~ 1 + Test * a1 * Sex + (0 + Test | Child)),
                      dat; contrasts=contr1b, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0233	0.0154	-1.51	0.1309
Test: Star-Run	-0.0263	0.0262	-1.00	0.3165
Test: S20-Star	-0.0195	0.0244	-0.80	0.4244
Test: SLJ-S20	0.0015	0.0220	0.07	0.9438
Test: BPT-SLJ	-0.0027	0.0250	-0.11	0.9131
a1	0.2601	0.0522	4.98	<1e-06
Sex: Boys	0.2133	0.0154	13.83	<1e-42
Test: Star-Run & a1	0.0929	0.0888	1.05	0.2952
Test: S20-Star & a1	0.0691	0.0825	0.84	0.4022
Test: SLJ-S20 & a1	-0.0871	0.0744	-1.17	0.2422
Test: BPT-SLJ & a1	0.2672	0.0846	3.16	0.0016
Test: Star-Run & Sex: Boys	-0.1565	0.0262	-5.96	<1e-08
Test: S20-Star & Sex: Boys	0.0378	0.0244	1.55	0.1215
Test: SLJ-S20 & Sex: Boys	0.0358	0.0220	1.63	0.1032
Test: BPT-SLJ & Sex: Boys	0.1534	0.0250	6.14	<1e-09
a1 & Sex: Boys	0.0342	0.0522	0.66	0.5119
Test: Star-Run & a1 & Sex: Boys	-0.0152	0.0888	-0.17	0.8637
Test: S20-Star & a1 & Sex: Boys	0.0445	0.0825	0.54	0.5892
Test: SLJ-S20 & a1 & Sex: Boys	0.0932	0.0744	1.25	0.2105
Test: BPT-SLJ & a1 & Sex: Boys	-0.0971	0.0846	-1.15	0.2515
Test: Run					0.7710
Test: Star_r					0.8121
Test: S20_r					0.7850
Test: SLJ					0.7233
Test: BPT					0.6865
Residual	0.5941

VarCorr(m_cpx_0_SeqDiff)

	Column	Variance	Std.Dev	Corr.
Child	Test: Run	0.594487	0.771030
	Test: Star_r	0.659509	0.812102	+0.56
	Test: S20_r	0.616155	0.784955	+0.61	+0.70
	Test: SLJ	0.523118	0.723269	+0.68	+0.70	+0.84
	Test: BPT	0.471300	0.686513	+0.27	+0.41	+0.46	+0.53
Residual		0.352937	0.594085

m_cpx_0_SeqDiff.PCA

(Child = 
Principal components based on correlation matrix
 Test: Run     1.0    .     .     .     .
 Test: Star_r  0.56  1.0    .     .     .
 Test: S20_r   0.61  0.7   1.0    .     .
 Test: SLJ     0.68  0.7   0.84  1.0    .
 Test: BPT     0.27  0.41  0.46  0.53  1.0

Normalized cumulative variances:
[0.6711, 0.8211, 0.9065, 0.9704, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5
 Test: Run     -0.42   0.47   0.69   0.31   0.18
 Test: Star_r  -0.46   0.1   -0.63   0.62  -0.0
 Test: S20_r   -0.49   0.07  -0.23  -0.6    0.58
 Test: SLJ     -0.51   0.02   0.03  -0.35  -0.78
 Test: BPT     -0.34  -0.87   0.28   0.18   0.12,)

f_cpx_1 = @formula(zScore ~ 1 + Test * a1 * Sex + (1 + Test | Child))

FormulaTerm
Response:
  zScore(unknown)
Predictors:
  1
  Test(unknown)
  a1(unknown)
  Sex(unknown)
  Test(unknown) & a1(unknown)
  Test(unknown) & Sex(unknown)
  a1(unknown) & Sex(unknown)
  Test(unknown) & a1(unknown) & Sex(unknown)
  (Test,Child)->(1 + Test) | Child

m_cpx_1_SeqDiff = lmm(f_cpx_1, dat; contrasts=contr1b, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0233	0.0154	-1.51	0.1309	0.6089
Test: Star-Run	-0.0263	0.0262	-1.00	0.3165	0.6814
Test: S20-Star	-0.0195	0.0244	-0.80	0.4244	0.5423
Test: SLJ-S20	0.0015	0.0220	0.07	0.9438	0.3164
Test: BPT-SLJ	-0.0027	0.0250	-0.11	0.9131	0.6141
a1	0.2601	0.0522	4.98	<1e-06
Sex: Boys	0.2133	0.0154	13.83	<1e-42
Test: Star-Run & a1	0.0929	0.0888	1.05	0.2952
Test: S20-Star & a1	0.0691	0.0825	0.84	0.4022
Test: SLJ-S20 & a1	-0.0871	0.0744	-1.17	0.2422
Test: BPT-SLJ & a1	0.2672	0.0846	3.16	0.0016
Test: Star-Run & Sex: Boys	-0.1565	0.0262	-5.96	<1e-08
Test: S20-Star & Sex: Boys	0.0378	0.0244	1.55	0.1215
Test: SLJ-S20 & Sex: Boys	0.0358	0.0220	1.63	0.1032
Test: BPT-SLJ & Sex: Boys	0.1534	0.0250	6.14	<1e-09
a1 & Sex: Boys	0.0342	0.0522	0.66	0.5119
Test: Star-Run & a1 & Sex: Boys	-0.0152	0.0888	-0.17	0.8637
Test: S20-Star & a1 & Sex: Boys	0.0445	0.0825	0.54	0.5892
Test: SLJ-S20 & a1 & Sex: Boys	0.0932	0.0744	1.25	0.2105
Test: BPT-SLJ & a1 & Sex: Boys	-0.0971	0.0846	-1.15	0.2515
Residual	0.6308

m_cpx_1_SeqDiff.PCA

(Child = 
Principal components based on correlation matrix
 (Intercept)      1.0     .      .      .      .
 Test: Star-Run   0.13   1.0     .      .      .
 Test: S20-Star   0.03  -0.5    1.0     .      .
 Test: SLJ-S20   -0.1   -0.21  -0.38   1.0     .
 Test: BPT-SLJ   -0.37   0.12  -0.13   0.07   1.0

Normalized cumulative variances:
[0.3251, 0.6022, 0.831, 0.9488, 1.0]

Component loadings
                   PC1    PC2    PC3    PC4    PC5
 (Intercept)     -0.24   0.63  -0.27   0.68  -0.11
 Test: Star-Run   0.46   0.55   0.35  -0.11   0.59
 Test: S20-Star  -0.67  -0.24   0.24   0.18   0.63
 Test: SLJ-S20    0.33  -0.29  -0.74   0.17   0.48
 Test: BPT-SLJ    0.41  -0.41   0.44   0.68  -0.09,)

2.4 PCA-based hypothesis coding: `contr4`

The fourth set of contrasts uses hypothesis coding to specify the set of contrasts implementing the loadings of the four principal components of the published LMM based on test scores, not test effects (contrasts); in other words, roughly according to their signs. This is actually a very interesting and plausible solution nobody had proposed a priori.

PC1: BPT - Run_r
PC2: (Star_r + S20_r + SLJ) - (BPT + Run_r)
PC3: Star_r - (S20_r + SLJ)
PC4: S20_r - SLJ

PC1 contrasts the worst and the best indicator of physical health; PC2 contrasts these two against the core indicators of physical fitness; PC3 contrasts the cognitive and the physical tests within the narrow set of physical fitness components; and PC4, finally, contrasts two types of lower muscular fitness differing in speed and power.

contr4 = Dict(
  :Sex => EffectsCoding(; levels=["Girls", "Boys"]),
  :Test => HypothesisCoding(
    [
      -1  0  0  0 +1
      -3 +2 +2 +2 -3
       0 +2 -1 -1  0
       0  0 +1 -1  0
    ];
    levels=["Run", "Star_r", "S20_r", "SLJ", "BPT"],
    labels=["c5.1", "c234.15", "c2.34", "c3.4"],
  ),
);

m_cpx_1_PC = lmm(f_cpx_1, dat; contrasts=contr4, progress)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0233	0.0154	-1.51	0.1309	0.6060
Test: c5.1	-0.0469	0.0282	-1.66	0.0964	0.8079
Test: c234.15	-0.0917	0.1001	-0.92	0.3596	2.4743
Test: c2.34	0.0374	0.0431	0.87	0.3853	0.9549
Test: c3.4	-0.0015	0.0220	-0.07	0.9438	0.2541
a1	0.2601	0.0522	4.98	<1e-06
Sex: Boys	0.2133	0.0154	13.83	<1e-42
Test: c5.1 & a1	0.3422	0.0957	3.58	0.0003
Test: c234.15 & a1	-0.3668	0.3392	-1.08	0.2795
Test: c2.34 & a1	-0.0511	0.1457	-0.35	0.7257
Test: c3.4 & a1	0.0871	0.0744	1.17	0.2422
Test: c5.1 & Sex: Boys	0.0704	0.0282	2.50	0.0126
Test: c234.15 & Sex: Boys	-0.9277	0.1001	-9.27	<1e-19
Test: c2.34 & Sex: Boys	-0.1113	0.0431	-2.58	0.0098
Test: c3.4 & Sex: Boys	-0.0358	0.0220	-1.63	0.1032
a1 & Sex: Boys	0.0342	0.0522	0.66	0.5119
Test: c5.1 & a1 & Sex: Boys	0.0255	0.0957	0.27	0.7902
Test: c234.15 & a1 & Sex: Boys	0.1968	0.3392	0.58	0.5618
Test: c2.34 & a1 & Sex: Boys	-0.1823	0.1457	-1.25	0.2107
Test: c3.4 & a1 & Sex: Boys	-0.0932	0.0744	-1.25	0.2105
Residual	0.6447

VarCorr(m_cpx_1_PC)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.367259	0.606019
	Test: c5.1	0.652653	0.807870	-0.20
	Test: c234.15	6.122295	2.474327	+0.41	-0.17
	Test: c2.34	0.911834	0.954900	-0.00	-0.01	+0.20
	Test: c3.4	0.064591	0.254148	+0.12	-0.02	+0.54	-0.12
Residual		0.415660	0.644717

m_cpx_1_PC.PCA

(Child = 
Principal components based on correlation matrix
 (Intercept)     1.0     .      .      .      .
 Test: c5.1     -0.2    1.0     .      .      .
 Test: c234.15   0.41  -0.17   1.0     .      .
 Test: c2.34    -0.0   -0.01   0.2    1.0     .
 Test: c3.4      0.12  -0.02   0.54  -0.12   1.0

Normalized cumulative variances:
[0.361, 0.5789, 0.7834, 0.938, 1.0]

Component loadings
                  PC1    PC2    PC3    PC4    PC5
 (Intercept)    -0.48  -0.16  -0.33  -0.74   0.3
 Test: c5.1      0.27   0.34   0.68  -0.58  -0.07
 Test: c234.15  -0.66  -0.03   0.23   0.04  -0.71
 Test: c2.34    -0.09  -0.77   0.54   0.11   0.3
 Test: c3.4     -0.51   0.51   0.27   0.32   0.55,)

There is a numerical interaction with a z-value > 2.0 for the first PCA (i.e., BPT - Run_r). This interaction would really need to be replicated to be taken seriously. It is probably due to larger “unfitness” gains in boys than girls (i.e., in BPT) relative to the slightly larger health-related “fitness” gains of girls than boys (i.e., in Run_r).

contr4b = Dict(
    :Sex => EffectsCoding(; levels=["Girls", "Boys"]),
    :Test => HypothesisCoding(
      [
        0.49 -0.04  0.20  0.03 -0.85
        0.70 -0.56 -0.21 -0.13  0.37
        0.31  0.68 -0.56 -0.35  0.00
        0.04  0.08  0.61 -0.78  0.13
      ];
      levels=["Run", "Star_r", "S20_r", "SLJ", "BPT"],
      labels=["c5.1", "c234.15", "c12.34", "c3.4"],
    )
);

# LN_NEWUOA does a poor job here
m_cpx_1_PC_2 = lmm(f_cpx_1, dat; contrasts=contr4b, progress, optimizer=:LN_BOBYQA)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0235	0.0157	-1.50	0.1338	0.6342
Test: c5.1	0.0249	0.0201	1.24	0.2163	0.6650
Test: c234.15	0.0183	0.0185	0.99	0.3233	0.5625
Test: c12.34	0.0259	0.0174	1.48	0.1380	0.4975
Test: c3.4	0.0010	0.0154	0.06	0.9486	0.3590
a1	0.2514	0.0531	4.74	<1e-05
Sex: Boys	0.2049	0.0157	13.07	<1e-38
Test: c5.1 & a1	-0.2386	0.0682	-3.50	0.0005
Test: c234.15 & a1	0.0094	0.0627	0.15	0.8804
Test: c12.34 & a1	-0.0638	0.0590	-1.08	0.2791
Test: c3.4 & a1	0.0822	0.0523	1.57	0.1159
Test: c5.1 & Sex: Boys	-0.0911	0.0201	-4.53	<1e-05
Test: c234.15 & Sex: Boys	0.1606	0.0185	8.67	<1e-17
Test: c12.34 & Sex: Boys	-0.0056	0.0174	-0.32	0.7475
Test: c3.4 & Sex: Boys	-0.0058	0.0154	-0.38	0.7060
a1 & Sex: Boys	0.0368	0.0531	0.69	0.4878
Test: c5.1 & a1 & Sex: Boys	-0.0055	0.0682	-0.08	0.9352
Test: c234.15 & a1 & Sex: Boys	-0.0101	0.0627	-0.16	0.8724
Test: c12.34 & a1 & Sex: Boys	-0.0724	0.0590	-1.23	0.2192
Test: c3.4 & a1 & Sex: Boys	-0.0784	0.0523	-1.50	0.1338
Residual	0.5587

VarCorr(m_cpx_1_PC_2)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.402164	0.634164
	Test: c5.1	0.442264	0.665029	+0.27
	Test: c234.15	0.316441	0.562531	-0.27	-0.04
	Test: c12.34	0.247512	0.497506	-0.06	+0.10	-0.04
	Test: c3.4	0.128915	0.359047	-0.03	-0.07	-0.13	-0.04
Residual		0.312154	0.558708

m_cpx_1_PC_2.PCA

(Child = 
Principal components based on correlation matrix
 (Intercept)     1.0     .      .      .      .
 Test: c5.1      0.27   1.0     .      .      .
 Test: c234.15  -0.27  -0.04   1.0     .      .
 Test: c12.34   -0.06   0.1   -0.04   1.0     .
 Test: c3.4     -0.03  -0.07  -0.13  -0.04   1.0

Normalized cumulative variances:
[0.2811, 0.5094, 0.7123, 0.8807, 1.0]

Component loadings
                  PC1    PC2    PC3    PC4    PC5
 (Intercept)    -0.68  -0.02  -0.28   0.07   0.68
 Test: c5.1     -0.51   0.44  -0.01  -0.6   -0.44
 Test: c234.15   0.52   0.43  -0.22  -0.49   0.5
 Test: c12.34   -0.08   0.37   0.87   0.13   0.27
 Test: c3.4     -0.02  -0.7    0.33  -0.62   0.16,)

# LN_BOBYQA does a poor job here
m_zcp_1_PC_2 = lmm(@formula(zScore ~ 1 + Test*a1*Sex + zerocorr(1 + Test | Child)),
                   dat; contrasts=contr4b, progress, optimizer=:LN_NEWUOA)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0231	0.0155	-1.49	0.1368	0.6035
Test: c5.1	0.0247	0.0200	1.24	0.2167	0.5405
Test: c234.15	0.0177	0.0184	0.96	0.3372	0.4135
Test: c12.34	0.0262	0.0174	1.50	0.1330	0.3318
Test: c3.4	0.0006	0.0156	0.04	0.9685	0.0000
a1	0.2522	0.0525	4.80	<1e-05
Sex: Boys	0.2048	0.0155	13.21	<1e-39
Test: c5.1 & a1	-0.2391	0.0678	-3.53	0.0004
Test: c234.15 & a1	0.0096	0.0624	0.15	0.8782
Test: c12.34 & a1	-0.0634	0.0590	-1.08	0.2819
Test: c3.4 & a1	0.0821	0.0529	1.55	0.1205
Test: c5.1 & Sex: Boys	-0.0907	0.0200	-4.53	<1e-05
Test: c234.15 & Sex: Boys	0.1611	0.0184	8.74	<1e-17
Test: c12.34 & Sex: Boys	-0.0060	0.0174	-0.35	0.7293
Test: c3.4 & Sex: Boys	-0.0059	0.0156	-0.38	0.7052
a1 & Sex: Boys	0.0357	0.0525	0.68	0.4959
Test: c5.1 & a1 & Sex: Boys	-0.0074	0.0678	-0.11	0.9134
Test: c234.15 & a1 & Sex: Boys	-0.0122	0.0624	-0.19	0.8456
Test: c12.34 & a1 & Sex: Boys	-0.0710	0.0590	-1.20	0.2282
Test: c3.4 & a1 & Sex: Boys	-0.0823	0.0529	-1.56	0.1198
Residual	0.6712

VarCorr(m_zcp_1_PC_2)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.364186	0.603478
	Test: c5.1	0.292170	0.540527	.
	Test: c234.15	0.170951	0.413462	.	.
	Test: c12.34	0.110087	0.331793	.	.	.
	Test: c3.4	0.000000	0.000000	.	.	.	.
Residual		0.450453	0.671158

likelihoodratiotest(m_zcp_1_PC_2, m_cpx_1_PC_2)

	model-dof	-2 logLik	χ²	χ²-dof	P(>χ²)
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + zerocorr(1 + Test \| Child)	26	-24978
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-24837	142	10	<1e-24

3 Other topics

3.1 Contrasts are re-parameterizations of the same model

The choice of contrast does not affect the model objective, in other words, they all yield the same goodness of fit. It does not matter whether a contrast is orthogonal or not.

[
  objective(m_ovi_SeqDiff),
  objective(m_ovi_Helmert),
  objective(m_ovi_Hypo),
]

3-element Vector{Float64}:
 25095.94518164001
 25095.94518164003
 25095.94518164004

3.2 VCs and CPs depend on contrast coding

Trivially, the meaning of a contrast depends on its definition. Consequently, the contrast specification has a big effect on the random-effect structure. As an illustration, we refit the LMMs with variance components (VCs) and correlation parameters (CPs) for Child-related contrasts of Test. Unfortunately, it is not easy, actually rather quite difficult, to grasp the meaning of correlations of contrast-based effects; they represent two-way interactions.

f_Child = @formula(zScore ~ 1 + Test * a1 * Sex + (1 + Test | Child))
m_Child_SDC = lmm(f_Child, dat; contrasts=contr1, progress, optimizer=:LN_BOBYQA)
m_Child_HeC = lmm(f_Child, dat; contrasts=contr2, progress, optimizer=:LN_NELDERMEAD)
m_Child_HyC = lmm(f_Child, dat; contrasts=contr3, progress, optimizer=:LN_NELDERMEAD)
m_Child_PCA = lmm(f_Child, dat; contrasts=contr4, progress, optimizer=:LN_NELDERMEAD);

VarCorr(m_Child_SDC)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.372579	0.610393
	Test: Star_r	0.482005	0.694265	+0.12
	Test: S20_r	0.311722	0.558321	+0.03	-0.50
	Test: SLJ	0.117813	0.343239	-0.09	-0.19	-0.38
	Test: BPT	0.394746	0.628289	-0.36	+0.12	-0.12	+0.02
Residual		0.389052	0.623740

VarCorr(m_Child_HeC)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.441401	0.664380
	Test: Star_r	0.292554	0.540883	+0.07
	Test: S20_r	0.083714	0.289333	+0.06	+0.05
	Test: SLJ	0.034662	0.186176	+0.01	-0.02	+0.10
	Test: BPT	0.034710	0.186306	-0.22	+0.04	+0.04	+0.09
Residual		0.044945	0.212002

VarCorr(m_Child_HyC)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.363751	0.603118
	Test: BPT-other	6.119207	2.473703	-0.36
	Test: Star-End	0.393735	0.627483	+0.14	+0.11
	Test: S20-Star	0.223449	0.472704	+0.04	+0.03	-0.51
	Test: SLJ-S20	0.029538	0.171865	-0.19	+0.42	-0.41	-0.36
Residual		0.433188	0.658170

VarCorr(m_Child_PCA)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.363796	0.603155
	Test: c5.1	0.618023	0.786144	-0.20
	Test: c234.15	5.602961	2.367058	+0.43	-0.18
	Test: c2.34	0.807921	0.898844	-0.00	-0.01	+0.22
	Test: c3.4	0.029970	0.173119	+0.18	-0.03	+0.82	-0.19
Residual		0.432973	0.658007

The CPs for the various contrasts are in line with expectations. For the SDC we observe substantial negative CPs between most neighboring contrasts. For the orthogonal HeC, all CPs are small; they are uncorrelated. HyC contains some of the SDC contrasts and we observe again the negative CPs. Most CPs for the (roughly) PCA-based contrasts are also small, but there are exceptions; for example, there is a sizeable CP of +.43 between GM and the core of adjusted physical fitness (c234.15).

Do these differences in CPs imply that we can move to zcpLMMs when we have orthogonal contrasts? We pursue this question with by refitting the four LMMs with zerocorr() and compare the goodness of fit.

f_Child0 = @formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Child))
m_Child_SDC0 = lmm(f_Child0, dat; contrasts=contr1, progress, optimizer=:LN_NEWUOA)
m_Child_HeC0 = lmm(f_Child0, dat; contrasts=contr2, progress, optimizer=:LN_NEWUOA)
m_Child_HyC0 = lmm(f_Child0, dat; contrasts=contr3, progress, optimizer=:LN_NEWUOA)
m_Child_PCA0 = lmm(f_Child0, dat; contrasts=contr4, progress, optimizer=:LN_NEWUOA)

	Est.	SE	z	p	σ_Child
(Intercept)	-0.0229	0.0154	-1.49	0.1373	0.6027
Test: c5.1	-0.0467	0.0282	-1.65	0.0980	0.7812
Test: c234.15	-0.0889	0.1002	-0.89	0.3751	2.3592
Test: c2.34	0.0380	0.0431	0.88	0.3775	0.8966
Test: c3.4	-0.0016	0.0220	-0.07	0.9425	0.1697
a1	0.2610	0.0522	5.00	<1e-06
Sex: Boys	0.2133	0.0154	13.84	<1e-42
Test: c5.1 & a1	0.3439	0.0956	3.60	0.0003
Test: c234.15 & a1	-0.3675	0.3394	-1.08	0.2790
Test: c2.34 & a1	-0.0518	0.1456	-0.36	0.7222
Test: c3.4 & a1	0.0866	0.0745	1.16	0.2448
Test: c5.1 & Sex: Boys	0.0700	0.0282	2.48	0.0131
Test: c234.15 & Sex: Boys	-0.9289	0.1002	-9.27	<1e-19
Test: c2.34 & Sex: Boys	-0.1125	0.0431	-2.61	0.0090
Test: c3.4 & Sex: Boys	-0.0359	0.0220	-1.63	0.1023
a1 & Sex: Boys	0.0336	0.0522	0.64	0.5199
Test: c5.1 & a1 & Sex: Boys	0.0283	0.0956	0.30	0.7672
Test: c234.15 & a1 & Sex: Boys	0.2031	0.3394	0.60	0.5496
Test: c2.34 & a1 & Sex: Boys	-0.1802	0.1456	-1.24	0.2158
Test: c3.4 & a1 & Sex: Boys	-0.0992	0.0745	-1.33	0.1828
Residual	0.6590

likelihoodratiotest(m_Child_SDC0, m_Child_SDC)

	model-dof	-2 logLik	χ²	χ²-dof	P(>χ²)
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + zerocorr(1 + Test \| Child)	26	-24988
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-24837	152	10	<1e-26

likelihoodratiotest(m_Child_HeC0, m_Child_HeC)

	model-dof	-2 logLik	χ²	χ²-dof	P(>χ²)
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + zerocorr(1 + Test \| Child)	26	-24982
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-24837	146	10	<1e-25

likelihoodratiotest(m_Child_HyC0, m_Child_HyC)

	model-dof	-2 logLik	χ²	χ²-dof	P(>χ²)
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + zerocorr(1 + Test \| Child)	26	-25000
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-24837	163	10	<1e-29

likelihoodratiotest(m_Child_PCA0, m_Child_PCA)

	model-dof	-2 logLik	χ²	χ²-dof	P(>χ²)
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + zerocorr(1 + Test \| Child)	26	-24974
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-24837	138	10	<1e-23

Obviously, we can not drop CPs from any of the LMMs. The full LMMs all have the same objective, but we can compare the goodness-of-fit statistics of zcpLMMs more directly.

zcpLMM = ["SDC0", "HeC0", "HyC0", "PCA0"]
mods = [m_Child_SDC0, m_Child_HeC0, m_Child_HyC0, m_Child_PCA0]
gof_summary = sort!(
  DataFrame(;
    zcpLMM=zcpLMM,
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    BIC=bic.(mods),
  ),
  :deviance,
)

4×5 DataFrame

Row	zcpLMM	dof	deviance	AIC	BIC
	String	Int64	Float64	Float64	Float64
1	PCA0	26	24974.2	25026.2	25212.9
2	HeC0	26	24982.5	25034.5	25221.2
3	SDC0	26	24988.4	25040.4	25227.0
4	HyC0	26	24999.6	25051.6	25238.3

The best fit was obtained for the PCA-based zcpLMM, followed by the orthogonal Helmert coding; the two non-orthogonal codings (SDC and HyC) fit worst. This ordering is consistent with the idea that orthogonal contrasts leave less information in the correlation parameters, so forcing the CPs to zero costs them less. Keep in mind, though, that the exact ordering depends on the subsample.

3.3 VCs and CPs depend on random factor

VCs and CPs resulting from a set of test contrasts can also be estimated for the random factor School. Of course, these VCs and CPs may look different from the ones we just estimated for Child.

The effect of age (i.e., developmental gain) varies within School. Therefore, we also include its VCs and CPs in this model; the school-related VC for Sex was not significant.

f_School = @formula(zScore ~ 1 + Test * a1 * Sex + (1 + Test + a1 | School))
m_School_SeqDiff = lmm(f_School, dat; contrasts=contr1, progress)
m_School_Helmert = lmm(f_School, dat; contrasts=contr2, progress)
m_School_Hypo = lmm(f_School, dat; contrasts=contr3, progress)
m_School_PCA = lmm(f_School, dat; contrasts=contr4, progress);

VarCorr(m_School_SeqDiff)

	Column	Variance	Std.Dev	Corr.
School	(Intercept)	0.0928218	0.3046667
	Test: Star_r	0.0849545	0.2914695	+0.14
	Test: S20_r	0.0310065	0.1760865	+0.15	-0.77
	Test: SLJ	0.0045536	0.0674804	-0.65	-0.40	-0.08
	Test: BPT	0.0270604	0.1645004	-0.60	-0.01	+0.30	+0.10
	a1	0.6320271	0.7950013	-0.11	+0.07	-0.12	+0.62	+0.12
Residual		0.7743279	0.8799590

VarCorr(m_School_Helmert)

	Column	Variance	Std.Dev	Corr.
School	(Intercept)	0.0927876	0.3046105
	Test: Star_r	0.0212455	0.1457583	+0.14
	Test: S20_r	0.0013936	0.0373308	+0.43	+0.09
	Test: SLJ	0.0002230	0.0149329	-0.20	-0.35	+0.51
	Test: BPT	0.0015703	0.0396266	-0.55	-0.09	+0.49	+0.80
	a1	0.6309028	0.7942939	-0.11	+0.07	-0.10	+0.56	+0.23
Residual		0.7743539	0.8799738

VarCorr(m_School_Hypo)

	Column	Variance	Std.Dev	Corr.
School	(Intercept)	0.0928339	0.3046865
	Test: BPT-other	0.6280309	0.7924840	-0.55
	Test: Star-End	0.0849194	0.2914094	+0.14	-0.09
	Test: S20-Star	0.0309923	0.1760462	+0.15	+0.39	-0.77
	Test: SLJ-S20	0.0045594	0.0675236	-0.65	+0.16	-0.40	-0.08
	a1	0.6318621	0.7948975	-0.11	+0.23	+0.07	-0.12	+0.62
Residual		0.7743262	0.8799581

VarCorr(m_School_PCA)

	Column	Variance	Std.Dev	Corr.
School	(Intercept)	0.0928313	0.3046823
	Test: c5.1	0.0690266	0.2627291	-0.29
	Test: c234.15	0.8210628	0.9061251	+0.54	+0.02
	Test: c2.34	0.1245499	0.3529163	-0.03	+0.05	+0.81
	Test: c3.4	0.0045554	0.0674939	+0.65	+0.18	+0.54	+0.11
	a1	0.6318655	0.7948997	-0.11	+0.23	-0.07	-0.00	-0.62
Residual		0.7743302	0.8799604

We compare again how much of the fit resides in the CPs.

f_School0 = @formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test + a1 | School))
m_School_SDC0 = lmm(f_School0, dat; contrasts=contr1, progress)
m_School_HeC0 = lmm(f_School0, dat; contrasts=contr2, progress)
m_School_HyC0 = lmm(f_School0, dat; contrasts=contr3, progress)
m_School_PCA0 = lmm(f_School0, dat; contrasts=contr4, progress);

zcpLMM2 = ["SDC0", "HeC0", "HyC0", "PCA0"]
mods2 = [
  m_School_SDC0, m_School_HeC0, m_School_HyC0, m_School_PCA0
]
gof_summary2 = sort!(
  DataFrame(;
    zcpLMM=zcpLMM2,
    dof=dof.(mods2),
    deviance=deviance.(mods2),
    AIC=aic.(mods2),
    BIC=bic.(mods2),
  ),
  :deviance,
)

4×5 DataFrame

Row	zcpLMM	dof	deviance	AIC	BIC
	String	Int64	Float64	Float64	Float64
1	HeC0	27	26078.4	26132.4	26326.2
2	PCA0	27	26085.2	26139.2	26333.1
3	SDC0	27	26088.7	26142.7	26336.5
4	HyC0	27	26089.2	26143.2	26337.1

For the random factor School the Helmert contrast, followed by PCA-based contrasts, has the least information in the CPs; the non-orthogonal SDC and HyC leave the largest contribution in the CPs. Interesting.

4 That’s it

That’s it for this tutorial. It is time to try your own contrast coding. You can use these data; there are many alternatives to set up hypotheses for the five tests. Of course and even better, code up some contrasts for data of your own.

Have fun!

Fühner, T., Granacher, U., Golle, K., & Kliegl, R. (2021). Age and sex effects in physical fitness components of 108,295 third graders including 515 primary schools and 9 cohorts. Scientific Reports, 11(1). https://doi.org/10.1038/s41598-021-97000-4

4.1 Exercises

Same fit, different coefficients. Fit the model with two different contrast schemes for Test (e.g. SeqDiffCoding and HelmertCoding). Confirm that the log-likelihood is identical while the fixed-effect coefficients differ, and explain why.

Solution

Different full-rank coding schemes for the same factor span the same model space, so the fitted values and log-likelihood are identical; only the parametrization (which comparisons the coefficients represent) changes. Choose the scheme whose coefficients answer your research question directly.

Match coding to a hypothesis. Suppose you want each coefficient to test a specific, theory-driven comparison among the levels of Test. Which coding scheme on this page supports that, and what do you have to supply?

Solution

HypothesisCoding lets you state the contrasts directly: you supply a matrix whose rows are the comparisons you care about (and, for interpretability, ensure they are sensibly scaled). The PCA-based HypothesisCoding variant on this page shows one way to derive such contrasts empirically.

This page was rendered from git revision 29e8d33 using Quarto 1.10.18 and Julia 1.12.7.

1 Setup

1.1 Packages and functions

1.2 Preprocessing

2 Contrast coding

2.1 Sequential difference coding: contr1

2.2 _Helmert coding: contr2

2.3 Hypothesis coding: contr3

2.4 PCA-based hypothesis coding: contr4

3 Other topics

3.1 Contrasts are re-parameterizations of the same model

3.2 VCs and CPs depend on contrast coding

3.3 VCs and CPs depend on random factor

4 That’s it

4.1 Exercises

2.1 Sequential difference coding: `contr1`

2.2 _Helmert coding: `contr2`

2.3 Hypothesis coding: `contr3`

2.4 PCA-based hypothesis coding: `contr4`