Better handling of constant response (#578)

* Better handling of constant response

For GLMM, this means allow constructing a model with a constant response.
This is useful for simulation.
For both GLMM and LMM, catch an initial PosDefException.
If the PosDefException is from a constant response, throw ArgumentError.

* formatter

* remove extra blank line

* broaden codecoverage

* always do a constant response check for GLMM

.github/workflows/Tier1.yml

@@ -32,4 +32,3 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
       - uses: julia-actions/[email protected]
-        if:  ${{ startsWith(matrix.os, 'Ubuntu') && startsWith(matrix.julia-version, '1.6') }}

NEWS.md

@@ -1,3 +1,8 @@
+MixedModels v4.5.0 Release Notes
+========================
+* Allow constructing a `GeneralizedLinearMixedModel` with constant response, but don't update the ``L`` matrix nor initialize its deviance. This allows for the model to still be used for simulation where the response will be changed before fitting. [#578]
+* Catch `PosDefException` during the first optimization step and throw a more informative `ArgumentError` if the response is constant. [#578]
+
 MixedModels v4.4.1 Release Notes
 ========================
 * Fix type parameterization in MixedModelsBootstrap to support models with a mixture of correlation structures (i.e. `zerocorr` in some but not all RE terms) [#577]
@@ -311,3 +316,4 @@ Package dependencies
 [#572]: https://github.com/JuliaStats/MixedModels.jl/issues/572
 [#573]: https://github.com/JuliaStats/MixedModels.jl/issues/573
 [#577]: https://github.com/JuliaStats/MixedModels.jl/issues/577
+[#578]: https://github.com/JuliaStats/MixedModels.jl/issues/578

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

src/generalizedlinearmixedmodel.jl

@@ -288,6 +288,10 @@ function fit!(
         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
+
     if !fast
         optsum.lowerbd = vcat(fill!(similar(β), T(-Inf)), optsum.lowerbd)
         optsum.initial = vcat(β, m.θ)
@@ -396,6 +400,7 @@ function GeneralizedLinearMixedModel(
 
     LMM = LinearMixedModel(f, tbl; contrasts=contrasts, wts=wts)
     y = copy(LMM.y)
+    constresponse = all(==(first(y)), y)
     # the sqrtwts field must be the correct length and type but we don't know those
     # until after the model is constructed if wt is empty.  Because a LinearMixedModel
     # type is immutable, another one must be created.
@@ -413,7 +418,11 @@ function GeneralizedLinearMixedModel(
             LMM.optsum,
         )
     end
-    updateL!(LMM)
+    # if the response is constant, there's no point (and this may even fail)
+    # we allow this instead of simply failing so that a constant response can
+    # be used as the starting point to simulation where the response will be
+    # overwritten before fitting
+    constresponse || updateL!(LMM)
     # fit a glm to the fixed-effects only
     T = eltype(LMM.Xymat)
     gl = glm(LMM.X, y, d, l; wts=convert(Vector{T}, wts), offset=convert(Vector{T}, offset))
@@ -439,7 +448,10 @@ function GeneralizedLinearMixedModel(
         similar(vv),
         similar(vv),
     )
-    deviance!(res, 1)
+
+    # if the response is constant, there's no point (and this may even fail)
+    constresponse || deviance!(res, 1)
+
     return res
 end
 

src/linearmixedmodel.jl

@@ -443,7 +443,17 @@ function fit!(
         return val
     end
     NLopt.min_objective!(opt, obj)
-    optsum.finitial = obj(optsum.initial, T[])
+    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
+    end
     empty!(fitlog)
     push!(fitlog, (copy(optsum.initial), optsum.finitial))
     fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial))

test/pirls.jl

@@ -122,6 +122,16 @@ end
         @test gm2r.β ≈ -gm2.β atol=1e-3
         @test gm2r.θ ≈ gm2.θ atol=1e-3
     end
+
+    @testset "constant response" begin
+        cbconst = DataFrame(cbpp)
+        cbconst.incid = zero(cbconst.incid)
+        # we do construction and fitting in two separate steps to make sure
+        # that construction succeeds and that the ArgumentError occurs in fitting.
+        mcbconst = GeneralizedLinearMixedModel(first(gfms[:cbpp]), cbconst, Binomial(); wts=float(cbpp.hsz))
+        @test mcbconst isa GeneralizedLinearMixedModel
+        @test_throws ArgumentError("The response is constant and thus model fitting has failed") fit!(mcbconst; progress=false)
+    end
 end
 
 @testset "verbagg" begin

test/pls.jl

@@ -375,6 +375,12 @@ end
     @test isa(b3tbl, NamedTuple)
     @test Tables.istable(only(b3tbl))
 
+    @testset "PosDefException from constant response" begin
+        slp = MixedModels.dataset(:sleepstudy)
+        @test_throws ArgumentError("The response is constant and thus model fitting has failed") refit!(fm, zero(slp.reaction); progress=false)
+        refit!(fm, slp.reaction; progress=false)
+    end
+
     simulate!(fm)  # to test one of the unscaledre methods
     # must restore state of fm as it is cached in the global fittedmodels
     slp = MixedModels.dataset(:sleepstudy)