Code
+using Arrow
+using AlgebraOfGraphics
+using CairoMakie # for displaying static plots
+using DataFrames
+using Statistics
+using StatsBase
+using SMLP2023: dataset
+
+activate!(; type="svg") # use SVG (other options include PNG) CairoMakie.
diff --git a/.nojekyll b/.nojekyll index befa032..233c9d0 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -7fe3b56d \ No newline at end of file +671831f7 \ No newline at end of file diff --git a/AoGPlots.html b/AoGPlots.html new file mode 100644 index 0000000..f130cd8 --- /dev/null +++ b/AoGPlots.html @@ -0,0 +1,1072 @@ + +
+ + + + + + + + + +Phillip Alday, Douglas Bates, and Reinhold Kliegl
+2023-08-21
+This notebook shows creating a multi-panel plot similar to Figure 2 of Fühner et al. (2021).
+The data are available from the SMLP2023 example datasets.
+Arrow.Table with 525126 rows, 7 columns, and schema:
+ :Cohort String
+ :School String
+ :Child String
+ :Sex String
+ :age Float64
+ :Test String
+ :score Float64
+The response to be plotted is the mean score by Test
and Sex
and age
, rounded to the nearest 0.1 years.
The first task is to round the age
to 1 digit after the decimal place, which can be done with select
applied to a DataFrame
. In some ways this is the most complicated expression in creating the plot so we will break it down. select
is applied to DataFrame(dat)
, which is the conversion of the Arrow.Table
, dat
, to a DataFrame
. This is necessary because an Arrow.Table
is immutable but a DataFrame
can be modified.
The arguments after the DataFrame
describe how to modify the contents. The first :
indicates that all the existing columns should be included. The other expression can be pairs (created with the =>
operator) of the form :col => function
or of the form :col => function => :newname
. (See the documentation of the DataFrames package for details.)
In this case the function is an anonymous function of the form round.(x, digits=1)
where “dot-broadcasting” is used to apply to the entire column (see this documentation for details).
transform!(df, :age, :age => (x -> x .- 8.5) => :a1) # centered age (linear)
+select!(groupby(df, :Test), :, :score => zscore => :zScore) # z-score
+tlabels = [ # establish order and labels of tbl.Test
+ "Run" => "Endurance",
+ "Star_r" => "Coordination",
+ "S20_r" => "Speed",
+ "SLJ" => "PowerLOW",
+ "BPT" => "PowerUP",
+];
The next stage is a group-apply-combine operation to group the rows by Sex
, Test
and rnd_age
then apply mean
to the zScore
and also apply length
to zScore
to record the number in each group.
df2 = combine(
+ groupby(
+ select(df, :, :age => ByRow(x -> round(x; digits=1)) => :age),
+ [:Sex, :Test, :age],
+ ),
+ :zScore => mean => :zScore,
+ :zScore => length => :n,
+)
Row | +Sex | +Test | +age | +zScore | +n | +
---|---|---|---|---|---|
+ | String | +String | +Float64 | +Float64 | +Int64 | +
1 | +male | +S20_r | +8.0 | +-0.0265138 | +1223 | +
2 | +male | +BPT | +8.0 | +0.026973 | +1227 | +
3 | +male | +SLJ | +8.0 | +0.121609 | +1227 | +
4 | +male | +Star_r | +8.0 | +-0.0571726 | +1186 | +
5 | +male | +Run | +8.0 | +0.292695 | +1210 | +
6 | +female | +S20_r | +8.0 | +-0.35164 | +1411 | +
7 | +female | +BPT | +8.0 | +-0.610355 | +1417 | +
8 | +female | +SLJ | +8.0 | +-0.279872 | +1418 | +
9 | +female | +Star_r | +8.0 | +-0.268221 | +1381 | +
10 | +female | +Run | +8.0 | +-0.245573 | +1387 | +
11 | +male | +S20_r | +8.1 | +0.0608397 | +3042 | +
12 | +male | +BPT | +8.1 | +0.0955413 | +3069 | +
13 | +male | +SLJ | +8.1 | +0.123099 | +3069 | +
⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +
109 | +male | +Star_r | +9.0 | +0.254973 | +4049 | +
110 | +male | +Run | +9.0 | +0.258082 | +4034 | +
111 | +female | +S20_r | +9.1 | +-0.0286172 | +1154 | +
112 | +female | +BPT | +9.1 | +-0.0752301 | +1186 | +
113 | +female | +SLJ | +9.1 | +-0.094587 | +1174 | +
114 | +female | +Star_r | +9.1 | +0.00276252 | +1162 | +
115 | +female | +Run | +9.1 | +-0.235591 | +1150 | +
116 | +male | +S20_r | +9.1 | +0.325745 | +1303 | +
117 | +male | +BPT | +9.1 | +0.616416 | +1320 | +
118 | +male | +SLJ | +9.1 | +0.267577 | +1310 | +
119 | +male | +Star_r | +9.1 | +0.254342 | +1297 | +
120 | +male | +Run | +9.1 | +0.251045 | +1294 | +
The AlgebraOfGraphics
package applies operators to the results of functions such as data
(specify the data table to be used), mapping
(designate the roles of columns), and visual
(type of visual presentation).
let
+ design = mapping(:age, :zScore; color=:Sex, col=:Test)
+ lines = design * linear()
+ means = design * visual(Scatter; markersize=5)
+ draw(data(df2) * means + data(df) * lines)
+end
Phillip Alday, Douglas Bates, and Reinhold Kliegl
+2023-08-21
+In Fühner et al. (2021) the original metric of two tasks (Star, S20) is time, but they were transformed to speed scores in the publication prior to computing z-scores. The critical result is the absence of evidence for the age x Sex x Test interaction. Is this interaction significant if we analyse all tasks in their original metric?
+Fitting the LMM of the publication takes time, roughly 1 hour. However, if you save the model parameters (and other relevant information), you can restore the fitted model object very quickly. The notebook also illustrates this procedure.
+dat = DataFrame(Arrow.Table(joinpath(datadir, "fggk21.arrow")))
+@transform!(dat, :a1 = :age - 8.5);
+select!(groupby(dat, :Test), :, :score => zscore => :zScore);
+describe(dat)
Row | +variable | +mean | +min | +median | +max | +nmissing | +eltype | +
---|---|---|---|---|---|---|---|
+ | Symbol | +Union… | +Any | +Union… | +Any | +Int64 | +DataType | +
1 | +Test | ++ | BPT | ++ | Star_r | +0 | +String | +
2 | +Cohort | ++ | 2011 | ++ | 2019 | +0 | +String | +
3 | +School | ++ | S100043 | ++ | S800200 | +0 | +String | +
4 | +Child | ++ | C002352 | ++ | C117966 | +0 | +String | +
5 | +Sex | ++ | female | ++ | male | +0 | +String | +
6 | +age | +8.56073 | +7.99452 | +8.55852 | +9.10609 | +0 | +Float64 | +
7 | +score | +226.141 | +1.14152 | +4.65116 | +1530.0 | +0 | +Float64 | +
8 | +a1 | +0.0607297 | +-0.505476 | +0.0585216 | +0.606092 | +0 | +Float64 | +
9 | +zScore | +-3.91914e-13 | +-3.1542 | +0.00031088 | +3.55078 | +0 | +Float64 | +
+ | Est. | +SE | +z | +p | +σ_Child | +σ_School | +σ_Cohort | +
---|---|---|---|---|---|---|---|
(Intercept) | +-0.0383 | +0.0108 | +-3.56 | +0.0004 | +0.5939 | +0.2024 | +0.0157 | +
Test: Run | +-0.0228 | +0.0274 | +-0.83 | +0.4052 | +0.8384 | +0.3588 | +0.0651 | +
Test: S20_r | +-0.0147 | +0.0405 | +-0.36 | +0.7171 | +0.5825 | +0.3596 | +0.1107 | +
Test: SLJ | +0.0328 | +0.0330 | +0.99 | +0.3198 | +0.4127 | +0.3027 | +0.0896 | +
Test: Star_r | +0.0006 | +0.0197 | +0.03 | +0.9763 | +0.5574 | +0.3620 | +0.0313 | +
a1 | +0.2713 | +0.0086 | +31.63 | +<1e-99 | ++ | 0.0966 | ++ |
Sex: male | +0.2064 | +0.0024 | +86.55 | +<1e-99 | ++ | 0.0245 | ++ |
Test: Run & a1 | +-0.4464 | +0.0131 | +-34.05 | +<1e-99 | ++ | + | + |
Test: S20_r & a1 | +0.1473 | +0.0114 | +12.97 | +<1e-37 | ++ | + | + |
Test: SLJ & a1 | +-0.0068 | +0.0103 | +-0.66 | +0.5116 | ++ | + | + |
Test: Star_r & a1 | +0.0761 | +0.0111 | +6.84 | +<1e-11 | ++ | + | + |
Test: Run & Sex: male | +-0.0900 | +0.0037 | +-24.10 | +<1e-99 | ++ | + | + |
Test: S20_r & Sex: male | +-0.0912 | +0.0032 | +-28.23 | +<1e-99 | ++ | + | + |
Test: SLJ & Sex: male | +0.0330 | +0.0029 | +11.24 | +<1e-28 | ++ | + | + |
Test: Star_r & Sex: male | +-0.0720 | +0.0032 | +-22.65 | +<1e-99 | ++ | + | + |
a1 & Sex: male | +0.0010 | +0.0069 | +0.14 | +0.8876 | ++ | + | + |
Test: Run & a1 & Sex: male | +-0.0154 | +0.0126 | +-1.22 | +0.2233 | ++ | + | + |
Test: S20_r & a1 & Sex: male | +0.0129 | +0.0109 | +1.18 | +0.2380 | ++ | + | + |
Test: SLJ & a1 & Sex: male | +-0.0098 | +0.0100 | +-0.98 | +0.3256 | ++ | + | + |
Test: Star_r & a1 & Sex: male | +0.0166 | +0.0108 | +1.54 | +0.1241 | ++ | + | + |
Residual | +0.5880 | ++ | + | + | + | + | + |
+ | Column | +Variance | +Std.Dev | +Corr. | ++ | + | + | + | + |
---|---|---|---|---|---|---|---|---|---|
Child | +(Intercept) | +0.3527294 | +0.5939103 | ++ | + | + | + | + | + |
+ | Test: Run | +0.7029003 | +0.8383915 | ++0.11 | ++ | + | + | + | + |
+ | Test: S20_r | +0.3393356 | +0.5825252 | ++0.19 | +-0.53 | ++ | + | + | + |
+ | Test: SLJ | +0.1702900 | +0.4126621 | ++0.05 | +-0.14 | +-0.29 | ++ | + | + |
+ | Test: Star_r | +0.3107227 | +0.5574251 | +-0.10 | ++0.01 | +-0.13 | +-0.42 | ++ | + |
School | +(Intercept) | +0.0409640 | +0.2023957 | ++ | + | + | + | + | + |
+ | Test: Run | +0.1287690 | +0.3588440 | ++0.26 | ++ | + | + | + | + |
+ | Test: S20_r | +0.1293351 | +0.3596319 | ++0.01 | +-0.57 | ++ | + | + | + |
+ | Test: SLJ | +0.0916522 | +0.3027411 | +-0.13 | ++0.01 | +-0.53 | ++ | + | + |
+ | Test: Star_r | +0.1310575 | +0.3620187 | ++0.26 | ++0.09 | +-0.06 | +-0.28 | ++ | + |
+ | a1 | +0.0093412 | +0.0966499 | ++0.48 | ++0.25 | +-0.15 | +-0.01 | ++0.12 | ++ |
+ | Sex: male | +0.0005999 | +0.0244934 | ++0.09 | ++0.13 | +-0.01 | ++0.05 | +-0.19 | ++0.25 | +
Cohort | +(Intercept) | +0.0002452 | +0.0156587 | ++ | + | + | + | + | + |
+ | Test: Run | +0.0042389 | +0.0651068 | +. | ++ | + | + | + | + |
+ | Test: S20_r | +0.0122535 | +0.1106954 | +. | +. | ++ | + | + | + |
+ | Test: SLJ | +0.0080210 | +0.0895599 | +. | +. | +. | ++ | + | + |
+ | Test: Star_r | +0.0009828 | +0.0313498 | +. | +. | +. | +. | ++ | + |
Residual | ++ | 0.3456872 | +0.5879517 | ++ | + | + | + | + | + |
Residual plots for published LMM
+ + +Residual plots for LMM with Star and Speed in original metric.
+ + +m_ovi
m_zcp
m_cpx
m1
to the reduced dataAge x Sex
nested in levels of Test
Phillip Alday, Douglas Bates, and Reinhold Kliegl
+2023-08-21
+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 four main sections:
+Test
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 ProgressMeter
+using Random
+using Statistics
+using StatsBase
+using SMLP2023: dataset
+
+ProgressMeter.ijulia_behavior(:clear)
+CairoMakie.activate!(; type="svg")
Number of scores: 525126 in dataset(:fggk21)
Cohort: 9 levels; 2011-2019
School: 515 levels
Child: 108295 levels; all children are between 8.0 and 8.99 years old
Sex: “Girls” (n=55,086), “Boys” (n= 53,209)
age: testdate - middle of month of birthdate
Test: 5 levels
+Run
): 6 minute endurance run [m]; to nearest 9m in 9x18m fieldStar_r
): star coordination run [m/s]; 9x9m field, 4 x diagonal = 50.912 mS20_r
): 20-meters sprint [m/s]SLJ
): standing long jump [cm]BPT
): 1-kg medicine ball push test [m]score - see units
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)
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 | +
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.
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 | +
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 (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 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))
Row | +Test | +Cohort | +School | +Child | +Sex | +age | +a1 | +zScore | +
---|---|---|---|---|---|---|---|---|
+ | Cat… | +String | +String | +String | +Cat… | +Float64 | +Float64 | +Float64 | +
1 | +Speed | +2013 | +S101540 | +C002482 | +Girls | +7.99452 | +-0.505476 | +0.0551481 | +
2 | +PowerUP | +2013 | +S101540 | +C002482 | +Girls | +7.99452 | +-0.505476 | +-0.342524 | +
3 | +PowerLOW | +2013 | +S101540 | +C002482 | +Girls | +7.99452 | +-0.505476 | +0.74162 | +
4 | +Coordination | +2013 | +S101540 | +C002482 | +Girls | +7.99452 | +-0.505476 | +-0.209186 | +
5 | +Endurance | +2013 | +S101540 | +C002482 | +Girls | +7.99452 | +-0.505476 | +-0.127938 | +
6 | +Speed | +2013 | +S102090 | +C002513 | +Girls | +7.99452 | +-0.505476 | +-0.422921 | +
7 | +PowerUP | +2013 | +S102090 | +C002513 | +Girls | +7.99452 | +-0.505476 | +-0.762963 | +
8 | +PowerLOW | +2013 | +S102090 | +C002513 | +Girls | +7.99452 | +-0.505476 | +-0.188712 | +
9 | +Coordination | +2013 | +S102090 | +C002513 | +Girls | +7.99452 | +-0.505476 | +0.81382 | +
10 | +Endurance | +2013 | +S102090 | +C002513 | +Girls | +7.99452 | +-0.505476 | +0.452557 | +
11 | +Speed | +2013 | +S103299 | +C002621 | +Girls | +7.99452 | +-0.505476 | +0.859704 | +
12 | +PowerUP | +2013 | +S103299 | +C002621 | +Girls | +7.99452 | +-0.505476 | +0.218061 | +
13 | +PowerLOW | +2013 | +S103299 | +C002621 | +Girls | +7.99452 | +-0.505476 | +0.74162 | +
⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +
9652 | +Coordination | +2018 | +S111508 | +C117805 | +Girls | +9.10609 | +0.606092 | +-0.95589 | +
9653 | +Endurance | +2018 | +S111508 | +C117805 | +Girls | +9.10609 | +0.606092 | +-1.67367 | +
9654 | +Speed | +2018 | +S111648 | +C117827 | +Boys | +9.10609 | +0.606092 | +0.859704 | +
9655 | +PowerUP | +2018 | +S111648 | +C117827 | +Boys | +9.10609 | +0.606092 | +0.6385 | +
9656 | +PowerLOW | +2018 | +S111648 | +C117827 | +Boys | +9.10609 | +0.606092 | +0.483194 | +
9657 | +Coordination | +2018 | +S111648 | +C117827 | +Boys | +9.10609 | +0.606092 | +0.708474 | +
9658 | +Endurance | +2018 | +S111648 | +C117827 | +Boys | +9.10609 | +0.606092 | +-0.337186 | +
9659 | +Speed | +2018 | +S112197 | +C117868 | +Boys | +9.10609 | +0.606092 | +0.578748 | +
9660 | +PowerUP | +2018 | +S112197 | +C117868 | +Boys | +9.10609 | +0.606092 | +0.498354 | +
9661 | +PowerLOW | +2018 | +S112197 | +C117868 | +Boys | +9.10609 | +0.606092 | +-0.550508 | +
9662 | +Coordination | +2018 | +S112197 | +C117868 | +Boys | +9.10609 | +0.606092 | +0.538978 | +
9663 | +Endurance | +2018 | +S112197 | +C117868 | +Boys | +9.10609 | +0.606092 | +0.938553 | +
Due to missing scores for some tests we have about 2% less than 10,000 observtions.
+age x Sex x Test
interactionThe 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,
+)
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 | +
The core results of the article are reported in Figure 2 of Fühner et al. (2021). In summary:
+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.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
.Test
and Sex
: The sex difference is test dependent, that is the difference between the slopes depends on Test
.age x Sex x Test
interaction, that is the slopes for boys and girls are statistically parallel for each of the five tests.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),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
.
Row | +Sex | +Test | +age | +zScore | +n | +
---|---|---|---|---|---|
+ | Cat… | +Cat… | +Float64 | +Float64 | +Int64 | +
1 | +Girls | +Speed | +8.0 | +-0.323114 | +28 | +
2 | +Girls | +PowerUP | +8.0 | +-0.590476 | +26 | +
3 | +Girls | +PowerLOW | +8.0 | +0.0677992 | +27 | +
4 | +Girls | +Coordination | +8.0 | +0.0273318 | +25 | +
5 | +Girls | +Endurance | +8.0 | +-0.17337 | +26 | +
6 | +Boys | +Speed | +8.0 | +0.394634 | +19 | +
7 | +Boys | +PowerUP | +8.0 | +0.328703 | +19 | +
8 | +Boys | +PowerLOW | +8.0 | +0.105077 | +19 | +
9 | +Boys | +Coordination | +8.0 | +-0.170018 | +19 | +
10 | +Boys | +Endurance | +8.0 | +0.407084 | +19 | +
11 | +Girls | +Speed | +8.1 | +-0.205934 | +54 | +
12 | +Girls | +PowerUP | +8.1 | +-0.653961 | +54 | +
13 | +Girls | +PowerLOW | +8.1 | +-0.157127 | +54 | +
⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +⋮ | +
109 | +Girls | +Coordination | +9.0 | +-0.000188842 | +68 | +
110 | +Girls | +Endurance | +9.0 | +-0.234448 | +68 | +
111 | +Girls | +Speed | +9.1 | +-0.285047 | +25 | +
112 | +Girls | +PowerUP | +9.1 | +-0.112684 | +25 | +
113 | +Girls | +PowerLOW | +9.1 | +-0.117645 | +24 | +
114 | +Girls | +Coordination | +9.1 | +-0.127844 | +25 | +
115 | +Girls | +Endurance | +9.1 | +-0.587496 | +24 | +
116 | +Boys | +Speed | +9.1 | +0.379808 | +17 | +
117 | +Boys | +PowerUP | +9.1 | +0.523085 | +17 | +
118 | +Boys | +PowerLOW | +9.1 | +0.294696 | +17 | +
119 | +Boys | +Coordination | +9.1 | +0.209309 | +16 | +
120 | +Boys | +Endurance | +9.1 | +0.266512 | +16 | +
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),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.
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) | +
Row | +Sex | +Test | +coef | +
---|---|---|---|
+ | Cat… | +Cat… | +Tuple… | +
1 | +Boys | +Endurance | +(0.39203, -0.0150694) | +
2 | +Boys | +Coordination | +(-3.33518, 0.401051) | +
3 | +Boys | +Speed | +(-1.75685, 0.228662) | +
4 | +Boys | +PowerLOW | +(-1.06646, 0.151546) | +
5 | +Boys | +PowerUP | +(-4.15536, 0.5245) | +
6 | +Girls | +Endurance | +(0.941712, -0.141158) | +
7 | +Girls | +Coordination | +(-0.681898, 0.0714891) | +
8 | +Girls | +Speed | +(-0.786382, 0.0725931) | +
9 | +Girls | +PowerLOW | +(-0.208472, 0.00150731) | +
10 | +Girls | +PowerUP | +(-5.23593, 0.570806) | +
Test
SeqDiffCoding was used in the publication. This specification tests pairwise differences between the five neighboring levels of Test
, that is:
Star_r
- Run
(2-1)S20_r
- Star_r
(3-2)SLJ
- S20_r
(4-3)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.
+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.
+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.
m_ovi
: allowing only varying intercepts (“Grand Means”);m_zcp
: adding variance components (VCs) for the four Test
contrasts, Sex
, and age
to LMM m_ovi
, yielding the zero-correlation parameters 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.
+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)
+end
+ | Est. | +SE | +z | +p | +σ_School | +
---|---|---|---|---|---|
(Intercept) | +-0.0174 | +0.0198 | +-0.88 | +0.3787 | +0.3417 | +
Test: Star_r | +-0.0026 | +0.0303 | +-0.09 | +0.9304 | ++ |
Test: S20_r | +0.0097 | +0.0302 | +0.32 | +0.7469 | ++ |
Test: SLJ | +0.0033 | +0.0300 | +0.11 | +0.9119 | ++ |
Test: BMT | +-0.0539 | +0.0299 | +-1.80 | +0.0711 | ++ |
Sex: Girls | +-0.4372 | +0.0205 | +-21.35 | +<1e-99 | ++ |
a1 | +0.1761 | +0.0354 | +4.97 | +<1e-06 | ++ |
Test: Star_r & Sex: Girls | +0.3802 | +0.0605 | +6.28 | +<1e-09 | ++ |
Test: S20_r & Sex: Girls | +-0.2038 | +0.0604 | +-3.38 | +0.0007 | ++ |
Test: SLJ & Sex: Girls | +-0.0669 | +0.0600 | +-1.11 | +0.2651 | ++ |
Test: BMT & Sex: Girls | +-0.2730 | +0.0598 | +-4.57 | +<1e-05 | ++ |
Test: Star_r & a1 | +0.3146 | +0.1038 | +3.03 | +0.0024 | ++ |
Test: S20_r & a1 | +-0.0767 | +0.1034 | +-0.74 | +0.4584 | ++ |
Test: SLJ & a1 | +-0.0745 | +0.1027 | +-0.73 | +0.4683 | ++ |
Test: BMT & a1 | +0.4724 | +0.1025 | +4.61 | +<1e-05 | ++ |
Sex: Girls & a1 | +-0.1044 | +0.0709 | +-1.47 | +0.1405 | ++ |
Test: Star_r & Sex: Girls & a1 | +-0.1912 | +0.2076 | +-0.92 | +0.3570 | ++ |
Test: S20_r & Sex: Girls & a1 | +0.1742 | +0.2068 | +0.84 | +0.3997 | ++ |
Test: SLJ & Sex: Girls & a1 | +0.0084 | +0.2054 | +0.04 | +0.9672 | ++ |
Test: BMT & Sex: Girls & a1 | +0.1923 | +0.2050 | +0.94 | +0.3482 | ++ |
Residual | +0.9152 | ++ | + | + | + |
Is the model singular (overparameterized, degenerate)? In other words: Is the model not supported by the data?
+ +Models varying only in intercepts are almost always supported by the data.
+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)
+end
Minimizing 148 Time: 0:00:00 ( 1.87 ms/it)
+ objective: 25923.460469750346
++ | Est. | +SE | +z | +p | +σ_School | +
---|---|---|---|---|---|
(Intercept) | +-0.0247 | +0.0198 | +-1.24 | +0.2136 | +0.3257 | +
Test: Star_r | +0.0029 | +0.0311 | +0.09 | +0.9254 | +0.2208 | +
Test: S20_r | +0.0102 | +0.0286 | +0.36 | +0.7209 | +0.0503 | +
Test: SLJ | +0.0049 | +0.0283 | +0.17 | +0.8634 | +0.0000 | +
Test: BMT | +-0.0460 | +0.0309 | +-1.49 | +0.1368 | +0.2286 | +
Sex: Girls | +-0.4399 | +0.0322 | +-13.64 | +<1e-41 | +0.4503 | +
a1 | +0.1840 | +0.0573 | +3.21 | +0.0013 | +0.7780 | +
Test: Star_r & Sex: Girls | +0.3784 | +0.0577 | +6.55 | +<1e-10 | ++ |
Test: S20_r & Sex: Girls | +-0.2019 | +0.0570 | +-3.54 | +0.0004 | ++ |
Test: SLJ & Sex: Girls | +-0.0666 | +0.0566 | +-1.18 | +0.2392 | ++ |
Test: BMT & Sex: Girls | +-0.2780 | +0.0570 | +-4.87 | +<1e-05 | ++ |
Test: Star_r & a1 | +0.3101 | +0.0992 | +3.13 | +0.0018 | ++ |
Test: S20_r & a1 | +-0.0728 | +0.0976 | +-0.75 | +0.4558 | ++ |
Test: SLJ & a1 | +-0.0784 | +0.0968 | +-0.81 | +0.4180 | ++ |
Test: BMT & a1 | +0.4741 | +0.0979 | +4.84 | +<1e-05 | ++ |
Sex: Girls & a1 | +-0.1453 | +0.0791 | +-1.84 | +0.0662 | ++ |
Test: Star_r & Sex: Girls & a1 | +-0.1570 | +0.1982 | +-0.79 | +0.4284 | ++ |
Test: S20_r & Sex: Girls & a1 | +0.1799 | +0.1952 | +0.92 | +0.3566 | ++ |
Test: SLJ & Sex: Girls & a1 | +0.0061 | +0.1936 | +0.03 | +0.9749 | ++ |
Test: BMT & Sex: Girls & a1 | +0.1657 | +0.1958 | +0.85 | +0.3974 | ++ |
Residual | +0.8623 | ++ | + | + | + |
Depending on sampling, this model estimating variance components for School
may or may not be supported by the data.
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)
+end
Minimizing 2515 Time: 0:00:02 ( 0.85 ms/it)
+ objective: 25864.389680661036
++ | Est. | +SE | +z | +p | +σ_School | +
---|---|---|---|---|---|
(Intercept) | +-0.0268 | +0.0198 | +-1.35 | +0.1757 | +0.3264 | +
Test: Star_r | +0.0056 | +0.0310 | +0.18 | +0.8575 | +0.2203 | +
Test: S20_r | +0.0104 | +0.0301 | +0.35 | +0.7294 | +0.1757 | +
Test: SLJ | +0.0024 | +0.0294 | +0.08 | +0.9357 | +0.1477 | +
Test: BMT | +-0.0401 | +0.0305 | +-1.32 | +0.1883 | +0.2245 | +
Sex: Girls | +-0.4356 | +0.0320 | +-13.61 | +<1e-41 | +0.4482 | +
a1 | +0.1877 | +0.0566 | +3.31 | +0.0009 | +0.7688 | +
Test: Star_r & Sex: Girls | +0.3759 | +0.0575 | +6.53 | +<1e-10 | ++ |
Test: S20_r & Sex: Girls | +-0.1969 | +0.0572 | +-3.44 | +0.0006 | ++ |
Test: SLJ & Sex: Girls | +-0.0673 | +0.0567 | +-1.19 | +0.2350 | ++ |
Test: BMT & Sex: Girls | +-0.2717 | +0.0566 | +-4.80 | +<1e-05 | ++ |
Test: Star_r & a1 | +0.3112 | +0.0988 | +3.15 | +0.0016 | ++ |
Test: S20_r & a1 | +-0.0707 | +0.0981 | +-0.72 | +0.4711 | ++ |
Test: SLJ & a1 | +-0.0642 | +0.0971 | +-0.66 | +0.5085 | ++ |
Test: BMT & a1 | +0.4678 | +0.0971 | +4.82 | +<1e-05 | ++ |
Sex: Girls & a1 | +-0.1405 | +0.0783 | +-1.79 | +0.0728 | ++ |
Test: Star_r & Sex: Girls & a1 | +-0.1632 | +0.1974 | +-0.83 | +0.4085 | ++ |
Test: S20_r & Sex: Girls & a1 | +0.1862 | +0.1961 | +0.95 | +0.3424 | ++ |
Test: SLJ & Sex: Girls & a1 | +0.0011 | +0.1941 | +0.01 | +0.9954 | ++ |
Test: BMT & Sex: Girls & a1 | +0.1638 | +0.1942 | +0.84 | +0.3989 | ++ |
Residual | +0.8598 | ++ | + | + | + |
We also need to see the VCs and CPs of the random-effect structure (RES).
++ | Column | +Variance | +Std.Dev | +Corr. | ++ | + | + | + | + |
---|---|---|---|---|---|---|---|---|---|
School | +(Intercept) | +0.106543 | +0.326409 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.048516 | +0.220263 | ++0.10 | ++ | + | + | + | + |
+ | Test: S20_r | +0.030859 | +0.175666 | +-0.19 | +-0.22 | ++ | + | + | + |
+ | Test: SLJ | +0.021828 | +0.147744 | +-0.10 | ++0.07 | +-0.83 | ++ | + | + |
+ | Test: BMT | +0.050409 | +0.224519 | +-0.75 | ++0.43 | +-0.13 | ++0.16 | ++ | + |
+ | Sex: Girls | +0.200924 | +0.448245 | +-0.04 | ++0.25 | ++0.20 | +-0.26 | +-0.10 | ++ |
+ | a1 | +0.591016 | +0.768776 | ++0.02 | +-0.13 | ++0.03 | +-0.49 | ++0.17 | +-0.06 | +
Residual | ++ | 0.739266 | +0.859806 | ++ | + | + | + | + | + |
The complex model may or may not be supported by the data.
+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?
+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.
++ | model-dof | +deviance | +χ² | +χ²-dof | +P(>χ²) | +
---|---|---|---|---|---|
zScore ~ 1 + Test + Sex + a1 + Test & Sex + Test & a1 + Sex & a1 + Test & Sex & a1 + (1 | School) | +22 | +26273 | ++ | + | + |
zScore ~ 1 + Test + Sex + a1 + Test & Sex + Test & a1 + Sex & a1 + Test & Sex & a1 + zerocorr(1 + Test + Sex + a1 | School) | +28 | +25923 | +349 | +6 | +<1e-71 | +
zScore ~ 1 + Test + Sex + a1 + Test & Sex + Test & a1 + Sex & a1 + Test & Sex & a1 + (1 + Test + Sex + a1 | School) | +49 | +25864 | +59 | +21 | +<1e-04 | +
Row | +name | +dof | +deviance | +AIC | +AICc | +BIC | +
---|---|---|---|---|---|---|
+ | Symbol | +Int64 | +Float64 | +Float64 | +Float64 | +Float64 | +
1 | +m_ovi | +22 | +26272.9 | +26316.9 | +26317.0 | +26474.8 | +
2 | +m_zcp | +28 | +25923.5 | +25979.5 | +25979.6 | +26180.4 | +
3 | +m_cpx | +49 | +25864.4 | +25962.4 | +25962.9 | +26314.0 | +
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. Basicly, 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.
+m_ovi
, m_zcp
, and m_cpx
We check whether enriching the RES changed the significance of fixed effects in the final model.
+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)
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.02 | +-0.03 | +-0.88 | +-1.24 | +-1.35 | +
2 | +Test: Star_r | +-0.0 | +0.0 | +0.01 | +-0.09 | +0.09 | +0.18 | +
3 | +Test: S20_r | +0.01 | +0.01 | +0.01 | +0.32 | +0.36 | +0.35 | +
4 | +Test: SLJ | +0.0 | +0.0 | +0.0 | +0.11 | +0.17 | +0.08 | +
5 | +Test: BMT | +-0.05 | +-0.05 | +-0.04 | +-1.8 | +-1.49 | +-1.32 | +
6 | +Sex: Girls | +-0.44 | +-0.44 | +-0.44 | +-21.35 | +-13.64 | +-13.61 | +
7 | +a1 | +0.18 | +0.18 | +0.19 | +4.97 | +3.21 | +3.31 | +
8 | +Test: Star_r & Sex: Girls | +0.38 | +0.38 | +0.38 | +6.28 | +6.55 | +6.53 | +
9 | +Test: S20_r & Sex: Girls | +-0.2 | +-0.2 | +-0.2 | +-3.38 | +-3.54 | +-3.44 | +
10 | +Test: SLJ & Sex: Girls | +-0.07 | +-0.07 | +-0.07 | +-1.11 | +-1.18 | +-1.19 | +
11 | +Test: BMT & Sex: Girls | +-0.27 | +-0.28 | +-0.27 | +-4.57 | +-4.87 | +-4.8 | +
12 | +Test: Star_r & a1 | +0.31 | +0.31 | +0.31 | +3.03 | +3.13 | +3.15 | +
13 | +Test: S20_r & a1 | +-0.08 | +-0.07 | +-0.07 | +-0.74 | +-0.75 | +-0.72 | +
14 | +Test: SLJ & a1 | +-0.07 | +-0.08 | +-0.06 | +-0.73 | +-0.81 | +-0.66 | +
15 | +Test: BMT & a1 | +0.47 | +0.47 | +0.47 | +4.61 | +4.84 | +4.82 | +
16 | +Sex: Girls & a1 | +-0.1 | +-0.15 | +-0.14 | +-1.47 | +-1.84 | +-1.79 | +
17 | +Test: Star_r & Sex: Girls & a1 | +-0.19 | +-0.16 | +-0.16 | +-0.92 | +-0.79 | +-0.83 | +
18 | +Test: S20_r & Sex: Girls & a1 | +0.17 | +0.18 | +0.19 | +0.84 | +0.92 | +0.95 | +
19 | +Test: SLJ & Sex: Girls & a1 | +0.01 | +0.01 | +0.0 | +0.04 | +0.03 | +0.01 | +
20 | +Test: BMT & Sex: Girls & a1 | +0.19 | +0.17 | +0.16 | +0.94 | +0.85 | +0.84 | +
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.
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)
+end
+ | Est. | +SE | +z | +p | +σ_Cohort | +
---|---|---|---|---|---|
(Intercept) | +-0.0009 | +0.0161 | +-0.06 | +0.9548 | +0.0378 | +
Test: Star_r | +-0.0073 | +0.0380 | +-0.19 | +0.8482 | +0.0614 | +
Test: S20_r | +0.0056 | +0.0383 | +0.15 | +0.8845 | +0.0636 | +
Test: SLJ | +0.0101 | +0.0451 | +0.22 | +0.8226 | +0.0959 | +
Test: BMT | +-0.0556 | +0.0335 | +-1.66 | +0.0968 | +0.0330 | +
a1 | +0.2055 | +0.0349 | +5.90 | +<1e-08 | ++ |
Sex: Girls | +-0.4238 | +0.0201 | +-21.07 | +<1e-97 | ++ |
Test: Star_r & a1 | +0.2849 | +0.1101 | +2.59 | +0.0097 | ++ |
Test: S20_r & a1 | +-0.1172 | +0.1096 | +-1.07 | +0.2851 | ++ |
Test: SLJ & a1 | +-0.0270 | +0.1094 | +-0.25 | +0.8049 | ++ |
Test: BMT & a1 | +0.4555 | +0.1084 | +4.20 | +<1e-04 | ++ |
Test: Star_r & Sex: Girls | +0.3700 | +0.0640 | +5.78 | +<1e-08 | ++ |
Test: S20_r & Sex: Girls | +-0.2116 | +0.0638 | +-3.32 | +0.0009 | ++ |
Test: SLJ & Sex: Girls | +-0.0552 | +0.0634 | +-0.87 | +0.3844 | ++ |
Test: BMT & Sex: Girls | +-0.2718 | +0.0632 | +-4.30 | +<1e-04 | ++ |
a1 & Sex: Girls | +-0.1368 | +0.0690 | +-1.98 | +0.0473 | ++ |
Test: Star_r & a1 & Sex: Girls | +-0.2099 | +0.2194 | +-0.96 | +0.3387 | ++ |
Test: S20_r & a1 & Sex: Girls | +0.1609 | +0.2186 | +0.74 | +0.4616 | ++ |
Test: SLJ & a1 & Sex: Girls | +0.0225 | +0.2172 | +0.10 | +0.9174 | ++ |
Test: BMT & a1 & Sex: Girls | +0.1901 | +0.2166 | +0.88 | +0.3801 | ++ |
Residual | +0.9671 | ++ | + | + | + |
+ | Column | +Variance | +Std.Dev | +Corr. | ++ | + | + |
---|---|---|---|---|---|---|---|
Cohort | +(Intercept) | +0.0014273 | +0.0377794 | ++ | + | + | + |
+ | Test: Star_r | +0.0037651 | +0.0613601 | +-0.90 | ++ | + | + |
+ | Test: S20_r | +0.0040449 | +0.0635995 | +-1.00 | ++0.87 | ++ | + |
+ | Test: SLJ | +0.0091885 | +0.0958567 | ++0.98 | +-0.97 | +-0.97 | ++ |
+ | Test: BMT | +0.0010882 | +0.0329872 | +-1.00 | ++0.91 | ++0.99 | +-0.99 | +
Residual | ++ | 0.9353139 | +0.9671163 | ++ | + | + | + |
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)
+end
+ | Est. | +SE | +z | +p | +σ_Cohort | +
---|---|---|---|---|---|
(Intercept) | +-0.0022 | +0.0152 | +-0.15 | +0.8837 | +0.0341 | +
Test: Star_r | +-0.0042 | +0.0339 | +-0.12 | +0.9023 | +0.0331 | +
Test: S20_r | +0.0088 | +0.0319 | +0.27 | +0.7837 | +0.0000 | +
Test: SLJ | +0.0045 | +0.0317 | +0.14 | +0.8876 | +0.0000 | +
Test: BMT | +-0.0536 | +0.0316 | +-1.69 | +0.0903 | +0.0000 | +
a1 | +0.1999 | +0.0351 | +5.70 | +<1e-07 | ++ |
Sex: Girls | +-0.4245 | +0.0201 | +-21.08 | +<1e-97 | ++ |
Test: Star_r & a1 | +0.3078 | +0.1101 | +2.80 | +0.0052 | ++ |
Test: S20_r & a1 | +-0.0849 | +0.1093 | +-0.78 | +0.4377 | ++ |
Test: SLJ & a1 | +-0.0748 | +0.1086 | +-0.69 | +0.4911 | ++ |
Test: BMT & a1 | +0.4717 | +0.1084 | +4.35 | +<1e-04 | ++ |
Test: Star_r & Sex: Girls | +0.3729 | +0.0640 | +5.82 | +<1e-08 | ++ |
Test: S20_r & Sex: Girls | +-0.2079 | +0.0638 | +-3.26 | +0.0011 | ++ |
Test: SLJ & Sex: Girls | +-0.0611 | +0.0634 | +-0.96 | +0.3354 | ++ |
Test: BMT & Sex: Girls | +-0.2697 | +0.0632 | +-4.27 | +<1e-04 | ++ |
a1 & Sex: Girls | +-0.1372 | +0.0691 | +-1.99 | +0.0470 | ++ |
Test: Star_r & a1 & Sex: Girls | +-0.2041 | +0.2195 | +-0.93 | +0.3525 | ++ |
Test: S20_r & a1 & Sex: Girls | +0.1733 | +0.2187 | +0.79 | +0.4280 | ++ |
Test: SLJ & a1 & Sex: Girls | +0.0051 | +0.2172 | +0.02 | +0.9811 | ++ |
Test: BMT & a1 & Sex: Girls | +0.1962 | +0.2168 | +0.90 | +0.3655 | ++ |
Residual | +0.9680 | ++ | + | + | + |
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)
+end
+ | Est. | +SE | +z | +p | +σ_Cohort | +
---|---|---|---|---|---|
(Intercept) | +-0.0022 | +0.0152 | +-0.15 | +0.8835 | +0.0340 | +
Test: Star_r | +-0.0032 | +0.0320 | +-0.10 | +0.9196 | ++ |
Test: S20_r | +0.0087 | +0.0319 | +0.27 | +0.7847 | ++ |
Test: SLJ | +0.0045 | +0.0317 | +0.14 | +0.8869 | ++ |
Test: BMT | +-0.0536 | +0.0316 | +-1.69 | +0.0903 | ++ |
a1 | +0.1999 | +0.0351 | +5.69 | +<1e-07 | ++ |
Sex: Girls | +-0.4245 | +0.0201 | +-21.08 | +<1e-97 | ++ |
Test: Star_r & a1 | +0.3135 | +0.1097 | +2.86 | +0.0043 | ++ |
Test: S20_r & a1 | +-0.0848 | +0.1093 | +-0.78 | +0.4378 | ++ |
Test: SLJ & a1 | +-0.0748 | +0.1086 | +-0.69 | +0.4910 | ++ |
Test: BMT & a1 | +0.4718 | +0.1084 | +4.35 | +<1e-04 | ++ |
Test: Star_r & Sex: Girls | +0.3736 | +0.0640 | +5.84 | +<1e-08 | ++ |
Test: S20_r & Sex: Girls | +-0.2079 | +0.0638 | +-3.26 | +0.0011 | ++ |
Test: SLJ & Sex: Girls | +-0.0611 | +0.0634 | +-0.96 | +0.3356 | ++ |
Test: BMT & Sex: Girls | +-0.2697 | +0.0632 | +-4.27 | +<1e-04 | ++ |
a1 & Sex: Girls | +-0.1371 | +0.0691 | +-1.99 | +0.0471 | ++ |
Test: Star_r & a1 & Sex: Girls | +-0.2022 | +0.2195 | +-0.92 | +0.3568 | ++ |
Test: S20_r & a1 & Sex: Girls | +0.1733 | +0.2187 | +0.79 | +0.4281 | ++ |
Test: SLJ & a1 & Sex: Girls | +0.0051 | +0.2172 | +0.02 | +0.9811 | ++ |
Test: BMT & a1 & Sex: Girls | +0.1961 | +0.2168 | +0.90 | +0.3657 | ++ |
Residual | +0.9681 | ++ | + | + | + |
This solves the problem with singularity, but does LMM m_zcpCohort
fit noise relative to the LMM m_oviCohort
?
+ | model-dof | +deviance | +χ² | +χ²-dof | +P(>χ²) | +
---|---|---|---|---|---|
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + (1 | Cohort) | +22 | +26803 | ++ | + | + |
zScore ~ 1 + Test + a1 + Sex + Test & a1 + Test & Sex + a1 & Sex + Test & a1 & Sex + zerocorr(1 + Test | Cohort) | +26 | +26803 | +0 | +4 | +0.9968 | +
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
Row | +dof | +deviance | +AIC | +AICc | +BIC | +
---|---|---|---|---|---|
+ | Int64 | +Float64 | +Float64 | +Float64 | +Float64 | +
1 | +22 | +26802.9 | +26846.9 | +26847.0 | +27004.7 | +
2 | +26 | +26802.7 | +26854.7 | +26854.8 | +27041.3 | +
3 | +36 | +26790.5 | +26862.5 | +26862.7 | +27120.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
?
m1
to the reduced dataThe 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 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 m1
are [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)
+end
Minimizing 2162 Time: 0:00:31 (14.80 ms/it)
+ objective: 24651.01322999849
++ | Est. | +SE | +z | +p | +σ_Child | +σ_School | +σ_Cohort | +
---|---|---|---|---|---|---|---|
(Intercept) | +-0.0133 | +0.0192 | +-0.69 | +0.4878 | +0.5908 | +0.2138 | +0.0132 | +
Test: Star_r | +0.0088 | +0.0375 | +0.24 | +0.8135 | +0.7705 | +0.3329 | +0.0637 | +
Test: S20_r | +0.0092 | +0.0293 | +0.31 | +0.7542 | +0.6656 | +0.3281 | +0.0068 | +
Test: SLJ | +0.0034 | +0.0299 | +0.11 | +0.9105 | +0.5786 | +0.3212 | +0.0339 | +
Test: BMT | +-0.0383 | +0.0296 | +-1.29 | +0.1954 | +0.7491 | +0.3130 | +0.0000 | +
a1 | +0.1955 | +0.0544 | +3.59 | +0.0003 | ++ | 0.2823 | ++ |
Sex: Girls | +-0.4287 | +0.0312 | +-13.75 | +<1e-42 | ++ | 0.1438 | ++ |
Test: Star_r & a1 | +0.2898 | +0.0887 | +3.27 | +0.0011 | ++ | + | + |
Test: S20_r & a1 | +-0.0784 | +0.0817 | +-0.96 | +0.3373 | ++ | + | + |
Test: SLJ & a1 | +-0.0631 | +0.0769 | +-0.82 | +0.4116 | ++ | + | + |
Test: BMT & a1 | +0.4713 | +0.0851 | +5.54 | +<1e-07 | ++ | + | + |
Test: Star_r & Sex: Girls | +0.3756 | +0.0511 | +7.35 | +<1e-12 | ++ | + | + |
Test: S20_r & Sex: Girls | +-0.1867 | +0.0475 | +-3.93 | +<1e-04 | ++ | + | + |
Test: SLJ & Sex: Girls | +-0.0695 | +0.0445 | +-1.56 | +0.1184 | ++ | + | + |
Test: BMT & Sex: Girls | +-0.2812 | +0.0495 | +-5.68 | +<1e-07 | ++ | + | + |
a1 & Sex: Girls | +-0.1267 | +0.1048 | +-1.21 | +0.2269 | ++ | + | + |
Test: Star_r & a1 & Sex: Girls | +-0.1382 | +0.1756 | +-0.79 | +0.4312 | ++ | + | + |
Test: S20_r & a1 & Sex: Girls | +0.1753 | +0.1633 | +1.07 | +0.2832 | ++ | + | + |
Test: SLJ & a1 & Sex: Girls | +-0.0154 | +0.1530 | +-0.10 | +0.9196 | ++ | + | + |
Test: BMT & a1 & Sex: Girls | +0.1548 | +0.1702 | +0.91 | +0.3630 | ++ | + | + |
Residual | +0.5136 | ++ | + | + | + | + | + |
+ | Column | +Variance | +Std.Dev | +Corr. | ++ | + | + | + | + |
---|---|---|---|---|---|---|---|---|---|
Child | +(Intercept) | +0.34903980 | +0.59079590 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.59364220 | +0.77048180 | ++0.14 | ++ | + | + | + | + |
+ | Test: S20_r | +0.44304305 | +0.66561479 | ++0.00 | +-0.52 | ++ | + | + | + |
+ | Test: SLJ | +0.33481089 | +0.57862846 | ++0.05 | +-0.04 | +-0.37 | ++ | + | + |
+ | Test: BMT | +0.56120773 | +0.74913799 | +-0.33 | ++0.13 | +-0.17 | +-0.25 | ++ | + |
School | +(Intercept) | +0.04569734 | +0.21376937 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.11081338 | +0.33288644 | +-0.06 | ++ | + | + | + | + |
+ | Test: S20_r | +0.10763442 | +0.32807686 | +-0.08 | +-0.39 | ++ | + | + | + |
+ | Test: SLJ | +0.10314745 | +0.32116576 | +-0.19 | ++0.21 | +-0.80 | ++ | + | + |
+ | Test: BMT | +0.09796213 | +0.31298902 | +-0.33 | +-0.02 | ++0.13 | +-0.38 | ++ | + |
+ | a1 | +0.07971464 | +0.28233780 | ++0.62 | +-0.20 | ++0.11 | +-0.61 | ++0.45 | ++ |
+ | Sex: Girls | +0.02066903 | +0.14376729 | +-0.45 | ++0.45 | ++0.31 | +-0.27 | +-0.15 | +-0.37 | +
Cohort | +(Intercept) | +0.00017336 | +0.01316656 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.00405975 | +0.06371614 | +. | ++ | + | + | + | + |
+ | Test: S20_r | +0.00004631 | +0.00680502 | +. | +. | ++ | + | + | + |
+ | Test: SLJ | +0.00114631 | +0.03385719 | +. | +. | +. | ++ | + | + |
+ | Test: BMT | +0.00000000 | +0.00000000 | +. | +. | +. | +. | ++ | + |
Residual | ++ | 0.26377271 | +0.51358807 | ++ | + | + | + | + | + |
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)
+end
Minimizing 1464 Time: 0:00:21 (14.46 ms/it)
+ objective: 24646.87190533457
++ | Est. | +SE | +z | +p | +σ_Child | +σ_School | +σ_Cohort | +
---|---|---|---|---|---|---|---|
(Intercept) | +-0.0136 | +0.0195 | +-0.70 | +0.4844 | ++ | + | + |
Test: Star_r | +0.0093 | +0.0365 | +0.26 | +0.7987 | +0.7972 | +0.3103 | +0.0208 | +
Test: S20_r | +0.0060 | +0.0354 | +0.17 | +0.8653 | +0.7678 | +0.3305 | +0.0563 | +
Test: SLJ | +0.0052 | +0.0339 | +0.15 | +0.8787 | +0.7806 | +0.2478 | +0.0157 | +
Test: BMT | +-0.0388 | +0.0300 | +-1.29 | +0.1967 | +0.6970 | +0.1969 | +0.0000 | +
a1 | +0.1935 | +0.0545 | +3.55 | +0.0004 | ++ | 0.2833 | ++ |
Sex: Girls | +-0.4290 | +0.0312 | +-13.75 | +<1e-42 | ++ | 0.1433 | ++ |
Test: Star_r & a1 | +0.2932 | +0.0887 | +3.31 | +0.0009 | ++ | + | + |
Test: S20_r & a1 | +-0.1007 | +0.0826 | +-1.22 | +0.2228 | ++ | + | + |
Test: SLJ & a1 | +-0.0456 | +0.0773 | +-0.59 | +0.5549 | ++ | + | + |
Test: BMT & a1 | +0.4685 | +0.0853 | +5.50 | +<1e-07 | ++ | + | + |
Test: Star_r & Sex: Girls | +0.3759 | +0.0511 | +7.36 | +<1e-12 | ++ | + | + |
Test: S20_r & Sex: Girls | +-0.1892 | +0.0475 | +-3.98 | +<1e-04 | ++ | + | + |
Test: SLJ & Sex: Girls | +-0.0680 | +0.0445 | +-1.53 | +0.1262 | ++ | + | + |
Test: BMT & Sex: Girls | +-0.2815 | +0.0495 | +-5.68 | +<1e-07 | ++ | + | + |
a1 & Sex: Girls | +-0.1262 | +0.1049 | +-1.20 | +0.2290 | ++ | + | + |
Test: Star_r & a1 & Sex: Girls | +-0.1371 | +0.1756 | +-0.78 | +0.4348 | ++ | + | + |
Test: S20_r & a1 & Sex: Girls | +0.1720 | +0.1634 | +1.05 | +0.2926 | ++ | + | + |
Test: SLJ & a1 & Sex: Girls | +-0.0125 | +0.1529 | +-0.08 | +0.9348 | ++ | + | + |
Test: BMT & a1 & Sex: Girls | +0.1542 | +0.1702 | +0.91 | +0.3648 | ++ | + | + |
Test: Run | ++ | + | + | + | 0.7501 | +0.3702 | +0.0548 | +
Residual | +0.5118 | ++ | + | + | + | + | + |
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).
+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.
+ | Column | +Variance | +Std.Dev | +Corr. | ++ | + | + | + | + |
---|---|---|---|---|---|---|---|---|---|
Child | +Test: Run | +0.5626315 | +0.7500877 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.6354871 | +0.7971744 | ++0.50 | ++ | + | + | + | + |
+ | Test: S20_r | +0.5894630 | +0.7677650 | ++0.56 | ++0.64 | ++ | + | + | + |
+ | Test: SLJ | +0.6093249 | +0.7805927 | ++0.55 | ++0.60 | ++0.72 | ++ | + | + |
+ | Test: BMT | +0.4858396 | +0.6970219 | ++0.19 | ++0.40 | ++0.37 | ++0.49 | ++ | + |
School | +Test: Run | +0.1370367 | +0.3701847 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.0962842 | +0.3102970 | ++0.53 | ++ | + | + | + | + |
+ | Test: S20_r | +0.1092459 | +0.3305237 | ++0.47 | ++0.48 | ++ | + | + | + |
+ | Test: SLJ | +0.0613916 | +0.2477733 | ++0.47 | ++0.76 | ++0.41 | ++ | + | + |
+ | Test: BMT | +0.0387555 | +0.1968641 | ++0.15 | ++0.38 | ++0.18 | ++0.02 | ++ | + |
+ | a1 | +0.0802430 | +0.2832720 | ++0.58 | ++0.47 | ++0.56 | +-0.05 | ++0.65 | ++ |
+ | Sex: Girls | +0.0205488 | +0.1433486 | +-0.63 | +-0.27 | ++0.06 | +-0.28 | +-0.58 | +-0.37 | +
Cohort | +Test: Run | +0.0030067 | +0.0548338 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.0004347 | +0.0208494 | +. | ++ | + | + | + | + |
+ | Test: S20_r | +0.0031725 | +0.0563249 | +. | +. | ++ | + | + | + |
+ | Test: SLJ | +0.0002452 | +0.0156604 | +. | +. | +. | ++ | + | + |
+ | Test: BMT | +0.0000000 | +0.0000000 | +. | +. | +. | +. | ++ | + |
Residual | ++ | 0.2619146 | +0.5117760 | ++ | + | + | + | + | + |
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 Fühner et al. (2021) for an application.
+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!
(Child =
+Principal components based on correlation matrix
+ (Intercept) 1.0 . . . .
+ Test: Star_r 0.14 1.0 . . .
+ Test: S20_r 0.0 -0.52 1.0 . .
+ Test: SLJ 0.05 -0.04 -0.37 1.0 .
+ Test: BMT -0.33 0.13 -0.17 -0.25 1.0
+
+Normalized cumulative variances:
+[0.3283, 0.6158, 0.8282, 0.9389, 1.0]
+
+Component loadings
+ PC1 PC2 PC3 PC4 PC5
+ (Intercept) -0.08 0.54 0.59 -0.59 -0.02
+ Test: Star_r -0.6 -0.11 0.45 0.42 0.5
+ Test: S20_r 0.71 0.04 0.15 0.06 0.68
+ Test: SLJ -0.33 0.44 -0.65 -0.22 0.48
+ Test: BMT -0.14 -0.71 -0.01 -0.65 0.24, School =
+Principal components based on correlation matrix
+ (Intercept) 1.0 . . . . . .
+ Test: Star_r -0.06 1.0 . . . . .
+ Test: S20_r -0.08 -0.39 1.0 . . . .
+ Test: SLJ -0.19 0.21 -0.8 1.0 . . .
+ Test: BMT -0.33 -0.02 0.13 -0.38 1.0 . .
+ a1 0.62 -0.2 0.11 -0.61 0.45 1.0 .
+ Sex: Girls -0.45 0.45 0.31 -0.27 -0.15 -0.37 1.0
+
+Normalized cumulative variances:
+[0.3569, 0.6351, 0.8058, 0.9687, 1.0, 1.0, 1.0]
+
+Component loadings
+ PC1 PC2 PC3 PC4 PC5 PC6 PC7
+ (Intercept) -0.25 -0.49 -0.42 -0.37 0.18 0.52 0.29
+ Test: Star_r 0.31 0.12 0.17 -0.75 0.47 -0.25 -0.1
+ Test: S20_r -0.4 0.44 -0.31 0.2 0.52 -0.29 0.39
+ Test: SLJ 0.55 -0.31 0.11 0.2 0.03 -0.21 0.71
+ Test: BMT -0.29 0.14 0.79 0.0 0.15 0.41 0.28
+ a1 -0.52 -0.27 0.16 -0.31 -0.41 -0.57 0.2
+ Sex: Girls 0.14 0.6 -0.2 -0.34 -0.53 0.23 0.36, 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 0.0
+
+Normalized cumulative variances:
+[0.25, 0.5, 0.75, 1.0, 1.0]
+
+Component loadings
+ PC1 PC2 PC3 PC4 PC5
+ (Intercept) 0.0 0.0 1.0 0.0 0.0
+ Test: Star_r 1.0 0.0 0.0 0.0 0.0
+ Test: S20_r 0.0 0.0 0.0 1.0 0.0
+ Test: SLJ 0.0 1.0 0.0 0.0 0.0
+ Test: BMT 0.0 0.0 0.0 0.0 NaN)
+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.
For the complete data, all PCs had unique variance associated with them.
+(Child =
+Principal components based on correlation matrix
+ Test: Run 1.0 . . . .
+ Test: Star_r 0.5 1.0 . . .
+ Test: S20_r 0.56 0.64 1.0 . .
+ Test: SLJ 0.55 0.6 0.72 1.0 .
+ Test: BMT 0.19 0.4 0.37 0.49 1.0
+
+Normalized cumulative variances:
+[0.6114, 0.7767, 0.8675, 0.9484, 1.0]
+
+Component loadings
+ PC1 PC2 PC3 PC4 PC5
+ Test: Run -0.42 0.53 0.63 0.37 0.08
+ Test: Star_r -0.47 0.03 -0.65 0.58 -0.16
+ Test: S20_r -0.49 0.15 -0.24 -0.49 0.66
+ Test: SLJ -0.5 -0.05 0.09 -0.5 -0.7
+ Test: BMT -0.34 -0.83 0.33 0.2 0.2, School =
+Principal components based on correlation matrix
+ Test: Run 1.0 . . . . . .
+ Test: Star_r 0.53 1.0 . . . . .
+ Test: S20_r 0.47 0.48 1.0 . . . .
+ Test: SLJ 0.47 0.76 0.41 1.0 . . .
+ Test: BMT 0.15 0.38 0.18 0.02 1.0 . .
+ a1 0.58 0.47 0.56 -0.05 0.65 1.0 .
+ Sex: Girls -0.63 -0.27 0.06 -0.28 -0.58 -0.37 1.0
+
+Normalized cumulative variances:
+[0.4822, 0.6938, 0.8508, 0.9553, 1.0, 1.0, 1.0]
+
+Component loadings
+ PC1 PC2 PC3 PC4 PC5 PC6 PC7
+ Test: Run -0.44 0.08 0.15 -0.64 -0.16 -0.51 0.29
+ Test: Star_r -0.44 0.28 0.02 0.41 -0.56 0.32 0.38
+ Test: S20_r -0.34 0.27 -0.58 -0.1 0.6 0.24 0.21
+ Test: SLJ -0.32 0.56 0.33 0.25 0.22 -0.23 -0.55
+ Test: BMT -0.32 -0.53 -0.02 0.55 0.25 -0.47 0.19
+ a1 -0.41 -0.35 -0.42 -0.14 -0.33 0.11 -0.62
+ Sex: Girls 0.34 0.35 -0.59 0.18 -0.29 -0.55 -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 0.0
+
+Normalized cumulative variances:
+[0.25, 0.5, 0.75, 1.0, 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 1.0 0.0 0.0
+ Test: SLJ 0.0 0.0 0.0 1.0 0.0
+ Test: BMT 0.0 0.0 0.0 0.0 NaN)
+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.
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)
+end
Minimizing 1781 Time: 0:00:25 (14.48 ms/it)
+ objective: 24646.871904034535
++ | Est. | +SE | +z | +p | +σ_Child | +σ_School | +σ_Cohort | +
---|---|---|---|---|---|---|---|
(Intercept) | +-0.0136 | +0.0195 | +-0.70 | +0.4844 | ++ | + | + |
Test: Star_r | +0.0093 | +0.0365 | +0.26 | +0.7987 | +0.7918 | +0.3103 | +0.0208 | +
Test: S20_r | +0.0060 | +0.0354 | +0.17 | +0.8653 | +0.7621 | +0.3305 | +0.0563 | +
Test: SLJ | +0.0052 | +0.0339 | +0.15 | +0.8787 | +0.7751 | +0.2478 | +0.0157 | +
Test: BMT | +-0.0388 | +0.0300 | +-1.29 | +0.1967 | +0.6908 | +0.1968 | +0.0000 | +
Test: Run & a1 | +-0.0561 | +0.0783 | +-0.72 | +0.4735 | ++ | + | + |
Test: Star_r & a1 | +0.2370 | +0.0801 | +2.96 | +0.0031 | ++ | + | + |
Test: S20_r & a1 | +0.1364 | +0.0786 | +1.74 | +0.0826 | ++ | + | + |
Test: SLJ & a1 | +0.0908 | +0.0774 | +1.17 | +0.2409 | ++ | + | + |
Test: BMT & a1 | +0.5593 | +0.0716 | +7.81 | +<1e-14 | ++ | + | + |
Test: Run & Sex: Girls | +-0.5327 | +0.0448 | +-11.90 | +<1e-31 | ++ | + | + |
Test: Star_r & Sex: Girls | +-0.1568 | +0.0461 | +-3.40 | +0.0007 | ++ | + | + |
Test: S20_r & Sex: Girls | +-0.3460 | +0.0449 | +-7.70 | +<1e-13 | ++ | + | + |
Test: SLJ & Sex: Girls | +-0.4140 | +0.0447 | +-9.26 | +<1e-19 | ++ | + | + |
Test: BMT & Sex: Girls | +-0.6955 | +0.0414 | +-16.80 | +<1e-62 | ++ | + | + |
Test: Run & a1 & Sex: Girls | +-0.1455 | +0.1523 | +-0.95 | +0.3396 | ++ | + | + |
Test: Star_r & a1 & Sex: Girls | +-0.2826 | +0.1572 | +-1.80 | +0.0722 | ++ | + | + |
Test: S20_r & a1 & Sex: Girls | +-0.1106 | +0.1530 | +-0.72 | +0.4695 | ++ | + | + |
Test: SLJ & a1 & Sex: Girls | +-0.1231 | +0.1516 | +-0.81 | +0.4167 | ++ | + | + |
Test: BMT & a1 & Sex: Girls | +0.0311 | +0.1404 | +0.22 | +0.8248 | ++ | + | + |
a1 | ++ | + | + | + | + | 0.2833 | ++ |
Sex: Girls | ++ | + | + | + | + | 0.1434 | ++ |
Test: Run | ++ | + | + | + | 0.7443 | +0.3702 | +0.0548 | +
Residual | +0.5201 | ++ | + | + | + | + | + |
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.
+Row | +name | +dof | +deviance | +AIC | +AICc | +BIC | +
---|---|---|---|---|---|---|
+ | Symbol | +Int64 | +Float64 | +Float64 | +Float64 | +Float64 | +
1 | +m1 | +69 | +24651.0 | +24789.0 | +24790.0 | +25284.2 | +
2 | +m2 | +69 | +24646.9 | +24784.9 | +24785.9 | +25280.0 | +
3 | +m2_nested | +69 | +24646.9 | +24784.9 | +24785.9 | +25280.0 | +
In prinicple, 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.
+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.
+ + + + + + +Here we introduce most of the commands available in the MixedModels.jl package that allow the immediated 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.
++ 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))
+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 | +25840.65375540981 | +
Optimizer settings | ++ |
Optimizer (from NLopt) | +LN_BOBYQA |
+
Lower bounds | +[0.0, -Inf, -Inf, -Inf, -Inf, 0.0, -Inf, -Inf, -Inf, 0.0, -Inf, -Inf, 0.0, -Inf, 0.0, 0.0, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, 0.0, -Inf, -Inf, -Inf, -Inf, -Inf, 0.0, -Inf, -Inf, -Inf, -Inf, 0.0, -Inf, -Inf, -Inf, 0.0, -Inf, -Inf, 0.0, -Inf, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] | +
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 |
+[0.75, 1.0, 1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 0.75, 1.0, 0.75, 0.75, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 0.75, 1.0, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75] | +
maxfeval |
+-1 | +
maxtime |
+-1.0 | +
Result | ++ |
Function evaluations | +2161 | +
Final parameter vector | +[1.1503, 0.2062, 0.0009, 0.0507, -0.4821, 1.486, -0.6797, -0.0491, 0.2518, 1.1034, -0.5191, -0.1336, 0.9974, -0.4436, 1.2717, 0.4162, -0.0382, -0.0488, -0.1159, -0.2032, 0.3417, -0.1254, 0.647, -0.2549, 0.1256, -0.0254, -0.092, 0.1178, 0.5837, -0.5055, 0.0597, 0.0542, 0.1363, 0.326, -0.4193, -0.4, -0.0256, 0.3874, 0.1186, -0.1718, 0.0037, -0.001, 0.0, 0.0256, 0.1241, 0.0132, 0.0659, 0.0] | +
Final objective value | +24651.0132 | +
Return code | +FTOL_REACHED |
+
-12325.506614999245
+24651.01322999849
+24651.01322999849
+5746.176138702353
++ julia> coeftable(m1) # StatsBase.coeftable: fixed-effects statiscs;
+ default level=0.95
++ julia> Arrow.write("./data/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
++ | Coef. | +Std. Error | +z | +Pr(> | +
---|---|---|---|---|
(Intercept) | +-0.0133171 | +0.0191925 | +-0.69 | +0.4878 | +
Test: Star_r | +0.0088395 | +0.0374641 | +0.24 | +0.8135 | +
Test: S20_r | +0.0091599 | +0.029254 | +0.31 | +0.7542 | +
Test: SLJ | +0.0033674 | +0.0299438 | +0.11 | +0.9105 | +
Test: BMT | +-0.0382907 | +0.0295746 | +-1.29 | +0.1954 | +
a1 | +0.195545 | +0.0544488 | +3.59 | +0.0003 | +
Sex: Girls | +-0.428744 | +0.0311858 | +-13.75 | +<1e-42 | +
Test: Star_r & a1 | +0.28979 | +0.0887481 | +3.27 | +0.0011 | +
Test: S20_r & a1 | +-0.0783964 | +0.0817053 | +-0.96 | +0.3373 | +
Test: SLJ & a1 | +-0.0631287 | +0.0768871 | +-0.82 | +0.4116 | +
Test: BMT & a1 | +0.471317 | +0.0851499 | +5.54 | +<1e-07 | +
Test: Star_r & Sex: Girls | +0.375601 | +0.0510844 | +7.35 | +<1e-12 | +
Test: S20_r & Sex: Girls | +-0.186727 | +0.0475087 | +-3.93 | +<1e-04 | +
Test: SLJ & Sex: Girls | +-0.0695118 | +0.0445172 | +-1.56 | +0.1184 | +
Test: BMT & Sex: Girls | +-0.281233 | +0.0495206 | +-5.68 | +<1e-07 | +
a1 & Sex: Girls | +-0.126657 | +0.104826 | +-1.21 | +0.2269 | +
Test: Star_r & a1 & Sex: Girls | +-0.138207 | +0.175584 | +-0.79 | +0.4312 | +
Test: S20_r & a1 & Sex: Girls | +0.175298 | +0.163336 | +1.07 | +0.2832 | +
Test: SLJ & a1 & Sex: Girls | +-0.0154409 | +0.152985 | +-0.10 | +0.9196 | +
Test: BMT & a1 & Sex: Girls | +0.154782 | +0.170157 | +0.91 | +0.3630 | +
20-element Vector{Float64}:
+ -0.01331706655216196
+ 0.00883950028956812
+ 0.009159901903096396
+ 0.0033673954396057993
+ -0.03829070122404062
+ 0.19554482244322863
+ -0.4287443626724418
+ 0.28979031920970977
+ -0.07839639356634587
+ -0.06312871281837554
+ 0.47131713732193803
+ 0.37560077961217664
+ -0.18672725220576128
+ -0.06951177262870464
+ -0.28123301997348354
+ -0.1266566750821929
+ -0.13820729140481497
+ 0.17529825709418737
+ -0.015440927030435017
+ 0.15478201529135271
+20-element Vector{Float64}:
+ -0.01331706655216196
+ 0.00883950028956812
+ 0.009159901903096396
+ 0.0033673954396057993
+ -0.03829070122404062
+ 0.19554482244322863
+ -0.4287443626724418
+ 0.28979031920970977
+ -0.07839639356634587
+ -0.06312871281837554
+ 0.47131713732193803
+ 0.37560077961217664
+ -0.18672725220576128
+ -0.06951177262870464
+ -0.28123301997348354
+ -0.1266566750821929
+ -0.13820729140481497
+ 0.17529825709418737
+ -0.015440927030435017
+ 0.15478201529135271
+20-element Vector{Float64}:
+ -0.01331706655216196
+ 0.00883950028956812
+ 0.009159901903096396
+ 0.0033673954396057993
+ -0.03829070122404062
+ 0.19554482244322863
+ -0.4287443626724418
+ 0.28979031920970977
+ -0.07839639356634587
+ -0.06312871281837554
+ 0.47131713732193803
+ 0.37560077961217664
+ -0.18672725220576128
+ -0.06951177262870464
+ -0.28123301997348354
+ -0.1266566750821929
+ -0.13820729140481497
+ 0.17529825709418737
+ -0.015440927030435017
+ 0.15478201529135271
+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"
+20×20 Matrix{Float64}:
+ 0.000368351 2.45468e-5 -1.72453e-5 … 3.89056e-6 2.86117e-6
+ 2.45468e-5 0.00140356 -0.000421444 -9.31972e-7 -1.04343e-7
+ -1.72453e-5 -0.000421444 0.000855795 2.05736e-5 -6.92183e-7
+ -2.68168e-5 5.45073e-5 -0.000468396 -2.44359e-5 1.37239e-6
+ -0.000143239 3.24513e-5 -7.84251e-6 1.42159e-6 -1.19387e-5
+ -4.42745e-5 -7.46407e-5 2.77403e-5 … -6.02716e-5 -5.09392e-5
+ -4.43171e-5 6.33573e-5 4.36499e-5 -7.38976e-5 0.000282215
+ -2.73546e-5 -0.000442542 0.000212864 8.70091e-6 2.62083e-6
+ 3.45871e-6 0.000213247 -0.000395414 -0.000104376 8.17941e-6
+ -2.00311e-5 7.81719e-8 0.000175323 0.000158956 -6.11846e-5
+ 7.30133e-5 -2.67094e-5 2.85589e-5 … -5.7891e-5 8.27717e-5
+ 6.89346e-6 8.80541e-7 -2.06783e-6 8.67367e-6 -0.000108136
+ 5.13451e-6 -1.96787e-6 7.69769e-6 0.000673738 0.000117904
+ -5.23405e-6 -1.56257e-6 -4.65589e-6 -0.00136392 0.000576753
+ -1.14233e-6 6.34682e-7 -9.06233e-7 0.000577419 -0.00171453
+ -4.58649e-6 -1.20769e-5 3.57658e-7 … 0.000152536 -0.00391907
+ -1.2815e-5 -2.83723e-5 2.85677e-5 -0.000139039 0.00179084
+ 4.96066e-7 2.86311e-5 -4.70739e-5 -0.0115807 -0.00195294
+ 3.89056e-6 -9.31972e-7 2.05736e-5 0.0234043 -0.00998529
+ 2.86117e-6 -1.04343e-7 -6.92183e-7 -0.00998529 0.0289533
+20×20 Matrix{Float64}:
+ 1.0 0.0341389 -0.0307153 … 0.00132505 0.000876119
+ 0.0341389 1.0 -0.384539 -0.000162607 -1.63681e-5
+ -0.0307153 -0.384539 1.0 0.00459702 -0.000139055
+ -0.0466626 0.0485884 -0.534714 -0.00533426 0.000269353
+ -0.252355 0.0292886 -0.00906467 0.000314202 -0.0023724
+ -0.0423677 -0.0365909 0.0174156 … -0.00723564 -0.00549812
+ -0.074043 0.0542282 0.0478456 -0.0154891 0.0531832
+ -0.0160598 -0.133101 0.0819894 0.000640852 0.000173552
+ 0.00220563 0.0696654 -0.165431 -0.00835027 0.000588332
+ -0.0135744 2.71383e-5 0.077947 0.0135137 -0.00467669
+ 0.0446773 -0.00837267 0.0114649 … -0.00444405 0.00571279
+ 0.00703102 0.000460094 -0.0013837 0.00110986 -0.0124404
+ 0.00563112 -0.00110563 0.00553863 0.0926979 0.014585
+ -0.00612604 -0.000936911 -0.00357512 -0.200269 0.0761401
+ -0.00120191 0.000342101 -0.00062556 0.0762178 -0.203474
+ -0.00227972 -0.00307519 0.000116631 … 0.00951161 -0.219717
+ -0.0038028 -0.00431315 0.00556167 -0.0051761 0.0599409
+ 0.000158243 0.00467887 -0.00985172 -0.463453 -0.0702678
+ 0.00132505 -0.000162607 0.00459702 1.0 -0.383587
+ 0.000876119 -1.63681e-5 -0.000139055 -0.383587 1.0
+20-element Vector{Float64}:
+ 0.019192465300358272
+ 0.03746407305661245
+ 0.02925397597009649
+ 0.029943795561659332
+ 0.02957456318772991
+ 0.054448757362360926
+ 0.031185791059940196
+ 0.08874814582057454
+ 0.08170531215279139
+ 0.07688713611571989
+ 0.0851499496378772
+ 0.05108438663261352
+ 0.04750873665393441
+ 0.044517155988056785
+ 0.0495206365951175
+ 0.10482604370023066
+ 0.17558419241258222
+ 0.16333633596447028
+ 0.1529846844619703
+ 0.17015682547428446
+(:formula, :reterms, :Xymat, :feterm, :sqrtwts, :parmap, :dims, :A, :L, :optsum, :θ, :theta, :β, :beta, :βs, :betas, :λ, :lambda, :stderror, :σ, :sigma, :σs, :sigmas, :σρs, :sigmarhos, :b, :u, :lowerbd, :X, :y, :corr, :vcov, :PCA, :rePCA, :objective, :pvalues)
+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
+
+
++ | Column | +Variance | +Std.Dev | +Corr. | ++ | + | + | + | + |
---|---|---|---|---|---|---|---|---|---|
Child | +(Intercept) | +0.34903980 | +0.59079590 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.59364220 | +0.77048180 | ++0.14 | ++ | + | + | + | + |
+ | Test: S20_r | +0.44304305 | +0.66561479 | ++0.00 | +-0.52 | ++ | + | + | + |
+ | Test: SLJ | +0.33481089 | +0.57862846 | ++0.05 | +-0.04 | +-0.37 | ++ | + | + |
+ | Test: BMT | +0.56120773 | +0.74913799 | +-0.33 | ++0.13 | +-0.17 | +-0.25 | ++ | + |
School | +(Intercept) | +0.04569734 | +0.21376937 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.11081338 | +0.33288644 | +-0.06 | ++ | + | + | + | + |
+ | Test: S20_r | +0.10763442 | +0.32807686 | +-0.08 | +-0.39 | ++ | + | + | + |
+ | Test: SLJ | +0.10314745 | +0.32116576 | +-0.19 | ++0.21 | +-0.80 | ++ | + | + |
+ | Test: BMT | +0.09796213 | +0.31298902 | +-0.33 | +-0.02 | ++0.13 | +-0.38 | ++ | + |
+ | a1 | +0.07971464 | +0.28233780 | ++0.62 | +-0.20 | ++0.11 | +-0.61 | ++0.45 | ++ |
+ | Sex: Girls | +0.02066903 | +0.14376729 | +-0.45 | ++0.45 | ++0.31 | +-0.27 | +-0.15 | +-0.37 | +
Cohort | +(Intercept) | +0.00017336 | +0.01316656 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.00405975 | +0.06371614 | +. | ++ | + | + | + | + |
+ | Test: S20_r | +0.00004631 | +0.00680502 | +. | +. | ++ | + | + | + |
+ | Test: SLJ | +0.00114631 | +0.03385719 | +. | +. | +. | ++ | + | + |
+ | Test: BMT | +0.00000000 | +0.00000000 | +. | +. | +. | +. | ++ | + |
Residual | ++ | 0.26377271 | +0.51358807 | ++ | + | + | + | + | + |
+ | Column | +Variance | +Std.Dev | +Corr. | ++ | + | + | + | + |
---|---|---|---|---|---|---|---|---|---|
Child | +Test: Run | +0.5626315 | +0.7500877 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.6354871 | +0.7971744 | ++0.50 | ++ | + | + | + | + |
+ | Test: S20_r | +0.5894630 | +0.7677650 | ++0.56 | ++0.64 | ++ | + | + | + |
+ | Test: SLJ | +0.6093249 | +0.7805927 | ++0.55 | ++0.60 | ++0.72 | ++ | + | + |
+ | Test: BMT | +0.4858396 | +0.6970219 | ++0.19 | ++0.40 | ++0.37 | ++0.49 | ++ | + |
School | +Test: Run | +0.1370367 | +0.3701847 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.0962842 | +0.3102970 | ++0.53 | ++ | + | + | + | + |
+ | Test: S20_r | +0.1092459 | +0.3305237 | ++0.47 | ++0.48 | ++ | + | + | + |
+ | Test: SLJ | +0.0613916 | +0.2477733 | ++0.47 | ++0.76 | ++0.41 | ++ | + | + |
+ | Test: BMT | +0.0387555 | +0.1968641 | ++0.15 | ++0.38 | ++0.18 | ++0.02 | ++ | + |
+ | a1 | +0.0802430 | +0.2832720 | ++0.58 | ++0.47 | ++0.56 | +-0.05 | ++0.65 | ++ |
+ | Sex: Girls | +0.0205488 | +0.1433486 | +-0.63 | +-0.27 | ++0.06 | +-0.28 | +-0.58 | +-0.37 | +
Cohort | +Test: Run | +0.0030067 | +0.0548338 | ++ | + | + | + | + | + |
+ | Test: Star_r | +0.0004347 | +0.0208494 | +. | ++ | + | + | + | + |
+ | Test: S20_r | +0.0031725 | +0.0563249 | +. | +. | ++ | + | + | + |
+ | Test: SLJ | +0.0002452 | +0.0156604 | +. | +. | +. | ++ | + | + |
+ | Test: BMT | +0.0000000 | +0.0000000 | +. | +. | +. | +. | ++ | + |
Residual | ++ | 0.2619146 | +0.5117760 | ++ | + | + | + | + | + |
(Child = (var"(Intercept)" = 0.5907959029225862, var"Test: Star_r" = 0.7704817994480245, var"Test: S20_r" = 0.6656147888960222, var"Test: SLJ" = 0.5786284552355224, var"Test: BMT" = 0.7491379937832613), School = (var"(Intercept)" = 0.2137693688364067, var"Test: Star_r" = 0.33288643566454384, var"Test: S20_r" = 0.3280768567505694, var"Test: SLJ" = 0.3211657642930787, var"Test: BMT" = 0.31298901749687913, a1 = 0.28233780381354795, var"Sex: Girls" = 0.143767290743394), Cohort = (var"(Intercept)" = 0.013166555574091778, var"Test: Star_r" = 0.06371614416899145, var"Test: S20_r" = 0.006805016452716729, var"Test: SLJ" = 0.03385719231820434, var"Test: BMT" = 0.0))
+48-element Vector{Float64}:
+ 1.1503302598343477
+ 0.2062190787609727
+ 0.0008666159469680149
+ 0.050734934180392224
+ -0.48211325954936624
+ 1.4859529039387687
+ -0.6797347055373657
+ -0.04906500731877601
+ 0.25177147495411384
+ 1.103448875092165
+ -0.5190786027370973
+ -0.13364340349638246
+ 0.9974426672137859
+ ⋮
+ -0.02563052910620881
+ 0.3873988146494795
+ 0.118601913906712
+ -0.1718026965130414
+ 0.0036579916237788343
+ -0.0010172580519802405
+ 0.0
+ 0.02563641220215602
+ 0.12406079377832567
+ 0.013249950288257796
+ 0.06592285532786576
+ 0.0
+rows | +Child | +School | +Cohort | +fixed | +
---|---|---|---|---|
9930 | +BlkDiag | ++ | + | + |
3101 | +Sparse | +BlkDiag | ++ | + |
45 | +Dense | +Dense | +BlkDiag/Dense | ++ |
21 | +Dense | +Dense | +Dense | +Dense | +
48-element Vector{Float64}:
+ 1.46565629578844
+ 0.7825331498310524
+ 0.8427270268123437
+ 0.843701714893452
+ 0.26423363532022687
+ 1.3468316700709517
+ 0.6159105518922656
+ 0.5685881881743816
+ 0.4806008412988357
+ 1.0775237645019875
+ 0.5454461221597267
+ 0.21582202447173554
+ 0.9968896467848608
+ ⋮
+ -0.01778801467231765
+ 0.31874007820888184
+ 0.13868543043388296
+ -0.18628651012633302
+ 0.0
+ -0.00019459500089828938
+ 0.0
+ 0.10714419938783641
+ 0.040739228366369604
+ 0.11005771477604791
+ 0.03060005959856037
+ 0.0
+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
.
Some plotting functions are currently available from the MixedModelsMakie
package or via custom functions.
3-element Vector{Array{Float64, 3}}:
+ [0.05575383274003557 0.00863915213999665 … 0.0031232443363528556 -0.021138218415038885; 0.00863915213999665 0.29001266660027136 … -0.0070474409935052055 0.022922985962176668; … ; 0.0031232443363528556 -0.0070474409935052055 … 0.20790206584505386 -0.0852354729075977; -0.021138218415038885 0.022922985962176668 … -0.0852354729075977 0.26330097819700116;;; 0.06213230185040765 0.010826308858088864 … 0.003929097771018621 -0.028409878110103502; 0.010826308858088864 0.29938324063036736 … -0.005483186816479629 0.023531240948500384; … ; 0.003929097771018621 -0.005483186816479629 … 0.21199432719260117 -0.08620203115574254; -0.028409878110103502 0.023531240948500384 … -0.08620203115574254 0.2731885630201976;;; 0.05928754515087498 0.009999345149668016 … 0.003214830309751238 -0.02513168722029122; 0.009999345149668016 0.29438206427015273 … -0.006437102215677196 0.022965997644015704; … ; 0.003214830309751238 -0.006437102215677196 … 0.21076926571106885 -0.08552674167274826; -0.02513168722029122 0.022965997644015704 … -0.08552674167274826 0.26861055510165105;;; … ;;; 0.0907210265033368 0.010497457003751263 … -0.0060407174687821434 -0.023131077664248153; 0.010497457003751263 0.2970544531740206 … -0.0069478640031876895 0.024285798673900034; … ; -0.0060407174687821434 -0.0069478640031876895 … 0.2099626452419368 -0.08344670225122619; -0.023131077664248153 0.024285798673900034 … -0.08344670225122619 0.2685457288945719;;; 0.09719108603781947 0.007906212183110295 … -0.001320199364202595 -0.022070542259081175; 0.007906212183110295 0.29347487103474784 … -0.007201306571266091 0.02468292835627069; … ; -0.001320199364202595 -0.007201306571266091 … 0.20924944238722887 -0.08407452186026657; -0.022070542259081175 0.02468292835627069 … -0.08407452186026657 0.2654490089917222;;; 0.08379719940391851 0.007847588198883811 … -0.00028220224056425754 -0.01854945143543425; 0.007847588198883811 0.29055593804014934 … -0.007421355885115292 0.02415358501664865; … ; -0.00028220224056425754 -0.007421355885115292 … 0.20794862679330195 -0.0842550154482498; -0.01854945143543425 0.02415358501664865 … -0.0842550154482498 0.2624831493252433]
+ [0.03738099438108044 -0.001778150236782388 … 0.030005959239663836 -0.012422815955098213; -0.001778150236782388 0.09046590526227734 … -0.012438199387792876 0.017400101424582894; … ; 0.030005959239663836 -0.012438199387792876 … 0.06539708247917235 -0.012780944067289027; -0.012422815955098213 0.017400101424582894 … -0.012780944067289027 0.01833845610541896;;; 0.032174800419166066 0.00021552464206599673 … 0.02631804193514885 -0.010643849286266155; 0.00021552464206599673 0.07636177622427075 … -0.007781645246191157 0.014216263238686322; … ; 0.02631804193514885 -0.007781645246191157 … 0.057233707504033296 -0.010780210452738961; -0.010643849286266155 0.014216263238686322 … -0.010780210452738961 0.015905591438488714;;; 0.029307729689901536 -0.0031500145971753725 … 0.02159243212220725 -0.01044773671375439; -0.0031500145971753725 0.0877447231001567 … -0.01330551483540402 0.018117822615212966; … ; 0.02159243212220725 -0.01330551483540402 … 0.05168567480761258 -0.010391773146284176; -0.01044773671375439 0.018117822615212966 … -0.010391773146284176 0.01678423124230707;;; … ;;; 0.026752237575879263 0.0029124911661105326 … 0.01908410221796392 -0.006843494296351528; 0.0029124911661105326 0.07610913276468977 … -0.005653623936266287 0.012914109119085188; … ; 0.01908410221796392 -0.005653623936266287 … 0.05072498368147368 -0.006134494719465586; -0.006843494296351528 0.012914109119085188 … -0.006134494719465586 0.013714332731263679;;; 0.02722117937531621 0.0035247368383708056 … 0.020363601389567814 -0.0061486862691171144; 0.0035247368383708056 0.07602932922491482 … -0.005124378717151021 0.012614623673192117; … ; 0.020363601389567814 -0.005124378717151021 … 0.05286205808611119 -0.00631660908074008; -0.0061486862691171144 0.012614623673192117 … -0.00631660908074008 0.013738866013923129;;; 0.04333514144377728 -0.0030526201099861485 … 0.0364639269042496 -0.013644784026381104; -0.0030526201099861485 0.09939148209978196 … -0.01520556786587638 0.0191560598024237; … ; 0.0364639269042496 -0.01520556786587638 … 0.07495035348527168 -0.014381992472318614; -0.013644784026381104 0.0191560598024237 … -0.014381992472318614 0.01951413892155994]
+ [0.00016287570134958826 6.856559889853957e-6 … -4.739951921609176e-6 0.0; 6.856559889853957e-6 0.0023548313779246824 … -0.00015351717222254004 0.0; … ; -4.739951921609176e-6 -0.00015351717222254004 … 0.0008614762660391023 0.0; 0.0 0.0 … 0.0 0.0;;; 0.00015950907262546907 8.701350327885633e-6 … -6.024504894509734e-6 0.0; 8.701350327885633e-6 0.0020423106983111475 … -0.0001673374995288987 0.0; … ; -6.024504894509734e-6 -0.0001673374995288987 … 0.0007931845117108109 0.0; 0.0 0.0 … 0.0 0.0;;; 0.00016145306769277504 7.573327674888944e-6 … -5.152636375745283e-6 0.0; 7.573327674888944e-6 0.0022393880860264387 … -0.00015896888705918904 0.0; … ; -5.152636375745283e-6 -0.00015896888705918904 … 0.0008333330451138073 0.0; 0.0 0.0 … 0.0 0.0;;; 0.00016084702419261226 8.35246981386616e-6 … -5.4742810949165876e-6 0.0; 8.35246981386616e-6 0.0021616288990247254 … -0.00016342062464208147 0.0; … ; -5.4742810949165876e-6 -0.00016342062464208147 … 0.0008213033177068655 0.0; 0.0 0.0 … 0.0 0.0;;; 0.00015845270342141037 8.619670728536786e-6 … -6.900073371582094e-6 0.0; 8.619670728536786e-6 0.0019601068318098897 … -0.0001690666309293983 0.0; … ; -6.900073371582094e-6 -0.0001690666309293983 … 0.0007719484491124605 0.0; 0.0 0.0 … 0.0 0.0;;; 0.00015768719984250653 7.385065727666733e-6 … -7.336092614509354e-6 0.0; 7.385065727666733e-6 0.0019316205603731977 … -0.00016776153068661235 0.0; … ; -7.336092614509354e-6 -0.00016776153068661235 … 0.0007548457845576088 0.0; 0.0 0.0 … 0.0 0.0;;; 0.00015975999081971895 7.948268432845527e-6 … -6.413248551251487e-6 0.0; 7.948268432845527e-6 0.002074914216290986 … -0.00016692868092659584 0.0; … ; -6.413248551251487e-6 -0.00016692868092659584 … 0.0007975022546437895 0.0; 0.0 0.0 … 0.0 0.0;;; 0.00015744703326980335 8.402243512073218e-6 … -7.90817159440362e-6 0.0; 8.402243512073218e-6 0.0018868038496581622 … -0.00017276721709532413 0.0; … ; -7.90817159440362e-6 -0.00017276721709532413 … 0.0007513634534431922 0.0; 0.0 0.0 … 0.0 0.0;;; 0.00015564965231731795 9.68631259477342e-6 … -7.697020113052748e-6 0.0; 9.68631259477342e-6 0.0017866573072894876 … -0.0001720968882090284 0.0; … ; -7.697020113052748e-6 -0.0001720968882090284 … 0.0007251764673322595 0.0; 0.0 0.0 … 0.0 0.0]
+3-element Vector{Array{Float64, 3}}:
+ [0.17116315582859484 0.027498862453680717 … 0.034104094305388286 -0.007244970029438276; 0.027498862453680717 0.17384103461076678 … 0.03748145360126994 0.018816089248565385; … ; 0.034104094305388286 0.03748145360126994 … 0.16031916639849703 0.028749779922481; -0.007244970029438276 0.018816089248565385 … 0.028749779922481 0.1607955523119319;;; 0.18238552619993495 0.03427097227300673 … 0.04141109934983114 -0.009148439621905625; 0.03427097227300673 0.1856359829925568 … 0.04863182300871331 0.021358489400357405; … ; 0.04141109934983114 0.04863182300871331 … 0.17268807116668955 0.031166483369606347; -0.009148439621905625 0.021358489400357405 … 0.031166483369606347 0.16324754394182825;;; 0.17621626015158762 0.031072180113851766 … 0.038125649002372385 -0.007995567564654946; 0.031072180113851766 0.1803469801138963 … 0.04349850263409247 0.020097760108707865; … ; 0.038125649002372385 0.04349850263409247 … 0.16701836305224055 0.030006687972602957; -0.007995567564654946 0.020097760108707865 … 0.030006687972602957 0.16197496984573145;;; … ;;; 0.21740912937655923 0.07090924481039461 … 0.06741379636911818 0.021486080664454548; 0.07090924481039461 0.22153412311552326 … 0.07279946779268884 0.050948373881341553; … ; 0.06741379636911818 0.07279946779268884 … 0.18373137691509117 0.05058489384545355; 0.021486080664454548 0.050948373881341553 … 0.05058489384545355 0.18633185886033954;;; 0.22426025206639164 0.0776902430712113 … 0.07687730598105631 0.03219132055270273; 0.0776902430712113 0.22462307925799707 … 0.07888460223906732 0.05867761080950808; … ; 0.07687730598105631 0.07888460223906732 … 0.19238983868372397 0.060896061536589774; 0.03219132055270273 0.05867761080950808 … 0.060896061536589774 0.19517714451586868;;; 0.2061169154182078 0.061463168353437944 … 0.062212999195302 0.022103225527677538; 0.061463168353437944 0.2073660464425763 … 0.0640017165502183 0.047830931525945246; … ; 0.062212999195302 0.0640017165502183 … 0.1792425246311685 0.05164496814401823; 0.022103225527677538 0.047830931525945246 … 0.05164496814401823 0.18684232972258574]
+ [0.10768178505825274 0.048400872742620095 … 0.04766684633075814 -0.027767663870375308; 0.048400872742620095 0.07949976129667788 … 0.034944549376798424 -0.010301353681294811; … ; 0.04766684633075814 0.034944549376798424 … 0.06575697359151865 -0.012753936623032553; -0.027767663870375308 -0.010301353681294811 … -0.012753936623032553 0.01822072725963938;;; 0.08790095967031233 0.04039487399830374 … 0.04006050669017243 -0.022694796915798992; 0.04039487399830374 0.06914919686791131 … 0.03207348888389284 -0.008419559170898776; … ; 0.04006050669017243 0.03207348888389284 … 0.05750637187929682 -0.010759600050894616; -0.022694796915798992 -0.008419559170898776 … -0.010759600050894616 0.01580107265173481;;; 0.0963598911033401 0.035885072470449086 … 0.038333060297426495 -0.02541792355752846; 0.035885072470449086 0.0630834010975428 … 0.02477182157087897 -0.007236741864081645; … ; 0.038333060297426495 0.02477182157087897 … 0.0519930767818036 -0.010389841513651113; -0.02541792355752846 -0.007236741864081645 … -0.010389841513651113 0.016683879213529238;;; … ;;; 0.07440697397002378 0.030529015222540118 … 0.029374057723458375 -0.017031449145314057; 0.030529015222540118 0.06267553791989461 … 0.023497129411826363 -0.004055529443914723; … ; 0.029374057723458375 0.023497129411826363 … 0.051035223116704755 -0.006132115409660271; -0.017031449145314057 -0.004055529443914723 … -0.006132115409660271 0.013629228762574913;;; 0.07355978268459702 0.03053833755418159 … 0.030108406727685198 -0.01606992180353426; 0.03053833755418159 0.06346698084161012 … 0.024774759577003626 -0.0033909215418409006; … ; 0.030108406727685198 0.024774759577003626 … 0.053155675557087244 -0.006318319173994548; -0.01606992180353426 -0.0033909215418409006 … -0.006318319173994548 0.013652196613227967;;; 0.12472878923414969 0.057966342588725205 … 0.05756315402857382 -0.03100070689366643; 0.057966342588725205 0.09051345642227474 … 0.042031396044692 -0.011772108325862396; … ; 0.05756315402857382 0.042031396044692 … 0.0753703800124635 -0.01436086970585532; -0.03100070689366643 -0.011772108325862396 … -0.01436086970585532 0.019396316425718988]
+ [0.0018739464805060115 3.099298859387236e-5 … 2.2116659335547737e-5 0.0; 3.099298859387236e-5 0.0003982925350932456 … 5.662327674244246e-6 0.0; … ; 2.2116659335547737e-5 5.662327674244246e-6 … 0.0002316450289147418 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0016451178801496406 3.732469788303523e-5 … 2.7075622912847006e-5 0.0; 3.732469788303523e-5 0.00038699606081539056 … 7.752861937193192e-6 0.0; … ; 2.7075622912847006e-5 7.752861937193192e-6 … 0.0002272553166207261 0.0; 0.0 0.0 … 0.0 0.0;;; 0.001788360817926057 3.246216882159844e-5 … 2.409075582915546e-5 0.0; 3.246216882159844e-5 0.0003941993142896808 … 6.47076235799776e-6 0.0; … ; 2.409075582915546e-5 6.47076235799776e-6 … 0.0002298593116535934 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0017314736634095084 3.439352491212058e-5 … 2.510100552477303e-5 0.0; 3.439352491212058e-5 0.00039210347603502567 … 6.686612028470298e-6 0.0; … ; 2.510100552477303e-5 6.686612028470298e-6 … 0.0002292053931970069 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0015869175336522604 3.8916114227008614e-5 … 2.7858034993426766e-5 0.0; 3.8916114227008614e-5 0.00038392382166180716 … 8.107240306589517e-6 0.0; … ; 2.7858034993426766e-5 8.107240306589517e-6 … 0.00022584769979068763 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0015689110909946945 4.003621092371261e-5 … 2.8522941942561645e-5 0.0; 4.003621092371261e-5 0.00038047593288191215 … 8.919058051767636e-6 0.0; … ; 2.8522941942561645e-5 8.919058051767636e-6 … 0.00022463172084040967 0.0; 0.0 0.0 … 0.0 0.0;;; 0.001670592740351977 3.6961837280171235e-5 … 2.6093380656973352e-5 0.0; 3.6961837280171235e-5 0.00038740854803262644 … 7.506517475684638e-6 0.0; … ; 2.6093380656973352e-5 7.506517475684638e-6 … 0.00022749979656451797 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0015337205158612683 4.093131889936314e-5 … 2.9075417115821492e-5 0.0; 4.093131889936314e-5 0.000380334357158905 … 8.821163224637285e-6 0.0; … ; 2.9075417115821492e-5 8.821163224637285e-6 … 0.00022435470926792088 0.0; 0.0 0.0 … 0.0 0.0;;; 0.001453753506796497 4.288550370752128e-5 … 3.0677599882156927e-5 0.0; 4.288550370752128e-5 0.00037427960644329807 … 9.99134452808877e-6 0.0; … ; 3.0677599882156927e-5 9.99134452808877e-6 … 0.00022244049186675522 0.0; 0.0 0.0 … 0.0 0.0]
+They are hard to look at. Let’s take pictures.
+These are just teasers. We will pick this up in a separate tutorial. Enjoy!
+ + +Phillip Alday, Douglas Bates, and Reinhold Kliegl
+2023-08-21
+Arrow.Table with 525126 rows, 7 columns, and schema:
+ :Cohort String
+ :School String
+ :Child String
+ :Sex String
+ :age Float64
+ :Test String
+ :score Float64
+Row | +variable | +mean | +min | +median | +max | +nmissing | +eltype | +
---|---|---|---|---|---|---|---|
+ | Symbol | +Union… | +Any | +Union… | +Any | +Int64 | +DataType | +
1 | +Cohort | ++ | 2011 | ++ | 2019 | +0 | +String | +
2 | +School | ++ | S100043 | ++ | S800200 | +0 | +String | +
3 | +Child | ++ | C002352 | ++ | C117966 | +0 | +String | +
4 | +Sex | ++ | female | ++ | male | +0 | +String | +
5 | +age | +8.56073 | +7.99452 | +8.55852 | +9.10609 | +0 | +Float64 | +
6 | +Test | ++ | BPT | ++ | Star_r | +0 | +String | +
7 | +score | +226.141 | +1.14152 | +4.65116 | +1530.0 | +0 | +Float64 | +
let
+ fdensity = Figure(; resolution=(1000, 500))
+ axs = Axis(fdensity[1, 1])
+ tdf = filter(:Test => ==(test), df)
+ colors = Makie.cgrad(:PuOr_4, 2; categorical=true, alpha=0.6)
+ if by_sex
+ density!(
+ axs,
+ filter(:Sex => ==("female"), tdf).score;
+ color=colors[1],
+ label="Girls",
+ )
+ density!(
+ axs,
+ filter(:Sex => ==("male"), tdf).score;
+ color=colors[2],
+ label="Boys",
+ )
+ axislegend(axs; position=:lt)
+ else
+ density!(axs, tdf.score)
+ end
+ fdensity
+end