The Emotikon Project: Model Complexification and PCA

Published

2026-04-10

This script uses a subset of data reported in @Fuehner2021. 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 four main sections:

  1. Setup with reading and examining the data, plotting the main results, and specifying the contrasts for the fixed factor Test
  2. 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,
  3. specification of nested fixed effects or interactions in the levels of another, superordinate factors,
  4. a Glossary of MixedModels.jl commands to inspect the information generated for a fitted model object.

1 Packages and functions

Code 
using AlgebraOfGraphics
using AlgebraOfGraphics: linear
using Arrow
using CairoMakie
using CategoricalArrays
using Chain
using DataFrameMacros
using DataFrames
using MixedModels
using MixedModelsMakie
using MixedModelsMakie: simplelinreg
using Random
using Statistics
using StatsBase
using SMLP2026: dataset, fit_or_restore, DATADIR

const progress = isinteractive()
false

2 Readme

Number of scores: 525126 in dataset(:fggk21)

  1. Cohort: 9 levels; 2011-2019

  2. School: 515 levels

  3. Child: 108295 levels; all children are between 8.0 and 8.99 years old

  4. Sex: “Girls” (n=55,086), “Boys” (n= 53,209)

  5. age: testdate - middle of month of birthdate

  6. 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]

  7. score - see units

3 Preprocessing

3.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

3.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

3.2 Extract a stratified subsample

For the prupose 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 [@Fuehner2021]. 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 generate a Child table with information about children. MersenneTwister(42) specifies 42 as the seed for the random number generator to ensure reproducibility of the stratification. For a different pattern of results choose, for example, 84. We randomly sample 1000 boys and 1000 girls from this table; they are stored in samp. Then, we extract the corresponding subset of these children’s test scores from df and store them dat.

Child = unique(select(df, :Cohort, :School, :Child, :Sex, :age))
sample = let
  rng = MersenneTwister(42)
  combine(
    groupby(Child, :Sex), x -> x[rand(rng, 1:nrow(x), 1000), :]
  )
end
insamp(x) = x ∈ sample.Child
dat = @subset(df, insamp(:Child))
9621×8 DataFrame
9596 rows omitted
Row Test Cohort School Child Sex age a1 zScore
Cat… String String String Cat… Float64 Float64 Float64
1 Speed 2013 S100547 C002378 Boys 7.99452 -0.505476 0.578748
2 PowerUP 2013 S100547 C002378 Boys 7.99452 -0.505476 -0.762963
3 PowerLOW 2013 S100547 C002378 Boys 7.99452 -0.505476 -0.808933
4 Coordination 2013 S100547 C002378 Boys 7.99452 -0.505476 -0.846558
5 Endurance 2013 S100547 C002378 Boys 7.99452 -0.505476 -0.397936
6 Speed 2013 S102489 C002562 Girls 7.99452 -0.505476 -0.189199
7 PowerUP 2013 S102489 C002562 Girls 7.99452 -0.505476 -0.622817
8 PowerLOW 2013 S102489 C002562 Girls 7.99452 -0.505476 -1.63589
9 Coordination 2013 S102489 C002562 Girls 7.99452 -0.505476 -1.62614
10 Endurance 2013 S102489 C002562 Girls 7.99452 -0.505476 0.12181
11 Speed 2013 S103226 C002618 Boys 7.99452 -0.505476 -0.422921
12 PowerUP 2013 S103226 C002618 Boys 7.99452 -0.505476 -0.482671
13 PowerLOW 2013 S103226 C002618 Boys 7.99452 -0.505476 -0.292082
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
9610 PowerLOW 2018 S111570 C117818 Girls 9.10609 0.606092 0.483194
9611 Coordination 2018 S111570 C117818 Girls 9.10609 0.606092 -0.31554
9612 Endurance 2018 S111570 C117818 Girls 9.10609 0.606092 -0.0941883
9613 PowerUP 2018 S111739 C117841 Girls 9.10609 0.606092 0.0779146
9614 PowerLOW 2018 S111739 C117841 Girls 9.10609 0.606092 -1.11904
9615 Coordination 2018 S111739 C117841 Girls 9.10609 0.606092 -0.95589
9616 Endurance 2018 S111739 C117841 Girls 9.10609 0.606092 -2.58492
9617 Speed 2018 S112999 C117902 Girls 9.10609 0.606092 1.7913
9618 PowerUP 2018 S112999 C117902 Girls 9.10609 0.606092 -0.0622317
9619 PowerLOW 2018 S112999 C117902 Girls 9.10609 0.606092 0.483194
9620 Coordination 2018 S112999 C117902 Girls 9.10609 0.606092 0.708474
9621 Endurance 2018 S112999 C117902 Girls 9.10609 0.606092 0.574056

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

3.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

3.3.1 Figure(s) of interaction

The core results of the article are reported in Figure 2 of @Fuehner2021. 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.0108467 20
2 Boys PowerUP 8.0 0.0779146 20
3 Boys PowerLOW 8.0 0.0800505 20
4 Boys Coordination 8.0 0.0843624 19
5 Boys Endurance 8.0 0.287361 19
6 Girls Speed 8.0 -0.430054 27
7 Girls PowerUP 8.0 -0.799298 27
8 Girls PowerLOW 8.0 -0.408374 28
9 Girls Coordination 8.0 -0.323287 28
10 Girls Endurance 8.0 -0.380386 25
11 Girls Speed 8.1 -0.461211 55
12 Girls PowerUP 8.1 -0.566758 55
13 Girls PowerLOW 8.1 -0.309937 55
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
109 Girls Coordination 9.0 -0.0750398 57
110 Girls Endurance 9.0 -0.42162 57
111 Boys Speed 9.1 0.119451 29
112 Boys PowerUP 9.1 0.618479 28
113 Boys PowerLOW 9.1 -0.131489 28
114 Boys Coordination 9.1 0.145517 28
115 Boys Endurance 9.1 -0.0891883 27
116 Girls Speed 9.1 -0.0435788 20
117 Girls PowerUP 9.1 0.230801 22
118 Girls PowerLOW 9.1 -0.0538813 23
119 Girls Coordination 9.1 0.282882 23
120 Girls Endurance 9.1 -0.0865179 22
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.

3.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.681676, 0.111322)
2 Boys Coordination (-1.11067, 0.141854)
3 Boys Speed (-1.37849, 0.178035)
4 Boys PowerLOW (-1.44439, 0.194219)
5 Boys PowerUP (-3.26371, 0.420349)
6 Girls Endurance (-0.350538, 0.0141756)
7 Girls Coordination (-2.08942, 0.234567)
8 Girls Speed (-1.93461, 0.213441)
9 Girls PowerLOW (0.0629858, -0.0240463)
10 Girls PowerUP (-4.30863, 0.469259)

3.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.

Note that random factors Child, School, and Cohort are declared as Grouping variables. Technically, this specification is required for variables with a very large number of levels (e.g., 100K+ children). We recommend the explicit specification for all random factors as a general coding style.

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" => "BMT",
)
contrasts = Dict{Symbol,Any}(nm => SeqDiffCoding() for nm in (:Test, :Sex))
Dict{Symbol, Any} with 2 entries:
  :Test => SeqDiffCoding(nothing)
  :Sex  => SeqDiffCoding(nothing)

The statistical disadvantage of SeqDiffCoding 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 limitations 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.

Various options for contrast coding are the topic of the MixedModelsTutorial_contrasts_emotikon.jl and MixedModelsTutorial_contrasts_kwdyz.jl notebooks.

4 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.

  1. LMM m_ovi: allowing only varying intercepts (“Grand Means”);
  2. 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;
  3. 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.

4.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.

4.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

4.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.1013606 0.3183718
Test: Star_r 0.0648041 0.2545665 +0.15
Test: S20_r 0.0257467 0.1604578 +0.03 -0.37
Test: SLJ 0.0020276 0.0450292 +0.04 -0.65 -0.35
Test: BMT 0.0446242 0.2112445 -0.57 -0.07 -0.34 +0.43
Sex: Girls 0.2424297 0.4923715 -0.42 -0.04 -0.06 +0.20 +0.36
a1 0.6364939 0.7978057 +0.18 +0.01 +0.09 -0.02 -0.14 +0.01
Residual 0.7444343 0.8628060 
issingular(m_cpx)
false

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

4.4 Model comparisons

