Mixed Models Tutorial: Contrast Coding

Author

Reinhold Kliegl

Published

2026-06-21

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 AlgebraOfGraphics
using CairoMakie
using Chain
using CategoricalArrays
using DataFrames
using DataFrameMacros
using MixedModels
using ProgressMeter
using StableRNGs
using Statistics
using StatsBase
using MixedModels: likelihoodratiotest
using SMLP2026: dataset

const progress = isinteractive()

1.2 Readme for `dataset("fggk21")`

Number of scores: 525126

Cohort: 9 levels; 2011-2019
School: 515 levels
Child: 108295 levels; all children are between 8.0 and 8.99 years old
Sex: “Girls” (n=55,086), “Boys” (n= 53,209)
age: testdate - middle of month of birthdate
Test: 5 levels
- Endurance (Run): 6 minute endurance run [m]; to nearest 9m in 9x18m field
- Coordination (Star_r): star coordination run [m/s]; 9x9m field, 4 x diagonal = 50.912 m
- Speed(S20_r): 20-meters sprint [m/s]
- Muscle power low (SLJ): standing long jump [cm]
- Muscle power up (BPT): 1-kg medicine ball push test [m]
score - see units

1.3 Preprocessing

1.3.1 Read data

tbl = dataset(:fggk21)

Arrow.Table with 525126 rows, 7 columns, and schema:
 :Cohort  String
 :School  String
 :Child   String
 :Sex     String
 :age     Float64
 :Test    String
 :score   Float64

df = DataFrame(tbl)
describe(df)

7×7 DataFrame

Row	variable	mean	min	median	max	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Int64	DataType
1	Cohort		2011		2019	0	String
2	School		S100043		S800200	0	String
3	Child		C002352		C117966	0	String
4	Sex		female		male	0	String
5	age	8.56073	7.99452	8.55852	9.10609	0	Float64
6	Test		BPT		Star_r	0	String
7	score	226.141	1.14152	4.65116	1530.0	0	Float64

1.3.2 Extract a stratified subsample

We extract a random sample of 500 children from the Sex (2) x Test (5) cells of the design. Cohort and School are random.

dat = @chain df begin
  @transform(:Sex = :Sex == "female" ? "Girls" : "Boys")
  @groupby(:Test, :Sex)
  combine(x -> x[sample(StableRNG(32123), 1:nrow(x), 500), :])
end

5000×7 DataFrame

4975 rows omitted

Row	Test	Sex	Cohort	School	Child	age	score
	String	String	String	String	String	Float64	Float64
1	S20_r	Boys	2014	S104954	C045199	8.45722	5.0
2	S20_r	Boys	2019	S100833	C095158	8.88433	5.0
3	S20_r	Boys	2016	S104474	C082616	8.77207	5.26316
4	S20_r	Boys	2013	S101710	C012744	8.16427	4.87805
5	S20_r	Boys	2017	S105041	C038822	8.39425	4.7619
6	S20_r	Boys	2013	S101308	C006351	8.0794	4.65116
7	S20_r	Boys	2015	S102350	C083988	8.79671	4.34783
8	S20_r	Boys	2014	S100572	C092059	8.87064	4.87805
9	S20_r	Boys	2011	S112057	C061342	8.58042	4.87805
10	S20_r	Boys	2016	S105077	C037054	8.3833	3.38983
11	S20_r	Boys	2013	S104863	C021572	8.24641	4.34783
12	S20_r	Boys	2011	S101576	C079489	8.75017	4.08163
13	S20_r	Boys	2015	S104917	C094302	8.87885	4.44444
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
4989	Run	Girls	2016	S104255	C036639	8.37782	765.0
4990	Run	Girls	2018	S104930	C090981	8.85421	1314.0
4991	Run	Girls	2017	S102532	C048831	8.4846	1107.0
4992	Run	Girls	2016	S100304	C111778	9.04586	738.0
4993	Run	Girls	2014	S100122	C091947	8.87064	911.0
4994	Run	Girls	2015	S102465	C037011	8.3833	828.0
4995	Run	Girls	2018	S106409	C082278	8.76934	918.0
4996	Run	Girls	2015	S104930	C034971	8.3614	1053.0
4997	Run	Girls	2018	S113220	C102192	8.93908	990.0
4998	Run	Girls	2013	S103548	C087265	8.82683	910.0
4999	Run	Girls	2014	S103421	C035759	8.37509	972.0
5000	Run	Girls	2019	S100924	C019069	8.22177	1080.0

1.3.3 Transformations

transform!(dat, :age, :age => (x -> x .- 8.5) => :a1) # centered age (linear)
select!(groupby(dat, :Test), :, :score => zscore => :zScore) # z-score

5000×9 DataFrame

4975 rows omitted

Row	Test	Sex	Cohort	School	Child	age	score	a1	zScore
	String	String	String	String	String	Float64	Float64	Float64	Float64
