GLMM models

.gitignore

@@ -0,0 +1,2 @@
+Project.toml
+*.png

Project.toml

@@ -0,0 +1,13 @@
+[deps]
+AlgebraOfGraphics = "cbdf2221-f076-402e-a563-3d30da359d67"
+Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Chain = "8be319e6-bccf-4806-a6f7-6fae938471bc"
+DataFrameMacros = "75880514-38bc-4a95-a458-c2aea5a3a702"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
+MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
+MixedModelsMakie = "b12ae82c-6730-437f-aff9-d2c38332a376"
+StandardizedPredictors = "5064a6a7-f8c2-40e2-8bdc-797ec6f1ae18"

scripts/glmm.jl

@@ -0,0 +1,76 @@
+using AlgebraOfGraphics
+using Arrow
+using BenchmarkTools
+using CairoMakie
+using DataFrameMacros
+using DataFrames
+using MixedModels
+using MixedModelsMakie
+using MKL
+using StandardizedPredictors
+
+items = DataFrame(Arrow.Table("./arrow/ldt_item.arrow"));
+subj = DataFrame(Arrow.Table("./arrow/ldt_subj.arrow"));
+trials = DataFrame(Arrow.Table("./arrow/ldt_trial.arrow"));
+
+dat = leftjoin(
+    leftjoin(
+        @subset(
+            trials,
+            !ismissing(:acc) && (200 < :rt < 4000) # bounds may need adjustment
+        ),
+        select(items, :item, :isword, :wrdlen),
+        on=:item,
+    ),
+    select(subj, :subj, :sex, :vocabAge),
+    on=:subj,
+)
+
+nobs_subj = combine(groupby(dat, :subj), nrow => :n);
+
+extrema(nobs_subj.n)
+
+data(nobs_subj) * mapping(:n => "Number of observations per subject") * AlgebraOfGraphics.density() |> draw
+
+contrasts = Dict(
+    :subj => Grouping(),
+    :item => Grouping(),
+    :isword => HelmertCoding(),
+    :sex => HelmertCoding(),
+    :wrdlen => Center(8),
+    :vocabAge => Center(17),
+);
+
+thin = 1;   # preserve all evaluations in optsum.fitlog
+
+formula = @formula(acc ~ 1 + wrdlen + vocabAge + (1 | item) + (1 | subj))
+
+@time fit(MixedModel, formula, dat, Bernoulli(); fast=true, contrasts, thin)
+
+@time fit(MixedModel, formula, dat, Bernoulli(); fast=true, contrasts, thin, lmminit=[:β, :θ])
+@time fit(MixedModel, formula, dat, Bernoulli(); fast=true, contrasts, thin, lmminit=[:θ])
+@time fit(MixedModel, formula, dat, Bernoulli(); fast=true, contrasts, thin, lmminit=[:β])
+
+
+@time m = fit(MixedModel, formula, dat, Bernoulli(); fast=false, contrasts, thin);
+@time mβθ = fit(MixedModel, formula, dat, Bernoulli(); fast=false, contrasts, thin, lmminit=[:β, :θ]);
+@time mθ = fit(MixedModel, formula, dat, Bernoulli(); fast=false, contrasts, thin, lmminit=[:θ]);
+@time mβ = fit(MixedModel, formula, dat, Bernoulli(); fast=false, contrasts, thin, lmminit=[:β]);
+
+
+df = DataFrame(; parameter=first.(m.optsum.fitlog), objective=last.(m.optsum.fitlog))
+df[!, :init] .= "glm"
+df[!, :iter] = 1:nrow(df)
+
+for (init, model) in ["lmm-βθ" => mβθ, "lmm-θ" => mθ, "lmm-β" => mβ]
+    df2 = DataFrame(; parameter=first.(model.optsum.fitlog), objective=last.(model.optsum.fitlog))
+    df2[!, :init] .= init
+    df2[!, :iter] = 1:nrow(df2)
+    append!(df, df2)
+end
+
+Arrow.write("glmm_fitlog_by_init.arrow", df; compress=:zstd)
+
+data(filter(:iter => <(50), df)) * mapping(:iter, :objective; color=:init) * visual(Lines) |> draw
+
+data(filter(:iter => >(50), df)) * mapping(:iter, :objective; color=:init) * visual(Lines) |> draw