The checks of model singularity indicate that the three models are supported by the data. Does model complexification also increase the goodness of fit or are we only fitting noise?

4.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 -26190
zScore ~ 1 + Test + Sex + a1 + Test & Sex + Test & a1 + Sex & a1 + Test & Sex & a1 + zerocorr(1 + Test + Sex + a1 | School) 28 -25862 328 6 <1e-67
zScore ~ 1 + Test + Sex + a1 + Test & Sex + Test & a1 + Sex & a1 + Test & Sex & a1 + (1 + Test + Sex + a1 | School) 49 -25814 48 21 0.0006
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 26189.8 26233.8 26234.0 26391.6
2 m_zcp 28 25862.2 25918.2 25918.4 26119.0
3 m_cpx 49 25813.8 25911.8 25912.4 26263.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.

4.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.01 0.0 -0.0 -0.4 0.01 -0.1
2 Test: Star_r -0.02 -0.01 -0.01 -0.68 -0.44 -0.46
3 Test: S20_r 0.01 0.01 0.01 0.32 0.39 0.33
4 Test: SLJ 0.02 0.02 0.02 0.81 0.84 0.78
5 Test: BMT -0.04 -0.03 -0.02 -1.21 -0.92 -0.75
6 Sex: Girls -0.38 -0.38 -0.37 -18.49 -11.17 -10.9
7 a1 0.19 0.19 0.17 5.32 3.13 2.94
8 Test: Star_r & Sex: Girls 0.3 0.31 0.31 5.02 5.37 5.35
9 Test: S20_r & Sex: Girls -0.07 -0.07 -0.07 -1.12 -1.2 -1.17
10 Test: SLJ & Sex: Girls -0.09 -0.09 -0.09 -1.49 -1.53 -1.53
11 Test: BMT & Sex: Girls -0.28 -0.29 -0.29 -4.67 -4.98 -5.03
12 Test: Star_r & a1 0.12 0.11 0.12 1.19 1.15 1.21
13 Test: S20_r & a1 0.01 0.02 0.02 0.11 0.18 0.18
14 Test: SLJ & a1 -0.11 -0.11 -0.11 -1.1 -1.15 -1.18
15 Test: BMT & a1 0.35 0.35 0.35 3.45 3.53 3.56
16 Sex: Girls & a1 0.06 0.1 0.09 0.87 1.21 1.12
17 Test: Star_r & Sex: Girls & a1 0.18 0.11 0.12 0.85 0.56 0.59
18 Test: S20_r & Sex: Girls & a1 -0.05 -0.05 -0.05 -0.22 -0.26 -0.27
19 Test: SLJ & Sex: Girls & a1 -0.26 -0.26 -0.25 -1.26 -1.34 -1.28
20 Test: BMT & Sex: Girls & a1 0.26 0.26 0.24 1.29 1.32 1.24

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., @Schielzeth2008]. 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.

4.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.0011062 0.0332591
Test: Star_r 0.0039335 0.0627180 +0.52
Test: S20_r 0.0041780 0.0646373 -1.00 -0.50
Test: SLJ 0.0021814 0.0467057 +0.58 -0.38 -0.60
Test: BMT 0.0011850 0.0344239 +0.68 +0.94 -0.67 -0.18
Residual 0.9349908 0.9669492 
issingular(m_cpxCohort)
true

The model is overparameterized with several CPs estimated between |.98| and |1.00|. 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. Three of the 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)
false

This solves the problem with singularity, but 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 -26674
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + zerocorr(1 + Test | Cohort) 26 -26674 0 4 0.9997 
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 26674.0 26718.0 26718.1 26875.8
2 26 26674.0 26726.0 26726.1 26912.4
3 36 26669.6 26741.6 26741.9 26999.8

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?

4.6 Fitting the published LMM m1 to the reduced data

Warning

The following LMMs m1, m2, etc. take a bit longer (e.g., close to 6 minutes in the Pluto notebook, close to 3 minutes in the REPL on a MacBook Pro).

LMM m1 reported in @Fuehner2021 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
Est. SE z p σ_Child σ_School σ_Cohort
(Intercept) 0.0005 0.0198 0.02 0.9812 0.6003 0.2027 0.0201
Test: Star_r -0.0213 0.0391 -0.55 0.5848 0.6595 0.3872 0.0685
Test: S20_r 0.0191 0.0323 0.59 0.5546 0.6453 0.3234 0.0409
Test: SLJ 0.0161 0.0260 0.62 0.5340 0.5441 0.2500 0.0086
Test: BMT -0.0187 0.0309 -0.60 0.5455 0.7110 0.2901 0.0314
a1 0.1805 0.0543 3.32 0.0009 0.1663
Sex: Girls -0.3743 0.0329 -11.38 <1e-29 0.2132
Test: Star_r & a1 0.1201 0.0858 1.40 0.1614
Test: S20_r & a1 0.0061 0.0832 0.07 0.9416
Test: SLJ & a1 -0.1207 0.0767 -1.58 0.1152
Test: BMT & a1 0.3397 0.0855 3.97 <1e-04
Test: Star_r & Sex: Girls 0.3113 0.0494 6.30 <1e-09
Test: S20_r & Sex: Girls -0.0634 0.0482 -1.31 0.1890
Test: SLJ & Sex: Girls -0.0892 0.0447 -1.99 0.0463
Test: BMT & Sex: Girls -0.2890 0.0498 -5.81 <1e-08
a1 & Sex: Girls 0.0233 0.1066 0.22 0.8270
Test: Star_r & a1 & Sex: Girls 0.0481 0.1700 0.28 0.7771
Test: S20_r & a1 & Sex: Girls -0.0716 0.1656 -0.43 0.6656
Test: SLJ & a1 & Sex: Girls -0.2182 0.1535 -1.42 0.1551
Test: BMT & a1 & Sex: Girls 0.2560 0.1705 1.50 0.1333
Residual 0.5394 
VarCorr(m1)
Column Variance Std.Dev Corr.
Child (Intercept) 0.36033239 0.60027693
Test: Star_r 0.43495553 0.65951159 +0.04
Test: S20_r 0.41638097 0.64527589 +0.11 -0.49
Test: SLJ 0.29608881 0.54414043 +0.05 -0.02 -0.41
Test: BMT 0.50554575 0.71101741 -0.31 +0.17 -0.17 -0.30
School (Intercept) 0.04109627 0.20272215
Test: Star_r 0.14992416 0.38720041 +0.16
Test: S20_r 0.10460004 0.32341929 -0.20 -0.48
Test: SLJ 0.06249852 0.24999704 -0.03 -0.14 -0.46
Test: BMT 0.08414285 0.29007387 -0.14 -0.09 -0.11 -0.28
a1 0.02766141 0.16631719 +0.81 -0.08 +0.39 -0.19 -0.19
Sex: Girls 0.04547115 0.21323966 -0.97 -0.18 +0.09 +0.23 +0.18 -0.81
Cohort (Intercept) 0.00040530 0.02013200
Test: Star_r 0.00469459 0.06851710 .
Test: S20_r 0.00167603 0.04093934 . .
Test: SLJ 0.00007426 0.00861717 . . .
Test: BMT 0.00098435 0.03137434 . . . .
Residual 0.29100268 0.53944665 
issingular(m1)
true

Depending on the random number for stratified samplign, LMM m1 may or may not be supported by the data.

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
Est. SE z p σ_Child σ_School σ_Cohort
(Intercept) 0.0000 0.0199 0.00 0.9984
Test: Star_r -0.0223 0.0424 -0.53 0.5991 0.7374 0.3453 0.0681
Test: S20_r 0.0188 0.0401 0.47 0.6392 0.7822 0.2832 0.0449
Test: SLJ 0.0212 0.0315 0.67 0.5003 0.7924 0.2497 0.0289
Test: BMT -0.0205 0.0321 -0.64 0.5233 0.7076 0.2450 0.0289
a1 0.1802 0.0542 3.32 0.0009 0.1653
Sex: Girls -0.3743 0.0329 -11.38 <1e-29 0.2132
Test: Star_r & a1 0.1185 0.0861 1.38 0.1686
Test: S20_r & a1 0.0112 0.0839 0.13 0.8942
Test: SLJ & a1 -0.1095 0.0774 -1.41 0.1571
Test: BMT & a1 0.3384 0.0857 3.95 <1e-04
Test: Star_r & Sex: Girls 0.3117 0.0494 6.31 <1e-09
Test: S20_r & Sex: Girls -0.0653 0.0482 -1.35 0.1757
Test: SLJ & Sex: Girls -0.0864 0.0447 -1.93 0.0535
Test: BMT & Sex: Girls -0.2900 0.0498 -5.83 <1e-08
a1 & Sex: Girls 0.0236 0.1066 0.22 0.8245
Test: Star_r & a1 & Sex: Girls 0.0510 0.1699 0.30 0.7641
Test: S20_r & a1 & Sex: Girls -0.0757 0.1655 -0.46 0.6472
Test: SLJ & a1 & Sex: Girls -0.2174 0.1534 -1.42 0.1563
Test: BMT & a1 & Sex: Girls 0.2571 0.1705 1.51 0.1315
Test: Run 0.7398 0.3070 0.0494
Residual 0.5234 
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 2.8 below).

