Eye–Voice Span in Rapid Automatized Naming: Pan et al. (2013)

Authors

Reinhold Kliegl

Jinger Pan

Published

2026-08-28

After working through this page you will be able to:

reproduce a covariate-based LMM analysis (gaze duration and eye–voice span as predictors of psychometric RAN) from a published study;
complexify a fixed-effect structure with continuous-by-categorical interactions and prune it with likelihood-ratio tests;
visualize higher-order interactions three ways — with partial effects (MixedModelsExtras.partial_fitted, the Julia analogue of R’s remef()), with marginal effects (Effects.jl), and with model diagnostics (MixedModelsMakie.jl).

Before you start

Prerequisites: Analysis of the sleepstudy data for the modeling vocabulary; Transformations and Effects for Effects.jl; Creating multi-panel plots for the plotting stack.

Data used: pylsk13.rda, distributed alongside this page (examples/eyevoicespan/). The original R script (pylsk13_v4.R) and the paper (Pan_etal.DevSci.2013.pdf) are in the same folder.

1 Background

Pan et al. (2013) measured the eye movements of Chinese dyslexic and control children during rapid automatized naming (RAN) of digits (alphanumeric) and dice surfaces (symbolic). Both stimulus types require identical oral responses, which controls for effects associated with speech production.

The question addressed on this page is a psychometric one: how well do two eye-movement measures obtained from the computerized RAN task — gaze duration (gaze, the summed fixation time on an item) and the eye–voice span (evs, the distance in characters between the currently fixated item and the item currently being named) — predict naming speed (ran, seconds to name a card of items), and does that prediction differ between the two Conditions (digit/dice) and the two Groups (control/dyslexic)?

The analysis follows pylsk13_v4.R:

Setup — read the data, check whether the response needs a transformation, center the covariates.
Model building — start from an additive covariate model, add the Condition × Group interaction, then the full Condition × Group × evs × gaze factorial, and prune it back with likelihood-ratio tests.
Visualization — partial-effect plots (Figures 1 and 2 of the paper), marginal-effect plots, and residual diagnostics.

2 Packages

Code

using AlgebraOfGraphics
using AlgebraOfGraphics: density, linear
using BoxCox
using CairoMakie
using CategoricalArrays
using DataFrames
using Effects
using MixedModels
using MixedModelsMakie
using MixedModelsExtras          # partial_fitted(): the Julia analogue of remef()
using RData                      # to read the .rda file
using Statistics
using StatsBase

using SMLP2026: fit_or_restore

const progress = false
set_aog_theme!()

3 Read and prepare the data

pylsk13.rda contains a single data frame dat with one row per child × condition (56 children × 2 conditions = 112 rows).

dat = DataFrame(load(joinpath(@__DIR__, "pylsk13.rda"))["dat"])

112×8 DataFrame

87 rows omitted

Row	subj	group	condition	ran	gaze	evs	accuracy	fluency
	Int32	Int32	Int32	Float64	Int32	Float64	Int32	Int32
1	101	1	1	13.79	242	1.33	130	1192
2	101	1	2	19.61	397	0.97	130	1192
3	102	1	1	13.28	318	0.84	119	1696
4	102	1	2	17.43	372	1.06	119	1696
5	103	1	1	22.66	403	0.83	132	1233
6	103	1	2	27.47	487	0.8	132	1233
7	104	1	1	13.85	295	1.28	127	1124
8	104	1	2	25.05	442	0.95	127	1124
9	105	1	1	12.75	281	1.15	103	922
10	105	1	2	26.17	412	0.89	103	922
11	107	1	1	16.36	303	1.13	124	1140
12	107	1	2	28.14	441	0.87	124	1140
13	108	1	1	15.33	323	1.25	139	1005
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
101	231	2	1	16.21	371	1.05	83	317
102	231	2	2	30.27	511	0.75	83	317
103	232	2	1	24.96	312	1.15	75	344
104	232	2	2	28.24	411	1.05	75	344
105	233	2	1	22.55	420	0.18	92	583
106	233	2	2	24.79	528	0.67	92	583
107	234	2	1	25.52	435	0.73	90	610
108	234	2	2	29.59	557	0.77	90	610
109	235	2	1	17.39	382	0.84	93	527
110	235	2	2	28.1	520	0.68	93	527
111	236	2	1	19.48	326	0.97	94	752
112	236	2	2	28.89	495	0.54	94	752

