fix lowerbd(::MixedModel) to work as singularity reference check

Author: palday

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <me@phillipalday.com>", "Douglas Bates <dmbates@gmail.com>"]
-version = "5.0.2"
+version = "5.0.3"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

src/generalizedlinearmixedmodel.jl

@@ -523,6 +523,14 @@ function StatsAPI.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T}
     return accum - (mapreduce(u -> sum(abs2, u), +, m.u) + logdet(m)) / 2
 end
 
+"""
+    lowerbd(m::GeneralizedLinearMixedModel)
+
+Return the vector of _canonical_ lower bounds on the parameters, `θ`.
+
+Note that this method does not distinguish between constrained optimization and
+unconstrained optimization with post-fit canonicalization.
+"""
 lowerbd(m::GeneralizedLinearMixedModel) = lowerbd(m.LMM)
 
 # Base.Fix1 doesn't forward kwargs

src/linearmixedmodel.jl

@@ -787,7 +787,15 @@ function StatsAPI.loglikelihood(m::LinearMixedModel)
     return -objective(m) / 2
 end
 
-lowerbd(m::LinearMixedModel) = foldl(vcat, lowerbd(c) for c in m.reterms)
+"""
+    lowerbd(m::LinearMixedModel)
+
+Return the vector of _canonical_ lower bounds on the parameters, `θ`.
+
+Note that this method does not distinguish between constrained optimization and
+unconstrained optimization with post-fit canonicalization.
+"""
+lowerbd(m::LinearMixedModel) = [(pm[2] == pm[3]) ? zero(T) : T(-Inf) for pm in m.parmap]
 
 function mkparmap(reterms::Vector{<:AbstractReMat{T}}) where {T}
     parmap = NTuple{3,Int}[]

src/mixedmodel.jl

@@ -71,7 +71,7 @@ Equality comparisons are used b/c small non-negative θ values are replaced by 0
 function issingular(
     m::MixedModel{T}, θ=m.θ; atol::Real=0, rtol::Real=atol > 0 ? 0 : √eps()
 ) where {T}
-    lb = [(pm[2] == pm[3]) ? zero(T) : T(-Inf) for pm in m.parmap]
+    lb = lowerbd(m)
     return _issingular(lb, θ; atol, rtol)
 end
 

test/pls.jl

@@ -43,6 +43,7 @@ end
     @test length(fm1.A) == 3
     @test size(fm1.reterms) == (1,)
     @test fm1.optsum.initial == ones(1)
+    @test lowerbd(fm1) == zeros(1)
     fm1.θ = ones(1)
     @test fm1.θ == ones(1)
     @test islinear(fm1)