Caution

RK: The order of RE terms is critical. In formula f2 the zerocorr() term must be placed last as shown. If it is placed first, School-related and Child-related CPs are estimated/reported (?) as zero. This was not the case for formula m1. Thus, it appears to be related to the 0-intercepts in School and Child terms. Need a reprex.

VarCorr(m2)
Column Variance Std.Dev Corr.
Child Test: Run 0.5473394 0.7398239
Test: Star_r 0.5437814 0.7374153 +0.57
Test: S20_r 0.6118591 0.7822142 +0.60 +0.61
Test: SLJ 0.6278340 0.7923598 +0.61 +0.62 +0.74
Test: BMT 0.5006439 0.7075620 +0.30 +0.45 +0.46 +0.53
School Test: Run 0.0942521 0.3070050
Test: Star_r 0.1191998 0.3452533 +0.29
Test: S20_r 0.0801860 0.2831713 +0.24 +0.47
Test: SLJ 0.0623318 0.2496633 +0.40 +0.50 +0.56
Test: BMT 0.0600153 0.2449801 +0.38 +0.36 +0.24 +0.31
a1 0.0273094 0.1652554 +0.46 +0.34 +0.84 +0.76 +0.56
Sex: Girls 0.0454408 0.2131685 -0.63 -0.76 -0.83 -0.70 -0.50 -0.82
Cohort Test: Run 0.0024418 0.0494145
Test: Star_r 0.0046420 0.0681322 .
Test: S20_r 0.0020181 0.0449231 . .
Test: SLJ 0.0008351 0.0288977 . . .
Test: BMT 0.0008351 0.0288978 . . . .
Residual 0.2739686 0.5234201

4.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). This topic will be picked uo in a separate vignette. See also @Fuehner2021 for an application.

4.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.04   1.0     .      .     .
 Test: S20_r    0.11  -0.49   1.0     .     .
 Test: SLJ      0.05  -0.02  -0.41   1.0    .
 Test: BMT     -0.31   0.17  -0.17  -0.3   1.0

