Emotikon: Comparing Original and Transformed Metrics

Published

2026-04-10

In @Fuehner2021 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.

1 Getting the packages and data

Code 
using AlgebraOfGraphics
using Arrow
using CairoMakie
using DataFrames
using DataFrameMacros
using MixedModels
using MixedModelsMakie
using RCall
using SMLP2026: dataset, FITSDIR
using Serialization
using StatsBase

1.1 Data and figure in publication

dat = DataFrame(dataset(:fggk21))
525126×7 DataFrame
525101 rows omitted
Row Cohort School Child Sex age Test score
String String String String Float64 String Float64
1 2013 S100067 C002352 male 7.99452 S20_r 5.26316
2 2013 S100067 C002352 male 7.99452 BPT 3.7
3 2013 S100067 C002352 male 7.99452 SLJ 125.0
4 2013 S100067 C002352 male 7.99452 Star_r 2.47146
5 2013 S100067 C002352 male 7.99452 Run 1053.0
6 2013 S100067 C002353 male 7.99452 S20_r 5.0
7 2013 S100067 C002353 male 7.99452 BPT 4.1
8 2013 S100067 C002353 male 7.99452 SLJ 116.0
9 2013 S100067 C002353 male 7.99452 Star_r 1.76778
10 2013 S100067 C002353 male 7.99452 Run 1089.0
11 2013 S100067 C002354 male 7.99452 S20_r 4.54545
12 2013 S100067 C002354 male 7.99452 BPT 3.9
13 2013 S100067 C002354 male 7.99452 SLJ 111.0
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
525115 2018 S401470 C117964 male 9.10609 Star_r 1.63704
525116 2018 S401470 C117964 male 9.10609 Run 864.0
525117 2018 S401470 C117965 female 9.10609 S20_r 4.65116
525118 2018 S401470 C117965 female 9.10609 BPT 3.8
525119 2018 S401470 C117965 female 9.10609 SLJ 123.0
525120 2018 S401470 C117965 female 9.10609 Star_r 1.52889
525121 2018 S401470 C117965 female 9.10609 Run 1080.0
525122 2018 S800200 C117966 male 9.10609 S20_r 4.54545
525123 2018 S800200 C117966 male 9.10609 BPT 3.8
525124 2018 S800200 C117966 male 9.10609 SLJ 100.0
525125 2018 S800200 C117966 male 9.10609 Star_r 2.18506
525126 2018 S800200 C117966 male 9.10609 Run 990.0 
@transform!(dat, :a1 = :age - 8.5);
select!(groupby(dat, :Test), :, :score => zscore => :zScore);
describe(dat)
9×7 DataFrame
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

1.2 Data and figure with z-scores based on original metric

2 LMMs

2.1 Contrasts

contrasts = Dict(
  :Test => SeqDiffCoding(),
  :Sex => HelmertCoding(),
);

2.2 Formula

f1 = @formula zScore ~
  1 +
  Test * a1 * Sex +
  (1 + Test + a1 + Sex | School) +
  (1 + Test | Child) +
  zerocorr(1 + Test | Cohort);

2.3 Restore LMM m1 from publication

  • Command for fitting LMM m1 = fit(MixedModel, f1, dat, contrasts=contr)
  • Fit statistics for LMM m1: Minimizing 5179 Time: 0 Time: 1:00:38 ( 0.70 s/it)
m1x = LinearMixedModel(f1, dat; contrasts)
restoreoptsum!(m1x, joinpath(FITSDIR, "fggk21_m1_optsum.json"))
┌ Warning: optsum was saved with an older version of MixedModels.jl: consider resaving.
└ @ MixedModels ~/.julia/packages/MixedModels/ACdLy/src/serialization.jl:91
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 
VarCorr(m1x)
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

2.4 Restore new LMM m1_om Star and S20 in original metric

  • Command for fitting LMM m1_om = fit(MixedModel, f1, dat_om, contrasts=contr)
  • Minimizing 10502 Time: 0 Time: 2:09:40 ( 0.74 s/it)
  • Store with: julia> saveoptsum(joinpath(FITSDIR, “fggk21_m1_om_optsum.json”), m1_om)
  • Only for short-term and when desperate: julia> serialize(joinpath(FITSDIR, “m1_om.jls”), m1_om);

2.4.1 … restoreoptsum!()

m1_om = LinearMixedModel(f1, dat; contrasts);
restoreoptsum!(m1_om, joinpath(FITSDIR, "fggk21_m1_om_optsum.json"));

2.4.2 … deserialize()

m1x_om = deserialize(joinpath(FITSDIR, "m1_om.jls"))
VarCorr(m1x_om)

2.5 Residual diagnostics for LMM m1

Residual plots for published LMM

#scatter(fitted(m1x), residuals(m1x)
#qqnorm(m1x)

2.6 Residual diagnostics for LMM m1_om

Residual plots for LMM with Star and Speed in original metric.

#scatter(fitted(m1_om_v2), residuals(m1_om_v2)
#qqnorm(m1_om_v2)

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

Back to top