1	S20_r	Boys	2014	S104954	C045199	8.45722	5.0	-0.0427789	1.1713
2	S20_r	Boys	2019	S100833	C095158	8.88433	5.0	0.384326	1.1713
3	S20_r	Boys	2016	S104474	C082616	8.77207	5.26316	0.272074	1.80951
4	S20_r	Boys	2013	S101710	C012744	8.16427	4.87805	-0.335729	0.875541
5	S20_r	Boys	2017	S105041	C038822	8.39425	4.7619	-0.105749	0.593867
6	S20_r	Boys	2013	S101308	C006351	8.0794	4.65116	-0.420602	0.325294
7	S20_r	Boys	2015	S102350	C083988	8.79671	4.34783	0.296715	-0.410363
8	S20_r	Boys	2014	S100572	C092059	8.87064	4.87805	0.370637	0.875541
9	S20_r	Boys	2011	S112057	C061342	8.58042	4.87805	0.0804244	0.875541
10	S20_r	Boys	2016	S105077	C037054	8.3833	3.38983	-0.116701	-2.73371
11	S20_r	Boys	2013	S104863	C021572	8.24641	4.34783	-0.253593	-0.410363
12	S20_r	Boys	2011	S101576	C079489	8.75017	4.08163	0.250171	-1.05594
13	S20_r	Boys	2015	S104917	C094302	8.87885	4.44444	0.37885	-0.176043
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
4989	Run	Girls	2016	S104255	C036639	8.37782	765.0	-0.122177	-1.49104
4990	Run	Girls	2018	S104930	C090981	8.85421	1314.0	0.354209	2.04567
4991	Run	Girls	2017	S102532	C048831	8.4846	1107.0	-0.0154004	0.712154
4992	Run	Girls	2016	S100304	C111778	9.04586	738.0	0.545859	-1.66498
4993	Run	Girls	2014	S100122	C091947	8.87064	911.0	0.370637	-0.550496
4994	Run	Girls	2015	S102465	C037011	8.3833	828.0	-0.116701	-1.08519
4995	Run	Girls	2018	S106409	C082278	8.76934	918.0	0.269336	-0.505401
4996	Run	Girls	2015	S104930	C034971	8.3614	1053.0	-0.138604	0.364281
4997	Run	Girls	2018	S113220	C102192	8.93908	990.0	0.439083	-0.0415708
4998	Run	Girls	2013	S103548	C087265	8.82683	910.0	0.326831	-0.556938
4999	Run	Girls	2014	S103421	C035759	8.37509	972.0	-0.124914	-0.157528
5000	Run	Girls	2019	S100924	C019069	8.22177	1080.0	-0.278234	0.538217

combine(
  groupby(dat, [:Test, :Sex]),
  :score => mean,
  :score => std,
  :zScore => mean,
  :zScore => std,
)

10×6 DataFrame

Row	Test	Sex	score_mean	score_std	zScore_mean	zScore_std
	String	String	Float64	Float64	Float64	Float64
1	S20_r	Boys	4.5841	0.418209	0.162654	1.01425
2	BPT	Boys	3.9752	0.722433	0.311603	0.99519
3	SLJ	Boys	129.996	18.466	0.199082	0.972295
4	Star_r	Boys	2.06142	0.302301	0.0727694	1.05774
5	Run	Boys	1036.57	166.222	0.258463	1.07082
6	S20_r	Girls	4.44997	0.39556	-0.162654	0.959317
7	BPT	Girls	3.5228	0.656156	-0.311603	0.903891
8	SLJ	Girls	122.434	18.7746	-0.199082	0.988544
9	Star_r	Girls	2.01982	0.266968	-0.0727694	0.934114
10	Run	Girls	956.332	131.851	-0.258463	0.849399

1.3.4 Figure of age x Sex x Test interactions

The main results of relevance here are shown in Figure 2 of Scientific Reports 11:17566.

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.0014	0.0143	0.10	0.9216	0.7593
Test: Star_r	-0.0044	0.0440	-0.10	0.9210
Test: S20_r	-0.0090	0.0442	-0.20	0.8385
Test: SLJ	-0.0062	0.0441	-0.14	0.8878
Test: BPT	0.0082	0.0440	0.19	0.8529
Residual	0.6512

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.0077	0.0144	-0.53	0.5960	0.7538
Test: Star_r	-0.0182	0.0442	-0.41	0.6814
Test: S20_r	-0.0019	0.0445	-0.04	0.9655
Test: SLJ	0.0010	0.0444	0.02	0.9820
Test: BPT	-0.0091	0.0442	-0.21	0.8373
a1	0.2029	0.0500	4.06	<1e-04
Test: Star_r & a1	0.3822	0.1519	2.52	0.0119
Test: S20_r & a1	-0.1670	0.1530	-1.09	0.2753
Test: SLJ & a1	-0.2320	0.1519	-1.53	0.1266
Test: BPT & a1	0.5125	0.1506	3.40	0.0007
Residual	0.6523