We follow the general recommendation to code the levels of grouping variables as strings (here S101, S102, …), to give factors meaningful level names, and to fix a sensible level order. ran is log-transformed (justified below) and the two covariates are centered (not scaled) so that the intercept and the lower-order terms are interpretable at the average gaze and evs.

transform!(dat,
  :subj => (x -> categorical(string.("S", x))) => :Subj,
  :group => (x -> categorical(recode(x, 1 => "control", 2 => "dyslexic"))) => :Group,
  :condition => (x -> categorical(recode(x, 1 => "digit", 2 => "dice"))) => :Condition,
  :ran => ByRow(log) => :lran,
)
levels!(dat.Group, ["control", "dyslexic"])
levels!(dat.Condition, ["digit", "dice"])
dat.evs_c = dat.evs .- mean(dat.evs)
dat.gaze_c = dat.gaze .- mean(dat.gaze)
describe(dat)

14×7 DataFrame

Row	variable	mean	min	median	max	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Int64	DataType
1	subj	171.696	101	202.5	236	0	Int32
2	group	1.53571	1	2.0	2	0	Int32
3	condition	1.5	1	1.5	2	0	Int32
4	ran	21.9127	10.06	22.58	42.42	0	Float64
5	gaze	402.429	242	406.0	557	0	Int32
6	evs	0.871607	0.18	0.84	1.73	0	Float64
7	accuracy	104.554	55	94.0	141	0	Int32
8	fluency	782.429	253	730.0	1696	0	Int32
9	Subj		S101		S236	0	CategoricalValue{String, UInt32}
10	Group		control		dyslexic	0	CategoricalValue{String, UInt32}
11	Condition		digit		dice	0	CategoricalValue{String, UInt32}
12	lran	3.04789	2.30857	3.11706	3.74762	0	Float64
13	evs_c	-2.83503e-16	-0.691607	-0.0316071	0.858393	0	Float64
14	gaze_c	-1.82711e-14	-160.429	3.57143	154.571	0	Float64

3.1 Does the response need a transformation?

The R script uses MASS::boxcox() to justify a log transform of ran (but not of gaze or evs). The Julia equivalent fits a BoxCoxTransformation to a mixed model of the untransformed response.

m_ran = fit(MixedModel,
  @formula(ran ~ 1 + Condition * Group + evs_c + gaze_c + (1 | Subj)),
  dat;
  contrasts=Dict(:Condition => DummyCoding(; base="digit"),
                 :Group => DummyCoding(; base="control")),
  progress,
)
bc = fit(BoxCoxTransformation, m_ran; progress)

Box-Cox transformation

estimated λ: 0.3804
resultant transformation:

 y^0.4 - 1
-----------
    0.4

Code

boxcoxplot(bc; conf_level=0.95)

Figure 1: Box–Cox profile for the RAN response. The optimum is close to λ = 0, i.e. a log transform.

The profile peaks near λ = 0, so log(ran) (= lran) is a reasonable choice and matches the published analysis.

3.2 Condition and Group means

cell_means = combine(
  groupby(dat, [:Group, :Condition]),
  :lran => mean => :lran_m,
  :lran => (x -> std(x) / sqrt(length(x))) => :lran_se,
  :evs => mean => :evs_m,
  :gaze => mean => :gaze_m,
  nrow => :n,
)
cell_means

4×7 DataFrame

Row	Group	Condition	lran_m	lran_se	evs_m	gaze_m	n
	Cat…	Cat…	Float64	Float64	Float64	Float64	Int64
1	control	digit	2.66863	0.0375581	1.09731	317.346	26
2	control	dice	3.15163	0.0304036	0.876923	426.038	26
3	dyslexic	digit	3.02093	0.0317818	0.810333	385.9	30
4	dyslexic	dice	3.31362	0.0266264	0.732667	472.233	30

Code

draw(
  data(cell_means) *
  mapping(:Condition, :lran_m => "log(RAN)"; color=:Group, group=:Group) *
  visual(ScatterLines; markersize=16),
)

Figure 2: Mean log(RAN) by Condition and Group. Naming dice is slower than naming digits in both groups, and the group difference is larger for digits.

Code

draw(
  data(dat) *
  mapping(:lran => "log(RAN)"; color=:Condition) *
  density(),
)

