9  Julia

9.1 Supplementary Material Chapter 6

Code 
using Pkg; #Pkg.activate(".")

using CairoMakie        # graphics back-end
using CategoricalArrays
using DataFrames
using MixedModels
using MixedModelsMakie  # diagnostic plots
using MixedModelsExtras # partial effects
using PrettyTables
using RCall             # call R from Julia
using Parquet

CairoMakie.activate!(; type="png");
Code 
R"""
suppressWarnings(suppressMessages(library(tidyverse)))
load(file="data/simon_raw.rda")
load(file="data/info.rda")
load(file="data/ages.rda")


data <- info |> 
  select(Child, Time, gender) |> 
  left_join(ages |> 
              select(Child, Time, age), multiple = "all") |> 
  left_join(simon_raw, multiple = "all") |> 
  subset(!is.na(RT)) |> 
  mutate(block = case_when(trialnumber %in% -5:0  ~ 0,
                           trialnumber %in%  1:20 ~ 1,
                           trialnumber %in% 21:40 ~ 2,
                           trialnumber %in% 41:60 ~ 3)) |> 
  unite(Condition, c(Colour, Side, Congruence), sep = "_", remove = FALSE)

dat <- data |> 
  subset(trialnumber > 0 & 
           Answer == "correct" & 
           between(RT, 500, 10000))

dat$speed <- 1000/dat$RT
""";
@rget dat;

transform!(dat,
    :Child => categorical => :Child,
    :Time => categorical => :Time,
    :gender => categorical => :gender,
    :Condition  => categorical => :Condition,
    :Colour  => categorical => :Colour,
    :Side  => categorical => :Side,
    :Congruence  => categorical => :Congruence,
    :age => (x -> x .- 9.9) => :a1);

dat = dat[completecases(dat), :];

dat2 = DataFrame(read_parquet("data/emo.parquet"))
dat2 = select!(dat2, Not(:Score))
dat2 = unstack(dat2, :Test, :zScore)
dat2 = transform!(dat2,
    :Child => categorical => :Child,
    :Time => categorical => :Time);
dat2 = dat2[completecases(dat2), :];
Data = leftjoin(dat, dat2, on = [:Child, :Time]);
Data = Data[completecases(Data, :OLB_l), :];
describe(Data)

contrasts = merge(
      Dict(:Congruence => EffectsCoding(base= "incon"; levels=["con", "incon"])),
      Dict(:Colour => EffectsCoding(base= "blue"; levels=["red", "blue"])),
      Dict(:Side => EffectsCoding(base= "left"; levels=["right", "left"])),
      Dict(:gender => EffectsCoding(base= "girl"; levels=["girl", "boy"])),
      Dict(nm => Grouping() for nm in (:Child, :School))
);

9.1.1 Main effects Simon Task

9.1.1.1 Model selection

  • m0: speed ~ 1 + Congruence + Colour + Side + a1 + gender + (1 | Child)
  • m0C1: speed ~ 1 + Congruence * a1 + gender + Colour + Side + (1 | Child)
  • m0C2: speed ~ 1 + Congruence * (a1 + gender) + Colour + Side + (1 | Child)
Code 
f0 = @formula(speed ~  1 + Congruence + Colour + Side + a1 + gender + (1  | Child));
m0 = fit(MixedModel, f0, dat; contrasts);

f0C1 = @formula(speed ~  1 + Congruence * a1 + gender + Colour + Side + (1  | Child));
m0C1 = fit(MixedModel, f0C1, dat; contrasts);

f0C2 = @formula(speed ~  1 + Congruence * (a1 + gender) + Colour + Side + (1  | Child));
m0C2 = fit(MixedModel, f0C2, dat; contrasts);
Code 
lrt0 = gof_summary = let
         nms = [:m0, :m0C1, :m0C2]
         mods = eval.(nms)
         lrt = MixedModels.likelihoodratiotest(m0, m0C1, m0C2)
         DataFrame(;
           name = nms, 
           dof=dof.(mods),
           deviance=deviance.(mods),
           AIC=aic.(mods),
           AICc=aicc.(mods),
           BIC=bic.(mods),
           χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
           χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
           pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
         )
       end
