The Emotikon Project: Model Complexification and PCA

Authors

Reinhold Kliegl

Phillip Alday

Published

2026-08-28

Learning objectives — this capstone pulls the course together on a large dataset. After working through it you will be able to:

carry out a complete analysis: data preparation, plotting the main results, contrast specification, and model complexification;
use PCA of the random effects to guide the random-effects structure;
interpret a parsimonious model for the Emotikon fitness data.

Before you start

Prerequisites: Mixed Models Tutorial: Contrast Coding; Creating multi-panel plots.

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

This script uses a subset of data reported in Fühner et al. (2021). To circumvent delays associated with model fitting we work with models that are less complex than those in the reference publication. All the data to reproduce the models in the publication are used here, too; the script requires only a few changes to specify the more complex models in the article.

The script is structured in three main sections:

Setup with reading and examining the data, plotting the main results, and specifying the contrasts for the fixed factor Test
a demonstration of model complexification to determine a parsimonious random-effect structure appropriate for and supported by the data, including also a quite elaborate demonstration of principle component analyses (PCAs) of levels (scores) and effects,
specification of nested fixed effects or interactions in the levels of another, superordinate factors.

1 Packages and functions

Code

using AlgebraOfGraphics
using AlgebraOfGraphics: linear
using CairoMakie
using CategoricalArrays
using DataFrames
using MixedModels
using MixedModelsMakie
using MixedModelsMakie: simplelinreg
using Statistics
using StatsBase
using MixedModelsDatasets: dataset
using SMLP2026: fit_or_restore

const progress = isinteractive()

false

2 Preprocessing

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

2.1 Read data

df = DataFrame(dataset(:fggk21))
transform!(df,
    :age => (x -> x .- 8.5) => :a1,
    :Sex => categorical => :Sex,
    :Test => categorical => :Test,
  )
levels!(df.Sex, ["male", "female"])
recode!(df.Sex, "male" => "Boys", "female" => "Girls")
levels!(df.Test, ["Run", "Star_r", "S20_r", "SLJ", "BPT"])
recode!(
  df.Test,
  "Run" => "Endurance",
  "Star_r" => "Coordination",
  "S20_r" => "Speed",
  "SLJ" => "PowerLOW",
  "BPT" => "PowerUP",
)
describe(df)

8×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		Boys		Girls	0	CategoricalValue{String, UInt32}
5	age	8.56073	7.99452	8.55852	9.10609	0	Float64
6	Test		Endurance		PowerUP	0	CategoricalValue{String, UInt32}
7	score	226.141	1.14152	4.65116	1530.0	0	Float64
8	a1	0.0607297	-0.505476	0.0585216	0.606092	0	Float64

2.1.1 Transformations

We center age at 8.5 years and compute z-scores for each Test. With these variables the data frame df contains all variables used for the final model in the original publication.

select!(groupby(df, :Test),  Not(:score), :score => zscore => :zScore)

525126×8 DataFrame

525101 rows omitted

Row	Test	Cohort	School	Child	Sex	age	a1	zScore
	Cat…	String	String	String	Cat…	Float64	Float64	Float64
1	Speed	2013	S100067	C002352	Boys	7.99452	-0.505476	1.7913
2	PowerUP	2013	S100067	C002352	Boys	7.99452	-0.505476	-0.0622317
3	PowerLOW	2013	S100067	C002352	Boys	7.99452	-0.505476	-0.0336567
4	Coordination	2013	S100067	C002352	Boys	7.99452	-0.505476	1.46874
5	Endurance	2013	S100067	C002352	Boys	7.99452	-0.505476	0.331058
6	Speed	2013	S100067	C002353	Boys	7.99452	-0.505476	1.15471
7	PowerUP	2013	S100067	C002353	Boys	7.99452	-0.505476	0.498354
8	PowerLOW	2013	S100067	C002353	Boys	7.99452	-0.505476	-0.498822
9	Coordination	2013	S100067	C002353	Boys	7.99452	-0.505476	-0.9773
10	Endurance	2013	S100067	C002353	Boys	7.99452	-0.505476	0.574056
11	Speed	2013	S100067	C002354	Boys	7.99452	-0.505476	0.0551481
12	PowerUP	2013	S100067	C002354	Boys	7.99452	-0.505476	0.218061
13	PowerLOW	2013	S100067	C002354	Boys	7.99452	-0.505476	-0.757248
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
525115	Coordination	2018	S401470	C117964	Boys	9.10609	0.606092	-1.43175
525116	Endurance	2018	S401470	C117964	Boys	9.10609	0.606092	-0.944681
525117	Speed	2018	S401470	C117965	Girls	9.10609	0.606092	0.31086
525118	PowerUP	2018	S401470	C117965	Girls	9.10609	0.606092	0.0779146
525119	PowerLOW	2018	S401470	C117965	Girls	9.10609	0.606092	-0.137027
525120	Coordination	2018	S401470	C117965	Girls	9.10609	0.606092	-1.8077
525121	Endurance	2018	S401470	C117965	Girls	9.10609	0.606092	0.513306
525122	Speed	2018	S800200	C117966	Boys	9.10609	0.606092	0.0551481
525123	PowerUP	2018	S800200	C117966	Boys	9.10609	0.606092	0.0779146
525124	PowerLOW	2018	S800200	C117966	Boys	9.10609	0.606092	-1.32578
525125	Coordination	2018	S800200	C117966	Boys	9.10609	0.606092	0.473217
525126	Endurance	2018	S800200	C117966	Boys	9.10609	0.606092	-0.0941883

Tip

The preprocessing steps above are also packaged as SMLP2026.fggk21_preprocessed(), so that other pages (e.g. the contrast-coding tutorial) can reuse them without repeating the code.

2.2 Extract a stratified subsample

For the purpose of the tutorial, we extract a random sample of 1000 boys and 1000 girls. Child, School, and Cohort are grouping variables. Traditionally, they are called random factors because the units (levels) of the factor are assumed to be a random sample from the population of their units (levels).

