After working through this page you will be able to:
remef / MixedModelsExtras.partial_fitted), and the predicted (marginal) view (Effects.jl);m09 off a consistent three-panel figure;Prerequisite: Eye–Voice Span in Rapid Automatized Naming — this page reuses its preferred model m09 and its cached fit. Read that page first for the background, the data, and the model-building steps.
Data used: pylsk13.rda (examples/eyevoicespan/).
Model m09 from the companion page has four reliable interaction terms. Each one describes how the slope of one predictor changes with another, but a raw scatterplot rarely shows that directly: the other twelve terms in the model, the by-child random intercept, and the residual noise are all superimposed on it.
Three views pull the term apart:
| View | What it shows | How it is computed |
|---|---|---|
| Observed | log(RAN) against the predictor, split by the moderator, with an ordinary least-squares line per group |
nothing removed — the raw data |
| Partial | the same, after subtracting every model term except the target interaction and its lower-order relatives, and the random intercept | partial_fitted(m09, keep; mode=:include) .+ residuals(m09) — the Julia analogue of R’s remef(..., keep=TRUE, grouping=TRUE, ran=NULL) |
| Predicted | the model’s fitted mean on a regular grid, with a 95% confidence band, holding the other predictors at typical values | Effects.effects(grid, m09; level=0.95) |
The observed panel is honest about scatter but confounded. The partial panel keeps the residual scatter (so you can still judge influence and fit) while removing the confounds. The predicted panel removes the scatter too and adds proper inferential uncertainty, but it averages over the predictors that are not on the grid — so it can look weaker than the partial view when a moderator you have collapsed also interacts with the x-axis variable.
using AlgebraOfGraphics
using AlgebraOfGraphics: linear
using CairoMakie
using CategoricalArrays
using DataFrames
using Effects
using MixedModels
using MixedModelsExtras
using RData
using Statistics
using StatsBase
using SMLP2026: fit_or_restore
const progress = false
set_aog_theme!()dat = DataFrame(load(joinpath(@__DIR__, "pylsk13.rda"))["dat"])
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)
contrasts = Dict(
:Condition => DummyCoding(; base="digit"),
:Group => DummyCoding(; base="control"),
)
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 |
The four interaction terms and their coefficients:
let
targets = ["Condition: dice & evs_c", "Group: dyslexic & evs_c", "evs_c & gaze_c",
"Condition: dice & Group: dyslexic & evs_c"]
filter(:Name => in(targets), DataFrame(coeftable(m09)))
end| Row | Name | Coef. | Std. Error | z | Pr(>|z|) |
|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Float64 | |
| 1 | Condition: dice & evs_c | 0.60853 | 0.199783 | 3.04595 | 0.00231946 |
| 2 | Group: dyslexic & evs_c | 0.616318 | 0.184165 | 3.34656 | 0.000818218 |
| 3 | evs_c & gaze_c | -0.00512145 | 0.0014823 | -3.45508 | 0.000550125 |
| 4 | Condition: dice & Group: dyslexic & evs_c | -0.551238 | 0.201037 | -2.74198 | 0.00610709 |
All four are reliable: Condition & evs_c (p ≈ .002), Group & evs_c (p ≈ .001), evs_c & gaze_c (p ≈ .001), and the three-way Condition & Group & evs_c (p ≈ .006).
partial(keep) returns the observed response with every term not in keep removed, plus the by-child random intercept removed — exactly what remef(m09, fix = keep, keep = TRUE, grouping = TRUE, ran = NULL) produces in R. grouping = TRUE means “also keep every lower-order term built from the same variables”, which here we spell out explicitly in each keep vector.
partial(keep) = partial_fitted(m09, keep, Dict(:Subj => String[]); mode=:include) .+ residuals(m09)
keep_C_evs = ["(Intercept)", "Condition: dice", "evs_c", "Condition: dice & evs_c"]
keep_G_evs = ["(Intercept)", "Group: dyslexic", "evs_c", "Group: dyslexic & evs_c"]
keep_evs_gz = ["(Intercept)", "evs_c", "gaze_c", "evs_c & gaze_c"]
keep_CG_evs = filter(c -> !occursin("gaze_c", c), coefnames(m09)) # everything without gaze
dp = select(dat, :Subj, :Condition, :Group, :evs, :gaze, :evs_c, :gaze_c, :lran)
dp.pf_C_evs = partial(keep_C_evs)
dp.pf_G_evs = partial(keep_G_evs)
dp.pf_evs_gz = partial(keep_evs_gz)
dp.pf_CG_evs = partial(keep_CG_evs)
dp.cg = categorical(string.(dp.Condition, " / ", dp.Group))
# tertiles of eye–voice span, for the EVS × Gaze figure
cut = quantile(dat.evs, [1/3, 2/3])
evsgrp(x) = x ≤ cut[1] ? "low" : x ≤ cut[2] ? "mid" : "high"
dp.evs_grp = categorical(evsgrp.(dp.evs); levels=["low", "mid", "high"])
tert_med = [median(dat.evs[evsgrp.(dat.evs) .== g]) for g in ["low", "mid", "high"]]
first(dp, 6)| Row | Subj | Condition | Group | evs | gaze | evs_c | gaze_c | lran | pf_C_evs | pf_G_evs | pf_evs_gz | pf_CG_evs | cg | evs_grp |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Cat… | Cat… | Cat… | Float64 | Int32 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Cat… | Cat… | |
| 1 | S101 | digit | control | 1.33 | 242 | 0.458393 | -160.429 | 2.62394 | 2.81295 | 2.81295 | 2.61468 | 2.81295 | digit / control | high |
| 2 | S101 | dice | control | 0.97 | 397 | 0.0983929 | -5.42857 | 2.97604 | 2.97662 | 2.83973 | 2.82301 | 2.97662 | dice / control | mid |
| 3 | S102 | digit | control | 0.84 | 318 | -0.0316071 | -84.4286 | 2.58626 | 3.00638 | 3.00638 | 2.69016 | 3.00638 | digit / control | mid |
| 4 | S102 | dice | control | 1.06 | 372 | 0.188393 | -30.4286 | 2.85819 | 3.00322 | 2.81156 | 2.73188 | 3.00322 | dice / control | high |
| 5 | S103 | digit | control | 0.83 | 403 | -0.0416071 | 0.571429 | 3.1206 | 3.08125 | 3.08125 | 3.08342 | 3.08125 | digit / control | mid |
| 6 | S103 | dice | control | 0.8 | 487 | -0.0716071 | 84.5714 | 3.31309 | 3.049 | 3.01555 | 3.34963 | 3.049 | dice / control | mid |
"""
three_panel(; obs_df, xcol, xlab, obs_y, par_y, pred_df, color, colorlab, title)
One figure row: observed | partial | predicted, sharing the y-axis, with one
shared colour legend. `pred_df` must carry `:yhat`, `:lower`, `:upper` and the
same `color` column as `obs_df`.
"""
function three_panel(; obs_df, xcol, xlab, obs_y, par_y, pred_df, color, colorlab, title)
fig = Figure(; size=(1050, 380))
ylims = (nothing, nothing, 2.0, 4.0)
pobs = data(obs_df) * mapping(xcol => xlab, obs_y => "log(RAN)"; color=color => colorlab) *
(visual(Scatter; alpha=0.35) + linear())
ppar = data(obs_df) * mapping(xcol => xlab, par_y => "adjusted log(RAN)"; color=color => colorlab) *
(visual(Scatter; alpha=0.35) + linear())
ppred = data(pred_df) * mapping(xcol => xlab, color=color => colorlab) *
(mapping(:yhat => "predicted log(RAN)") * visual(Lines) +
mapping(:lower, :upper) * visual(Band; alpha=0.3))
g = draw!(fig[1, 1], pobs; axis=(; limits=ylims, title="Observed"))
draw!(fig[1, 2], ppar; axis=(; limits=ylims, title="Partial (remef)"))
draw!(fig[1, 3], ppred; axis=(; limits=ylims, title="Predicted (Effects)"))
legend!(fig[1, 4], g)
Label(fig[0, 1:4], title; fontsize=16, font=:bold)
return fig
endMain.Notebook.three_panellet
grid = Dict(:evs_c => range(extrema(dat.evs_c)...; length=50), :Condition => levels(dat.Condition))
e = effects(grid, m09; level=0.95)
rename!(e, :lran => :yhat)
e.evs = e.evs_c .+ mean(dat.evs)
three_panel(; obs_df=dp, xcol=:evs, xlab="Eye–voice span", obs_y=:lran, par_y=:pf_C_evs,
pred_df=e, color=:Condition, colorlab="Condition", title="Condition × EVS")
end
digit retains a slope once Group and gaze are removed. Predicted: same pattern, averaged over Group (hence weaker) and with model uncertainty.
A larger eye–voice span goes with faster naming, and the partial panel shows this benefit is carried almost entirely by digit naming: once the (also reliable) Group × EVS term and everything involving gaze are stripped out, the dice slope is essentially flat. The predicted panel collapses over Group, so its digit/dice lines separate less — the interaction is real but the marginal prediction dilutes it.
let
grid = Dict(:evs_c => range(extrema(dat.evs_c)...; length=50), :Group => levels(dat.Group))
e = effects(grid, m09; level=0.95)
rename!(e, :lran => :yhat)
e.evs = e.evs_c .+ mean(dat.evs)
three_panel(; obs_df=dp, xcol=:evs, xlab="Eye–voice span", obs_y=:lran, par_y=:pf_G_evs,
pred_df=e, color=:Group, colorlab="Group", title="Group × EVS")
end
The partial panel isolates the mirror-image of Figure 1: the steep negative EVS slope belongs to the control group, while the dyslexic slope is shallow. Substantively, a wide eye–voice span — reading ahead of the voice — pays off more for the children who are already fluent.
let
grid = Dict(:gaze_c => range(extrema(dat.gaze_c)...; length=50),
:evs_c => tert_med .- mean(dat.evs))
e = effects(grid, m09; level=0.95)
rename!(e, :lran => :yhat)
e.gaze = e.gaze_c .+ mean(dat.gaze)
pts = tert_med .- mean(dat.evs)
e.evs_grp = categorical(["low", "mid", "high"][[findmin(abs.(pts .- v))[2] for v in e.evs_c]];
levels=["low", "mid", "high"])
three_panel(; obs_df=dp, xcol=:gaze, xlab="Gaze duration [ms]", obs_y=:lran, par_y=:pf_evs_gz,
pred_df=e, color=:evs_grp, colorlab="EVS tertile", title="EVS × Gaze")
end
This is the clearest demonstration of what the partial view buys you. In the raw data the gaze-duration slope looks the same at every eye–voice span. After removing the Condition, Group, and random-intercept contributions, a fan opens up: longer gaze durations cost more (steeper positive slope) when the reader is not looking far ahead. The predicted panel confirms the crossing pattern with confidence bands.
let
grid = Dict(:evs_c => range(extrema(dat.evs_c)...; length=50),
:Condition => levels(dat.Condition), :Group => levels(dat.Group))
e = effects(grid, m09; level=0.95)
rename!(e, :lran => :yhat)
e.evs = e.evs_c .+ mean(dat.evs)
e.cg = categorical(string.(e.Condition, " / ", e.Group))
three_panel(; obs_df=dp, xcol=:evs, xlab="Eye–voice span", obs_y=:lran, par_y=:pf_CG_evs,
pred_df=e, color=:cg, colorlab="Condition / Group", title="Condition × Group × EVS")
end
Putting Condition and Group back together: the partial panel shows three flat cells and one steeply declining one — digit naming in control children. The EVS benefit found in Figures 1 and 2 is not additive; it is a property of that single cell. Because the three-way grid here retains all the non-gaze structure, the predicted panel now agrees closely with the partial panel — nothing has been averaged away.
pf_CG_evs with Dict(:Subj => ["(Intercept)"]) instead of Dict(:Subj => String[]). What changes in the partial panel, and which R remef argument does this correspond to?Each point moves vertically by that child’s estimated random intercept (± ~0.1 on the log scale), so the within-cell scatter grows and the cells separate a little more by their child composition. It corresponds to remef(..., ran = list("(Intercept)")) (keep the by-subject intercept) rather than ran = NULL.
The conclusion is unchanged — the gaze slope still decreases monotonically as EVS grows. The cost is only visual: five overlapping bands are harder to read than three, and the extreme quantiles are supported by fewer observations, so their bands are wider.
Condition × gaze_c (p ≈ .16 in m09). What do you expect the partial panel to look like, and why is that the right null result to show learners?The two partial slopes should be nearly parallel — the interaction coefficient is small and uncertain. Showing it next to the significant interactions makes the point that the partial-effect plot is not a device for manufacturing patterns: when the term is null, the plot looks null.
This page was rendered from git revision 29e8d33 using Quarto 1.10.18 and Julia 1.12.7.