diff --git a/Project.toml b/Project.toml index d34a2fe..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" +MixedModelsMakie = "0.4.17" MixedModelsSerialization = "0.2" MixedModelsSim = "0.2.7" PrettyTables = "3" diff --git a/glmm/glmm.qmd b/glmm/glmm.qmd index 8e4fc85..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() ``` @@ -264,29 +265,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 +326,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 +342,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 ``` @@ -372,6 +371,40 @@ 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} @@ -386,12 +419,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 +435,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 ``` @@ -443,7 +476,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)))) ``` 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} 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 new file mode 100644 index 0000000..3b10d1e --- /dev/null +++ b/scripts/plotting.jl @@ -0,0 +1,108 @@ +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 +using MixedModelsExtras # extra things that people ask for but we don't really endorse +using Random # for the random number generator + +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(base="maintain"), + :load => EffectsCoding())) + +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)) + +nestingplot(mkb07) +nestingtable(mkb07) +filter(:count => iszero, nestingtable(mkb07)) +nestingstructure(mkb07) +nestingplot(mi) +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), 5000, 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) + +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)) + + +DataFrame(ictable(mkb07, mkb07_small)) + +show(stdout, MIME("text/latex"), DataFrame(ictable(mkb07, mkb07_small)))