use a named tuple with coefnames for β in simulate (#262)

* use a named tuple with fixefnames for β in simulate

* handle rank-deficiency in parametricbootstrap and simulate!

Co-authored-by: Phillip Alday <[email protected]>

Author: kleinschmidt

src/simulate.jl

@@ -58,7 +58,7 @@ function parametricbootstrap(
     rng::AbstractRNG,
     n::Integer,
     morig::LinearMixedModel{T};
-    β = morig.β,
+    β = coef(morig),
     σ = morig.σ,
     θ = morig.θ,
     use_threads = false,
@@ -68,6 +68,14 @@ function parametricbootstrap(
     θ = convert(Vector{T},θ)
     βsc, θsc, p, k, m = similar(β), similar(θ), length(β), length(θ), deepcopy(morig)
     y₀ = copy(response(m))
+
+    β_names = (Symbol.(fixefnames(morig))..., )
+    rank = length(β_names)
+    # fixef! requires that we take all coefs, even for pivoted terms
+    if rank ≠ length(βsc)
+        resize!(βsc, rank)
+    end
+
     # we need to do for in-place operations to work across threads
     m_threads = [m]
     βsc_threads = [βsc]
@@ -93,7 +101,10 @@ function parametricbootstrap(
         (
          objective = mod.objective,
          σ = mod.σ,
-         β = SVector{p,T}(fixef!(βsc, mod)),
+         # fixef! does the pivoted, but not truncated coefs
+         # coef does the non-pivoted
+         # fixef does either pivoted+truncated or unpivoted
+         β = NamedTuple{β_names}(fixef!(βsc, mod)[1:rank]),
          θ = SVector{k,T}(getθ!(θsc, mod)),
         )
     end
@@ -154,11 +165,22 @@ function simulate!(
     σ = m.σ,
     θ = T[],
 ) where {T}
+    length(β) == length(fixef(m)) ||
+        length(β) == length(coef(m)) ||
+            throw(ArgumentError("You must specify all (non-singular) βs"))
+
     β = convert(Vector{T},β)
     σ = T(σ)
     θ = convert(Vector{T},θ)
     isempty(θ) || setθ!(m, θ)
 
+    if length(β) ≠ length(coef(m))
+        padding = length(coef(m)) - length(β)
+        for ii in 1:padding
+            push!(β, -0.0)
+        end
+    end
+
     y = randn!(rng, response(m))      # initialize y to standard normal
 
     for trm in m.reterms              # add the unscaled random effects
@@ -169,7 +191,7 @@ function simulate!(
     m
 end
 
-function simulate!(m::LinearMixedModel{T}; β = m.β, σ = m.σ, θ = T[]) where {T}
+function simulate!(m::LinearMixedModel{T}; β = coef(m), σ = m.σ, θ = T[]) where {T}
     simulate!(Random.GLOBAL_RNG, m, β = β, σ = σ, θ = θ)
 end