Figure 3: Comparative density of log(RAN) by Condition.

4 Linear mixed models

The two categorical predictors use treatment (dummy) contrasts with the same reference levels as the published analysis (digit and control), so the coefficients reproduce Table 3 of Pan et al. (2013).

Tip

If the intercept should estimate the grand mean rather than the digit/control cell — usually preferable in factorial designs — swap DummyCoding for EffectsCoding. The omnibus tests below are unchanged; only the meaning of the lower-order coefficients changes.

contrasts = Dict(
  :Condition => DummyCoding(; base="digit"),
  :Group => DummyCoding(; base="control"),
)

Dict{Symbol, DummyCoding} with 2 entries:
  :Group     => DummyCoding("control", nothing)
  :Condition => DummyCoding("digit", nothing)

All models keep the maximal grouping structure supported by this design: a single by-child random intercept ((1 | Subj)). With only two observations per child there is no room for by-child random slopes.

4.1 Model building

m00 = fit_or_restore("eyevoicespan_m00.json", MixedModel,
  @formula(lran ~ 1 + Condition + Group + evs_c + gaze_c + (1 | Subj)),
  dat; contrasts, progress)

	Est.	SE	z	p	σ_Subj
(Intercept)	2.8926	0.0354	81.80	<1e-99	0.0857
Condition: dice	0.1861	0.0378	4.92	<1e-06
Group: dyslexic	0.1163	0.0362	3.22	0.0013
evs_c	-0.1934	0.0713	-2.71	0.0067
gaze_c	0.0017	0.0004	4.69	<1e-05
Residual	0.1091

m05 = fit_or_restore("eyevoicespan_m05.json", MixedModel,
  @formula(lran ~ 1 + Condition * Group + evs_c + gaze_c + (1 | Subj)),
  dat; contrasts, progress)

	Est.	SE	z	p	σ_Subj
(Intercept)	2.8436	0.0374	76.04	<1e-99	0.0903
Condition: dice	0.2710	0.0438	6.19	<1e-09
Group: dyslexic	0.1936	0.0427	4.53	<1e-05
evs_c	-0.1688	0.0696	-2.43	0.0153
gaze_c	0.0016	0.0004	4.54	<1e-05
Condition: dice & Group: dyslexic	-0.1303	0.0390	-3.34	0.0008
Residual	0.1000

m10 = fit_or_restore("eyevoicespan_m10.json", MixedModel,
  @formula(lran ~ 1 + Condition * Group * evs_c * gaze_c + (1 | Subj)),
  dat; contrasts, progress)

	Est.	SE	z	p	σ_Subj
(Intercept)	3.0252	0.0788	38.39	<1e-99	0.0948
Condition: dice	0.0758	0.0867	0.87	0.3817
Group: dyslexic	-0.0243	0.0848	-0.29	0.7747
evs_c	-0.8446	0.2443	-3.46	0.0005
gaze_c	0.0040	0.0011	3.73	0.0002
Condition: dice & Group: dyslexic	0.1388	0.1001	1.39	0.1656
Condition: dice & evs_c	0.6823	0.2797	2.44	0.0147
Group: dyslexic & evs_c	0.8354	0.2687	3.11	0.0019
Condition: dice & gaze_c	-0.0018	0.0011	-1.66	0.0970
Group: dyslexic & gaze_c	-0.0024	0.0012	-2.02	0.0437
evs_c & gaze_c	-0.0075	0.0029	-2.57	0.0103
Condition: dice & Group: dyslexic & evs_c	-0.6611	0.3634	-1.82	0.0689
Condition: dice & Group: dyslexic & gaze_c	0.0008	0.0013	0.61	0.5421
Condition: dice & evs_c & gaze_c	0.0076	0.0040	1.89	0.0588
Group: dyslexic & evs_c & gaze_c	0.0011	0.0035	0.31	0.7556
Condition: dice & Group: dyslexic & evs_c & gaze_c	-0.0058	0.0050	-1.17	0.2416
Residual	0.0851

m10 is the full Condition × Group × evs × gaze factorial. The three highest-order terms that involve the evs × gaze product are not significant, so we drop them: the Condition × Group × evs × gaze four-way, the Group × evs × gaze three-way, and the Condition × evs × gaze three-way. What remains (m09) is the model whose coefficients correspond to Table 3 of Pan et al. (2013). (The z statistics here differ slightly from the published t values because of software differences — ML vs. REML, contrast parameterization — but the substantive pattern is the same.)

