From e2e3d5d9e0163ecaab285576901571701cb504d7 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Wed, 15 Sep 2021 11:03:52 -0500 Subject: [PATCH] specify re by name (#53) * specify re by name * Apply suggestions from code review Co-authored-by: Douglas Bates --- Project.toml | 2 +- src/MixedModelsSim.jl | 2 ++ src/utilities.jl | 79 ++++++++++++++++++++++++++++++++++++++++--- test/power.jl | 4 +-- test/utilities.jl | 34 ++++++++++++++++--- 5 files changed, 109 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index d0cb6c3..ebd0bcb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MixedModelsSim" uuid = "d5ae56c5-23ca-4a1f-b505-9fc4796fc1fe" authors = ["Phillip Alday", "Douglas Bates", "Lisa DeBruine", "Reinhold Kliegl"] -version = "0.2.4" +version = "0.2.5" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/MixedModelsSim.jl b/src/MixedModelsSim.jl index 376b3ec..71e891e 100644 --- a/src/MixedModelsSim.jl +++ b/src/MixedModelsSim.jl @@ -12,6 +12,8 @@ using MixedModels: replicate export create_re, + create_theta, + createθ, cyclicshift, factorproduct, flatlowertri, diff --git a/src/utilities.jl b/src/utilities.jl index 98bdb8c..73b5b84 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -155,7 +155,7 @@ end """ - update!(m::MixedModel; θ) + _update!(m::MixedModel, θ) Update the mixed model to use θ as its new parameter vector. @@ -167,8 +167,8 @@ Update the mixed model to use θ as its new parameter vector. !!! note For GLMMs, this only sets θ and not β, even for `fast=false` fits. """ -update!(m::LinearMixedModel; θ) = updateL!(MixedModels.setθ!(m, θ)) -update!(m::GeneralizedLinearMixedModel; θ) = pirls!(MixedModels.setθ!(m, θ), false) +_update!(m::LinearMixedModel, θ) = updateL!(MixedModels.setθ!(m, θ)) +_update!(m::GeneralizedLinearMixedModel, θ) = pirls!(MixedModels.setθ!(m, θ), false) # arguably type piracy, but we're all the same developers.... @@ -182,16 +182,87 @@ The `re` can be created using [`create_re`](@ref). They should be specified in the order specified in `VarCorr(m)`. +!!! warning + We recommend against calling this method directly. Instead, use the method + with keyword arguments to specify the different `re` by name. + +!!! note + This is a convenience function for installing a particular parameter vector + and the resulting model fit. It does not actually perform any type of + optimization. + Details ======== The `re` used as the λ fields of the model's `ReTerm`s and should be specified as the lower Cholesky factor of covariance matrices. """ function update!(m::MixedModel, re...) + Base.depwarn("Specifying the random effects by position instead of name is deprecated", + :update!vararg) + θ = vcat((flatlowertri(rr) for rr in re)...) - update!(m; θ=θ) + return _update!(m, θ) end +""" + update!(m::MixedModel; namedre...) + update!(m::MixedModel; θ) + +Update the mixed model to use the random-effects covariance matrices. + +The `namedre` can be created using [`create_re`](@ref). The `namedre` are specified +by the name of the blocking variable, e.g. `subj=create_re(...)`. + +!!! warning + Setting θ directly as a keyword-argument is deprecated. + +!!! note + This is a convenience function for installing a particular parameter vector + and the resulting model fit. It does not actually perform any type of + optimization. + +Details +======== +The `re` is used as the λ fields of the model's `ReTerm`s and should be specified +as the lower Cholesky factor of covariance matrices. +""" +function update!(m::MixedModel; θ=nothing, namedre...) + if θ === nothing + length(namedre) == 0 && throw(ArgumentError("no random effects specified")) + θ = createθ(m; namedre...) + elseif θ !== nothing && length(namedre) != 0 + throw(ArgumentError("θ must not be specified when named random effects are")) + end + return _update!(m, θ) +end + +""" + createθ(m::MixedModel; named_re) + create_theta(m::MixedModel; named_re) + +Create the parameter vector θ corresponding to the random effects. + +The `named_re` can be created using [`create_re`](@ref). The `named_re` are specified +by the name of the blocking variable, e.g. `subj=create_re(...)`. + +The model must be specified because the parameters are sorted internally for computational +efficiency. +""" +function createθ(m::MixedModel; named_re...) + newre = Dict(named_re...) + fns = fnames(m) + + if Set(fns) != keys(newre) + throw(ArgumentError("Specified blocking variables do not match model blocking variables.")) + end + + return mapfoldl(vcat, fns) do fname + flatlowertri(newre[fname]) + end +end + +const create_theta = createθ + """ create_re(sigmas...; corrmat=Matrix{Float64}(I, length(sigmas), length(sigmas)) diff --git a/test/power.jl b/test/power.jl index 11b2849..3732931 100644 --- a/test/power.jl +++ b/test/power.jl @@ -15,8 +15,8 @@ dat = simdat_crossed(StableRNG(42), subj_n, item_n, both_win = both_win) form = @formula(dv ~ 1 + age + pet + cond + time + (1|subj) + (1|item)); cont = Dict(nm => HelmertCoding() for nm in (:age, :pet, :cond, :time)); -fm1 = fit(MixedModel, form, dat; contrasts=cont, REML=false); -sim = parametricbootstrap(StableRNG(42), 10, fm1; use_threads=false); +fm1 = fit(MixedModel, form, dat; contrasts=cont, REML=false, progress=false); +sim = parametricbootstrap(StableRNG(42), 10, fm1; use_threads=false, hide_progress=true); @testset "power_table" begin pt = power_table(sim) diff --git a/test/utilities.jl b/test/utilities.jl index 2c5fb66..3c6c9e3 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -60,7 +60,7 @@ end @testset "update!" begin @testset "LMM" begin fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), - MixedModels.dataset(:sleepstudy)) + MixedModels.dataset(:sleepstudy); progress=false) update!(fm1; θ=[1.0, 0.0, 1.0]) @test all(values(first(fm1.σs)) .== fm1.σ) @@ -70,7 +70,7 @@ end @testset "GLMM" begin cbpp = MixedModels.dataset(:cbpp) gm1 = fit(MixedModel, @formula((incid/hsz) ~ 1 + period + (1|herd)), cbpp, Binomial(); - wts=cbpp.hsz, fast=true) + wts=cbpp.hsz, fast=true, progress=false) β = repeat([-1.], 4) θ = [0.5] @@ -80,7 +80,7 @@ end @test all(values(first(gm1.σs)) .== θ) @test all(gm1.β .== β₀) - refit!(gm1, fast=false) + refit!(gm1, fast=false, progress=false) β₀ = copy(fixef(gm1)) update!(gm1; θ=θ) @test all(values(first(gm1.σs)) .== θ) @@ -89,9 +89,9 @@ end end -@testset "create_re" begin +@testset "create_re|θ" begin fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), - MixedModels.dataset(:sleepstudy)) + MixedModels.dataset(:sleepstudy); progress=false) σs = values(first(fm1.σs)) ./ fm1.σ ρ = only(first(VarCorr(fm1).σρ).ρ) @@ -105,10 +105,34 @@ end fm1unfitted = LinearMixedModel(@formula(reaction ~ 1 + days + (1 + days|subj)), MixedModels.dataset(:sleepstudy)) + @test_logs((:warn, + "Specifying the random effects by position instead of name is deprecated"), + update!(fm1unfitted, subjre)) + update!(fm1unfitted, subjre) vcu, vc = VarCorr(fm1unfitted), VarCorr(fm1) @test all(values(vcu.σρ.subj.σ) .≈ values(vc.σρ.subj.σ)) @test all(values(vcu.σρ.subj.ρ) .≈ values(vc.σρ.subj.ρ)) + + @test_throws ArgumentError update!(fm1unfitted; subjitem=subjre, θ=fm1.θ) + @test_throws ArgumentError update!(fm1unfitted) + @test_throws ArgumentError update!(fm1unfitted; item=subjre) + + fmkb = fit(MixedModel, + @formula(rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1 + prec|item)), + MixedModels.dataset(:kb07); progress=false) + + fmkb2 = LinearMixedModel(@formula(rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1 + prec|item)), + MixedModels.dataset(:kb07)) + update!(fmkb2; subj=fmkb.reterms[2].λ, item=fmkb.reterms[1].λ) + + vcu, vc = VarCorr(fmkb2), VarCorr(fmkb2) + + @test all(values(vcu.σρ.subj.σ) .≈ values(vc.σρ.subj.σ)) + @test all(values(vcu.σρ.subj.ρ) .≈ values(vc.σρ.subj.ρ)) + @test all(values(vcu.σρ.item.σ) .≈ values(vc.σρ.item.σ)) + @test all(values(vcu.σρ.item.ρ) .≈ values(vc.σρ.item.ρ)) + end