From f2caaf461dd66d6645cc7c8fcb2e4bbe25b2173b Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 27 Aug 2026 09:14:49 +0200 Subject: [PATCH 01/10] initial plotting script --- Project.toml | 2 +- scripts/plotting.jl | 31 +++++++++++++++++++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) create mode 100644 scripts/plotting.jl diff --git a/Project.toml b/Project.toml index d34a2fe..ddc2508 100644 --- a/Project.toml +++ b/Project.toml @@ -67,7 +67,7 @@ MacroTools = "0.5" MixedModels = "5" MixedModelsDatasets = "0.2" MixedModelsExtras = "2" -MixedModelsMakie = "0.4" +MixedModelsMakie = "0.4.15" MixedModelsSerialization = "0.2" MixedModelsSim = "0.2.7" PrettyTables = "3" diff --git a/scripts/plotting.jl b/scripts/plotting.jl new file mode 100644 index 0000000..15d6146 --- /dev/null +++ b/scripts/plotting.jl @@ -0,0 +1,31 @@ +using CairoMakie # plotting backend needed for the other plotting stuff +using AlgebraOfGraphics # ggplot2 type interface +using MixedModelsMakie # mixed models specials for plotting +using DataFrames +using Effects # like effects or emmeans +using MixedModels # you know this one +using MixedModelsDatasets # for the data + +kb07 = dataset(:kb07) +insteval = dataset(:insteval) +ml1m = dataset(:ml1m) + +mkb07 = lmm(@formula(rt_trunc ~ 1 + spkr * prec * load + + (1 + spkr + prec + load | subj) + + (1 + spkr + prec + load | item)), + kb07; + contrasts=Dict(:spkr => EffectsCoding(), + :prec => EffectsCoding(), + :load => EffectsCoding())) + +mm = lmm(@formula(Y ~ 1 + (1|G) + (1|H)), dataset(:ml1m)) +mi = lmm(@formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept)), dataset(:insteval)) + +upsetplot(kb07; cols=Not([:subj, :item])) +upsetplot(mkb07, :subj) +upsetplot(mkb07, :item) +nestingplot(mkb07) +nestingtable(mkb07) +filter(:count => ==(0), nestingtable(mkb07)) +nestingplot(mi) +nestingplot(mm) From 22f89029dcee7a4b7065677fa3955b5bb450ae93 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 27 Aug 2026 09:32:24 +0200 Subject: [PATCH 02/10] notes --- scripts/plotting.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/scripts/plotting.jl b/scripts/plotting.jl index 15d6146..8ecb82a 100644 --- a/scripts/plotting.jl +++ b/scripts/plotting.jl @@ -18,14 +18,17 @@ mkb07 = lmm(@formula(rt_trunc ~ 1 + spkr * prec * load :prec => EffectsCoding(), :load => EffectsCoding())) -mm = lmm(@formula(Y ~ 1 + (1|G) + (1|H)), dataset(:ml1m)) mi = lmm(@formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept)), dataset(:insteval)) +mm = lmm(@formula(Y ~ 1 + (1|G) + (1|H)), dataset(:ml1m)) -upsetplot(kb07; cols=Not([:subj, :item])) -upsetplot(mkb07, :subj) -upsetplot(mkb07, :item) nestingplot(mkb07) nestingtable(mkb07) -filter(:count => ==(0), nestingtable(mkb07)) +filter(:count => iszero, nestingtable(mkb07)) +nestingstructure(mkb07) nestingplot(mi) nestingplot(mm) + +# BUG HERE. I WILL FIX +upsetplot(kb07; cols=Not([:subj, :item])) +upsetplot(mkb07, :subj) +upsetplot(mkb07, :item) From 2bf6c7f41555464ec1b8b3f7ef6773c74948e986 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 27 Aug 2026 09:45:55 +0200 Subject: [PATCH 03/10] wip --- scripts/plotting.jl | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/scripts/plotting.jl b/scripts/plotting.jl index 8ecb82a..83ab615 100644 --- a/scripts/plotting.jl +++ b/scripts/plotting.jl @@ -5,6 +5,7 @@ using DataFrames using Effects # like effects or emmeans using MixedModels # you know this one using MixedModelsDatasets # for the data +using Random # for the random number generator kb07 = dataset(:kb07) insteval = dataset(:insteval) @@ -15,7 +16,7 @@ mkb07 = lmm(@formula(rt_trunc ~ 1 + spkr * prec * load + (1 + spkr + prec + load | item)), kb07; contrasts=Dict(:spkr => EffectsCoding(), - :prec => EffectsCoding(), + :prec => EffectsCoding(base="maintain"), :load => EffectsCoding())) mi = lmm(@formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept)), dataset(:insteval)) @@ -32,3 +33,31 @@ nestingplot(mm) upsetplot(kb07; cols=Not([:subj, :item])) upsetplot(mkb07, :subj) upsetplot(mkb07, :item) + +mkb07_small = lmm(@formula(rt_trunc ~ 1 + spkr * prec * load + + (1 + spkr | subj) + + (1 + load | item)), + kb07; + contrasts=Dict(:spkr => EffectsCoding(), + :prec => EffectsCoding(base="maintain"), + :load => EffectsCoding())) + +bkb07 = parametricbootstrap(MersenneTwister(2708), 500, mkb07_small) + +coefplot(mkb07) + +coefplot(mkb07, mkb07_small; + show_intercept=false, + labels=["big", "small"]) + +coefplot(mkb07_small, bkb07; + show_intercept=false, + labels=["wald", "boot"]) + +ridgeplot(bkb07; show_intercept=false) + +ridgeplot(bkb07; ptype=:σ) + +ridgeplot(bkb07; + ptype=:sigma, + group=:subj) From a1bff81d11fc56da516d2ab3cbdd67b7ee23e9ef Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 27 Aug 2026 10:14:14 +0200 Subject: [PATCH 04/10] wip --- .../hypothesis_drop_col.jl | 0 scripts/plotting.jl | 43 ++++++++++++++++++- 2 files changed, 42 insertions(+), 1 deletion(-) rename hypothesis_drop_col.jl => scripts/hypothesis_drop_col.jl (100%) diff --git a/hypothesis_drop_col.jl b/scripts/hypothesis_drop_col.jl similarity index 100% rename from hypothesis_drop_col.jl rename to scripts/hypothesis_drop_col.jl diff --git a/scripts/plotting.jl b/scripts/plotting.jl index 83ab615..515e57e 100644 --- a/scripts/plotting.jl +++ b/scripts/plotting.jl @@ -42,7 +42,7 @@ mkb07_small = lmm(@formula(rt_trunc ~ 1 + spkr * prec * load :prec => EffectsCoding(base="maintain"), :load => EffectsCoding())) -bkb07 = parametricbootstrap(MersenneTwister(2708), 500, mkb07_small) +bkb07 = parametricbootstrap(MersenneTwister(2708), 5000, mkb07_small) coefplot(mkb07) @@ -61,3 +61,44 @@ ridgeplot(bkb07; ptype=:σ) ridgeplot(bkb07; ptype=:sigma, group=:subj) + +ridgeplot(bkb07; + ptype=:rho) + +ridgeplot(bkb07; + ptype=:rho, + histogram=true) + +ridgeplot(bkb07; + ptype=:rho, + histogram=true, + bins=20) + +ridgeplot(bkb07; + ptype=:θ) + +ridgeplot(bkb07; + ptype=:θ, + histogram=true) + +eff = effects(Dict(:spkr => ["old", "new"], + :prec => ["break", "maintain"], + :load => ["yes", "no"]), + mkb07) + +plt = data(eff) * mapping(:spkr, :rt_trunc; + color=:load, + col=:prec) * visual(ScatterLines) + +plt = data(eff) * mapping(:spkr; color=:load, col=:prec) * + (mapping(:rt_trunc) * visual(ScatterLines) + + mapping(:lower, :upper) * visual(Band, alpha=0.3)) + +draw(plt; + figure=(; title="kb07 model"), + axis=(; ylabel="Reaction Time (ms)", + xlabel="Speaker"), + legend=(; position=:bottom, + titleposition=:left, + framevisible=false)) + From 8a397245b84875e12f8ee1bc2c21b86074e87edb Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 27 Aug 2026 10:19:54 +0200 Subject: [PATCH 05/10] wip --- scripts/plotting.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/scripts/plotting.jl b/scripts/plotting.jl index 515e57e..e1a354b 100644 --- a/scripts/plotting.jl +++ b/scripts/plotting.jl @@ -5,6 +5,7 @@ using DataFrames using Effects # like effects or emmeans using MixedModels # you know this one using MixedModelsDatasets # for the data +using MixedModelsExtras # extra things that people ask for but we don't really endorse using Random # for the random number generator kb07 = dataset(:kb07) @@ -102,3 +103,7 @@ draw(plt; titleposition=:left, framevisible=false)) + +DataFrame(ictable(mkb07, mkb07_small)) + +show(stdout, MIME("text/latex"), DataFrame(ictable(mkb07, mkb07_small))) From 89f0be89e81c155d44e60d3cf264554d9a846345 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Fri, 28 Aug 2026 08:46:48 +0200 Subject: [PATCH 06/10] labeled shrinkageplot --- Project.toml | 2 +- lmm-intro/sleepstudy.qmd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index ddc2508..c7fdc7f 100644 --- a/Project.toml +++ b/Project.toml @@ -67,7 +67,7 @@ MacroTools = "0.5" MixedModels = "5" MixedModelsDatasets = "0.2" MixedModelsExtras = "2" -MixedModelsMakie = "0.4.15" +MixedModelsMakie = "0.4.17" MixedModelsSerialization = "0.2" MixedModelsSim = "0.2.7" PrettyTables = "3" diff --git a/lmm-intro/sleepstudy.qmd b/lmm-intro/sleepstudy.qmd index cbaae0a..4f34284 100644 --- a/lmm-intro/sleepstudy.qmd +++ b/lmm-intro/sleepstudy.qmd @@ -262,7 +262,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(; size=(500, 500)), m1) +shrinkageplot!(Figure(; size=(500, 500)), m1; labels=:auto) ``` ::: {.callout-note} From 73593b49173f582fa19448573f5e107f0e2a5871 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Fri, 28 Aug 2026 09:14:24 +0200 Subject: [PATCH 07/10] bug fixed --- scripts/plotting.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/plotting.jl b/scripts/plotting.jl index e1a354b..3b10d1e 100644 --- a/scripts/plotting.jl +++ b/scripts/plotting.jl @@ -30,7 +30,6 @@ nestingstructure(mkb07) nestingplot(mi) nestingplot(mm) -# BUG HERE. I WILL FIX upsetplot(kb07; cols=Not([:subj, :item])) upsetplot(mkb07, :subj) upsetplot(mkb07, :item) From 683cc3c8106ba66a4a736ff9580d3308ed078720 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Fri, 28 Aug 2026 09:58:02 +0200 Subject: [PATCH 08/10] enable the correct chunks --- glmm/glmm.qmd | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/glmm/glmm.qmd b/glmm/glmm.qmd index 8e4fc85..8e1e721 100644 --- a/glmm/glmm.qmd +++ b/glmm/glmm.qmd @@ -264,29 +264,27 @@ contrasts = Dict( :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(MixedModel, form, contra, dist; nAGQ, contrasts, progress) + fit(MixedModel, form, contra, dist; contrasts, progress) end ``` ```{julia} -#| include: false +#| echo: false contrasts = Dict( :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) + fit_or_restore("glmm_gm1.json", MixedModel, form, contra, dist; contrasts, progress) end ``` @@ -327,12 +325,12 @@ gm2 = let urban + (1 | dist) ) - fit(MixedModel, form, contra, dist; nAGQ, contrasts, progress) + fit(MixedModel, form, contra, dist; contrasts, progress) end ``` ```{julia} -#| include: false +#| echo: false gm2 = let form = @formula( use ~ @@ -343,7 +341,7 @@ gm2 = let urban + (1 | dist) ) - fit_or_restore("glmm_gm2.json", MixedModel, form, contra, dist; nAGQ, contrasts, progress) + fit_or_restore("glmm_gm2.json", MixedModel, form, contra, dist; contrasts, progress) end ``` @@ -371,6 +369,8 @@ At present the calculation of the `geomdof` as `sum(influence(m))` is not correc ### Using `urban&dist` as a grouping factor + + It turns out that there can be more difference between urban and rural settings within the same political district than there is between districts. To model this difference we build a model with `urban&dist` as a grouping factor. @@ -386,12 +386,12 @@ gm3 = let urban + (1 | urban & dist) ) - fit(MixedModel, form, contra, dist; nAGQ, contrasts, progress) + fit(MixedModel, form, contra, dist; contrasts, progress) end ``` ```{julia} -#| include: false +#| echo: false gm3 = let form = @formula( use ~ @@ -402,7 +402,7 @@ gm3 = let urban + (1 | urban & dist) ) - fit_or_restore("glmm_gm3.json", MixedModel, form, contra, dist; nAGQ, contrasts, progress) + fit_or_restore("glmm_gm3.json", MixedModel, form, contra, dist; contrasts, progress) end ``` From 560e1ec9a557f70673687bd759553c44880e7d69 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Fri, 28 Aug 2026 09:58:13 +0200 Subject: [PATCH 09/10] additional effects plots --- glmm/glmm.qmd | 48 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 47 insertions(+), 1 deletion(-) diff --git a/glmm/glmm.qmd b/glmm/glmm.qmd index 8e1e721..cb36e2d 100644 --- a/glmm/glmm.qmd +++ b/glmm/glmm.qmd @@ -443,7 +443,53 @@ using Effects design = Dict( :children => ["Y", "N"], :urban => ["Y", "N"], :age => [0.0] ) -preds = effects(design, gm3; invlink=AutoInvLink()) +preds = effects(design, gm3) +``` + +We can plot this with a few more values for age: + +```{julia} +design = Dict( + :children => ["Y", "N"], + :urban => ["Y", "N"], + :age => -10:10 +) +preds = effects(design, gm3; level=0.95) +base = data(preds) * mapping(:age; + color=:children, + col=:urban => renamer(["N" => "rural", "Y" => "urban"])) + + +lines = mapping("use: Y") * visual(Lines) +bands = mapping(:lower, :upper) * visual(Band; alpha=0.3) + + +draw(base * (lines + bands), + legend = (; position = :top, + framevisible=false), + axis=(; ylabel="Log odds of contraception use", + xlabel="Centered age")) +``` + +We can also plot this on the response scale, i.e. the probability scale: + +```{julia} +preds = effects(design, gm3; invlink=AutoInvLink(), level=0.95) +base = data(preds) * mapping(:age; + color=:children, + col=:urban => renamer(["N" => "rural", "Y" => "urban"])) + + +lines = mapping("use: Y") * visual(Lines) +bands = mapping(:lower, :upper) * visual(Band; alpha=0.3) + + +draw(base * (lines + bands); + legend = (; position = :top, + framevisible=false), + axis=(; ylabel="Probability of contraception use", + xlabel="Centered age", + limits=(nothing, (0, 1)))) ``` From 496419b7351922ea96197d2629b0e032b77b621a Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Fri, 28 Aug 2026 10:48:49 +0200 Subject: [PATCH 10/10] district urban vs rural distribution --- glmm/glmm.qmd | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/glmm/glmm.qmd b/glmm/glmm.qmd index cb36e2d..ce6339c 100644 --- a/glmm/glmm.qmd +++ b/glmm/glmm.qmd @@ -35,6 +35,7 @@ using MixedModels using MixedModelsMakie using MixedModelsDatasets: dataset using SMLP2026: fit_or_restore +using Statistics const progress = isinteractive() ``` @@ -369,9 +370,41 @@ At present the calculation of the `geomdof` as `sum(influence(m))` is not correc ### Using `urban&dist` as a grouping factor - - It turns out that there can be more difference between urban and rural settings within the same political district than there is between districts. + +```{julia} +dist_mean = combine(groupby(contra, :dist), + :use => (x -> mean(x .== "Y")) => "dist_mean") +plt = data(dist_mean) * mapping(:dist_mean => "Distribution of district means") * AlgebraOfGraphics.density() +draw(plt) +``` + +```{julia} +dist_urban_mean = combine(groupby(contra, [:dist, :urban]), + :use => (x -> mean(x .== "Y")) => "dist_urban_mean") +plt = data(dist_urban_mean) * mapping(:dist_urban_mean => "Distribution of district × urban means"; color=:urban) * AlgebraOfGraphics.density() +draw(plt) +``` + +```{julia} +dum_sorter = combine(groupby(dist_urban_mean, :dist), + :dist_urban_mean => diff => :urban_rural_diff) +transform!(dum_sorter, :urban_rural_diff => ByRow(abs); renamecols=false) +all_dists = DataFrame(; dist=unique(dist_urban_mean.dist)) +dum_sorter = leftjoin(all_dists, dum_sorter; on=:dist) +transform!(dum_sorter, + :urban_rural_diff => ByRow(x -> coalesce(x, 0)); + renamecols=false) +sort!(dum_sorter, :urban_rural_diff) + +plt = data(dist_urban_mean) * + mapping(:dist_urban_mean => "Proportion contraception use", + :dist => sorter(dum_sorter.dist) => "District") * + (mapping(; color=:urban) * visual(Scatter) + + mapping(; group=:dist) * visual(Lines)) +draw(plt; figure=(;size=(500, 950), title="Distribution of district × urban means"), legend=(; position=:top)) +``` + To model this difference we build a model with `urban&dist` as a grouping factor. ```{julia}