specify re by name (#53)

* specify re by name

* Apply suggestions from code review

Co-authored-by: Douglas Bates <[email protected]>

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"

src/MixedModelsSim.jl

@@ -12,6 +12,8 @@ using MixedModels: replicate
 
 export
     create_re,
+    create_theta,
+    createθ,
     cyclicshift,
     factorproduct,
     flatlowertri,

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))
 

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)

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