From 4833c11a6b3077418b80384763c34cb758df5a84 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Wed, 9 Aug 2023 15:16:52 +0000 Subject: [PATCH] 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 | 5 ++++- Project.toml | 2 +- src/generalizedlinearmixedmodel.jl | 10 +++++++++- test/pirls.jl | 7 ++++++- 4 files changed, 20 insertions(+), 4 deletions(-) diff --git a/NEWS.md b/NEWS.md index 67cbb1e66..c3019459f 100644 --- a/NEWS.md +++ b/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 diff --git a/Project.toml b/Project.toml index 498deacd3..0aceb5fa4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MixedModels" uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" author = ["Phillip Alday ", "Douglas Bates ", "Jose Bayoan Santiago Calderon "] -version = "4.16.0" +version = "4.17.0" [deps] Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index dacefefc6..568740f4e 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/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)]) diff --git a/test/pirls.jl b/test/pirls.jl index b50f0fa1a..96a005033 100644 --- a/test/pirls.jl +++ b/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