The difference between older and younger children is larger for Star_r than for Run (0.2473). S20_r did not differ significantly from Star_r (-0.0377) and SLJ (-0.0113) The largest difference in developmental gain was between BPT and SLJ (0.3355).

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.0073	0.0141	-0.52	0.6033	0.7310
Test: Star_r	-0.0216	0.0432	-0.50	0.6177
Test: S20_r	0.0015	0.0435	0.04	0.9720
Test: SLJ	0.0006	0.0433	0.01	0.9887
Test: BPT	-0.0089	0.0432	-0.21	0.8362
a1	0.1853	0.0488	3.80	0.0001
Sex: Boys	0.2004	0.0139	14.39	<1e-46
Test: Star_r & a1	0.4133	0.1486	2.78	0.0054
Test: S20_r & a1	-0.1919	0.1496	-1.28	0.1998
Test: SLJ & a1	-0.2320	0.1484	-1.56	0.1180
Test: BPT & a1	0.5010	0.1473	3.40	0.0007
Test: Star_r & Sex: Boys	-0.1857	0.0429	-4.33	<1e-04
Test: S20_r & Sex: Boys	0.1003	0.0431	2.33	0.0198
Test: SLJ & Sex: Boys	0.0319	0.0430	0.74	0.4582
Test: BPT & Sex: Boys	0.1110	0.0429	2.59	0.0097
Residual	0.6424

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.0069	0.0141	-0.49	0.6233	0.7283
Test: Star_r	-0.0261	0.0432	-0.60	0.5458
Test: S20_r	0.0054	0.0435	0.12	0.9018
Test: SLJ	0.0017	0.0433	0.04	0.9687
Test: BPT	-0.0104	0.0432	-0.24	0.8096
a1	0.1850	0.0488	3.79	0.0001
Sex: Boys	0.2026	0.0141	14.38	<1e-46
Test: Star_r & a1	0.4198	0.1486	2.83	0.0047
Test: S20_r & a1	-0.1955	0.1495	-1.31	0.1911
Test: SLJ & a1	-0.2317	0.1484	-1.56	0.1183
Test: BPT & a1	0.5024	0.1472	3.41	0.0006
Test: Star_r & Sex: Boys	-0.1969	0.0432	-4.55	<1e-05
Test: S20_r & Sex: Boys	0.1093	0.0435	2.51	0.0119
Test: SLJ & Sex: Boys	0.0390	0.0433	0.90	0.3680
Test: BPT & Sex: Boys	0.1078	0.0432	2.50	0.0125
a1 & Sex: Boys	-0.0457	0.0488	-0.94	0.3484
Test: Star_r & a1 & Sex: Boys	0.2786	0.1486	1.87	0.0609
Test: S20_r & a1 & Sex: Boys	-0.2118	0.1495	-1.42	0.1566
Test: SLJ & a1 & Sex: Boys	-0.1757	0.1484	-1.18	0.2364
Test: BPT & a1 & Sex: Boys	0.0865	0.1472	0.59	0.5570
Residual	0.6441

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.0069	0.0141	-0.49	0.6233	0.7283
Test: Star_r	-0.0131	0.0216	-0.60	0.5458
Test: S20_r	-0.0026	0.0126	-0.20	0.8381
Test: SLJ	-0.0009	0.0088	-0.10	0.9227
Test: BPT	-0.0026	0.0069	-0.38	0.7051
a1	0.1850	0.0488	3.79	0.0001
Sex: Boys	0.2026	0.0141	14.38	<1e-46
Test: Star_r & a1	0.2099	0.0743	2.83	0.0047
Test: S20_r & a1	0.0048	0.0432	0.11	0.9116
Test: SLJ & a1	-0.0555	0.0303	-1.83	0.0669
Test: BPT & a1	0.0672	0.0235	2.86	0.0042
Test: Star_r & Sex: Boys	-0.0985	0.0216	-4.55	<1e-05
Test: S20_r & Sex: Boys	0.0036	0.0126	0.29	0.7735
Test: SLJ & Sex: Boys	0.0116	0.0088	1.31	0.1910
Test: BPT & Sex: Boys	0.0285	0.0069	4.16	<1e-04
a1 & Sex: Boys	-0.0457	0.0488	-0.94	0.3484
Test: Star_r & a1 & Sex: Boys	0.1393	0.0743	1.87	0.0609
Test: S20_r & a1 & Sex: Boys	-0.0242	0.0432	-0.56	0.5760
Test: SLJ & a1 & Sex: Boys	-0.0560	0.0303	-1.85	0.0646
Test: BPT & a1 & Sex: Boys	-0.0163	0.0235	-0.70	0.4869
Residual	0.6441

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.0069	0.0141	-0.49	0.6233	0.7283
Test: BPT-other	-0.0519	0.1371	-0.38	0.7051
Test: Star-End	-0.0261	0.0432	-0.60	0.5458
Test: S20-Star	0.0054	0.0435	0.12	0.9018
Test: SLJ-S20	0.0017	0.0433	0.04	0.9687
a1	0.1850	0.0488	3.79	0.0001
Sex: Boys	0.2026	0.0141	14.38	<1e-46
Test: BPT-other & a1	1.3433	0.4692	2.86	0.0042
Test: Star-End & a1	0.4198	0.1486	2.83	0.0047
Test: S20-Star & a1	-0.1955	0.1495	-1.31	0.1911
Test: SLJ-S20 & a1	-0.2317	0.1484	-1.56	0.1183
Test: BPT-other & Sex: Boys	0.5699	0.1371	4.16	<1e-04
Test: Star-End & Sex: Boys	-0.1969	0.0432	-4.55	<1e-05
Test: S20-Star & Sex: Boys	0.1093	0.0435	2.51	0.0119
Test: SLJ-S20 & Sex: Boys	0.0390	0.0433	0.90	0.3680
a1 & Sex: Boys	-0.0457	0.0488	-0.94	0.3484
Test: BPT-other & a1 & Sex: Boys	-0.3262	0.4692	-0.70	0.4869
Test: Star-End & a1 & Sex: Boys	0.2786	0.1486	1.87	0.0609
Test: S20-Star & a1 & Sex: Boys	-0.2118	0.1495	-1.42	0.1566
Test: SLJ-S20 & a1 & Sex: Boys	-0.1757	0.1484	-1.18	0.2364
Residual	0.6441

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.0069	0.0141	-0.49	0.6233	0.7283
Test: Star-Run	-0.0261	0.0432	-0.60	0.5458
Test: S20-Star	0.0054	0.0435	0.12	0.9018
Test: SLJ-S20	0.0017	0.0433	0.04	0.9687
Test: BPT-SLJ	-0.0104	0.0432	-0.24	0.8096
a1	0.1850	0.0488	3.79	0.0001
Sex: Boys	0.2026	0.0141	14.38	<1e-46
Test: Star-Run & a1	0.4198	0.1486	2.83	0.0047
Test: S20-Star & a1	-0.1955	0.1495	-1.31	0.1911
Test: SLJ-S20 & a1	-0.2317	0.1484	-1.56	0.1183
Test: BPT-SLJ & a1	0.5024	0.1472	3.41	0.0006
Test: Star-Run & Sex: Boys	-0.1969	0.0432	-4.55	<1e-05
Test: S20-Star & Sex: Boys	0.1093	0.0435	2.51	0.0119
Test: SLJ-S20 & Sex: Boys	0.0390	0.0433	0.90	0.3680
Test: BPT-SLJ & Sex: Boys	0.1078	0.0432	2.50	0.0125
a1 & Sex: Boys	-0.0457	0.0488	-0.94	0.3484
Test: Star-Run & a1 & Sex: Boys	0.2786	0.1486	1.87	0.0609
Test: S20-Star & a1 & Sex: Boys	-0.2118	0.1495	-1.42	0.1566
Test: SLJ-S20 & a1 & Sex: Boys	-0.1757	0.1484	-1.18	0.2364
Test: BPT-SLJ & a1 & Sex: Boys	0.0865	0.1472	0.59	0.5570
Residual	0.6441

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.0065	0.0141	-0.46	0.6432	0.7248
Test: Star-Run	-0.0239	0.0434	-0.55	0.5815	0.0000
Test: S20-Star	0.0063	0.0433	0.15	0.8843	0.5276
Test: SLJ-S20	-0.0024	0.0429	-0.06	0.9548	0.4473
Test: BPT-SLJ	0.0007	0.0430	0.02	0.9865	0.0941
a1	0.1852	0.0488	3.79	0.0001
Sex: Boys	0.2027	0.0141	14.38	<1e-46
Test: Star-Run & a1	0.4300	0.1487	2.89	0.0038
Test: S20-Star & a1	-0.2070	0.1488	-1.39	0.1642
Test: SLJ-S20 & a1	-0.2208	0.1468	-1.50	0.1326
Test: BPT-SLJ & a1	0.4694	0.1459	3.22	0.0013
Test: Star-Run & Sex: Boys	-0.1960	0.0434	-4.51	<1e-05
Test: S20-Star & Sex: Boys	0.1059	0.0433	2.45	0.0144
Test: SLJ-S20 & Sex: Boys	0.0332	0.0429	0.77	0.4398
Test: BPT-SLJ & Sex: Boys	0.1164	0.0430	2.71	0.0068
a1 & Sex: Boys	-0.0454	0.0488	-0.93	0.3522
Test: Star-Run & a1 & Sex: Boys	0.2626	0.1487	1.77	0.0774
Test: S20-Star & a1 & Sex: Boys	-0.2006	0.1488	-1.35	0.1776
Test: SLJ-S20 & a1 & Sex: Boys	-0.1702	0.1468	-1.16	0.2462
Test: BPT-SLJ & a1 & Sex: Boys	0.0834	0.1459	0.57	0.5678
Residual	0.5521

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.0062	0.0140	-0.45	0.6559
Test: Star-Run	-0.0245	0.0437	-0.56	0.5748
Test: S20-Star	0.0048	0.0443	0.11	0.9137
Test: SLJ-S20	0.0131	0.0443	0.30	0.7674
Test: BPT-SLJ	-0.0287	0.0432	-0.66	0.5073
a1	0.1844	0.0483	3.82	0.0001
Sex: Boys	0.2018	0.0140	14.44	<1e-46
Test: Star-Run & a1	0.4064	0.1507	2.70	0.0070
Test: S20-Star & a1	-0.1807	0.1527	-1.18	0.2366
Test: SLJ-S20 & a1	-0.2782	0.1524	-1.82	0.0680
Test: BPT-SLJ & a1	0.5586	0.1489	3.75	0.0002
Test: Star-Run & Sex: Boys	-0.1971	0.0437	-4.51	<1e-05
Test: S20-Star & Sex: Boys	0.1014	0.0443	2.29	0.0220
Test: SLJ-S20 & Sex: Boys	0.0451	0.0443	1.02	0.3089
Test: BPT-SLJ & Sex: Boys	0.1019	0.0432	2.36	0.0185
a1 & Sex: Boys	-0.0448	0.0483	-0.93	0.3536
Test: Star-Run & a1 & Sex: Boys	0.3097	0.1507	2.06	0.0398
Test: S20-Star & a1 & Sex: Boys	-0.2051	0.1527	-1.34	0.1792
Test: SLJ-S20 & a1 & Sex: Boys	-0.2030	0.1524	-1.33	0.1831
Test: BPT-SLJ & a1 & Sex: Boys	0.0979	0.1489	0.66	0.5107
Test: Run					0.9577
Test: Star_r					0.9839
Test: S20_r					0.9855
Test: SLJ					0.9716
Test: BPT					0.9394
Residual	0.0001

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.0062	0.0140	-0.45	0.6559
Test: Star-Run	-0.0245	0.0437	-0.56	0.5748
Test: S20-Star	0.0048	0.0443	0.11	0.9137
Test: SLJ-S20	0.0131	0.0443	0.30	0.7674
Test: BPT-SLJ	-0.0287	0.0432	-0.66	0.5073
a1	0.1844	0.0483	3.82	0.0001
Sex: Boys	0.2018	0.0140	14.44	<1e-46
Test: Star-Run & a1	0.4064	0.1507	2.70	0.0070
Test: S20-Star & a1	-0.1807	0.1527	-1.18	0.2366
Test: SLJ-S20 & a1	-0.2782	0.1524	-1.82	0.0680
Test: BPT-SLJ & a1	0.5586	0.1489	3.75	0.0002
Test: Star-Run & Sex: Boys	-0.1971	0.0437	-4.51	<1e-05
Test: S20-Star & Sex: Boys	0.1014	0.0443	2.29	0.0220
Test: SLJ-S20 & Sex: Boys	0.0451	0.0443	1.02	0.3089
Test: BPT-SLJ & Sex: Boys	0.1019	0.0432	2.36	0.0185
a1 & Sex: Boys	-0.0448	0.0483	-0.93	0.3536
Test: Star-Run & a1 & Sex: Boys	0.3097	0.1507	2.06	0.0398
Test: S20-Star & a1 & Sex: Boys	-0.2051	0.1527	-1.34	0.1792
Test: SLJ-S20 & a1 & Sex: Boys	-0.2030	0.1524	-1.33	0.1831
Test: BPT-SLJ & a1 & Sex: Boys	0.0979	0.1489	0.66	0.5107
Test: Run					0.9577
Test: Star_r					0.9839
Test: S20_r					0.9855
Test: SLJ					0.9716
Test: BPT					0.9394
Residual	0.0001