Normalized cumulative variances:
[0.3321, 0.6249, 0.828, 0.941, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5
 (Intercept)   -0.19  -0.44   0.7   -0.52  -0.06
 Test: Star_r   0.57   0.04   0.54   0.44   0.44
 Test: S20_r   -0.7    0.17   0.07   0.13   0.68
 Test: SLJ      0.27  -0.61  -0.46  -0.26   0.53
 Test: BMT      0.29   0.64  -0.0   -0.67   0.25, School = 
Principal components based on correlation matrix
 (Intercept)    1.0     .      .      .      .      .      .
 Test: Star_r   0.16   1.0     .      .      .      .      .
 Test: S20_r   -0.2   -0.48   1.0     .      .      .      .
 Test: SLJ     -0.03  -0.14  -0.46   1.0     .      .      .
 Test: BMT     -0.14  -0.09  -0.11  -0.28   1.0     .      .
 a1             0.81  -0.08   0.39  -0.19  -0.19   1.0     .
 Sex: Girls    -0.97  -0.18   0.09   0.23   0.18  -0.81   1.0

Normalized cumulative variances:
[0.4037, 0.6536, 0.8347, 0.9629, 1.0, 1.0, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5    PC6    PC7
 (Intercept)   -0.56   0.18  -0.01   0.22   0.07  -0.51   0.58
 Test: Star_r  -0.09   0.48   0.43  -0.56  -0.49   0.09   0.13
 Test: S20_r   -0.05  -0.71  -0.12  -0.25  -0.34   0.22   0.5
 Test: SLJ      0.14   0.41  -0.64   0.32  -0.45   0.27   0.17
 Test: BMT      0.14  -0.1    0.62   0.68  -0.32   0.13   0.08
 a1            -0.54  -0.22  -0.12   0.07  -0.49  -0.19  -0.6
 Sex: Girls     0.58  -0.09  -0.07  -0.07  -0.3   -0.74   0.0, Cohort = 
Principal components based on correlation matrix
 (Intercept)   1.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: BMT     0.0  0.0  0.0  0.0  1.0

Normalized cumulative variances:
[0.2, 0.4, 0.6, 0.8, 1.0]

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

4.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.57  1.0    .     .     .
 Test: S20_r   0.6   0.61  1.0    .     .
 Test: SLJ     0.61  0.62  0.74  1.0    .
 Test: BMT     0.3   0.45  0.46  0.53  1.0

Normalized cumulative variances:
[0.6443, 0.788, 0.8725, 0.9489, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5
 Test: Run     -0.43   0.52   0.04  -0.73  -0.12
 Test: Star_r  -0.45   0.09  -0.84   0.27   0.06
 Test: S20_r   -0.48   0.1    0.39   0.48  -0.61
 Test: SLJ     -0.49  -0.01   0.38   0.18   0.76
 Test: BMT     -0.37  -0.84  -0.02  -0.36  -0.15, School = 
Principal components based on correlation matrix
 Test: Run      1.0     .      .      .      .      .      .
 Test: Star_r   0.29   1.0     .      .      .      .      .
 Test: S20_r    0.24   0.47   1.0     .      .      .      .
 Test: SLJ      0.4    0.5    0.56   1.0     .      .      .
 Test: BMT      0.38   0.36   0.24   0.31   1.0     .      .
 a1             0.46   0.34   0.84   0.76   0.56   1.0     .
 Sex: Girls    -0.63  -0.76  -0.83  -0.7   -0.5   -0.82   1.0

Normalized cumulative variances:
[0.6069, 0.7374, 0.8431, 0.9369, 1.0, 1.0, 1.0]

Component loadings
                 PC1    PC2    PC3    PC4    PC5    PC6    PC7
 Test: Run     -0.29   0.58  -0.09   0.68  -0.19  -0.26   0.08
 Test: Star_r  -0.33  -0.05   0.84  -0.07   0.04  -0.33  -0.24
 Test: S20_r   -0.39  -0.49  -0.16  -0.08  -0.46  -0.33   0.49
 Test: SLJ     -0.39  -0.18  -0.1    0.17   0.83   0.0    0.29
 Test: BMT     -0.29   0.61  -0.07  -0.69   0.04  -0.03   0.26
 a1            -0.44  -0.11  -0.46  -0.13  -0.0   -0.15  -0.74
 Sex: Girls     0.47   0.03  -0.16  -0.1    0.23  -0.83   0.0, Cohort = 
Principal components based on correlation matrix
 Test: Run     1.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: BMT     0.0  0.0  0.0  0.0  1.0

Normalized cumulative variances:
[0.2, 0.4, 0.6, 0.8, 1.0]

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

4.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, the interactions between age and c1 and c4 contrasts and the interactions between Sex and three test contrasts (c1, c2, c4) are replicated. Obviously, the subset of data is much noisier than the full set.

5 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
Est. SE z p σ_Child σ_School σ_Cohort
(Intercept) 0.0000 0.0199 0.00 0.9982
Test: Star_r -0.0223 0.0424 -0.53 0.5991 0.7388 0.3452 0.0681
Test: S20_r 0.0188 0.0401 0.47 0.6391 0.7835 0.2830 0.0449
Test: SLJ 0.0212 0.0314 0.67 0.5003 0.7937 0.2497 0.0289
Test: BMT -0.0205 0.0321 -0.64 0.5233 0.7090 0.2449 0.0289
Test: Run & a1 0.0549 0.0768 0.71 0.4751
Test: Star_r & a1 0.1734 0.0770 2.25 0.0243
Test: S20_r & a1 0.1846 0.0784 2.35 0.0186
Test: SLJ & a1 0.0750 0.0782 0.96 0.3376
Test: BMT & a1 0.4134 0.0731 5.65 <1e-07
Test: Run & Sex: Girls -0.4920 0.0453 -10.86 <1e-26
Test: Star_r & Sex: Girls -0.1802 0.0455 -3.97 <1e-04
Test: S20_r & Sex: Girls -0.2456 0.0463 -5.30 <1e-06
Test: SLJ & Sex: Girls -0.3319 0.0465 -7.14 <1e-12
Test: BMT & Sex: Girls -0.6219 0.0434 -14.32 <1e-45
Test: Run & a1 & Sex: Girls 0.0639 0.1514 0.42 0.6732
Test: Star_r & a1 & Sex: Girls 0.1149 0.1512 0.76 0.4476
Test: S20_r & a1 & Sex: Girls 0.0391 0.1547 0.25 0.8004
Test: SLJ & a1 & Sex: Girls -0.1783 0.1549 -1.15 0.2496
Test: BMT & a1 & Sex: Girls 0.0788 0.1445 0.55 0.5856
Test: Run 0.7412 0.3070 0.0495
a1 0.1652
Sex: Girls 0.2131
Residual 0.5215

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.

5.0.1 CONSTRUCTION SITE: More 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 24467.7 24605.7 24606.7 25100.5
2 m2 69 24463.1 24601.1 24602.1 25096.0
3 m2_nested 69 24463.1 24601.1 24602.1 25096.0 
n, p, q, k = size(m1)  # nobs, fe params, VCs+CPs, re terms
(9621, 20, 13072, 3)

In principle, the models should yield the save 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.

5.0.2 Geometric degrees of freedom

From MixedModels documentation: “The sum of the leverage values is the rank of the model matrix and n - sum(leverage(m)) is the degrees of freedom for residuals. The sum of the leverage values is also the trace of the so-called”hat” matrixH.”

New term: geometric degrees of freedom.

m1_geomdf = sum(leverage(m1))  # geom_dof
5280.732584533316
sum(leverage(m2))
5530.923279710569
sum(leverage(m2_nested))
5561.244293311253
n - m1_geomdf
4340.267415466684
m1.feterm.rank
20
dof(m1)
69

6 Glossary of MixedModels.jl commands

Here we introduce most of the commands available in the MixedModels.jl package that allow the immediate inspection and analysis of results returned in a fitted linear mixed-effect model.

Postprocessing related to conditional modes will be dealt with in a different tutorial.

6.1 Overall summary statistics

+ julia> m1.optsum         # MixedModels.OptSummary:  gets all info
+ julia> loglikelihood(m1) # StatsBase.loglikelihood: return loglikelihood
                             of the model
+ julia> deviance(m1)      # StatsBase.deviance: negative twice the log-likelihood
                             relative to saturated model
+ julia> objective(m1)     # MixedModels.objective: saturated model not clear:
                             negative twice the log-likelihood
+ julia> nobs(m1)          # n of observations; they are not independent
+ julia> dof(m1)           # n of degrees of freedom is number of model parameters
+ julia> aic(m1)           # objective(m1) + 2*dof(m1)
+ julia> bic(m1)           # objective(m1) + dof(m1)*log(nobs(m1))
m1.optsum            # MixedModels.OptSummary:  gets all info
Initialization
Initial parameter vector [1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
Initial objective value 25653.67035388993
Optimizer settings
Optimizer LN_NEWUOA
Backend nlopt
ftol_rel 1.0e-12
ftol_abs 1.0e-8
xtol_rel 0.0
xtol_abs [1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10]
initial_step [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
maxfeval -1
maxtime -1.0
xtol_zero_abs 0.001
ftol_zero_abs 1.0e-5
Result
Function evaluations 2235
Final parameter vector [1.1128, 0.054, 0.1311, 0.0477, -0.4061, 1.2214, -0.5933, -0.027, 0.2404, 1.0304, -0.5068, -0.064, 0.8704, -0.4629, 1.1385, 0.3758, 0.118, -0.12, -0.0142, -0.0754, 0.2484, -0.3841, 0.708, -0.2743, -0.0645, -0.0365, -0.0673, -0.0083, 0.5194, -0.2822, -0.1037, 0.16, -0.0537, 0.3616, -0.2833, 0.0466, 0.0602, 0.4371, 0.0328, 0.0465, 0.0, 0.0, 0.0, 0.0373, 0.127, 0.0759, 0.016, 0.0582]
Final objective value 24467.6721
Return code FTOL_REACHED 
loglikelihood(m1) # StatsBase.loglikelihood: return loglikelihood of the model
-12233.83605956268
deviance(m1)      # StatsBase.deviance: negative twice the log-likelihood relative to saturated mode`
24467.67211912536
objective(m1)    # MixedModels.objective: saturated model not clear: negative twice the log-likelihood
24467.67211912536
nobs(m1) # n of observations; they are not independent
9621
n_, p_, q_, k_ = size(m1)
(9621, 20, 13072, 3)
dof(m1)  # n of degrees of freedom is number of model parameters
69
geom_df = sum(leverage(m1)) # trace of hat / rank of model matrix / geom dof
5280.732584533316
resid_df = nobs(m1) - geom_df  # eff. residual degrees of freedom
4340.267415466684
aic(m1)  # objective(m1) + 2*dof(m1)
24605.67211912536
bic(m1)  # objective(m1) + dof(m1)*log(nobs(m1))
25100.519659822286

6.2 Fixed-effect statistics

+ julia> coeftable(m1)     # StatsBase.coeftable: fixed-effects statiscs;
                             default level=0.95
+ julia> Arrow.write(joinpath(DATADIR, "m_cpx_fe.arrow"), DataFrame(coeftable(m1)));
+ julia> coef(m1)          # StatsBase.coef - parts of the table
+ julia> fixef(m1)         # MixedModels.fixef: not the same as coef()
                             for rank-deficient case
+ julia> m1.beta           # alternative extractor
+ julia> fixefnames(m1)    # works also for coefnames(m1)
+ julia> vcov(m1)          # StatsBase.vcov: var-cov matrix of fixed-effects coef.
+ julia> stderror(m1)      # StatsBase.stderror: SE for fixed-effects coefficients
+ julia> propertynames(m1) # names of available extractors
coeftable(m1) # StatsBase.coeftable: fixed-effects statiscs; default level=0.95
Coef. Std. Error z Pr(>
(Intercept) 0.000467515 0.0198427 0.02 0.9812
Test: Star_r -0.021334 0.0390506 -0.55 0.5848
Test: S20_r 0.019095 0.0323144 0.59 0.5546
Test: SLJ 0.016147 0.0259626 0.62 0.5340
Test: BMT -0.0186619 0.0308689 -0.60 0.5455
a1 0.180452 0.0542785 3.32 0.0009
Sex: Girls -0.3743 0.0328964 -11.38 <1e-29
Test: Star_r & a1 0.120144 0.0857883 1.40 0.1614
Test: S20_r & a1 0.00609658 0.083157 0.07 0.9416
Test: SLJ & a1 -0.120749 0.0766585 -1.58 0.1152
Test: BMT & a1 0.339748 0.0855176 3.97 <1e-04
Test: Star_r & Sex: Girls 0.311303 0.0493993 6.30 <1e-09
Test: S20_r & Sex: Girls -0.0633753 0.0482497 -1.31 0.1890
Test: SLJ & Sex: Girls -0.0891757 0.0447464 -1.99 0.0463
Test: BMT & Sex: Girls -0.289001 0.0497775 -5.81 <1e-08
a1 & Sex: Girls 0.0233092 0.106647 0.22 0.8270
Test: Star_r & a1 & Sex: Girls 0.0481288 0.169996 0.28 0.7771
Test: S20_r & a1 & Sex: Girls -0.0715778 0.165604 -0.43 0.6656
Test: SLJ & a1 & Sex: Girls -0.218184 0.153476 -1.42 0.1551
Test: BMT & a1 & Sex: Girls 0.256011 0.170514 1.50 0.1333 
#Arrow.write(joinpath(DATADIR, "m_cpx_fe.arrow"), DataFrame(coeftable(m1)));
coef(m1)              # StatsBase.coef; parts of the table
20-element Vector{Float64}:
  0.00046751547688959346
 -0.021334027008666566
  0.0190949541906168
  0.01614699476981148
 -0.018661932603837496
  0.18045212242925016
 -0.37430001142073505
  0.12014440464563587
  0.006096582607952155
 -0.12074929905558543
  0.3397475366907411
  0.31130293165696654
 -0.06337528930387396
 -0.08917566161354086
 -0.28900143535349726
  0.023309215568065384
  0.04812879198855719
 -0.0715777986505176
 -0.2181844069972714
  0.2560107198300933
fixef(m1)    # MixedModels.fixef: not the same as coef() for rank-deficient case
20-element Vector{Float64}:
  0.00046751547688959346
 -0.021334027008666566
  0.0190949541906168
  0.01614699476981148
 -0.018661932603837496
  0.18045212242925016
 -0.37430001142073505
  0.12014440464563587
  0.006096582607952155
 -0.12074929905558543
  0.3397475366907411
  0.31130293165696654
 -0.06337528930387396
 -0.08917566161354086
 -0.28900143535349726
  0.023309215568065384
  0.04812879198855719
 -0.0715777986505176
 -0.2181844069972714
  0.2560107198300933
m1.β                  # alternative extractor
20-element Vector{Float64}:
  0.00046751547688959346
 -0.021334027008666566
  0.0190949541906168
  0.01614699476981148
 -0.018661932603837496
  0.18045212242925016
 -0.37430001142073505
  0.12014440464563587
  0.006096582607952155
 -0.12074929905558543
  0.3397475366907411
  0.31130293165696654
 -0.06337528930387396
 -0.08917566161354086
 -0.28900143535349726
  0.023309215568065384
  0.04812879198855719
 -0.0715777986505176
 -0.2181844069972714
  0.2560107198300933
fixefnames(m1)        # works also for coefnames(m1)
20-element Vector{String}:
 "(Intercept)"
 "Test: Star_r"
 "Test: S20_r"
 "Test: SLJ"
 "Test: BMT"
 "a1"
 "Sex: Girls"
 "Test: Star_r & a1"
 "Test: S20_r & a1"
 "Test: SLJ & a1"
 "Test: BMT & a1"
 "Test: Star_r & Sex: Girls"
 "Test: S20_r & Sex: Girls"
 "Test: SLJ & Sex: Girls"
 "Test: BMT & Sex: Girls"
 "a1 & Sex: Girls"
 "Test: Star_r & a1 & Sex: Girls"
 "Test: S20_r & a1 & Sex: Girls"
 "Test: SLJ & a1 & Sex: Girls"
 "Test: BMT & a1 & Sex: Girls"
vcov(m1)   # StatsBase.vcov: var-cov matrix of fixed-effects coefficients
20×20 Matrix{Float64}:
  0.000393733   4.46622e-5   -1.51428e-5   …  -4.12578e-6   -1.85046e-5
  4.46622e-5    0.00152495   -0.000453875     -1.94887e-6    7.35688e-6
 -1.51428e-5   -0.000453875   0.00104422      -5.94438e-5   -1.0784e-5
  4.64753e-6   -4.25236e-5   -0.00034755       0.000108165  -4.5075e-5
 -9.4094e-5     1.59917e-5   -7.05281e-5      -4.47355e-5    0.000144059
 -7.23977e-5   -2.69605e-5    5.44624e-5   …   4.26813e-5   -1.4437e-5
 -0.000146004  -4.41554e-5    1.7248e-5       -3.62769e-5    0.000180536
 -1.17356e-5   -0.000403606   0.000200281     -6.20066e-6    3.32417e-5
 -1.02559e-6    0.000199515  -0.00038039      -0.000334583  -4.06852e-5
 -9.74107e-6    6.11005e-6    0.000160751      0.000562632  -0.000208963
  4.42149e-5   -2.57755e-5    2.95454e-5   …  -0.000208583   0.000564775
 -7.12781e-6   -3.05672e-5    2.04268e-5       2.29973e-5   -0.000103516
  3.03941e-6    2.04248e-5   -3.07057e-5       0.000655821   0.000123651
  4.76105e-6    8.62571e-7    8.65223e-6      -0.00132596    0.00061476
  5.90695e-6   -4.7203e-7     2.01707e-6       0.000616024  -0.00169323
  6.03417e-5    1.31915e-5   -6.25537e-6   …   0.000357947  -0.00338097
  1.15575e-5    0.000165275  -9.46663e-5      -0.000389212   0.00187177
 -2.13162e-6   -9.46511e-5    0.000161572     -0.0118661    -0.00206608
 -4.12578e-6   -1.94887e-6   -5.94438e-5       0.023555     -0.0107439
 -1.85046e-5    7.35688e-6   -1.0784e-5       -0.0107439     0.0290752
vcov(m1; corr=true) # StatsBase.vcov: correlation matrix of fixed-effects coefficients
20×20 Matrix{Float64}:
  1.0           0.0576384    -0.0236162   …  -0.00135477   -0.00546913
  0.0576384     1.0          -0.359677       -0.000325173   0.00110486
 -0.0236162    -0.359677      1.0            -0.0119859    -0.00195714
  0.00902139   -0.0419426    -0.41426         0.0271456    -0.0101819
 -0.153617      0.0132662    -0.0707041      -0.00944256    0.027369
 -0.0672196    -0.0127196     0.0310508   …   0.00512352   -0.00155987
 -0.223673     -0.0343723     0.0162254      -0.00718524    0.0321851
 -0.00689411   -0.120476      0.0722463      -0.000470943   0.00227245
 -0.000621548   0.0614398    -0.141558       -0.0262159    -0.00286931
 -0.00640392    0.00204107    0.0648931       0.0478215    -0.0159863
  0.0260563    -0.00771836    0.0106915   …  -0.0158922     0.038731
 -0.00727168   -0.0158456     0.0127963       0.0030333    -0.0122893
  0.00317463    0.0108402    -0.0196938       0.0885625     0.0150294
  0.00536221    0.00049364    0.00598376     -0.193077      0.0805725
  0.00598039   -0.000242834   0.00125398      0.080635     -0.19949
  0.0285147     0.00316753   -0.00181514  …   0.021869     -0.185923
  0.00342628    0.0248966    -0.017233       -0.0149178     0.0645732
 -0.000648693  -0.0146362     0.0301924      -0.46687      -0.0731672
 -0.00135477   -0.000325173  -0.0119859       1.0          -0.410543
 -0.00546913    0.00110486   -0.00195714     -0.410543      1.0
stderror(m1)       # StatsBase.stderror: SE for fixed-effects coefficients
20-element Vector{Float64}:
 0.01984269558668302
 0.03905055880425583
 0.03231440836412586
 0.025962560472886434
 0.030868907838576043
 0.054278508198501385
 0.03289639720659251
 0.08578825556180032
 0.08315695866689764
 0.07665847281866368
 0.08551757730309334
 0.04939928626019333
 0.04824967140335991
 0.04474637583840272
 0.04977749123952791
 0.10664666150971132
 0.16999630177557448
 0.16560395167985534
 0.15347624563699208
 0.17051441074910598
propertynames(m1)  # names of available extractors
(:formula, :reterms, :Xymat, :feterm, :sqrtwts, :parmap, :dims, :A, :L, :optsum, :θ, :theta, :β, :beta, :βs, :betas, :λ, :lambda, :stderror, :σ, :sigma, :σs, :sigmas, :σρs, :sigmarhos, :b, :u, :X, :y, :corr, :vcov, :PCA, :rePCA, :objective, :pvalues)

6.3 Covariance parameter estimates

These commands inform us about the model parameters associated with the RES.

+ julia> issingular(m1)        # Test singularity for param. vector m1.theta
+ julia> VarCorr(m1)           # MixedModels.VarCorr: est. of RES
+ julia> propertynames(m1)
+ julia> m1.σ                  # residual; or: m1.sigma
+ julia> m1.σs                 # VCs; m1.sigmas
+ julia> m1.θ                  # Parameter vector for RES (w/o residual); m1.theta
+ julia> MixedModels.sdest(m1) #  prsqrt(MixedModels.varest(m1))
+ julia> BlockDescription(m1)  #  Description of blocks of A and L in an LMM
issingular(m1) # Test if model is singular for parameter vector m1.theta (default)
true
issingular(m2)
true
VarCorr(m1) # MixedModels.VarCorr: estimates of random-effect structure (RES)
Column Variance Std.Dev Corr.
Child (Intercept) 0.36033239 0.60027693
Test: Star_r 0.43495553 0.65951159 +0.04
Test: S20_r 0.41638097 0.64527589 +0.11 -0.49
Test: SLJ 0.29608881 0.54414043 +0.05 -0.02 -0.41
Test: BMT 0.50554575 0.71101741 -0.31 +0.17 -0.17 -0.30
School (Intercept) 0.04109627 0.20272215
Test: Star_r 0.14992416 0.38720041 +0.16
Test: S20_r 0.10460004 0.32341929 -0.20 -0.48
Test: SLJ 0.06249852 0.24999704 -0.03 -0.14 -0.46
Test: BMT 0.08414285 0.29007387 -0.14 -0.09 -0.11 -0.28
a1 0.02766141 0.16631719 +0.81 -0.08 +0.39 -0.19 -0.19
Sex: Girls 0.04547115 0.21323966 -0.97 -0.18 +0.09 +0.23 +0.18 -0.81
Cohort (Intercept) 0.00040530 0.02013200
Test: Star_r 0.00469459 0.06851710 .
Test: S20_r 0.00167603 0.04093934 . .
Test: SLJ 0.00007426 0.00861717 . . .
Test: BMT 0.00098435 0.03137434 . . . .
Residual 0.29100268 0.53944665 
VarCorr(m2)
Column Variance Std.Dev Corr.
Child Test: Run 0.5473394 0.7398239
Test: Star_r 0.5437814 0.7374153 +0.57
Test: S20_r 0.6118591 0.7822142 +0.60 +0.61
Test: SLJ 0.6278340 0.7923598 +0.61 +0.62 +0.74
Test: BMT 0.5006439 0.7075620 +0.30 +0.45 +0.46 +0.53
School Test: Run 0.0942521 0.3070050
Test: Star_r 0.1191998 0.3452533 +0.29
Test: S20_r 0.0801860 0.2831713 +0.24 +0.47
Test: SLJ 0.0623318 0.2496633 +0.40 +0.50 +0.56
Test: BMT 0.0600153 0.2449801 +0.38 +0.36 +0.24 +0.31
a1 0.0273094 0.1652554 +0.46 +0.34 +0.84 +0.76 +0.56
Sex: Girls 0.0454408 0.2131685 -0.63 -0.76 -0.83 -0.70 -0.50 -0.82
Cohort Test: Run 0.0024418 0.0494145
Test: Star_r 0.0046420 0.0681322 .
Test: S20_r 0.0020181 0.0449231 . .
Test: SLJ 0.0008351 0.0288977 . . .
Test: BMT 0.0008351 0.0288978 . . . .
Residual 0.2739686 0.5234201 
m1.σs      # VCs; m1.sigmas
(Child = (var"(Intercept)" = 0.6002769314424566, var"Test: Star_r" = 0.6595115868256303, var"Test: S20_r" = 0.6452758896278719, var"Test: SLJ" = 0.544140434624562, var"Test: BMT" = 0.7110174077935371), School = (var"(Intercept)" = 0.20272214895029594, var"Test: Star_r" = 0.38720041335022076, var"Test: S20_r" = 0.3234192913062064, var"Test: SLJ" = 0.24999703867581158, var"Test: BMT" = 0.2900738686545868, a1 = 0.16631718564740655, var"Sex: Girls" = 0.21323966128257757), Cohort = (var"(Intercept)" = 0.020132001586338804, var"Test: Star_r" = 0.06851709691438108, var"Test: S20_r" = 0.040939340943969844, var"Test: SLJ" = 0.008617170012474924, var"Test: BMT" = 0.031374342939087745))
m1.θ       # Parameter vector for RES (w/o residual); m1.theta
48-element Vector{Float64}:
  1.1127642275182457
  0.05395891146748544
  0.13112205732461346
  0.04774927196860498
 -0.40605539183625283
  1.2213792198768771
 -0.5932956280182006
 -0.02701029907573906
  0.24044471057717898
  1.0303671940424057
  ⋮
  0.046475648212836077
  0.0
  2.22044225634205e-5
  0.0
  0.03731972697965639
  0.1270136672360744
  0.07589136232697247
  0.01597409134028042
  0.05816023347859276
BlockDescription(m1) #  Description of blocks of A and L in a LinearMixedModel
rows Child School Cohort fixed
9905 BlkDiag
3122 Sparse BlkDiag
45 Dense Dense BlkDiag/Dense
21 Dense Dense Dense Dense 
m2.θ
48-element Vector{Float64}:
  1.41344188645469
  0.8064340384256612
  0.892792599664444
  0.9296770066479014
  0.4037800334352948
  1.1552034695499855
  0.49755286787692177
  0.4892487589421097
  0.4583871856940061
  1.0902663878832757
  ⋮
 -0.03799157772996963
  0.0
  8.187819876586734e-5
  0.0
  0.09440702228233282
  0.1301673073369134
  0.08582599342290523
  0.05520935784694587
  0.05520952065294796
BlockDescription(m2)
rows Child School Cohort fixed
9905 BlkDiag
3122 Sparse BlkDiag
45 Dense Dense BlkDiag/Dense
21 Dense Dense Dense Dense

6.4 Model “predictions”

These commands inform us about extracion of conditional modes/means and (co-)variances, that using the model parameters to improve the predictions for units (levels) of the grouping (random) factors. We need this information, e.g., for partial-effect response profiles (e.g., facet plot) or effect profiles (e.g., caterpillar plot), or visualizing the borrowing-strength effect for correlation parameters (e.g., shrinkage plots). We are using the fit of LMM m2.

julia> condVar(m2)

Some plotting functions are currently available from the MixedModelsMakie package or via custom functions.

+ julia> caterpillar!(m2)
+ julia> shrinkage!(m2)

6.4.1 Conditional covariances

condVar(m1)
3-element Vector{Array{Float64, 3}}:
 [0.07420163921893243 0.006246631778172274 … 0.002029974984645099 -0.018290105835260095; 0.006246631778172274 0.26043146788912047 … -0.005987019207423662 0.026504216716524438; … ; 0.002029974984645099 -0.005987019207423662 … 0.20011536387900336 -0.08601644897554256; -0.018290105835260095 0.026504216716524438 … -0.08601644897554256 0.26191506839588446;;; 0.05144496548725417 0.003694319926726746 … 0.0033961011026631942 -0.010357627316085376; 0.003694319926726746 0.2574677448144253 … -0.005683277790987534 0.026938674362014314; … ; 0.0033961011026631942 -0.005683277790987534 … 0.19952797213889106 -0.08609034682307835; -0.010357627316085376 0.026938674362014314 … -0.08609034682307835 0.257879510631094;;; 0.0753764989530672 0.006570848768885153 … 0.0019709649686385543 -0.01879606522494719; 0.006570848768885153 0.26160064870616845 … -0.006100648757492659 0.026489728945137434; … ; 0.0019709649686385543 -0.006100648757492659 … 0.20040025807171596 -0.08599494635458116; -0.01879606522494719 0.026489728945137434 … -0.08599494635458116 0.26276039307469234;;; … ;;; 0.06619880394307251 0.0028635468566462915 … 0.002739217336934494 -0.016052700133800477; 0.0028635468566462915 0.26128470414043836 … -0.006196270147891165 0.02728909944872443; … ; 0.002739217336934494 -0.006196270147891165 … 0.20046587194493404 -0.08613038325444577; -0.016052700133800477 0.02728909944872443 … -0.08613038325444577 0.2620010285506084;;; 0.08471643183252144 0.0024597727271779197 … -0.031144245119574573 -0.026177577727926237; 0.0024597727271779197 0.2613407748148688 … -0.005364849865599808 0.027521641515535492; … ; -0.031144245119574573 -0.005364849865599808 … 0.2643156889036811 -0.06737386544295736; -0.026177577727926237 0.027521641515535492 … -0.06737386544295736 0.26758298721340557;;; 0.0645514376778271 0.0020081139409852005 … 0.0028725198439486267 -0.015376994673724964; 0.0020081139409852005 0.26036753622472153 … -0.00609932885097756 0.027443859404963595; … ; 0.0028725198439486267 -0.00609932885097756 … 0.19995986608139413 -0.0861423441336708; -0.015376994673724964 0.027443859404963595 … -0.0861423441336708 0.2605800115340165]
 [0.027811603793077684 0.007395790924316572 … 0.019039480001574013 -0.02822880611822501; 0.007395790924316572 0.11441775923746687 … -0.006148291052989946 -0.008363339110063591; … ; 0.019039480001574013 -0.006148291052989946 … 0.021137242510252122 -0.020166463503234925; -0.02822880611822501 -0.008363339110063591 … -0.020166463503234925 0.030844469281447663;;; 0.02868078727049485 0.007319064305100705 … 0.019823392359319778 -0.029134679773970967; 0.007319064305100705 0.1143812655003623 … -0.006176882626156598 -0.008286707847374413; … ; 0.019823392359319778 -0.006176882626156598 … 0.02180846593512761 -0.0209805457084117; -0.029134679773970967 -0.008286707847374413 … -0.0209805457084117 0.03178834538508992;;; 0.038947689195447485 0.008341508790463355 … 0.027004809862487246 -0.039701420740468854; 0.008341508790463355 0.10289838690461442 … -0.0044550294925543845 -0.009095139923081256; … ; 0.027004809862487246 -0.0044550294925543845 … 0.02611641812385842 -0.02822587191307009; -0.039701420740468854 -0.009095139923081256 … -0.02822587191307009 0.042497578355228216;;; … ;;; 0.030598791396820296 0.009646157087130803 … 0.020093390341393725 -0.03109460326672016; 0.009646157087130803 0.1298034468672959 … -0.005987298303009213 -0.011018485660283227; … ; 0.020093390341393725 -0.005987298303009213 … 0.022012186152513758 -0.02136061340602325; -0.03109460326672016 -0.011018485660283227 … -0.02136061340602325 0.03394957110692568;;; 0.029703900102630822 0.008199735889556858 … 0.020146116313597844 -0.030245778061352847; 0.008199735889556858 0.104311848330787 … -0.004646662898071603 -0.008851636578050008; … ; 0.020146116313597844 -0.004646662898071603 … 0.020937813097195196 -0.021188221813950404; -0.030245778061352847 -0.008851636578050008 … -0.021188221813950404 0.03271774263789534;;; 0.029819299059823258 0.007835885042132843 … 0.019750617610598174 -0.03029789835952142; 0.007835885042132843 0.12307137838241447 … -0.005604242239536485 -0.008854426461858472; … ; 0.019750617610598174 -0.005604242239536485 … 0.02103331926706263 -0.020826379826384786; -0.03029789835952142 -0.008854426461858472 … -0.020826379826384786 0.03283125944668525]
 [0.00035704410831846416 1.1212052697180321e-5 … -2.1942430054311961e-7 -1.3288405553777801e-5; 1.1212052697180321e-5 0.002658670484631714 … -1.0434384526615128e-5 -3.4354028367291965e-5; … ; -2.1942430054311961e-7 -1.0434384526615128e-5 … 7.293192949610833e-5 -6.260181209342795e-6; -1.3288405553777801e-5 -3.4354028367291965e-5 … -6.260181209342795e-6 0.0008450015327209852;;; 0.00034794189429652913 1.266005013244718e-5 … -2.427911356759967e-7 -1.5087455257657452e-5; 1.266005013244718e-5 0.002421370820688978 … -1.0909989442265396e-5 -3.23976395471745e-5; … ; -2.427911356759967e-7 -1.0909989442265396e-5 … 7.267482543124194e-5 -7.1743923802881e-6; -1.5087455257657452e-5 -3.23976395471745e-5 … -7.1743923802881e-6 0.0008206252884776522;;; 0.0003402934393195578 1.4982542049551411e-5 … -2.625576477517161e-7 -1.6775125293400984e-5; 1.4982542049551411e-5 0.002225268841481532 … -1.1136439969200303e-5 -3.13152825361867e-5; … ; -2.625576477517161e-7 -1.1136439969200303e-5 … 7.24416211205863e-5 -8.143291337648845e-6; -1.6775125293400984e-5 -3.13152825361867e-5 … -8.143291337648845e-6 0.0007976234587123457;;; 0.00034882738856309393 1.3842807974195054e-5 … -2.5449397333559213e-7 -1.5075973310685458e-5; 1.3842807974195054e-5 0.0024234748666635427 … -1.0918330040410884e-5 -3.334030538094929e-5; … ; -2.5449397333559213e-7 -1.0918330040410884e-5 … 7.268739896328334e-5 -7.201100251499824e-6; -1.5075973310685458e-5 -3.334030538094929e-5 … -7.201100251499824e-6 0.0008215628344350112;;; 0.0003412330571835013 1.4823676070814796e-5 … -2.9123108025869365e-7 -1.6829310786925892e-5; 1.4823676070814796e-5 0.0022591987984080276 … -1.1310499965088666e-5 -3.211970030518062e-5; … ; -2.9123108025869365e-7 -1.1310499965088666e-5 … 7.2476504747266e-5 -7.968379172209978e-6; -1.6829310786925892e-5 -3.211970030518062e-5 … -7.968379172209978e-6 0.000800092298364764;;; 0.0003391889467495915 1.5193724915651492e-5 … -2.781787542249665e-7 -1.68866712626123e-5; 1.5193724915651492e-5 0.0022195318187253235 … -1.1248308096525554e-5 -3.260676280234238e-5; … ; -2.781787542249665e-7 -1.1248308096525554e-5 … 7.241280922972752e-5 -8.202870499260977e-6; -1.68866712626123e-5 -3.260676280234238e-5 … -8.202870499260977e-6 0.0007949422455376611;;; 0.0003361090575844702 1.3707337112862039e-5 … -3.6842198296176886e-7 -1.776155976899562e-5; 1.3707337112862039e-5 0.0021638978611889744 … -1.1487671655525522e-5 -3.151598600376191e-5; … ; -3.6842198296176886e-7 -1.1487671655525522e-5 … 7.229485640743411e-5 -8.47917466700245e-6; -1.776155976899562e-5 -3.151598600376191e-5 … -8.47917466700245e-6 0.0007840481029990294;;; 0.00032913230669072825 1.426259418977282e-5 … -3.296329226388193e-7 -1.9023843284094874e-5; 1.426259418977282e-5 0.0020421421447891176 … -1.1386299620017928e-5 -3.0247998699223892e-5; … ; -3.296329226388193e-7 -1.1386299620017928e-5 … 7.211872930867235e-5 -9.185273148195754e-6; -1.9023843284094874e-5 -3.0247998699223892e-5 … -9.185273148195754e-6 0.0007672810913842903;;; 0.00031684584216499535 1.3355567397232205e-5 … -4.029066541633233e-7 -2.1493287429899596e-5; 1.3355567397232205e-5 0.0018858548786504388 … -1.130597890783126e-5 -2.764848712954746e-5; … ; -4.029066541633233e-7 -1.130597890783126e-5 … 7.17740770975573e-5 -1.0280687427020646e-5; -2.1493287429899596e-5 -2.764848712954746e-5 … -1.0280687427020646e-5 0.0007361891142169222]
condVar(m2)
3-element Vector{Array{Float64, 3}}:
 [0.1834791428369631 0.05174737225101482 … 0.05261696011012374 0.01404003715888325; 0.05174737225101482 0.18630739402998828 … 0.051946731289680456 0.03660971225113556; … ; 0.05261696011012374 0.051946731289680456 … 0.17231734708116886 0.04183013762478513; 0.01404003715888325 0.03660971225113556 … 0.04183013762478513 0.17838984359022692;;; 0.15980391179377965 0.02720873837569012 … 0.028491756732142177 -0.0022054210124423665; 0.02720873837569012 0.15750228845387587 … 0.02525807667821121 0.01821799795731268; … ; 0.028491756732142177 0.02525807667821121 … 0.14588684572930063 0.02454642361508862; -0.0022054210124423665 0.01821799795731268 … 0.02454642361508862 0.16608795132658638;;; 0.18503798532084798 0.05291886499198959 … 0.053770608279896576 0.014670749710339107; 0.05291886499198959 0.18845601356096237 … 0.053432608173227654 0.03754568156770724; … ; 0.053770608279896576 0.053432608173227654 … 0.17384706106001702 0.04255782512863491; 0.014670749710339107 0.03754568156770724 … 0.04255782512863491 0.1792514647513263;;; … ;;; 0.17317836267933334 0.037634552579831865 … 0.044333919520989834 0.00833296891166478; 0.037634552579831865 0.16945162166087177 … 0.04018105657000757 0.0281210365654034; … ; 0.044333919520989834 0.04018105657000757 … 0.16782133959500503 0.03923170566557307; 0.00833296891166478 0.0281210365654034 … 0.03923170566557307 0.17783782794987718;;; 0.17860593954705672 0.04287942523344259 … 0.05257564759714597 0.011535010795276177; 0.04287942523344259 0.17459623987980966 … 0.04819770495307883 0.031262619605128765; … ; 0.05257564759714597 0.04819770495307883 … 0.18049406026016496 0.04406528271747994; 0.011535010795276177 0.031262619605128765 … 0.04406528271747994 0.17977147300729193;;; 0.17225965233499474 0.03594900139945862 … 0.04287506313438603 0.007507595071675819; 0.03594900139945862 0.16595033134316514 … 0.03804567167310316 0.026785223842648642; … ; 0.04287506313438603 0.03804567167310316 … 0.16549729351146875 0.03810923157602961; 0.007507595071675819 0.026785223842648642 … 0.03810923157602961 0.17633162284517606]
 [0.06935668085706434 0.018481973397403654 … 0.01565968329448124 -0.027793659524471534; 0.018481973397403654 0.0826871356417127 … 0.010530640116219552 -0.035890748826687684; … ; 0.01565968329448124 0.010530640116219552 … 0.020769734385725113 -0.02028680094944356; -0.027793659524471534 -0.035890748826687684 … -0.02028680094944356 0.030808606228451847;;; 0.07011654765239465 0.019170532910597247 … 0.01641046233424289 -0.028651092880002952; 0.019170532910597247 0.08327505672349995 … 0.011268453945737536 -0.03669726685222868; … ; 0.01641046233424289 0.011268453945737536 … 0.021431736409315786 -0.021095321273622052; -0.028651092880002952 -0.03669726685222868 … -0.021095321273622052 0.03175794660811435;;; 0.07817151237540611 0.0343623875143453 … 0.024485988033850343 -0.04035255476152567; 0.0343623875143453 0.09401197851467276 … 0.0209665432739896 -0.04919257609903791; … ; 0.024485988033850343 0.0209665432739896 … 0.025830752058229284 -0.028369632976778667; -0.04035255476152567 -0.04919257609903791 … -0.028369632976778667 0.04244355720620964;;; … ;;; 0.07750248386416694 0.02122757646140273 … 0.01668293642642119 -0.030576843614035013; 0.02122757646140273 0.0955623500419963 … 0.011763219903347076 -0.04123387172620302; … ; 0.01668293642642119 0.011763219903347076 … 0.021665131157143615 -0.021523476524858958; -0.030576843614035013 -0.04123387172620302 … -0.021523476524858958 0.03389235057072086;;; 0.06777449943166543 0.022831506945702126 … 0.01722260045116368 -0.03014855257740921; 0.022831506945702126 0.08282763529957686 … 0.013500187082547027 -0.03869830786036218; … ; 0.01722260045116368 0.013500187082547027 … 0.02064856673400251 -0.02132328274557488; -0.03014855257740921 -0.03869830786036218 … -0.02132328274557488 0.032644659950718954;;; 0.07746534999777382 0.02107342582926149 … 0.01686630049312837 -0.030813310109198154; 0.02107342582926149 0.08840261207099885 … 0.012312851155844396 -0.03932872077841236; … ; 0.01686630049312837 0.012312851155844396 … 0.020717777383905844 -0.02097042017642257; -0.030813310109198154 -0.03932872077841236 … -0.02097042017642257 0.03276270896719467]
 [0.0016565649252985144 0.00016050848355271747 … 6.262221923322494e-5 1.0124287356148408e-5; 0.00016050848355271747 0.0023895743217316008 … 8.118463776180893e-5 5.431295226442004e-5; … ; 6.262221923322494e-5 8.118463776180893e-5 … 0.0007076908963497845 2.4071539718671977e-5; 1.0124287356148408e-5 5.431295226442004e-5 … 2.4071539718671977e-5 0.0007142976843305117;;; 0.0015390715131331555 0.0001725747015696921 … 7.239690106357643e-5 1.2424453140698236e-5; 0.0001725747015696921 0.0021338679607886407 … 9.115237981142446e-5 5.9582468762130715e-5; … ; 7.239690106357643e-5 9.115237981142446e-5 … 0.000684543232603804 2.837459875753945e-5; 1.2424453140698236e-5 5.9582468762130715e-5 … 2.837459875753945e-5 0.0006931403140804265;;; 0.0014349945412343407 0.00018155721012576398 … 7.962493985415916e-5 1.590229264818653e-5; 0.00018155721012576398 0.001967922088781173 … 9.64450968129515e-5 6.477264665392986e-5; … ; 7.962493985415916e-5 9.64450968129515e-5 … 0.0006664962829353928 3.1834104878137954e-5; 1.590229264818653e-5 6.477264665392986e-5 … 3.1834104878137954e-5 0.0006732836828646648;;; 0.0015387250780588787 0.0001739761997975363 … 7.16423816185037e-5 1.3018108701845242e-5; 0.0001739761997975363 0.002161420191039463 … 8.951405112555336e-5 5.9963847128235556e-5; … ; 7.16423816185037e-5 8.951405112555336e-5 … 0.000686625352562016 2.7929885351398195e-5; 1.3018108701845242e-5 5.9963847128235556e-5 … 2.7929885351398195e-5 0.0006940008438064902;;; 0.0014504904161016922 0.00017640111957789784 … 7.91782105423148e-5 1.5235674377077875e-5; 0.00017640111957789784 0.002007881862892075 … 9.544420082310591e-5 6.296921318909321e-5; … ; 7.91782105423148e-5 9.544420082310591e-5 … 0.0006680681117790557 3.181315980194258e-5; 1.5235674377077875e-5 6.296921318909321e-5 … 3.181315980194258e-5 0.0006752863386655104;;; 0.00142878589505435 0.00017637743004036936 … 7.951838882253712e-5 1.63838135583607e-5; 0.00017637743004036936 0.001957854565490459 … 9.712270281491905e-5 6.375803397253176e-5; … ; 7.951838882253712e-5 9.712270281491905e-5 … 0.0006635813769283605 3.249978763646735e-5; 1.63838135583607e-5 6.375803397253176e-5 … 3.249978763646735e-5 0.000670998947884876;;; 0.0014013081103788402 0.00017836309560561477 … 8.381663429958933e-5 1.77259805203796e-5; 0.00017836309560561477 0.001891982408659417 … 0.00010189356736687529 6.599072853632253e-5; … ; 8.381663429958933e-5 0.00010189356736687529 … 0.0006519970532222221 3.47505456283251e-5; 1.77259805203796e-5 6.599072853632253e-5 … 3.47505456283251e-5 0.0006616506237241126;;; 0.0013333105826927895 0.00017794082592507943 … 8.745195622868091e-5 1.9497887182972066e-5; 0.00017794082592507943 0.0017604464123340464 … 0.00010494324071326891 6.788628914306501e-5; … ; 8.745195622868091e-5 0.00010494324071326891 … 0.0006396773495634702 3.687143530149366e-5; 1.9497887182972066e-5 6.788628914306501e-5 … 3.687143530149366e-5 0.0006467079489776267;;; 0.0012449400420847963 0.00017747389102892603 … 9.255396183293524e-5 2.171182279942239e-5; 0.00017747389102892603 0.0015951237959942373 … 0.00011026456766466698 7.129656373709772e-5; … ; 9.255396183293524e-5 0.00011026456766466698 … 0.000613716350052961 4.104981014548895e-5; 2.171182279942239e-5 7.129656373709772e-5 … 4.104981014548895e-5 0.0006188680965184902]

They are hard to look at. Let’s take pictures.

6.4.2 Caterpillar plots

Code 
caterpillar!(
  Figure(; size=(800, 400)), ranefinfo(m1, :Cohort)
)
Figure 3: Prediction intervals of the random effects for Cohort in model m1

6.4.3 Shrinkage plots

These are just teasers. We will pick this up in a separate tutorial. Enjoy!

Code 
shrinkageplot!(Figure(; size=(800, 800)), m1, :Cohort)
Figure 4: Shrinkage plot of the random effects for Cohort in model m1
Code 
shrinkageplot!(Figure(; size=(800, 800)), m2, :Cohort)
Figure 5: Shrinkage plot of the random effects for Cohort in model m2

6.5 See Also

This page was rendered from git revision 0e399c4 using Quarto 1.9.36.

Back to top