diff --git a/Project.toml b/Project.toml
index 968249c..0c1d3ac 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,6 +22,7 @@ GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
+MixedModelsExtras = "781a26e1-49f4-409a-8f4c-c3159d78c17e"
 MixedModelsMakie = "b12ae82c-6730-437f-aff9-d2c38332a376"
 NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -56,7 +57,8 @@ Effects = "0.1,1"
 FreqTables = "0.4"
 GLM = "1"
 MixedModels = "4,5"
-MixedModelsMakie = "0.3.13"
+MixedModelsExtras = "1.1"
+MixedModelsMakie = "0.3.20"
 NLopt = "0.6"
 ProgressMeter = "1.7"
 RCall = "0.13"
diff --git a/_quarto.yml b/_quarto.yml
index 0a4bf16..acc1983 100644
--- a/_quarto.yml
+++ b/_quarto.yml
@@ -20,6 +20,7 @@ book:
     - datatables.qmd
     - linalg.qmd
     - aGHQ.qmd
+    - viz.qmd
 
 execute-dir: project
 
diff --git a/viz.jl b/viz.jl
new file mode 100644
index 0000000..3a64044
--- /dev/null
+++ b/viz.jl
@@ -0,0 +1,51 @@
+const hide_progress = true
+
+using AlgebraOfGraphics
+using CairoMakie
+using DataFrames
+if contains(first(Sys.cpu_info()).model, "Intel")
+  using MKL
+end
+using MixedModels
+using MixedModelsMakie
+using ProgressMeter
+using Random
+
+kb07 = DataFrame(MixedModels.dataset("kb07"))
+contrasts = Dict(:item => Grouping(),
+                 :subj => Grouping(),
+                 :spkr => EffectsCoding(),
+                 :prec => EffectsCoding(),
+                 :load => EffectsCoding())
+
+fm1 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1|item) + (1|subj)), kb07; contrasts, progress=true, thin=1)
+plot_blups(args...; kwargs...) = plot_blups!(Figure(; resolution=(1000, 600)), args...;kwargs...)
+function plot_blups!(f::Makie.FigureLike, m::MixedModel)
+
+    re = ranefinfo(m)
+    for (idx, grp) in enumerate(propertynames(re))
+      gl = f[1,idx] = GridLayout()
+      qqcaterpillar!(gl, re[grp])
+      Label(gl[end+1, :], string(grp); font=:bold, tellwidth=false)
+    end
+    Label(f[0, :], "Conditional Modes", tellwidth=false)
+    for idx in 2:length(re)
+      ratio = length(re[idx].cnames) / length(re[1].cnames)
+      colsize!(f.layout, idx, Auto(ratio))
+    end
+    return f
+end
+
+function eff_dim(m::MixedModel, threshold::Real=0.95)
+  re = MixedModels.rePCA(m)
+  return NamedTuple{propertynames(re)}(findfirst(>=(threshold), val) for val in values(re))
+end
+
+fm_max = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + spkr * prec * load |item) + (1 + spkr * prec * load|subj)), kb07; contrasts, progress=true, thin=1)
+fm2 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + prec + load|subj)), kb07; contrasts, progress=true, thin=1)
+fm3 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + prec + load|subj)), kb07; contrasts, progress=true, thin=1)
+fm4 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + load|subj)), kb07; contrasts, progress=true, thin=1)
+
+function plot_objective(m::MixedModel)
+    fitlog = optsum.fitlog
+end
diff --git a/viz.qmd b/viz.qmd
new file mode 100644
index 0000000..d29be3b
--- /dev/null
+++ b/viz.qmd
@@ -0,0 +1,73 @@
+---
+jupyter: julia-1.9
+---
+
+# Vizualization of Mixed Models {#sec-viz}
+
+
+```{julia}
+#| code-fold: true
+const hide_progress = true
+
+using AlgebraOfGraphics
+using CairoMakie
+using DataFrames
+if contains(first(Sys.cpu_info()).model, "Intel")
+  using MKL
+end
+using MixedModels
+using MixedModelsMakie
+using ProgressMeter
+using Random
+
+ProgressMeter.ijulia_behavior(:clear);
+CairoMakie.activate!(; type="svg");
+```
+
+```{julia}
+kb07 = DataFrame(MixedModels.dataset("kb07"))
+contrasts = Dict(:item => Grouping(),
+                 :subj => Grouping(),
+                 :spkr => EffectsCoding(),
+                 :prec => EffectsCoding(),
+                 :load => EffectsCoding())
+
+fm1 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1|item) + (1|subj)), kb07; contrasts, progress=true, thin=1)
+```
+
+
+```{julia}
+plot_blups(args...; kwargs...) = plot_blups!(Figure(; resolution=(1000, 600)), args...;kwargs...)
+function plot_blups!(f::Makie.FigureLike, m::MixedModel)
+
+    re = ranefinfo(m)
+    for (idx, grp) in enumerate(propertynames(re))
+      gl = f[1,idx] = GridLayout()
+      qqcaterpillar!(gl, re[grp])
+      Label(gl[end+1, :], string(grp); font=:bold, tellwidth=false)
+    end
+    Label(f[0, :], "Conditional Modes", tellwidth=false)
+    for idx in 2:length(re)
+      ratio = length(re[idx].cnames) / length(re[1].cnames)
+      colsize!(f.layout, idx, Auto(ratio))
+    end
+    return f
+end
+```
+
+```{julia}
+function eff_dim(m::MixedModel, threshold::Real=0.95)
+  re = MixedModels.rePCA(m)
+  return NamedTuple{propertynames(re)}(findfirst(>=(threshold), val) for val in values(re))
+end
+```
+
+```{julia}
+fm_max = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + spkr * prec * load |item) + (1 + spkr * prec * load|subj)), kb07; contrasts, progress=true, thin=1)
+```
+
+fm2 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + prec + load|subj)), kb07; contrasts, progress=true, thin=1)
+
+fm3 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + prec + load|subj)), kb07; contrasts, progress=true, thin=1)
+
+fm4 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + load|subj)), kb07; contrasts, progress=true, thin=1)