VarCorr(m_cpx_0_SeqDiff)

	Column	Variance	Std.Dev	Corr.
Child	Test: Run	0.91718958	0.95770015
	Test: Star_r	0.96803107	0.98388570	+0.28
	Test: S20_r	0.97111832	0.98545336	+0.55	+0.51
	Test: SLJ	0.94393655	0.97156397	+0.41	+0.32	-0.05
	Test: BPT	0.88244094	0.93938328	+0.25	+0.67	+0.43	-0.02
Residual		0.00000001	0.00010119

m_cpx_0_SeqDiff.PCA

(Child = 
Principal components based on correlation matrix
 Test: Run      1.0     .      .      .      .
 Test: Star_r   0.28   1.0     .      .      .
 Test: S20_r    0.55   0.51   1.0     .      .
 Test: SLJ      0.41   0.32  -0.05   1.0     .
 Test: BPT      0.25   0.67   0.43  -0.02   1.0

Normalized cumulative variances:
[0.4867, 0.7246, 0.8955, 0.9725, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5
 Test: Run     -0.44   0.39   0.54  -0.44  -0.41
 Test: Star_r  -0.54  -0.09  -0.46   0.42  -0.56
 Test: S20_r   -0.49  -0.23   0.51   0.49   0.45
 Test: SLJ     -0.22   0.79  -0.34   0.1    0.44
 Test: BPT     -0.47  -0.4   -0.35  -0.62   0.34,)

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.0061	0.0140	-0.43	0.6657	0.7230
Test: Star-Run	-0.0236	0.0431	-0.55	0.5839	0.8364
Test: S20-Star	0.0039	0.0441	0.09	0.9299	0.8995
Test: SLJ-S20	0.0091	0.0438	0.21	0.8355	0.8356
Test: BPT-SLJ	-0.0174	0.0429	-0.41	0.6846	1.0836
a1	0.1853	0.0486	3.81	0.0001
Sex: Boys	0.2027	0.0140	14.43	<1e-46
Test: Star-Run & a1	0.4134	0.1470	2.81	0.0049
Test: S20-Star & a1	-0.1890	0.1516	-1.25	0.2124
Test: SLJ-S20 & a1	-0.2609	0.1501	-1.74	0.0822
Test: BPT-SLJ & a1	0.5195	0.1467	3.54	0.0004
Test: Star-Run & Sex: Boys	-0.1987	0.0431	-4.61	<1e-05
Test: S20-Star & Sex: Boys	0.1104	0.0441	2.50	0.0124
Test: SLJ-S20 & Sex: Boys	0.0259	0.0438	0.59	0.5540
Test: BPT-SLJ & Sex: Boys	0.1143	0.0429	2.67	0.0077
a1 & Sex: Boys	-0.0445	0.0486	-0.91	0.3604
Test: Star-Run & a1 & Sex: Boys	0.2960	0.1470	2.01	0.0440
Test: S20-Star & a1 & Sex: Boys	-0.2087	0.1516	-1.38	0.1686
Test: SLJ-S20 & a1 & Sex: Boys	-0.1540	0.1501	-1.03	0.3050
Test: BPT-SLJ & a1 & Sex: Boys	0.0603	0.1467	0.41	0.6810
Residual	0.0001

m_cpx_1_SeqDiff.PCA

(Child = 
Principal components based on correlation matrix
 (Intercept)      1.0     .      .      .      .
 Test: Star-Run   0.17   1.0     .      .      .
 Test: S20-Star   0.06  -0.6    1.0     .      .
 Test: SLJ-S20   -0.22   0.11  -0.12   1.0     .
 Test: BPT-SLJ   -0.12   0.61  -0.54  -0.36   1.0

Normalized cumulative variances:
[0.4343, 0.6908, 0.8911, 0.9703, 1.0]

Component loadings
                   PC1    PC2    PC3    PC4    PC5
 (Intercept)     -0.01   0.54   0.78   0.12   0.3
 Test: Star-Run  -0.58  -0.06   0.32  -0.53  -0.52
 Test: S20-Star   0.56   0.23  -0.07  -0.79   0.1
 Test: SLJ-S20    0.06  -0.78   0.41  -0.17   0.44
 Test: BPT-SLJ   -0.58   0.19  -0.36  -0.24   0.66,)

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.0060	0.0141	-0.43	0.6686	0.7446
Test: c5.1	-0.0285	0.0427	-0.67	0.5042	1.1465
Test: c234.15	-0.0023	0.1647	-0.01	0.9887	0.2807
Test: c2.34	-0.0159	0.0759	-0.21	0.8338	1.5034
Test: c3.4	-0.0112	0.0441	-0.25	0.7989	1.3408
a1	0.1838	0.0486	3.78	0.0002
Sex: Boys	0.2025	0.0141	14.41	<1e-46
Test: c5.1 & a1	0.4904	0.1470	3.33	0.0009
Test: c234.15 & a1	-0.3516	0.5587	-0.63	0.5291
Test: c2.34 & a1	0.6574	0.2605	2.52	0.0116
Test: c3.4 & a1	0.2749	0.1517	1.81	0.0700
Test: c5.1 & Sex: Boys	0.0547	0.0427	1.28	0.2000
Test: c234.15 & Sex: Boys	-0.8430	0.1647	-5.12	<1e-06
Test: c2.34 & Sex: Boys	-0.2469	0.0759	-3.25	0.0011
Test: c3.4 & Sex: Boys	-0.0451	0.0441	-1.02	0.3075
a1 & Sex: Boys	-0.0467	0.0486	-0.96	0.3365
Test: c5.1 & a1 & Sex: Boys	-0.0078	0.1470	-0.05	0.9575
Test: c234.15 & a1 & Sex: Boys	0.5470	0.5587	0.98	0.3275
Test: c2.34 & a1 & Sex: Boys	0.5878	0.2605	2.26	0.0240
Test: c3.4 & a1 & Sex: Boys	0.1878	0.1517	1.24	0.2159
Residual	0.0001

VarCorr(m_cpx_1_PC)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.554445410	0.744610912
	Test: c5.1	1.314478795	1.146507216	-0.04
	Test: c234.15	0.078781275	0.280680023	-0.61	-0.39
	Test: c2.34	2.260246069	1.503411477	+0.25	+0.10	-0.33
	Test: c3.4	1.797793071	1.340818060	+0.08	+0.17	-0.69	+0.13
Residual		0.000000012	0.000108827

m_cpx_1_PC.PCA

(Child = 
Principal components based on correlation matrix
 (Intercept)     1.0     .      .      .      .
 Test: c5.1     -0.04   1.0     .      .      .
 Test: c234.15  -0.61  -0.39   1.0     .      .
 Test: c2.34     0.25   0.1   -0.33   1.0     .
 Test: c3.4      0.08   0.17  -0.69   0.13   1.0

Normalized cumulative variances:
[0.4455, 0.6656, 0.8394, 0.9872, 1.0]

Component loadings
                  PC1    PC2    PC3    PC4    PC5
 (Intercept)    -0.41   0.61  -0.1    0.5   -0.44
 Test: c5.1     -0.28  -0.63   0.52   0.44  -0.23
 Test: c234.15   0.65   0.03   0.15  -0.11  -0.74
 Test: c2.34    -0.33   0.32   0.67  -0.57  -0.06
 Test: c3.4     -0.47  -0.35  -0.49  -0.46  -0.46,)

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.0055	0.0142	-0.39	0.6971	0.7581
Test: c5.1	0.0331	0.0301	1.10	0.2709	0.6112
Test: c234.15	0.0109	0.0302	0.36	0.7179	0.5295
Test: c12.34	0.0040	0.0308	0.13	0.8965	0.7168
Test: c3.4	-0.0121	0.0314	-0.38	0.7003	0.8100
a1	0.1704	0.0491	3.47	0.0005
Sex: Boys	0.1954	0.0142	13.75	<1e-42
Test: c5.1 & a1	-0.4004	0.1031	-3.89	0.0001
Test: c234.15 & a1	-0.1340	0.1031	-1.30	0.1939
Test: c12.34 & a1	0.1402	0.1062	1.32	0.1868
Test: c3.4 & a1	0.2557	0.1079	2.37	0.0178
Test: c5.1 & Sex: Boys	-0.0686	0.0301	-2.28	0.0226
Test: c234.15 & Sex: Boys	0.1630	0.0302	5.39	<1e-07
Test: c12.34 & Sex: Boys	-0.0576	0.0308	-1.87	0.0617
Test: c3.4 & Sex: Boys	-0.0159	0.0314	-0.51	0.6133
a1 & Sex: Boys	-0.0413	0.0491	-0.84	0.4010
Test: c5.1 & a1 & Sex: Boys	0.0392	0.1031	0.38	0.7037
Test: c234.15 & a1 & Sex: Boys	-0.1682	0.1031	-1.63	0.1030
Test: c12.34 & a1 & Sex: Boys	0.1777	0.1062	1.67	0.0942
Test: c3.4 & a1 & Sex: Boys	0.1415	0.1079	1.31	0.1899
Residual	0.0000

VarCorr(m_cpx_1_PC_2)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.57466736	0.75806818
	Test: c5.1	0.37351796	0.61116116	+0.10
	Test: c234.15	0.28035780	0.52948825	+0.04	+0.01
	Test: c12.34	0.51374850	0.71676251	+0.09	-0.57	+0.27
	Test: c3.4	0.65604892	0.80996847	+0.24	-0.20	+0.18	+0.46
Residual		0.00000000	0.00000692

m_cpx_1_PC_2.PCA

(Child = 
Principal components based on correlation matrix
 (Intercept)     1.0     .      .      .      .
 Test: c5.1      0.1    1.0     .      .      .
 Test: c234.15   0.04   0.01   1.0     .      .
 Test: c12.34    0.09  -0.57   0.27   1.0     .
 Test: c3.4      0.24  -0.2    0.18   0.46   1.0

Normalized cumulative variances:
[0.3866, 0.6217, 0.8113, 0.9351, 1.0]

Component loadings
                  PC1    PC2    PC3    PC4    PC5
 (Intercept)    -0.15  -0.72   0.44   0.51   0.08
 Test: c5.1      0.48  -0.54  -0.21  -0.35  -0.56
 Test: c234.15  -0.28  -0.29  -0.86   0.22   0.23
 Test: c12.34   -0.64   0.14  -0.0    0.12  -0.75
 Test: c3.4     -0.51  -0.3    0.15  -0.75   0.26,)

# 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.0059	0.0143	-0.41	0.6822	0.7290
Test: c5.1	0.0239	0.0291	0.82	0.4114	0.4594
Test: c234.15	0.0134	0.0314	0.43	0.6684	0.8128
Test: c12.34	0.0004	0.0318	0.01	0.9890	0.7930
Test: c3.4	-0.0046	0.0314	-0.15	0.8831	0.7566
a1	0.1709	0.0494	3.46	0.0005
Sex: Boys	0.1952	0.0143	13.66	<1e-41
Test: c5.1 & a1	-0.3756	0.0997	-3.77	0.0002
Test: c234.15 & a1	-0.1347	0.1081	-1.25	0.2128
Test: c12.34 & a1	0.1509	0.1099	1.37	0.1695
Test: c3.4 & a1	0.2364	0.1076	2.20	0.0280
Test: c5.1 & Sex: Boys	-0.0739	0.0291	-2.54	0.0110
Test: c234.15 & Sex: Boys	0.1678	0.0314	5.34	<1e-07
Test: c12.34 & Sex: Boys	-0.0571	0.0318	-1.79	0.0728
Test: c3.4 & Sex: Boys	-0.0181	0.0314	-0.58	0.5649
a1 & Sex: Boys	-0.0440	0.0494	-0.89	0.3734
Test: c5.1 & a1 & Sex: Boys	0.0349	0.0997	0.35	0.7266
Test: c234.15 & a1 & Sex: Boys	-0.1766	0.1081	-1.63	0.1023
Test: c12.34 & a1 & Sex: Boys	0.1888	0.1099	1.72	0.0857
Test: c3.4 & a1 & Sex: Boys	0.1440	0.1076	1.34	0.1806
Residual	0.0000

VarCorr(m_zcp_1_PC_2)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.53143023	0.72899261
	Test: c5.1	0.21104625	0.45939771	.
	Test: c234.15	0.66071258	0.81284229	.	.
	Test: c12.34	0.62892085	0.79304530	.	.	.
	Test: c3.4	0.57244504	0.75660098	.	.	.	.
Residual		0.00000000	0.00000703

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	-12940
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-12928	12	10	0.2856

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}:
 13847.543243300042
 13847.54324330003
 13847.543243299973

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.4815585	0.6939441
	Test: Star_r	0.4225459	0.6500353	+0.21
	Test: S20_r	0.6091305	0.7804681	-0.30	-0.13
	Test: SLJ	1.2739964	1.1287145	+0.02	+0.16	-0.57
	Test: BPT	1.0858291	1.0420313	-0.22	-0.47	-0.21	-0.08