m09 = fit_or_restore("eyevoicespan_m09.json", MixedModel,
  @formula(lran ~ 1 + Condition + Group + evs_c + gaze_c +
    Condition & Group + Condition & evs_c + Condition & gaze_c +
    Group & evs_c + Group & gaze_c + evs_c & gaze_c +
    Condition & Group & evs_c + Condition & Group & gaze_c +
    (1 | Subj)),
  dat; contrasts, progress)

	Est.	SE	z	p	σ_Subj
(Intercept)	2.9984	0.0683	43.88	<1e-99	0.0982
Condition: dice	0.0770	0.0760	1.01	0.3108
Group: dyslexic	0.0117	0.0789	0.15	0.8825
evs_c	-0.6308	0.1440	-4.38	<1e-04
gaze_c	0.0036	0.0009	4.17	<1e-04
Condition: dice & Group: dyslexic	0.1299	0.0971	1.34	0.1807
Condition: dice & evs_c	0.6085	0.1998	3.05	0.0023
Condition: dice & gaze_c	-0.0013	0.0009	-1.41	0.1598
Group: dyslexic & evs_c	0.6163	0.1842	3.35	0.0008
Group: dyslexic & gaze_c	-0.0020	0.0011	-1.90	0.0569
evs_c & gaze_c	-0.0051	0.0015	-3.46	0.0006
Condition: dice & Group: dyslexic & evs_c	-0.5512	0.2010	-2.74	0.0061
Condition: dice & Group: dyslexic & gaze_c	0.0002	0.0011	0.18	0.8570
Residual	0.0861

MixedModels.likelihoodratiotest(m00, m05, m09, m10)

	model-dof	-2 logLik	χ²	χ²-dof	P(>χ²)
lran ~ 1 + Condition + Group + evs_c + gaze_c + (1 \| Subj)	7	134
lran ~ 1 + Condition + Group + evs_c + gaze_c + Condition & Group + (1 \| Subj)	8	144	10	1	0.0013
lran ~ 1 + Condition + Group + evs_c + gaze_c + Condition & Group + Condition & evs_c + Condition & gaze_c + Group & evs_c + Group & gaze_c + evs_c & gaze_c + Condition & Group & evs_c + Condition & Group & gaze_c + (1 \| Subj)	15	160	16	7	0.0258
lran ~ 1 + Condition + Group + evs_c + gaze_c + Condition & Group + Condition & evs_c + Group & evs_c + Condition & gaze_c + Group & gaze_c + evs_c & gaze_c + Condition & Group & evs_c + Condition & Group & gaze_c + Condition & evs_c & gaze_c + Group & evs_c & gaze_c + Condition & Group & evs_c & gaze_c + (1 \| Subj)	18	164	4	3	0.2197

Adding Condition × Group (m05) and the covariate interactions retained in m09 each improve the fit reliably; the extra terms in m10 do not. m09 is the preferred model.

4.2 The preferred model

VarCorr(m09)

	Column	Variance	Std.Dev
Subj	(Intercept)	0.0096367	0.0981667
Residual		0.0074142	0.0861058

m09

	Est.	SE	z	p	σ_Subj
(Intercept)	2.9984	0.0683	43.88	<1e-99	0.0982
Condition: dice	0.0770	0.0760	1.01	0.3108
Group: dyslexic	0.0117	0.0789	0.15	0.8825
evs_c	-0.6308	0.1440	-4.38	<1e-04
gaze_c	0.0036	0.0009	4.17	<1e-04
Condition: dice & Group: dyslexic	0.1299	0.0971	1.34	0.1807
Condition: dice & evs_c	0.6085	0.1998	3.05	0.0023
Condition: dice & gaze_c	-0.0013	0.0009	-1.41	0.1598
Group: dyslexic & evs_c	0.6163	0.1842	3.35	0.0008
Group: dyslexic & gaze_c	-0.0020	0.0011	-1.90	0.0569
evs_c & gaze_c	-0.0051	0.0015	-3.46	0.0006
Condition: dice & Group: dyslexic & evs_c	-0.5512	0.2010	-2.74	0.0061
Condition: dice & Group: dyslexic & gaze_c	0.0002	0.0011	0.18	0.8570
Residual	0.0861

