recover from PIRLS failure by returning finitial

src/generalizedlinearmixedmodel.jl

@@ -282,7 +282,12 @@ function 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 PosDefException
+            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)])

src/linearmixedmodel.jl

@@ -451,13 +451,7 @@ function fit!(
                 ArgumentError("The response is constant and thus model fitting has failed")
             )
         else
-            # give it one more try with a massive change in scaling
-            @info "Initial step failed, rescaling initial guess and trying again."
-            @warn """Failure of the initial step is often indicative of a model specification
-                     that is not well supported by the data and/or a poorly scaled model.
-                  """
-            optsum.initial ./= maximum(m.sqrtwts)^2 * maximum(response(m))
-            optsum.finitial = obj(optsum.initial, T[])
+            rethrow()
         end
     end
     empty!(fitlog)

test/pirls.jl

@@ -148,7 +148,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)
@@ -159,6 +159,10 @@ 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.θ atol=0.05
+    @test gm4.β[2:end] ≈ gm4slow.β[2:end] atol=0.1
 end
 
 @testset "goldstein" begin # from a 2020-04-22 msg by Ben Goldstein to R-SIG-Mixed-Models