We take the kwdyz11 dataset (Kliegl et al., 2011) from an experiment looking at three effects of visual cueing under four different cue-target relations (CTRs). Two horizontal rectangles are displayed above and below a central fixation point or they displayed in vertical orientation to the left and right of the fixation point. Subjects react to the onset of a small visual target occurring at one of the four ends of the two rectangles. The target is cued validly on 70% of trials by a brief flash of the corner of the rectangle at which it appears; it is cued invalidly at the three other locations 10% of the trials each.
We specify three contrasts for the four-level factor CTR that are derived from spatial, object-based, and attractor-like features of attention. They map onto sequential differences between appropriately ordered factor levels. At the level of fixed effects, there is the noteworthy result, that the attraction effect was estimated at 2 ms, that is clearly not significant. Nevertheless, there was a highly reliable variance component (VC) estimated for this effect. Moreover, the reliable individual differences in the attraction effect were negatively correlated with those in the spatial effect.
Unfortunately, a few years after the publication, we determined that the reported LMM is actually singular and that the singularity is linked to a theoretically critical correlation parameter (CP) between the spatial effect and the attraction effect. Fortunately, there is also a larger dataset kkl15 from a replication and extension of this study (Kliegl et al., 2015), analyzed with kkl15.jl notebook. The critical CP (along with other fixed effects and CPs) was replicated in this study.
A more comprehensive analysis was reported in the parsimonious mixed-model paper (Bates et al., 2015). Data and R scripts are also available in R-package RePsychLing. In this and the complementary kkl15.jl scripts, we provide some corresponding analyses with MixedModels.jl.
using AlgebraOfGraphics
using AlgebraOfGraphics: density
using CairoMakie
using CategoricalArrays
using Chain
using DataFrames
using DataFrameMacros
using Distributions
using MixedModels
using MixedModelsMakie
using Random
using SMLP2023: dataset
using StatsBase
CairoMakie.activate!(; type="svg")
import ProgressMeter
ProgressMeter.ijulia_behavior(:clear)dat = DataFrame(dataset(:kwdyz11))
describe(dat)| Row | variable | mean | min | median | max | nmissing | eltype |
|---|---|---|---|---|---|---|---|
| Symbol | Union… | Any | Union… | Any | Int64 | DataType | |
| 1 | Item | I002 | I603 | 0 | String | ||
| 2 | Subj | S01 | S61 | 0 | String | ||
| 3 | CTR | dod | val | 0 | String | ||
| 4 | dir | hor | ver | 0 | String | ||
| 5 | rt | 370.426 | 150.1 | 358.6 | 705.7 | 0 | Float32 |
We recommend to code the levels/units of random factor / grouping variable not as a number, but as a string starting with a letter and of the same length for all levels/units.
We also recommend to sort levels of factors into a meaningful order, that is overwrite the default alphabetic ordering. This is also a good place to choose alternative names for variables in the context of the present analysis.
The LMM analysis is based on log-transformed reaction times lrt, indicated by a boxcox() check of model residuals. With the exception of diagnostic plots of model residuals, the analysis of untransformed reaction times did not lead to different results and exhibited the same problems of model identification (see Kliegl et al., 2011).
Comparative density plots of all response times by cue-target relation, Figure 1, show the times for valid cues to be faster than for the other conditions.
draw(
data(dat) *
mapping(
:rt => log => "log(Reaction time [ms])";
color=:CTR =>
renamer("val" => "valid cue", "sod" => "some obj/diff pos", "dos" => "diff obj/same pos", "dod" => "diff obj/diff pos") => "Cue-target relation",
) *
density(),
)
An alternative visualization without overlap of the conditions can be accomplished with ridge plots.
To be done
For the next set of plots we average subjects’ data within the four experimental conditions. This table could be used as input for a repeated-measures ANOVA.
dat_subj = combine(
groupby(dat, [:Subj, :CTR]),
:rt => length => :n,
:rt => mean => :rt_m,
:rt => (r -> mean(log, r)) => :lrt_m,
)
dat_subj.CTR = categorical(dat_subj.CTR, levels=levels(dat.CTR))
dat_subj| Row | Subj | CTR | n | rt_m | lrt_m |
|---|---|---|---|---|---|
| String | Cat… | Int64 | Float32 | Float32 | |
| 1 | S01 | val | 330 | 413.332 | 6.01272 |
| 2 | S01 | sod | 48 | 437.721 | 6.0682 |
| 3 | S01 | dos | 47 | 443.366 | 6.08269 |
| 4 | S01 | dod | 45 | 434.316 | 6.06079 |
| 5 | S02 | val | 333 | 365.899 | 5.87511 |
| 6 | S02 | sod | 47 | 396.149 | 5.95916 |
| 7 | S02 | dos | 46 | 439.487 | 6.06949 |
| 8 | S02 | dod | 48 | 441.042 | 6.06883 |
| 9 | S03 | val | 336 | 371.446 | 5.90489 |
| 10 | S03 | sod | 46 | 446.854 | 6.09464 |
| 11 | S03 | dos | 48 | 471.302 | 6.14287 |
| 12 | S03 | dod | 47 | 476.532 | 6.15706 |
| 13 | S04 | val | 336 | 403.181 | 5.99124 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 233 | S59 | val | 330 | 313.032 | 5.72433 |
| 234 | S59 | sod | 46 | 304.852 | 5.70367 |
| 235 | S59 | dos | 48 | 342.612 | 5.80927 |
| 236 | S59 | dod | 47 | 308.191 | 5.71206 |
| 237 | S60 | val | 330 | 372.835 | 5.90832 |
| 238 | S60 | sod | 48 | 382.087 | 5.93183 |
| 239 | S60 | dos | 44 | 395.109 | 5.96918 |
| 240 | S60 | dod | 47 | 387.689 | 5.95208 |
| 241 | S61 | val | 320 | 312.374 | 5.72542 |
| 242 | S61 | sod | 45 | 383.002 | 5.92798 |
| 243 | S61 | dos | 45 | 357.998 | 5.85563 |
| 244 | S61 | dod | 46 | 391.087 | 5.95344 |
boxplot(
dat_subj.CTR.refs,
dat_subj.lrt_m;
orientation=:horizontal,
show_notch=true,
axis=(;
yticks=(
1:4,
[
"valid cue",
"same obj/diff pos",
"diff obj/same pos",
"diff obj/diff pos",
],
),
),
figure=(; resolution=(800, 300)),
)
Mean of log reaction times for four cue-target relations. Targets appeared at (a) the cued position (valid) in a rectangle, (b) in the same rectangle cue, but at its other end, (c) on the second rectangle, but at a corresponding horizontal/vertical physical distance, or (d) at the other end of the second rectangle, that is \(\sqrt{2}\) of horizontal/vertical distance diagonally across from the cue, that is also at larger physical distance compared to (c).
A better alternative to the boxplot is a dotplot. It also displays subjects’ condition means.
To be done
contrasts = Dict(
:CTR => SeqDiffCoding(; levels=["val", "sod", "dos", "dod"]),
:Subj => Grouping(),
)
m1 = let
form = @formula(log(rt) ~ 1 + CTR + (1 + CTR | Subj))
fit(MixedModel, form, dat; contrasts)
endMinimizing 240 Time: 0:00:00 ( 1.40 ms/it)
objective: -12782.373697791405| Est. | SE | z | p | σ_Subj | |
|---|---|---|---|---|---|
| (Intercept) | 5.9358 | 0.0185 | 320.54 | <1e-99 | 0.1441 |
| CTR: sod | 0.0878 | 0.0084 | 10.48 | <1e-24 | 0.0582 |
| CTR: dos | 0.0366 | 0.0062 | 5.92 | <1e-08 | 0.0274 |
| CTR: dod | -0.0086 | 0.0060 | -1.43 | 0.1515 | 0.0249 |
| Residual | 0.1920 |
VarCorr(m1)| Column | Variance | Std.Dev | Corr. | |||
|---|---|---|---|---|---|---|
| Subj | (Intercept) | 0.0207636 | 0.1440957 | |||
| CTR: sod | 0.0033848 | 0.0581791 | +0.48 | |||
| CTR: dos | 0.0007528 | 0.0274366 | -0.24 | -0.15 | ||
| CTR: dod | 0.0006220 | 0.0249400 | +0.30 | +0.93 | -0.43 | |
| Residual | 0.0368544 | 0.1919749 |
issingular(m1)trueLMM m1 is not fully supported by the data; it is overparameterized. This is also visible in the PCA: only three, not four PCS are needed to account for all the variance and covariance in the random-effect structure. The problem is the +.93 CP for spatial sod and attraction dod effects.
first(MixedModels.PCA(m1))
Principal components based on correlation matrix
(Intercept) 1.0 . . .
CTR: sod 0.48 1.0 . .
CTR: dos -0.24 -0.15 1.0 .
CTR: dod 0.3 0.93 -0.43 1.0
Normalized cumulative variances:
[0.5887, 0.8095, 1.0, 1.0]
Component loadings
PC1 PC2 PC3 PC4
(Intercept) -0.4 0.04 0.9 0.17
CTR: sod -0.6 0.4 -0.16 -0.68
CTR: dos 0.33 0.91 0.06 0.23
CTR: dod -0.61 0.08 -0.41 0.68Do model residuals meet LMM assumptions? Classic plots are
The slant in residuals show a lower and upper boundary of reaction times, that is we have have too few short and too few long residuals. Not ideal, but at least width of the residual band looks similar across the fitted values, that is there is no evidence for heteroskedasticity.
CairoMakie.activate!(; type="png")
set_aog_theme!()
draw(
data((; f=fitted(m1), r=residuals(m1))) *
mapping(:f => "Fitted values", :r => "Residual from model m1") *
visual(Scatter);
)
With many observations the scatterplot is not that informative. Contour plots or heatmaps may be an alternative.
CairoMakie.activate!(; type="png")
draw(
data((; f=fitted(m1), r=residuals(m1))) *
mapping(
:f => "Fitted log response time", :r => "Residual from model m1"
) *
density();
)
The plot of quantiles of model residuals over corresponding quantiles of the normal distribution should yield a straight line along the main diagonal.
qqnorm(residuals(m1); qqline=:none)
The violation of expectation is again due to the fact that the distribution of residuals is much narrower than expected from a normal distribution, as shown in Figure 5. Overall, it does not look too bad.
let
n = nrow(dat)
dat_rz = DataFrame(;
value=vcat(residuals(m1) ./ std(residuals(m1)), randn(n)),
curve=vcat(fill.("residual", n), fill.("normal", n)),
)
draw(
data(dat_rz) *
mapping(:value => "Standardized residuals"; color=:curve) *
density(; bandwidth=0.1);
)
end
Now we move on to visualizations that are based on model parameters and subjects’ data, that is “predictions” of the LMM for subject’s GM and experimental effects. Three important plots are
The first plot overlays shrinkage-corrected conditional modes of the random effects with within-subject-based and pooled GMs and experimental effects.
To be done
The caterpillar plot, Figure 6, also reveals the high correlation between spatial sod and attraction dod effects.
caterpillar!(
Figure(; resolution=(800, 1000)), ranefinfo(m1, :Subj); orderby=2
)
Figure 7 provides more evidence for a problem with the visualization of the spatial sod and attraction dod CP. The corresponding panel illustrates an implosion of conditional modes.
shrinkageplot!(Figure(; resolution=(1000, 1000)), m1)
Here we
We generate 2500 samples for the 15 model parameters (4 fixed effect, 4 VCs, 6 CPs, and 1 residual).
Random.seed!(1234321)
samp = parametricbootstrap(2500, m1)
tbl = samp.tblTable with 26 columns and 2500 rows:
obj β1 β2 β3 β4 σ ⋯
┌───────────────────────────────────────────────────────────────────
1 │ -12802.3 5.93392 0.0864217 0.048887 -0.0121824 0.192002 ⋯
2 │ -12794.3 5.95776 0.0989839 0.0318752 -0.00666179 0.191977 ⋯
3 │ -12930.6 5.93669 0.088561 0.0209803 -0.000658234 0.191496 ⋯
4 │ -12946.5 5.93165 0.0767128 0.0343559 -0.00499475 0.191438 ⋯
5 │ -12677.5 5.91938 0.0984591 0.0272961 -0.00555834 0.19247 ⋯
6 │ -12869.7 5.96313 0.104457 0.0321629 -0.00719225 0.191751 ⋯
7 │ -12478.1 5.93281 0.0947864 0.0409705 -0.0046158 0.193076 ⋯
8 │ -12463.7 5.94439 0.0977164 0.0321099 -0.0062546 0.193115 ⋯
9 │ -12836.2 5.9728 0.0971351 0.0416612 6.2017e-5 0.191843 ⋯
10 │ -12932.3 5.94709 0.0742454 0.0505263 -0.0131451 0.191583 ⋯
11 │ -12736.2 5.93625 0.0728173 0.047942 -0.0183152 0.192193 ⋯
12 │ -13048.0 5.89421 0.0680602 0.0481648 -0.0214149 0.191101 ⋯
13 │ -13313.8 5.92608 0.0892317 0.047389 -0.0141679 0.190142 ⋯
14 │ -12575.0 5.89235 0.0846168 0.0341529 -0.00613056 0.192728 ⋯
15 │ -12833.1 5.94217 0.0774876 0.0422338 -0.00278894 0.191772 ⋯
16 │ -12504.5 5.94385 0.0898018 0.0429783 -0.00705609 0.192833 ⋯
17 │ -12875.5 5.96624 0.0976098 0.0370804 -0.000718399 0.191577 ⋯
18 │ -13287.4 5.97492 0.101651 0.0407503 -0.0163215 0.190317 ⋯
19 │ -13088.2 5.95426 0.101522 0.0342398 0.0022965 0.190991 ⋯
20 │ -12820.6 5.92995 0.0775968 0.0394342 -0.0112547 0.191861 ⋯
21 │ -12634.4 5.92989 0.0836736 0.035749 -0.00206045 0.192456 ⋯
22 │ -13097.6 5.92438 0.088869 0.0332515 -0.00701304 0.191018 ⋯
23 │ -12837.4 5.90782 0.0861051 0.0385006 -0.00752588 0.191891 ⋯
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱The upper limit of the interval for the critical CP CTR: sod, CTR: dod is hitting the upper wall of a perfect correlation. This is evidence of singularity. The other intervals do not exhibit such pathologies; they appear to be ok.
confint(samp)DictTable with 2 columns and 15 rows:
par lower upper
────┬───────────────────────
β1 │ 5.89917 5.97256
β2 │ 0.0719305 0.104588
β3 │ 0.0251181 0.0491664
β4 │ -0.0207184 0.00268322
ρ1 │ 0.243806 0.712967
ρ2 │ -0.926998 0.231027
ρ3 │ -0.687501 0.536383
ρ4 │ -0.132528 0.739327
ρ5 │ 0.573732 0.999995
ρ6 │ -0.896244 0.44397
σ │ 0.19051 0.19359
σ1 │ 0.116558 0.168559
σ2 │ 0.0454979 0.0708042
σ3 │ 0.00939896 0.041265
σ4 │ 0.0133637 0.0375904draw(
data(tbl) *
mapping(:σ => "Residual standard deviation") *
density();
)
The shortest coverage interval for the GM ranges from 376 to 404 ms. To keep the plot range small we do not include its density here.
labels = [
"CTR: sod" => "spatial effect",
"CTR: dos" => "object effect",
"CTR: dod" => "attraction effect",
"(Intercept)" => "grand mean",
]
draw(
data(tbl) *
mapping(
[:β2, :β3, :β4] .=> "Experimental effect size [ms]";
color=dims(1) => renamer(["spatial", "object", "attraction"] .* " effect") =>
"Experimental effects",
) *
density();
)
The densitiies correspond nicely with the shortest coverage intervals.
draw(
data(tbl) *
mapping(
[:σ1, :σ2, :σ3, :σ4] .=> "Standard deviations [ms]";
color=dims(1) =>
renamer(append!(["Grand mean"],["spatial", "object", "attraction"] .* " effect")) =>
"Variance components",
) *
density();
)
The VC are all very nicely defined.
let
labels = [
"(Intercept), CTR: sod" => "GM, spatial",
"(Intercept), CTR: dos" => "GM, object",
"CTR: sod, CTR: dos" => "spatial, object",
"(Intercept), CTR: dod" => "GM, attraction",
"CTR: sod, CTR: dod" => "spatial, attraction",
"CTR: dos, CTR: dod" => "object, attraction",
]
draw(
data(tbl) *
mapping(
[:ρ1, :ρ2, :ρ3, :ρ4, :ρ5, :ρ6] .=> "Correlation";
color=dims(1) => renamer(last.(labels)) => "Correlation parameters",
) *
density();
)
end
Two of the CPs stand out positively. First, the correlation between GM and the spatial effect is well defined. Second, as discussed throughout this script, the CP between spatial and attraction effect is close to the 1.0 border and clearly not well defined. Therefore, this CP will be replicated with a larger sample in script kkl15.jl (Kliegl et al., 2015).