tolerances for issingular (#783)

* tolerances for issingular

* version bump

* NEWS

Author: palday

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

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>", "Jose Bayoan Santiago Calderon <[email protected]>"]
-version = "4.25.4"
+version = "4.26.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

src/bootstrap.jl

@@ -408,15 +408,22 @@ 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.
 
 Equality comparisons are used b/c small non-negative θ values are replaced by 0 in `fit!`.
 
 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)
 

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