pretty_table(lrt0, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.1: LRT of m0, m0C1, and m0C2
name dof devianc AIC AICc BIC χ² χ²_dof p
m0 8 882.824 898.824 898.831 962.41 . . .
m0C1 9 881.472 899.472 899.481 971.007 1.351 1.0 0.245
m0C2 10 880.731 900.731 900.742 980.214 0.741 1.0 0.389

No improvements of the models fit due to interactions between the Congruence with age and gender.

  • m1: speed ~ 1 + Congruence + Side + Colour + a1 + gender + (1 + Congruence + a1 | Child)
  • m2: speed ~ 1 + Congruence + Side + Colour + a1 + gender + (1 + Congruence + Side + Colour + a1 | Child)
Code 
f1 = @formula(speed ~  1 + Congruence + Side + Colour + a1 + gender +
                      (1 + Congruence + a1 | Child));
m1 = fit(MixedModel, f1, dat; contrasts);

f2 = @formula(speed ~  1 + Congruence + Side + Colour + a1 + gender +
                      (1 + Congruence + Side + Colour + a1 | Child));
m2 = fit(MixedModel, f2, dat; contrasts);
Code 
lrt1 = gof_summary = let
  nms = [:m0, :m1, :m2]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m0, m1, m2)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrt1, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.2: LRT of m0, m1, and m2
name dof devianc AIC AICc BIC χ² χ²_dof p
m0 8 882.824 898.824 898.831 962.41 . . .
m1 13 500.56 526.56 526.577 629.887 382.264 5.0 0.0
m2 22 381.669 425.669 425.717 600.531 118.891 9.0 0.0

Including Congruence, Colour, Side, and age as variance components in the random effect structure of Child improves the models fit.

  • m3: speed ~ 1 + Congruence + Side + Colour + a1 + gender + zerocorr(1 + Congruence + Side + Colour + a1 | Child)
Code 
f3 = @formula(speed ~  1 + Congruence + Side + Colour + a1 + gender +
              zerocorr(1 + Congruence + Side + Colour + a1 | Child));
m3 = fit(MixedModel, f3, dat; contrasts);
Code 
lrt2 = gof_summary = let
  nms = [:m3, :m2]
  mods = eval.(nms)
   lrt = MixedModels.likelihoodratiotest(m3, m2)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrt2, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.3: LRT of m2 and m3
name dof devianc AIC AICc BIC χ² χ²_dof p
m3 12 408.198 432.198 432.212 527.577 . . .
m2 22 381.669 425.669 425.717 600.531 26.529 10.0 0.003

Keeping all correlation parameters from the random effect structure of Child has the better fit.

  • m4: speed ~ 1 + Congruence + Side + Colour + a1 + gender + (1 + a1 | Child) + zerocorr(0 + Congruence + Side + Colour | Child)
Code 
f4 = @formula(speed ~  1 + Congruence + Side + Colour + a1 + gender +
                       (1 + a1 | Child) +  
                       zerocorr(0 + Congruence + Side + Colour | Child));
m4 = fit(MixedModel, f4, dat; contrasts);
Code 
lrt3 = gof_summary = let
  nms = [:m3, :m4, :m2]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m3, m4, m2)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrt3, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.4: LRT of m3, m4, and m2
name dof devianc AIC AICc BIC χ² χ²_dof p
m3 12 408.198 432.198 432.212 527.577 . . .
m4 13 403.936 429.936 429.954 533.264 4.261 1.0 0.039
m2 22 381.669 425.669 425.717 600.531 22.267 9.0 0.008

Including the individual correlation parameters Congruence, Colour, Side, and age gives the best fit. Accordingly, m2 is reported in Section 6.3.2.

9.1.1.2 Means and standard deviations of model variabels

