diff --git a/Project.toml b/Project.toml
index 414a15d..68e3a67 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,6 +25,7 @@ Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 MixedModelsExtras = "781a26e1-49f4-409a-8f4c-c3159d78c17e"
 MixedModelsMakie = "b12ae82c-6730-437f-aff9-d2c38332a376"
+MixedModelsSim = "d5ae56c5-23ca-4a1f-b505-9fc4796fc1fe"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
 RegressionFormulae = "545c379f-4ec2-4339-9aea-38f2fb6a8ba2"
@@ -49,6 +50,7 @@ Loess = "0.5.4, 0.6.2"
 MacroTools = "0.5"
 MixedModels = "4.20"
 MixedModelsMakie = "0.3.26"
+MixedModelsSim = "0.2.7"
 RCall = "0.13"
 RegressionFormulae = "0.1.3"
 StandardizedPredictors = "1"
diff --git a/partial_within.qmd b/partial_within.qmd
new file mode 100644
index 0000000..9251650
--- /dev/null
+++ b/partial_within.qmd
@@ -0,0 +1,124 @@
+---
+title: "Partially-within subjects designs"
+jupyter: julia-1.9
+---
+
+Begin by loading the packages to be used.
+
+```{julia}
+#| code-fold: true
+#| output: false
+using AlgebraOfGraphics
+using CairoMakie
+using DataFrames
+using MixedModels
+using MixedModlesMakie
+using MixedModelsSim
+using ProgressMeter
+using Random
+
+CairoMakie.activate!(; type="svg")
+
+ProgressMeter.ijulia_behavior(:clear)
+```
+
+```{julia}
+#| code-fold: true
+n_subj = 40
+n_item = 3
+# things are expressed as "between", so "within subjects" is "between items"
+item_btwn = Dict(:frequency => ["high", "medium", "low"])
+design = simdat_crossed(MersenneTwister(42), n_subj, n_item;
+                        item_btwn = item_btwn)
+design = DataFrame(design)
+```
+
+```{julia}
+#| code-fold: true
+unique!(select(design, :item, :frequency))
+```
+
+```{julia}
+#| code-fold: true
+m0 = let contrasts, form
+    contrasts = Dict(:frequency => HelmertCoding(base="high"))
+    form = @formula(dv ~ 1 + frequency +
+                    (1 + frequency | subj))
+    fit(MixedModel, form, design; contrasts)
+end
+```
+
+```{julia}
+#| code-fold: true
+corrmat = [ 1    0.1 -0.2
+            0.1  1    0.1
+           -0.2  0.1  1 ]
+re_subj = create_re(1.2, 1.5, -1.5; corrmat)
+```
+
+```{julia}
+#| code-fold: true
+θ = createθ(m0; subj=re_subj)
+```
+
+```{julia}
+#| code-fold: true
+σ = 1;
+β = [1.0, -3, -2];
+```
+
+```{julia}
+#| code-fold: true
+fit!(simulate!(m0; θ, β, σ))
+```
+
+
+```{julia}
+#| code-fold: true
+shrinkageplot(m0)
+```
+
+```{julia}
+#| code-fold: true
+caterpillar(m0; orderby=3)
+```
+
+
+```{julia}
+#| code-fold: true
+design[!, :dv] .= response(m0)
+```
+
+```{julia}
+#| code-fold: true
+design_partial = filter(design) do row
+    subj = parse(Int, row.subj[2:end])
+    item = parse(Int, row.item[2:end])
+    # for even-numbered subjects, we keep all conditions
+    # for odd-numbered subjects, we keep only the two "odd" items,
+    # i.e. the first and last conditions
+    return iseven(subj) || isodd(item)
+end
+sort!(unique!(select(design_partial, :subj, :frequency)), :subj)
+```
+
+```{julia}
+#| code-fold: true
+
+m1 = let contrasts, form
+    contrasts = Dict(:frequency => HelmertCoding(base="high"))
+    form = @formula(dv ~ 1 + frequency +
+                    (1 + frequency | subj))
+    fit(MixedModel, form, design_partial; contrasts)
+end
+```
+
+```{julia}
+#| code-fold: true
+shrinkageplot(m1)
+```
+
+```{julia}
+#| code-fold: true
+caterpillar(m1; orderby=3)
+```