Residual		0.0000000	0.0000175

VarCorr(m_Child_HeC)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.23132450	0.48096205
	Test: Star_r	0.38241758	0.61839921	+0.03
	Test: S20_r	0.14694408	0.38333286	-0.10	+0.15
	Test: SLJ	0.05866651	0.24221171	-0.29	+0.07	+0.60
	Test: BPT	0.06296274	0.25092378	-0.20	-0.27	-0.60	-0.72
Residual		0.00000000	0.00000037

VarCorr(m_Child_HyC)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.45648953	0.67564009
	Test: BPT-other	1.43810923	1.19921192	+0.98
	Test: Star-End	2.20582996	1.48520368	+0.03	+0.11
	Test: S20-Star	1.38484315	1.17679359	+0.21	+0.04	-0.78
	Test: SLJ-S20	1.72665474	1.31402235	-0.20	+0.00	+0.26	-0.55
Residual		0.00000000	0.00000036

VarCorr(m_Child_PCA)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.54829777	0.74047131
	Test: c5.1	0.39609331	0.62935944	-0.24
	Test: c234.15	3.90793528	1.97684984	-0.59	-0.51
	Test: c2.34	3.45080286	1.85763367	+0.12	-0.64	+0.56
	Test: c3.4	2.22678560	1.49224180	+0.09	+0.22	+0.19	+0.21
