attempt rescaling of poorly scaled problems

src/linearmixedmodel.jl

@@ -428,6 +428,11 @@ function fit!(
     if optsum.feval > 0
         throw(ArgumentError("This model has already been fitted. Use refit!() instead."))
     end
+    if all(==(first(m.y)), m.y)
+        throw(
+            ArgumentError("The response is constant and thus model fitting has failed")
+        )
+    end
     opt = Opt(optsum)
     optsum.REML = REML
     prog = ProgressUnknown("Minimizing"; showspeed=true)
@@ -436,7 +441,13 @@ function fit!(
     fitlog = optsum.fitlog
     function obj(x, g)
         isempty(g) || throw(ArgumentError("g should be empty for this objective"))
-        val = objective(updateL!(setθ!(m, x)))
+        val = try
+            objective(updateL!(setθ!(m, x)))
+        catch ex
+            ex isa PosDefException || rethrow()
+            iter == 0 && rethrow()
+            m.optsum.finitial
+        end
         iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val))
         progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
         iter += 1
@@ -445,14 +456,17 @@ function fit!(
     NLopt.min_objective!(opt, obj)
     try
         optsum.finitial = obj(optsum.initial, T[])
-    catch PosDefException
-        if all(==(first(m.y)), m.y)
-            throw(
-                ArgumentError("The response is constant and thus model fitting has failed")
-            )
-        else
-            rethrow()
-        end
+    catch ex
+        ex isa PosDefException || rethrow()
+        # 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 ./= (isempty(m.sqrtwts) ? 1.0 : maximum(m.sqrtwts)^2) *
+            maximum(response(m))
+        optsum.finitial = obj(optsum.initial, T[])
     end
     empty!(fitlog)
     push!(fitlog, (copy(optsum.initial), optsum.finitial))

test/pls.jl

@@ -605,3 +605,10 @@ end
     # that we're starting to support a non trivial type hierarchy
     @test typeof(re) == Vector{AbstractReMat{Float64}}
 end
+
+@testset "recovery from misscaling" begin
+    fm1 = MixedModels.unfit!(deepcopy(last(models(:insteval))))
+    fm1.optsum.initial .*= 1e6
+    @test_logs (:info, r"Initial step failed") (:warn, r"Failure of the initial step") fit!(fm1; progress=false)
+    @test objective(fm1) ≈ objective(last(models(:insteval)))
+end