Code 
R"""
suppressWarnings(suppressMessages(library(flextable)))
dat_table <- suppressWarnings(suppressMessages(
rbind(
dat |>
  group_by(Time,Child) |>
  summarise(n=n(),
            age=mean(age),
            gender=unique(gender)) |>
  ungroup(Child) |>
  summarise(n=sum(n),
            age=paste0(round(mean(age),1),"±",round(sd(age),1)),
            gender=(paste0(sum(gender=="girl"),"/",sum(gender=="boy")))) |>
  t() |>
  janitor::row_to_names(row_number = 1) |>
  as.data.frame() |>
  rownames_to_column(var="Condition") |>
  mutate(Condition=case_when(Condition=="n"~"number of reaction times [n]",
                             Condition=="gender"~"gender [girls/boys]",
                             Condition=="age"~"age [years]")),
dat |>
  group_by(Time, Condition) |>
  summarize(m=paste0(round(mean(RT),0),"±",round(sd(RT),0))) |>
  pivot_wider(names_from = Time, values_from = c(m)) |>
  mutate(Condition=case_when(Condition=="blue_left_incon"~"blue left (incongruent)",
                             Condition=="blue_right_con"~"blue right (congruent)",
                             Condition=="red_left_con"~"red left (congruent)",
                             Condition=="red_right_incon"~"red right (incongruent)"))) 
))
"""
@rget dat_table;

pretty_table(dat_table, 
             header = (["Variable","t0","t1","t2","t3","t6","t8"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r ],
             backend = Val(:html))
Variable t0 t1 t2 t3 t6 t8
number of reaction times [n] 4388 4059 3904 3838 2586 2141
age [years] 9.2±0.5 9.4±0.5 9.6±0.5 10.1±0.5 10.7±0.5 11.6±0.5
gender [girls/boys] 35/41 30/39 28/40 28/37 18/26 16/21
blue left (incongruent) 1318±455 1206±427 1172±501 1090±434 1088±478 917±299
blue right (congruent) 1265±573 1200±575 1093±350 1047±425 1013±494 859±317
red left (congruent) 1231±505 1168±514 1092±514 1018±425 1049±499 897±398
red right (incongruent) 1303±481 1279±635 1204±679 1079±442 1119±524 896±290

9.1.2 Simon Task x EMOTIKON

9.1.2.1 Model selection

A minimal model m4_m0 was fitted according to the structure of m2. All EMOTIKON tests were included in m4_mall. Six additional models were fitted (m4_mRUN, m4_mStar, m4_mS20, m4_mSLJ, m4_mBPT, & m4_mOLB), with each dropping one EMOTIKON test.

  • m4_m0: speed ~ 1 + Congruence + Side + Colour + a1 + gender + (1 + Congruence + Side + Colour + a1 | Child)
  • m4_mRun: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Star_r + S20_r + SLJ + BPT + OLB_l + (1 + Congruence + Side + Colour + a1 | Child)
  • m4_mStar: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + S20_r + SLJ + BPT + OLB_l + (1 + Congruence + Side + Colour + a1 | Child)
  • m4_mS20: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + Star_r + SLJ + BPT + OLB_l + (1 + Congruence + Side + Colour + a1 | Child)
  • m4_mSLJ: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + Star_r + S20_r + BPT + OLB_l + (1 + Congruence + Side + Colour + a1 | Child)
  • m4_mBPT: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + Star_r + S20_r + SLJ + OLB_l + (1 + Congruence + Side + Colour + a1 | Child)
  • m4_mOLB: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + Star_r + S20_r + SLJ + BPT + (1 + Congruence + Side + Colour + a1 | Child)
  • m4_mall: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + Star_r + S20_r + SLJ + BPT + OLB_l + (1 + Congruence + Side + Colour + a1 | Child)
Code 
m4_m0 = fit(MixedModel, f2, Data; contrasts);

f4_mRun = @formula(speed ~ 1 + Congruence + Side + Colour + a1 + gender +
                                 Star_r + S20_r + SLJ + BPT + OLB_l +
                           (1 + a1 + Congruence + Side + Colour | Child));
m4_mRun = fit(MixedModel, f4_mRun, Data; contrasts);

f4_mStar = @formula(speed ~ 1 + Congruence + Side + Colour + a1 + gender +
                           Run +          S20_r + SLJ + BPT + OLB_l +
                           (1 + a1 + Congruence + Side + Colour | Child));
m4_mStar = fit(MixedModel, f4_mStar, Data; contrasts);

f4_mS20 = @formula(speed ~ 1 + Congruence + Side + Colour + a1 + gender +
                           Run + Star_r +         SLJ + BPT + OLB_l +
                           (1 + a1 + Congruence + Side + Colour | Child));
m4_mS20 = fit(MixedModel, f4_mS20, Data; contrasts);

f4_mSLJ = @formula(speed ~ 1 + Congruence + Side + Colour + a1 + gender +
                           Run + Star_r + S20_r +       BPT + OLB_l +
                           (1 + a1 + Congruence + Side + Colour | Child));
m4_mSLJ = fit(MixedModel, f4_mSLJ, Data; contrasts);

f4_mBPT = @formula(speed ~ 1 + Congruence + Side + Colour + a1 + gender +
                           Run + Star_r + S20_r + SLJ +       OLB_l +
                           (1 + a1 + Congruence + Side + Colour | Child));
m4_mBPT = fit(MixedModel, f4_mBPT, Data; contrasts);

f4_mOLB = @formula(speed ~ 1 + Congruence + Side + Colour + a1 + gender +
                           Run + Star_r + S20_r + SLJ + BPT +
                           (1 + a1 +  Congruence + Side + Colour | Child));
m4_mOLB = fit(MixedModel, f4_mOLB, Data; contrasts);

f4_mall = @formula(speed ~  1 + Congruence + Side + Colour + a1 + gender +
                           Run + Star_r + S20_r + SLJ + BPT + OLB_l +
                      (1 + Congruence + Side + Colour + a1 | Child));
m4_mall = fit(MixedModel, f4_mall, Data; contrasts);
Code 
lrtRun = gof_summary = let
  nms = [:m4_m0, :m4_mRun, :m4_mall]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m4_m0, m4_mRun, m4_mall)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
  end

