From 515286816bd384dfb745345658a284d957fb68d5 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Fri, 26 Dec 2025 23:12:51 -0600
Subject: [PATCH 1/9] more quarto ignore

---
 .gitignore | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.gitignore b/.gitignore
index e9fdf3e..7ec82e1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -49,3 +49,5 @@ site_libs/
 *.ipynb
 .cache/
 
+
+**/*.quarto_ipynb

From 09af25fec949dae8c6b43287b0e83c97d723b10b Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 29 Dec 2025 16:03:33 -0600
Subject: [PATCH 2/9] basic caching code

---
 src/SMLP2026.jl | 35 ++++++++++++++++++++++++++++++++++-
 1 file changed, 34 insertions(+), 1 deletion(-)

diff --git a/src/SMLP2026.jl b/src/SMLP2026.jl
index 08ce54b..de21321 100644
--- a/src/SMLP2026.jl
+++ b/src/SMLP2026.jl
@@ -6,6 +6,7 @@ using DataFrames
 using Dates
 using Downloads
 using Markdown
+using MixedModels
 using MixedModelsDatasets
 using PooledArrays
 using Scratch
@@ -39,6 +40,38 @@ end
 
 export GENRES,
     age_at_event,
-    tagpad
+    tagpad,
+    fit_or_restore,
+    fit_or_restore!
 
+function _normalize_path(path)
+    dir, file = splitdir(path)
+    if isempty(dir) 
+        return joinpath(dirname(@__DIR__), "fits", file)
+    else 
+        return path
+    end
+end
+
+function fit_or_restore(fname, args...;
+                         force_fit=false, fit_kwargs=(;), restore_kwargs=(;))
+    model = MixedModel(args...)
+    return fit_or_restore!(model, fname; force_fit, fit_kwargs, restore_kwargs)
+end
+    
+function fit_or_restore!(model::MixedModel, fname;
+                         force_fit=false, fit_kwargs=(;), restore_kwargs=(;))
+    fname = _normalize_path(fname)
+    if !isfile(fname) || force_fit
+        fit!(model; fitlog=true, fit_kwargs...)
+        saveoptsum(fname, model)
+    else
+        @info "Restoring from cache"
+        restoreoptsum!(model, fname; restore_kwargs...)
+    end
+
+    return model
+end
+
+    
 end # module EmbraceUncertainty

From 93745d34cde0109e71034e6b0cb4fe1025aef483 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 29 Dec 2025 17:01:13 -0600
Subject: [PATCH 3/9] model caching!

---
 fits/.gitattributes            |  1 +
 fits/kkl15_m_cpx.json.zip      |  3 ++
 fits/kkl15_m_max.json.zip      |  3 ++
 fits/kkl15_m_prm1.json.zip     |  3 ++
 fits/kkl15_m_prm2.json.zip     |  3 ++
 fits/kkl15_m_zcp1.json.zip     |  3 ++
 fits/kkl15_m_zcp1_rdc.json.zip |  3 ++
 fits/kkl15_m_zcp2.json.zip     |  3 ++
 kkl15.qmd                      | 93 +++++++++++++++++++++++++++++++---
 src/SMLP2026.jl                | 63 ++++++++++++++++-------
 10 files changed, 154 insertions(+), 24 deletions(-)
 create mode 100644 fits/.gitattributes
 create mode 100644 fits/kkl15_m_cpx.json.zip
 create mode 100644 fits/kkl15_m_max.json.zip
 create mode 100644 fits/kkl15_m_prm1.json.zip
 create mode 100644 fits/kkl15_m_prm2.json.zip
 create mode 100644 fits/kkl15_m_zcp1.json.zip
 create mode 100644 fits/kkl15_m_zcp1_rdc.json.zip
 create mode 100644 fits/kkl15_m_zcp2.json.zip

diff --git a/fits/.gitattributes b/fits/.gitattributes
new file mode 100644
index 0000000..486a232
--- /dev/null
+++ b/fits/.gitattributes
@@ -0,0 +1 @@
+*.zip filter=lfs diff=lfs merge=lfs -text
diff --git a/fits/kkl15_m_cpx.json.zip b/fits/kkl15_m_cpx.json.zip
new file mode 100644
index 0000000..80bc77a
--- /dev/null
+++ b/fits/kkl15_m_cpx.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bf87c12148db68ed590eae8bdd211f115b3d8ba7dad746c11e0046f124f5eecf
+size 58708
diff --git a/fits/kkl15_m_max.json.zip b/fits/kkl15_m_max.json.zip
new file mode 100644
index 0000000..54da652
--- /dev/null
+++ b/fits/kkl15_m_max.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4d529820ab07e70f5957da9e0efcd81fda911d02cf60c87718d232ff178679fb
+size 439320
diff --git a/fits/kkl15_m_prm1.json.zip b/fits/kkl15_m_prm1.json.zip
new file mode 100644
index 0000000..8e7bc15
--- /dev/null
+++ b/fits/kkl15_m_prm1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bc69caa8cb92eeede04f5d7b5430b8a59c867151f86561e635f9435efe942eac
+size 58708
diff --git a/fits/kkl15_m_prm2.json.zip b/fits/kkl15_m_prm2.json.zip
new file mode 100644
index 0000000..36ae6ae
--- /dev/null
+++ b/fits/kkl15_m_prm2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f1bf861cf57aa46e4fe8e681383ecb88afcf3f5e7cc47f2f1ce8beb6dd71741e
+size 34920
diff --git a/fits/kkl15_m_zcp1.json.zip b/fits/kkl15_m_zcp1.json.zip
new file mode 100644
index 0000000..e9030fa
--- /dev/null
+++ b/fits/kkl15_m_zcp1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8fe9590f918f440728206863aac10a1d0fe4e2070a9230c3d70bd983b96c4f3e
+size 16106
diff --git a/fits/kkl15_m_zcp1_rdc.json.zip b/fits/kkl15_m_zcp1_rdc.json.zip
new file mode 100644
index 0000000..0ab0a67
--- /dev/null
+++ b/fits/kkl15_m_zcp1_rdc.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3cdd880fde7ebf3c3799b97a39ed6be79406fafc4ef00d77ccd0c5729b29a4db
+size 7191
diff --git a/fits/kkl15_m_zcp2.json.zip b/fits/kkl15_m_zcp2.json.zip
new file mode 100644
index 0000000..43557ce
--- /dev/null
+++ b/fits/kkl15_m_zcp2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7d33b6162680cf7088d891ce21812826c9ed618dec05036fbde72098a5d44b4f
+size 7191
diff --git a/kkl15.qmd b/kkl15.qmd
index 534f0b4..c0c7969 100644
--- a/kkl15.qmd
+++ b/kkl15.qmd
@@ -48,6 +48,11 @@ using SMLP2026: dataset
 progress = isinteractive()
 ```
 
+```{julia}
+#| echo: false
+using SMLP2026: fit_or_restore
+```
+
 # Read data, compute and plot means
 
 ```{julia}
@@ -86,7 +91,7 @@ boxplot(
       ]
     )
   ),
-  figure=(; resolution=(800, 300)),
+  figure=(; size=(800, 300)),
 )
 ```
 
@@ -108,6 +113,7 @@ This is the maximum LMM for the design; `size` is a between-subject factor,
 ignoring other information such as trial number, age and gender of subjects.
 
 ```{julia}
+#| eval: false
 m_max = let
   form = @formula rt ~ 1 + CTR * cardinal * size +
                            (1 + CTR * cardinal | Subj)
@@ -115,6 +121,15 @@ m_max = let
 end;
 ```
 
+```{julia}
+#| include: false
+m_max = let
+  form = @formula rt ~ 1 + CTR * cardinal * size +
+                           (1 + CTR * cardinal | Subj)
+  fit_or_restore("kkl15_m_max.json", MixedModel, form, dat; contrasts)
+end
+```
+
 ```{julia}
 issingular(m_max)
 ```
@@ -136,6 +151,7 @@ The LMM `m_max` is overparameterized but it is not immediately apparent why.
 Force CPs to zero. _Reduction strategy 1_ is more suited for reducing model w/o theoretical expectations about CPs. The better reduction strategy for the present experiment with an _a priori_ interest in CPs is described as _Reduction strategy 2_.
 
 ```{julia}
+#| eval: false
 m_zcp1 = let
   form = @formula rt ~ 1 + CTR * cardinal * size +
                    zerocorr(1 + CTR * cardinal | Subj)
@@ -143,6 +159,16 @@ m_zcp1 = let
 end;
 ```
 
+```{julia}
+#| include: false
+m_zcp1 = let
+  form = @formula rt ~ 1 + CTR * cardinal * size +
+                  zerocorr(1 + CTR * cardinal | Subj)
+  fit_or_restore("kkl15_m_zcp1.json", MixedModel, form, dat; contrasts)
+end
+```
+
+
 ```{julia}
 issingular(m_zcp1)
 ```
@@ -162,6 +188,7 @@ The LMM `m_zcp1` is also overparameterized, but now there is clear evidence for
 Take out VCs for interactions.
 
 ```{julia}
+#| eval: false
 m_zcp1_rdc = let
   form = @formula rt ~ 1 + CTR * cardinal * size +
                    zerocorr(1 + CTR + cardinal | Subj)
@@ -169,6 +196,16 @@ m_zcp1_rdc = let
 end;
 ```
 
+```{julia}
+#| include: false
+m_zcp1_rdc = let
+  form = @formula rt ~ 1 + CTR * cardinal * size +
+                   zerocorr(1 + CTR + cardinal | Subj)
+  fit_or_restore("kkl15_m_zcp1_rdc.json", MixedModel, form, dat; contrasts)
+end
+```
+
+
 ```{julia}
 issingular(m_zcp1_rdc)
 ```
@@ -188,6 +225,7 @@ LMM `m_zcp_rdc` is ok . We add in CPs.
 Extend zcp-reduced LMM with CPs
 
 ```{julia}
+#| eval: false
 m_prm1 = let
   form = @formula rt ~ 1 + CTR * cardinal * size +
                            (1 + CTR + cardinal | Subj)
@@ -195,6 +233,16 @@ m_prm1 = let
 end;
 ```
 
+```{julia}
+#| include: false
+m_prm1 = let
+  form = @formula rt ~ 1 + CTR * cardinal * size +
+                           (1 + CTR + cardinal | Subj)
+  fit_or_restore("kkl15_m_prm1.json", MixedModel, form, dat; contrasts)
+end
+```
+
+
 ```{julia}
 issingular(m_prm1)
 ```
@@ -242,6 +290,7 @@ As the CPs were one reason for conducting this experiment, AIC is the criterion
 Relative to LMM `m_max`, first we take out interaction VCs and associated CPs, because these VCs are very small. This is the same as LMM `m_prm1` above.
 
 ```{julia}
+#| eval: false
 m_cpx = let
   form = @formula rt ~ 1 + CTR * cardinal * size +
                       (1 + CTR + cardinal | Subj)
@@ -249,6 +298,15 @@ m_cpx = let
 end;
 ```
 
+```{julia}
+#| include: false
+m_cpx = let
+  form = @formula rt ~ 1 + CTR * cardinal * size +
+                      (1 + CTR + cardinal | Subj)
+  fit_or_restore("kkl15_m_cpx.json", MixedModel, form, dat; contrasts)
+end
+```
+
 ## Zero-correlation parameter LMM (2)
 Now we check the significance of ensemble of CPs.
 
@@ -260,6 +318,14 @@ m_zcp2 = let
 end;
 ```
 
+```{julia}
+#| include: false
+m_zcp2 = let
+  form = @formula rt ~ 1 + CTR * cardinal * size  +
+              zerocorr(1 + CTR + cardinal | Subj)
+  fit_or_restore("kkl15_m_zcp2.json", MixedModel, form, dat; contrasts)
+end
+```
 
 ```{julia}
 VarCorr(m_zcp2)
@@ -270,6 +336,7 @@ VarCorr(m_zcp2)
 The cardinal-related CPs are quite small. Do we need them?
 
 ```{julia}
+#| eval: false
 m_prm2 = let
   form = @formula(rt ~ 1 + CTR * cardinal * size  +
                       (1 + CTR | Subj) + (0 + cardinal | Subj))
@@ -277,6 +344,15 @@ m_prm2 = let
 end;
 ```
 
+```{julia}
+#| include: false
+m_prm2 = let
+  form = @formula(rt ~ 1 + CTR * cardinal * size  +
+                      (1 + CTR | Subj) + (0 + cardinal | Subj))
+  fit_or_restore("kkl15_m_prm2.json", MixedModel, form, dat; contrasts)
+end
+```
+
 ```{julia}
 VarCorr(m_prm2)
 ```
@@ -371,7 +447,7 @@ qqnorm(
 #| label: fig-caterpillarm1
 #| fig-cap: Prediction intervals of the subject random effects in model m1
 cm1 = only(ranefinfo(m_prm1))
-caterpillar!(Figure(; resolution=(800, 1200)), cm1; orderby=2)
+caterpillar!(Figure(; size=(800, 1200)), cm1; orderby=2)
 ```
 
 ## Shrinkage plot
@@ -381,7 +457,7 @@ caterpillar!(Figure(; resolution=(800, 1200)), cm1; orderby=2)
 #| code-fold: true
 #| label: fig-caterpillarm1L
 #| fig-cap: Shrinkage plots of the subject random effects in model m1L
-shrinkageplot!(Figure(; resolution=(1000, 1200)), m_prm1)
+shrinkageplot!(Figure(; size=(1000, 1200)), m_prm1)
 ```
 
 # Parametric bootstrap
@@ -396,6 +472,9 @@ Here we
 
 We generate 2500 samples for the 15 model parameters (4 fixed effect, 7 VCs, 15 CPs, and 1 residual).
 
+<!--
+FIXME: cache bootstrap!
+-->
 ```{julia}
 samp = parametricbootstrap(MersenneTwister(1234321), 2500, m_prm1;
                            optsum_overrides=(; ftol_rel=1e-8));
@@ -432,7 +511,7 @@ draw(
   data(tbl) *
   mapping(:σ => "Residual") *
   density();
-  figure=(; resolution=(800, 400)),
+  figure=(; size=(800, 400)),
 )
 ```
 
@@ -465,7 +544,7 @@ draw(
     "Experimental effects",
   ) *
   density();
-  figure=(; resolution=(800, 350)),
+  figure=(; size=(800, 350)),
 )
 ```
 
@@ -484,7 +563,7 @@ draw(
     "Variance components",
   ) *
   density();
-  figure=(; resolution=(800, 350)),
+  figure=(; size=(800, 350)),
 )
 ```
 
@@ -506,7 +585,7 @@ draw(
     "Correlation parameters",
   ) *
   density();
-  figure=(; resolution=(800, 350)),
+  figure=(; size=(800, 350)),
 )
 ```
 
diff --git a/src/SMLP2026.jl b/src/SMLP2026.jl
index de21321..e71fb70 100644
--- a/src/SMLP2026.jl
+++ b/src/SMLP2026.jl
@@ -44,34 +44,63 @@ export GENRES,
     fit_or_restore,
     fit_or_restore!
 
-function _normalize_path(path)
+function _normalize_model_cache_path(path)
     dir, file = splitdir(path)
-    if isempty(dir) 
-        return joinpath(dirname(@__DIR__), "fits", file)
-    else 
-        return path
+    if isempty(dir)
+        path = joinpath(dirname(@__DIR__), "fits", file)
     end
+
+    _, ext = splitext(path)
+    if ext != ".zip"
+        path = path * ".zip"
+    end
+
+    return path
 end
 
-function fit_or_restore(fname, args...;
-                         force_fit=false, fit_kwargs=(;), restore_kwargs=(;))
-    model = MixedModel(args...)
-    return fit_or_restore!(model, fname; force_fit, fit_kwargs, restore_kwargs)
+function fit_or_restore(fname, ::Type{<:MixedModel}, args...; kwargs...)
+    return fit_or_restore(fname, args...; kwargs...)
 end
-    
+
+function fit_or_restore(fname, args...; contrasts=Dict{Symbol}(), kwargs...)
+    model = MixedModel(args...; contrasts)
+    return fit_or_restore!(model, fname; kwargs...)
+end
+
 function fit_or_restore!(model::MixedModel, fname;
-                         force_fit=false, fit_kwargs=(;), restore_kwargs=(;))
-    fname = _normalize_path(fname)
+                         force_fit=false, restore_kwargs=(; atol=1e-8), fit_kwargs...)
+    fname = _normalize_model_cache_path(fname)
+    @debug "cache path: $(fname)"
     if !isfile(fname) || force_fit
-        fit!(model; fitlog=true, fit_kwargs...)
-        saveoptsum(fname, model)
+        @debug "fitting model"
+        fit!(model; fit_kwargs...)
+        zip = ZipFile.Writer(fname)
+        try
+            mktempdir() do dir
+                f = ZipFile.addfile(zip, "model.json"; method=ZipFile.Deflate)
+                saveoptsum(f, model)
+                return nothing
+            end
+        catch
+            error("Something went wrong in saving the model cache")
+        finally
+            close(zip)
+        end
     else
-        @info "Restoring from cache"
-        restoreoptsum!(model, fname; restore_kwargs...)
+        @debug "restoring from cache"
+        zip = ZipFile.Reader(fname)
+        try
+            restoreoptsum!(model, only(zip.files); restore_kwargs...)
+        catch
+            error("Something went wrong in reading the model cache")
+        finally
+            close(zip)
+        end
     end
 
     return model
 end
 
-    
+# TODO: wrapper function for bootstrap, which crosschecks the PRNG
+
 end # module EmbraceUncertainty

From 00ee9fee81916cd7758ac7ef068d4ac729666110 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 29 Dec 2025 23:44:06 -0600
Subject: [PATCH 4/9] bootstrap caching

---
 fits/.gitattributes     |  1 +
 fits/kkl15_m_prm1.arrow |  3 +++
 kkl15.qmd               | 12 ++++++++----
 src/SMLP2026.jl         | 32 ++++++++++++++++++++++++++++----
 4 files changed, 40 insertions(+), 8 deletions(-)
 create mode 100644 fits/kkl15_m_prm1.arrow

diff --git a/fits/.gitattributes b/fits/.gitattributes
index 486a232..24df8a2 100644
--- a/fits/.gitattributes
+++ b/fits/.gitattributes
@@ -1 +1,2 @@
 *.zip filter=lfs diff=lfs merge=lfs -text
+*.arrow filter=lfs diff=lfs merge=lfs -text
diff --git a/fits/kkl15_m_prm1.arrow b/fits/kkl15_m_prm1.arrow
new file mode 100644
index 0000000..68df100
--- /dev/null
+++ b/fits/kkl15_m_prm1.arrow
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d0cfdda1d6739ad26f6070df190027a5c3df0b14589d4c44718eb9b9e72f8f1a
+size 1004138
diff --git a/kkl15.qmd b/kkl15.qmd
index c0c7969..b269498 100644
--- a/kkl15.qmd
+++ b/kkl15.qmd
@@ -50,7 +50,7 @@ progress = isinteractive()
 
 ```{julia}
 #| echo: false
-using SMLP2026: fit_or_restore
+using SMLP2026: fit_or_restore, bootstrap_or_restore
 ```
 
 # Read data, compute and plot means
@@ -472,14 +472,18 @@ Here we
 
 We generate 2500 samples for the 15 model parameters (4 fixed effect, 7 VCs, 15 CPs, and 1 residual).
 
-<!--
-FIXME: cache bootstrap!
--->
 ```{julia}
+#| eval: false
 samp = parametricbootstrap(MersenneTwister(1234321), 2500, m_prm1;
                            optsum_overrides=(; ftol_rel=1e-8));
 ```
 
+```{julia}
+#| include: false
+samp = bootstrap_or_restore("kkl15_m_prm1.arrow", MersenneTwister(1234321), 2500, m_prm1;
+                            optsum_overrides=(; ftol_rel=1e-8))
+```
+
 ```{julia}
 tbl = samp.tbl
 ```
diff --git a/src/SMLP2026.jl b/src/SMLP2026.jl
index e71fb70..7f28e3c 100644
--- a/src/SMLP2026.jl
+++ b/src/SMLP2026.jl
@@ -9,6 +9,7 @@ using Markdown
 using MixedModels
 using MixedModelsDatasets
 using PooledArrays
+using Random
 using Scratch
 using SHA
 using TypedTables
@@ -44,12 +45,17 @@ export GENRES,
     fit_or_restore,
     fit_or_restore!
 
-function _normalize_model_cache_path(path)
+function _normalize_cache_path(path)
     dir, file = splitdir(path)
     if isempty(dir)
         path = joinpath(dirname(@__DIR__), "fits", file)
     end
 
+    return path
+end
+
+function _normalize_model_cache_path(path)
+    path = _normalize_cache_path(path)
     _, ext = splitext(path)
     if ext != ".zip"
         path = path * ".zip"
@@ -68,10 +74,10 @@ function fit_or_restore(fname, args...; contrasts=Dict{Symbol}(), kwargs...)
 end
 
 function fit_or_restore!(model::MixedModel, fname;
-                         force_fit=false, restore_kwargs=(; atol=1e-8), fit_kwargs...)
+                         force=false, restore_kwargs=(; atol=1e-8), fit_kwargs...)
     fname = _normalize_model_cache_path(fname)
     @debug "cache path: $(fname)"
-    if !isfile(fname) || force_fit
+    if !isfile(fname) || force
         @debug "fitting model"
         fit!(model; fit_kwargs...)
         zip = ZipFile.Writer(fname)
@@ -101,6 +107,24 @@ function fit_or_restore!(model::MixedModel, fname;
     return model
 end
 
-# TODO: wrapper function for bootstrap, which crosschecks the PRNG
+# TODO: cache invalidation if PRNG / replicates don't match
+function bootstrap_or_restore(fname,  args...; kwargs...)
+    return bootstrap_or_restore(fname, Random.default_rng(), args...; kwargs...)
+end
+function bootstrap_or_restore(fname, rng::AbstractRNG, n::Integer, model::MixedModel, args...;
+                              force=false, bootstrap_kwargs...)
+    fname = _normalize_cache_path(fname)
+    @debug "cache path: $(fname)"
+    if !isfile(fname) || force
+        @debug "performing bootstrap"
+        boot = parametricbootstrap(rng, n, model, args...; bootstrap_kwargs...)
+        savereplicates(fname, boot)
+    else
+        @debug "restoring from cache"
+        boot = restorereplicates(fname, model)
+    end
+
+    return boot
+end
 
 end # module EmbraceUncertainty

From c270506059d3ac235688f2af8e8cba8c2bdd9ec5 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 29 Dec 2025 23:44:31 -0600
Subject: [PATCH 5/9] s/resolution/size/g

---
 bootstrap.qmd          |  2 +-
 fggk21.qmd             |  6 +++---
 kb07.qmd               |  6 +++---
 kwdyz11.qmd            |  6 +++---
 largescaledesigned.qmd | 10 +++++-----
 profiling.qmd          | 12 ++++++------
 selection.qmd          |  2 +-
 sleepstudy.qmd         | 22 +++++++++++-----------
 sleepstudy_speed.qmd   |  4 ++--
 9 files changed, 35 insertions(+), 35 deletions(-)

diff --git a/bootstrap.qmd b/bootstrap.qmd
index 438ae9f..b56a90d 100644
--- a/bootstrap.qmd
+++ b/bootstrap.qmd
@@ -150,7 +150,7 @@ confint(samp; method=:shortest)
 ```{julia}
 draw(
   data(samp.β) * mapping(:β; color=:coefname) * AoG.density();
-  figure=(; resolution=(800, 450)),
+  figure=(; size=(800, 450)),
 )
 ```
 
diff --git a/fggk21.qmd b/fggk21.qmd
index b984153..cde3a8a 100644
--- a/fggk21.qmd
+++ b/fggk21.qmd
@@ -926,7 +926,7 @@ Let's take pictures.
 #| label: fig-caterpillarCohort
 #| fig-cap: Prediction intervals of the random effects for Cohort in model m1
 caterpillar!(
-  Figure(; resolution=(800, 400)), ranefinfo(m1, :Cohort)
+  Figure(; size=(800, 400)), ranefinfo(m1, :Cohort)
 )
 ```
 
@@ -938,14 +938,14 @@ These are just teasers. We will pick this up in a separate tutorial. Enjoy!
 #| code-fold: true
 #| label: fig-m1shrinkageCohort
 #| fig-cap: Shrinkage plot of the random effects for Cohort in model m1
-shrinkageplot!(Figure(; resolution=(800, 800)), m1, :Cohort)
+shrinkageplot!(Figure(; size=(800, 800)), m1, :Cohort)
 ```
 
 ```{julia}
 #| code-fold: true
 #| label: fig-m2shrinkageCohort
 #| fig-cap: Shrinkage plot of the random effects for Cohort in model m2
-shrinkageplot!(Figure(; resolution=(800, 800)), m2, :Cohort)
+shrinkageplot!(Figure(; size=(800, 800)), m2, :Cohort)
 ```
 
 ::: {#refs}
diff --git a/kb07.qmd b/kb07.qmd
index e4c69d7..74abf72 100644
--- a/kb07.qmd
+++ b/kb07.qmd
@@ -116,7 +116,7 @@ For the fixed-effects parameters the estimators are close to being normally dist
 #| label: fig-kbm02fedens
 draw(
   data(kbm02samp.β) * mapping(:β; color=:coefname) * AlgebraOfGraphics.density();
-  figure=(; resolution=(800, 450)),
+  figure=(; size=(800, 450)),
 )
 ```
 
@@ -129,7 +129,7 @@ let pars = ["σ1", "σ2", "σ3"]
     data(kbm02tbl) *
     mapping(pars .=> "σ"; color=dims(1) => renamer(pars)) *
     AlgebraOfGraphics.density();
-    figure=(; resolution=(800, 450)),
+    figure=(; size=(800, 450)),
   )
 end
 ```
@@ -142,7 +142,7 @@ draw(
   data(kbm02tbl) *
   mapping(:ρ1 => "Correlation") *
   AlgebraOfGraphics.density();
-  figure=(; resolution=(800, 450)),
+  figure=(; size=(800, 450)),
 )
 ```
 
diff --git a/kwdyz11.qmd b/kwdyz11.qmd
index 89c43d3..197221c 100644
--- a/kwdyz11.qmd
+++ b/kwdyz11.qmd
@@ -116,7 +116,7 @@ boxplot(
       ],
     ),
   ),
-  figure=(; resolution=(800, 300)),
+  figure=(; size=(800, 300)),
 )
 ```
 
@@ -247,7 +247,7 @@ The caterpillar plot, @fig-m1caterpillar, also reveals the high correlation betw
 #| fig-cap: "Prediction intervals on the random effects for Subj in model m1"
 #| label: fig-m1caterpillar
 caterpillar!(
-  Figure(; resolution=(800, 1000)), ranefinfo(m1, :Subj); orderby=2
+  Figure(; size=(800, 1000)), ranefinfo(m1, :Subj); orderby=2
 )
 ```
 
@@ -260,7 +260,7 @@ The corresponding panel illustrates an *implosion* of conditional modes.
 #| code-fold: true
 #| fig-cap: "Shrinkage plot of the conditional means of the random effects for model m1"
 #| label: fig-m1shrinkage
-shrinkageplot!(Figure(; resolution=(1000, 1000)), m1)
+shrinkageplot!(Figure(; size=(1000, 1000)), m1)
 ```
 
 # Parametric bootstrap
diff --git a/largescaledesigned.qmd b/largescaledesigned.qmd
index 2a709e9..b9cf67d 100644
--- a/largescaledesigned.qmd
+++ b/largescaledesigned.qmd
@@ -193,7 +193,7 @@ let
       color=:isword,
       markersize=:msz,
     );
-    figure=(; resolution=(800, 450)),
+    figure=(; size=(800, 450)),
   )
 end
 ```
@@ -297,7 +297,7 @@ draw(
   data(pruned) *
   mapping(:rt => "Response time (ms.) for correct responses") *
   AlgebraOfGraphics.density();
-  figure=(; resolution=(800, 450)),
+  figure=(; size=(800, 450)),
 )
 ```
 
@@ -314,7 +314,7 @@ draw(
   mapping(
     :rt => (x -> 1000 / x) => "Response speed (s⁻¹) for correct responses") *
   AlgebraOfGraphics.density();
-  figure=(; resolution=(800, 450)),
+  figure=(; size=(800, 450)),
 )
 ```
 
@@ -562,7 +562,7 @@ A caterpillar plot, @fig-elm01caterpillarsubj,
 #| fig-cap: Conditional means and 95% prediction intervals for subject random effects in elm01.
 #| label: fig-elm01caterpillarsubj
 qqcaterpillar!(
-  Figure(resolution=(800, 650)),
+  Figure(size=(800, 650)),
   ranefinfo(elm01, :subj),
 )
 ```
@@ -577,7 +577,7 @@ Also, there is at least one outlier with a conditional mode over 1.0.
 #| fig-cap: Conditional means and 95% prediction intervals for subject random effects in elm02.
 #| label: fig-elm02caterpillarsubj
 qqcaterpillar!(
-  Figure(resolution=(800, 650)),
+  Figure(size=(800, 650)),
   ranefinfo(elm02, :subj),
 )
 ```
diff --git a/profiling.qmd b/profiling.qmd
index 4c4c84b..d0bd0ca 100644
--- a/profiling.qmd
+++ b/profiling.qmd
@@ -117,7 +117,7 @@ To review:
 #| code-fold: true
 #| fig-cap: "ζ versus the value of the coefficient for the fixed-effects parameters in a model of response speed for the kb07 data."
 #| label: fig-kb07zetabeta
-zetaplot!(Figure(; resolution=(1200, 350)), pr01; ptyp='β')
+zetaplot!(Figure(; size=(1200, 350)), pr01; ptyp='β')
 ```
 
 The lines on these panels are read like normal probability plots, i.e. QQ plots against a standard normal distribution.
@@ -137,7 +137,7 @@ as shown in @fig-kb07abszetabeta, which shows the absolute value of $\zeta$, whi
 #| code-fold: true
 #| fig-cap: "Absolute value of ζ versus value of the coefficient for the fixed-effects parameters in a model of response speed for the kb07 data.  The horizontal lines are confidence intervals with nominal 50%, 80%, 90%, 95% and 99% confidence."
 #| label: fig-kb07abszetabeta
-zetaplot!(Figure(; resolution=(1200, 330)), pr01; ptyp='β', absv=true)
+zetaplot!(Figure(; size=(1200, 330)), pr01; ptyp='β', absv=true)
 ```
 
 The 95% confidence intervals are the second horizontal lines from the top in each panel, at 1.96 on the vertical scale.
@@ -146,7 +146,7 @@ The 95% confidence intervals are the second horizontal lines from the top in eac
 #| code-fold: true
 #| fig-cap: "Absolute value of ζ versus value of the coefficient for the variance component parameters in a model of response speed for the kb07 data.  The horizontal lines are confidence intervals with nominal 50%, 80%, 90%, 95% and 99% confidence."
 #| label: fig-kb07abszetasigma
-zetaplot!(Figure(; resolution=(1200, 330)), pr01; ptyp='σ', absv=true)
+zetaplot!(Figure(; size=(1200, 330)), pr01; ptyp='σ', absv=true)
 ```
 
 @fig-kb07abszetasigma shows similar confidence intervals on the parameters representing standard deviations as does @fig-kb07abszetatheta for the $\theta$ parameters.
@@ -155,7 +155,7 @@ zetaplot!(Figure(; resolution=(1200, 330)), pr01; ptyp='σ', absv=true)
 #| code-fold: true
 #| fig-cap: "Absolute value of ζ versus parameter value for the θ parameters in a model of response speed for the kb07 data.  The horizontal lines are confidence intervals with nominal 50%, 80%, 90%, 95% and 99% confidence."
 #| label: fig-kb07abszetatheta
-zetaplot!(Figure(; resolution=(1200, 330)), pr01; ptyp='θ', absv=true)
+zetaplot!(Figure(; size=(1200, 330)), pr01; ptyp='θ', absv=true)
 ```
 
 # Comparisons with the parametric bootstrap
@@ -192,7 +192,7 @@ The profile-$\zeta$ function can be transformed to an equivalent density functio
 #| code-fold: true
 #| fig-cap: "Density functions of the marginal distribution of the estimators of variance component parameters in a model of response speed in the kb07 data."
 #| label: fig-kb07densigma
-profiledensity!(Figure(; resolution=(1200, 300)), pr01; share_y_scale=false)
+profiledensity!(Figure(; size=(1200, 300)), pr01; share_y_scale=false)
 ```
 
 Alternatively we can see the skewness in the plots of $\zeta$
@@ -201,7 +201,7 @@ Alternatively we can see the skewness in the plots of $\zeta$
 #| code-fold: true
 #| label: fig-kb07zetasigma
 #| fig-cap: "Plot of ζ scores for the variance component parameters in a model of response speed for the kb07 data"
-zetaplot!(Figure(; resolution=(1100, 330)), pr01; ptyp='σ')
+zetaplot!(Figure(; size=(1100, 330)), pr01; ptyp='σ')
 ```
 
 The skewness is less obvious in @fig-kb07bskdesigma because of the stochastic nature of the bootstrap sample and the kernel density estimators.
diff --git a/selection.qmd b/selection.qmd
index f5869b0..df3e6e0 100644
--- a/selection.qmd
+++ b/selection.qmd
@@ -25,7 +25,7 @@ describe(df)
 ```{julia}
 #| eval: false
 let
-  fdensity = Figure(; resolution=(1000, 500))
+  fdensity = Figure(; size=(1000, 500))
   axs = Axis(fdensity[1, 1])
   tdf = filter(:Test => ==(test), df)
   colors = Makie.cgrad(:PuOr_4, 2; categorical=true, alpha=0.6)
diff --git a/sleepstudy.qmd b/sleepstudy.qmd
index 72224be..243de32 100644
--- a/sleepstudy.qmd
+++ b/sleepstudy.qmd
@@ -53,8 +53,8 @@ sleepstudy = DataFrame(dataset("sleepstudy"))
 let f = Figure(; size=(700, 400))
   yrange = maximum(sleepstudy.reaction) - minimum(sleepstudy.reaction)
   xrange = maximum(sleepstudy.days) - minimum(sleepstudy.days)
-  
-  reg = combine(groupby(sleepstudy, :subj), 
+
+  reg = combine(groupby(sleepstudy, :subj),
                 [:days, :reaction] => NamedTuple{(:intercept, :slope)} ∘ simplelinreg => AsTable)
   sort!(reg, :intercept)
 
@@ -70,7 +70,7 @@ let f = Figure(; size=(700, 400))
   # set up all the axes and plot the simple regression lines
   for row in eachrow(reg)
     pos = gridpos[row.subj]
-    ax = Axis(f[pos...]; title=row.subj, 
+    ax = Axis(f[pos...]; title=row.subj,
               autolimitaspect=xrange/yrange)
     if pos[1] == 1
       hidexdecorations!(ax; grid=false, ticks=false)
@@ -87,11 +87,11 @@ let f = Figure(; size=(700, 400))
     pos = gridpos[grouping.subj]
     scatter!(f[pos...], gdf.days, gdf.reaction)
   end
-  Label(f[end+1, :], "Days of sleep deprivation"; 
+  Label(f[end+1, :], "Days of sleep deprivation";
         tellwidth=false, tellheight=true)
-  Label(f[:, 0], "Average reaction time (ms)"; 
+  Label(f[:, 0], "Average reaction time (ms)";
         tellwidth=true, tellheight=false, rotation=pi/2)
-  
+
   linkaxes!(axes...)
 
   # tweak the layout a little
@@ -227,7 +227,7 @@ If the BLUPs are strongly shrunk towards zero then the additional complexity in
 #| code-fold: true
 #| fig-cap: Shrinkage plot of means of the random effects in model m1
 #| label: fig-m1shrinkage
-shrinkageplot!(Figure(; resolution=(500, 500)), m1)
+shrinkageplot!(Figure(; size=(500, 500)), m1)
 ```
 
 ::: {.callout-note}
@@ -269,7 +269,7 @@ An empirical density plot of the estimates for the fixed-effects coefficients, @
 #| fig-cap: 'Empirical density plots of bootstrap replications of fixed-effects parameter estimates'
 #| label: fig-bsbetadensity
 begin
-  f1 = Figure(; resolution=(1000, 400))
+  f1 = Figure(; size=(1000, 400))
   CairoMakie.density!(
     Axis(f1[1, 1]; xlabel="Intercept [ms]"), tbl.β1
   )
@@ -301,7 +301,7 @@ hist(
   tbl.ρ1;
   bins=40,
   axis=(; xlabel="Estimated correlation of the random effects"),
-  figure=(; resolution=(500, 500)),
+  figure=(; size=(500, 500)),
 )
 ```
 
@@ -312,7 +312,7 @@ Finally, density plots for the variance components (but on the scale of the stan
 #| fig-cap: 'Empirical density plots of bootstrap replicates of standard deviation estimates'
 #| label: fig-bssigmadensity
 begin
-  f2 = Figure(; resolution=(1000, 300))
+  f2 = Figure(; size=(1000, 300))
   CairoMakie.density!(
     Axis(f2[1, 1]; xlabel="Residual σ"),
     tbl.σ,
@@ -344,7 +344,7 @@ let
     tbl.β1, tbl.β2,
     color=(:blue, 0.20),
     axis=(; xlabel="Intercept", ylabel="Coefficient of days"),
-    figure=(; resolution=(500, 500)),
+    figure=(; size=(500, 500)),
   )
   contour!(kde((tbl.β1, tbl.β2)))
   current_figure()
diff --git a/sleepstudy_speed.qmd b/sleepstudy_speed.qmd
index 7d9b5a4..420322d 100644
--- a/sleepstudy_speed.qmd
+++ b/sleepstudy_speed.qmd
@@ -95,7 +95,7 @@ The panels are arranged by increasing initial reaction speed starting at the low
 let
   ord = sortperm(first.(within.coef))
   labs = values(only.(keys(gdf)))[ord]       # labels for panels
-  f = clevelandaxes!(Figure(; resolution=(1000, 750)), labs, (2, 9))
+  f = clevelandaxes!(Figure(; size=(1000, 750)), labs, (2, 9))
   for (axs, sdf) in zip(f.content, gdf[ord]) # iterate over the panels and groups
     scatter!(axs, sdf.day, sdf.speed)      # add the points
     coef = simplelinreg(sdf.day, sdf.speed)
@@ -176,7 +176,7 @@ caterpillar(m2)
 #| code-fold: true
 #| fig-cap: "Shrinkage plot of the means of the random effects in model m2"
 #| label: fig-shrinkagem2
-shrinkageplot!(Figure(; resolution=(500, 500)), m2)
+shrinkageplot!(Figure(; size=(500, 500)), m2)
 ```
 
 # References

From 2375b275c21d12b7aea2388cfe3adf673dbad627 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 30 Mar 2026 16:48:54 -0500
Subject: [PATCH 6/9] use claude to update the remaining files

---
 Project.toml                               |   2 +-
 bootstrap.qmd                              |  36 ++-
 contrasts_fggk21.qmd                       |   3 +-
 contrasts_kwdyz11.qmd                      | 252 ++++++++++++++++++++-
 fggk21.qmd                                 |  63 +++++-
 fits/bootstrap_m0.json.zip                 |   3 +
 fits/bootstrap_m1.json.zip                 |   3 +
 fits/bootstrap_mdye.json.zip               |   3 +
 fits/bootstrap_msleep.json.zip             |   3 +
 fits/contrasts_kwdyz11_m1.json.zip         |   3 +
 fits/contrasts_kwdyz11_m1b.json.zip        |   3 +
 fits/contrasts_kwdyz11_m2.json.zip         |   3 +
 fits/contrasts_kwdyz11_m2b.json.zip        |   3 +
 fits/contrasts_kwdyz11_m2c.json.zip        |   3 +
 fits/contrasts_kwdyz11_m2d.json.zip        |   3 +
 fits/contrasts_kwdyz11_m3.json.zip         |   3 +
 fits/contrasts_kwdyz11_m3b.json.zip        |   3 +
 fits/contrasts_kwdyz11_m3c.json.zip        |   3 +
 fits/contrasts_kwdyz11_m4.json.zip         |   3 +
 fits/contrasts_kwdyz11_m4b.json.zip        |   3 +
 fits/contrasts_kwdyz11_m4c.json.zip        |   3 +
 fits/contrasts_kwdyz11_m4d.json.zip        |   3 +
 fits/contrasts_kwdyz11_m5.json.zip         |   3 +
 fits/contrasts_kwdyz11_m6.json.zip         |   3 +
 fits/contrasts_kwdyz11_m6_cpx.json.zip     |   3 +
 fits/contrasts_kwdyz11_m6_ovi.json.zip     |   3 +
 fits/contrasts_kwdyz11_m6_zcp.json.zip     |   3 +
 fits/fggk21_m_cpx.json.zip                 |   3 +
 fits/fggk21_m_cpxCohort.json.zip           |   3 +
 fits/fggk21_m_ovi.json.zip                 |   3 +
 fits/fggk21_m_oviCohort.json.zip           |   3 +
 fits/fggk21_m_zcp.json.zip                 |   3 +
 fits/fggk21_m_zcpCohort.json.zip           |   3 +
 fits/glmm_gm1.json.zip                     |   3 +
 fits/glmm_gm2.json.zip                     |   3 +
 fits/glmm_gm3.json.zip                     |   3 +
 fits/glmm_m1.json.zip                      |   3 +
 fits/kwdyz11_m1.json.zip                   |   3 +
 fits/largescaledesigned_elm01.json.zip     |   3 +
 fits/largescaledesigned_elm02.json.zip     |   3 +
 fits/mrk17_m_cpx.json.zip                  |   3 +
 fits/mrk17_m_mrk17_crossed.json.zip        |   3 +
 fits/mrk17_m_mrk17_nested.json.zip         |   3 +
 fits/mrk17_m_prm.json.zip                  |   3 +
 fits/mrk17_m_prm3.json.zip                 |   3 +
 fits/mrk17_m_prm3_posthoc.json.zip         |   3 +
 fits/mrk17_m_prm_1.json.zip                |   3 +
 fits/mrk17_m_prm_2.json.zip                |   3 +
 fits/mrk17_m_prm_3.json.zip                |   3 +
 fits/mrk17_m_prm_4.json.zip                |   3 +
 fits/mrk17_m_prm_5.json.zip                |   3 +
 fits/mrk17_m_zcp.json.zip                  |   3 +
 fits/mrk17_m_zcp2.json.zip                 |   3 +
 fits/mrk17_m_zcp3.json.zip                 |   3 +
 fits/mrk17_m_zcp4.json.zip                 |   3 +
 fits/partial_within_m0.json.zip            |   3 +
 fits/partial_within_m1.json.zip            |   3 +
 fits/profiling_pr01.json.zip               |   3 +
 fits/shrinkageplot_m1.json.zip             |   3 +
 fits/shrinkageplot_m2.json.zip             |   3 +
 fits/shrinkageplot_m3.json.zip             |   3 +
 fits/shrinkageplot_m4.json.zip             |   3 +
 fits/sleepstudy_m1.json.zip                |   3 +
 fits/sleepstudy_m2.json.zip                |   3 +
 fits/sleepstudy_speed_m1.json.zip          |   3 +
 fits/sleepstudy_speed_m2.json.zip          |   3 +
 fits/transformation_days_centered.json.zip |   3 +
 fits/transformation_model_bc.json.zip      |   3 +
 fits/transformation_slp.json.zip           |   3 +
 glmm.qmd                                   |  64 +++++-
 kwdyz11.qmd                                |  14 +-
 largescaledesigned.qmd                     |  25 +-
 mrk17.qmd                                  | 161 ++++++++++++-
 partial_within.qmd                         |  24 +-
 power_simulation.qmd                       |  16 +-
 profiling.qmd                              |  12 +-
 shrinkageplot.qmd                          |  57 +++++
 sleepstudy.qmd                             |  19 +-
 sleepstudy_speed.qmd                       |  20 +-
 src/SMLP2026.jl                            |   2 +-
 transformation.qmd                         |  33 ++-
 81 files changed, 970 insertions(+), 25 deletions(-)
 create mode 100644 fits/bootstrap_m0.json.zip
 create mode 100644 fits/bootstrap_m1.json.zip
 create mode 100644 fits/bootstrap_mdye.json.zip
 create mode 100644 fits/bootstrap_msleep.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m1.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m1b.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m2.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m2b.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m2c.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m2d.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m3.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m3b.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m3c.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m4.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m4b.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m4c.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m4d.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m5.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m6.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m6_cpx.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m6_ovi.json.zip
 create mode 100644 fits/contrasts_kwdyz11_m6_zcp.json.zip
 create mode 100644 fits/fggk21_m_cpx.json.zip
 create mode 100644 fits/fggk21_m_cpxCohort.json.zip
 create mode 100644 fits/fggk21_m_ovi.json.zip
 create mode 100644 fits/fggk21_m_oviCohort.json.zip
 create mode 100644 fits/fggk21_m_zcp.json.zip
 create mode 100644 fits/fggk21_m_zcpCohort.json.zip
 create mode 100644 fits/glmm_gm1.json.zip
 create mode 100644 fits/glmm_gm2.json.zip
 create mode 100644 fits/glmm_gm3.json.zip
 create mode 100644 fits/glmm_m1.json.zip
 create mode 100644 fits/kwdyz11_m1.json.zip
 create mode 100644 fits/largescaledesigned_elm01.json.zip
 create mode 100644 fits/largescaledesigned_elm02.json.zip
 create mode 100644 fits/mrk17_m_cpx.json.zip
 create mode 100644 fits/mrk17_m_mrk17_crossed.json.zip
 create mode 100644 fits/mrk17_m_mrk17_nested.json.zip
 create mode 100644 fits/mrk17_m_prm.json.zip
 create mode 100644 fits/mrk17_m_prm3.json.zip
 create mode 100644 fits/mrk17_m_prm3_posthoc.json.zip
 create mode 100644 fits/mrk17_m_prm_1.json.zip
 create mode 100644 fits/mrk17_m_prm_2.json.zip
 create mode 100644 fits/mrk17_m_prm_3.json.zip
 create mode 100644 fits/mrk17_m_prm_4.json.zip
 create mode 100644 fits/mrk17_m_prm_5.json.zip
 create mode 100644 fits/mrk17_m_zcp.json.zip
 create mode 100644 fits/mrk17_m_zcp2.json.zip
 create mode 100644 fits/mrk17_m_zcp3.json.zip
 create mode 100644 fits/mrk17_m_zcp4.json.zip
 create mode 100644 fits/partial_within_m0.json.zip
 create mode 100644 fits/partial_within_m1.json.zip
 create mode 100644 fits/profiling_pr01.json.zip
 create mode 100644 fits/shrinkageplot_m1.json.zip
 create mode 100644 fits/shrinkageplot_m2.json.zip
 create mode 100644 fits/shrinkageplot_m3.json.zip
 create mode 100644 fits/shrinkageplot_m4.json.zip
 create mode 100644 fits/sleepstudy_m1.json.zip
 create mode 100644 fits/sleepstudy_m2.json.zip
 create mode 100644 fits/sleepstudy_speed_m1.json.zip
 create mode 100644 fits/sleepstudy_speed_m2.json.zip
 create mode 100644 fits/transformation_days_centered.json.zip
 create mode 100644 fits/transformation_model_bc.json.zip
 create mode 100644 fits/transformation_slp.json.zip

diff --git a/Project.toml b/Project.toml
index a585c45..8624414 100644
--- a/Project.toml
+++ b/Project.toml
@@ -70,5 +70,5 @@ Random = "1"
 RegressionFormulae = "0.1.3"
 StableRNGs = "1.0.3"
 StandardizedPredictors = "1"
-StatsModels = "0.7 - 0.7.4, 0.7.6"
+StatsModels = "0.7.6"
 julia = "1.10"
diff --git a/bootstrap.qmd b/bootstrap.qmd
index b56a90d..fca1686 100644
--- a/bootstrap.qmd
+++ b/bootstrap.qmd
@@ -38,7 +38,7 @@ using DataFrames
 using MixedModels
 using MixedModelsMakie
 using Random
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 
 using AlgebraOfGraphics: AlgebraOfGraphics as AoG
 const progress=false
@@ -77,12 +77,21 @@ The `EffectsCoding` contrast is used with these to create a ±1 encoding.
 We can look at an initial fit of moderate complexity:
 
 ```{julia}
+#| eval: false
 form = @formula(rt_trunc ~ 1 + spkr * prec * load +
                           (1 + spkr + prec + load | subj) +
                           (1 + spkr + prec + load | item))
 m0 = fit(MixedModel, form, kb07; contrasts, progress)
 ```
 
+```{julia}
+#| include: false
+form = @formula(rt_trunc ~ 1 + spkr * prec * load +
+                          (1 + spkr + prec + load | subj) +
+                          (1 + spkr + prec + load | item))
+m0 = fit_or_restore("bootstrap_m0.json", MixedModel, form, kb07; contrasts, progress)
+```
+
 The default display in Quarto uses the [pretty MIME show method](https://juliastats.org/MixedModels.jl/v4.7.1/mime/) for the model and omits the estimated correlations of the random effects.
 
 The `VarCorr` extractor displays these.
@@ -95,12 +104,21 @@ None of the two-factor or three-factor interaction terms in the fixed-effects ar
 In the random-effects terms only the scalar random effects and the `prec` random effect for `item` appear to be warranted, leading to the reduced formula
 
 ```{julia}
+#| eval: false
 # formula f4 from https://doi.org/10.33016/nextjournal.100002
 form = @formula(rt_trunc ~ 1 + spkr * prec * load + (1 | subj) + (1 + prec | item))
 
 m1 = fit(MixedModel, form, kb07; contrasts, progress)
 ```
 
+```{julia}
+#| include: false
+# formula f4 from https://doi.org/10.33016/nextjournal.100002
+form = @formula(rt_trunc ~ 1 + spkr * prec * load + (1 | subj) + (1 + prec | item))
+
+m1 = fit_or_restore("bootstrap_m1.json", MixedModel, form, kb07; contrasts, progress)
+```
+
 ```{julia}
 VarCorr(m1)
 ```
@@ -173,10 +191,17 @@ ridgeplot(samp; show_intercept=false)
 Let's consider the classic dysetuff dataset:
 
 ```{julia}
+#| eval: false
 dyestuff = dataset(:dyestuff)
 mdye = fit(MixedModel, @formula(yield ~ 1 + (1 | batch)), dyestuff)
 ```
 
+```{julia}
+#| include: false
+dyestuff = dataset(:dyestuff)
+mdye = fit_or_restore("bootstrap_mdye.json", MixedModel, @formula(yield ~ 1 + (1 | batch)), dyestuff)
+```
+
 
 ```{julia}
 #| warning: false
@@ -204,11 +229,20 @@ However, it is not as straightforward to detect singularity in vector-valued ran
 
 For example, if we bootstrap a model fit to the `sleepstudy` data
 ```{julia}
+#| eval: false
 sleepstudy = dataset(:sleepstudy)
 msleep = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)),
              sleepstudy)
 ```
 
+```{julia}
+#| include: false
+sleepstudy = dataset(:sleepstudy)
+msleep = fit_or_restore("bootstrap_msleep.json", MixedModel,
+                        @formula(reaction ~ 1 + days + (1 + days | subj)),
+                        sleepstudy)
+```
+
 ```{julia}
 sampsleep = parametricbootstrap(MersenneTwister(666), 10_000, msleep)
 tblsleep = sampsleep.tbl
diff --git a/contrasts_fggk21.qmd b/contrasts_fggk21.qmd
index ae6dfe5..1b997ec 100644
--- a/contrasts_fggk21.qmd
+++ b/contrasts_fggk21.qmd
@@ -407,7 +407,8 @@ VarCorr(m_zcp_1_PC_2)
 ```
 
 ```{julia}
-likelihoodratiotest(m_zcp_1_PC_2, m_cpx_1_PC_2)
+# XXX LRT is non monotonic in model complexity!
+# likelihoodratiotest(m_zcp_1_PC_2, m_cpx_1_PC_2)
 ```
 
 # Other topics
diff --git a/contrasts_kwdyz11.qmd b/contrasts_kwdyz11.qmd
index 38caa0b..d2f44ad 100644
--- a/contrasts_kwdyz11.qmd
+++ b/contrasts_kwdyz11.qmd
@@ -15,7 +15,7 @@ using CairoMakie
 using Chain
 using DataFrames
 using MixedModels
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 using StatsBase
 using StatsModels
 
@@ -89,6 +89,7 @@ The `SeqDiffCoding` contrast corresponds to `MASS::contr.sdif()` in R.
 The assignment of random factors such as `Subj` to `Grouping()` is necessary when the sample size is very large. We recommend to include it always, but in this tutorial we do so only in the first example.
 
 ```{julia}
+#| eval: false
 m1 = let levels = ["val", "sod", "dos", "dod"]
   contrasts = Dict(
     :CTR => SeqDiffCoding(; levels),
@@ -97,6 +98,16 @@ m1 = let levels = ["val", "sod", "dos", "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m1 = let levels = ["val", "sod", "dos", "dod"]
+  contrasts = Dict(
+    :CTR => SeqDiffCoding(; levels),
+  )
+  fit_or_restore("contrasts_kwdyz11_m1.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 
 What does the intercept represent?
 
@@ -112,6 +123,7 @@ Grand Mean is mean of condition means.
 `HypothesisCoding` is the most general option available. We can implement all "canned" contrasts ourselves. The next example reproduces the test statistics from `SeqDiffCoding` - with a minor modification illustrating the flexibility of going beyond the default version.
 
 ```{julia}
+#| eval: false
 m1b = let levels = ["val", "sod", "dos", "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
@@ -128,6 +140,24 @@ m1b = let levels = ["val", "sod", "dos", "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m1b = let levels = ["val", "sod", "dos", "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+        -1  1 0  0
+         0 -1 1  0
+         0  0 1 -1
+      ];
+      levels,
+      labels=["spt", "obj", "grv"],
+    ),
+  )
+  fit_or_restore("contrasts_kwdyz11_m1b.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 The difference to the preprogrammed `SeqDiffCoding` is that for the third contrast we changed the direction of the contrast such that the sign of the effect is positive when the result is in agreement with theoretical expectation, that is we subtract the fourth level from the third, not the third level from the fourth.
 
 ## DummyCoding
@@ -135,12 +165,21 @@ The difference to the preprogrammed `SeqDiffCoding` is that for the third contra
 This contrast corresponds to `contr.treatment()` in R
 
 ```{julia}
+#| eval: false
 m2 = let
   contrasts = Dict(:CTR => DummyCoding(; base="val"))
   fit(MixedModel, form, dat1; contrasts, progress)
 end
 ```
 
+```{julia}
+#| include: false
+m2 = let
+  contrasts = Dict(:CTR => DummyCoding(; base="val"))
+  fit_or_restore("contrasts_kwdyz11_m2.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 The `DummyCoding` contrast has the disadvantage that the intercept returns the mean of the level specified as `base`, default is the first level, not the GM.
 
 ## YchycaeitCoding
@@ -149,6 +188,7 @@ The contrasts returned by `DummyCoding` may be exactly what we want.
 Can't we have them, but also have the intercept estimate the GM, rather than the mean of the base level? Yes, we can! We call this "You can have your cake and it eat, too"-Coding (YchycaeitCoding). And we use `HypothesisCoding` to achieve this outcome.
 
 ```{julia}
+#| eval: false
 m2b = let levels = ["val", "sod", "dos", "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
@@ -165,9 +205,28 @@ m2b = let levels = ["val", "sod", "dos", "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m2b = let levels = ["val", "sod", "dos", "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+        -1 1 0 0
+        -1 0 1 0
+        -1 0 0 1
+      ];
+      levels,
+      labels=levels[2:end],
+    )
+  )
+  fit_or_restore("contrasts_kwdyz11_m2b.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 We can simply move the column with -1s for a different base.
 
 ```{julia}
+#| eval: false
 m2c = let levels = ["val", "sod", "dos", "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
@@ -184,9 +243,28 @@ m2c = let levels = ["val", "sod", "dos", "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m2c = let levels = ["val", "sod", "dos", "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+       1 -1  0  0
+       0 -1  1  0
+       0 -1  0  1
+      ];
+      levels,
+      labels=["val", "dos", "dod"],
+    )
+  )
+  fit_or_restore("contrasts_kwdyz11_m2c.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 We can simply relevel the factor with a different base.
 
 ```{julia}
+#| eval: false
 m2d = let levels = ["sod", "val", "dos", "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
@@ -203,29 +281,66 @@ m2d = let levels = ["sod", "val", "dos", "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m2d = let levels = ["sod", "val", "dos", "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+        -1 1 0 0
+        -1 0 1 0
+        -1 0 0 1
+      ];
+      levels,
+      labels=levels[2:end],
+    )
+  )
+  fit_or_restore("contrasts_kwdyz11_m2d.json", MixedModel, form, dat1; contrasts)
+end
+```
+
 ## EffectsCoding
 
 `EffectsCoding` estimates the difference between the Grand Mean and three of the four levels. The difference of the fourth levels can be computed from the Grand Mean and these three differences. 
 
 ```{julia}
+#| eval: false
 m3 = let levels = ["val", "sod", "dos",  "dod"]
   contrasts = Dict(:CTR => EffectsCoding(; levels, base="val"))
   fit(MixedModel, form, dat1; contrasts, progress)
 end
 ```
 
+```{julia}
+#| include: false
+m3 = let levels = ["val", "sod", "dos",  "dod"]
+  contrasts = Dict(:CTR => EffectsCoding(; levels, base="val"))
+  fit_or_restore("contrasts_kwdyz11_m3.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 This contrast corresponds almost to `contr.sum()` in R. The "almost" qualification refers to the fact that `EffectsCoding` uses the first level as default base; `contr.sum()` uses the last factor level.
 
 ```{julia}
+#| eval: false
 m3b = let levels = ["val", "sod", "dos",  "dod"]
   contrasts = Dict(:CTR => EffectsCoding(; levels, base = "dod"))
   fit(MixedModel, form, dat1; contrasts, progress)
 end
 ```
 
+```{julia}
+#| include: false
+m3b = let levels = ["val", "sod", "dos",  "dod"]
+  contrasts = Dict(:CTR => EffectsCoding(; levels, base = "dod"))
+  fit_or_restore("contrasts_kwdyz11_m3b.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 How could we achieve the default result with HypothesisCoding?
 
 ```{julia}
+#| eval: false
 m3c = let levels = ["val", "sod", "dos",  "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
@@ -242,18 +357,45 @@ m3c = let levels = ["val", "sod", "dos",  "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m3c = let levels = ["val", "sod", "dos",  "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+         -1/4   3/4 -1/4  -1/4   # b - GM = b - (a+b+c+d)/4 =>  -1/4*a + 3/4*b - 1/4*c - 1/4*d
+         -1/4  -1/4  3/4  -1/4   # c - GM = c - (a+b+c+d)/4 =>  -1/4*a - 1/4*b + 3/4*c - 1/4*d
+         -1/4  -1/4 -1/4   3/4   # d - GM = d - (a+b+c+d)/4 =>  -1/4*a - 1/4*b - 1/4*c + 3/4*d
+      ];
+      levels,
+      labels=levels[2:end],
+    )
+  )
+  fit_or_restore("contrasts_kwdyz11_m3c.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 
 ## HelmertCoding
 
 `HelmertCoding` codes each level as the difference from the average of the lower levels. With the default order of `CTR` levels we get the following test statistics. These contrasts are orthogonal.
 
 ```{julia}
+#| eval: false
 m4 = let levels = ["val", "sod", "dos",  "dod"]
   contrasts = Dict(:CTR => HelmertCoding(; levels))
   fit(MixedModel, form, dat1; contrasts, progress)
 end
 ```
 
+```{julia}
+#| include: false
+m4 = let levels = ["val", "sod", "dos",  "dod"]
+  contrasts = Dict(:CTR => HelmertCoding(; levels))
+  fit_or_restore("contrasts_kwdyz11_m4.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 ```sh
 + HeC1: (b - a)/2           # (391 - 358)/2 = 16.5
 + HeC2: (c - (b+a)/2)/3     # (405 - (391 + 358)/2)/3 = 10.17 
@@ -263,14 +405,15 @@ end
 We can reconstruct the estimates, but they are scaled by the number of levels involved. With `HypothesisCoding` we can estimate the "unscaled" differences. Also the labeling of the contrasts is not as informative as they could be.
 
 ```{julia}
+#| eval: false
 m4b = let levels = ["val", "sod", "dos",  "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
       [
-          -1    1    0   0 
+          -1    1    0   0
          -1/2 -1/2   1   0
          -1/3 -1/3 -1/3  1
-        
+
       ];
       levels,
       labels= ["2-1", "3-21", "4-321"]
@@ -280,17 +423,45 @@ m4b = let levels = ["val", "sod", "dos",  "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m4b = let levels = ["val", "sod", "dos",  "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+          -1    1    0   0
+         -1/2 -1/2   1   0
+         -1/3 -1/3 -1/3  1
+
+      ];
+      levels,
+      labels= ["2-1", "3-21", "4-321"]
+    )
+  )
+  fit_or_restore("contrasts_kwdyz11_m4b.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 ## Reverse HelmertCoding
 
 `Reverse HelmertCoding` codes each level as the difference from the average of the higher levels. To estimate these effects we simply reverse the order of factor levels. Of course, the contrasts are also orthogonal.
 
 ```{julia}
+#| eval: false
 m4c = let levels = reverse(StatsModels.levels(dat1.CTR))
   contrasts = Dict(:CTR => HelmertCoding(; levels))
   fit(MixedModel, form, dat1; contrasts, progress)
 end
 ```
 
+```{julia}
+#| include: false
+m4c = let levels = reverse(StatsModels.levels(dat1.CTR))
+  contrasts = Dict(:CTR => HelmertCoding(; levels))
+  fit_or_restore("contrasts_kwdyz11_m4c.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 ```sh
 + HeC1:(c - d)/2            # (405 - 402)/2 = 1.5
 + HeC2:(b - (c+d)/2)/3      # (391 - (405 + 402)/2)/3 = -4.17
@@ -300,11 +471,12 @@ end
 ... and the unscaled-by-number-of-levels estimates.
 
 ```{julia}
+#| eval: false
 m4d = let levels = ["val", "sod", "dos",  "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
       [
-        0    0     1   -1 
+        0    0     1   -1
         0    1   -1/2 -1/2
         1  -1/3  -1/3 -1/3
       ];
@@ -316,6 +488,24 @@ m4d = let levels = ["val", "sod", "dos",  "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m4d = let levels = ["val", "sod", "dos",  "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+        0    0     1   -1
+        0    1   -1/2 -1/2
+        1  -1/3  -1/3 -1/3
+      ];
+      levels,
+      labels= ["3-4", "2-34", "1-234"]
+    )
+  )
+  fit_or_restore("contrasts_kwdyz11_m4d.json", MixedModel, form, dat1; contrasts)
+end
+```
+
 # Other orthogonal contrasts
 
 For factors with more than four levels there are many options for specifying orthogonal contrasts as long as one proceeds in a top-down strictly hierarchical fashion.
@@ -362,6 +552,7 @@ In a figure With A on the x-axis and the levels of B shown as two lines, the int
 A positive coefficient implies overadditivity (diverging lines toward the right) and a negative coefficient underadditivity (converging lines).
 
 ```{julia}
+#| eval: false
 m5 = let levels = ["val", "sod", "dos", "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
@@ -378,6 +569,24 @@ m5 = let levels = ["val", "sod", "dos", "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m5 = let levels = ["val", "sod", "dos", "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+        -1 -1 +1 +1          # A
+        -1 +1 -1 +1          # B
+        +1 -1 -1 +1          # A x B
+      ];
+      levels,
+      labels=["A", "B", "AxB"],
+    ),
+  )
+  fit_or_restore("contrasts_kwdyz11_m5.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 It is also helpful to see the corresponding layout of the four means for the interaction of A and B (i.e., the third contrast)
 
 ```
@@ -408,7 +617,7 @@ The following contrast specification returns an estimate for the main effect of
 In a figure With A on the x-axis and the levels of B shown as two lines, the second contrast tests whether A1-B1 is different from A1-B2 and the third contrast tests whether A2-B1 is different from A2-B2.
 
 ```{julia}
-
+#| eval: false
 m6 = let levels = ["val", "sod", "dos", "dod"]
   contrasts = Dict(
     :CTR => HypothesisCoding(
@@ -425,6 +634,24 @@ m6 = let levels = ["val", "sod", "dos", "dod"]
 end
 ```
 
+```{julia}
+#| include: false
+m6 = let levels = ["val", "sod", "dos", "dod"]
+  contrasts = Dict(
+    :CTR => HypothesisCoding(
+      [
+        -1 -1 +1 +1
+        -1 +1  0  0
+         0  0 +1 -1
+      ];
+      levels,
+      labels=["do_so", "spt", "grv"],
+    ),
+  )
+  fit_or_restore("contrasts_kwdyz11_m6.json", MixedModel, form, dat1; contrasts, progress)
+end
+```
+
 The three contrasts for one main effect and two nested contrasts are orthogonal. There is no test of the interaction (parallelism).
 
 # An Example for a complex example: Group(2) x A(2) x B(2)
@@ -451,17 +678,28 @@ cntrst2 = Dict(
   )
 
 f6_cpx = @formula dv ~ 1 + Group*A*B   + (1 + A+B | Subj);
+f6_zcp = @formula dv ~ 1 + Group*A*B   + zerocorr(1 + A+B | Subj);
+f6_ovi = @formula dv ~ 1 + Group*A*B   +  (1 | Subj);
+```
+
+```{julia}
+#| eval: false
 m6_cpx = fit(MixedModel, f6_cpx, dat2; contrasts=cntrst2)
 issingular(m6_cpx)
 
-f6_zcp = @formula dv ~ 1 + Group*A*B   + zerocorr(1 + A+B | Subj);
 m6_zcp = fit(MixedModel, f6_zcp, dat2; contrasts=cntrst2)
 issingular(m6_zcp)
 
-f6_ovi = @formula dv ~ 1 + Group*A*B   +  (1 | Subj);
 m6_ovi = fit(MixedModel, f6_ovi, dat2; contrasts=cntrst2)
 ```
 
+```{julia}
+#| include: false
+m6_cpx = fit_or_restore("contrasts_kwdyz11_m6_cpx.json", MixedModel, f6_cpx, dat2; contrasts=cntrst2)
+m6_zcp = fit_or_restore("contrasts_kwdyz11_m6_zcp.json", MixedModel, f6_zcp, dat2; contrasts=cntrst2)
+m6_ovi = fit_or_restore("contrasts_kwdyz11_m6_ovi.json", MixedModel, f6_ovi, dat2; contrasts=cntrst2)
+```
+
 ```{julia}
 lrtest(m6_ovi, m6_zcp, m6_cpx)
 ```
diff --git a/fggk21.qmd b/fggk21.qmd
index cde3a8a..c5d09d5 100644
--- a/fggk21.qmd
+++ b/fggk21.qmd
@@ -34,7 +34,7 @@ using MixedModelsMakie: simplelinreg
 using Random
 using Statistics
 using StatsBase
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 
 const progress=false
 ```
@@ -287,12 +287,21 @@ In its random-effect structure (RES) we only vary intercepts (i.e., Grand Means)
 It is well known that such a simple RES is likely to be anti-conservative with respect to fixed-effect test statistics.
 
 ```{julia}
+#| eval: false
 m_ovi = let
   f = @formula zScore ~ 1 + Test * Sex * a1 + (1 | School)
   fit(MixedModel, f, dat; contrasts, progress)
 end
 ```
 
+```{julia}
+#| include: false
+m_ovi = let
+  f = @formula zScore ~ 1 + Test * Sex * a1 + (1 | School)
+  fit_or_restore("fggk21_m_ovi.json", MixedModel, f, dat; contrasts, progress)
+end
+```
+
 Is the model singular (overparameterized, degenerate)?
 In other words: Is the model not supported by the data?
 
@@ -309,6 +318,7 @@ In this LMM we allow that schools differ not only in `GM`, but also in the size
 We assume that there is covariance associated with these CPs beyond residual noise, that is we assume that there is no detectable evidence in the data that the CPs are different from zero.
 
 ```{julia}
+#| eval: false
 m_zcp = let
   f = @formula(
     zScore ~
@@ -318,6 +328,17 @@ m_zcp = let
 end
 ```
 
+```{julia}
+#| include: false
+m_zcp = let
+  f = @formula(
+    zScore ~
+      1 + Test * Sex * a1 + zerocorr(1 + Test + Sex + a1 | School)
+  )
+  fit_or_restore("fggk21_m_zcp.json", MixedModel, f, dat; contrasts, progress)
+end
+```
+
 Depending on sampling, this model estimating variance components for `School`  may or may not be supported by the data.
 
 ```{julia}
@@ -329,6 +350,7 @@ issingular(m_zcp)
 In the complex LMM investigated in this sequence we give up the assumption of zero-correlation between VCs.
 
 ```{julia}
+#| eval: false
 m_cpx = let
   f = @formula(
     zScore ~ 1 + Test * Sex * a1 + (1 + Test + Sex + a1 | School)
@@ -337,6 +359,16 @@ m_cpx = let
 end
 ```
 
+```{julia}
+#| include: false
+m_cpx = let
+  f = @formula(
+    zScore ~ 1 + Test * Sex * a1 + (1 + Test + Sex + a1 | School)
+  )
+  fit_or_restore("fggk21_m_cpx.json", MixedModel, f, dat; contrasts, progress)
+end
+```
+
 We also need to see the VCs and CPs of the random-effect structure (RES).
 
 ```{julia}
@@ -438,12 +470,21 @@ As an illustration, we fit a full CP matrix for the `Cohort`.
 As there are only nine cohorts in the data, we may be asking too much to estimate 5*6/2 = 15 VC/CP parameters.
 
 ```{julia}
+#| eval: false
 m_cpxCohort = let
   f = @formula zScore ~ 1 + Test * a1 * Sex + (1 + Test | Cohort)
   fit(MixedModel, f, dat; contrasts, progress)
 end
 ```
 
+```{julia}
+#| include: false
+m_cpxCohort = let
+  f = @formula zScore ~ 1 + Test * a1 * Sex + (1 + Test | Cohort)
+  fit_or_restore("fggk21_m_cpxCohort.json", MixedModel, f, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 VarCorr(m_cpxCohort)
 ```
@@ -456,6 +497,7 @@ The model is overparameterized with several CPs estimated between |.98| and |1.0
 How about the **zero-correlation parameter** (zcp) version of this LMM?
 
 ```{julia}
+#| eval: false
 m_zcpCohort = let
   f = @formula(
     zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Cohort)
@@ -464,6 +506,16 @@ m_zcpCohort = let
 end
 ```
 
+```{julia}
+#| include: false
+m_zcpCohort = let
+  f = @formula(
+    zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Cohort)
+  )
+  fit_or_restore("fggk21_m_zcpCohort.json", MixedModel, f, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 issingular(m_zcpCohort)
 ```
@@ -472,12 +524,21 @@ This `zcpLMM`  is also singular. Three of the five VCs are estimated as zero.
 This raises the possibility that LMM `m_oviCohort`  might fit as well as LMM `m_zcpCohort`.
 
 ```{julia}
+#| eval: false
 m_oviCohort = let
   f = @formula zScore ~ 1 + Test * a1 * Sex + (1 | Cohort)
   fit(MixedModel, f, dat; contrasts, progress)
 end
 ```
 
+```{julia}
+#| include: false
+m_oviCohort = let
+  f = @formula zScore ~ 1 + Test * a1 * Sex + (1 | Cohort)
+  fit_or_restore("fggk21_m_oviCohort.json", MixedModel, f, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 issingular(m_oviCohort)
 ```
diff --git a/fits/bootstrap_m0.json.zip b/fits/bootstrap_m0.json.zip
new file mode 100644
index 0000000..56f08b0
--- /dev/null
+++ b/fits/bootstrap_m0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3c886179fe1f23d009d6a9cbdbbd19d329903ed9cdeed3cde657a56cea6bb520
+size 148842
diff --git a/fits/bootstrap_m1.json.zip b/fits/bootstrap_m1.json.zip
new file mode 100644
index 0000000..921b18f
--- /dev/null
+++ b/fits/bootstrap_m1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8443413cc5d12636f73942b690b05d9f74c664a87d525fa43dc3d9415aab1663
+size 4529
diff --git a/fits/bootstrap_mdye.json.zip b/fits/bootstrap_mdye.json.zip
new file mode 100644
index 0000000..e4b32d8
--- /dev/null
+++ b/fits/bootstrap_mdye.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c4f679871514817e0cd0d909626c42a166e567b3183f648170074f29cb500b01
+size 651
diff --git a/fits/bootstrap_msleep.json.zip b/fits/bootstrap_msleep.json.zip
new file mode 100644
index 0000000..2abfda7
--- /dev/null
+++ b/fits/bootstrap_msleep.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fa657f5b28250fb16711f72840a6bbbb999a85660f5b6ad1360911ea83dbe804
+size 3071
diff --git a/fits/contrasts_kwdyz11_m1.json.zip b/fits/contrasts_kwdyz11_m1.json.zip
new file mode 100644
index 0000000..d80227f
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:597a694906259dff46bc4a569b6334c73d5ee8d9fb721ee23af2be91277ef445
+size 28665
diff --git a/fits/contrasts_kwdyz11_m1b.json.zip b/fits/contrasts_kwdyz11_m1b.json.zip
new file mode 100644
index 0000000..69f4074
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m1b.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:734dfa4e0be4eea5c4f180f363513f153435fb02d9c4c92d699728e67da0efbf
+size 26848
diff --git a/fits/contrasts_kwdyz11_m2.json.zip b/fits/contrasts_kwdyz11_m2.json.zip
new file mode 100644
index 0000000..57f2ed2
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ab324e302c03f598f8805d793defed5d02a2edbc1f31dd224903838780b8c1c0
+size 37124
diff --git a/fits/contrasts_kwdyz11_m2b.json.zip b/fits/contrasts_kwdyz11_m2b.json.zip
new file mode 100644
index 0000000..8e21428
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m2b.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1d9d6fa7c07d2b45b6576e8ebc2a7e59a142b5f81c9ba1b1391378bbaa020a61
+size 32202
diff --git a/fits/contrasts_kwdyz11_m2c.json.zip b/fits/contrasts_kwdyz11_m2c.json.zip
new file mode 100644
index 0000000..bea8f43
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m2c.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:562bc151c1fb97ad002e2a5a4b14e156ae762a903c18ef285bbf5567e7762cf6
+size 24219
diff --git a/fits/contrasts_kwdyz11_m2d.json.zip b/fits/contrasts_kwdyz11_m2d.json.zip
new file mode 100644
index 0000000..45995a7
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m2d.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d426e56bd4f85bc900c0e27438c195e5546bb50ae653c5f3c9455c711a50a297
+size 27517
diff --git a/fits/contrasts_kwdyz11_m3.json.zip b/fits/contrasts_kwdyz11_m3.json.zip
new file mode 100644
index 0000000..d8952ee
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m3.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bd4ce20ed7991a756c137dd1db86f9d95c4c503789ebb46dc5565aa7aed46e06
+size 18099
diff --git a/fits/contrasts_kwdyz11_m3b.json.zip b/fits/contrasts_kwdyz11_m3b.json.zip
new file mode 100644
index 0000000..ff7e0d1
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m3b.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:295ff5a768f747833570a1462dd9ac331ffc452c78b12cddf6b5ca04b3f57319
+size 21855
diff --git a/fits/contrasts_kwdyz11_m3c.json.zip b/fits/contrasts_kwdyz11_m3c.json.zip
new file mode 100644
index 0000000..d46ccd7
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m3c.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c09e90d44e455910dfdfbfe8c019e71310e5a6612ab389e999c3548e195e4b07
+size 25154
diff --git a/fits/contrasts_kwdyz11_m4.json.zip b/fits/contrasts_kwdyz11_m4.json.zip
new file mode 100644
index 0000000..0e4dc26
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m4.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:edb2acdab6f843516c6bbcbbf502c063aec7ed4dc6abc2d0175c512d358135d3
+size 49551
diff --git a/fits/contrasts_kwdyz11_m4b.json.zip b/fits/contrasts_kwdyz11_m4b.json.zip
new file mode 100644
index 0000000..0676991
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m4b.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ceb9c62ed4be7e000103a84c0534cc2e80c323f7c24505ddf88049910ab08822
+size 26971
diff --git a/fits/contrasts_kwdyz11_m4c.json.zip b/fits/contrasts_kwdyz11_m4c.json.zip
new file mode 100644
index 0000000..df93cd3
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m4c.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:21c16b2e324c984d608a4fa8b02515fbf49b00323803cd222c9d5aaf2dd466e1
+size 43089
diff --git a/fits/contrasts_kwdyz11_m4d.json.zip b/fits/contrasts_kwdyz11_m4d.json.zip
new file mode 100644
index 0000000..1528ff0
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m4d.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9075f5467115022e2f3fd677fcf020644aae5e97fba1b2e02cae477725e8d380
+size 22737
diff --git a/fits/contrasts_kwdyz11_m5.json.zip b/fits/contrasts_kwdyz11_m5.json.zip
new file mode 100644
index 0000000..dc50130
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m5.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bf425cf1433f9c2330d764c7505bff96f4e7a87a34226cb34c976b04c1f60263
+size 18090
diff --git a/fits/contrasts_kwdyz11_m6.json.zip b/fits/contrasts_kwdyz11_m6.json.zip
new file mode 100644
index 0000000..142ed67
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m6.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d35f49e4df733fbc24c8f71ab0d4e4917488858229fd725a5a2f8aa6f7f94f3a
+size 22340
diff --git a/fits/contrasts_kwdyz11_m6_cpx.json.zip b/fits/contrasts_kwdyz11_m6_cpx.json.zip
new file mode 100644
index 0000000..fc57b1a
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m6_cpx.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2ba5270c9d275ce04f6e6565e97a5a5c63b86ab4661548699846d56168182424
+size 19982
diff --git a/fits/contrasts_kwdyz11_m6_ovi.json.zip b/fits/contrasts_kwdyz11_m6_ovi.json.zip
new file mode 100644
index 0000000..a9d9ff7
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m6_ovi.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fcb6a68b4c8e51b976a38ab594d7e682ee404e235b0410d4db5ed0b797fe9636
+size 650
diff --git a/fits/contrasts_kwdyz11_m6_zcp.json.zip b/fits/contrasts_kwdyz11_m6_zcp.json.zip
new file mode 100644
index 0000000..5c71185
--- /dev/null
+++ b/fits/contrasts_kwdyz11_m6_zcp.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fcebdfb465ba42c72100a62672a876b98cd0adc70327a9b33bd5dfb95520baf0
+size 3568
diff --git a/fits/fggk21_m_cpx.json.zip b/fits/fggk21_m_cpx.json.zip
new file mode 100644
index 0000000..075b3c8
--- /dev/null
+++ b/fits/fggk21_m_cpx.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f7fe50c219f6db3b634b04924103f11bf8baf4d8a3bb01e141b941dd45df2641
+size 751140
diff --git a/fits/fggk21_m_cpxCohort.json.zip b/fits/fggk21_m_cpxCohort.json.zip
new file mode 100644
index 0000000..aae48aa
--- /dev/null
+++ b/fits/fggk21_m_cpxCohort.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6257c784c577266fcb28d3a22a71bcf3e88226aaae929ce4e0a8bc7e634c21db
+size 101416
diff --git a/fits/fggk21_m_ovi.json.zip b/fits/fggk21_m_ovi.json.zip
new file mode 100644
index 0000000..5c41347
--- /dev/null
+++ b/fits/fggk21_m_ovi.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0f398a0ebfb86e85b86daccd87bf2a819f3af2ba03f51a17e41ba6dd59cae3d8
+size 749
diff --git a/fits/fggk21_m_oviCohort.json.zip b/fits/fggk21_m_oviCohort.json.zip
new file mode 100644
index 0000000..2d60952
--- /dev/null
+++ b/fits/fggk21_m_oviCohort.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:27482f2a9f01f626e3d0bb24c6db890f0971dfb9d4aba91b3a89c40f8ce7c3f8
+size 648
diff --git a/fits/fggk21_m_zcp.json.zip b/fits/fggk21_m_zcp.json.zip
new file mode 100644
index 0000000..6d1784a
--- /dev/null
+++ b/fits/fggk21_m_zcp.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8472d8378cdf6786b7d8bf030b1238e30328b4f91222b7d17e552e7d842827b9
+size 8575
diff --git a/fits/fggk21_m_zcpCohort.json.zip b/fits/fggk21_m_zcpCohort.json.zip
new file mode 100644
index 0000000..98f1e85
--- /dev/null
+++ b/fits/fggk21_m_zcpCohort.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a4ee2b6b195b428eed421e6f8607d3195736d2e1073d7c33ac44a096e727ac32
+size 9928
diff --git a/fits/glmm_gm1.json.zip b/fits/glmm_gm1.json.zip
new file mode 100644
index 0000000..517f489
--- /dev/null
+++ b/fits/glmm_gm1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fd62ed53c2f58bf05f236023c3592ed7d368eba5ef3a6abcb303d55384c539e7
+size 92159
diff --git a/fits/glmm_gm2.json.zip b/fits/glmm_gm2.json.zip
new file mode 100644
index 0000000..b2dec71
--- /dev/null
+++ b/fits/glmm_gm2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:141980d8a3d6aee415fd6efa60bd94f636157960ccd3107230fa03b772e5f73e
+size 75904
diff --git a/fits/glmm_gm3.json.zip b/fits/glmm_gm3.json.zip
new file mode 100644
index 0000000..92e0b93
--- /dev/null
+++ b/fits/glmm_gm3.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8ed5cd7b2c7f6f862e8a931079f176baf7353cecc626c8a9bc0491a210a4b1b5
+size 69518
diff --git a/fits/glmm_m1.json.zip b/fits/glmm_m1.json.zip
new file mode 100644
index 0000000..482e0b4
--- /dev/null
+++ b/fits/glmm_m1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8661adb70cd4eeb3233b86ec7aa3cd90312805947f4213057e238ce30b182fb8
+size 3071
diff --git a/fits/kwdyz11_m1.json.zip b/fits/kwdyz11_m1.json.zip
new file mode 100644
index 0000000..069be85
--- /dev/null
+++ b/fits/kwdyz11_m1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:922cd70f85dcada83ef53a2f67eb500898475b04aabd0c1cf682efe1533aa123
+size 27608
diff --git a/fits/largescaledesigned_elm01.json.zip b/fits/largescaledesigned_elm01.json.zip
new file mode 100644
index 0000000..69dd445
--- /dev/null
+++ b/fits/largescaledesigned_elm01.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6bcbe2b24bb2ff1048fb7ec69d97238cd752c8c4beda1dba0cc7a1c9faad09ea
+size 1859
diff --git a/fits/largescaledesigned_elm02.json.zip b/fits/largescaledesigned_elm02.json.zip
new file mode 100644
index 0000000..8bc7c38
--- /dev/null
+++ b/fits/largescaledesigned_elm02.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:966473968a9bb14a0b6fe9d521409ec2bc1379e69aefbcda7f9e2c0cded29475
+size 2021
diff --git a/fits/mrk17_m_cpx.json.zip b/fits/mrk17_m_cpx.json.zip
new file mode 100644
index 0000000..2bc42c8
--- /dev/null
+++ b/fits/mrk17_m_cpx.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:34243da587e0b03dff8823c4b795475cdb540964fa3511f5068e2613771907fb
+size 651633
diff --git a/fits/mrk17_m_mrk17_crossed.json.zip b/fits/mrk17_m_mrk17_crossed.json.zip
new file mode 100644
index 0000000..b6c5ab8
--- /dev/null
+++ b/fits/mrk17_m_mrk17_crossed.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:063a6ecee3953a24e302e7ff2cbc173654beb1784b5c4d5b6f9d89c0877c7fe5
+size 9088
diff --git a/fits/mrk17_m_mrk17_nested.json.zip b/fits/mrk17_m_mrk17_nested.json.zip
new file mode 100644
index 0000000..c163a18
--- /dev/null
+++ b/fits/mrk17_m_mrk17_nested.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4c8e3fbf95c5bbd85d577ef731082403ede29ec450689c126875be3cbff7fcfb
+size 12207
diff --git a/fits/mrk17_m_prm.json.zip b/fits/mrk17_m_prm.json.zip
new file mode 100644
index 0000000..91ebb94
--- /dev/null
+++ b/fits/mrk17_m_prm.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b33f814f7ab36333d29df6838abc8e605bcfbfbf91c8df196303b00bce002ac2
+size 15108
diff --git a/fits/mrk17_m_prm3.json.zip b/fits/mrk17_m_prm3.json.zip
new file mode 100644
index 0000000..5acd107
--- /dev/null
+++ b/fits/mrk17_m_prm3.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3cdd19579e3da2111c40681df42db71e596033c8ca8ef6d19bf69584f5378111
+size 11982
diff --git a/fits/mrk17_m_prm3_posthoc.json.zip b/fits/mrk17_m_prm3_posthoc.json.zip
new file mode 100644
index 0000000..475ca44
--- /dev/null
+++ b/fits/mrk17_m_prm3_posthoc.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:302d817208cccbc46f892771689828014d6b4cb3336aa7acd9d556625a45eac1
+size 13463
diff --git a/fits/mrk17_m_prm_1.json.zip b/fits/mrk17_m_prm_1.json.zip
new file mode 100644
index 0000000..9592c44
--- /dev/null
+++ b/fits/mrk17_m_prm_1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4ecd2de3ae3d4790aa3240144ce775116aae9e27dfabdd3c02bf358c0a8b726e
+size 15124
diff --git a/fits/mrk17_m_prm_2.json.zip b/fits/mrk17_m_prm_2.json.zip
new file mode 100644
index 0000000..56fe3d6
--- /dev/null
+++ b/fits/mrk17_m_prm_2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ab6b8de0b146adab7d95a4d524b0c979ceaa876f44799ce80b604e76c3ed857a
+size 17425
diff --git a/fits/mrk17_m_prm_3.json.zip b/fits/mrk17_m_prm_3.json.zip
new file mode 100644
index 0000000..89cc64f
--- /dev/null
+++ b/fits/mrk17_m_prm_3.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9b4d7fec1d93c80e3fb34bc0fd7c94fba986e89279a2cebf3d58fb7f7e58e823
+size 14989
diff --git a/fits/mrk17_m_prm_4.json.zip b/fits/mrk17_m_prm_4.json.zip
new file mode 100644
index 0000000..3a99e97
--- /dev/null
+++ b/fits/mrk17_m_prm_4.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:60e835ef721f2dbfa3db0c5c41a85de28a3eb37273de5e3e50fb85a4614b0f76
+size 14324
diff --git a/fits/mrk17_m_prm_5.json.zip b/fits/mrk17_m_prm_5.json.zip
new file mode 100644
index 0000000..814b5e8
--- /dev/null
+++ b/fits/mrk17_m_prm_5.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6a0020e799d21e65ca600161304146c80fa0c9276f88f4193d5c8b722568f234
+size 15108
diff --git a/fits/mrk17_m_zcp.json.zip b/fits/mrk17_m_zcp.json.zip
new file mode 100644
index 0000000..e2267b0
--- /dev/null
+++ b/fits/mrk17_m_zcp.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:22c67ee27eded51588988aef6a154994ee7fdbe8b8f7c2bb8fd6d57b07744d8e
+size 63067
diff --git a/fits/mrk17_m_zcp2.json.zip b/fits/mrk17_m_zcp2.json.zip
new file mode 100644
index 0000000..26cd6ba
--- /dev/null
+++ b/fits/mrk17_m_zcp2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0bfd5de1cf1feac7fc4ed17252c5ef2018811766330a90295cf0746b5d194cd9
+size 22863
diff --git a/fits/mrk17_m_zcp3.json.zip b/fits/mrk17_m_zcp3.json.zip
new file mode 100644
index 0000000..719eaab
--- /dev/null
+++ b/fits/mrk17_m_zcp3.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e80ec8e690f33198dc0e04dc9776618ca9411a22f70d1d59e41b32b0c085308a
+size 5692
diff --git a/fits/mrk17_m_zcp4.json.zip b/fits/mrk17_m_zcp4.json.zip
new file mode 100644
index 0000000..a491b82
--- /dev/null
+++ b/fits/mrk17_m_zcp4.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8d4384d00055801922636d1c55a4caf32e2955a8c839560091706fc497213505
+size 3654
diff --git a/fits/partial_within_m0.json.zip b/fits/partial_within_m0.json.zip
new file mode 100644
index 0000000..9d635fb
--- /dev/null
+++ b/fits/partial_within_m0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:eef25195246d8cd5d2f740b5c7bbb20a1e63e966e798efaf4b44be7abd04a816
+size 4287
diff --git a/fits/partial_within_m1.json.zip b/fits/partial_within_m1.json.zip
new file mode 100644
index 0000000..ab9c42e
--- /dev/null
+++ b/fits/partial_within_m1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:62b350608be138a6e9a2633c2f0e5d98f229f98bd2b6d317277754240367f933
+size 4883
diff --git a/fits/profiling_pr01.json.zip b/fits/profiling_pr01.json.zip
new file mode 100644
index 0000000..4b588fc
--- /dev/null
+++ b/fits/profiling_pr01.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:da61ca28cc25ffe731a7f194cdae96aaf9189017ba55ae08bc603308cbcd964b
+size 4087
diff --git a/fits/shrinkageplot_m1.json.zip b/fits/shrinkageplot_m1.json.zip
new file mode 100644
index 0000000..dc1a93b
--- /dev/null
+++ b/fits/shrinkageplot_m1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f7188d13bc3e84fb17fc86199efe44a1ec49c7ac1b7c7493ec4a0681595346a9
+size 148842
diff --git a/fits/shrinkageplot_m2.json.zip b/fits/shrinkageplot_m2.json.zip
new file mode 100644
index 0000000..1f42aa6
--- /dev/null
+++ b/fits/shrinkageplot_m2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d7b9f067dd1b32c3c96d47b1884e675a34f8becbacb04c8c37cee24d4e6e0b97
+size 4532
diff --git a/fits/shrinkageplot_m3.json.zip b/fits/shrinkageplot_m3.json.zip
new file mode 100644
index 0000000..7aa1964
--- /dev/null
+++ b/fits/shrinkageplot_m3.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a4b81f4a09409a8256b01ba307bed9c4d49396921d64c074937f9928f609b9d7
+size 4312
diff --git a/fits/shrinkageplot_m4.json.zip b/fits/shrinkageplot_m4.json.zip
new file mode 100644
index 0000000..d1ea4df
--- /dev/null
+++ b/fits/shrinkageplot_m4.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:062367a7abb3e60431d4b993bcbcbdca63a8db2a10882ddfbb8d4b163e736789
+size 4312
diff --git a/fits/sleepstudy_m1.json.zip b/fits/sleepstudy_m1.json.zip
new file mode 100644
index 0000000..fc41064
--- /dev/null
+++ b/fits/sleepstudy_m1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:df021e0d991fbc176e2b646aaaab367edbaa5107cf7ece3bd283c862eb77b179
+size 3071
diff --git a/fits/sleepstudy_m2.json.zip b/fits/sleepstudy_m2.json.zip
new file mode 100644
index 0000000..6bdc162
--- /dev/null
+++ b/fits/sleepstudy_m2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a4d53f0fc60cc1a52596fc5e13c87b9b2523a781206a4a839673c7bb41e7579f
+size 1828
diff --git a/fits/sleepstudy_speed_m1.json.zip b/fits/sleepstudy_speed_m1.json.zip
new file mode 100644
index 0000000..6e74137
--- /dev/null
+++ b/fits/sleepstudy_speed_m1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0a3c5d24c194266a3d76f031af5dbdd65b880f08a6cfd963c7f58a87b11ed518
+size 2820
diff --git a/fits/sleepstudy_speed_m2.json.zip b/fits/sleepstudy_speed_m2.json.zip
new file mode 100644
index 0000000..bc96f84
--- /dev/null
+++ b/fits/sleepstudy_speed_m2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a4a51314b1936db2950076b7810654563af653b479c050382b1f1b02fb178beb
+size 1601
diff --git a/fits/transformation_days_centered.json.zip b/fits/transformation_days_centered.json.zip
new file mode 100644
index 0000000..741d799
--- /dev/null
+++ b/fits/transformation_days_centered.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:34e32b8ea5c7436583986e37073dbae33266446391c803d35845e1c2cee034d1
+size 1983
diff --git a/fits/transformation_model_bc.json.zip b/fits/transformation_model_bc.json.zip
new file mode 100644
index 0000000..a1b3be1
--- /dev/null
+++ b/fits/transformation_model_bc.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7bb2894fc4cb129cfbed80dbfe67419bed6676200c4e5688db4b548dbda48f8b
+size 2820
diff --git a/fits/transformation_slp.json.zip b/fits/transformation_slp.json.zip
new file mode 100644
index 0000000..f6375e2
--- /dev/null
+++ b/fits/transformation_slp.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:431b3e612dd7899243883f32e768320a1d5b3101b24b9494ff8234a525bad7de
+size 3071
diff --git a/glmm.qmd b/glmm.qmd
index 1ba72e9..ecb8467 100644
--- a/glmm.qmd
+++ b/glmm.qmd
@@ -19,7 +19,7 @@ using DataFrameMacros
 using DataFrames
 using MixedModels
 using MixedModelsMakie
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 
 const progress=false
 ```
@@ -36,6 +36,7 @@ sleepstudy = DataFrame(dataset(:sleepstudy))
 ```
 
 ```{julia}
+#| eval: false
 contrasts = Dict(:subj => Grouping())
 m1 = let f = @formula reaction ~ 1 + days + (1 + days | subj)
   fit(MixedModel, f, sleepstudy; contrasts, progress)
@@ -43,6 +44,15 @@ end
 println(m1)
 ```
 
+```{julia}
+#| include: false
+contrasts = Dict(:subj => Grouping())
+m1 = let f = @formula reaction ~ 1 + days + (1 + days | subj)
+  fit_or_restore("glmm_m1.json", MixedModel, f, sleepstudy; contrasts, progress)
+end
+println(m1)
+```
+
 The response vector, y, has 180 elements.  The fixed-effects coefficient vector, β, has 2 elements and the fixed-effects model matrix, X, is of size 180 × 2.
 
 ```{julia}
@@ -235,6 +245,7 @@ draw(
 ```
 
 ```{julia}
+#| eval: false
 contrasts = Dict(
   :dist => Grouping(),
   :urban => HelmertCoding(),
@@ -250,6 +261,23 @@ gm1 = let
 end
 ```
 
+```{julia}
+#| include: false
+contrasts = Dict(
+  :dist => Grouping(),
+  :urban => HelmertCoding(),
+  :livch => DummyCoding(), # default, but no harm in being explicit
+)
+nAGQ = 9
+dist = Bernoulli()
+gm1 = let
+  form = @formula(
+    use ~ 1 + age + abs2(age) + urban + livch + (1 | dist)
+  )
+  fit_or_restore("glmm_gm1.json", MixedModel, form, contra, dist; nAGQ, contrasts, progress)
+end
+```
+
 Notice that the linear term for `age` is not significant but the quadratic term for `age` is highly significant.
 We usually retain the lower order term, even if it is not significant, if the higher order term is significant.
 
@@ -276,6 +304,7 @@ contrasts[:children] = EffectsCoding()
 ```
 
 ```{julia}
+#| eval: false
 gm2 = let
   form = @formula(
     use ~
@@ -290,6 +319,22 @@ gm2 = let
 end
 ```
 
+```{julia}
+#| include: false
+gm2 = let
+  form = @formula(
+    use ~
+      1 +
+      age * children +
+      abs2(age) +
+      children +
+      urban +
+      (1 | dist)
+  )
+  fit_or_restore("glmm_gm2.json", MixedModel, form, contra, dist; nAGQ, contrasts, progress)
+end
+```
+
 ```{julia}
 #| code-fold: true
 let
@@ -318,6 +363,7 @@ It turns out that there can be more difference between urban and rural settings
 To model this difference we build a model with `urban&dist` as a grouping factor.
 
 ```{julia}
+#| eval: false
 gm3 = let
   form = @formula(
     use ~
@@ -332,6 +378,22 @@ gm3 = let
 end
 ```
 
+```{julia}
+#| include: false
+gm3 = let
+  form = @formula(
+    use ~
+      1 +
+      age * children +
+      abs2(age) +
+      children +
+      urban +
+      (1 | urban & dist)
+  )
+  fit_or_restore("glmm_gm3.json", MixedModel, form, contra, dist; nAGQ, contrasts, progress)
+end
+```
+
 ```{julia}
 #| code-fold: true
 let
diff --git a/kwdyz11.qmd b/kwdyz11.qmd
index 197221c..beb6d23 100644
--- a/kwdyz11.qmd
+++ b/kwdyz11.qmd
@@ -42,7 +42,7 @@ using Distributions
 using MixedModels
 using MixedModelsMakie
 using Random
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 using StatsBase
 ```
 
@@ -129,6 +129,7 @@ A better alternative to the boxplot is a dotplot. It also displays subjects' con
 # Linear mixed model
 
 ```{julia}
+#| eval: false
 contrasts = Dict(
   :CTR => SeqDiffCoding(; levels=["val", "sod", "dos", "dod"]),
 )
@@ -138,6 +139,17 @@ m1 = let
 end
 ```
 
+```{julia}
+#| include: false
+contrasts = Dict(
+  :CTR => SeqDiffCoding(; levels=["val", "sod", "dos", "dod"]),
+)
+m1 = let
+  form = @formula(log(rt) ~ 1 + CTR + (1 + CTR | Subj))
+  fit_or_restore("kwdyz11_m1.json", MixedModel, form, dat; contrasts)
+end
+```
+
 ```{julia}
 VarCorr(m1)
 ```
diff --git a/largescaledesigned.qmd b/largescaledesigned.qmd
index b9cf67d..8d80f1f 100644
--- a/largescaledesigned.qmd
+++ b/largescaledesigned.qmd
@@ -19,7 +19,7 @@ using DataFrames
 using Effects
 using MixedModels
 using MixedModelsMakie
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 using StandardizedPredictors
 using StatsBase
 ```
@@ -427,6 +427,7 @@ contrasts = Dict(
 and fit a first model with simple, scalar, random effects for `subj` and `item`.
 
 ```{julia}
+#| eval: false
 elm01 = let
   form = @formula(
     1000 / rt ~ 1 + isword * wrdlen + (1 | item) + (1 | subj)
@@ -435,6 +436,16 @@ elm01 = let
 end
 ```
 
+```{julia}
+#| include: false
+elm01 = let
+  form = @formula(
+    1000 / rt ~ 1 + isword * wrdlen + (1 | item) + (1 | subj)
+  )
+  fit_or_restore("largescaledesigned_elm01.json", MixedModel, form, pruned; contrasts, progress=false)
+end
+```
+
 The predicted response speed by word length and word/nonword status can be summarized as
 
 ```{julia}
@@ -445,6 +456,7 @@ effects(Dict(:isword => [false, true], :wrdlen => 4:2:12), elm01)
 If we restrict to only those subjects with 80% accuracy or greater the model becomes
 
 ```{julia}
+#| eval: false
 elm02 = let
   form = @formula(
     1000 / rt ~ 1 + isword * wrdlen + (1 | item) + (1 | subj)
@@ -454,6 +466,17 @@ elm02 = let
 end
 ```
 
+```{julia}
+#| include: false
+elm02 = let
+  form = @formula(
+    1000 / rt ~ 1 + isword * wrdlen + (1 | item) + (1 | subj)
+  )
+  dat = @subset(pruned, :spropacc > 0.8)
+  fit_or_restore("largescaledesigned_elm02.json", MixedModel, form, dat; contrasts, progress=false)
+end
+```
+
 ```{julia}
 effects(Dict(:isword => [false, true], :wrdlen => 4:2:12), elm02)
 ```
diff --git a/mrk17.qmd b/mrk17.qmd
index f540a18..b38a799 100644
--- a/mrk17.qmd
+++ b/mrk17.qmd
@@ -15,7 +15,7 @@ using MixedModels
 using MixedModelsMakie
 using StatsBase
 
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 using Statistics: mean, std
 
 const progress = isinteractive()
@@ -138,7 +138,10 @@ contrasts =
           :lQ =>EffectsCoding(; levels=["deg", "clr"]),
           :lT =>EffectsCoding(; levels=["NW", "WD"])
           );
+```
 
+```{julia}
+#| eval: false
 m_cpx = let
     form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
                                (1+F+P+Q+lQ+lT | subj) +
@@ -149,6 +152,16 @@ end
 VarCorr(m_cpx)
 ```
 
+```{julia}
+#| include: false
+m_cpx = let
+    form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
+                               (1+F+P+Q+lQ+lT | subj) +
+                               (1  +P+Q+lQ+lT | item);
+    fit_or_restore("mrk17_m_cpx.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 issingular(m_cpx)
 ```
@@ -183,6 +196,7 @@ VP is zero for fourth diagonal entry; not supported by data.
 We take out correlation parameters.
 
 ```{julia}
+#| eval: false
 m_zcp = let
     form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
                        zerocorr(1+F+P+Q+lQ+lT | subj) +
@@ -195,6 +209,16 @@ issingular(m_zcp)
 MixedModels.PCA(m_zcp)
 ```
 
+```{julia}
+#| include: false
+m_zcp = let
+    form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
+                       zerocorr(1+F+P+Q+lQ+lT | subj) +
+                       zerocorr(1  +P+Q+lQ+lT | item);
+    fit_or_restore("mrk17_m_zcp.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 MixedModels.likelihoodratiotest(m_zcp, m_cpx)
 ```
@@ -202,6 +226,7 @@ MixedModels.likelihoodratiotest(m_zcp, m_cpx)
 Looks ok. It might be a good idea to prune the LMM by removing small VCs.
 
 ```{julia}
+#| eval: false
 m_zcp2 = let
     form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
                        zerocorr(1  +P+Q+lQ+lT | subj) +
@@ -211,6 +236,16 @@ end
 VarCorr(m_zcp2)
 ```
 
+```{julia}
+#| include: false
+m_zcp2 = let
+    form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
+                       zerocorr(1  +P+Q+lQ+lT | subj) +
+                       zerocorr(1  +P+Q   +lT | item);
+    fit_or_restore("mrk17_m_zcp2.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 MixedModels.likelihoodratiotest(m_zcp2, m_zcp, m_cpx)
 ```
@@ -218,6 +253,7 @@ MixedModels.likelihoodratiotest(m_zcp2, m_zcp, m_cpx)
 We can perhaps remove some more.
 
 ```{julia}
+#| eval: false
 m_zcp3 = let
     form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
                        zerocorr(1    +Q   +lT | subj) + (1 | item);
@@ -226,6 +262,15 @@ end
 VarCorr(m_zcp3)
 ```
 
+```{julia}
+#| include: false
+m_zcp3 = let
+    form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
+                       zerocorr(1    +Q   +lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_zcp3.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 MixedModels.likelihoodratiotest(m_zcp3, m_zcp2, m_zcp, m_cpx)
 ```
@@ -233,6 +278,7 @@ MixedModels.likelihoodratiotest(m_zcp3, m_zcp2, m_zcp, m_cpx)
 And another iteration.
 
 ```{julia}
+#| eval: false
 m_zcp4 = let
     form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
                        zerocorr(1         +lT | subj) + (1 | item);
@@ -241,6 +287,15 @@ end
 VarCorr(m_zcp4)
 ```
 
+```{julia}
+#| include: false
+m_zcp4 = let
+    form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
+                       zerocorr(1         +lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_zcp4.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 MixedModels.likelihoodratiotest(m_zcp4, m_zcp3, m_zcp2, m_zcp, m_cpx)
 ```
@@ -248,6 +303,7 @@ MixedModels.likelihoodratiotest(m_zcp4, m_zcp3, m_zcp2, m_zcp, m_cpx)
 Too much removed. Stay with `m_zcp3`, but extend with CPs.
 
 ```{julia}
+#| eval: false
 m_prm = let
     form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
                                (1+    Q+lT | subj) + (1 | item);
@@ -255,8 +311,20 @@ m_prm = let
 end
 VarCorr(m_prm)
 ```
+
+```{julia}
+#| include: false
+m_prm = let
+    form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
+                               (1+    Q+lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_prm.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 ### post-hoc LMM
+
 ```{julia}
+#| eval: false
 m_prm = let
     form = @formula (1000/rt) ~ 1+F*P*Q*lQ*lT +
                                (1+    Q+lT | subj) + (1 | item);
@@ -300,6 +368,7 @@ The LMM is not nested in the previous sequence.
 ## Crossed fixed effects
 
 ```{julia}
+#| eval: false
 m_mrk17_crossed =let
    form = @formula (1000/rt) ~ 1 + F*P*Q*lQ*lT +
         (1+Q | subj) + zerocorr(0+lT | subj) + zerocorr(1 + P | item) ;
@@ -309,6 +378,15 @@ end
 VarCorr(m_prm)
 ```
 
+```{julia}
+#| include: false
+m_mrk17_crossed =let
+   form = @formula (1000/rt) ~ 1 + F*P*Q*lQ*lT +
+        (1+Q | subj) + zerocorr(0+lT | subj) + zerocorr(1 + P | item) ;
+    fit_or_restore("mrk17_m_mrk17_crossed.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 ```{julia}
 show(m_mrk17_crossed)
 ```
@@ -318,6 +396,7 @@ Finally, a look at the fixed effects. The four-factor interaction reported in Ma
 ## Nested fixed effects
 
 ```{julia}
+#| eval: false
 m_mrk17_nested =let
    form = @formula (1000/rt) ~ 1 + Q/(F/P) +
         (1+Q | subj) + zerocorr(0+lT | subj) + zerocorr(1 + P | item) ;
@@ -325,6 +404,15 @@ m_mrk17_nested =let
 end
 ```
 
+```{julia}
+#| include: false
+m_mrk17_nested =let
+   form = @formula (1000/rt) ~ 1 + Q/(F/P) +
+        (1+Q | subj) + zerocorr(0+lT | subj) + zerocorr(1 + P | item) ;
+    fit_or_restore("mrk17_m_mrk17_nested.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 # Questions from the Discussion forum
 
 ## Nesting within products of factors
@@ -332,6 +420,7 @@ end
 Include parenthesis
 
 ```{julia}
+#| eval: false
 m_mrk17_nested =let
    form = @formula (1000/rt) ~ 1 + Q/(F/P) +
         (1+Q | subj) + zerocorr(0+lT | subj) + zerocorr(1 + P | item) ;
@@ -342,6 +431,7 @@ end
 ## Selection in fixed effects
 
 ```{julia}
+#| eval: false
 using RegressionFormulae
 # m_prm_5 is equivalent to m_prm
 m_prm_5 = let
@@ -389,6 +479,55 @@ gof_summary = let
 end
 ```
 
+```{julia}
+#| include: false
+using RegressionFormulae
+# m_prm_5 is equivalent to m_prm
+m_prm_5 = let
+    form = @formula (1000/rt) ~ 1+(F+P+Q+lQ+lT)^5 + (1+Q+lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_prm_5.json", MixedModel, form, dat; contrasts, progress)
+end
+
+m_prm_4 = let
+    form = @formula (1000/rt) ~ 1+(F+P+Q+lQ+lT)^4 + (1+Q+lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_prm_4.json", MixedModel, form, dat; contrasts, progress)
+end
+
+m_prm_3 = let
+    form = @formula (1000/rt) ~ 1+(F+P+Q+lQ+lT)^3 + (1+Q+lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_prm_3.json", MixedModel, form, dat; contrasts, progress)
+end
+
+m_prm_2 = let
+    form = @formula (1000/rt) ~ 1+(F+P+Q+lQ+lT)^2 + (1+Q+lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_prm_2.json", MixedModel, form, dat; contrasts, progress)
+end
+
+m_prm_1 = let
+    form = @formula (1000/rt) ~ 1+ F+P+Q+lQ+lT +   (1+Q+lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_prm_1.json", MixedModel, form, dat; contrasts, progress)
+end
+
+# Compare the fits
+gof_summary = let
+  nms = [:m_prm_1, :m_prm_2, :m_prm_3, :m_prm_4, :m_prm_5]
+  mods = eval.(nms)
+  lrt = MixedModels.likelihoodratiotest(m_prm_1, m_prm_2, m_prm_3, m_prm_4, m_prm_5)
+  DataFrame(;
+    name = nms,
+    dof=dof.(mods),
+    deviance=deviance.(mods),
+    AIC=aic.(mods),
+    AICc=aicc.(mods),
+    BIC=bic.(mods),
+    χ²=vcat(:., round.(Int, diff(collect(lrt.lrt.deviance)))),
+    χ²_dof=vcat(:., diff(collect(lrt.lrt.dof))),
+    # StatsBase.PValue includes some pretty-printing methods
+    pvalue=vcat(:., StatsBase.PValue.(collect(lrt.lrt.pval[2:end])))
+  )
+end
+```
+
 Depending on the level of precision of your hypothesis. You could stay with main effect (BIC), include 2-factor interactions (AIC; also called _simple_ interactions) or include 3-factor interactions [χ² < 2*(χ²-dof); also called _2-way_ interactions].
 
 ## Posthoc LMM
@@ -396,6 +535,7 @@ Depending on the level of precision of your hypothesis. You could stay with main
 We are using only three factors for the illustruation.
 
 ```{julia}
+#| eval: false
 m_prm3 = let
     form = @formula (1000/rt) ~ 1 + lT*lQ*Q +
                                (1+    Q+lT | subj) + (1 | item);
@@ -403,9 +543,19 @@ m_prm3 = let
 end
 ```
 
+```{julia}
+#| include: false
+m_prm3 = let
+    form = @formula (1000/rt) ~ 1 + lT*lQ*Q +
+                               (1+    Q+lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_prm3.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 The `lT & lQ & Q` interactions is significant. Let's follow it up with a post-hoc LMM, that is looking at the `lQ & Q` interaction in the two levels of whether the last word was a target or not.
 
 ```{julia}
+#| eval: false
 m_prm3_posthoc = let
     form = @formula (1000/rt) ~ 1 + lT/(lQ*Q) +
                                (1+    Q+lT | subj) + (1 | item);
@@ -413,6 +563,15 @@ m_prm3_posthoc = let
 end
 ```
 
+```{julia}
+#| include: false
+m_prm3_posthoc = let
+    form = @formula (1000/rt) ~ 1 + lT/(lQ*Q) +
+                               (1+    Q+lT | subj) + (1 | item);
+    fit_or_restore("mrk17_m_prm3_posthoc.json", MixedModel, form, dat; contrasts, progress)
+end
+```
+
 The source of the interaction are trials where the last trial was a word target; there is no evidence for the interaction when the last trial was a nonword target.
 
 The original and post-hoc LMM have the same goodness of fit.
diff --git a/partial_within.qmd b/partial_within.qmd
index ccb49d2..20113a3 100644
--- a/partial_within.qmd
+++ b/partial_within.qmd
@@ -18,6 +18,7 @@ using MixedModelsMakie
 using MixedModelsSim
 using ProgressMeter
 using Random
+using SMLP2026: fit_or_restore
 
 const progress=false
 ```
@@ -39,6 +40,7 @@ unique!(select(design, :item, :frequency))
 ```
 
 ```{julia}
+#| eval: false
 #| code-fold: true
 m0 = let contrasts, form
     contrasts = Dict(:frequency => HelmertCoding(base="high"))
@@ -48,6 +50,16 @@ m0 = let contrasts, form
 end
 ```
 
+```{julia}
+#| include: false
+m0 = let contrasts, form
+    contrasts = Dict(:frequency => HelmertCoding(base="high"))
+    form = @formula(dv ~ 1 + frequency +
+                    (1 + frequency | subj))
+    fit_or_restore("partial_within_m0.json", MixedModel, form, design; contrasts, progress)
+end
+```
+
 ```{julia}
 #| code-fold: true
 corrmat = [ 1    0.1 -0.2
@@ -103,8 +115,8 @@ sort!(unique!(select(design_partial, :subj, :frequency)), :subj)
 ```
 
 ```{julia}
+#| eval: false
 #| code-fold: true
-
 m1 = let contrasts, form
     contrasts = Dict(:frequency => HelmertCoding(base="high"))
     form = @formula(dv ~ 1 + frequency +
@@ -113,6 +125,16 @@ m1 = let contrasts, form
 end
 ```
 
+```{julia}
+#| include: false
+m1 = let contrasts, form
+    contrasts = Dict(:frequency => HelmertCoding(base="high"))
+    form = @formula(dv ~ 1 + frequency +
+                    (1 + frequency | subj))
+    fit_or_restore("partial_within_m1.json", MixedModel, form, design_partial; contrasts, progress)
+end
+```
+
 ```{julia}
 #| code-fold: true
 shrinkageplot(m1)
diff --git a/power_simulation.qmd b/power_simulation.qmd
index 614452d..a45b4f0 100644
--- a/power_simulation.qmd
+++ b/power_simulation.qmd
@@ -20,7 +20,7 @@ using MixedModels
 using MixedModelsMakie
 using MixedModelsSim
 using StatsBase
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 using Random
 
 const progress=false
@@ -31,6 +31,7 @@ const progress=false
 Let us consider the `kb07` dataset.
 
 ```{julia}
+#| eval: false
 kb07 = dataset(:kb07)
 contrasts = Dict(:spkr => EffectsCoding(),
                  :prec => EffectsCoding(),
@@ -42,6 +43,19 @@ fm1 = fit(MixedModel,
           kb07; contrasts, progress)
 ```
 
+```{julia}
+#| include: false
+kb07 = dataset(:kb07)
+contrasts = Dict(:spkr => EffectsCoding(),
+                 :prec => EffectsCoding(),
+                 :load => EffectsCoding())
+fm1 = fit_or_restore("power_simulation_fm1.json", MixedModel,
+                     @formula(rt_trunc ~ 1 * spkr * prec * load +
+                                        (1 | subj) +
+                                        (1 | item)),
+                     kb07; contrasts, progress)
+```
+
 We can perform a *parametric bootstrap* on the model to get estimates of our uncertainty.
 In the parametric bootstrap, we use the *parameters* we estimated to simulate new data.
 If we repeat this process many times, we are able to "pick ourselves up by our bootstraps"
diff --git a/profiling.qmd b/profiling.qmd
index d0bd0ca..cbaddc3 100644
--- a/profiling.qmd
+++ b/profiling.qmd
@@ -34,7 +34,7 @@ using CairoMakie
 using MixedModels
 using MixedModelsMakie
 using Random
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 
 const progress=false
 ```
@@ -55,6 +55,7 @@ The response is defined as `1000 / rt_raw` where `rt_raw` is measured in millise
 Thus the response being modeled is the speed measured in responses per second.
 
 ```{julia}
+#| eval: false
 pr01 = let f = @formula 1000 / rt_raw ~
     1 + load + spkr + prec + (1 + prec | item) + (1 | subj)
   profile(fit(MixedModel, f, kb07; contrasts, progress))
@@ -62,6 +63,15 @@ end
 println(pr01.m) # model is a property of the profile object
 ```
 
+```{julia}
+#| include: false
+pr01 = let f = @formula 1000 / rt_raw ~
+    1 + load + spkr + prec + (1 + prec | item) + (1 | subj)
+  profile(fit_or_restore("profiling_pr01.json", MixedModel, f, kb07; contrasts, progress))
+end
+println(pr01.m) # model is a property of the profile object
+```
+
 Evaluation of `pr01` is similar to other model fits in these notes except that the call to `fit` is wrapped in a call to `profile`.
 Because the object returned from `profile` includes the original model fit as its `m` property, it is not necessary to save the original model fit separately.
 
diff --git a/shrinkageplot.qmd b/shrinkageplot.qmd
index 9cb32e2..6edfe6e 100644
--- a/shrinkageplot.qmd
+++ b/shrinkageplot.qmd
@@ -20,6 +20,7 @@ using MixedModels
 using MixedModelsMakie
 using Random
 using ProgressMeter
+using SMLP2026: fit_or_restore
 
 const progress = false
 ```
@@ -31,6 +32,7 @@ kb07 = MixedModels.dataset(:kb07)
 ```
 
 ```{julia}
+#| eval: false
 contrasts = Dict(
   :spkr => HelmertCoding(),
   :prec => HelmertCoding(),
@@ -48,6 +50,25 @@ m1 = let
 end
 ```
 
+```{julia}
+#| include: false
+contrasts = Dict(
+  :spkr => HelmertCoding(),
+  :prec => HelmertCoding(),
+  :load => HelmertCoding(),
+)
+m1 = let
+  form = @formula(
+    rt_trunc ~
+      1 +
+      spkr * prec * load +
+      (1 + spkr + prec + load | subj) +
+      (1 + spkr + prec + load | item)
+  )
+  fit_or_restore("shrinkageplot_m1.json", MixedModel, form, kb07; contrasts, progress)
+end
+```
+
 ```{julia}
 VarCorr(m1)
 ```
@@ -125,6 +146,7 @@ X1 * X1'
   - For the `item` both the intercept and the `spkr` appear to be contributing
 
 ```{julia}
+#| eval: false
 m2 = let
   form = @formula(
     rt_trunc ~
@@ -134,6 +156,17 @@ m2 = let
 end
 ```
 
+```{julia}
+#| include: false
+m2 = let
+  form = @formula(
+    rt_trunc ~
+      1 + prec * spkr * load + (1 | subj) + (1 + prec | item)
+  )
+  fit_or_restore("shrinkageplot_m2.json", MixedModel, form, kb07; contrasts, progress)
+end
+```
+
 ```{julia}
 VarCorr(m2)
 ```
@@ -146,6 +179,7 @@ shrinkageplot(m2)
 ```
 
 ```{julia}
+#| eval: false
 #| lst-label: m1def
 m3 = let
   form = @formula(
@@ -156,6 +190,17 @@ m3 = let
 end
 ```
 
+```{julia}
+#| include: false
+m3 = let
+  form = @formula(
+    rt_trunc ~
+      1 + prec + spkr + load + (1 | subj) + (1 + prec | item)
+  )
+  fit_or_restore("shrinkageplot_m3.json", MixedModel, form, kb07; contrasts, progress)
+end
+```
+
 ```{julia}
 VarCorr(m3)
 ```
@@ -185,6 +230,7 @@ ridgeplot(m3btstrp; show_intercept=false)
 ```
 
 ```{julia}
+#| eval: false
 m4 = let
   form = @formula(
     rt_trunc ~
@@ -194,6 +240,17 @@ m4 = let
 end
 ```
 
+```{julia}
+#| include: false
+m4 = let
+  form = @formula(
+    rt_trunc ~
+      1 + prec + spkr + load + (1 + prec | item) + (1 | subj)
+  )
+  fit_or_restore("shrinkageplot_m4.json", MixedModel, form, kb07; contrasts, progress)
+end
+```
+
 ```{julia}
 m4bstrp = parametricbootstrap(rng, 2000, m4);
 ```
diff --git a/sleepstudy.qmd b/sleepstudy.qmd
index 243de32..8a5b7ea 100644
--- a/sleepstudy.qmd
+++ b/sleepstudy.qmd
@@ -34,7 +34,7 @@ using MixedModelsMakie # diagnostic plots
 using Random           # random number generators
 using RCall            # call R from Julia
 using MixedModelsMakie: simplelinreg
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 ```
 
 The `sleepstudy` data are one of the datasets available with the `MixedModels` package.  It is re-exported by the `SMLP2026` package's `dataset` function.
@@ -114,6 +114,7 @@ There are some deviations from linearity within the panels but the deviations ar
 # Fitting an initial model
 
 ```{julia}
+#| eval: false
 #| lst-label: m1
 contrasts = Dict{Symbol,Any}(:subj => Grouping())
 m1 = let f = @formula(reaction ~ 1 + days + (1 + days | subj))
@@ -121,6 +122,14 @@ m1 = let f = @formula(reaction ~ 1 + days + (1 + days | subj))
 end
 ```
 
+```{julia}
+#| include: false
+contrasts = Dict{Symbol,Any}(:subj => Grouping())
+m1 = let f = @formula(reaction ~ 1 + days + (1 + days | subj))
+  fit_or_restore("sleepstudy_m1.json", MixedModel, f, sleepstudy; contrasts)
+end
+```
+
 This model includes fixed effects for the intercept, representing the typical reaction time at the beginning of the experiment with zero days of sleep deprivation, and the slope w.r.t. days of sleep deprivation.
 The parameter estimates are about 250 ms. typical reaction time without deprivation and a typical increase of 10.5 ms. per day of sleep deprivation.
 
@@ -158,12 +167,20 @@ caterpillar(m1; vline_at_zero=true)
 The `zerocorr` function applied to a random-effects term creates uncorrelated vector-valued per-subject random effects.
 
 ```{julia}
+#| eval: false
 #| lst-label: m2
 m2 = let f = @formula reaction ~ 1 + days + zerocorr(1 + days | subj)
   fit(MixedModel, f, sleepstudy; contrasts)
 end
 ```
 
+```{julia}
+#| include: false
+m2 = let f = @formula reaction ~ 1 + days + zerocorr(1 + days | subj)
+  fit_or_restore("sleepstudy_m2.json", MixedModel, f, sleepstudy; contrasts)
+end
+```
+
 Again, the default display doesn't show that there is no correlation parameter to be estimated in this model, but the `VarCorr` display does.
 
 ```{julia}
diff --git a/sleepstudy_speed.qmd b/sleepstudy_speed.qmd
index 420322d..f20b85e 100644
--- a/sleepstudy_speed.qmd
+++ b/sleepstudy_speed.qmd
@@ -39,7 +39,7 @@ using DataFrameMacros
 using DataFrames
 using MixedModels
 using MixedModelsMakie   # plots specific to mixed-effects models using Makie
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 ```
 
 # Preprocessing
@@ -108,12 +108,21 @@ end
 # Basic LMM
 
 ```{julia}
+#| eval: false
 m1 = let
   form = @formula speed ~ 1 + day + (1 + day | Subj)
   fit(MixedModel, form, df)
 end
 ```
 
+```{julia}
+#| include: false
+m1 = let
+  form = @formula speed ~ 1 + day + (1 + day | Subj)
+  fit_or_restore("sleepstudy_speed_m1.json", MixedModel, form, df)
+end
+```
+
 This model includes fixed effects for the intercept which estimates the average speed on the baseline day of the experiment prior to sleep deprivation, and the slowing per day of sleep deprivation. In this case about -0.11/second.
 
 The random effects represent shifts from the typical behavior for each subject.The shift in the intercept has a standard deviation of about 0.42/s.
@@ -125,12 +134,21 @@ The within-subject correlation of the random effects for intercept and slope is
 The `zerocorr` function applied to a random-effects term estimates one parameter less than LMM `m1`-- the CP is now fixed to zero.
 
 ```{julia}
+#| eval: false
 m2 = let
   form = @formula speed ~ 1 + day + zerocorr(1 + day | Subj)
   fit(MixedModel, form, df)
 end
 ```
 
+```{julia}
+#| include: false
+m2 = let
+  form = @formula speed ~ 1 + day + zerocorr(1 + day | Subj)
+  fit_or_restore("sleepstudy_speed_m2.json", MixedModel, form, df)
+end
+```
+
 LMM `m2` has a slghtly  lower log-likelihood than LMM `m1` but also one fewer parameters.
 A likelihood-ratio test is used to compare these nested models.
 
diff --git a/src/SMLP2026.jl b/src/SMLP2026.jl
index 7f28e3c..8ef1a84 100644
--- a/src/SMLP2026.jl
+++ b/src/SMLP2026.jl
@@ -68,7 +68,7 @@ function fit_or_restore(fname, ::Type{<:MixedModel}, args...; kwargs...)
     return fit_or_restore(fname, args...; kwargs...)
 end
 
-function fit_or_restore(fname, args...; contrasts=Dict{Symbol}(), kwargs...)
+function fit_or_restore(fname, args...; contrasts=Dict{Symbol,Any}(), kwargs...)
     model = MixedModel(args...; contrasts)
     return fit_or_restore!(model, fname; kwargs...)
 end
diff --git a/transformation.qmd b/transformation.qmd
index 68d8809..5e183b2 100644
--- a/transformation.qmd
+++ b/transformation.qmd
@@ -47,22 +47,39 @@ using DataFrames
 using Effects
 using MixedModels
 using StandardizedPredictors
-using SMLP2026: dataset
+using SMLP2026: dataset, fit_or_restore
 ```
 
 ```{julia}
-slp = fit(MixedModel, 
+#| eval: false
+slp = fit(MixedModel,
           @formula(reaction ~ 1 + days + (1 + days |subj)),
           dataset(:sleepstudy))
 ```
 
 ```{julia}
-days_centered = fit(MixedModel, 
+#| include: false
+slp = fit_or_restore("transformation_slp.json", MixedModel,
+                     @formula(reaction ~ 1 + days + (1 + days |subj)),
+                     dataset(:sleepstudy))
+```
+
+```{julia}
+#| eval: false
+days_centered = fit(MixedModel,
                     @formula(reaction ~ 1 + days + (1 + days |subj)),
                     dataset(:sleepstudy);
                     contrasts=Dict(:days => Center()))
 ```
 
+```{julia}
+#| include: false
+days_centered = fit_or_restore("transformation_days_centered.json", MixedModel,
+                               @formula(reaction ~ 1 + days + (1 + days |subj)),
+                               dataset(:sleepstudy);
+                               contrasts=Dict(:days => Center()))
+```
+
 If we look at the log-likelihood, AIC, BIC, etc. of these two models, we see that they are the same:
 
 ```{julia}
@@ -155,9 +172,17 @@ While useful at times, speed has a natural interpretation and so we instead use
 Because `reaction` is stored in milliseconds, we use `1000 / reaction` instead of `1 / reaction` so that our speed units are responses per second.
 
 ```{julia}
+#| eval: false
 model_bc = fit(MixedModel,
                @formula(1000 / reaction ~ 1 + days + (1 + days | subj)),
-                dataset(:sleepstudy))
+               dataset(:sleepstudy))
+```
+
+```{julia}
+#| include: false
+model_bc = fit_or_restore("transformation_model_bc.json", MixedModel,
+                          @formula(1000 / reaction ~ 1 + days + (1 + days | subj)),
+                          dataset(:sleepstudy))
 ```
 
 For our original model on the untransformed scale, the intercept was approximately 250, which means that the average response time was about 250 milliseconds.

From 89a519274b91c21d99c6f35173dfad8f963a827b Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 30 Mar 2026 18:42:02 -0500
Subject: [PATCH 7/9] call to lmm

---
 contrasts_fggk21.qmd                          | 154 +++++++++++++++++-
 fits/contrasts_fggk21_m_Child_HeC.json.zip    |   3 +
 fits/contrasts_fggk21_m_Child_HeC0.json.zip   |   3 +
 fits/contrasts_fggk21_m_Child_HyC.json.zip    |   3 +
 fits/contrasts_fggk21_m_Child_HyC0.json.zip   |   3 +
 fits/contrasts_fggk21_m_Child_PCA.json.zip    |   3 +
 fits/contrasts_fggk21_m_Child_PCA0.json.zip   |   3 +
 fits/contrasts_fggk21_m_Child_SDC.json.zip    |   3 +
 fits/contrasts_fggk21_m_Child_SDC0.json.zip   |   3 +
 fits/contrasts_fggk21_m_School_HeC0.json.zip  |   3 +
 ...contrasts_fggk21_m_School_Helmert.json.zip |   3 +
 fits/contrasts_fggk21_m_School_HyC0.json.zip  |   3 +
 fits/contrasts_fggk21_m_School_Hypo.json.zip  |   3 +
 fits/contrasts_fggk21_m_School_PCA.json.zip   |   3 +
 fits/contrasts_fggk21_m_School_PCA0.json.zip  |   3 +
 fits/contrasts_fggk21_m_School_SDC0.json.zip  |   3 +
 ...contrasts_fggk21_m_School_SeqDiff.json.zip |   3 +
 .../contrasts_fggk21_m_cpx_0_SeqDiff.json.zip |   3 +
 fits/contrasts_fggk21_m_cpx_1_PC.json.zip     |   3 +
 fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip   |   3 +
 .../contrasts_fggk21_m_cpx_1_SeqDiff.json.zip |   3 +
 fits/contrasts_fggk21_m_ovi_Helmert.json.zip  |   3 +
 fits/contrasts_fggk21_m_ovi_Hypo.json.zip     |   3 +
 fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip  |   3 +
 .../contrasts_fggk21_m_ovi_SeqDiff_1.json.zip |   3 +
 .../contrasts_fggk21_m_ovi_SeqDiff_2.json.zip |   3 +
 .../contrasts_fggk21_m_ovi_SeqDiff_3.json.zip |   3 +
 ...contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip |   3 +
 fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip   |   3 +
 fits/contrasts_fggk21_m_zcp_SeqD.json.zip     |   3 +
 fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip   |   3 +
 fits/singularity_m01.json.zip                 |   3 +
 singularity.qmd                               |  13 ++
 src/SMLP2026.jl                               |  10 +-
 34 files changed, 263 insertions(+), 7 deletions(-)
 create mode 100644 fits/contrasts_fggk21_m_Child_HeC.json.zip
 create mode 100644 fits/contrasts_fggk21_m_Child_HeC0.json.zip
 create mode 100644 fits/contrasts_fggk21_m_Child_HyC.json.zip
 create mode 100644 fits/contrasts_fggk21_m_Child_HyC0.json.zip
 create mode 100644 fits/contrasts_fggk21_m_Child_PCA.json.zip
 create mode 100644 fits/contrasts_fggk21_m_Child_PCA0.json.zip
 create mode 100644 fits/contrasts_fggk21_m_Child_SDC.json.zip
 create mode 100644 fits/contrasts_fggk21_m_Child_SDC0.json.zip
 create mode 100644 fits/contrasts_fggk21_m_School_HeC0.json.zip
 create mode 100644 fits/contrasts_fggk21_m_School_Helmert.json.zip
 create mode 100644 fits/contrasts_fggk21_m_School_HyC0.json.zip
 create mode 100644 fits/contrasts_fggk21_m_School_Hypo.json.zip
 create mode 100644 fits/contrasts_fggk21_m_School_PCA.json.zip
 create mode 100644 fits/contrasts_fggk21_m_School_PCA0.json.zip
 create mode 100644 fits/contrasts_fggk21_m_School_SDC0.json.zip
 create mode 100644 fits/contrasts_fggk21_m_School_SeqDiff.json.zip
 create mode 100644 fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip
 create mode 100644 fits/contrasts_fggk21_m_cpx_1_PC.json.zip
 create mode 100644 fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip
 create mode 100644 fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip
 create mode 100644 fits/contrasts_fggk21_m_ovi_Helmert.json.zip
 create mode 100644 fits/contrasts_fggk21_m_ovi_Hypo.json.zip
 create mode 100644 fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip
 create mode 100644 fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip
 create mode 100644 fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip
 create mode 100644 fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip
 create mode 100644 fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip
 create mode 100644 fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip
 create mode 100644 fits/contrasts_fggk21_m_zcp_SeqD.json.zip
 create mode 100644 fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip
 create mode 100644 fits/singularity_m01.json.zip

diff --git a/contrasts_fggk21.qmd b/contrasts_fggk21.qmd
index 1b997ec..485345b 100644
--- a/contrasts_fggk21.qmd
+++ b/contrasts_fggk21.qmd
@@ -56,6 +56,11 @@ using SMLP2026: dataset
 progress = isinteractive()
 ```
 
+```{julia}
+#| echo: false
+using SMLP2026: fit_or_restore
+```
+
 ## Readme for `dataset("fggk21")`
 
 Number of scores: 525126
@@ -160,19 +165,35 @@ contr1 = Dict(
 ```
 
 ```{julia}
+#| eval: false
 m_ovi_SeqDiff_1 = lmm(@formula(zScore ~ 1 + Test + (1 | Child)),
                       dat; contrasts=contr1, progress)
 ```
 
+```{julia}
+#| include: false
+m_ovi_SeqDiff_1 = fit_or_restore("contrasts_fggk21_m_ovi_SeqDiff_1.json",
+                                  @formula(zScore ~ 1 + Test + (1 | Child)),
+                                  dat; contrasts=contr1, progress)
+```
+
 In this case, any differences between tests identified by the contrasts would be  spurious because each test was standardized (i.e., _M_=0, $SD$=1). The differences could also be due to an imbalance in the number of boys and girls or in the number of missing observations for each test.
 
 The primary interest in this study related to interactions of the test contrasts with and `age` and `Sex`. We start with age (linear) and its interaction with the four test contrasts.
 
 ```{julia}
+#| eval: false
 m_ovi_SeqDiff_2 = lmm(@formula(zScore ~ 1 + Test * a1 + (1 | Child)),
                       dat; contrasts=contr1, progress)
 ```
 
+```{julia}
+#| include: false
+m_ovi_SeqDiff_2 = fit_or_restore("contrasts_fggk21_m_ovi_SeqDiff_2.json",
+                                  @formula(zScore ~ 1 + Test * a1 + (1 | Child)),
+                                  dat; contrasts=contr1, progress)
+```
+
 The difference between older and younger children is larger for `Star_r` than for `Run` (0.2473).  `S20_r` did not differ significantly from `Star_r` (-0.0377) and `SLJ` (-0.0113) The largest difference in developmental gain was between `BPT` and `SLJ` (0.3355).
 
 **Please note that standard errors of this LMM are anti-conservative because the LMM is missing a lot of information in the RES (e..g., contrast-related VCs snd CPs for `Child`, `School`, and `Cohort`.**
@@ -180,10 +201,18 @@ The difference between older and younger children is larger for `Star_r` than fo
 Next we add the main effect of `Sex` and its interaction with the four test contrasts.
 
 ```{julia}
+#| eval: false
 m_ovi_SeqDiff_3 = lmm(@formula(zScore ~ 1 + Test * (a1 + Sex) + (1 | Child)),
                       dat; contrasts=contr1, progress)
 ```
 
+```{julia}
+#| include: false
+m_ovi_SeqDiff_3 = fit_or_restore("contrasts_fggk21_m_ovi_SeqDiff_3.json",
+                                  @formula(zScore ~ 1 + Test * (a1 + Sex) + (1 | Child)),
+                                  dat; contrasts=contr1, progress)
+```
+
 The significant interactions with `Sex` reflect mostly differences related to muscle power, where the physiological constitution gives boys an advantage. The sex difference is smaller when coordination and cognition play a role -- as in the `Star_r` test. (Caveat: SEs are estimated with an underspecified RES.)
 
 The final step in this first series is to add the interactions between the three covariates. A significant interaction between any of the four `Test` contrasts and age (linear) x `Sex` was hypothesized to reflect a prepubertal signal (i.e., hormones start to rise in girls' ninth year of life). However, this hypothesis is linked to a specific shape of the interaction: Girls would need to gain more than boys in tests of muscular power.
@@ -192,9 +221,18 @@ The final step in this first series is to add the interactions between the three
 # we'll be re-using this formula with other contrast specifications,
 # so we assign it to its own variable
 f_ovi = @formula(zScore ~ 1 + Test * a1 * Sex + (1 | Child))
+```
+
+```{julia}
+#| eval: false
 m_ovi_SeqDiff = lmm(f_ovi, dat; contrasts=contr1, progress)
 ```
 
+```{julia}
+#| include: false
+m_ovi_SeqDiff = fit_or_restore("contrasts_fggk21_m_ovi_SeqDiff.json", f_ovi, dat; contrasts=contr1, progress)
+```
+
 The results are very clear: Despite an abundance of statistical power there is no evidence for the differences between boys and girls in how much they gain in the ninth year of life in these five tests. The authors argue that, in this case, absence of evidence looks very much like evidence of absence of a hypothesized interaction.
 
 In the next two sections we use different contrasts. Does this have a bearing on this result?  We still ignore for now that we are looking at anti-conservative test statistics.
@@ -230,9 +268,15 @@ contr2 = Dict(
 ```
 
 ```{julia}
+#| eval: false
 m_ovi_Helmert = lmm(f_ovi, dat; contrasts=contr2, progress)
 ```
 
+```{julia}
+#| include: false
+m_ovi_Helmert = fit_or_restore("contrasts_fggk21_m_ovi_Helmert.json", f_ovi, dat; contrasts=contr2, progress)
+```
+
 We forego a detailed discussion of the effects, but note that again none of the interactions between `age x Sex` with the four test contrasts was significant.
 
 The default labeling of Helmert contrasts may lead to confusions with other contrasts. Therefore, we could provide our own labels:
@@ -267,9 +311,15 @@ contr3 = Dict(
 ```
 
 ```{julia}
+#| eval: false
 m_ovi_Hypo = lmm(f_ovi, dat; contrasts=contr3, progress)
 ```
 
+```{julia}
+#| include: false
+m_ovi_Hypo = fit_or_restore("contrasts_fggk21_m_ovi_Hypo.json", f_ovi, dat; contrasts=contr3, progress)
+```
+
 With _HypothesisCoding_ we must generate our own labels for the contrasts. The default labeling of contrasts is usually not interpretable. Therefore, we provide our own.
 
 Anyway, none of the interactions between `age` x `Sex` with the four `Test` contrasts was significant for these contrasts.
@@ -291,24 +341,54 @@ contr1b = Dict(
 ```
 
 ```{julia}
+#| eval: false
 m_ovi_SeqDiff_v2 = lmm(f_ovi, dat; contrasts=contr1b, progress)
 ```
 
 ```{julia}
+#| include: false
+m_ovi_SeqDiff_v2 = fit_or_restore("contrasts_fggk21_m_ovi_SeqDiff_v2.json", f_ovi, dat; contrasts=contr1b, progress)
+```
+
+```{julia}
+#| eval: false
 m_zcp_SeqD = lmm(@formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Child)),
                  dat; contrasts=contr1b, progress)
 ```
 
 ```{julia}
+#| include: false
+m_zcp_SeqD = fit_or_restore("contrasts_fggk21_m_zcp_SeqD.json",
+                             @formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Child)),
+                             dat; contrasts=contr1b, progress)
+```
+
+```{julia}
+#| eval: false
 m_zcp_SeqD_2 = lmm(@formula(zScore ~ 1 + Test * a1 * Sex + (0 + Test | Child)),
                    dat; contrasts=contr1b, progress)
 ```
 
 ```{julia}
+#| include: false
+m_zcp_SeqD_2 = fit_or_restore("contrasts_fggk21_m_zcp_SeqD_2.json",
+                               @formula(zScore ~ 1 + Test * a1 * Sex + (0 + Test | Child)),
+                               dat; contrasts=contr1b, progress)
+```
+
+```{julia}
+#| eval: false
 m_cpx_0_SeqDiff = lmm(@formula(zScore ~ 1 + Test * a1 * Sex + (0 + Test | Child)),
                       dat; contrasts=contr1b, progress)
 ```
 
+```{julia}
+#| include: false
+m_cpx_0_SeqDiff = fit_or_restore("contrasts_fggk21_m_cpx_0_SeqDiff.json",
+                                  @formula(zScore ~ 1 + Test * a1 * Sex + (0 + Test | Child)),
+                                  dat; contrasts=contr1b, progress)
+```
+
 ```{julia}
 VarCorr(m_cpx_0_SeqDiff)
 ```
@@ -319,9 +399,18 @@ m_cpx_0_SeqDiff.PCA
 
 ```{julia}
 f_cpx_1 = @formula(zScore ~ 1 + Test * a1 * Sex + (1 + Test | Child))
+```
+
+```{julia}
+#| eval: false
 m_cpx_1_SeqDiff = lmm(f_cpx_1, dat; contrasts=contr1b, progress)
 ```
 
+```{julia}
+#| include: false
+m_cpx_1_SeqDiff = fit_or_restore("contrasts_fggk21_m_cpx_1_SeqDiff.json", f_cpx_1, dat; contrasts=contr1b, progress)
+```
+
 ```{julia}
 m_cpx_1_SeqDiff.PCA
 ```
@@ -354,9 +443,15 @@ contr4 = Dict(
 ```
 
 ```{julia}
+#| eval: false
 m_cpx_1_PC = lmm(f_cpx_1, dat; contrasts=contr4, progress)
 ```
 
+```{julia}
+#| include: false
+m_cpx_1_PC = fit_or_restore("contrasts_fggk21_m_cpx_1_PC.json", f_cpx_1, dat; contrasts=contr4, progress)
+```
+
 ```{julia}
 VarCorr(m_cpx_1_PC)
 ```
@@ -384,10 +479,16 @@ contr4b = Dict(
 ```
 
 ```{julia}
+#| eval: false
 # LN_NEWUOA does a poor job here
 m_cpx_1_PC_2 = lmm(f_cpx_1, dat; contrasts=contr4b, progress, optimizer=:LN_BOBYQA)
 ```
 
+```{julia}
+#| include: false
+m_cpx_1_PC_2 = fit_or_restore("contrasts_fggk21_m_cpx_1_PC_2.json", f_cpx_1, dat; contrasts=contr4b, progress, optimizer=:LN_BOBYQA)
+```
+
 ```{julia}
 VarCorr(m_cpx_1_PC_2)
 ```
@@ -397,11 +498,19 @@ m_cpx_1_PC_2.PCA
 ```
 
 ```{julia}
+#| eval: false
 # LN_BOBYQA does a poor job here
 m_zcp_1_PC_2 = lmm(@formula(zScore ~ 1 + Test*a1*Sex + zerocorr(1 + Test | Child)),
                    dat; contrasts=contr4b, progress, optimizer=:LN_NEWUOA)
 ```
 
+```{julia}
+#| include: false
+m_zcp_1_PC_2 = fit_or_restore("contrasts_fggk21_m_zcp_1_PC_2.json",
+                               @formula(zScore ~ 1 + Test*a1*Sex + zerocorr(1 + Test | Child)),
+                               dat; contrasts=contr4b, progress, optimizer=:LN_NEWUOA)
+```
+
 ```{julia}
 VarCorr(m_zcp_1_PC_2)
 ```
@@ -431,7 +540,7 @@ Trivially, the meaning of a contrast depends on its definition. Consequently, th
 
 <!-- These optima are really poorly defined! The choice of optimizer makes a big difference -->
 ```{julia}
-#| output: false
+#| eval: false
 f_Child = @formula(zScore ~ 1 + Test * a1 * Sex + (1 + Test | Child))
 m_Child_SDC = lmm(f_Child, dat; contrasts=contr1, progress, optimizer=:LN_BOBYQA)
 m_Child_HeC = lmm(f_Child, dat; contrasts=contr2, progress, optimizer=:LN_NELDERMEAD)
@@ -439,6 +548,15 @@ m_Child_HyC = lmm(f_Child, dat; contrasts=contr3, progress, optimizer=:LN_NELDER
 m_Child_PCA = lmm(f_Child, dat; contrasts=contr4, progress, optimizer=:LN_NELDERMEAD)
 ```
 
+```{julia}
+#| include: false
+f_Child = @formula(zScore ~ 1 + Test * a1 * Sex + (1 + Test | Child))
+m_Child_SDC = fit_or_restore("contrasts_fggk21_m_Child_SDC.json", f_Child, dat; contrasts=contr1, progress, optimizer=:LN_BOBYQA)
+m_Child_HeC = fit_or_restore("contrasts_fggk21_m_Child_HeC.json", f_Child, dat; contrasts=contr2, progress, optimizer=:LN_NELDERMEAD)
+m_Child_HyC = fit_or_restore("contrasts_fggk21_m_Child_HyC.json", f_Child, dat; contrasts=contr3, progress, optimizer=:LN_NELDERMEAD)
+m_Child_PCA = fit_or_restore("contrasts_fggk21_m_Child_PCA.json", f_Child, dat; contrasts=contr4, progress, optimizer=:LN_NELDERMEAD)
+```
+
 ```{julia}
 VarCorr(m_Child_SDC)
 ```
@@ -460,7 +578,7 @@ The CPs for the various contrasts are in line with expectations. For the SDC we
 Do these differences in CPs imply that we can move to zcpLMMs when we have orthogonal contrasts? We pursue this question with by refitting the four LMMs with zerocorr() and compare the goodness of fit.
 
 ```{julia}
-#| output: false
+#| eval: false
 f_Child0 = @formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Child))
 m_Child_SDC0 = lmm(f_Child0, dat; contrasts=contr1, progress, optimizer=:LN_NEWUOA)
 m_Child_HeC0 = lmm(f_Child0, dat; contrasts=contr2, progress, optimizer=:LN_NEWUOA)
@@ -468,6 +586,15 @@ m_Child_HyC0 = lmm(f_Child0, dat; contrasts=contr3, progress, optimizer=:LN_NEWU
 m_Child_PCA0 = lmm(f_Child0, dat; contrasts=contr4, progress, optimizer=:LN_NEWUOA)
 ```
 
+```{julia}
+#| include: false
+f_Child0 = @formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test | Child))
+m_Child_SDC0 = fit_or_restore("contrasts_fggk21_m_Child_SDC0.json", f_Child0, dat; contrasts=contr1, progress, optimizer=:LN_NEWUOA)
+m_Child_HeC0 = fit_or_restore("contrasts_fggk21_m_Child_HeC0.json", f_Child0, dat; contrasts=contr2, progress, optimizer=:LN_NEWUOA)
+m_Child_HyC0 = fit_or_restore("contrasts_fggk21_m_Child_HyC0.json", f_Child0, dat; contrasts=contr3, progress, optimizer=:LN_NEWUOA)
+m_Child_PCA0 = fit_or_restore("contrasts_fggk21_m_Child_PCA0.json", f_Child0, dat; contrasts=contr4, progress, optimizer=:LN_NEWUOA)
+```
+
 ```{julia}
 likelihoodratiotest(m_Child_SDC0, m_Child_SDC)
 ```
@@ -512,7 +639,7 @@ The effect of `age` (i.e., developmental gain) varies within `School`.
 Therefore, we also include its VCs and CPs in this model; the school-related VC for `Sex` was not significant.
 
 ```{julia}
-#| output: false
+#| eval: false
 f_School = @formula(zScore ~ 1 + Test * a1 * Sex + (1 + Test + a1 | School))
 m_School_SeqDiff = lmm(f_School, dat; contrasts=contr1, progress)
 m_School_Helmert = lmm(f_School, dat; contrasts=contr2, progress)
@@ -520,6 +647,15 @@ m_School_Hypo = lmm(f_School, dat; contrasts=contr3, progress)
 m_School_PCA = lmm(f_School, dat; contrasts=contr4, progress)
 ```
 
+```{julia}
+#| include: false
+f_School = @formula(zScore ~ 1 + Test * a1 * Sex + (1 + Test + a1 | School))
+m_School_SeqDiff = fit_or_restore("contrasts_fggk21_m_School_SeqDiff.json", f_School, dat; contrasts=contr1)
+m_School_Helmert = fit_or_restore("contrasts_fggk21_m_School_Helmert.json", f_School, dat; contrasts=contr2)
+m_School_Hypo = fit_or_restore("contrasts_fggk21_m_School_Hypo.json", f_School, dat; contrasts=contr3)
+m_School_PCA = fit_or_restore("contrasts_fggk21_m_School_PCA.json", f_School, dat; contrasts=contr4)
+```
+
 ```{julia}
 VarCorr(m_School_SeqDiff)
 ```
@@ -539,12 +675,24 @@ VarCorr(m_School_PCA)
 We compare again how much of the fit resides in the CPs.
 
 ```{julia}
+#| eval: false
 f_School0 = @formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test + a1 | School))
 m_School_SDC0 = lmm(f_School0, dat; contrasts=contr1, progress)
 m_School_HeC0 = lmm(f_School0, dat; contrasts=contr2, progress)
 m_School_HyC0 = lmm(f_School0, dat; contrasts=contr3, progress)
 m_School_PCA0 = lmm(f_School0, dat; contrasts=contr4, progress)
+```
 
+```{julia}
+#| include: false
+f_School0 = @formula(zScore ~ 1 + Test * a1 * Sex + zerocorr(1 + Test + a1 | School))
+m_School_SDC0 = fit_or_restore("contrasts_fggk21_m_School_SDC0.json", f_School0, dat; contrasts=contr1, progress)
+m_School_HeC0 = fit_or_restore("contrasts_fggk21_m_School_HeC0.json", f_School0, dat; contrasts=contr2, progress)
+m_School_HyC0 = fit_or_restore("contrasts_fggk21_m_School_HyC0.json", f_School0, dat; contrasts=contr3, progress)
+m_School_PCA0 = fit_or_restore("contrasts_fggk21_m_School_PCA0.json", f_School0, dat; contrasts=contr4, progress)
+```
+
+```{julia}
 zcpLMM2 = ["SDC0", "HeC0", "HyC0", "PCA0"]
 mods2 = [
   m_School_SDC0, m_School_HeC0, m_School_HyC0, m_School_PCA0
diff --git a/fits/contrasts_fggk21_m_Child_HeC.json.zip b/fits/contrasts_fggk21_m_Child_HeC.json.zip
new file mode 100644
index 0000000..4ae8d8e
--- /dev/null
+++ b/fits/contrasts_fggk21_m_Child_HeC.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c57f271a20cfaae2b55456ffa762b89378280a4f1c74993eef6df39addf42a08
+size 197261
diff --git a/fits/contrasts_fggk21_m_Child_HeC0.json.zip b/fits/contrasts_fggk21_m_Child_HeC0.json.zip
new file mode 100644
index 0000000..43c4377
--- /dev/null
+++ b/fits/contrasts_fggk21_m_Child_HeC0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ba371128e8820b0ee4a30678439b8623675301822ba6b5321c3fe3ade281ddbd
+size 81800
diff --git a/fits/contrasts_fggk21_m_Child_HyC.json.zip b/fits/contrasts_fggk21_m_Child_HyC.json.zip
new file mode 100644
index 0000000..28403ba
--- /dev/null
+++ b/fits/contrasts_fggk21_m_Child_HyC.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8ab3b8388d1732c543f130e55c4f94a43b9001e40d7158885890548bbf735687
+size 195212
diff --git a/fits/contrasts_fggk21_m_Child_HyC0.json.zip b/fits/contrasts_fggk21_m_Child_HyC0.json.zip
new file mode 100644
index 0000000..538da33
--- /dev/null
+++ b/fits/contrasts_fggk21_m_Child_HyC0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a633388dc857f8ccd317db5bbac12b4496643333ac2360ea30067ed9a3d73457
+size 87120
diff --git a/fits/contrasts_fggk21_m_Child_PCA.json.zip b/fits/contrasts_fggk21_m_Child_PCA.json.zip
new file mode 100644
index 0000000..e83cadf
--- /dev/null
+++ b/fits/contrasts_fggk21_m_Child_PCA.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:593ba911af64f6add8b37062ecd34b168a5f4d0cecbaa87709783debad4ccdfb
+size 204103
diff --git a/fits/contrasts_fggk21_m_Child_PCA0.json.zip b/fits/contrasts_fggk21_m_Child_PCA0.json.zip
new file mode 100644
index 0000000..bcfc5a4
--- /dev/null
+++ b/fits/contrasts_fggk21_m_Child_PCA0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:15fc7b9ec5ac5491a3fe331c3a4e560522af8f8c4c15ba28d28669058f629da0
+size 84188
diff --git a/fits/contrasts_fggk21_m_Child_SDC.json.zip b/fits/contrasts_fggk21_m_Child_SDC.json.zip
new file mode 100644
index 0000000..6f4aed9
--- /dev/null
+++ b/fits/contrasts_fggk21_m_Child_SDC.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3c9b850d9350b3f934c1032e0fb0c709915f5efce2ba6a4719d725ab5ff61eb3
+size 355181
diff --git a/fits/contrasts_fggk21_m_Child_SDC0.json.zip b/fits/contrasts_fggk21_m_Child_SDC0.json.zip
new file mode 100644
index 0000000..3bb4f28
--- /dev/null
+++ b/fits/contrasts_fggk21_m_Child_SDC0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a562ccc7820fd733c114b7d287334977ee9e6964188241fcc0e27fda09af8422
+size 8970
diff --git a/fits/contrasts_fggk21_m_School_HeC0.json.zip b/fits/contrasts_fggk21_m_School_HeC0.json.zip
new file mode 100644
index 0000000..83740c3
--- /dev/null
+++ b/fits/contrasts_fggk21_m_School_HeC0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:543b4c9e2ea78f5d2e0542fb725dc4ceb2e746a23d22b2bf3e5f132b4559bed1
+size 9544
diff --git a/fits/contrasts_fggk21_m_School_Helmert.json.zip b/fits/contrasts_fggk21_m_School_Helmert.json.zip
new file mode 100644
index 0000000..401b9fa
--- /dev/null
+++ b/fits/contrasts_fggk21_m_School_Helmert.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a7866e0645a862ee22da9d9a39c28cee887804b35826e667922a29bf05c08ce9
+size 258638
diff --git a/fits/contrasts_fggk21_m_School_HyC0.json.zip b/fits/contrasts_fggk21_m_School_HyC0.json.zip
new file mode 100644
index 0000000..e687192
--- /dev/null
+++ b/fits/contrasts_fggk21_m_School_HyC0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:aef8f8b368e76fa3049a46fc538d747323fdf1cbe4c48776c25924e76ff30172
+size 16065
diff --git a/fits/contrasts_fggk21_m_School_Hypo.json.zip b/fits/contrasts_fggk21_m_School_Hypo.json.zip
new file mode 100644
index 0000000..b819608
--- /dev/null
+++ b/fits/contrasts_fggk21_m_School_Hypo.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:481b3405af291b9064adc87740ad43d04c96e983d26218088fbdc306c8245e75
+size 160761
diff --git a/fits/contrasts_fggk21_m_School_PCA.json.zip b/fits/contrasts_fggk21_m_School_PCA.json.zip
new file mode 100644
index 0000000..d03e80b
--- /dev/null
+++ b/fits/contrasts_fggk21_m_School_PCA.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5b205e2ea4f750d5c4060695a57944e1a1cada278604e6fda4fec0d56b709f40
+size 215820
diff --git a/fits/contrasts_fggk21_m_School_PCA0.json.zip b/fits/contrasts_fggk21_m_School_PCA0.json.zip
new file mode 100644
index 0000000..4a1172d
--- /dev/null
+++ b/fits/contrasts_fggk21_m_School_PCA0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:13fbcb7619bfb67821b4ab1491b63d605e97db5d677419abf260cc53f6af7342
+size 18401
diff --git a/fits/contrasts_fggk21_m_School_SDC0.json.zip b/fits/contrasts_fggk21_m_School_SDC0.json.zip
new file mode 100644
index 0000000..846fee6
--- /dev/null
+++ b/fits/contrasts_fggk21_m_School_SDC0.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e97fde51cf874bde961e62fdea8b47ca6cfec5f2f5304d7227930fd0216b96b3
+size 12897
diff --git a/fits/contrasts_fggk21_m_School_SeqDiff.json.zip b/fits/contrasts_fggk21_m_School_SeqDiff.json.zip
new file mode 100644
index 0000000..b525f40
--- /dev/null
+++ b/fits/contrasts_fggk21_m_School_SeqDiff.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c00e177e669866bc7a0b6635ebb94fc6baea7bf0130690044eeee41efeed7b00
+size 169672
diff --git a/fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip b/fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip
new file mode 100644
index 0000000..08c31ad
--- /dev/null
+++ b/fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:79d81cbe3772d79c48efc3daf4d37f787376fa72e288898e4f9c726ce106eb87
+size 1547457
diff --git a/fits/contrasts_fggk21_m_cpx_1_PC.json.zip b/fits/contrasts_fggk21_m_cpx_1_PC.json.zip
new file mode 100644
index 0000000..b505732
--- /dev/null
+++ b/fits/contrasts_fggk21_m_cpx_1_PC.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c0cc419df65187f8b096f788c4905e221058595762558635396706444246b66a
+size 1903048
diff --git a/fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip b/fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip
new file mode 100644
index 0000000..1ea2d1f
--- /dev/null
+++ b/fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1db75667e854240f4d4fd4b3984c176867124593c662bbafb6b1f1f9ba259e81
+size 268934
diff --git a/fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip b/fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip
new file mode 100644
index 0000000..488003f
--- /dev/null
+++ b/fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bfc777892c686adcdde256cea19f32035d81f4a03889e2de5c560e8583230baf
+size 975764
diff --git a/fits/contrasts_fggk21_m_ovi_Helmert.json.zip b/fits/contrasts_fggk21_m_ovi_Helmert.json.zip
new file mode 100644
index 0000000..42079da
--- /dev/null
+++ b/fits/contrasts_fggk21_m_ovi_Helmert.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0c23af10a03265a3f43cf676bd6f7f2716106777613ae3a2c282cad368eaba57
+size 584
diff --git a/fits/contrasts_fggk21_m_ovi_Hypo.json.zip b/fits/contrasts_fggk21_m_ovi_Hypo.json.zip
new file mode 100644
index 0000000..db9afb3
--- /dev/null
+++ b/fits/contrasts_fggk21_m_ovi_Hypo.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ab55693e280fe55d01e6ad5d3d2b293fe7375208c0ecd67ec2a5feaaa0869b5b
+size 576
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip
new file mode 100644
index 0000000..7f08de6
--- /dev/null
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f9a3d89b601312b8d083ae3b3563509734573e169b87c1a85c45173275e04359
+size 597
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip
new file mode 100644
index 0000000..44c563e
--- /dev/null
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0cb209d4ea7c6f6063249c7a0940b90ebeb6e8aaa6766350b2290435328456c1
+size 582
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip
new file mode 100644
index 0000000..73444f6
--- /dev/null
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6eea77485b3f2c006ce134764abf3f5f9b7785db54ea2c9a59478bb2034f018b
+size 591
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip
new file mode 100644
index 0000000..774d159
--- /dev/null
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7e5a5fce5234f4c60bb5f81b5d67a12e63ada43434de65d910359c0d165a68b1
+size 575
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip
new file mode 100644
index 0000000..de16259
--- /dev/null
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8add5816393826f5417d490cb73e8753b0a877f9a2730c6f826c8c75eda0d410
+size 578
diff --git a/fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip b/fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip
new file mode 100644
index 0000000..6b877ba
--- /dev/null
+++ b/fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fe1ecb81fd2f531d5f1a88813745f7fa5580e6fbc9d4d4cc45c2cbb60a613ea1
+size 70767
diff --git a/fits/contrasts_fggk21_m_zcp_SeqD.json.zip b/fits/contrasts_fggk21_m_zcp_SeqD.json.zip
new file mode 100644
index 0000000..fee3122
--- /dev/null
+++ b/fits/contrasts_fggk21_m_zcp_SeqD.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ef04048bc526ed5956608248af91d6715e3537454b5f53ed8645c999e2b8163c
+size 8277
diff --git a/fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip b/fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip
new file mode 100644
index 0000000..f88bd11
--- /dev/null
+++ b/fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c235a44fd503722478c56ebcdf62ab4275f2c9e81cef1401838f547ccbd3a8fa
+size 1547457
diff --git a/fits/singularity_m01.json.zip b/fits/singularity_m01.json.zip
new file mode 100644
index 0000000..8e9d5c0
--- /dev/null
+++ b/fits/singularity_m01.json.zip
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:83d643ddf4d8d69dd9dc82ad88319ea37e8ceaba8be12931a07f4199db4cf193
+size 3071
diff --git a/singularity.qmd b/singularity.qmd
index 86ec373..2c05b98 100644
--- a/singularity.qmd
+++ b/singularity.qmd
@@ -22,11 +22,24 @@ using MixedModelsMakie
 const progress = isinteractive()
 ```
 
+```{julia}
+#| echo: false
+using SMLP2026: fit_or_restore
+```
+
 Fit a model for reaction time in the sleepstudy example.
 
 ```{julia}
+#| eval: false
 m01 = lmm(@formula(reaction ~ 1 + days + (1 + days|subj)), MixedModels.dataset(:sleepstudy); progress)
+```
 
+```{julia}
+#| include: false
+m01 = fit_or_restore("singularity_m01.json", @formula(reaction ~ 1 + days + (1 + days|subj)), MixedModels.dataset(:sleepstudy); progress)
+```
+
+```{julia}
 print(m01)
 ```
 
diff --git a/src/SMLP2026.jl b/src/SMLP2026.jl
index 8ef1a84..769406b 100644
--- a/src/SMLP2026.jl
+++ b/src/SMLP2026.jl
@@ -87,8 +87,9 @@ function fit_or_restore!(model::MixedModel, fname;
                 saveoptsum(f, model)
                 return nothing
             end
-        catch
-            error("Something went wrong in saving the model cache")
+        catch ex
+            @error "Something went wrong in saving the model cache"
+            rethrow(ex)
         finally
             close(zip)
         end
@@ -97,8 +98,9 @@ function fit_or_restore!(model::MixedModel, fname;
         zip = ZipFile.Reader(fname)
         try
             restoreoptsum!(model, only(zip.files); restore_kwargs...)
-        catch
-            error("Something went wrong in reading the model cache")
+        catch ex
+            @error "Something went wrong in reading the model cache"
+            rethrow(ex)
         finally
             close(zip)
         end

From 47d744097cb0e211ad6e8bb4f1143ca03c577060 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 30 Mar 2026 23:54:20 -0500
Subject: [PATCH 8/9] stable sampling for stable caching

---
 contrasts_fggk21.qmd                            |  5 +++--
 fits/contrasts_fggk21_m_Child_HeC.json.zip      |  4 ++--
 fits/contrasts_fggk21_m_Child_HeC0.json.zip     |  4 ++--
 fits/contrasts_fggk21_m_Child_HyC.json.zip      |  4 ++--
 fits/contrasts_fggk21_m_Child_HyC0.json.zip     |  4 ++--
 fits/contrasts_fggk21_m_Child_PCA.json.zip      |  4 ++--
 fits/contrasts_fggk21_m_Child_PCA0.json.zip     |  4 ++--
 fits/contrasts_fggk21_m_Child_SDC.json.zip      |  4 ++--
 fits/contrasts_fggk21_m_Child_SDC0.json.zip     |  4 ++--
 fits/contrasts_fggk21_m_School_HeC0.json.zip    |  4 ++--
 fits/contrasts_fggk21_m_School_Helmert.json.zip |  4 ++--
 fits/contrasts_fggk21_m_School_HyC0.json.zip    |  4 ++--
 fits/contrasts_fggk21_m_School_Hypo.json.zip    |  4 ++--
 fits/contrasts_fggk21_m_School_PCA.json.zip     |  4 ++--
 fits/contrasts_fggk21_m_School_PCA0.json.zip    |  4 ++--
 fits/contrasts_fggk21_m_School_SDC0.json.zip    |  4 ++--
 fits/contrasts_fggk21_m_School_SeqDiff.json.zip |  4 ++--
 fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip  |  4 ++--
 fits/contrasts_fggk21_m_cpx_1_PC.json.zip       |  4 ++--
 fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip     |  4 ++--
 fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip  |  4 ++--
 fits/contrasts_fggk21_m_ovi_Helmert.json.zip    |  4 ++--
 fits/contrasts_fggk21_m_ovi_Hypo.json.zip       |  4 ++--
 fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip    |  4 ++--
 fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip  |  4 ++--
 fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip  |  4 ++--
 fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip  |  4 ++--
 fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip |  4 ++--
 fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip     |  4 ++--
 fits/contrasts_fggk21_m_zcp_SeqD.json.zip       |  4 ++--
 fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip     |  4 ++--
 src/SMLP2026.jl                                 | 11 ++++-------
 32 files changed, 67 insertions(+), 69 deletions(-)

diff --git a/contrasts_fggk21.qmd b/contrasts_fggk21.qmd
index 485345b..863aba1 100644
--- a/contrasts_fggk21.qmd
+++ b/contrasts_fggk21.qmd
@@ -47,6 +47,7 @@ using DataFrames
 using DataFrameMacros
 using MixedModels
 using ProgressMeter
+using StableRNGs
 using Statistics
 using StatsBase
 
@@ -101,7 +102,7 @@ We extract a random sample of 500 children from the Sex (2)  x Test (5) cells of
 dat = @chain df begin
   @transform(:Sex = :Sex == "female" ? "Girls" : "Boys")
   @groupby(:Test, :Sex)
-  combine(x -> x[sample(1:nrow(x), 500), :])
+  combine(x -> x[sample(StableRNG(32123), 1:nrow(x), 500), :])
 end
 ```
 
@@ -113,7 +114,7 @@ select!(groupby(dat, :Test), :, :score => zscore => :zScore) # z-score
 ```
 
 ```{julia}
-dat2 = combine(
+combine(
   groupby(dat, [:Test, :Sex]),
   :score => mean,
   :score => std,
diff --git a/fits/contrasts_fggk21_m_Child_HeC.json.zip b/fits/contrasts_fggk21_m_Child_HeC.json.zip
index 4ae8d8e..5dd7472 100644
--- a/fits/contrasts_fggk21_m_Child_HeC.json.zip
+++ b/fits/contrasts_fggk21_m_Child_HeC.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:c57f271a20cfaae2b55456ffa762b89378280a4f1c74993eef6df39addf42a08
-size 197261
+oid sha256:337ef9bdc5a306b36da21f64f1387ae3569d21d849f5f88bf9b079c6b60b11c7
+size 193318
diff --git a/fits/contrasts_fggk21_m_Child_HeC0.json.zip b/fits/contrasts_fggk21_m_Child_HeC0.json.zip
index 43c4377..ef461ea 100644
--- a/fits/contrasts_fggk21_m_Child_HeC0.json.zip
+++ b/fits/contrasts_fggk21_m_Child_HeC0.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:ba371128e8820b0ee4a30678439b8623675301822ba6b5321c3fe3ade281ddbd
-size 81800
+oid sha256:60d50ee06fc8d64a057307fc33e5bf3ced866b0dcede7810bbffd5289d8f1571
+size 89740
diff --git a/fits/contrasts_fggk21_m_Child_HyC.json.zip b/fits/contrasts_fggk21_m_Child_HyC.json.zip
index 28403ba..22da382 100644
--- a/fits/contrasts_fggk21_m_Child_HyC.json.zip
+++ b/fits/contrasts_fggk21_m_Child_HyC.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:8ab3b8388d1732c543f130e55c4f94a43b9001e40d7158885890548bbf735687
-size 195212
+oid sha256:cf8d9c2b3829615ef9cf671ca42ec06d811a6f9649c81a0cf37b628eccada1d3
+size 244423
diff --git a/fits/contrasts_fggk21_m_Child_HyC0.json.zip b/fits/contrasts_fggk21_m_Child_HyC0.json.zip
index 538da33..84da38c 100644
--- a/fits/contrasts_fggk21_m_Child_HyC0.json.zip
+++ b/fits/contrasts_fggk21_m_Child_HyC0.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:a633388dc857f8ccd317db5bbac12b4496643333ac2360ea30067ed9a3d73457
-size 87120
+oid sha256:77a0b06b82a442368ef31d729e9ba0e4ff6a2962e406c93fbe2c37bb503afa61
+size 88113
diff --git a/fits/contrasts_fggk21_m_Child_PCA.json.zip b/fits/contrasts_fggk21_m_Child_PCA.json.zip
index e83cadf..ffbbc8c 100644
--- a/fits/contrasts_fggk21_m_Child_PCA.json.zip
+++ b/fits/contrasts_fggk21_m_Child_PCA.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:593ba911af64f6add8b37062ecd34b168a5f4d0cecbaa87709783debad4ccdfb
-size 204103
+oid sha256:9485f547ab2478b991fd90bab60cc660b51d7b779905c9ce749454f681e58497
+size 208357
diff --git a/fits/contrasts_fggk21_m_Child_PCA0.json.zip b/fits/contrasts_fggk21_m_Child_PCA0.json.zip
index bcfc5a4..289e50c 100644
--- a/fits/contrasts_fggk21_m_Child_PCA0.json.zip
+++ b/fits/contrasts_fggk21_m_Child_PCA0.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:15fc7b9ec5ac5491a3fe331c3a4e560522af8f8c4c15ba28d28669058f629da0
-size 84188
+oid sha256:a1e76577e0ad093402dcdfac433a13b80d51a7a037f52c44107a05b8aa1c72c3
+size 90934
diff --git a/fits/contrasts_fggk21_m_Child_SDC.json.zip b/fits/contrasts_fggk21_m_Child_SDC.json.zip
index 6f4aed9..ef1f2bf 100644
--- a/fits/contrasts_fggk21_m_Child_SDC.json.zip
+++ b/fits/contrasts_fggk21_m_Child_SDC.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:3c9b850d9350b3f934c1032e0fb0c709915f5efce2ba6a4719d725ab5ff61eb3
-size 355181
+oid sha256:b6858d8bcef238075e6624034d288710e8b7903cbfef0e282ac4e0f4daefdc66
+size 189828
diff --git a/fits/contrasts_fggk21_m_Child_SDC0.json.zip b/fits/contrasts_fggk21_m_Child_SDC0.json.zip
index 3bb4f28..3d866cb 100644
--- a/fits/contrasts_fggk21_m_Child_SDC0.json.zip
+++ b/fits/contrasts_fggk21_m_Child_SDC0.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:a562ccc7820fd733c114b7d287334977ee9e6964188241fcc0e27fda09af8422
-size 8970
+oid sha256:55fd71a76e521b2f18dd90c6167dcb8086aa090f2a84459e00809cea9c8f9424
+size 6484
diff --git a/fits/contrasts_fggk21_m_School_HeC0.json.zip b/fits/contrasts_fggk21_m_School_HeC0.json.zip
index 83740c3..2383c91 100644
--- a/fits/contrasts_fggk21_m_School_HeC0.json.zip
+++ b/fits/contrasts_fggk21_m_School_HeC0.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:543b4c9e2ea78f5d2e0542fb725dc4ceb2e746a23d22b2bf3e5f132b4559bed1
-size 9544
+oid sha256:1bbc8a4d83b72c2813bcceeac79068eb1ea327e012b7f6c299117ca3705b31b5
+size 10405
diff --git a/fits/contrasts_fggk21_m_School_Helmert.json.zip b/fits/contrasts_fggk21_m_School_Helmert.json.zip
index 401b9fa..b0e0a5b 100644
--- a/fits/contrasts_fggk21_m_School_Helmert.json.zip
+++ b/fits/contrasts_fggk21_m_School_Helmert.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:a7866e0645a862ee22da9d9a39c28cee887804b35826e667922a29bf05c08ce9
-size 258638
+oid sha256:2e86a1ef8251fd9ea506c9d9c88d546fa419fbc68aadfada3e1f8737dabb7a7a
+size 219000
diff --git a/fits/contrasts_fggk21_m_School_HyC0.json.zip b/fits/contrasts_fggk21_m_School_HyC0.json.zip
index e687192..df2823c 100644
--- a/fits/contrasts_fggk21_m_School_HyC0.json.zip
+++ b/fits/contrasts_fggk21_m_School_HyC0.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:aef8f8b368e76fa3049a46fc538d747323fdf1cbe4c48776c25924e76ff30172
-size 16065
+oid sha256:71c0491444386bc1056d0c8ae6fb8eaddfe02197c655a7066a254060770f308b
+size 15974
diff --git a/fits/contrasts_fggk21_m_School_Hypo.json.zip b/fits/contrasts_fggk21_m_School_Hypo.json.zip
index b819608..680c745 100644
--- a/fits/contrasts_fggk21_m_School_Hypo.json.zip
+++ b/fits/contrasts_fggk21_m_School_Hypo.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:481b3405af291b9064adc87740ad43d04c96e983d26218088fbdc306c8245e75
-size 160761
+oid sha256:a47144e2d15898f4401283368304c44a70d3f915da5c09e206a1201fca0f9210
+size 187600
diff --git a/fits/contrasts_fggk21_m_School_PCA.json.zip b/fits/contrasts_fggk21_m_School_PCA.json.zip
index d03e80b..e118bcd 100644
--- a/fits/contrasts_fggk21_m_School_PCA.json.zip
+++ b/fits/contrasts_fggk21_m_School_PCA.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:5b205e2ea4f750d5c4060695a57944e1a1cada278604e6fda4fec0d56b709f40
-size 215820
+oid sha256:597ed223ba9b78ba44698ddd4d6edbd00b6afe1b574c792e4e33b36ee79cd497
+size 243597
diff --git a/fits/contrasts_fggk21_m_School_PCA0.json.zip b/fits/contrasts_fggk21_m_School_PCA0.json.zip
index 4a1172d..a00f996 100644
--- a/fits/contrasts_fggk21_m_School_PCA0.json.zip
+++ b/fits/contrasts_fggk21_m_School_PCA0.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:13fbcb7619bfb67821b4ab1491b63d605e97db5d677419abf260cc53f6af7342
-size 18401
+oid sha256:ebab64e7ddbf1b480f77440c11286692d6610eed678732fa9ad7f5179cd65529
+size 13600
diff --git a/fits/contrasts_fggk21_m_School_SDC0.json.zip b/fits/contrasts_fggk21_m_School_SDC0.json.zip
index 846fee6..a378595 100644
--- a/fits/contrasts_fggk21_m_School_SDC0.json.zip
+++ b/fits/contrasts_fggk21_m_School_SDC0.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:e97fde51cf874bde961e62fdea8b47ca6cfec5f2f5304d7227930fd0216b96b3
-size 12897
+oid sha256:c6d31dff56162d14f71413f79048064017bf5e20b13e7c844d764adfba4a6eb4
+size 15038
diff --git a/fits/contrasts_fggk21_m_School_SeqDiff.json.zip b/fits/contrasts_fggk21_m_School_SeqDiff.json.zip
index b525f40..83d21a2 100644
--- a/fits/contrasts_fggk21_m_School_SeqDiff.json.zip
+++ b/fits/contrasts_fggk21_m_School_SeqDiff.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:c00e177e669866bc7a0b6635ebb94fc6baea7bf0130690044eeee41efeed7b00
-size 169672
+oid sha256:5830042847b696db03d222692eaa20c1a5ee0e682372bd1c8895bc5831ddf606
+size 130024
diff --git a/fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip b/fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip
index 08c31ad..345e303 100644
--- a/fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip
+++ b/fits/contrasts_fggk21_m_cpx_0_SeqDiff.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:79d81cbe3772d79c48efc3daf4d37f787376fa72e288898e4f9c726ce106eb87
-size 1547457
+oid sha256:f30efa109dd4fc0d79b7b2ff9e132053c3a35442c185313db2fc687210b19c31
+size 2094579
diff --git a/fits/contrasts_fggk21_m_cpx_1_PC.json.zip b/fits/contrasts_fggk21_m_cpx_1_PC.json.zip
index b505732..e40a5d2 100644
--- a/fits/contrasts_fggk21_m_cpx_1_PC.json.zip
+++ b/fits/contrasts_fggk21_m_cpx_1_PC.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:c0cc419df65187f8b096f788c4905e221058595762558635396706444246b66a
-size 1903048
+oid sha256:f0785a1757a8816aa38fe99d59e1f6917bc974e4b12fb39474ea2e64e9acfa6b
+size 1529606
diff --git a/fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip b/fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip
index 1ea2d1f..093bfb6 100644
--- a/fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip
+++ b/fits/contrasts_fggk21_m_cpx_1_PC_2.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:1db75667e854240f4d4fd4b3984c176867124593c662bbafb6b1f1f9ba259e81
-size 268934
+oid sha256:506e0c5cddfc3f42da9641981d1eb93655c6ab5794d4c216eadb5e269a0a02f9
+size 221686
diff --git a/fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip b/fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip
index 488003f..c27ab42 100644
--- a/fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip
+++ b/fits/contrasts_fggk21_m_cpx_1_SeqDiff.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:bfc777892c686adcdde256cea19f32035d81f4a03889e2de5c560e8583230baf
-size 975764
+oid sha256:4c0062b2cc826655653fde0f16cbd2a1d58f8d4d8ee75f17043f94af2bc90d86
+size 1104205
diff --git a/fits/contrasts_fggk21_m_ovi_Helmert.json.zip b/fits/contrasts_fggk21_m_ovi_Helmert.json.zip
index 42079da..bc1a7a6 100644
--- a/fits/contrasts_fggk21_m_ovi_Helmert.json.zip
+++ b/fits/contrasts_fggk21_m_ovi_Helmert.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:0c23af10a03265a3f43cf676bd6f7f2716106777613ae3a2c282cad368eaba57
-size 584
+oid sha256:2b613700bd6e3c012620999847ee1aa6cd012d065d94bf53328b4ac6a03a48a1
+size 645
diff --git a/fits/contrasts_fggk21_m_ovi_Hypo.json.zip b/fits/contrasts_fggk21_m_ovi_Hypo.json.zip
index db9afb3..7ccb6ed 100644
--- a/fits/contrasts_fggk21_m_ovi_Hypo.json.zip
+++ b/fits/contrasts_fggk21_m_ovi_Hypo.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:ab55693e280fe55d01e6ad5d3d2b293fe7375208c0ecd67ec2a5feaaa0869b5b
-size 576
+oid sha256:9587666ad09bfffd0f0c35e42c8274391c3f7a3f5d8d4541ff01e83d7686f6af
+size 632
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip
index 7f08de6..6bb4462 100644
--- a/fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:f9a3d89b601312b8d083ae3b3563509734573e169b87c1a85c45173275e04359
-size 597
+oid sha256:032309dbb8a3e3779f63f07b7d9532651de4a345a525c75e717ccbd347f04c5f
+size 632
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip
index 44c563e..9b02795 100644
--- a/fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff_1.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:0cb209d4ea7c6f6063249c7a0940b90ebeb6e8aaa6766350b2290435328456c1
-size 582
+oid sha256:d318464347ddfb961602429645e3c0f7cb5500ed7f6ba20e1e3a7f472b219af8
+size 605
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip
index 73444f6..cacc29a 100644
--- a/fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff_2.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:6eea77485b3f2c006ce134764abf3f5f9b7785db54ea2c9a59478bb2034f018b
-size 591
+oid sha256:3ca67d178b77ae47ca2c323ddfff4da3d2e28e022421626fdb8882689d0adb9a
+size 602
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip
index 774d159..6db8beb 100644
--- a/fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff_3.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:7e5a5fce5234f4c60bb5f81b5d67a12e63ada43434de65d910359c0d165a68b1
-size 575
+oid sha256:11f824015dcb2f9dd86ceb195edbf45fdbc28216d5443aea405a9fdd0bd217fb
+size 641
diff --git a/fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip b/fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip
index de16259..d30f6bd 100644
--- a/fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip
+++ b/fits/contrasts_fggk21_m_ovi_SeqDiff_v2.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:8add5816393826f5417d490cb73e8753b0a877f9a2730c6f826c8c75eda0d410
-size 578
+oid sha256:f9b1e14126bdaa8833b09b816a57834f853e47e4cf7a70d5d22cf846b25238dc
+size 632
diff --git a/fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip b/fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip
index 6b877ba..cb042af 100644
--- a/fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip
+++ b/fits/contrasts_fggk21_m_zcp_1_PC_2.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:fe1ecb81fd2f531d5f1a88813745f7fa5580e6fbc9d4d4cc45c2cbb60a613ea1
-size 70767
+oid sha256:07a8234d03731e588a267f8716647631f02542cd2fdb0e1a0dd7f231d831a49c
+size 78683
diff --git a/fits/contrasts_fggk21_m_zcp_SeqD.json.zip b/fits/contrasts_fggk21_m_zcp_SeqD.json.zip
index fee3122..e74efff 100644
--- a/fits/contrasts_fggk21_m_zcp_SeqD.json.zip
+++ b/fits/contrasts_fggk21_m_zcp_SeqD.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:ef04048bc526ed5956608248af91d6715e3537454b5f53ed8645c999e2b8163c
-size 8277
+oid sha256:a8eab7cae79603fb1d43059fa796ca61f8a1635b704042a5bb18aa54319a5649
+size 6754
diff --git a/fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip b/fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip
index f88bd11..f719677 100644
--- a/fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip
+++ b/fits/contrasts_fggk21_m_zcp_SeqD_2.json.zip
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:c235a44fd503722478c56ebcdf62ab4275f2c9e81cef1401838f547ccbd3a8fa
-size 1547457
+oid sha256:4d56537735fdb8f1c78127d5ead46113f957f5f2408f8db6753cb44492ef74c3
+size 2094579
diff --git a/src/SMLP2026.jl b/src/SMLP2026.jl
index 769406b..c8a2c92 100644
--- a/src/SMLP2026.jl
+++ b/src/SMLP2026.jl
@@ -82,13 +82,10 @@ function fit_or_restore!(model::MixedModel, fname;
         fit!(model; fit_kwargs...)
         zip = ZipFile.Writer(fname)
         try
-            mktempdir() do dir
-                f = ZipFile.addfile(zip, "model.json"; method=ZipFile.Deflate)
-                saveoptsum(f, model)
-                return nothing
-            end
+            f = ZipFile.addfile(zip, "model.json"; method=ZipFile.Deflate)
+            saveoptsum(f, model)
         catch ex
-            @error "Something went wrong in saving the model cache"
+            @error "Something went wrong in saving the model cache to $(fname)"
             rethrow(ex)
         finally
             close(zip)
@@ -99,7 +96,7 @@ function fit_or_restore!(model::MixedModel, fname;
         try
             restoreoptsum!(model, only(zip.files); restore_kwargs...)
         catch ex
-            @error "Something went wrong in reading the model cache"
+            @error "Something went wrong in reading the model cache from $(fname)"
             rethrow(ex)
         finally
             close(zip)

From 04b6eeb7765f92cd93a6245c36ea74d1c0406000 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 30 Mar 2026 23:58:59 -0500
Subject: [PATCH 9/9] lrt again

---
 contrasts_fggk21.qmd | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/contrasts_fggk21.qmd b/contrasts_fggk21.qmd
index 863aba1..cfc5305 100644
--- a/contrasts_fggk21.qmd
+++ b/contrasts_fggk21.qmd
@@ -517,8 +517,7 @@ VarCorr(m_zcp_1_PC_2)
 ```
 
 ```{julia}
-# XXX LRT is non monotonic in model complexity!
-# likelihoodratiotest(m_zcp_1_PC_2, m_cpx_1_PC_2)
+likelihoodratiotest(m_zcp_1_PC_2, m_cpx_1_PC_2)
 ```
 
 # Other topics