Residual		0.00000000	0.00000032

The CPs for the various contrasts are in line with expectations. For the SDC we observe substantial negative CPs between 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. The (roughly) PCA-based contrasts are small with one exception; there is a sizeable CP of +.41 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.0060	0.0141	-0.42	0.6722	0.7514
Test: c5.1	-0.0221	0.0429	-0.52	0.6056	1.1831
Test: c234.15	-0.0434	0.1649	-0.26	0.7926	0.0299
Test: c2.34	-0.0100	0.0768	-0.13	0.8966	1.9567
Test: c3.4	-0.0053	0.0435	-0.12	0.9035	1.0495
a1	0.1843	0.0487	3.78	0.0002
Sex: Boys	0.2027	0.0141	14.41	<1e-46
Test: c5.1 & a1	0.4655	0.1477	3.15	0.0016
Test: c234.15 & a1	-0.2459	0.5596	-0.44	0.6603
Test: c2.34 & a1	0.6414	0.2642	2.43	0.0152
Test: c3.4 & a1	0.2499	0.1491	1.68	0.0937
Test: c5.1 & Sex: Boys	0.0583	0.0429	1.36	0.1743
Test: c234.15 & Sex: Boys	-0.8636	0.1649	-5.24	<1e-06
Test: c2.34 & Sex: Boys	-0.2425	0.0768	-3.16	0.0016
Test: c3.4 & Sex: Boys	-0.0375	0.0435	-0.86	0.3885
a1 & Sex: Boys	-0.0470	0.0487	-0.96	0.3348
Test: c5.1 & a1 & Sex: Boys	-0.0028	0.1477	-0.02	0.9848
Test: c234.15 & a1 & Sex: Boys	0.6425	0.5596	1.15	0.2509
Test: c2.34 & a1 & Sex: Boys	0.5886	0.2642	2.23	0.0259
Test: c3.4 & a1 & Sex: Boys	0.1689	0.1491	1.13	0.2573
Residual	0.0000

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	-13842
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-13008	835	10	<1e-99

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	-13018
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-12760	258	10	<1e-48

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	-13097
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-12707	390	10	<1e-76

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	-12958
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 + Test \| Child)	36	-12702	255	10	<1e-48

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	12957.5	13009.5	13179.0
2	HeC0	26	13017.8	13069.8	13239.3
3	HyC0	26	13097.5	13149.5	13318.9
4	SDC0	26	13842.5	13894.5	14063.9