pretty_table(lrtRun, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.5: LRT of m4_m0, m4_mRun, and m4_mall
name dof devianc AIC AICc BIC χ² χ²_dof p
m4_m0 22 410.105 454.105 454.155 628.218 . . .
m4_mRun 27 367.127 421.127 421.202 634.811 42.978 5.0 0.0
m4_mall 28 371.32 427.32 427.401 648.919 -4.193 1.0 NaN

LRT reveals that including Run improves the models fit.

Code 
lrtStar = gof_summary = let
  nms = [:m4_m0, :m4_mStar, :m4_mall]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m4_m0, m4_mStar, m4_mall)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrtStar, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.6: LRT of m4_m0, m4_mStar, and m4_mall
name dof devianc AIC AICc BIC χ² χ²_dof p
m4_m0 22 410.105 454.105 454.155 628.218 . . .
m4_mStar 27 380.002 434.002 434.077 647.686 30.103 5.0 0.0
m4_mall 28 371.32 427.32 427.401 648.919 8.682 1.0 0.003

LRT reveals that including Star improves the models fit.

Code 
lrtS20 = gof_summary = let
  nms = [:m4_m0, :m4_mS20, :m4_mall]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m4_m0, m4_mS20, m4_mall)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrtS20, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.7: LRT of m4_m0, m4_mS20, and m4_mall
name dof devianc AIC AICc BIC χ² χ²_dof p
m4_m0 22 410.105 454.105 454.155 628.218 . . .
m4_mS20 27 364.421 418.421 418.496 632.105 45.684 5.0 0.0
m4_mall 28 371.32 427.32 427.401 648.919 -6.899 1.0 NaN

LRT reveals that including S20_r does not improve the models fit.

Code 
lrtSLJ = gof_summary = let
  nms = [:m4_m0, :m4_mSLJ, :m4_mall]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m4_m0, m4_mSLJ, m4_mall)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrtSLJ, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.8: LRT of m4_m0, m4_mSLJ, and m4_mall
name dof devianc AIC AICc BIC χ² χ²_dof p
m4_m0 22 410.105 454.105 454.155 628.218 . . .
m4_mSLJ 27 361.089 415.089 415.164 628.773 49.016 5.0 0.0
m4_mall 28 371.32 427.32 427.401 648.919 -10.231 1.0 NaN

LRT reveals that including SLJ does not improve the models fit.

