Recover from PIRLS failure by returning finitial (#616)

* attempt rescaling of poorly scaled problems

* recover from PIRLS failure by returning finitial

* fix try-catch

* catch DomainError

* NEWS

* version bump

* try a different test

NEWS.md

@@ -1,6 +1,9 @@
+MixedModels v4.17.0 Release Notes
+==============================
 * New kwarg `amalgamate` can be used to disable amalgation of random0effects terms sharing a single grouping variable. Generally, `amalgamate=false` will result in a slower fit, but may improve convergence in some pathological cases. Note that this feature is experimental and changes to it are **not** considered breakin. [#673]
 * More informative error messages when passing a `Distribution` or `Link` type instead of the desired instance. [#698]
 * More informative error message on the intentional decision not to define methods for the coefficient of determination. [#698]
+* **EXPERIMENTAL** Return `finitial` when PIRLS drifts into a portion of the parameter space that yields a (numerically) invalid covariance matrix. This recovery strategy may be removed in a future release. [#616]
 
 MixedModels v4.16.0 Release Notes
 ==============================
@@ -80,7 +83,6 @@ MixedModels v4.6.5 Release Notes
 * Attempt recovery when the initial parameter values lead to an invalid covariance matrix by rescaling [#615]
 * Return `finitial` when the optimizer drifts into a portion of the parameter space that yields a (numerically) invalid covariance matrix [#615]
 
-
 MixedModels v4.6.4 Release Notes
 ========================
 * Support transformed responses in `predict` [#614]
@@ -431,6 +433,7 @@ Package dependencies
 [#608]: https://github.com/JuliaStats/MixedModels.jl/issues/608
 [#614]: https://github.com/JuliaStats/MixedModels.jl/issues/614
 [#615]: https://github.com/JuliaStats/MixedModels.jl/issues/615
+[#616]: https://github.com/JuliaStats/MixedModels.jl/issues/616
 [#628]: https://github.com/JuliaStats/MixedModels.jl/issues/628
 [#634]: https://github.com/JuliaStats/MixedModels.jl/issues/634
 [#637]: https://github.com/JuliaStats/MixedModels.jl/issues/637

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.16.0"
+version = "4.17.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

src/generalizedlinearmixedmodel.jl

@@ -284,7 +284,15 @@ function StatsAPI.fit!(
     fitlog = optsum.fitlog
     function obj(x, g)
         isempty(g) || throw(ArgumentError("g should be empty for this objective"))
-        val = deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ)
+        val = try
+            deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ)
+        catch ex
+            # this allows us to recover from models where e.g. the link isn't
+            # as constraining as it should be
+            ex isa Union{PosDefException,DomainError} || rethrow()
+            iter == 1 && rethrow()
+            m.optsum.finitial
+        end
         iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val))
         verbose && println(round(val; digits=5), " ", x)
         progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])

test/pirls.jl

@@ -162,7 +162,7 @@ end
     center(v::AbstractVector) = v .- (sum(v) / length(v))
     grouseticks = DataFrame(dataset(:grouseticks))
     grouseticks.ch = center(grouseticks.height)
-    gm4 = fit(MixedModel, only(gfms[:grouseticks]), grouseticks, Poisson(); fast=true, progress=false)  # fails in pirls! with fast=false
+    gm4 = fit(MixedModel, only(gfms[:grouseticks]), grouseticks, Poisson(); fast=true, progress=false)
     @test isapprox(deviance(gm4), 851.4046, atol=0.001)
     # these two values are not well defined at the optimum
     #@test isapprox(sum(x -> sum(abs2, x), gm4.u), 196.8695297987013, atol=0.1)
@@ -173,6 +173,11 @@ end
     @test sdest(gm4) === missing
     @test varest(gm4) === missing
     @test gm4.σ === missing
+    gm4slow = fit(MixedModel, only(gfms[:grouseticks]), grouseticks, Poisson(); fast=false, progress=false)
+    # this tolerance isn't great, but then again the optimum isn't well defined
+    # @test gm4.θ ≈ gm4slow.θ rtol=0.05
+    # @test gm4.β[2:end] ≈ gm4slow.β[2:end] atol=0.1
+    @test isapprox(deviance(gm4), deviance(gm4slow); rtol=0.1)
 end
 
 @testset "goldstein" begin # from a 2020-04-22 msg by Ben Goldstein to R-SIG-Mixed-Models