From 3f495319d75f55beb80dcb353114ecdd7dafc097 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 28 May 2023 12:25:50 -0500 Subject: [PATCH 1/2] start work on visualization appendix --- _quarto.yml | 1 + viz.qmd | 60 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+) create mode 100644 viz.qmd diff --git a/_quarto.yml b/_quarto.yml index 0a4bf16..acc1983 100644 --- a/_quarto.yml +++ b/_quarto.yml @@ -20,6 +20,7 @@ book: - datatables.qmd - linalg.qmd - aGHQ.qmd + - viz.qmd execute-dir: project diff --git a/viz.qmd b/viz.qmd new file mode 100644 index 0000000..59a2a7f --- /dev/null +++ b/viz.qmd @@ -0,0 +1,60 @@ +--- +jupyter: julia-1.9 +--- + +# Vizualization of Mixed Models {#sec-viz} + + +```{julia} +#| code-fold: true +const hide_progress = true + +using AlgebraOfGraphics +using CairoMakie +using DataFrames +if contains(first(Sys.cpu_info()).model, "Intel") + using MKL +end +using MixedModels +using MixedModelsMakie +using ProgressMeter +using Random + +ProgressMeter.ijulia_behavior(:clear); +CairoMakie.activate!(; type="svg"); +``` + +```{julia} +kb07 = DataFrame(MixedModels.dataset("kb07")) +contrasts = Dict(:item => Grouping(), + :subj => Grouping(), + :spkr => EffectsCoding(), + :prec => EffectsCoding(), + :load => EffectsCoding()) + +fm1 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1|item) + (1|subj)), kb07; contrasts, progress=true, thin=1) +``` + + +```{julia} +let f = Figure() + ax1 = Axis(f[1,1], title="Subject") + caterpillar!(ax1, fm1, :subj) + # caterpillar!(ax2, fm1, :subj) + f +end + + let f = Figure(; resolution=(1000, 600)) + gl = f[1,1] = GridLayout() + re = ranefinfo(m2) + qqcaterpillar!(gl, re[:item]) + Label(gl[end+1, :], "Item"; font=:bold) + gl = f[1,2] = GridLayout() + qqcaterpillar!(gl, re[:subj]) + Label(gl[end+1, :], "Subject"; font=:bold) + Label(f[0, :], "Conditional Modes") + colsize!(f.layout, 1, Auto(0.5)) + save(joinpath(OUTDIR, "qqcat_kb07_joint.png"), f) + f + end +``` From d01c2398ec9a3ff9bca4edfe69576896eaeb36f5 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 29 May 2023 00:04:38 -0500 Subject: [PATCH 2/2] notes --- Project.toml | 4 +++- viz.jl | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++ viz.qmd | 49 +++++++++++++++++++++++++++++++------------------ 3 files changed, 85 insertions(+), 19 deletions(-) create mode 100644 viz.jl diff --git a/Project.toml b/Project.toml index 968249c..0c1d3ac 100644 --- a/Project.toml +++ b/Project.toml @@ -22,6 +22,7 @@ GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" +MixedModelsExtras = "781a26e1-49f4-409a-8f4c-c3159d78c17e" MixedModelsMakie = "b12ae82c-6730-437f-aff9-d2c38332a376" NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" @@ -56,7 +57,8 @@ Effects = "0.1,1" FreqTables = "0.4" GLM = "1" MixedModels = "4,5" -MixedModelsMakie = "0.3.13" +MixedModelsExtras = "1.1" +MixedModelsMakie = "0.3.20" NLopt = "0.6" ProgressMeter = "1.7" RCall = "0.13" diff --git a/viz.jl b/viz.jl new file mode 100644 index 0000000..3a64044 --- /dev/null +++ b/viz.jl @@ -0,0 +1,51 @@ +const hide_progress = true + +using AlgebraOfGraphics +using CairoMakie +using DataFrames +if contains(first(Sys.cpu_info()).model, "Intel") + using MKL +end +using MixedModels +using MixedModelsMakie +using ProgressMeter +using Random + +kb07 = DataFrame(MixedModels.dataset("kb07")) +contrasts = Dict(:item => Grouping(), + :subj => Grouping(), + :spkr => EffectsCoding(), + :prec => EffectsCoding(), + :load => EffectsCoding()) + +fm1 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1|item) + (1|subj)), kb07; contrasts, progress=true, thin=1) +plot_blups(args...; kwargs...) = plot_blups!(Figure(; resolution=(1000, 600)), args...;kwargs...) +function plot_blups!(f::Makie.FigureLike, m::MixedModel) + + re = ranefinfo(m) + for (idx, grp) in enumerate(propertynames(re)) + gl = f[1,idx] = GridLayout() + qqcaterpillar!(gl, re[grp]) + Label(gl[end+1, :], string(grp); font=:bold, tellwidth=false) + end + Label(f[0, :], "Conditional Modes", tellwidth=false) + for idx in 2:length(re) + ratio = length(re[idx].cnames) / length(re[1].cnames) + colsize!(f.layout, idx, Auto(ratio)) + end + return f +end + +function eff_dim(m::MixedModel, threshold::Real=0.95) + re = MixedModels.rePCA(m) + return NamedTuple{propertynames(re)}(findfirst(>=(threshold), val) for val in values(re)) +end + +fm_max = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + spkr * prec * load |item) + (1 + spkr * prec * load|subj)), kb07; contrasts, progress=true, thin=1) +fm2 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + prec + load|subj)), kb07; contrasts, progress=true, thin=1) +fm3 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + prec + load|subj)), kb07; contrasts, progress=true, thin=1) +fm4 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + load|subj)), kb07; contrasts, progress=true, thin=1) + +function plot_objective(m::MixedModel) + fitlog = optsum.fitlog +end diff --git a/viz.qmd b/viz.qmd index 59a2a7f..d29be3b 100644 --- a/viz.qmd +++ b/viz.qmd @@ -37,24 +37,37 @@ fm1 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1|item) + ( ```{julia} -let f = Figure() - ax1 = Axis(f[1,1], title="Subject") - caterpillar!(ax1, fm1, :subj) - # caterpillar!(ax2, fm1, :subj) - f -end +plot_blups(args...; kwargs...) = plot_blups!(Figure(; resolution=(1000, 600)), args...;kwargs...) +function plot_blups!(f::Makie.FigureLike, m::MixedModel) - let f = Figure(; resolution=(1000, 600)) - gl = f[1,1] = GridLayout() - re = ranefinfo(m2) - qqcaterpillar!(gl, re[:item]) - Label(gl[end+1, :], "Item"; font=:bold) - gl = f[1,2] = GridLayout() - qqcaterpillar!(gl, re[:subj]) - Label(gl[end+1, :], "Subject"; font=:bold) - Label(f[0, :], "Conditional Modes") - colsize!(f.layout, 1, Auto(0.5)) - save(joinpath(OUTDIR, "qqcat_kb07_joint.png"), f) - f + re = ranefinfo(m) + for (idx, grp) in enumerate(propertynames(re)) + gl = f[1,idx] = GridLayout() + qqcaterpillar!(gl, re[grp]) + Label(gl[end+1, :], string(grp); font=:bold, tellwidth=false) + end + Label(f[0, :], "Conditional Modes", tellwidth=false) + for idx in 2:length(re) + ratio = length(re[idx].cnames) / length(re[1].cnames) + colsize!(f.layout, idx, Auto(ratio)) end + return f +end +``` + +```{julia} +function eff_dim(m::MixedModel, threshold::Real=0.95) + re = MixedModels.rePCA(m) + return NamedTuple{propertynames(re)}(findfirst(>=(threshold), val) for val in values(re)) +end ``` + +```{julia} +fm_max = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + spkr * prec * load |item) + (1 + spkr * prec * load|subj)), kb07; contrasts, progress=true, thin=1) +``` + +fm2 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + prec + load|subj)), kb07; contrasts, progress=true, thin=1) + +fm3 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + prec + load|subj)), kb07; contrasts, progress=true, thin=1) + +fm4 = fit(MixedModel, @formula(1000 / rt_raw ~ spkr * prec * load + (1 + prec + load |item) + (1 + load|subj)), kb07; contrasts, progress=true, thin=1)