The best fit was obtained for the PCA-based zcpLMM. Somewhat surprisingly the second best fit was obtained for the SDC. The relatively poor performance of HeC-based zcpLMM is puzzling to me. I thought it might be related to imbalance in design in the present data, but this does not appear to be the case. The same comparison of sequential difference and Helmert coding also showed a worse fit for the zcp-HeC LMM than the zcp-SDC LMM.

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.046804	0.216343
	Test: Star_r	0.160446	0.400557	+0.05
	Test: S20_r	0.083645	0.289214	-0.09	-0.65
	Test: SLJ	0.056538	0.237778	-0.23	-0.16	-0.31
	Test: BPT	0.085374	0.292188	-0.29	-0.24	+0.13	-0.17
	a1	0.130157	0.360773	+0.40	+0.57	-0.38	+0.36	-0.29
Residual		0.842496	0.917876

VarCorr(m_School_Helmert)

	Column	Variance	Std.Dev	Corr.
School	(Intercept)	0.0468016	0.2163367
	Test: Star_r	0.0400951	0.2002377	+0.05
	Test: S20_r	0.0053896	0.0734139	-0.07	+0.06
	Test: SLJ	0.0024695	0.0496939	-0.33	-0.15	+0.08
	Test: BPT	0.0034667	0.0588790	-0.45	-0.32	-0.02	+0.27
	a1	0.1302342	0.3608798	+0.40	+0.57	+0.02	+0.44	-0.07
