diff --git a/.github/workflows/current.yml b/.github/workflows/current.yml
index 5eaa42d52..d7fcc13a4 100644
--- a/.github/workflows/current.yml
+++ b/.github/workflows/current.yml
@@ -19,12 +19,13 @@ jobs:
   ci:
     runs-on: ${{ matrix.os }}
     strategy:
+      fail-fast: false
       matrix:
         julia-version: [1]
         # julia-arch: [x64]
         os:
-         - ubuntu-22.04
-         - macOS-14 # apple silicon!
+         - ubuntu-24.04
+         - macOS-15 # apple silicon!
     steps:
       - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
diff --git a/NEWS.md b/NEWS.md
index a16bbc05f..28dad2736 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@ MixedModels v5.0.0 Release Notes
 ==============================
 - Optimization is now performed _without constraints_. In a post-fitting step, the Cholesky factor is canonicalized to have non-negative diagonal elements. [#840]
 - The default optimizer has changed to NLopt's implementation of NEWUOA where possible. NLopt's implementation fails on 1-dimensional problems, so in the case of a single, scalar random effect, BOBYQA is used instead. In the future, the default optimizer backend will likely change to PRIMA and NLopt support will be moved to an extension. Blocking this change in backend is an issue with PRIMA.jl when running in VSCode's built-in REPL on Linux. [#840]
-- [BREAKING] Support for constrained optimization has been completely removed, i.e. the field `lowerbd` has been removed from `OptSummary`.
+- [BREAKING] Support for constrained optimization has been completely removed, i.e. the field `lowerbd` has been removed from `OptSummary`. [#849]
 - [BREAKING] A fitlog is always kept -- the deprecated keyword argument `thin` has been removed as has the `fitlog` keyword argument. [#850]
 - The fitlog is now stored as Tables.jl-compatible column table. [#850]
 - Internal code around the default optimizer has been restructured. In particular, the NLopt backend has been moved to a submodule, which will make it easier to move it to an extension if we promote another backend to the default. [#853]
diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl
index d8383da05..3131267c5 100644
--- a/ext/MixedModelsPRIMAExt.jl
+++ b/ext/MixedModelsPRIMAExt.jl
@@ -105,7 +105,7 @@ function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend;
         sc
     end
     info = _optimizer!(Val(optsum.optimizer), obj, optsum.final;
-        xl=optsum.lowerbd, maxfun,
+        maxfun,
         optsum.rhoend, optsum.rhobeg,
         scale)
     ProgressMeter.finish!(prog)
diff --git a/src/MixedModelsNLoptExt.jl b/src/MixedModelsNLoptExt.jl
index 6452e4fe5..a82527e46 100644
--- a/src/MixedModelsNLoptExt.jl
+++ b/src/MixedModelsNLoptExt.jl
@@ -97,8 +97,7 @@ function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend;
 end
 
 function NLopt.Opt(optsum::OptSummary)
-    lb = optsum.lowerbd
-    n = length(lb)
+    n = length(optsum.initial)
 
     if optsum.optimizer == :LN_NEWUOA && isone(n) # :LN_NEWUOA doesn't allow n == 1
         optsum.optimizer = :LN_BOBYQA
@@ -110,11 +109,10 @@ function NLopt.Opt(optsum::OptSummary)
     if length(optsum.xtol_abs) == n  # not true for fast=false optimization in GLMM
         NLopt.xtol_abs!(opt, optsum.xtol_abs) # absolute criterion on parameter values
     end
-    #    NLopt.lower_bounds!(opt, lb)   # use unconstrained optimization even for :LN_BOBYQA
     NLopt.maxeval!(opt, optsum.maxfeval)
     NLopt.maxtime!(opt, optsum.maxtime)
     if isempty(optsum.initial_step)
-        optsum.initial_step = NLopt.initial_step(opt, optsum.initial, similar(lb))
+        optsum.initial_step = NLopt.initial_step(opt, optsum.initial)
     else
         NLopt.initial_step!(opt, optsum.initial_step)
     end
diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index 330f995b3..e9e0354f4 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -12,7 +12,7 @@ Object returned by `parametericbootstrap` with fields
 - `fits`: the parameter estimates from the bootstrap replicates as a vector of named tuples.
 - `λ`: `Vector{LowerTriangular{T,Matrix{T}}}` containing copies of the λ field from `ReMat` model terms
 - `inds`: `Vector{Vector{Int}}` containing copies of the `inds` field from `ReMat` model terms
-- `lowerbd`: `Vector{T}` containing the vector of lower bounds (corresponds to the identically named field of [`OptSummary`](@ref))
+- `lowerbd`: `Vector{T}` containing the vector of lower bounds (corresponds to the `lowerbd(model)` of the original model)
 - `fcnames`: NamedTuple whose keys are the grouping factor names and whose values are the column names
 
 The schema of `fits` is, by default,
@@ -34,7 +34,7 @@ struct MixedModelBootstrap{T<:AbstractFloat} <: MixedModelFitCollection{T}
     fits::Vector
     λ::Vector{Union{LowerTriangular{T},Diagonal{T}}}
     inds::Vector{Vector{Int}}
-    lowerbd::Vector{T}
+    lowerbd::Vector{T} # we need to store this explicitly because we no longer have access to the ReMats
     fcnames::NamedTuple
 end
 
@@ -112,7 +112,7 @@ function restorereplicates(
         samp,
         map(vv -> T.(vv), m.λ), # also does a deepcopy if no type conversion is necessary
         getfield.(m.reterms, :inds),
-        T.(m.optsum.lowerbd[1:length(first(samp).θ)]),
+        T.(lowerbd(m)),
         NamedTuple{Symbol.(fnames(m))}(map(t -> Tuple(t.cnames), m.reterms)),
     )
 end
@@ -175,6 +175,8 @@ function Base.show(io::IO, mime::MIME"text/plain", x::MixedModelBootstrap)
     return nothing
 end
 
+lowerbd(x::MixedModelFitCollection) = x.lowerbd
+
 """
     parametricbootstrap([rng::AbstractRNG], nsamp::Integer, m::MixedModel{T}, ftype=T;
         β = fixef(m), σ = m.σ, θ = m.θ, progress=true, optsum_overrides=(;))
@@ -243,8 +245,6 @@ function parametricbootstrap(
     for (key, val) in pairs(optsum_overrides)
         setfield!(m.optsum, key, val)
     end
-    # this seemed to slow things down?!
-    # _copy_away_from_lowerbd!(m.optsum.initial, morig.optsum.final, m.lowerbd; incr=0.05)
 
     β_names = Tuple(Symbol.(coefnames(morig)))
 
@@ -267,7 +267,7 @@ function parametricbootstrap(
         samp,
         map(vv -> ftype.(vv), morig.λ), # also does a deepcopy if no type conversion is necessary
         getfield.(morig.reterms, :inds),
-        ftype.(morig.optsum.lowerbd[1:length(first(samp).θ)]),
+        ftype.(lowerbd(morig)),
         NamedTuple{Symbol.(fnames(morig))}(map(t -> Tuple(t.cnames), morig.reterms)),
     )
 end
@@ -421,7 +421,7 @@ function issingular(
     bsamp::MixedModelFitCollection; atol::Real=0, rtol::Real=atol > 0 ? 0 : √eps()
 )
     return map(bsamp.θ) do θ
-        return _issingular(bsamp.lowerbd, θ; atol, rtol)
+        return _issingular(lowerbd(bsamp), θ; atol, rtol)
     end
 end
 
diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl
index 3b41619ee..36f67f97c 100644
--- a/src/generalizedlinearmixedmodel.jl
+++ b/src/generalizedlinearmixedmodel.jl
@@ -30,7 +30,6 @@ In addition to the fieldnames, the following names are also accessible through t
 - `theta`: synonym for `θ`
 - `beta`: synonym for `β`
 - `σ` or `sigma`: common scale parameter (value is `NaN` for distributions without a scale parameter)
-- `lowerbd`: vector of lower bounds on the combined elements of `β` and `θ`
 - `formula`, `trms`, `A`, `L`, and `optsum`: fields of the `LMM` field
 - `X`: fixed-effects model matrix
 - `y`: response vector
@@ -280,7 +279,6 @@ function StatsAPI.fit!(
     end
 
     if !fast
-        optsum.lowerbd = vcat(fill!(similar(β), T(-Inf)), optsum.lowerbd)
         optsum.initial = vcat(β, lm.optsum.final)
         optsum.final = copy(optsum.initial)
     end
@@ -296,8 +294,9 @@ function StatsAPI.fit!(
 
     ## check if very small parameter values bounded below by zero can be set to zero
     xmin_ = copy(xmin)
+    lb = fast ? lowerbd(m) : vcat(zero(β), lowerbd(m))
     for i in eachindex(xmin_)
-        if iszero(optsum.lowerbd[i]) && zero(T) < xmin_[i] < optsum.xtol_zero_abs
+        if iszero(lb[i]) && zero(T) < xmin_[i] < optsum.xtol_zero_abs
             xmin_[i] = zero(T)
         end
     end
@@ -524,6 +523,8 @@ function StatsAPI.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T}
     return accum - (mapreduce(u -> sum(abs2, u), +, m.u) + logdet(m)) / 2
 end
 
+lowerbd(m::GeneralizedLinearMixedModel) = lowerbd(m.LMM)
+
 # Base.Fix1 doesn't forward kwargs
 function objective!(m::GeneralizedLinearMixedModel; fast=false, kwargs...)
     return x -> _objective!(m, x, Val(fast); kwargs...)
@@ -560,7 +561,6 @@ function Base.propertynames(m::GeneralizedLinearMixedModel, private::Bool=false)
         :sigma,
         :X,
         :y,
-        :lowerbd,
         :objective,
         :σρs,
         :σs,
@@ -785,11 +785,6 @@ function unfit!(model::GeneralizedLinearMixedModel{T}) where {T}
 
     reterms = model.LMM.reterms
     optsum = model.LMM.optsum
-    # we need to reset optsum so that it
-    # plays nice with the modifications fit!() does
-    optsum.lowerbd = mapfoldl(lowerbd, vcat, reterms)      # probably don't need this anymore - now trivial with all elements = -Inf
-    # for variances (bounded at zero), we have ones, while
-    # for everything else (bounded at -Inf), we have zeros
     optsum.initial = map(x -> T(x[2] == x[3]), model.LMM.parmap)
     optsum.final = copy(optsum.initial)
     optsum.xtol_abs = fill!(copy(optsum.initial), 1.0e-10)
@@ -839,7 +834,7 @@ function StatsAPI.weights(m::GeneralizedLinearMixedModel{T}) where {T}
 end
 
 # delegate GLMM method to LMM field
-for f in (:feL, :fetrm, :fixefnames, :(LinearAlgebra.logdet), :lowerbd, :PCA, :rePCA)
+for f in (:feL, :fetrm, :fixefnames, :(LinearAlgebra.logdet), :PCA, :rePCA)
     @eval begin
         $f(m::GeneralizedLinearMixedModel) = $f(m.LMM)
     end
diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl
index ce597fbcf..2f5ab3a20 100644
--- a/src/linearmixedmodel.jl
+++ b/src/linearmixedmodel.jl
@@ -24,7 +24,6 @@ Linear mixed-effects model representation
 * `σ` or `sigma`: current value of the standard deviation of the per-observation noise
 * `b`: random effects on the original scale, as a vector of matrices
 * `u`: random effects on the orthogonal scale, as a vector of matrices
-* `lowerbd`: lower bounds on the elements of θ
 * `X`: the fixed-effects model matrix
 * `y`: the response vector
 """
@@ -173,9 +172,8 @@ function LinearMixedModel(
     reweight!.(reterms, Ref(sqrtwts))
     reweight!(Xy, sqrtwts)
     A, L = createAL(reterms, Xy)
-    lbd = foldl(vcat, lowerbd(c) for c in reterms)
     θ = foldl(vcat, getθ(c) for c in reterms)
-    optsum = OptSummary(θ, lbd)
+    optsum = OptSummary(θ)
     optsum.sigma = isnothing(σ) ? nothing : T(σ)
     return LinearMixedModel(
         form,
@@ -662,8 +660,6 @@ function Base.getproperty(m::LinearMixedModel{T}, s::Symbol) where {T}
         stderror(m)
     elseif s == :u
         ranef(m; uscale=true)
-    elseif s == :lowerbd
-        m.optsum.lowerbd
     elseif s == :X
         modelmatrix(m)
     elseif s == :y
@@ -791,7 +787,7 @@ function StatsAPI.loglikelihood(m::LinearMixedModel)
     return -objective(m) / 2
 end
 
-lowerbd(m::LinearMixedModel) = m.optsum.lowerbd
+lowerbd(m::LinearMixedModel) = foldl(vcat, lowerbd(c) for c in m.reterms)
 
 function mkparmap(reterms::Vector{<:AbstractReMat{T}}) where {T}
     parmap = NTuple{3,Int}[]
@@ -856,7 +852,6 @@ function Base.propertynames(m::LinearMixedModel, private::Bool=false)
         :sigmarhos,
         :b,
         :u,
-        :lowerbd,
         :X,
         :y,
         :corr,
diff --git a/src/mimeshow.jl b/src/mimeshow.jl
index e5ef705a2..88bc6b451 100644
--- a/src/mimeshow.jl
+++ b/src/mimeshow.jl
@@ -184,7 +184,7 @@ end
 function _markdown(s::OptSummary)
     optimizer_settings = [["Optimizer", "`$(s.optimizer)`"],
         ["Backend", "`$(s.backend)`"],
-        ["Lower bounds", string(s.lowerbd)]]
+    ]
 
     for param in opt_params(Val(s.backend))
         push!(optimizer_settings, [string(param), string(getfield(s, param))])
diff --git a/src/optsummary.jl b/src/optsummary.jl
index 4bf284aa4..929f3b503 100644
--- a/src/optsummary.jl
+++ b/src/optsummary.jl
@@ -8,7 +8,6 @@ Summary of an optimization
 ## Tolerances, initial and final values
 * `initial`: a copy of the initial parameter values in the optimization
 * `finitial`: the initial value of the objective
-* `lowerbd`: lower bounds on the parameter values
 * `final`: a copy of the final parameter values from the optimization
 * `fmin`: the final value of the objective
 * `feval`: the number of function evaluations
@@ -65,7 +64,6 @@ Similarly, the list of applicable optimization parameters can be inspected with
 Base.@kwdef mutable struct OptSummary{T<:AbstractFloat}
     initial::Vector{T}
     finitial::T = Inf * one(eltype(initial))
-    lowerbd::Vector{T}
     final::Vector{T} = copy(initial)
     fmin::T = Inf * one(eltype(initial))
     feval::Int = -1
@@ -100,11 +98,9 @@ end
 
 function OptSummary(
     initial::Vector{T},
-    lowerbd::Vector{S},
     optimizer::Symbol=:LN_NEWUOA; kwargs...,
-) where {T<:AbstractFloat,S<:AbstractFloat}
-    TS = promote_type(T, S)
-    return OptSummary{TS}(; initial, lowerbd, optimizer, kwargs...)
+) where {T<:AbstractFloat}
+    return OptSummary{T}(; initial, optimizer, kwargs...)
 end
 
 """
@@ -141,7 +137,6 @@ function Base.show(io::IO, ::MIME"text/plain", s::OptSummary)
     println(io)
     println(io, "Backend:                  ", s.backend)
     println(io, "Optimizer:                ", s.optimizer)
-    println(io, "Lower bounds:             ", s.lowerbd)
 
     for param in opt_params(Val(s.backend))
         println(io, rpad(string(param, ":"), length("Initial parameter vector: ")),
diff --git a/src/profile/fixefpr.jl b/src/profile/fixefpr.jl
index 98ab98166..2e2c35c9d 100644
--- a/src/profile/fixefpr.jl
+++ b/src/profile/fixefpr.jl
@@ -44,9 +44,7 @@ function FeProfile(m::LinearMixedModel, tc::TableColumns, j::Integer)
         LinearMixedModel(y₀ - xⱼ * m.β[j], feterm, reterms, m.formula); progress=false
     )
     # not sure this next call makes sense - should the second argument be m.optsum.final?
-    _copy_away_from_lowerbd!(
-        mnew.optsum.initial, mnew.optsum.final, mnew.lowerbd; incr=0.05
-    )
+    copyto!(mnew.optsum.initial, mnew.optsum.final)
     return FeProfile(mnew, tc, y₀, xⱼ, j)
 end
 
diff --git a/src/profile/sigmapr.jl b/src/profile/sigmapr.jl
index 2fbf9cf22..3336ebe52 100644
--- a/src/profile/sigmapr.jl
+++ b/src/profile/sigmapr.jl
@@ -43,7 +43,7 @@ function profileσ(m::LinearMixedModel{T}, tc::TableColumns{T}; threshold=4) whe
         throw(ArgumentError("Can't profile σ, which is fixed at $(optsum.sigma)"))
     θ = copy(optsum.final)
     θinitial = copy(optsum.initial)
-    _copy_away_from_lowerbd!(optsum.initial, optsum.final, optsum.lowerbd)
+    copyto!(optsum.initial, optsum.final)
     obj = optsum.fmin
     σ = m.σ
     pnm = (p=:σ,)
diff --git a/src/profile/thetapr.jl b/src/profile/thetapr.jl
index a89902504..6f60dd7f3 100644
--- a/src/profile/thetapr.jl
+++ b/src/profile/thetapr.jl
@@ -9,7 +9,6 @@ Return an `OptSummary` with the `j`'th component of the parameter omitted.
 function optsumj(os::OptSummary, j::Integer)
     return OptSummary(
         deleteat!(copy(os.final), j),
-        deleteat!(copy(os.lowerbd), j),
         os.optimizer;
         os.backend,
     )
diff --git a/src/profile/utilities.jl b/src/profile/utilities.jl
index f5a6963ac..f973c984d 100644
--- a/src/profile/utilities.jl
+++ b/src/profile/utilities.jl
@@ -145,35 +145,6 @@ function ρvals!(
     return v
 end
 
-"""
-    _copy_away_from_lowerbd!(sink, source, bd; incr=0.01)
-
-Replace `sink[i]` by `max(source[i], bd[i] + incr)`.  When `bd[i] == -Inf` this simply copies `source[i]`.
-"""
-function _copy_away_from_lowerbd!(sink, source, bd; incr=0.01)
-    for i in eachindex(sink, source, bd)
-        @inbounds sink[i] = max(source[i], bd[i] + incr)
-    end
-    return sink
-end
-
-#=  # It appears that this method is not used
-"""
-    stepsize(tbl::Vector{NamedTuple}, resp::Symbol, pred::Symbol; rev::Bool=false)
-
-Return the stepsize from the last value of `tbl.pred` to increase `resp` by approximately 0.5
-"""
-function stepsize(tbl::Vector{<:NamedTuple}, resp::Symbol, pred::Symbol)
-    ntbl = length(tbl)
-    lm1tbl = tbl[ntbl - 1]
-    x1 = getproperty(lm1tbl, pred)
-    y1 = getproperty(lm1tbl, resp)
-    x2 = getproperty(last(tbl), pred)
-    y2 = getproperty(last(tbl), resp)
-    return (x2 - x1) / (2 * (y2 - y1))
-end
-=#
-
 function isnondecreasing(spl::SplineInterpolation)
     return all(≥(0), (Derivative(1) * spl).(spl.x))
 end
diff --git a/src/profile/vcpr.jl b/src/profile/vcpr.jl
index 087996ecb..fca28ef55 100644
--- a/src/profile/vcpr.jl
+++ b/src/profile/vcpr.jl
@@ -36,10 +36,11 @@ function profileσs!(val::NamedTuple, tc::TableColumns{T}; nzlb=1.0e-8) where {T
     m = val.m
     (; λ, σ, β, optsum, parmap, reterms) = m
     isnothing(optsum.sigma) || throw(ArgumentError("Can't profile vc's when σ is fixed"))
-    (; initial, final, fmin, lowerbd) = optsum
-    lowerbd .+= T(nzlb)                       # lower bounds must be > 0 b/c θ's occur in denominators
+    (; initial, final, fmin) = optsum
+    # lowerbd .+= T(nzlb)                       # lower bounds must be > 0 b/c θ's occur in denominators
     saveinitial = copy(initial)
-    copyto!(initial, max.(final, lowerbd))
+    # copyto!(initial, max.(final, lowerbd))
+    copyto!(initial, final)
     zetazero = mkrow!(tc, m, zero(T))         # parameter estimates
     vcnms = filter(keys(first(val.tbl))) do sym
         str = string(sym)
@@ -84,7 +85,6 @@ function profileσs!(val::NamedTuple, tc::TableColumns{T}; nzlb=1.0e-8) where {T
     end
     copyto!(final, initial)
     copyto!(initial, saveinitial)
-    lowerbd .-= T(nzlb)
     optsum.sigma = nothing
     updateL!(setθ!(m, final))
     return val
diff --git a/src/serialization.jl b/src/serialization.jl
index 12884d617..217edd59f 100644
--- a/src/serialization.jl
+++ b/src/serialization.jl
@@ -44,17 +44,11 @@ function restoreoptsum!(
 
     theta_beta_len = length(m.θ) + length(m.β)
     if length(dict.initial) == theta_beta_len # fast=false
-        if length(ops.lowerbd) == length(m.θ)
-            prepend!(ops.lowerbd, fill(-Inf, length(m.β)))
-        end
         setpar! = setβθ!
         varyβ = false
     else # fast=true
         setpar! = setθ!
         varyβ = true
-        if length(ops.lowerbd) != length(m.θ)
-            deleteat!(ops.lowerbd, 1:length(m.β))
-        end
     end
     restoreoptsum!(ops, dict)
     for (par, obj_at_par) in (:initial => :finitial, :final => :fmin)
@@ -75,7 +69,7 @@ end
 function restoreoptsum!(ops::OptSummary{T}, dict::AbstractDict) where {T}
     warn_old_version = true
     allowed_missing = (
-        :lowerbd,       # never saved, -Inf not allowed in JSON
+        :lowerbd,       # never saved, -Inf not allowed in JSON, not used in v5
         :xtol_zero_abs, # added in v4.25.0
         :ftol_zero_abs, # added in v4.25.0
         :sigma,         # added in v4.1.0
@@ -98,10 +92,6 @@ function restoreoptsum!(ops::OptSummary{T}, dict::AbstractDict) where {T}
         warn_old_version = false
     end
 
-    if any(ops.lowerbd .> dict.initial) || any(ops.lowerbd .> dict.final)
-        @debug "" ops.lowerbd dict.initial dict.final
-        throw(ArgumentError("initial or final parameters in io do not satisfy lowerbd"))
-    end
     for fld in (:feval, :finitial, :fmin, :ftol_rel, :ftol_abs, :maxfeval, :nAGQ, :REML)
         setproperty!(ops, fld, getproperty(dict, fld))
     end
@@ -153,16 +143,12 @@ function _deserialize_fitlog(
 end
 
 StructTypes.StructType(::Type{<:OptSummary}) = StructTypes.Mutable()
-StructTypes.excludes(::Type{<:OptSummary}) = (:lowerbd,)
 
 """
     saveoptsum(io::IO, m::MixedModel)
     saveoptsum(filename, m::MixedModel)
 
-Save `m.optsum` (w/o the `lowerbd` field) in JSON format to an IO stream or a file
-
-The reason for omitting the `lowerbd` field is because it often contains `-Inf`
-values that are not allowed in JSON.
+Save `m.optsum` in JSON format to an IO stream or a file
 """
 saveoptsum(io::IO, m::MixedModel) = JSON3.write(io, m.optsum)
 function saveoptsum(filename, m::MixedModel)
diff --git a/test/FactorReTerm.jl b/test/FactorReTerm.jl
index 516042e3f..d4d529865 100644
--- a/test/FactorReTerm.jl
+++ b/test/FactorReTerm.jl
@@ -47,7 +47,6 @@ const LMM = LinearMixedModel
         @test MixedModels.nθ(sf) == 1
         @test MixedModels.getθ(sf) == ones(1)
         @test MixedModels.getθ!(Vector{Float64}(undef, 1), sf) == ones(1)
-        @test lowerbd(sf) == [-Inf]
         @test MixedModels.getθ(setθ!(sf, [0.5])) == [0.5]
         MixedModels.unscaledre!(Vector{Float64}(undef, 30), sf)
         @test_throws DimensionMismatch MixedModels.getθ!(Float64[], sf)
diff --git a/test/mime.jl b/test/mime.jl
index 0f380a340..7f740dbe9 100644
--- a/test/mime.jl
+++ b/test/mime.jl
@@ -134,7 +134,6 @@ lrt = likelihoodratiotest(fm0, fm1)
         | **Optimizer settings**   |                             |
         | Optimizer                | `LN_NEWUOA`                 |
         | Backend                  | `nlopt`                     |
-        | Lower bounds             | [-Inf, -Inf, -Inf]          |
         | ftol_rel                 | 1.0e-12                     |
         | ftol_abs                 | 1.0e-8                      |
         | xtol_rel                 | 0.0                         |
diff --git a/test/pirls.jl b/test/pirls.jl
index aeb16f834..274dce8d1 100644
--- a/test/pirls.jl
+++ b/test/pirls.jl
@@ -62,7 +62,6 @@ end
     # but not when run via Pkg.test(). I have no idea why.
     @test fitlog.θ[end] ≈ gm0.optsum.final
     @test fitlog.objective[end] ≈ gm0.optsum.fmin
-    @test gm0.lowerbd == [-Inf]
     @test isapprox(gm0.θ, [0.5720746212924732], atol=0.001)
     @test !issingular(gm0)
     @test issingular(gm0, [0])
@@ -87,7 +86,6 @@ end
 
     gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); progress=false)
     @test isapprox(gm1.θ, [0.5730523416716424], atol=0.005)
-    @test lowerbd(gm1) == fill(-Inf, 8)
     @test isapprox(deviance(gm1), 2361.545768866505, rtol=0.00001)
     @test isapprox(loglikelihood(gm1), -1180.772884433253, rtol=0.00001)
 
@@ -187,7 +185,6 @@ end
         MixedModel, only(gfms[:verbagg]), dataset(:verbagg), Bernoulli(); progress=false
     )
     @test deviance(gm3) ≈ 8151.40 rtol = 1e-5
-    @test lowerbd(gm3) == fill(-Inf, 8)
     @test fitted(gm3) == predict(gm3)
     # these two values are not well defined at the optimum
     @test isapprox(sum(x -> sum(abs2, x), gm3.u), 273.29646346940785, rtol=1e-3)
diff --git a/test/pls.jl b/test/pls.jl
index 03ab0faa4..f5a4ba87e 100644
--- a/test/pls.jl
+++ b/test/pls.jl
@@ -42,8 +42,6 @@ end
 
     @test length(fm1.A) == 3
     @test size(fm1.reterms) == (1,)
-    @test lowerbd(fm1) == [-Inf]
-    @test fm1.lowerbd == [-Inf]
     @test fm1.optsum.initial == ones(1)
     fm1.θ = ones(1)
     @test fm1.θ == ones(1)
@@ -171,7 +169,6 @@ end
 
 @testset "Dyestuff2" begin
     fm = only(models(:dyestuff2))
-    @test lowerbd(fm) == [-Inf]
     show(IOBuffer(), fm)
     @test fm.θ ≈ zeros(1)
     @test objective(fm) ≈ 162.87303665382575
@@ -198,7 +195,6 @@ end
     fm = only(models(:penicillin))
     @test size(fm) == (144, 1, 30, 2)
     @test fm.optsum.initial == ones(2)
-    @test lowerbd(fm) == fill(-Inf, 2)
 
     @test objective(fm) ≈ 332.1883486700085 atol = 0.001
     @test coef(fm) ≈ [22.97222222222222] atol = 0.001
@@ -241,7 +237,6 @@ end
     fm = last(models(:pastes))
     @test size(fm) == (60, 1, 40, 2)
     @test fm.optsum.initial == ones(2)
-    @test lowerbd(fm) == fill(-Inf, 2)
 
     @test objective(fm) ≈ 247.9944658624955 atol = 0.001
     @test coef(fm) ≈ [60.0533333333333] atol = 0.001
@@ -285,7 +280,6 @@ end
     fm1 = models(:insteval)[2]              # at one time this was the fist of the :insteval models
     @test size(fm1) == (73421, 2, 4114, 3)
     @test fm1.optsum.initial == ones(3)
-    @test lowerbd(fm1) == fill(-Inf, 3)
 
     spL = sparseL(fm1)
     @test size(spL) == (4114, 4114)
@@ -326,7 +320,6 @@ end
 
 @testset "sleep" begin
     fm = last(models(:sleepstudy))
-    @test all(==(-Inf), lowerbd(fm))
     A11 = first(fm.A)
     @test isa(A11, UniformBlockDiagonal{Float64})
     @test isa(first(fm.L), UniformBlockDiagonal{Float64})
@@ -417,7 +410,6 @@ end
     fmnc = models(:sleepstudy)[2]
     @test size(fmnc) == (180, 2, 36, 1)
     @test fmnc.optsum.initial == ones(2)
-    @test lowerbd(fmnc) == fill(-Inf, 2)
     sigmas = fmnc.σs
     @test length(only(sigmas)) == 2
     @test first(only(sigmas)) ≈ 24.171361283849798 atol = 1e-4
diff --git a/test/prima.jl b/test/prima.jl
index dcd0ab198..605ebb713 100644
--- a/test/prima.jl
+++ b/test/prima.jl
@@ -69,7 +69,6 @@ end
 
     Backend:                  prima
     Optimizer:                bobyqa
-    Lower bounds:             [-Inf]
     rhobeg:                   1.0
     rhoend:                   1.0e-6
     maxfeval:                 -1
@@ -94,7 +93,6 @@ end
     | **Optimizer settings**   |                   |
     | Optimizer                | `bobyqa`          |
     | Backend                  | `prima`           |
-    | Lower bounds             | [-Inf]            |
     | rhobeg                   | 1.0               |
     | rhoend                   | 1.0e-6            |
     | maxfeval                 | -1                |