diff --git a/NEWS.md b/NEWS.md
index 4e2304383..8fc431a7d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+MixedModels v4.26.0 Release Notes
+==============================
+- `issingular` now accepts comparison tolerances through the keyword arguments `atol` and `rtol`. [#783]
+
 MixedModels v4.25.4 Release Notes
 ==============================
 - Added additional precompilation for rePCA. [#749]
@@ -559,3 +563,4 @@ Package dependencies
 [#775]: https://github.com/JuliaStats/MixedModels.jl/issues/775
 [#776]: https://github.com/JuliaStats/MixedModels.jl/issues/776
 [#778]: https://github.com/JuliaStats/MixedModels.jl/issues/778
+[#783]: https://github.com/JuliaStats/MixedModels.jl/issues/783
diff --git a/Project.toml b/Project.toml
index bc051ae95..a49022ce7 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>", "Jose Bayoan Santiago Calderon <jbs3hp@virginia.edu>"]
-version = "4.25.4"
+version = "4.26.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index c90dc33a7..e59809bab 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -408,7 +408,8 @@ function Base.getproperty(bsamp::MixedModelFitCollection, s::Symbol)
 end
 
 """
-    issingular(bsamp::MixedModelFitCollection)
+    issingular(bsamp::MixedModelFitCollection;
+               atol::Real=0, rtol::Real=atol>0 ? 0 : √eps))
 
 Test each bootstrap sample for singularity of the corresponding fit.
 
@@ -416,7 +417,13 @@ Equality comparisons are used b/c small non-negative θ values are replaced by 0
 
 See also [`issingular(::MixedModel)`](@ref).
 """
-issingular(bsamp::MixedModelFitCollection) = map(θ -> any(θ .== bsamp.lowerbd), bsamp.θ)
+function issingular(
+    bsamp::MixedModelFitCollection; atol::Real=0, rtol::Real=atol > 0 ? 0 : √eps()
+)
+    return map(bsamp.θ) do θ
+        return _issingular(bsamp.lowerbd, θ; atol, rtol)
+    end
+end
 
 Base.length(x::MixedModelFitCollection) = length(x.fits)
 
diff --git a/src/mixedmodel.jl b/src/mixedmodel.jl
index 3b145f959..36e11fc81 100644
--- a/src/mixedmodel.jl
+++ b/src/mixedmodel.jl
@@ -58,7 +58,7 @@ function StatsAPI.dof_residual(m::MixedModel)
 end
 
 """
-    issingular(m::MixedModel, θ=m.θ)
+    issingular(m::MixedModel, θ=m.θ; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps)
 
 Test whether the model `m` is singular if the parameter vector is `θ`.
 
@@ -68,8 +68,15 @@ Equality comparisons are used b/c small non-negative θ values are replaced by 0
     For `GeneralizedLinearMixedModel`, the entire parameter vector (including
     β in the case `fast=false`) must be specified if the default is not used.
 """
-issingular(m::MixedModel, θ=m.θ) = any(lowerbd(m) .== θ)
-issingular(m::GeneralizedLinearMixedModel, θ=m.optsum.final) = any(lowerbd(m) .== θ)
+function issingular(m::MixedModel, θ=m.θ; atol::Real=0, rtol::Real=atol > 0 ? 0 : √eps())
+    return _issingular(m.lowerbd, θ; atol, rtol)
+end
+
+function _issingular(v, w; atol, rtol)
+    return any(zip(v, w)) do (x, y)
+        return isapprox(x, y; atol, rtol)
+    end
+end
 
 # FIXME: better to base this on m.optsum.returnvalue
 StatsAPI.isfitted(m::MixedModel) = m.optsum.feval > 0