diff --git a/NEWS.md b/NEWS.md
index 8ce4cb1bd..5aa022cfb 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+MixedModels v5.0.3 Release Notes
+==============================
+- `lowerbd(::MixedModel)` returns the _canonical_ lower bounds of a model's parameters, i.e. the expected bounds after rectification in unconstrained optimization. [#864]
+
 MixedModels v5.0.2 Release Notes
 ==============================
 - The default display and `confint` methods for bootstrap results from models without dispersion parameters has been fixed. [#861]
@@ -695,3 +699,4 @@ Package dependencies
 [#858]: https://github.com/JuliaStats/MixedModels.jl/issues/858
 [#860]: https://github.com/JuliaStats/MixedModels.jl/issues/860
 [#861]: https://github.com/JuliaStats/MixedModels.jl/issues/861
+[#864]: https://github.com/JuliaStats/MixedModels.jl/issues/864
diff --git a/Project.toml b/Project.toml
index 0ed61fc22..ed744c8e8 100644
--- a/Project.toml
+++ b/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"
diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl
index 36f67f97c..60b113e9f 100644
--- a/src/generalizedlinearmixedmodel.jl
+++ b/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
diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl
index 2f5ab3a20..de2b3bcc6 100644
--- a/src/linearmixedmodel.jl
+++ b/src/linearmixedmodel.jl
@@ -787,7 +787,16 @@ 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{T}) where {T} =
+    [(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}[]
diff --git a/src/mixedmodel.jl b/src/mixedmodel.jl
index efd498af1..036e65253 100644
--- a/src/mixedmodel.jl
+++ b/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
 
diff --git a/test/pls.jl b/test/pls.jl
index d5d3af870..d86ceadcf 100644
--- a/test/pls.jl
+++ b/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)