Cohort has only nine “groups” and could have been included as a set of polynomical fixed-effect contrasts rather than a random factor. This choice warrants a short excursion: The secular trends are very different for different tests and require the inclusion of interaction terms with Test contrasts (see Figure 4 in (Fühner et al., 2021). The authors opted to absorb these effects in cohort-related variance components for the Test contrasts and plan to address the details of secular changes in a separate analysis.

For complex designs, when they are in the theoretical focus of an article, factors and covariates should be specified as part of the fixed effects. If they are not in the theoretical focus, but serve as statistical control variables, they could be put in the RES - if supported by the data.

Stratified sampling: We build a table with one row per child, sample 1000 boys and 1000 girls from it without replacement, and then extract all of the sampled children’s test scores from df. Sampling whole children – rather than individual scores – preserves the within-child structure that by-Child random slopes rely on. This procedure is packaged as SMLP2026.fggk21_teaching_sample(); it uses a StableRNG with a fixed seed so that the stratification is reproducible across Julia versions (pass a different seed for a different pattern of results). The contrast-coding tutorial uses the same subsample, so estimates are directly comparable between the two pages.

using SMLP2026: fggk21_teaching_sample
dat = fggk21_teaching_sample(df)

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

Due to missing scores for some tests we have about 3% less than 10,000 observations.

2.3 No evidence for `age x Sex x Test` interaction

The main results are captured in the figure constructed in this section. We build it both for the full data and the stratified subset.

df2 = combine(
  groupby(
    select(df, :, :age => ByRow(x -> round(x; digits=1)) => :age),
    [:Sex, :Test, :age],
  ),
  :zScore => mean => :zScore,
  :zScore => length => :n,
)

120×5 DataFrame

95 rows omitted

Row	Sex	Test	age	zScore	n
	Cat…	Cat…	Float64	Float64	Int64
1	Boys	Speed	8.0	-0.0265138	1223
2	Boys	PowerUP	8.0	0.026973	1227
3	Boys	PowerLOW	8.0	0.121609	1227
4	Boys	Coordination	8.0	-0.0571726	1186
5	Boys	Endurance	8.0	0.292695	1210
6	Girls	Speed	8.0	-0.35164	1411
7	Girls	PowerUP	8.0	-0.610355	1417
8	Girls	PowerLOW	8.0	-0.279872	1418
9	Girls	Coordination	8.0	-0.268221	1381
10	Girls	Endurance	8.0	-0.245573	1387
11	Boys	Speed	8.1	0.0608397	3042
12	Boys	PowerUP	8.1	0.0955413	3069
13	Boys	PowerLOW	8.1	0.123099	3069
⋮	⋮	⋮	⋮	⋮	⋮
109	Boys	Coordination	9.0	0.254973	4049
110	Boys	Endurance	9.0	0.258082	4034
111	Girls	Speed	9.1	-0.0286172	1154
112	Girls	PowerUP	9.1	-0.0752301	1186
113	Girls	PowerLOW	9.1	-0.094587	1174
114	Girls	Coordination	9.1	0.00276252	1162
115	Girls	Endurance	9.1	-0.235591	1150
116	Boys	Speed	9.1	0.325745	1303
117	Boys	PowerUP	9.1	0.616416	1320
118	Boys	PowerLOW	9.1	0.267577	1310
119	Boys	Coordination	9.1	0.254342	1297
120	Boys	Endurance	9.1	0.251045	1294

2.3.1 Figure(s) of interaction

The core results of the article are reported in Figure 2 of Fühner et al. (2021). In summary:

Main effects of age and Sex: There are developmental gains in the ninth year of life; boys outperform girls. There is no main effect of Test because of z-scoring.
Interactions of Test and age: Tests differ in how much children improve during the year (i.e., the magnitude of developmental gain), that is slopes depend on Test.
Interactions of Test and Sex: The sex difference is test dependent, that is the difference between the slopes depends on Test.
The most distinctive result is the absence of evidence for an age x Sex x Test interaction, that is the slopes for boys and girls are statistically parallel for each of the five tests.

Code

let
  design1 = mapping(:age, :zScore; color=:Sex, col=:Test)
  lines1 = design1 * linear()
  means1 = design1 * visual(Scatter; markersize=5)
  draw(data(df2) * means1 + data(df) * lines1;)
end

Figure 1: Age trends by sex for each Test for the full data set

Figure 1 shows performance differences for the full set of data between 8.0 and 9.2 years by sex in the five physical fitness tests presented as z-transformed data computed separately for each test.

Endurance = cardiorespiratory endurance (i.e., 6-min-run test),
Coordination = star-run test,
Speed = 20-m linear sprint test,
PowerLOW = power of lower limbs (i.e., standing long jump test),
PowerUP = power of upper limbs (i.e., ball push test),
SD = standard deviation. Points are binned observed child means; lines are simple regression fits to the observations.

What do the results look like for the stratified subsample? Here the parallelism is much less clear. In the final LMM we test whether the two regression lines in each of the five panels are statistically parallel for this subset of data. That is, we test the interaction of Sex and age as nested within the levels of Test. Most people want to know the signficance of these five Sex x age interactions.

The theoretical focus of the article, however, is on comparisons between tests displayed next to each other. We ask whether the degree of parallelism is statistically the same for Endurance and Coordination (H1), Coordination and Speed (H2), Speed and PowerLOW (H3), and PowerLow and PowerUP (H4). Hypotheses H1 to H4 require Sequential Difference contrasts c1 to c4 for Test; they are tested as fixed effects for`H1 x age x Sex, H2 x age x Sex, H3 x age x Sex, and H4 x age x Sex.

Code

dat2 = combine(
  groupby(
    select(dat, :, :age => ByRow(x -> round(x; digits=1)) => :age),
    [:Sex, :Test, :age],
  ),
  :zScore => mean => :zScore,
  :zScore => length => :n,
)

120×5 DataFrame

95 rows omitted

Row	Sex	Test	age	zScore	n
	Cat…	Cat…	Float64	Float64	Int64
1	Boys	Speed	8.0	0.0789166	22
2	Boys	PowerUP	8.0	0.64487	22
3	Boys	PowerLOW	8.0	-0.0195608	22
4	Boys	Coordination	8.0	0.411284	21
5	Boys	Endurance	8.0	0.485693	22
6	Girls	Speed	8.0	-0.439071	31
7	Girls	PowerUP	8.0	-0.688219	30
8	Girls	PowerLOW	8.0	-0.343767	31
9	Girls	Coordination	8.0	-0.410926	31
10	Girls	Endurance	8.0	-0.42921	30
11	Girls	Speed	8.1	-0.346915	64
12	Girls	PowerUP	8.1	-0.541795	64
13	Girls	PowerLOW	8.1	-0.198557	63
⋮	⋮	⋮	⋮	⋮	⋮
109	Girls	Coordination	9.0	0.0919877	75
110	Girls	Endurance	9.0	-0.240405	71
111	Boys	Speed	9.1	0.403568	28
112	Boys	PowerUP	9.1	0.527349	29
113	Boys	PowerLOW	9.1	0.445292	30
114	Boys	Coordination	9.1	0.309464	29
115	Boys	Endurance	9.1	0.304958	30
116	Girls	Speed	9.1	-0.0798607	17
117	Girls	PowerUP	9.1	-0.482671	17
118	Girls	PowerLOW	9.1	-0.430986	16
119	Girls	Coordination	9.1	-0.57964	17
120	Girls	Endurance	9.1	-0.226234	16

Code

let
  design2 = mapping(:age, :zScore; color=:Sex, col=:Test)
  lines2 = design2 * linear()
  means2 = design2 * visual(Scatter; markersize=5)
  draw(data(dat2) * means2 + data(dat) * lines2;)
end

Figure 2: Age trends by sex for each Test for the stratified sample

Figure 2 Performance differences for subset of data between 8.0 and 9.2 years by sex in the five physical fitness tests presented as z-transformed data computed separately for each test.

Endurance = cardiorespiratory endurance (i.e., 6-min-run test),
Coordination = star-run test,
Speed = 20-m linear sprint test,
PowerLOW = power of lower limbs (i.e., standing long jump test),
PowerUP = power of upper limbs (i.e., ball push test),
SD = standard deviation. Points are binned observed child means; lines are simple regression fits to the observations.

2.3.2 Regression on `age` by `Sex` for each `Test`

Another set of relevant statistics are the slopes for the regression of performance on age for boys and girls in each of the five tests. The lines in Figures 1 and 2, however, are computed directly from the raw data with the linear() command.

combine(
  groupby(df, [:Sex, :Test]),
  [:age, :zScore] => simplelinreg => :coef,
)

10×3 DataFrame

Row	Sex	Test	coef
	Cat…	Cat…	Tuple…
1	Boys	Endurance	(0.00256718, 0.0291899)
2	Boys	Coordination	(-2.47279, 0.302819)
3	Boys	Speed	(-2.12689, 0.267153)
4	Boys	PowerLOW	(-1.4307, 0.189659)
5	Boys	PowerUP	(-4.35864, 0.549005)
6	Girls	Endurance	(-0.692022, 0.0523217)
7	Girls	Coordination	(-2.50524, 0.279119)
8	Girls	Speed	(-2.34431, 0.255687)
9	Girls	PowerLOW	(-1.87241, 0.196917)
10	Girls	PowerUP	(-4.82271, 0.524799)

combine(
  groupby(dat, [:Sex, :Test]),
  [:age, :zScore] => simplelinreg => :coef,
)

10×3 DataFrame

Row	Sex	Test	coef
	Cat…	Cat…	Tuple…
1	Boys	Endurance	(-0.809742, 0.128104)
2	Boys	Coordination	(-1.79221, 0.222836)
3	Boys	Speed	(-2.4715, 0.304701)
4	Boys	PowerLOW	(-2.76169, 0.343944)
5	Boys	PowerUP	(-3.88521, 0.492747)
6	Girls	Endurance	(-1.42132, 0.137163)
7	Girls	Coordination	(-2.14706, 0.237589)
8	Girls	Speed	(-2.62018, 0.286215)
9	Girls	PowerLOW	(-0.830218, 0.0721176)
10	Girls	PowerUP	(-4.08892, 0.436415)

2.4 SeqDiffCoding of `Test`

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

H1: Star_r - Run (2-1)
H2: S20_r - Star_r (3-2)
H3: SLJ - S20_r (4-3)
H4: BPT - SLJ (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.

The first command recodes names indicating the physical fitness components used in the above figures and tables back to the shorter actual test names. This reduces clutter in LMM outputs.

recode!(
  dat.Test,
  "Endurance" => "Run",
  "Coordination" => "Star_r",
  "Speed" => "S20_r",
  "PowerLOW" => "SLJ",
  "PowerUP" => "BPT",
)
contrasts = Dict{Symbol,Any}(nm => SeqDiffCoding() for nm in (:Test, :Sex))

Dict{Symbol, Any} with 2 entries:
  :Test => SeqDiffCoding(nothing)
  :Sex  => SeqDiffCoding(nothing)

See the contrast-coding tutorial for why sequential-difference contrasts are not orthogonal and what that implies for the interpretation of correlation parameters between neighboring contrasts, as well as for alternative contrast schemes for these data; contrast coding of visual-attention effects works through the general options on a four-level factor.

3 Model complexification

We fit and compare three LMMs with the same fixed-effect structure but increasing complexity of the random-effect structure for School. We ignore the other two random factors Child and Cohort to avoid undue delays when fitting the models.

LMM m_ovi: allowing only varying intercepts (“Grand Means”);
LMM m_zcp: adding variance components (VCs) for the four Test contrasts, Sex, and age to LMM m_ovi, yielding the zero-correlation parameters LMM;
LMM m_cpx: adding correlation parameters (CPs) to LMM m_zcp; yielding a complex LMM.

In a final part illustrate how to check whether the complex model is supported by the data, rather than leading to a singular fit and, if supported by the data, whether there is an increase in goodness of fit associated with the model complexification.

3.1 LMM `m_ovi`

In its random-effect structure (RES) we only vary intercepts (i.e., Grand Means) for School (LMM m_ovi), that is we allow that the schools differ in the average fitness of its children, average over the five tests.

It is well known that such a simple RES is likely to be anti-conservative with respect to fixed-effect test statistics.

m_ovi = let
  f = @formula zScore ~ 1 + Test * Sex * a1 + (1 | School)
  fit(MixedModel, f, dat; contrasts, progress)
end

Is the model singular (overparameterized, degenerate)? In other words: Is the model not supported by the data?

issingular(m_ovi)

false

Models varying only in intercepts are almost always supported by the data.

3.2 LMM `m_zcp`

In this LMM we allow that schools differ not only in GM, but also in the size of the four contrasts defined for Test, in the difference between boys and girls (Sex) and the developmental gain children achieve within the third grade (age).

We assume that there is covariance associated with these CPs beyond residual noise, that is we assume that there is no detectable evidence in the data that the CPs are different from zero.

m_zcp = let
  f = @formula(
    zScore ~
      1 + Test * Sex * a1 + zerocorr(1 + Test + Sex + a1 | School)
  )
  fit(MixedModel, f, dat; contrasts, progress)
end

Depending on sampling, this model estimating variance components for School may or may not be supported by the data.

issingular(m_zcp)

true

3.3 LMM `m_cpx`

In the complex LMM investigated in this sequence we give up the assumption of zero-correlation between VCs.

m_cpx = let
  f = @formula(
    zScore ~ 1 + Test * Sex * a1 + (1 + Test + Sex + a1 | School)
  )
  fit(MixedModel, f, dat; contrasts, progress)
end

We also need to see the VCs and CPs of the random-effect structure (RES).

VarCorr(m_cpx)

	Column	Variance	Std.Dev	Corr.
School	(Intercept)	0.0876384	0.2960378
	Test: Star_r	0.0971150	0.3116328	+0.11
	Test: S20_r	0.0356897	0.1889171	+0.14	-0.77
	Test: SLJ	0.0059103	0.0768787	-0.59	-0.35	-0.01
	Test: BPT	0.0314018	0.1772055	-0.54	+0.01	+0.30	-0.02
	Sex: Girls	0.2341024	0.4838413	-0.16	+0.06	+0.12	+0.45	-0.04
	a1	0.6565138	0.8102554	-0.10	+0.05	-0.10	+0.54	+0.13	-0.03
Residual		0.7306303	0.8547691

issingular(m_cpx)

false

The complex model may or may not be supported by the data.

3.4 Model comparisons

The singularity checks depend on the subsample: with the present one, m_ovi and m_cpx are supported by the data, while the zero-correlation model m_zcp is singular (one of its variance components is estimated as zero). Does model complexification also increase the goodness of fit or are we only fitting noise?

3.4.1 LRT and goodness-of-fit statistics

As the thee models are strictly hierarchically nested, we compare them with a likelihood-ratio tests (LRT) and AIC and BIC goodness-of-fit statistics derived from them.

MixedModels.likelihoodratiotest(m_ovi, m_zcp, m_cpx)

	model-dof	-2 logLik	χ²	χ²-dof	P(>χ²)
zScore ~ 1 + Test + Sex + a1 + Test & Sex + Test & a1 + Sex & a1 + Test & Sex & a1 + (1 \| School)	22	-26250
zScore ~ 1 + Test + Sex + a1 + Test & Sex + Test & a1 + Sex & a1 + Test & Sex & a1 + zerocorr(1 + Test + Sex + a1 \| School)	28	-25915	336	6	<1e-68
zScore ~ 1 + Test + Sex + a1 + Test & Sex + Test & a1 + Sex & a1 + Test & Sex & a1 + (1 + Test + Sex + a1 \| School)	49	-25867	47	21	0.0009

Code

gof_summary = let
  nms = [:m_ovi, :m_zcp, :m_cpx]
  mods = eval.(nms)
  DataFrame(;
    name=nms,
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
  )
end

3×6 DataFrame

Row	name	dof	deviance	AIC	AICc	BIC
	Symbol	Int64	Float64	Float64	Float64	Float64
1	m_ovi	22	26250.4	26294.4	26294.5	26452.3
2	m_zcp	28	25914.7	25970.7	25970.9	26171.7
3	m_cpx	49	25867.5	25965.5	25966.0	26317.3

These statistics will depend on sampling. In general, smaller deviance, AIC, and BIC indicate an improvement in goodness of fit. Usually, χ² should be larger than the associated degrees of freedom; for AIC and BIC the decrease should amount to more than 5, according to some literature. Severity of meeting these criteria increases from deviance to AIC to BIC. Therefore, it is not always the case that the criteria are unanimous in their verdict. Basically, the more confirmatory the analysis, the more one may go with deviance and AIC; for exploratory analyses the BIC is probably a better guide. There are grey zones here.

3.4.2 Comparing fixed effects of `m_ovi`, `m_zcp`, and `m_cpx`

We check whether enriching the RES changed the significance of fixed effects in the final model.

Code

m_ovi_fe = DataFrame(coeftable(m_ovi));
m_zcp_fe = DataFrame(coeftable(m_zcp));
m_cpx_fe = DataFrame(coeftable(m_cpx));
m_all = hcat(
  m_ovi_fe[:, [1, 2, 4]],
  leftjoin(
    m_zcp_fe[:, [1, 2, 4]],
    m_cpx_fe[:, [1, 2, 4]];
    on=:Name,
    makeunique=true,
  );
  makeunique=true,
)
rename!(
  m_all,
  "Coef." => "b_ovi",
  "Coef._2" => "b_zcp",
  "Coef._1" => "b_cpx",
  "z" => "z_ovi",
  "z_2" => "z_zcp",
  "z_1" => "z_cpx",
)
m_all2 =
  round.(
    m_all[:, [:b_ovi, :b_zcp, :b_cpx, :z_ovi, :z_zcp, :z_cpx]],
    digits=2,
  )
m_all3 = hcat(m_all.Name, m_all2)

20×7 DataFrame

Row	x1	b_ovi	b_zcp	b_cpx	z_ovi	z_zcp	z_cpx
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	(Intercept)	-0.02	-0.03	-0.03	-1.31	-1.45	-1.58
2	Test: Star_r	-0.03	-0.02	-0.01	-0.87	-0.58	-0.45
3	Test: S20_r	-0.02	-0.02	-0.02	-0.6	-0.61	-0.62
4	Test: SLJ	0.01	0.01	0.01	0.22	0.23	0.18
5	Test: BPT	-0.01	-0.0	-0.0	-0.34	-0.01	-0.05
6	Sex: Girls	-0.42	-0.4	-0.39	-20.6	-11.93	-11.78
7	a1	0.23	0.17	0.18	6.52	3.01	3.15
8	Test: Star_r & Sex: Girls	0.31	0.31	0.31	5.07	5.39	5.31
9	Test: S20_r & Sex: Girls	-0.08	-0.08	-0.08	-1.26	-1.35	-1.38
10	Test: SLJ & Sex: Girls	-0.08	-0.07	-0.07	-1.27	-1.31	-1.34
11	Test: BPT & Sex: Girls	-0.3	-0.31	-0.31	-5.12	-5.46	-5.46
12	Test: Star_r & a1	0.1	0.1	0.11	0.96	1.06	1.12
13	Test: S20_r & a1	0.06	0.07	0.05	0.64	0.72	0.53
14	Test: SLJ & a1	-0.09	-0.09	-0.09	-0.89	-0.95	-0.96
15	Test: BPT & a1	0.27	0.26	0.28	2.66	2.73	2.9
16	Sex: Girls & a1	-0.07	-0.08	-0.07	-1.06	-1.05	-0.94
17	Test: Star_r & Sex: Girls & a1	0.01	0.04	0.05	0.04	0.19	0.24
18	Test: S20_r & Sex: Girls & a1	-0.05	-0.05	-0.07	-0.24	-0.28	-0.38
19	Test: SLJ & Sex: Girls & a1	-0.23	-0.23	-0.23	-1.12	-1.2	-1.2
20	Test: BPT & Sex: Girls & a1	0.21	0.19	0.19	1.06	0.97	1.0

The three models usually do not differ in fixed-effect estimates. For main effects of age and Sex, z-values decrease strongly with the complexity of the model (i.e., standard errors are larger). For other coefficients, the changes are not very large and not consistent.

In general, dropping significant variance components and/or correlation parameters may lead to anti-conservative estimates of fixed effects (e.g., Schielzeth & Forstmeier, 2008). Basically, some of the variance allocated to age and Sex in LMM m_ovi could also be due to differences between schools. This ambiguity increased the uncertainty of the respective fixed effects in the other two LMMs.

3.5 Fitting an overparameterized LMM

The complex LMM was not overparameterized with respect to School, because there are over 400 schools in the data. When the number of units (levels) of a grouping factor is small relative to the number of parameters we are trying to estimate, we often end up with an overparameterized / degenerate random-effect structure.

As an illustration, we fit a full CP matrix for the Cohort. As there are only nine cohorts in the data, we may be asking too much to estimate 5*6/2 = 15 VC/CP parameters.

m_cpxCohort = let
  f = @formula zScore ~ 1 + Test * a1 * Sex + (1 + Test | Cohort)
  fit(MixedModel, f, dat; contrasts, progress)
end

VarCorr(m_cpxCohort)

	Column	Variance	Std.Dev	Corr.
Cohort	(Intercept)	0.0002056	0.0143393
	Test: Star_r	0.0015819	0.0397730	-0.54
	Test: S20_r	0.0046825	0.0684287	-0.51	-0.45
	Test: SLJ	0.0025294	0.0502930	+0.97	-0.31	-0.71
	Test: BPT	0.0022986	0.0479434	-0.64	+0.99	-0.33	-0.44
Residual		0.9219180	0.9601656

issingular(m_cpxCohort)

true

The model is overparameterized, with correlation parameters estimated close to ±1. How about the zero-correlation parameter (zcp) version of this LMM?

m_zcpCohort = let
  f = @formula(
    zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Cohort)
  )
  fit(MixedModel, f, dat; contrasts, progress)
end

issingular(m_zcpCohort)

true

This zcpLMM is also singular: with the present subsample, all five VCs are estimated as zero. This raises the possibility that LMM m_oviCohort might fit as well as LMM m_zcpCohort.

m_oviCohort = let
  f = @formula zScore ~ 1 + Test * a1 * Sex + (1 | Cohort)
  fit(MixedModel, f, dat; contrasts, progress)
end

issingular(m_oviCohort)

true

With the present subsample even this varying-intercepts model estimates the Cohort VC as zero — there is essentially no cohort-related variance left in the reduced data. Does LMM m_zcpCohort fit noise relative to the LMM m_oviCohort?

MixedModels.likelihoodratiotest(m_oviCohort, m_zcpCohort)

	model-dof	-2 logLik	χ²	χ²-dof	P(>χ²)
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 \| Cohort)	22	-26761
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + zerocorr(1 + Test \| Cohort)	26	-26761	0	4	1.0000

gof_summary2 = let
  mods = [m_oviCohort, m_zcpCohort, m_cpxCohort]
  DataFrame(;
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
  )
end

3×5 DataFrame

Row	dof	deviance	AIC	AICc	BIC
	Int64	Float64	Float64	Float64	Float64
1	22	26761.1	26805.1	26805.2	26963.0
2	26	26761.1	26813.1	26813.2	26999.8
3	36	26758.7	26830.7	26831.0	27089.2

Indeed, adding VCs is fitting noise. Again, the goodness of fit statistics unanimously favor the selection of the LMM m_oviCohort.

Not shown here, but the Cohort-related VCs for the Test contrasts could be estimated reliably for the full data. Thus, the small number of cohorts does not necessarily prevent the determination of reliable differences between tests across cohorts. What if we include VCs and CPs related to random factors Child and School?

3.6 Fitting the published LMM `m1` to the reduced data

Warning

The following LMMs m1, m2, etc. take a bit longer to fit (a few minutes each on a recent laptop). As elsewhere in these materials, the rendered page restores the fits from a cache via fit_or_restore, so re-rendering is fast; delete the corresponding files under fits/ to force a fresh fit.

LMM m1 reported in Fühner et al. (2021) included random factors for School, Child, and Cohort. The RES for School was specified like in LMM m_cpx. The RES for Child included VCs and CPs for Test, but not for linear developmental gain in the ninth year of life a1 or Sex; they are between-Child effects.

The RES for Cohort included only VCs, no CPs for Test. The parsimony was due to the small number of nine levels for this grouping factor.

Here we fit this LMM m1 for the reduced data. For a different subset of similar size on MacBook Pro [13 | 15 | 16] this took [303 | 250 | 244 ] s; for LMM m1a (i.e., dropping 1 school-relate VC for Sex), times are [212 | 165 | 160] s. The corresponding lme4 times for LMM m1 are [397 | 348 | 195].

Finally, times for fitting the full set of data –not in this script–, for LMM m1are [60 | 62 | 85] minutes (!); for LMM m1a the times were [46 | 48 | 34] minutes. It was not possible to fit the full set of data with lme4; after about 13 to 18 minutes the program stopped with: Error in eval_f(x, ...) : Downdated VtV is not positive definite.

m1 = let
  f = @formula(
    zScore ~
      1 +
      Test * a1 * Sex +
      (1 + Test + a1 + Sex | School) +
      (1 + Test | Child) +
      zerocorr(1 + Test | Cohort)
  )
  fit(MixedModel, f, dat; contrasts, progress)
end

m1

	Est.	SE	z	p	σ_Child	σ_School	σ_Cohort
(Intercept)	-0.0274	0.0186	-1.47	0.1405	0.5780	0.2141	0.0000
Test: Star_r	-0.0104	0.0323	-0.32	0.7474	0.6964	0.3888	0.0000
Test: S20_r	-0.0277	0.0317	-0.87	0.3818	0.5991	0.3225	0.0377
Test: SLJ	0.0038	0.0270	0.14	0.8879	0.4485	0.2712	0.0242
Test: BPT	0.0067	0.0285	0.23	0.8152	0.7018	0.2425	0.0187
a1	0.2359	0.0520	4.54	<1e-05		0.1752
Sex: Girls	-0.4194	0.0308	-13.62	<1e-41		0.1264
Test: Star_r & a1	0.1165	0.0866	1.35	0.1786
Test: S20_r & a1	0.0393	0.0814	0.48	0.6295
Test: SLJ & a1	-0.0911	0.0738	-1.23	0.2172
Test: BPT & a1	0.2680	0.0844	3.18	0.0015
Test: Star_r & Sex: Girls	0.3112	0.0510	6.10	<1e-08
Test: S20_r & Sex: Girls	-0.0674	0.0477	-1.41	0.1578
Test: SLJ & Sex: Girls	-0.0748	0.0433	-1.73	0.0841
Test: BPT & Sex: Girls	-0.3128	0.0496	-6.31	<1e-09
a1 & Sex: Girls	-0.0599	0.1023	-0.59	0.5581
Test: Star_r & a1 & Sex: Girls	0.0486	0.1731	0.28	0.7791
Test: S20_r & a1 & Sex: Girls	-0.1097	0.1617	-0.68	0.4974
Test: SLJ & a1 & Sex: Girls	-0.1900	0.1469	-1.29	0.1959
Test: BPT & a1 & Sex: Girls	0.1881	0.1683	1.12	0.2636
Residual	0.5564

VarCorr(m1)

	Column	Variance	Std.Dev	Corr.
Child	(Intercept)	0.3340612	0.5779802
	Test: Star_r	0.4849610	0.6963914	+0.13
	Test: S20_r	0.3589109	0.5990917	+0.00	-0.46
	Test: SLJ	0.2011919	0.4485442	-0.04	-0.10	-0.42
	Test: BPT	0.4925463	0.7018164	-0.33	+0.11	-0.13	-0.15
School	(Intercept)	0.0458521	0.2141311
	Test: Star_r	0.1511322	0.3887573	+0.02
	Test: S20_r	0.1040291	0.3225354	+0.11	-0.61
	Test: SLJ	0.0735391	0.2711809	-0.15	-0.13	-0.42
	Test: BPT	0.0588117	0.2425113	-0.13	+0.01	+0.12	-0.41
	a1	0.0306948	0.1751994	+0.29	-0.19	+0.03	+0.61	+0.18
	Sex: Girls	0.0159679	0.1263641	-0.58	+0.33	+0.02	+0.30	-0.09	+0.18
Cohort	(Intercept)	0.0000000	0.0000000
	Test: Star_r	0.0000000	0.0000000	.
	Test: S20_r	0.0014220	0.0377094	.	.
	Test: SLJ	0.0005858	0.0242029	.	.	.
	Test: BPT	0.0003508	0.0187294	.	.	.	.
Residual		0.3095329	0.5563568

issingular(m1)

true

Depending on the random number seed for the stratified sampling, LMM m1 may or may not be supported by the data; with the present subsample it is singular.

We also fit an alternative parameterization, estimating VCs and CPs for Test scores rather than Test effects by replacing the 1 + ... in the RE terms with 0 + ....

m2 = let
  f = @formula(
    zScore ~
      1 +
      Test * a1 * Sex +
      (0 + Test + a1 + Sex | School) +
      (0 + Test | Child) +
      zerocorr(0 + Test | Cohort)
  )
  fit(MixedModel, f, dat; contrasts, progress)
end

m2

	Est.	SE	z	p	σ_Child	σ_School	σ_Cohort
(Intercept)	-0.0286	0.0192	-1.49	0.1373
Test: Star_r	-0.0103	0.0338	-0.31	0.7596	0.8483	0.3255	0.0294
Test: S20_r	-0.0312	0.0370	-0.84	0.3987	0.8275	0.3074	0.0617
Test: SLJ	0.0082	0.0342	0.24	0.8092	0.7861	0.2606	0.0258
Test: BPT	0.0062	0.0291	0.21	0.8305	0.7546	0.2542	0.0000
a1	0.2327	0.0521	4.47	<1e-05		0.1717
Sex: Girls	-0.4198	0.0308	-13.62	<1e-41		0.1275
Test: Star_r & a1	0.1156	0.0870	1.33	0.1839
Test: S20_r & a1	0.0250	0.0821	0.30	0.7612
Test: SLJ & a1	-0.0670	0.0747	-0.90	0.3695
Test: BPT & a1	0.2635	0.0845	3.12	0.0018
Test: Star_r & Sex: Girls	0.3110	0.0510	6.10	<1e-08
Test: S20_r & Sex: Girls	-0.0692	0.0477	-1.45	0.1467
Test: SLJ & Sex: Girls	-0.0726	0.0432	-1.68	0.0931
Test: BPT & Sex: Girls	-0.3134	0.0496	-6.32	<1e-09
a1 & Sex: Girls	-0.0579	0.1023	-0.57	0.5714
Test: Star_r & a1 & Sex: Girls	0.0486	0.1731	0.28	0.7788
Test: S20_r & a1 & Sex: Girls	-0.0988	0.1616	-0.61	0.5409
Test: SLJ & a1 & Sex: Girls	-0.2008	0.1467	-1.37	0.1710
Test: BPT & a1 & Sex: Girls	0.1894	0.1682	1.13	0.2602
Test: Run					0.8177	0.3007	0.0000
Residual	0.4232

issingular(m2)

true

Depending on the random number generator seed, the model may or may not be supported in the alternative parameterization of scores. The fixed-effects profile is not affected (see the model comparisons below).

VarCorr(m2)

	Column	Variance	Std.Dev	Corr.
Child	Test: Run	0.6685769	0.8176655
	Test: Star_r	0.7196973	0.8483498	+0.46
	Test: S20_r	0.6848056	0.8275298	+0.47	+0.56
	Test: SLJ	0.6179891	0.7861228	+0.51	+0.56	+0.65
	Test: BPT	0.5694333	0.7546081	+0.17	+0.32	+0.32	+0.37
School	Test: Run	0.0904482	0.3007460
	Test: Star_r	0.1059671	0.3255260	+0.23
	Test: S20_r	0.0944839	0.3073823	+0.45	+0.48
	Test: SLJ	0.0679284	0.2606308	+0.54	+0.41	+0.55
	Test: BPT	0.0646129	0.2541907	+0.37	+0.26	+0.51	+0.56
	a1	0.0294933	0.1717361	+0.14	-0.10	-0.06	+0.55	+0.74
	Sex: Girls	0.0162487	0.1274703	-0.85	-0.40	-0.42	-0.16	-0.25	+0.19
Cohort	Test: Run	0.0000000	0.0000000
	Test: Star_r	0.0008625	0.0293682	.
	Test: S20_r	0.0038116	0.0617382	.	.
	Test: SLJ	0.0006674	0.0258334	.	.	.
	Test: BPT	0.0000000	0.0000000	.	.	.	.
Residual		0.1790652	0.4231610

3.7 Principle Component Analysis of Random Effect Structure (rePCA)

The ìssingular() command is sort of a shortcut for a quick inspection of the principle components (PCs) of the variance-covariance matrix of the RES. With the MixedModels.PCA() command, we also obtain information about the amount of cumulative variance accounted for as we add PCs.

The output also provides PC loadings which may facilitate interpretation of the CP matrices (if estimated). See Fühner et al. (2021) for an application.

3.8 Effects in RES

For every random factor, MixedModels.PCA() extracts as many PCs as there are VCs. Therefore, the cumulation of variance across PCs within a random factor will always add up to 100% – at the latest with the last VC, but, in the case of overparameterized LMMs, the ceiling will be reached earlier. The final PCs are usually quite small.

PCs are extracted in the order of the amount of unique variance they account for. The first PC accounts for the largest and the final PC for the least amount of variance. The number the PCs with percent variance above a certain threshold indicates the number of weighted composites needed and reflects the dimensionality of the orthogonal space within which (almost) all the variance can be accounted for. The weights for forming composite scores are the listed loadings. For ease of interpretation it is often useful to change the sign of some composite scores.

The PCA for LMM m1 shows that each of the five PCs for Child accounts for a non-zero percent of unique variance.

For School fewer than seven PCs have unique variance. The exact number depends on sampling. The overparameterization of School might be resolved when the CPs for Sex are dropped from the LMM.

Cohort was estimated with CPs forced to zero. Therefore, the VCs were forced to be orthogonal; they already represent the PCA solution. However, depending on sampling, not all PCs may be identified for this random factor either.

Importantly, again depending on sampling, a non-singular fit does not imply that unique variance is associated with all PCs (i.e., not for last PC for School). Embrace uncertainty!

MixedModels.PCA(m1)

(Child = 
Principal components based on correlation matrix
 (Intercept)    1.0     .      .      .      .
 Test: Star_r   0.13   1.0     .      .      .
 Test: S20_r    0.0   -0.46   1.0     .      .
 Test: SLJ     -0.04  -0.1   -0.42   1.0     .
 Test: BPT     -0.33   0.11  -0.13  -0.15   1.0

Normalized cumulative variances:
[0.3183, 0.5873, 0.8223, 0.9387, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4   PC5
 (Intercept)    0.0   -0.65  -0.41  -0.64  0.01
 Test: Star_r   0.53   0.01  -0.58   0.37  0.49
 Test: S20_r   -0.73   0.02  -0.07   0.04  0.68
 Test: SLJ      0.4   -0.26   0.69  -0.17  0.52
 Test: BPT      0.15   0.72  -0.11  -0.65  0.17, School = 
Principal components based on correlation matrix
 (Intercept)    1.0     .      .      .      .      .      .
 Test: Star_r   0.02   1.0     .      .      .      .      .
 Test: S20_r    0.11  -0.61   1.0     .      .      .      .
 Test: SLJ     -0.15  -0.13  -0.42   1.0     .      .      .
 Test: BPT     -0.13   0.01   0.12  -0.41   1.0     .      .
 a1             0.29  -0.19   0.03   0.61   0.18   1.0     .
 Sex: Girls    -0.58   0.33   0.02   0.3   -0.09   0.18   1.0

Normalized cumulative variances:
[0.2898, 0.5304, 0.7346, 0.8948, 1.0, 1.0, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5    PC6    PC7
 (Intercept)   -0.3   -0.3    0.59  -0.14   0.44   0.44  -0.25
 Test: Star_r   0.3    0.51   0.33  -0.31   0.4   -0.5   -0.19
 Test: S20_r   -0.45  -0.23  -0.48   0.15   0.46  -0.42  -0.34
 Test: SLJ      0.54  -0.45   0.06   0.11  -0.24  -0.1   -0.65
 Test: BPT     -0.23   0.11  -0.28  -0.79  -0.3    0.17  -0.35
 a1             0.22  -0.61  -0.03  -0.48   0.16  -0.27   0.5
 Sex: Girls     0.47   0.12  -0.48  -0.04   0.51   0.52   0.0, Cohort = 
Principal components based on correlation matrix
 (Intercept)   0.0  .    .    .    .
 Test: Star_r  0.0  0.0  .    .    .
 Test: S20_r   0.0  0.0  1.0  .    .
 Test: SLJ     0.0  0.0  0.0  1.0  .
 Test: BPT     0.0  0.0  0.0  0.0  1.0

Normalized cumulative variances:
[0.3333, 0.6667, 1.0, 1.0, 1.0]

Component loadings
                PC1   PC2   PC3   PC4   PC5
 (Intercept)   0.0   0.0   0.0   1.0   0.0
 Test: Star_r  0.0   0.0   0.0   0.0   1.0
 Test: S20_r   1.0   0.0   0.0   0.0   0.0
 Test: SLJ     0.0   1.0   0.0   0.0   0.0
 Test: BPT     0.0   0.0   1.0   0.0   0.0)

3.8.1 Scores in RES

Now lets looks at the PCA results for the alternative parameterization of LMM m2. It is important to note that the reparameterization to base estimates of VCs and CPs on scores rather than effects applies only to the Test factor (i.e., the first factor in the formula); VCs for Sex and age refer to the associated effects.

Depending on sampling, the difference between LMM m1 and LMM m2 may show that overparameterization according to PCs may depend on the specification chosen for the other the random-effect structure.

Note

For the complete data, all PCs had unique variance associated with them.

MixedModels.PCA(m2)

(Child = 
Principal components based on correlation matrix
 Test: Run     1.0    .     .     .     .
 Test: Star_r  0.46  1.0    .     .     .
 Test: S20_r   0.47  0.56  1.0    .     .
 Test: SLJ     0.51  0.56  0.65  1.0    .
 Test: BPT     0.17  0.32  0.32  0.37  1.0

Normalized cumulative variances:
[0.561, 0.7307, 0.8375, 0.9313, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5
 Test: Run     -0.42   0.45   0.78  -0.03  -0.13
 Test: Star_r  -0.47   0.07  -0.32  -0.82   0.05
 Test: S20_r   -0.49   0.08  -0.41   0.41  -0.64
 Test: SLJ     -0.51   0.03  -0.15   0.4    0.75
 Test: BPT     -0.31  -0.89   0.33  -0.03  -0.1, School = 
Principal components based on correlation matrix
 Test: Run      1.0     .      .      .      .      .      .
 Test: Star_r   0.23   1.0     .      .      .      .      .
 Test: S20_r    0.45   0.48   1.0     .      .      .      .
 Test: SLJ      0.54   0.41   0.55   1.0     .      .      .
 Test: BPT      0.37   0.26   0.51   0.56   1.0     .      .
 a1             0.14  -0.1   -0.06   0.55   0.74   1.0     .
 Sex: Girls    -0.85  -0.4   -0.42  -0.16  -0.25   0.19   1.0

Normalized cumulative variances:
[0.467, 0.7177, 0.8557, 0.93, 1.0, 1.0, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5    PC6    PC7
 Test: Run     -0.43   0.22   0.52   0.05   0.27   0.24   0.6
 Test: Star_r  -0.31   0.25  -0.6    0.62  -0.17  -0.02   0.24
 Test: S20_r   -0.41   0.17  -0.37  -0.71  -0.0   -0.37   0.16
 Test: SLJ     -0.44  -0.25  -0.17   0.08   0.68   0.17  -0.47
 Test: BPT     -0.42  -0.37   0.01  -0.16  -0.59   0.55  -0.11
 a1            -0.22  -0.67   0.19   0.24  -0.09  -0.61   0.18
 Sex: Girls     0.35  -0.47  -0.41  -0.15   0.27   0.32   0.54, Cohort = 
Principal components based on correlation matrix
 Test: Run     0.0  .    .    .    .
 Test: Star_r  0.0  1.0  .    .    .
 Test: S20_r   0.0  0.0  1.0  .    .
 Test: SLJ     0.0  0.0  0.0  1.0  .
 Test: BPT     0.0  0.0  0.0  0.0  0.0

Normalized cumulative variances:
[0.3333, 0.6667, 1.0, 1.0, 1.0]

Component loadings
                PC1   PC2   PC3   PC4   PC5
 Test: Run     0.0   0.0   0.0   1.0   0.0
 Test: Star_r  1.0   0.0   0.0   0.0   0.0
 Test: S20_r   0.0   1.0   0.0   0.0   0.0
 Test: SLJ     0.0   0.0   1.0   0.0   0.0
 Test: BPT     0.0   0.0   0.0   0.0   1.0)

3.9 Summary of results for stratified subset of data

Returning to the theoretical focus of the article, the significant main effects of age and Sex are replicated, as are the interaction between age and the c4 contrast and the interactions between Sex and two test contrasts (c1, c4). The remaining contrast interactions reported for the full data do not reach significance in this subsample. Obviously, the subset of data is much noisier than the full set.

4 `Age x Sex` nested in levels of `Test`

In this final LMM, we test post-hoc five age x Sex interactions by nesting the interaction in the levels of Test. As this LMM m2_nested is a reparameterization of LMM m2.

m2_nested = let
  f = @formula(
    zScore ~
      1 +
      Test +
      Test & (a1 * Sex) +
      (0 + Test + a1 + Sex | School) +
      (0 + Test | Child) +
      zerocorr(0 + Test | Cohort)
  )
  fit(MixedModel, f, dat; contrasts, progress)
end

m2_nested

	Est.	SE	z	p	σ_Child	σ_School	σ_Cohort
(Intercept)	-0.0286	0.0192	-1.49	0.1372
Test: Star_r	-0.0103	0.0338	-0.31	0.7595	0.8510	0.3256	0.0294
Test: S20_r	-0.0312	0.0370	-0.84	0.3988	0.8303	0.3074	0.0618
Test: SLJ	0.0082	0.0342	0.24	0.8093	0.7890	0.2606	0.0259
Test: BPT	0.0062	0.0291	0.21	0.8305	0.7576	0.2542	0.0000
Test: Run & a1	0.0993	0.0757	1.31	0.1893
Test: Star_r & a1	0.2149	0.0784	2.74	0.0061
Test: S20_r & a1	0.2399	0.0771	3.11	0.0019
Test: SLJ & a1	0.1729	0.0728	2.37	0.0176
Test: BPT & a1	0.4364	0.0704	6.20	<1e-09
Test: Run & Sex: Girls	-0.5354	0.0448	-11.95	<1e-32
Test: Star_r & Sex: Girls	-0.2244	0.0462	-4.86	<1e-05
Test: S20_r & Sex: Girls	-0.2936	0.0450	-6.52	<1e-10
Test: SLJ & Sex: Girls	-0.3662	0.0428	-8.56	<1e-16
Test: BPT & Sex: Girls	-0.6795	0.0415	-16.37	<1e-59
Test: Run & a1 & Sex: Girls	0.0049	0.1503	0.03	0.9740
Test: Star_r & a1 & Sex: Girls	0.0535	0.1548	0.35	0.7295
Test: S20_r & a1 & Sex: Girls	-0.0453	0.1511	-0.30	0.7645
Test: SLJ & a1 & Sex: Girls	-0.2460	0.1437	-1.71	0.0869
Test: BPT & a1 & Sex: Girls	-0.0566	0.1396	-0.41	0.6850
Test: Run					0.8204	0.3008	0.0000
a1						0.1719
Sex: Girls						0.1274
Residual	0.4177

The results show that none of the interactions in the panels of Figure 2 is significant. The size and direction of interaction effects correspond with what is shown in Figure 2.

4.0.1 Model comparisons

Code

gof_summary3 = let
  nms = [:m1, :m2, :m2_nested]
  mods = eval.(nms)
  DataFrame(;
    name=nms,
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
  )
end

3×6 DataFrame

Row	name	dof	deviance	AIC	AICc	BIC
	Symbol	Int64	Float64	Float64	Float64	Float64
1	m1	69	24641.9	24779.9	24780.9	25275.4
2	m2	69	24637.5	24775.5	24776.5	25270.9
3	m2_nested	69	24637.5	24775.5	24776.5	25270.9

In principle, the models should yield the same deviance. When models are not supported by the data, that is for singular models, there may be small differences between deviances for these reparameterizations. During optimization such models search for the absolute minimum in a very shallow surface and may end up in a local minimum instead.

4.1 See Also

Mixed Models Tutorial: Contrast Coding — applies a full range of contrast schemes (effect, Helmert, hypothesis, PCA-based) to the same Emotikon dataset.
Emotikon: Comparing Original and Transformed Metrics — investigates whether the age × Sex × Test interaction changes when tasks are analyzed in their original time metric rather than speed.
Creating Multi-Panel Plots — reproduces the publication figure from Fühner et al. (2021) using AlgebraOfGraphics.
MixedModels.jl cheatsheet — a tour of the extractor functions for inspecting a fitted model, formerly a section of this page.

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

Schielzeth, H., & Forstmeier, W. (2008). Conclusions beyond support: Overconfident estimates in mixed models. Behavioral Ecology, 20(2), 416–420. https://doi.org/10.1093/beheco/arn145

4.2 Exercises

Read the PCA. Inspect the principal-components decomposition of the random effects. How many components carry most of the variance, and how does that inform which random-effects terms to keep?

Solution

If the first few components account for nearly all of the random-effects variance, the remaining dimensions are effectively redundant — the data do not support a full-rank random-effects covariance. That is the quantitative basis for reducing the random-effects structure (e.g. dropping near-zero components or correlation parameters) rather than guessing.

Add a comparison. Add a contrast for the Test factor that encodes a comparison you find substantively interesting, refit, and check whether it changes the conclusions about the age-by-sex interaction.

Solution

Because re-coding a factor does not change the model fit, the omnibus conclusions (and the age-by-sex interaction) are unaffected; what changes is which specific comparison each Test coefficient reports. Use this to confirm that your substantive comparison is or is not significant without disturbing the rest of the model. See Mixed Models Tutorial: Contrast Coding.

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

1 Packages and functions

2 Preprocessing

2.1 Read data

2.1.1 Transformations

2.2 Extract a stratified subsample

2.3 No evidence for age x Sex x Test interaction

2.3.1 Figure(s) of interaction

2.3.2 Regression on age by Sex for each Test

2.4 SeqDiffCoding of Test

3 Model complexification

3.1 LMM m_ovi

3.2 LMM m_zcp

3.3 LMM m_cpx

3.4 Model comparisons

3.4.1 LRT and goodness-of-fit statistics

3.4.2 Comparing fixed effects of m_ovi, m_zcp, and m_cpx

3.5 Fitting an overparameterized LMM

3.6 Fitting the published LMM m1 to the reduced data

3.7 Principle Component Analysis of Random Effect Structure (rePCA)

3.8 Effects in RES

3.8.1 Scores in RES

3.9 Summary of results for stratified subset of data

4 Age x Sex nested in levels of Test

4.0.1 Model comparisons

4.1 See Also

4.2 Exercises

2.3 No evidence for `age x Sex x Test` interaction

2.3.2 Regression on `age` by `Sex` for each `Test`

2.4 SeqDiffCoding of `Test`

3.1 LMM `m_ovi`

3.2 LMM `m_zcp`

3.3 LMM `m_cpx`

3.4.2 Comparing fixed effects of `m_ovi`, `m_zcp`, and `m_cpx`

3.6 Fitting the published LMM `m1` to the reduced data

4 `Age x Sex` nested in levels of `Test`