The coefficients of interest:

evs_c (−0.63) — within digit/control, a larger eye–voice span predicts faster (lower log) naming.
Condition: dice & Group: dyslexic & evs_c (−0.55, p ≈ .006) — the three-way interaction: the EVS benefit is distributed differently across cells (Figure 1 of the paper).
evs_c & gaze_c (−0.005, p ≈ .001) — the two covariates interact.
gaze_c (+0.0036) — longer gaze durations predict slower naming, as expected.

5 Diagnostics

Code

draw(
  data((; f=fitted(m09), r=residuals(m09))) *
  mapping(:f => "Fitted log(RAN)", :r => "Residual (m09)") *
  visual(Scatter),
)

Figure 4: Residuals vs. fitted values for m09.

Code

qqnorm(residuals(m09); qqline=:fitrobust)

Figure 5: Normal quantile plot of the m09 residuals.

Code

caterpillar(m09, :Subj)

Figure 6: Prediction intervals on the by-child random intercepts (m09).

6 Partial-effect plots (the `remef()` analogue)

The R script visualizes the higher-order interactions with partial effects computed by remef(): it takes the observed response, removes the contribution of the nuisance terms (here everything involving gaze and the random intercept), and keeps the terms of interest plus the residual. The result is an “adjusted” response that isolates one slice of the model.

MixedModelsExtras.partial_fitted is the Julia counterpart. It returns the fitted values for a chosen set of coefficients; adding the residuals back reproduces remef(..., keep=TRUE).

# keep every fixed effect that does NOT involve gaze; drop the by-child intercept
keep_evs = filter(c -> !occursin("gaze_c", c), coefnames(m09))
# keep every fixed effect that does NOT involve evs; drop the by-child intercept
keep_gaze = filter(c -> !occursin("evs_c", c), coefnames(m09))

dp = select(dat, :Subj, :Group, :Condition, :evs, :gaze, :lran)
dp.partial_evs =
  partial_fitted(m09, keep_evs, Dict(:Subj => String[]); mode=:include) .+ residuals(m09)
dp.partial_gaze =
  partial_fitted(m09, keep_gaze, Dict(:Subj => String[]); mode=:include) .+ residuals(m09)
dp.cg = categorical(string.(dp.Condition, " / ", dp.Group))
first(dp, 6)

6×9 DataFrame

Row	Subj	Group	Condition	evs	gaze	lran	partial_evs	partial_gaze	cg
	Cat…	Cat…	Cat…	Float64	Int32	Float64	Float64	Float64	Cat…
1	S101	control	digit	1.33	242	2.62394	2.81295	2.52721	digit / control
2	S101	control	dice	0.97	397	2.97604	2.97662	2.96624	dice / control
3	S102	control	digit	0.84	318	2.58626	3.00638	2.68389	digit / control
4	S102	control	dice	1.06	372	2.85819	3.00322	2.93693	dice / control
5	S103	control	digit	0.83	403	3.1206	3.08125	3.05705	digit / control
6	S103	control	dice	0.8	487	3.31309	3.049	3.2433	dice / control

For plotting we stack the observed response and each partial response into long form, so that “observed vs. adjusted” becomes a facet column with a single shared legend.

stack_partial(col, label) =
  insertcols!(select(dp, :evs, :gaze, :cg, col => :y), :kind => label)

dp_evs = vcat(stack_partial(:lran, "observed"), stack_partial(:partial_evs, "adjusted"))
dp_gaze = vcat(stack_partial(:lran, "observed"), stack_partial(:partial_gaze, "adjusted"))
first(dp_evs, 4)

4×5 DataFrame

Row	evs	gaze	cg	y	kind
	Float64	Int32	Cat…	Float64	String
1	1.33	242	digit / control	2.62394	observed
2	0.97	397	dice / control	2.97604	observed
3	0.84	318	digit / control	2.58626	observed
4	1.06	372	dice / control	2.85819	observed

How the specification maps onto R’s remef()

remef(m09, fix = c(1, "condition:group:evs_c"), keep = TRUE, grouping = TRUE, ran = NULL) keeps the intercept, the three-way condition:group:evs_c term, and (because grouping = TRUE) every lower-order term built from condition, group, and evs_c; it removes everything involving gaze_c and all random effects. filter(c -> !occursin("gaze_c", c), coefnames(m09)) selects exactly that set of fixed-effect coefficients, and Dict(:Subj => String[]) drops the random intercept.