Code 
lrtBPT = gof_summary = let
  nms = [:m4_m0, :m4_mBPT, :m4_mall]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m4_m0, m4_mBPT, m4_mall)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrtBPT, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.9: LRT of m4_m0, m4_mBPT, and m4_mall
name dof devianc AIC AICc BIC χ² χ²_dof p
m4_m0 22 410.105 454.105 454.155 628.218 . . .
m4_mBPT 27 361.124 415.124 415.199 628.808 48.981 5.0 0.0
m4_mall 28 371.32 427.32 427.401 648.919 -10.196 1.0 NaN

LRT reveals that including BPT does not improve the models fit.

Code 
lrtOLB = gof_summary = let
  nms = [:m4_m0, :m4_mOLB, :m4_mall]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m4_m0, m4_mOLB, m4_mall)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrtOLB, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.10: LRT of m4_m0, m4_mOLB, and m4_mall
name dof devianc AIC AICc BIC χ² χ²_dof p
m4_m0 22 410.105 454.105 454.155 628.218 . . .
m4_mOLB 27 381.117 435.117 435.192 648.801 28.988 5.0 0.0
m4_mall 28 371.32 427.32 427.401 648.919 9.797 1.0 0.002

LRT reveals that including OLB improves the models fit.

Accordingly, a optimal model m5 using Run, Star, and OBL was fitted.

  • m5: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + Star_r + OLB_l + (1 + a1 + Congruence + Side + Colour | Child)
Code 
f5 = @formula(speed ~ 1 + Congruence + Side + Colour + a1  + gender +
                           Run + Star_r + OLB_l +
                           (1 + a1 + Congruence + Side + Colour | Child));
m5 = fit(MixedModel, f5, Data; contrasts);
Code 
lrt5 = gof_summary = let
  nms = [:m4_m0, :m5, :m4_mall]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m4_m0, m5, m4_mall)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrt5, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.11: LRT of m4_m0, m5, and m4_mall
name dof devianc AIC AICc BIC χ² χ²_dof p
m4_m0 22 410.105 454.105 454.155 628.218 . . .
m5 25 364.552 414.552 414.617 612.408 45.553 3.0 0.0
m4_mall 28 371.32 427.32 427.401 648.919 -6.768 3.0 NaN

Verification that dropping S20_r, SLJ, and BPT does not affect the models overall fit.

  • m6: speed ~ 1 + (Congruence + Side + Colour) * (Run + Star_r + OLB_l) + a1 + gender + (1 + a1 + Congruence + Side + Colour | Child)
Code 
f6 = @formula(speed ~ 1 + (Congruence + Side + Colour) * (Run + Star_r + OLB_l) +
                          a1  + gender +
                          (1 + a1 + Congruence + Side + Colour | Child));
m6 = fit(MixedModel, f6, Data; contrasts);
Code 
lrt6 = gof_summary = let
  nms = [:m5, :m6]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m5, m6)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrt6, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.12: LRT of m5 and m6
name dof devianc AIC AICc BIC χ² χ²_dof p
m5 25 364.552 414.552 414.617 612.408 . . .
m6 34 358.744 426.744 426.862 695.828 5.808 9.0 0.759

LRT shows no improvement of model fit for including interactions between Simon task parameters (i.e., Congruence, Colour, and Side) and EMOTIKON tests (Run, Star, and OLB).

Next, Run, Star and, OLB are included in the random effects without CP (i.e., m7) and with (i.e., m8).

  • m7: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + Star_r + OLB_l + (1 + a1 + Congruence + Side + Colour | Child) + zerocorr(0 + Run + Star_r + OLB_l | Child)
  • m8: speed ~ 1 + Congruence + Side + Colour + a1 + gender + Run + Star_r + OLB_l + (1 + a1 + Congruence + Side + Colour + Run + Star_r + OLB_l | Child)
Code 
f7 = @formula(speed ~ 1 + Congruence + Side + Colour + a1 + gender +
                           Run + Star_r + OLB_l + 
                           (1 + a1 + Congruence + Side + Colour | Child) +  
                           zerocorr(0 + Run + Star_r + OLB_l | Child));
m7 = fit(MixedModel, f7, Data; contrasts);

