diff --git a/.gitattributes b/.gitattributes
index 03d92dd..db95ba5 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -3,3 +3,5 @@
 *.RDs filter=lfs diff=lfs merge=lfs -text
 *.RDS filter=lfs diff=lfs merge=lfs -text
 *.rds filter=lfs diff=lfs merge=lfs -text
+*.pdf filter=lfs diff=lfs merge=lfs -text
+*.rda filter=lfs diff=lfs merge=lfs -text
diff --git a/_quarto.yml b/_quarto.yml
index 54d0dab..0516852 100644
--- a/_quarto.yml
+++ b/_quarto.yml
@@ -110,6 +110,10 @@ website:
           contents:
           - examples/visual-attention/kwdyz11.qmd
           - examples/visual-attention/kkl15.qmd
+        - section: "Eye–Voice Span in Rapid Automatized Naming"
+          contents:
+          - examples/eyevoicespan/pylsk13.qmd
+          - examples/eyevoicespan/pylsk13_interactions.qmd
       - section: "Reference & Appendix"
         contents:
         - reference/useful_packages.qmd
diff --git a/contrasts/contrasts_kwdyz11.qmd b/contrasts/contrasts_kwdyz11.qmd
index 06bf306..f51b51b 100644
--- a/contrasts/contrasts_kwdyz11.qmd
+++ b/contrasts/contrasts_kwdyz11.qmd
@@ -38,6 +38,7 @@ progress = false
 ```{julia}
 #| echo: false
 using SMLP2026: fit_or_restore
+using SMLP2026.RKUtils
 ```
 
 # A word of caution {#sec-caution}
diff --git a/examples/emotikon/fggk21.qmd b/examples/emotikon/fggk21.qmd
index 3de604a..41a9e0d 100644
--- a/examples/emotikon/fggk21.qmd
+++ b/examples/emotikon/fggk21.qmd
@@ -659,10 +659,6 @@ issingular(m2)
 
 Depending on the random number generator seed, the model may or may not be supported in the alternative parameterization of scores. The fixed-effects profile is not affected (see the model comparisons below).
 
-:::{.callout-caution}
-RK: The order of RE terms is critical. In formula `f2` the `zerocorr()` term must be placed last as shown. If it is placed first, School-related and Child-related CPs are estimated/reported (?) as zero. This was not the case for formula `m1`. Thus, it appears to be related to the `0`-intercepts in School and Child terms. Need a reprex.
-:::
-
 ```{julia}
 VarCorr(m2)
 ```
