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 438ae9f..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)
 ```
@@ -150,7 +168,7 @@ confint(samp; method=:shortest)
 ```{julia}
 draw(
   data(samp.β) * mapping(:β; color=:coefname) * AoG.density();
-  figure=(; resolution=(800, 450)),
+  figure=(; size=(800, 450)),
 )
 ```
 
@@ -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..cfc5305 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
 
@@ -56,6 +57,11 @@ using SMLP2026: dataset
 progress = isinteractive()
 ```
 
+```{julia}
+#| echo: false
+using SMLP2026: fit_or_restore
+```
+
 ## Readme for `dataset("fggk21")`
 
 Number of scores: 525126
@@ -96,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
 ```
 
@@ -108,7 +114,7 @@ select!(groupby(dat, :Test), :, :score => zscore => :zScore) # z-score
 ```
 
 ```{julia}
-dat2 = combine(
+combine(
   groupby(dat, [:Test, :Sex]),
   :score => mean,
   :score => std,
@@ -160,19 +166,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 +202,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 +222,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 +269,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 +312,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 +342,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 +400,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 +444,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 +480,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 +499,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)
 ```
@@ -430,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)
@@ -438,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)
 ```
@@ -459,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)
@@ -467,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)
 ```
@@ -511,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)
@@ -519,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)
 ```
@@ -538,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/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 b984153..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)
 ```
@@ -926,7 +987,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 +999,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/fits/.gitattributes b/fits/.gitattributes
new file mode 100644
index 0000000..24df8a2
--- /dev/null
+++ b/fits/.gitattributes
@@ -0,0 +1,2 @@
+*.zip filter=lfs diff=lfs merge=lfs -text
+*.arrow filter=lfs diff=lfs merge=lfs -text
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_fggk21_m_Child_HeC.json.zip b/fits/contrasts_fggk21_m_Child_HeC.json.zip
new file mode 100644
index 0000000..5dd7472
--- /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:337ef9bdc5a306b36da21f64f1387ae3569d21d849f5f88bf9b079c6b60b11c7
+size 193318
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..ef461ea
--- /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:60d50ee06fc8d64a057307fc33e5bf3ced866b0dcede7810bbffd5289d8f1571
+size 89740
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..22da382
--- /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:cf8d9c2b3829615ef9cf671ca42ec06d811a6f9649c81a0cf37b628eccada1d3
+size 244423
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..84da38c
--- /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:77a0b06b82a442368ef31d729e9ba0e4ff6a2962e406c93fbe2c37bb503afa61
+size 88113
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..ffbbc8c
--- /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:9485f547ab2478b991fd90bab60cc660b51d7b779905c9ce749454f681e58497
+size 208357
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..289e50c
--- /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:a1e76577e0ad093402dcdfac433a13b80d51a7a037f52c44107a05b8aa1c72c3
+size 90934
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..ef1f2bf
--- /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:b6858d8bcef238075e6624034d288710e8b7903cbfef0e282ac4e0f4daefdc66
+size 189828
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..3d866cb
--- /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:55fd71a76e521b2f18dd90c6167dcb8086aa090f2a84459e00809cea9c8f9424
+size 6484
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..2383c91
--- /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:1bbc8a4d83b72c2813bcceeac79068eb1ea327e012b7f6c299117ca3705b31b5
+size 10405
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..b0e0a5b
--- /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:2e86a1ef8251fd9ea506c9d9c88d546fa419fbc68aadfada3e1f8737dabb7a7a
+size 219000
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..df2823c
--- /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:71c0491444386bc1056d0c8ae6fb8eaddfe02197c655a7066a254060770f308b
+size 15974
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..680c745
--- /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:a47144e2d15898f4401283368304c44a70d3f915da5c09e206a1201fca0f9210
+size 187600
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..e118bcd
--- /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:597ed223ba9b78ba44698ddd4d6edbd00b6afe1b574c792e4e33b36ee79cd497
+size 243597
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..a00f996
--- /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:ebab64e7ddbf1b480f77440c11286692d6610eed678732fa9ad7f5179cd65529
+size 13600
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..a378595
--- /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:c6d31dff56162d14f71413f79048064017bf5e20b13e7c844d764adfba4a6eb4
+size 15038
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..83d21a2
--- /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: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
new file mode 100644
index 0000000..345e303
--- /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: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
new file mode 100644
index 0000000..e40a5d2
--- /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: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
new file mode 100644
index 0000000..093bfb6
--- /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: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
new file mode 100644
index 0000000..c27ab42
--- /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:4c0062b2cc826655653fde0f16cbd2a1d58f8d4d8ee75f17043f94af2bc90d86
+size 1104205
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..bc1a7a6
--- /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:2b613700bd6e3c012620999847ee1aa6cd012d065d94bf53328b4ac6a03a48a1
+size 645
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..7ccb6ed
--- /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:9587666ad09bfffd0f0c35e42c8274391c3f7a3f5d8d4541ff01e83d7686f6af
+size 632
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..6bb4462
--- /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: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
new file mode 100644
index 0000000..9b02795
--- /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: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
new file mode 100644
index 0000000..cacc29a
--- /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: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
new file mode 100644
index 0000000..6db8beb
--- /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: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
new file mode 100644
index 0000000..d30f6bd
--- /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: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
new file mode 100644
index 0000000..cb042af
--- /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:07a8234d03731e588a267f8716647631f02542cd2fdb0e1a0dd7f231d831a49c
+size 78683
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..e74efff
--- /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: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
new file mode 100644
index 0000000..f719677
--- /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:4d56537735fdb8f1c78127d5ead46113f957f5f2408f8db6753cb44492ef74c3
+size 2094579
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/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.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/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/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/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/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/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/kkl15.qmd b/kkl15.qmd
index 534f0b4..b269498 100644
--- a/kkl15.qmd
+++ b/kkl15.qmd
@@ -48,6 +48,11 @@ using SMLP2026: dataset
 progress = isinteractive()
 ```
 
+```{julia}
+#| echo: false
+using SMLP2026: fit_or_restore, bootstrap_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
@@ -397,10 +473,17 @@ Here we
 We generate 2500 samples for the 15 model parameters (4 fixed effect, 7 VCs, 15 CPs, and 1 residual).
 
 ```{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
 ```