f8  = @formula(speed ~ 1 + Congruence + Side + Colour + a1 + gender +
                           Run + Star_r + OLB_l + 
                           (1 + a1 + Congruence + Side + Colour + Run + Star_r + OLB_l | Child));
m8 = fit(MixedModel, f8, Data; contrasts);
Code 
lrt7 = gof_summary = let
  nms = [:m5, :m7, :m8]
  mods = eval.(nms)
  lrt = MixedModels.likelihoodratiotest(m5, m7, m8)
  DataFrame(;
    name = nms, 
    dof=dof.(mods),
    deviance=deviance.(mods),
    AIC=aic.(mods),
    AICc=aicc.(mods),
    BIC=bic.(mods),
    χ²=vcat(:., round.(lrt.tests.deviancediff, digits=3)),
    χ²_dof=vcat(:., round.(lrt.tests.dofdiff, digits=3)),
    pvalue=vcat(:., round.(lrt.tests.pvalues, digits=3)),
  )
end

pretty_table(lrt7, 
             header = (["name", "dof", "devianc", "AIC", "AICc", "BIC", "χ²", "χ²_dof", "p"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.13: LRT of m5, m7, and m8
name dof devianc AIC AICc BIC χ² χ²_dof p
m5 25 364.552 414.552 414.617 612.408 . . .
m7 28 -350.44 -294.44 -294.36 -72.8415 714.992 3.0 0.0
m8 46 -414.837 -322.837 -322.623 41.2173 64.397 18.0 0.0

Including Run, Star, and OLB as variance components with associated correlation parameters in the random effects structure has the best fit.

Accordingly, m8 was reported in Section 6.4.2.

9.1.2.2 Fixed and random effects of m5

Code 
coefm5 = transform!(DataFrame(coeftable(m5)), 
              [:z] .=> (x -> round.(x, digits=1)) .=> [:z],
              [:"Coef." :"Std. Error"] .=> (x -> round.(x, digits=3)) .=> [:Est :SE],
              [:"Pr(>|z|)"] .=> (x -> ifelse.(x .< 0.001, "<0.001", round.(x, digits=3))) .=> [:p]);

pretty_table(select(coefm5, :Name, :Est, :SE, :z, :p), 
             header = (["Name","Est.","SE","z","p"]),
             alignment=[:l, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.14: Fixed effects of the Simon Task main effects, 6 min run, star run, one leg balance, age, and gender
Name Est. SE z p
(Intercept) 0.987 0.018 54.5 <0.001
Congruence: con 0.029 0.003 11.2 <0.001
Side: right -0.004 0.002 -1.8 0.076
Colour: red 0.005 0.003 1.6 0.102
a1 0.125 0.009 13.5 <0.001
gender: boy 0.045 0.016 2.9 0.004
Run -0.011 0.005 -2.4 0.018
Star_r 0.023 0.005 4.7 <0.001
OLB_l 0.013 0.003 4.7 <0.001
Code 
vcm5 = VarCorr(m5);
c = round.(vcm5.σρ.Child.ρ, digits = 3);
v = round.(values(vcm5.σρ.Child.σ), digits = 3);
sd = round.(abs2.(values(vcm5.σρ.Child.σ)), digits = 3);
n = keys(vcm5.σρ.Child.σ);
vcm5_tbl =DataFrame(Groups=["Child"," "," "," "," ","Residual"],
          Name=["Grand mean",n[2],n[3],n[4],n[5]," "],
          Var=[v[1],v[2],v[3],v[4],v[5],round.(vcm5.s, digits=3)],
          SD=[sd[1],sd[2],sd[3],sd[4],sd[5],round.(abs2.(vcm5.s), digits=3)],
          c1=[" ", c[1], c[2], c[3], c[4], " "],
          c2=[" ", " ",  c[5], c[6], c[7], " "],
          c3=[" ", " ",  " ",  c[8], c[9], " "],
          c4=[" ", " ",  " ",  " ",  c[10]," "]);

pretty_table(vcm5_tbl, 
             header = (["Groups","Name","Var","SD","Corr"," "," "," ",]),
             alignment=[:l, :l, :r, :r, :r, :r, :r, :r],
             backend = Val(:html))
Table 9.15: Random effects of the Simon Task main effects, age, and gender
Groups Name Var SD Corr
Child Grand mean 0.148 0.022
a1 0.07 0.005 0.361
Congruence: con 0.016 0.0 0.066 0.068
Side: right 0.012 0.0 -0.09 0.22 -0.056
Colour: red 0.024 0.001 -0.43 0.113 -0.431 -0.617
Residual 0.24 0.058

9.1.2.3 Means and standard deviations of model variabels

Code 
@rput Data;
R"""
load(file="data/emo.rda")
Data_table <- suppressWarnings(suppressMessages(
rbind(
Data |>
  group_by(Time,Child) |>
  summarise(n=n(),
            age=mean(age),
            gender=unique(gender)) |>
  ungroup(Child) |>
  left_join(emo |>
              select(Child, Time, Test, Score) |>
              pivot_wider(names_from = Test, values_from = Score), multiple = "all") |>
  summarise(n=sum(n),
            gender=(paste0(sum(gender=="girl"),"/",sum(gender=="boy"))),
            across(c(age,S20_r,Star_r,BPT,OLB_l), function(x) paste0(round(mean(x),1), "±", round(sd(x),1))),
            across(c(Run,SLJ), function(x) paste0(round(mean(x),0), "±", round(sd(x),0)))) |>
  t() |>
  janitor::row_to_names(row_number = 1) |>
  as.data.frame() |>
  rownames_to_column(var="Condition") |>
  mutate(Condition=case_when(Condition=="n"~"number of reaction times [n]",
                             Condition=="gender"~"gender [girls/boys]",
                             Condition=="age"~"age [years]",
                             Condition=="Run"~"6 min run",
                             Condition=="S20_r"~"20 m sprint",
                             Condition=="SLJ"~"standing long jump",
                             Condition=="Star_r"~"star run",
                             Condition=="BPT"~"ball push test",
                             Condition=="OLB_l"~"one leg balance")),
Data |>
  group_by(Time, Condition) |>
  summarize(m=paste0(round(mean(RT),0),"±",round(sd(RT),0))) |>
  pivot_wider(names_from = Time, values_from = c(m)) |>
  mutate(Condition=case_when(Condition=="blue_left_incon"~"blue left (incongruent)",
                             Condition=="blue_right_con"~"blue right (congruent)",
                             Condition=="red_left_con"~"red left (congruent)",
                             Condition=="red_right_incon"~"red right (incongruent)")))
))
"""

@rget Data_table;

pretty_table(Data_table,              
             header = (["Variable","t0","t1","t2","t3","t6","t8"]),
             alignment=[:l, :r, :r, :r, :r, :r, :r ],
             backend = Val(:html))
Table 9.16: Means and standard deviations of Simon task conditions, EMOTIKON tests, and age
Variable t0 t1 t2 t3 t6 t8
number of reaction times [n] 4388 3940 3724 3611 2586 1967
gender [girls/boys] 35/41 29/38 28/37 28/33 18/26 15/19
age [years] 9.2±0.5 9.4±0.5 9.6±0.5 10.1±0.6 10.7±0.5 11.6±0.5
20 m sprint 4.1±0.4 4.1±0.4 4.1±0.4 4.1±0.4 4.2±0.4 4.3±0.4
star run 1.8±0.2 1.8±0.2 1.8±0.2 1.9±0.2 2±0.2 2±0.3
ball push test 3.6±0.7 3.7±0.7 3.7±0.6 3.8±0.7 4.1±0.7 4.8±0.8
one leg balance 2±0.8 2.2±0.7 2.1±0.8 2.1±0.8 2.4±0.6 2.5±0.7
6 min run 859±125 879±135 846±139 868±136 864±146 844±177
standing long jump 108±19 108±20 110±20 112±21 114±21 124±21
blue left (incongruent) 1318±455 1206±426 1177±510 1094±442 1088±478 922±307
blue right (congruent) 1265±573 1202±580 1095±355 1051±434 1013±494 861±319
red left (congruent) 1231±505 1170±516 1095±524 1016±395 1049±499 902±392
red right (incongruent) 1303±481 1274±627 1210±693 1083±449 1119±524 900±293