Residual		0.8424998	0.9178779

VarCorr(m_School_Hypo)

	Column	Variance	Std.Dev	Corr.
School	(Intercept)	0.046802	0.216338
	Test: BPT-other	1.387336	1.177852	-0.45
	Test: Star-End	0.160443	0.400553	+0.05	-0.32
	Test: S20-Star	0.083685	0.289284	-0.09	+0.21	-0.65
	Test: SLJ-S20	0.056562	0.237829	-0.23	+0.23	-0.16	-0.31
	a1	0.130139	0.360747	+0.40	-0.07	+0.57	-0.38	+0.36
Residual		0.842497	0.917876

VarCorr(m_School_PCA)

	Column	Variance	Std.Dev	Corr.
School	(Intercept)	0.046809	0.216354
	Test: c5.1	0.103784	0.322155	-0.45
	Test: c234.15	2.410929	1.552717	+0.22	-0.11
	Test: c2.34	0.304914	0.552190	+0.19	-0.02	+0.46
	Test: c3.4	0.056572	0.237849	+0.23	-0.11	+0.24	+0.10
	a1	0.130164	0.360782	+0.40	+0.36	+0.48	+0.24	-0.36
Residual		0.842522	0.917890

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	13821.6	13875.6	14051.5
2	PCA0	27	13822.6	13876.6	14052.5
3	HyC0	27	13827.8	13881.8	14057.8
4	SDC0	27	13828.9	13882.9	14058.8

For the random factor School the Helmert contrast, followed by PCA-based contrasts have least information in the CPs; SDC has the largest contribution from 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 e563d55 using Quarto 1.9.38 and Julia 1.12.6.

1 Setup

1.1 Packages and functions

1.2 Readme for dataset("fggk21")

1.3 Preprocessing

1.3.1 Read data

1.3.2 Extract a stratified subsample

1.3.3 Transformations

1.3.4 Figure of age x Sex x Test interactions

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

1.2 Readme for `dataset("fggk21")`

2.1 Sequential difference coding: `contr1`

2.2 _Helmert coding: `contr2`

2.3 Hypothesis coding: `contr3`

2.4 PCA-based hypothesis coding: `contr4`