6.1 Figure 1 — `Condition × Group × EVS`

Code

draw(
  data(dp_evs) *
  mapping(
    :evs => "Eye–voice span",
    :y => "log(RAN)";
    color=:cg => "Condition / Group",
    col=:kind,
  ) *
  (visual(Scatter; alpha=0.5) + linear());
  facet=(; linkyaxes=:all),
)

Figure 7: Left: observed log(RAN) vs. eye–voice span. Right: partial effect from m09 (gaze contribution and random intercept removed).

The partial plot sharpens the pattern noted in the paper: the EVS benefit (negative slope) is present for digit naming in the control group but is weak or absent in the dice / dyslexic cell.

6.2 Figure 2 — `Condition × Group × Gaze`

Code

draw(
  data(dp_gaze) *
  mapping(
    :gaze => "Gaze duration [ms]",
    :y => "log(RAN)";
    color=:cg => "Condition / Group",
    col=:kind,
  ) *
  (visual(Scatter; alpha=0.5) + linear());
  facet=(; linkyaxes=:all),
)

Figure 8: Left: observed log(RAN) vs. gaze duration. Right: partial effect from m09 (EVS contribution and random intercept removed).

7 Marginal-effect plots with `Effects.jl`

partial_fitted adjusts observed data points. Effects.jl instead evaluates the model’s prediction on a regular grid, holding the variables that are not on the grid at a typical value (the mean; for the centered covariates that is 0). This is the Julia analogue of R’s effects / emmeans.

7.1 EVS effect by Condition and Group

evsgrid = Dict(
  :evs_c => range(extrema(dat.evs_c)...; length=50),
  :Condition => levels(dat.Condition),
  :Group => levels(dat.Group),
)
eff_evs = effects(evsgrid, m09)
eff_evs.evs = eff_evs.evs_c .+ mean(dat.evs)
eff_evs.cg = categorical(string.(eff_evs.Condition, " / ", eff_evs.Group))
first(eff_evs, 6)

6×9 DataFrame

Row	evs_c	Group	Condition	lran	err	lower	upper	evs	cg
	Float64	Cat…	Cat…	Float64	Float64	Float64	Float64	Float64	Cat…
1	-0.691607	control	digit	3.43466	0.134987	3.29967	3.56965	0.18	digit / control
2	-0.659974	control	digit	3.4147	0.131029	3.28368	3.54573	0.211633	digit / control
3	-0.628342	control	digit	3.39475	0.12711	3.26764	3.52186	0.243265	digit / control
4	-0.596709	control	digit	3.3748	0.123235	3.25156	3.49803	0.274898	digit / control
5	-0.565077	control	digit	3.35484	0.119408	3.23543	3.47425	0.306531	digit / control
6	-0.533444	control	digit	3.33489	0.115635	3.21925	3.45052	0.338163	digit / control

Code

draw(
  data(eff_evs) *
  mapping(
    :evs => "Eye–voice span",
    color=:cg => "Condition / Group",
  ) * (
    mapping(:lran => "log(RAN)") * visual(Lines) + 
    mapping(:lower, :upper) * visual(Band; alpha=0.3)
    )
  )

Figure 9: Model-predicted log(RAN) as a function of eye–voice span, with 95% confidence bands, by Condition and Group.

7.2 Gaze effect by Condition and Group

Code

let
  gazegrid = Dict(
    :gaze_c => range(extrema(dat.gaze_c)...; length=50),
    :Condition => levels(dat.Condition),
    :Group => levels(dat.Group),
  )
  eff_gaze = effects(gazegrid, m09)
  eff_gaze.gaze = eff_gaze.gaze_c .+ mean(dat.gaze)
  eff_gaze.cg = categorical(string.(eff_gaze.Condition, " / ", eff_gaze.Group))
  draw(
    data(eff_gaze) *
    mapping(:gaze => "Gaze duration [ms]"; color=:cg => "Condition / Group") * (
      mapping(:lran => "log(RAN)") * visual(Lines) +
      mapping(:lower, :upper) * visual(Band; alpha=0.3)
  ))
end

Figure 10: Model-predicted log(RAN) as a function of gaze duration, with 95% confidence bands, by Condition and Group.