@@ -432,7 +515,7 @@ draw(
   data(tbl) *
   mapping(:σ => "Residual") *
   density();
-  figure=(; resolution=(800, 400)),
+  figure=(; size=(800, 400)),
 )
 ```
 
@@ -465,7 +548,7 @@ draw(
     "Experimental effects",
   ) *
   density();
-  figure=(; resolution=(800, 350)),
+  figure=(; size=(800, 350)),
 )
 ```
 
@@ -484,7 +567,7 @@ draw(
     "Variance components",
   ) *
   density();
-  figure=(; resolution=(800, 350)),
+  figure=(; size=(800, 350)),
 )
 ```
 
@@ -506,7 +589,7 @@ draw(
     "Correlation parameters",
   ) *
   density();
-  figure=(; resolution=(800, 350)),
+  figure=(; size=(800, 350)),
 )
 ```
 
diff --git a/kwdyz11.qmd b/kwdyz11.qmd
index 89c43d3..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
 ```
 
@@ -116,7 +116,7 @@ boxplot(
       ],
     ),
   ),
-  figure=(; resolution=(800, 300)),
+  figure=(; size=(800, 300)),
 )
 ```
 
@@ -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)
 ```
@@ -247,7 +259,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 +272,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..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
 ```
@@ -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)),
 )
 ```
 
@@ -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)
 ```
@@ -562,7 +585,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 +600,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/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 4c4c84b..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.
 
@@ -117,7 +127,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 +147,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 +156,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 +165,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 +202,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 +211,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/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/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/sleepstudy.qmd b/sleepstudy.qmd
index 72224be..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.
@@ -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
@@ -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}
@@ -227,7 +244,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 +286,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 +318,7 @@ hist(
   tbl.ρ1;
   bins=40,
   axis=(; xlabel="Estimated correlation of the random effects"),
-  figure=(; resolution=(500, 500)),
+  figure=(; size=(500, 500)),
 )
 ```
 
@@ -312,7 +329,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 +361,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..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
@@ -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)
@@ -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.
 
@@ -176,7 +194,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
diff --git a/src/SMLP2026.jl b/src/SMLP2026.jl
index 08ce54b..c8a2c92 100644
--- a/src/SMLP2026.jl
+++ b/src/SMLP2026.jl
@@ -6,8 +6,10 @@ using DataFrames
 using Dates
 using Downloads
 using Markdown
+using MixedModels
 using MixedModelsDatasets
 using PooledArrays
+using Random
 using Scratch
 using SHA
 using TypedTables
@@ -39,6 +41,89 @@ end
 
 export GENRES,
     age_at_event,
-    tagpad
+    tagpad,
+    fit_or_restore,
+    fit_or_restore!
+
+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"
+    end
+
+    return path
+end
+
+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,Any}(), kwargs...)
+    model = MixedModel(args...; contrasts)
+    return fit_or_restore!(model, fname; kwargs...)
+end
+
+function fit_or_restore!(model::MixedModel, fname;
+                         force=false, restore_kwargs=(; atol=1e-8), fit_kwargs...)
+    fname = _normalize_model_cache_path(fname)
+    @debug "cache path: $(fname)"
+    if !isfile(fname) || force
+        @debug "fitting model"
+        fit!(model; fit_kwargs...)
+        zip = ZipFile.Writer(fname)
+        try
+            f = ZipFile.addfile(zip, "model.json"; method=ZipFile.Deflate)
+            saveoptsum(f, model)
+        catch ex
+            @error "Something went wrong in saving the model cache to $(fname)"
+            rethrow(ex)
+        finally
+            close(zip)
+        end
+    else
+        @debug "restoring from cache"
+        zip = ZipFile.Reader(fname)
+        try
+            restoreoptsum!(model, only(zip.files); restore_kwargs...)
+        catch ex
+            @error "Something went wrong in reading the model cache from $(fname)"
+            rethrow(ex)
+        finally
+            close(zip)
+        end
+    end
+
+    return model
+end
+
+# 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
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.