diff --git a/examples/eyevoicespan/Pan_etal.DevSci.2013.pdf b/examples/eyevoicespan/Pan_etal.DevSci.2013.pdf
new file mode 100644
index 0000000..d4e8f26
Binary files /dev/null and b/examples/eyevoicespan/Pan_etal.DevSci.2013.pdf differ
diff --git a/examples/eyevoicespan/pylsk13.qmd b/examples/eyevoicespan/pylsk13.qmd
new file mode 100644
index 0000000..44e6903
--- /dev/null
+++ b/examples/eyevoicespan/pylsk13.qmd
@@ -0,0 +1,493 @@
+---
+title: "Eye–Voice Span in Rapid Automatized Naming: Pan et al. (2013)"
+engine: julia
+author:
+- Reinhold Kliegl
+- Jinger Pan
+julia:
+  exeflags: ["--project", "--threads=auto"]
+---
+
+::: {.objectives title="Learning objectives"}
+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`).
+:::
+
+::: {.callout-note title="Before you start"}
+**Prerequisites:** [Analysis of the sleepstudy data](../../lmm-intro/sleepstudy.qmd) for the modeling vocabulary; [Transformations and Effects](../../contrasts/transformation.qmd) for `Effects.jl`; [Creating multi-panel plots](../../visualization/AoGPlots.qmd) 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.
+:::
+
+# Background
+
+@Pan2013 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 `Condition`s (digit/dice) and the two `Group`s (control/dyslexic)?
+
+The analysis follows `pylsk13_v4.R`:
+
+1. **Setup** — read the data, check whether the response needs a transformation, center the covariates.
+2. **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.
+3. **Visualization** — partial-effect plots (Figures 1 and 2 of the paper), marginal-effect plots, and residual diagnostics.
+
+# Packages
+
+```{julia}
+#| code-fold: true
+#| output: false
+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!()
+```
+
+# 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).
+
+```{julia}
+dat = DataFrame(load(joinpath(@__DIR__, "pylsk13.rda"))["dat"])
+```
+
+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`.
+
+```{julia}
+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)
+```
+
+## 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.
+
+```{julia}
+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)
+```
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Box–Cox profile for the RAN response. The optimum is close to λ = 0, i.e. a log transform."
+#| label: fig-boxcox
+boxcoxplot(bc; conf_level=0.95)
+```
+
+The profile peaks near `λ = 0`, so `log(ran)` (= `lran`) is a reasonable choice and
+matches the published analysis.
+
+## Condition and Group means
+
+```{julia}
+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
+```
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "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."
+#| label: fig-cellmeans
+draw(
+  data(cell_means) *
+  mapping(:Condition, :lran_m => "log(RAN)"; color=:Group, group=:Group) *
+  visual(ScatterLines; markersize=16),
+)
+```
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Comparative density of log(RAN) by Condition."
+#| label: fig-density
+draw(
+  data(dat) *
+  mapping(:lran => "log(RAN)"; color=:Condition) *
+  density(),
+)
+```
+
+# 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 @Pan2013.
+
+::: {.callout-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.
+:::
+
+```{julia}
+contrasts = Dict(
+  :Condition => DummyCoding(; base="digit"),
+  :Group => DummyCoding(; base="control"),
+)
+```
+
+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.
+
+## Model building
+
+```{julia}
+m00 = fit_or_restore("eyevoicespan_m00.json", MixedModel,
+  @formula(lran ~ 1 + Condition + Group + evs_c + gaze_c + (1 | Subj)),
+  dat; contrasts, progress)
+```
+
+```{julia}
+m05 = fit_or_restore("eyevoicespan_m05.json", MixedModel,
+  @formula(lran ~ 1 + Condition * Group + evs_c + gaze_c + (1 | Subj)),
+  dat; contrasts, progress)
+```
+
+```{julia}
+m10 = fit_or_restore("eyevoicespan_m10.json", MixedModel,
+  @formula(lran ~ 1 + Condition * Group * evs_c * gaze_c + (1 | Subj)),
+  dat; contrasts, progress)
+```
+
+`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 @Pan2013. (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.)
+
+```{julia}
+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)
+```
+
+```{julia}
+MixedModels.likelihoodratiotest(m00, m05, m09, m10)
+```
+
+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.
+
+## The preferred model
+
+```{julia}
+VarCorr(m09)
+```
+
+```{julia}
+m09
+```
+
+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.
+
+# Diagnostics
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Residuals vs. fitted values for m09."
+#| label: fig-resvsfitted
+draw(
+  data((; f=fitted(m09), r=residuals(m09))) *
+  mapping(:f => "Fitted log(RAN)", :r => "Residual (m09)") *
+  visual(Scatter),
+)
+```
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Normal quantile plot of the m09 residuals."
+#| label: fig-qq
+qqnorm(residuals(m09); qqline=:fitrobust)
+```
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Prediction intervals on the by-child random intercepts (m09)."
+#| label: fig-caterpillar
+caterpillar(m09, :Subj)
+```
+
+# 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)`.
+
+```{julia}
+# 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)
+```
+
+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.
+
+```{julia}
+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)
+```
+
+::: {.callout-note collapse="true" title="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.
+:::
+
+## Figure 1 — `Condition × Group × EVS`
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Left: observed log(RAN) vs. eye–voice span. Right: partial effect from m09 (gaze contribution and random intercept removed)."
+#| label: fig-partial-evs
+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),
+)
+```
+
+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.
+
+## Figure 2 — `Condition × Group × Gaze`
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Left: observed log(RAN) vs. gaze duration. Right: partial effect from m09 (EVS contribution and random intercept removed)."
+#| label: fig-partial-gaze
+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),
+)
+```
+
+# 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`.
+
+## EVS effect by Condition and Group
+
+```{julia}
+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)
+```
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Model-predicted log(RAN) as a function of eye–voice span, with 95% confidence bands, by Condition and Group."
+#| label: fig-effects-evs
+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)
+    )
+  )
+```
+
+## Gaze effect by Condition and Group
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "Model-predicted log(RAN) as a function of gaze duration, with 95% confidence bands, by Condition and Group."
+#| label: fig-effects-gaze
+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
+```
+
+## `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).
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "The gaze-duration effect at small, medium, and large eye–voice span (model prediction, 95% bands)."
+#| label: fig-effects-evsxgaze
+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
+```
+
+# See also
+
+- [Transformations and Effects](../../contrasts/transformation.qmd) — the `Effects.jl` workflow and Box–Cox transformations in more detail.
+- [The Emotikon Project](../emotikon/fggk21.qmd) — model complexification and PCA of the random effects on a large dataset.
+- R's [`remef`](https://github.com/hohenstein/remef) and the `partial_fitted` docstring in `MixedModelsExtras`.
+
+# References
+
+::: {#refs}
+:::
+
+## Exercises
+
+1. **Contrast coding.** Refit `m09` with `EffectsCoding` for `Condition` and `Group`. Which coefficients change, which stay the same, and what does the intercept estimate now?
+
+::: {.callout-note collapse="true" title="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.
+:::
+
+2. **Partial vs. marginal.** @fig-partial-evs (partial effects) and @fig-effects-evs (marginal effects) both show the EVS slopes by cell. What does each one add that the other does not?
+
+::: {.callout-note collapse="true" title="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.
+:::
+
+3. **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?
+
+::: {.callout-note collapse="true" title="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 {{< git-rev short=true >}} using Quarto {{< version >}} and Julia {{< julia-version >}}.*
diff --git a/examples/eyevoicespan/pylsk13.rda b/examples/eyevoicespan/pylsk13.rda
new file mode 100755
index 0000000..ad39a45
--- /dev/null
+++ b/examples/eyevoicespan/pylsk13.rda
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e7d8ed826b031bcd62706999f9e994b8a963e41f03e9d8148e7f732fc77d62d5
+size 1549
diff --git a/examples/eyevoicespan/pylsk13_interactions.qmd b/examples/eyevoicespan/pylsk13_interactions.qmd
new file mode 100644
index 0000000..91447a8
--- /dev/null
+++ b/examples/eyevoicespan/pylsk13_interactions.qmd
@@ -0,0 +1,305 @@
+---
+title: "Interactions Three Ways: Observed, Partial, and Predicted (Pan et al., 2013)"
+engine: julia
+author:
+- Reinhold Kliegl
+- Jinger Pan
+julia:
+  exeflags: ["--project", "--threads=auto"]
+---
+
+::: {.objectives title="Learning objectives"}
+After working through this page you will be able to:
+
+- distinguish three complementary views of a model term — the **observed** relationship, the **partial-effect** view (`remef` / `MixedModelsExtras.partial_fitted`), and the **predicted** (marginal) view (`Effects.jl`);
+- read each of the significant interactions in `m09` off a consistent three-panel figure;
+- explain *why* an interaction that is invisible in the raw data can be obvious in the partial-effect plot.
+:::
+
+::: {.callout-note title="Before you start"}
+**Prerequisite:** [Eye–Voice Span in Rapid Automatized Naming](pylsk13.qmd) — 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/`).
+:::
+
+# Why three panels?
+
+Model `m09` from the [companion page](pylsk13.qmd) 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.
+
+# Setup
+
+```{julia}
+#| code-fold: true
+#| output: false
+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!()
+```
+
+```{julia}
+#| code-fold: true
+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)
+```
+
+The four interaction terms and their coefficients:
+
+```{julia}
+#| code-fold: true
+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
+```
+
+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 responses
+
+`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.
+
+```{julia}
+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)
+```
+
+# The three-panel helper
+
+```{julia}
+#| code-fold: true
+"""
+    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
+end
+```
+
+# (1) Condition × EVS
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "The eye–voice-span slope by naming Condition. Observed: both slopes negative. Partial: only `digit` retains a slope once Group and gaze are removed. Predicted: same pattern, averaged over Group (hence weaker) and with model uncertainty."
+#| label: fig-c-evs
+let
+  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
+```
+
+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.
+
+# (2) Group × EVS
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "The eye–voice-span slope by Group. The EVS benefit is concentrated in the control children."
+#| label: fig-g-evs
+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.
+
+# (3) EVS × Gaze
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "The gaze-duration slope at low, mid, and high eye–voice span (tertiles). Observed: the three slopes are nearly identical. Partial and Predicted: the gaze slope is steeper when the eye–voice span is short."
+#| label: fig-evs-gaze
+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.
+
+# (4) Condition × Group × EVS
+
+```{julia}
+#| code-fold: true
+#| fig-cap: "The three-way interaction (Figure 1 of Pan et al., 2013). The negative EVS slope is specific to the digit / control cell."
+#| label: fig-cg-evs
+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.
+
+# Takeaways
+
+- Read an interaction from the **partial** panel: it removes the confounds that hide
+  the term in raw data while keeping the scatter that tells you how well it is
+  determined.
+- Use the **predicted** panel for inference (confidence bands) and for communicating a
+  clean model summary — but remember it averages over whatever you left off the grid.
+- When the partial and predicted panels *disagree* (Figures 1–3) it is usually because
+  a collapsed moderator also interacts with the x-axis variable; when they *agree*
+  (Figure 4) the grid already contains everything that matters.
+
+# References
+
+::: {#refs}
+:::
+
+## Exercises
+
+1. **Keep the random intercept.** Recompute `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?
+
+::: {.callout-note collapse="true" title="Solution"}
+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`.
+:::
+
+2. **Effects grid resolution.** In Figure 3 the predicted lines are evaluated at three EVS values (the tertile medians). Replace them with five quantiles. Does the qualitative "fan" conclusion change? What is the cost?
+
+::: {.callout-note collapse="true" title="Solution"}
+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.
+:::
+
+3. **A non-significant term.** Build a three-panel figure for `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?
+
+::: {.callout-note collapse="true" title="Solution"}
+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 {{< git-rev short=true >}} using Quarto {{< version >}} and Julia {{< julia-version >}}.*
diff --git a/examples/eyevoicespan/pylsk13_v3.R b/examples/eyevoicespan/pylsk13_v3.R
new file mode 100644
index 0000000..25b46b9
--- /dev/null
+++ b/examples/eyevoicespan/pylsk13_v3.R
@@ -0,0 +1,304 @@
+# Predicting psychometric digit-RAN and dice-RAN with gaze duration and 
+# eye-voice span (EVS) from computerized assessment of digit-RAN and dice-RAN 
+# in Chinese control and dyslexic children. These results are reported in:
+
+# Pan et al. (2013). Eye–voice span during rapid automatized naming of digits  
+# and dice in Chinese normal and dyslexic children. Developmental Science. 
+
+# June 2013, Reinhold Kliegl & Jinger Pan
+# 7 & 14 December 2016, Reinhold Kliegl
+
+library(ggplot2)
+library(grid)
+library(reshape2)
+library(plyr)
+library(MASS)
+library(lme4)
+library("LMERConvenienceFunctions" ) 
+library(gtools)
+library(latticeExtra)
+
+library(RePsychLing)
+
+source("remef.v0.6.9.R")
+
+RKstats <- function(x) c(N=length(x), M=mean(x), SD=sd(x), SE=sd(x)/sqrt(length(x)) )
+
+vplayout <- function(x, y) {
+viewport(layout.pos.row = x, layout.pos.col = y) 
+}
+
+theme_set(theme_bw())
+
+# OVERVIEW 
+# Part 1: Setup
+# Part 2: Preliminary analyses
+# Part 3: Main LMM analyses: Using gaze and EVS as covariates
+# Part 4: Figures of paper -- including some new figures (3-group  visualisation)
+
+# PART 1: SETUP
+
+load("pylsk13.rda")
+
+data <- dat
+data$Subj <- factor(data$subj)
+data$Group <- factor(data$group, labels=c("control", "dyslexic"))
+data$Condition <- factor(data$condition, labels=c("digit", "dice"))
+
+
+# Check distributions of continuous variables
+boxcox(ran ~ subj*condition, data=data)
+boxcox(gaze ~ subj*condition, data=data)
+boxcox(evs ~ subj*condition, data=data)
+
+# Justifies log-transform of ran, not of gaze and evs
+data$lran <- log(data$ran)
+
+# Center covariates
+data$gaze.c <- scale(data$gaze,scale=FALSE,center=TRUE)
+data$evs.c <- scale(data$evs,scale=FALSE,center=TRUE)
+
+# PART 2: PRELIMINARY ANALYSES
+
+# Overall and within-group correlations
+data_w <- cbind(data[data$Condition=="digit", ], 
+                data[data$Condition=="dice", ])[ , c(9:10, 4, 19, 6, 21, 5, 20)]
+names(data_w) <- c("Subj", "Group", "g_ran", "c_ran", "g_evs", "c_evs", "g_gaze", "c_gaze")
+options(digits=2)
+
+#  Psychometric RAN (pmRAN)
+M_RAN <- ddply(data, c("Group", "Condition"), with, each(RKstats) (lran) )
+# Note: SEs only valid for between-subject comparisons (Group)
+
+ggplot(data=M_RAN, aes(x=Group, y=M, group=Condition, shape=Condition)) +
+       geom_line() + geom_point() +
+	     scale_y_continuous("log(RAN)") + 
+       geom_errorbar(aes(ymax=M+2*SE, ymin=M-2*SE, width=0.03)) +
+       coord_cartesian(ylim=c(2,4))
+
+r1  <- lmer(lran ~ group*condition + (1 | subj), REML=FALSE, data=data)
+print(summary(r1), cor=FALSE)
+
+r1b  <- lmer(lran ~ group*condition + (1 + condition || subj), REML=FALSE, data=data)
+print(summary(r1b), cor=FALSE)
+summary(rePCA(r1b))
+
+anova(r1, r1b) # not significant
+
+## max model not possible
+#r1c  <- lmer(lran ~ group*condition + (1 + condition | subj), REML=FALSE, data=data,
+#             control=lmerControl(check.nobs.vs.nRE = "ignore"))
+#print(summary(r1c), cor=FALSE)
+
+
+# Gaze -- replicates pattern for pmRAN
+M_gaze <- ddply(data, c("Group", "Condition"),  with, each(RKstats) (gaze) )
+# Note: SEs only valid for between-subject comparisons (Group)
+
+qplot(data=M_gaze, x=Group, y=M, group=Condition, 
+      shape=Condition, geom=c("line", "point"), ylab="Gaze duration (ms)") + 
+      geom_errorbar(aes(max=M+2*SE, min=M-2*SE, width=0.03)) 
+
+g1 <- lmer(gaze ~ group*condition + (1 | subj), REML=FALSE, data=data)
+print(summary(g1), cor=FALSE)
+
+g1b <- lmer(gaze ~ group*condition + (1 + condition || subj), REML=FALSE, data=data)
+print(summary(g1b), cor=FALSE)
+summary(rePCA(g1b))
+
+anova(g1, g1b)
+
+# So if we include gaze as covariate for pmRAN
+r1 <- lmer(lran ~ group*condition + (1 | subj), data=data, REML=FALSE)
+r2 <- lmer(lran ~ group*condition+gaze + (1 | subj), data=data, REML=FALSE)
+print(summary(r2), cor=FALSE)
+
+r2a <- lmer(lran ~ (group+condition+gaze.c)^2 + (1 | subj), data=data, REML=FALSE)
+print(summary(r2a), cor=FALSE)
+
+r2b <- lmer(lran ~ (group+condition+gaze.c)^3 + (1 | subj), data=data, REML=FALSE)
+print(summary(r2b), cor=FALSE)
+
+anova(r1, r2, r2a, r2b) 
+
+# Add variance components for within-subject effects
+r3a <- lmer(lran ~ group*condition+gaze.c + (1 + condition || subj), data=data, REML=FALSE)
+print(summary(r3a), cor=FALSE)
+summary(rePCA(r3a))
+anova(r2, r3a)
+
+r3b <- lmer(lran ~ group*condition+gaze.c + (1 + gaze.c || subj), data=data, REML=FALSE)
+print(summary(r3b), cor=FALSE)
+summary(rePCA(r3b))
+anova(r2, r3b)
+
+# r2 looks like best model
+
+# Three main effects (group, condition, gaze.c), plus significant group * condition interaction
+r2 <- lmer(lran ~ group*condition+gaze.c + (1 | subj), REML=FALSE, data=data)
+print(summary(r2))
+
+# EVS -- replicates pattern for pmRAN
+M_EVS <- ddply(data, c("Group", "Condition"), with, each(RKstats) (evs))
+# Note: SEs only valid for between-subject comparisons (Group)
+
+qplot(data=M_EVS, x=Group, y=M, group=Condition, shape=Condition, geom=c("line", "point"), ylab="Eye-voice span (char)") + 
+      geom_errorbar(aes(max=M+1*SE, min=M-1*SE, width=0.03)) 
+
+e1  <- lmer(evs ~ group*condition + (1 | subj), REML=FALSE, data=data)
+print(summary(e1), cor=FALSE)
+
+# So if we include evs as covariate for pmRAN,
+data$evs.c <- scale(data$evs,scale=FALSE,center=TRUE)
+
+r3 <- lmer(lran ~ group*condition+evs.c + (1 | subj), data=data, REML=FALSE)
+r3a <- lmer(lran ~ (group+condition+evs.c)^2 + (1 | subj), data=data, REML=FALSE)
+r3b <- lmer(lran ~ (group+condition+evs.c)^3 + (1 | subj), data=data, REML=FALSE)
+
+anova(r1, r3, r3a, r3b)  # evs.c is significant covariate
+# Three main effects (group, condition, evs.c), plus significant condition*group interaction
+
+# PART 3: Main LMM analyses: Using gaze and EVS as covariates
+
+m00 <- lmer(lran ~ condition+group+evs.c+gaze.c + (1  | subj), data=data, REML=FALSE)
+print(summary(m00), cor=FALSE)
+
+# Minimum model given previous results
+m05 <- lmer(lran ~ condition*group+evs.c+gaze.c + (1  | subj), data=data, REML=FALSE)
+print(summary(m05), cor=FALSE)
+
+# Full factorial 
+m10 <- lmer(lran ~ condition*group*evs.c*gaze.c + (1  | subj), data=data, REML=FALSE)
+print(summary(m10), cor=FALSE)
+
+# Remove 3 non-significant higher-order interactions involving evs.c:gaze.c
+m09 <- lmer(lran ~ condition*group*evs.c*gaze.c - condition:group:evs.c:gaze.c 
+	                     - group:evs.c:gaze.c - condition:gaze.c:evs.c
+	                     + (1  | subj), data=data, REML=FALSE)
+print(summary(m09), corr=FALSE)
+anova(m00, m05, m09, m10)
+# Best model; reported in Table 3
+
+# Check effect of log-transformation
+m09b <- lmer(ran ~ condition*group*evs.c*gaze.c - condition:group:evs.c:gaze.c 
+            - group:evs.c:gaze.c - condition:gaze.c:evs.c
+            + (1  | subj), data=data, REML=FALSE)
+print(summary(m09b), corr=FALSE)
+
+# Check residuals
+qqmath(resid(m09))
+
+qplot(x=fitted(m09), y=resid(m09), geom="point", 	
+	xlab="Fitted values", ylab="Standardized residuals") + 
+	geom_hline(yintercept=0) + 
+	geom_density2d(size=1) 
+
+# ... for log-transformed values
+qqmath(resid(m09b))      
+qplot(x=fitted(m09b), y=resid(m09), geom="point", 	
+      xlab="Fitted values", ylab="Standardized residuals") + 
+  geom_hline(yintercept=0) + 
+  geom_density2d(size=1) 
+
+# PART 4: Figures -- unadjusted and partial plots (model m09) of interactions
+
+
+# Figure 1: Condition x  Group x EVS (t=-2.5)
+
+# Unadjusted observed scores
+p1 <- qplot(data=data, y=lran, x=evs,  group=Condition:Group, color=Condition:Group, 
+            geom=c("point", "smooth"), method=lm,
+            xlab = "Eye-voice span", ylab="log(RAN)", ylim=c(2, 4)) + 
+  theme(legend.position = "none", panel.background=element_rect(fill = "white"))
+p1
+
+# -- evs for dice better predictor for dyslexic, 
+# -- evs for digits better predictor for control
+
+# Partial effects
+data$CndGrpEVS.m09 <- remef(m09, keep=TRUE, grouping=TRUE, fix = c(1, "condition:group:evs.c"), 
+                            ran = NULL, plot=FALSE) 
+
+p2 <- 
+  ggplot(data=data, aes(y=CndGrpEVS.m09, x=evs, 
+                        group=Condition:Group, color=Condition:Group)) +
+  geom_point() + 
+  geom_smooth(method="lm") +
+  scale_x_continuous("Eye-voice span") +
+  scale_y_continuous("Adjusted log(RAN)") + 
+  coord_cartesian(ylim=c(2,4)) +
+  theme(legend.position = c(.30, .20), 
+        legend.text = element_text(size=8), 
+        panel.background=element_rect(fill = "white")) 
+p2
+
+#  -- evs only predictor for control in digit-RAN!
+
+grid.newpage() # Figure 1
+pushViewport(viewport(layout = grid.layout(1,2))) 
+print(p1, vp=vplayout(1,1))
+print(p2, vp=vplayout(1,2))
+
+
+# Figure 2: EVSgroup x GAZE 
+
+# Partial effect for Condition x Group x GD
+data$CndGrpGAZE.m09 <- 
+  remef(m09, keep=TRUE, grouping=TRUE, fix = c(1, "condition:group:gaze.c"), ran = NULL) 
+
+p3 <- 
+  ggplot(data=data, aes(y=CndGrpGAZE.m09, x=gaze, 
+                        group=Condition:Group, color=Condition:Group)) + 
+  geom_point() + 
+  geom_smooth(method="lm") +
+  scale_x_continuous("Gaze duration") +
+  scale_y_continuous("Adjusted log(RAN)") + 
+  coord_cartesian(ylim=c(2,4)) +
+  theme(legend.position = c(.75, .25), panel.background=element_rect(fill = "white")) 
+p3
+
+# -- dice control slope looks different, but this difference is not strong enough?
+
+# Partial effect for EVS x GD, visualized for small and large EVS groups
+
+# ... form two EVS groups
+idEVS <- ddply(data, "Subj", with, each(M=mean)(evs))
+idEVS$EVSgroup.2 <- quantcut(idEVS$M, seq(0, 1, by=1/2), label=c("small", "large"))
+data1 <- merge(data, idEVS, by="Subj")
+
+data1$EVS_GAZE.m09  <- remef(m09, keep=TRUE, grouping=TRUE, fix = c(1, "evs.c:gaze.c"), ran = NULL, plot=FALSE)
+
+p4 <- qplot(data=data1, y=EVS_GAZE.m09, x=gaze, group=EVSgroup.2, color=EVSgroup.2, geom=c("point", "smooth"), method=lm, 
+		xlab = "Gaze duration", ylab="Adjusted log(RAN)", ylim=c(2, 4)) + scale_colour_hue("EVS group") +
+	theme(legend.position = c(.75, .25),panel.background=element_rect(fill = "white")) 
+
+p4b <- qplot(data=data1, y=lran, x=gaze, group=EVSgroup.2, color=EVSgroup.2, geom=c("point", "smooth"), method=lm, 
+            xlab = "Gaze duration", ylab="log(RAN)", ylim=c(2, 4)) + scale_colour_hue("EVS group") +
+  theme(legend.position = c(.75, .25),panel.background=element_rect(fill = "white")) 
+
+grid.newpage()  #  Figure 2
+pushViewport(viewport(layout = grid.layout(1,2))) 
+print(p4b, vp=vplayout(1,1))
+print(p4, vp=vplayout(1,2))
+
+# Partial effect for EVS x GD, visualized for small, medium, and large EVS groups
+
+# ... form three EVS groups
+idEVS2 <- ddply(data, "Subj", with, each(M=mean)(evs))
+idEVS2$EVSgroup.3 <- quantcut(idEVS2$M, seq(0, 1, by=1/3), label=c("small", "medium", "large"))
+data2 <- merge(data1, idEVS2, by="Subj")
+
+# ... ... partial efffects
+p5 <- qplot(data=data2, y=EVS_GAZE.m09, x=gaze, group=EVSgroup.3, color=EVSgroup.3, geom=c("point", "smooth"), method=lm, 
+            xlab = "Gaze duration", ylab="Adjusted log(RAN)", ylim=c(2, 4)) + scale_colour_hue("EVS group") +
+  theme(legend.position = c(.75, .25),panel.background=element_rect(fill = "white")) 
+
+# ... ... zero-order relations
+p5b <- qplot(data=data2, y=lran, x=gaze, group=EVSgroup.3, color=EVSgroup.3, geom=c("point", "smooth"), method=lm, 
+             xlab = "Gaze duration", ylab="log(RAN)", ylim=c(2, 4)) + scale_colour_hue("EVS group") +
+  theme(legend.position = c(.75, .25),panel.background=element_rect(fill = "white")) 
+
+grid.newpage()  # New Figure (Dec 2016)
+pushViewport(viewport(layout = grid.layout(1,2))) 
+print(p5b, vp=vplayout(1,1))
+print(p5, vp=vplayout(1,2))
diff --git a/examples/eyevoicespan/pylsk13_v4.R b/examples/eyevoicespan/pylsk13_v4.R
new file mode 100644
index 0000000..a00b49d
--- /dev/null
+++ b/examples/eyevoicespan/pylsk13_v4.R
@@ -0,0 +1,187 @@
+# Predicting psychometric digit-RAN and dice-RAN with gaze duration and 
+# eye-voice span (EVS) from computerized assessment of digit-RAN and dice-RAN 
+# in Chinese control and dyslexic children. These results are reported in:
+
+# Pan et al. (2013). Eye–voice span during rapid automatized naming of digits  
+# and dice in Chinese normal and dyslexic children. Developmental Science. 
+
+# June 2013, Reinhold Kliegl & Jinger Pan
+# 7 & 14 December 2016, Reinhold Kliegl
+
+library(ggplot2)
+library(grid)
+library(reshape2)
+library(plyr)
+library(MASS)
+library(lme4)
+library("LMERConvenienceFunctions" ) 
+library(gtools)
+library(latticeExtra)
+
+library(RePsychLing)
+
+source("remef.v0.6.9.R")
+
+RKstats <- function(x) c(N=length(x), M=mean(x), SD=sd(x), SE=sd(x)/sqrt(length(x)) )
+
+vplayout <- function(x, y) {
+viewport(layout.pos.row = x, layout.pos.col = y) 
+}
+
+theme_set(theme_bw())
+
+# OVERVIEW 
+# Part 1: Setup
+# Part 2: Preliminary analyses
+# Part 3: Main LMM analyses: Using gaze and EVS as covariates
+# Part 4: Figures of paper -- including some new figures (3-group  visualisation)
+
+# PART 1: SETUP
+
+load("pylsk13.rda")
+
+data <- dat
+data$Subj <- factor(data$subj)
+data$Group <- factor(data$group, labels=c("control", "dyslexic"))
+data$Condition <- factor(data$condition, labels=c("digit", "dice"))
+
+
+# Check distributions of continuous variables
+boxcox(ran ~ subj*condition, data=data)
+boxcox(gaze ~ subj*condition, data=data)
+boxcox(evs ~ subj*condition, data=data)
+
+# Justifies log-transform of ran, not of gaze and evs
+data$lran <- log(data$ran)
+
+# Center covariates
+data$gaze_c <- scale(data$gaze,scale=FALSE,center=TRUE)
+data$evs_c <- scale(data$evs,scale=FALSE,center=TRUE)
+
+# PART 2: Preliminary analysis: <Skipped see other script>
+
+# PART 3: Main LMM analyses: Using gaze and EVS as covariates
+
+m00 <- lmer(lran ~ condition+group+evs_c+gaze_c + (1  | subj), data=data, REML=FALSE)
+print(summary(m00), cor=FALSE)
+
+# Minimum model given previous results
+m05 <- lmer(lran ~ condition*group+evs_c+gaze_c + (1  | subj), data=data, REML=FALSE)
+print(summary(m05), cor=FALSE)
+
+# Full factorial 
+m10 <- lmer(lran ~ condition*group*evs_c*gaze_c + (1  | subj), data=data, REML=FALSE)
+print(summary(m10), cor=FALSE)
+
+# Remove 3 non-significant higher-order interactions involving evs_c:gaze_c
+m09 <- lmer(lran ~ condition*group*evs_c*gaze_c - condition:group:evs_c:gaze_c 
+	                     - group:evs_c:gaze_c - condition:gaze_c:evs_c
+	                     + (1  | subj), data=data, REML=FALSE)
+print(summary(m09), corr=FALSE)
+anova(m00, m05, m09, m10)
+# Best model; reported in Table 3
+
+# Check effect of log-transformation
+m09b <- lmer(ran ~ condition*group*evs_c*gaze_c - condition:group:evs_c:gaze_c 
+            - group:evs_c:gaze_c - condition:gaze_c:evs_c
+            + (1  | subj), data=data, REML=FALSE)
+print(summary(m09b), corr=FALSE)
+
+# Check residuals
+qqmath(resid(m09))
+
+qplot(x=fitted(m09), y=resid(m09), geom="point", 	
+	xlab="Fitted values", ylab="Standardized residuals") + 
+	geom_hline(yintercept=0) + 
+	geom_density2d(size=1) 
+
+# ... for log-transformed values
+qqmath(resid(m09b))      
+qplot(x=fitted(m09b), y=resid(m09), geom="point", 	
+      xlab="Fitted values", ylab="Standardized residuals") + 
+  geom_hline(yintercept=0) + 
+  geom_density2d(linewidth=1) 
+
+# PART 4: Figures -- unadjusted and partial plots (model m09) of interactions
+
+
+# Figure 1: Condition x  Group x EVS (t=-2.5)
+
+## Unadjusted observed scores
+p1 <- 
+  ggplot(data=data, aes(y=lran, x=evs,  group=Condition:Group, color=Condition:Group)) +
+         geom_point() + geom_smooth(method="lm") +
+         scale_x_continuous("Eye-voice span") +
+         scale_y_continuous("log(RAN)") + 
+  coord_cartesian(ylim=c(2,4)) +
+  theme(legend.position = "none", 
+        panel.background=element_rect(fill = "white"))
+p1
+
+# -- evs for dice better predictor for dyslexic, 
+# -- evs for digits better predictor for control
+
+## Partial effects
+data$CndGrpEVS.m09 <- 
+  remef(m09, keep=TRUE, grouping=TRUE, fix = c(1, "condition:group:evs_c"), 
+        ran = NULL, plot=FALSE) 
+
+p2 <- 
+  ggplot(data=data, aes(y=CndGrpEVS.m09, x=evs, 
+                          group=Condition:Group, color=Condition:Group)) +
+  geom_point() + 
+  geom_smooth(method="lm") +
+  scale_x_continuous("Eye-voice span") +
+  scale_y_continuous("Adjusted log(RAN)") + 
+  coord_cartesian(ylim=c(2,4)) +
+  theme(legend.position = c(.30, .20), 
+        legend.text = element_text(size=8), 
+        panel.background=element_rect(fill = "white")) 
+p2
+
+## Combine in two-panel figure
+
+grid.newpage() # Figure 1
+pushViewport(viewport(layout = grid.layout(1,2))) 
+print(p1, vp=vplayout(1,1))
+print(p2, vp=vplayout(1,2))
+
+#-- evs only predictor for control in digit-RAN!
+  
+# Figure 2: EVSgroup x GAZE 
+
+## Unadjusted observed scores
+p3 <- 
+  ggplot(data=data, aes(y=lran, x=gaze, 
+                        group=Condition:Group, color=Condition:Group)) + 
+  geom_point() + 
+  geom_smooth(method="lm") +
+  scale_x_continuous("Gaze duration") +
+  scale_y_continuous("log(RAN)") + 
+  coord_cartesian(ylim=c(2,4)) +
+  theme(legend.position = "none") 
+p3
+
+## Partial effect for Condition x Group x GD
+data$CndGrpGAZE.m09 <- 
+  remef(m09, keep=TRUE, grouping=TRUE, fix = c(1, "condition:group:gaze_c"), ran = NULL) 
+
+p4 <- 
+  ggplot(data=data, aes(y=CndGrpGAZE.m09, x=gaze, 
+                        group=Condition:Group, color=Condition:Group)) + 
+  geom_point() + 
+  geom_smooth(method="lm") +
+  scale_x_continuous("Gaze duration") +
+  scale_y_continuous("Adjusted log(RAN)") + 
+  coord_cartesian(ylim=c(2,4)) +
+	theme(legend.position = c(.75, .25), 
+	      panel.background=element_rect(fill = "white")) 
+p4
+
+## Combine in two-panel figure
+
+grid.newpage() # Figure 2
+pushViewport(viewport(layout = grid.layout(1,2))) 
+print(p3, vp=vplayout(1,1))
+print(p4, vp=vplayout(1,2))
+
diff --git a/examples/eyevoicespan/remef.v0.6.9.R b/examples/eyevoicespan/remef.v0.6.9.R
new file mode 100644
index 0000000..b79baa1
--- /dev/null
+++ b/examples/eyevoicespan/remef.v0.6.9.R
@@ -0,0 +1,226 @@
+#############################################################
+### function remef (REMove EFfects)
+# remove random factors variance and fixed effects from the dependent variable of an LMM analysis
+
+# by Sven Hohenstein, Reinhold Kliegl, 2011, 2012, 2013
+
+# v0.6.9, July 2013
+
+## VERSION HISTORY:
+# changes to last version (v0.6.8, June 2013):
+# - function now works with glmerMod models created with glmer
+# v0.6.7, December 2012:
+# - function now uses the new lme4 function for 'lmerMod' objects
+#   (lme4Eigen is no longer used)
+# - the plot parameter is deprecated
+# - minor code changes
+# v0.6.6, February 2012:
+# - the argument "all" (string) for the parameter 'ran'
+#   selects all random effects
+# v0.6.5, February 2012:
+# - the function now works with objects created by the library
+#  	lme4Eigen too
+# changes to last version (v0.6.2, August 2011):
+# - former parameter 'link' is now called 'family'
+# - family parameter is specified as object (not string),
+#	like, e.g., in lmer()
+# - functions are called from the correct package
+# - new parameter 'plot': if TRUE, a plot of the uncorrected
+# 	and corrected dependent variable is displayed
+# v0.6.2 August 2011
+# - a bug occuring when NULL was part of the random effects list
+#   was removed
+# v0.6, June 2011
+# - a bug was removed
+# v0.54, June 2011:
+# - new parameter 'link' specifying the logit link function; either
+#   "identity" or "logit"
+# - redundant entries are removed from the vectors in the list 'ran'
+# v0.54, June 2011:
+# - Fixed effects (in the vector 'fix') can be integers or strings
+#	(e.g., c(2:4, "Effect1", 6, "Effect2") )
+# v0.52, May 2011
+# - Correction: If grouping was TRUE and a specified fixed main effect of a
+#     factor was not present in any interaction, the function failed
+# v0.51, May 2011:
+# - parameter 'grouping' also works if keep = FALSE, if grouping is TRUE,
+#	the effect and all adssociated ones of higher order are chosen to be
+#	removed
+# v0.5, May 2011:
+# - new boolean parameter: 'grouping'; if both keep and grouping are TRUE,
+#	the effects in fix and all subeffects are chosen to be kept
+#	Note: In the present version, if keep is FALSE, grouping will be ignored
+# v0.4a, May 2011:
+# - if keep is TRUE, the to be kept effects are not removed from the
+#	dependent variable but are actually added to the residuals
+# - more efficient processing and less operations (if keep is set properly)
+# v0.4, April 2011:
+# - new parameter 'keep' allows the remove the effects not specified
+# - more efficient processing
+# v0.32, April 2011:
+# - Correction: If random effects of random factors with order numbers < 1
+#				were specified, the output was not correct.
+# v0.3a, April 2011:
+# changes to last version (v0.3):
+# - more efficient processing
+# v0.3, April 2011:
+# changes to last version (v0.2):
+# - now it is possible to remove all random effects (even random slopes)
+# - the parameter 'ran' must be a list of vectors
+# - switched order position of parameters 'fix' and 'ran'
+# v0.2, February 2011:
+# changes to last version (v0.1):
+# - Correction: If just one fixed effect was specified, the function failed.
+#				Now it is possible to specify a single number as parameter 'fix'
+
+
+## INPUT:
+# model:	an object of class 'mer' (lme4 packge), 'lmerMod' (lme4Eigen package)
+# ran:		list of vectors of natural numbers according to random effects of the model;
+#			the maximum length of this list is the number of random factors;
+#			each list member includes the numbers of the random effects (of the current random factor)
+#			(e.g., list(1:2, NULL, c(1, 3)) )
+#			default value: NULL (no random-factor related variance is removed)
+# fix:		vector of natural numbers or strings according to fixed effects of the model
+#			(e.g., c(2:4, 6, 8:10), c("Effect1", "Effect2"), or c(2:4, "Effect1", 6, "Effect2"))
+#			Vectors can consist of both integers and strings (R will transform all intergers to
+#			strings if at least one string is entered);
+#			redundant values will be ignored (one effect can't be removed more than once)
+#			default value: NULL (no fixed effects are removed)
+# keep:		logical value; if TRUE, the specified effects are not removed but kept
+#			(and all other effects are removed)
+# grouping:	logical value; if FALSE the specified effects in fix are used;
+#			if TRUE, effects are grouped:
+#				- if keep=TRUE, the specified effect (and all effects associated with the
+#					variables) as well as all effects of lower order consisting solitary of
+#					a subset of the variables of the specified effet are kept
+#					(e.g., if the interaction A:B:C is specified, the following effects are
+#					chosen:	A:B:C, A:B, A:C, B:C, A, B, and C)
+#				- if keep=FALSE, the specified effect (and all effects associated with the
+#					variables) as well as all effects of higher order consisting of
+#					all variables of the specified effet (and other variables) are removed
+#					(e.g., if the interaction A:B:C is specified, the following effects are
+#					chosen:	A:B:C and A:B:C:x)
+# family: a family object including a link function; 
+#         either "identity" (default) or "logit" link function; this must be the
+#         link function which was used for generating the model;
+#       - gaussian(link = "identity") is the right option for most purposes, e.g. Gaussian
+#         distributions; the dependent variable is not transformed
+#       - binomial(link = "logit") is used for binomial dependent variables; the output vector
+#         will contain probabilities which can be rounded to obtain data in the
+#         original metric (zeroes and ones); note that even if keep=FALSE, the
+#         base for the calculation are not the input values but the residuals
+
+## OUTPUT:	a numerical vector
+
+## How does it work?
+# Random factor variance and fixed effects are removed from the dependent variable of an LMM analysis.
+# The order of effects in the fixed-effects input vector is the same as in the model output.
+# The returned vector includes the dependent variable corrected by the specified effects.
+# Note. The currect version of this function works for the "identity" and "logit" model link functions
+# only.
+
+remef <- function(model, fix = NULL, ran = NULL, keep = FALSE, grouping = FALSE, family = gaussian(link = "identity"), plot = FALSE) {
+  # name of link function
+  mclass <- class(model)
+  sup_classes <- c("mer", "lmerMod", "glmerMod")	# supported object classes
+  if (!mclass %in% c(sup_classes)) stop("This class is not supported yet.")
+  fixef.fnc <- lme4::fixef
+  ranef.fnc <- lme4::ranef
+  moma <- model.matrix(model)
+  	# model matrix
+  if (is.function(family)) family <- family()
+  link <- family$link
+  if (!(link %in% c("identity", "logit"))) stop("Unknown link function specified.")
+  if (keep || link == "logit") {
+  	#DV <- lme4::residuals(model)			# use residuals as base and add effects
+  	DV <- residuals(model)
+  } else {
+  	# use actual data as base and remove effects (this is not possible with the logit link function)
+  	DV <- model@frame[ , 1]
+  }	
+  fix <- unique(fix)				# remove redundant values
+  if (any(is.character(fix))) {		# strings were entered as fixed effects
+  	suppressWarnings( fix.num <- as.numeric(fix) )
+  		# fix.num, numerical fix vector
+  	ef.names.idx <- is.na(fix.num)
+  		# logical vector indicating the positions of effect names in the fix vector
+  	for (ef in fix[ef.names.idx]) {
+  		if(!any(ef == names(fixef.fnc(model)))) stop(paste("The fixed effect", ef, "is not present in the model."))
+  			# stops the function if any fixed effect string is not present in the model
+  	}	
+  	fix.num[ef.names.idx] <- which( names(fixef.fnc(model)) %in% fix[ef.names.idx] )
+  	fix <- unique(fix.num)
+  }
+  if (grouping & length(fix) > 0) {			
+  	new.fix <- NULL				# help variable for the construction of a new 'fix' vector
+  	for (ef in fix) {			# choose the specified effects and all lower-order/higher-order effects of the relevant variables
+  		# var.str <- attr(lme4::model.matrix(model), "assign")
+  		var.str <- attr(moma, "assign")
+  			# var.str, structure of the variables
+  		if (!var.str[ef]) {		# the specified effect is the intercept
+  			new.fix <- unique(c(new.fix, ef))
+  		} else {				# the specified effect is NOT the intercept
+  			var.mat <- attr(terms(model), "factors")
+  				# var.mat, variable matrix
+  			ef.order <- attr(terms(model), "order")
+  				# ef.order, order of effects (0: intercept, 1: main effect, 2: two-factor interaction, etc.)
+  			if (keep) {			# choose specified effect an all associated ones of lower order
+  				ass.ef <- which( var.str %in% which( colSums( matrix( var.mat[ var.mat[ , var.str[ef] ] == 1, ] , ncol = length(ef.order)) ) == ef.order ) )
+  					# ass.ef, effects associated with the input effect (integer(s))
+  			} else {			# choose specified effect an all associated ones of higher order
+  				ass.ef <- which( var.str %in% which( colSums( matrix( var.mat[ var.mat[ , var.str[ef] ] == 1, ] , ncol = length(ef.order)) ) == ef.order[var.str[ef]] ) )
+  			}	
+  			new.fix <- unique(c(new.fix, ass.ef))
+  		}  		
+  	}	  
+  fix <- new.fix				  		
+  }
+  if (!keep && link == "logit") { # keep effects!
+  	fix <- setdiff(seq_along(fixef.fnc(model)), fix)
+  }  	
+  # remove random factor variance
+  if (identical(ran, "all")) {
+    ran <- lapply(ranef(model), seq_along)    
+  }
+  if (length(ran) > 0) {
+  	rf_before.end.at <- -1
+	# for (rf in 1 : length(ran)) {			
+    for (rf in seq_along(ranef.fnc(model))) {	# rf, random factor
+  		if (rf > length(ran)) {				# non-specified random factors
+  			if (!keep && link == "logit") { # keep effects!
+  				ran[[rf]] <- seq.int(ncol(ranef.fnc(model)[[rf]]))
+  			} else {
+  				ran[[rf]] <- NULL
+  			}		
+  		} else {  			
+  			if(!is.null(ran[[rf]])) ran[[rf]] <- unique(ran[[rf]])
+  			if (!keep && link == "logit") { # keep effects!
+  				ran[[rf]] <- setdiff(seq.int(ncol(ranef.fnc(model)[[rf]])), ran[[rf]])
+  			} 
+  		}
+  		n.rf.levels	<- nrow(ranef.fnc(model)[[rf]])
+			# n.rf.levels, number of random-factor levels
+      	for (re in ran[[rf]]) {				# re, random effect (e.g., intercept, slope)
+        	idx.re <- ( ((re - 1) * n.rf.levels + 1) : (re * n.rf.levels) ) + (rf_before.end.at + 1)
+          		#idx.re, numbers of lines of the (transpose) random effect matrix for the currect random effect
+        	if (mclass == "lmerMod") re.matrix <- model@pp$Zt else re.matrix <- model@Zt        	  
+        	DV <- DV + sign((keep || link == "logit") - 0.5) * ( as.vector( ranef.fnc(model)[[rf]][ , re] %*% re.matrix[idx.re, ] ) )
+          		# Note. re.matrix is the transpose sparse random-effect matrix (class dgCMatrix)
+				# Note. sign(keep - 0.5) is +1 if keep==TRUE and -1 otherwise
+		}
+      rf_before.end.at <- rf_before.end.at + prod(dim(ranef.fnc(model)[[rf]]))
+        # at which line of the random-effect matrix did the current random effect end?
+  	}
+  }
+  # remove fixed effects  
+  DV <- DV + sign((keep || link == "logit") - 0.5) * ( matrix(moma[ , fix], nrow = length(DV)) %*% fixef.fnc(model)[fix] )
+  if (link == "logit") DV <- ( 1 / (1 + exp(- DV)) )
+  	# if DV is rounded, zeros and ones will result
+  if (plot) {
+    warning("The 'plot' parameter is deprecated.")
+  }	
+  return(as.vector(DV))	
+}
+
+#############################################################
diff --git a/examples/visual-attention/kkl15.qmd b/examples/visual-attention/kkl15.qmd
index ba7bc5f..d07d4d3 100644
--- a/examples/visual-attention/kkl15.qmd
+++ b/examples/visual-attention/kkl15.qmd
@@ -32,7 +32,7 @@ We specify three contrasts for the four-level factor CTR that are derived from s
 
 This comparison is of interest because a few years after the publication of @Kliegl2011, the theoretically critical correlation parameter (CP) between the spatial effect and the attraction effect was determined as the source of a non-singular LMM in that paper. The present study served the purpose to estimate this parameter with a larger sample and a wider variety of experimental conditions.
 
-Here we also include two additional experimental manipulations of target size and orientation of cue rectangle. A similar analysis was reported in the parsimonious mixed-model paper [@Bates2015]; it was also used in a paper of GAMEMs [@Baayen2017]. Data and R scripts of those analyses are also available in [R-package RePsychLing](https://github.com/dmbates/RePsychLing/tree/master/data/).
+Here we also include two additional experimental manipulations of target size and orientation of cue rectangle. A similar analysis was reported in the parsimonious mixed-model paper [@Bates2015]; it was also used in a paper of GAMMs [@Baayen2017]. Data and R scripts of those analyses are also available in [R-package RePsychLing](https://github.com/dmbates/RePsychLing/tree/master/data/).
 
 The analysis is based on reaction times `rt` to maintain compatibility with @Kliegl2011.
 
@@ -40,7 +40,7 @@ In this vignette we focus on the reduction of model complexity. And we start wit
 
 “Neither the [maximal] nor the [minimal] linear mixed models are appropriate for most repeated measures analysis. Using the [maximal] model is generally wasteful and costly in terms of statiscal power for testing hypotheses. On the other hand, the [minimal] model fails to account for nontrivial correlation among repeated measurements. This results in inflated [T]ype I error rates when non-negligible correlation does in fact exist. We can usually find middle ground, a covariance model that adequately accounts for correlation but is more parsimonious than the [maximal] model. Doing so allows us full control over [T]ype I error rates without needlessly sacrificing power.”
 
-Stroup, W. W. (2012, p. 185). _Generalized linear mixed models: Modern concepts, methods and applica?ons._ CRC Press, Boca Raton.
+Stroup, W. W. (2012, p. 185). _Generalized linear mixed models: Modern concepts, methods and applications._ CRC Press, Boca Raton.
 
 # Packages
 
@@ -68,6 +68,8 @@ progress = isinteractive()
 ```{julia}
 #| echo: false
 using SMLP2026: fit_or_restore, bootstrap_or_restore
+
+using SMLP2026.RKUtils
 ```
 
 # Read data, compute and plot means
diff --git a/examples/visual-attention/kwdyz11.qmd b/examples/visual-attention/kwdyz11.qmd
index 735a9c2..7181486 100644
--- a/examples/visual-attention/kwdyz11.qmd
+++ b/examples/visual-attention/kwdyz11.qmd
@@ -59,8 +59,11 @@ using MixedModels
 using MixedModelsMakie
 using Random
 using MixedModelsDatasets: dataset
-using SMLP2026: fit_or_restore
+
 using StatsBase
+
+using SMLP2026: fit_or_restore, bootstrap_or_restore
+using SMLP2026.RKUtils
 ```
 
 # Read data, compute and plot densities and means
@@ -163,7 +166,7 @@ contrasts = Dict(
 )
 m1 = let
   form = @formula(log(rt) ~ 1 + CTR + (1 + CTR | Subj))
-  fit_or_restore("kwdyz11_m1.json", MixedModel, form, dat; contrasts)
+  fit_or_restore("kwdyz11_m1.json", MixedModel, form, dat; contrasts, REML=true)
 end
 ```
 
diff --git a/fits/eyevoicespan_m00.json.zip b/fits/eyevoicespan_m00.json.zip
new file mode 100644
index 0000000..b0453cd
--- /dev/null
+++ b/fits/eyevoicespan_m00.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:02618ff85db9afa15e436cc41dbd2e2babeacdb61904993d15b89003fcdb61ab
+size 646
diff --git a/fits/eyevoicespan_m05.json.zip b/fits/eyevoicespan_m05.json.zip
new file mode 100644
index 0000000..fa162b6
--- /dev/null
+++ b/fits/eyevoicespan_m05.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a49a374ece734172ebd0650b4b93ef87604ebb38121c058cbb1f4d4756d1773e
+size 640
diff --git a/fits/eyevoicespan_m09.json.zip b/fits/eyevoicespan_m09.json.zip
new file mode 100644
index 0000000..ea0c5ce
--- /dev/null
+++ b/fits/eyevoicespan_m09.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a9f4d85839984d1b6d6f83fe280998f8a962718b43ad3883bd816c92fe7dd5b6
+size 664
diff --git a/fits/eyevoicespan_m10.json.zip b/fits/eyevoicespan_m10.json.zip
new file mode 100644
index 0000000..db4a59f
--- /dev/null
+++ b/fits/eyevoicespan_m10.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4ddbbf0840639d0e3b8e6ffa969028522df983eea9dd9beb2eb0e1d809fb86ff
+size 616
diff --git a/references.bib b/references.bib
index 293b87d..1fe078d 100644
--- a/references.bib
+++ b/references.bib
@@ -80,6 +80,18 @@ @Article{Kliegl2011
   publisher = {Frontiers Media {SA}},
 }
 
+@Article{Pan2013,
+  author    = {Jinger Pan and Ming Yan and Jochen Laubrock and Hua Shu and Reinhold Kliegl},
+  journal   = {Developmental Science},
+  title     = {Eye--voice span during rapid automatized naming of digits and dice in {Chinese} normal and dyslexic children},
+  year      = {2013},
+  volume    = {16},
+  number    = {6},
+  pages     = {967--979},
+  doi       = {10.1111/desc.12075},
+  publisher = {Wiley},
+}
+
 @Misc{Bates2015,
   author    = {Bates, Douglas and Kliegl, Reinhold and Vasishth, Shravan and Baayen, Harald},
   title     = {Parsimonious Mixed Models},