7.3 `evs × gaze` interaction

The reliable evs_c & gaze_c term means the gaze slope changes with EVS. We show it by evaluating the model at three EVS levels (roughly the tertiles).

Code

let
  qs = quantile(dat.evs, [1/6, 1/2, 5/6])
  grid = Dict(
    :gaze_c => range(extrema(dat.gaze_c)...; length=50),
    :evs_c => qs .- mean(dat.evs),
  )
  eff = effects(grid, m09)
  eff.gaze = eff.gaze_c .+ mean(dat.gaze)
  eff.evs_level = categorical(
    round.(eff.evs_c .+ mean(dat.evs); digits=2);
  )
  draw(
    data(eff) *
    mapping(:gaze => "Gaze duration [ms]";
            color=:evs_level => "Eye–voice span") *
    (mapping(:lran => "log(RAN)") * visual(Lines) + 
     mapping(:lower, :upper) * visual(Band; alpha=0.3))
  )
end

Figure 11: The gaze-duration effect at small, medium, and large eye–voice span (model prediction, 95% bands).

8 See also

Transformations and Effects — the Effects.jl workflow and Box–Cox transformations in more detail.
The Emotikon Project — model complexification and PCA of the random effects on a large dataset.
R’s remef and the partial_fitted docstring in MixedModelsExtras.

9 References

Pan, J., Yan, M., Laubrock, J., Shu, H., & Kliegl, R. (2013). Eye–voice span during rapid automatized naming of digits and dice in Chinese normal and dyslexic children. Developmental Science, 16(6), 967–979. https://doi.org/10.1111/desc.12075

9.1 Exercises

Contrast coding. Refit m09 with EffectsCoding for Condition and Group. Which coefficients change, which stay the same, and what does the intercept estimate now?

Solution

The omnibus likelihood-ratio tests, VarCorr, the residual standard deviation, and the fitted values are all unchanged — the model space is the same. The intercept now estimates the grand mean (the unweighted average over the four Condition × Group cells) instead of the digit/control cell, and every lower-order Condition/Group coefficient becomes a deviation from that grand mean rather than a simple difference from a reference level. The highest-order interaction coefficient is unchanged up to a scale factor.

Partial vs. marginal. Figure 7 (partial effects) and Figure 9 (marginal effects) both show the EVS slopes by cell. What does each one add that the other does not?

Solution

The partial-effect plot keeps the residuals, so it shows the scatter of individual observations around the adjusted regression line — useful for spotting influential points and for judging how much of the variance the term actually explains. The marginal-effect plot shows the model’s prediction with a proper confidence band and no leftover noise from other terms, which makes the estimated slopes and their uncertainty directly comparable across cells. Use partial effects to inspect the data given the model; use marginal effects to communicate the model.

Why center? The covariates were centered before fitting. What would change in the m09 coefficient table if evs and gaze entered on their raw scales, and would any of the omnibus tests change?

Solution

Only the interpretation of the lower-order terms shifts. With raw covariates, Condition: dice would estimate the digit-vs-dice difference at evs = 0 and gaze = 0 — an extrapolation far outside the data — inflating its standard error and making it hard to interpret. The intercept would likewise refer to evs = gaze = 0. The higher-order interaction coefficients, the variance components, the fitted values, and every likelihood-ratio test are invariant to centering.

This page was rendered from git revision 29e8d33 using Quarto 1.10.18 and Julia 1.12.7.

1 Background

2 Packages

3 Read and prepare the data

3.1 Does the response need a transformation?

3.2 Condition and Group means

4 Linear mixed models

4.1 Model building

4.2 The preferred model

5 Diagnostics

6 Partial-effect plots (the remef() analogue)

6.1 Figure 1 — Condition × Group × EVS

6.2 Figure 2 — Condition × Group × Gaze

7 Marginal-effect plots with Effects.jl

7.1 EVS effect by Condition and Group

7.2 Gaze effect by Condition and Group

7.3 evs × gaze interaction

8 See also

9 References

9.1 Exercises

6 Partial-effect plots (the `remef()` analogue)

6.1 Figure 1 — `Condition × Group × EVS`

6.2 Figure 2 — `Condition × Group × Gaze`

7 Marginal-effect plots with `Effects.jl`

7.3 `evs × gaze` interaction