Repository: JuliaStats/MixedModels.jl (julia)

Issue:
`issingular` failed for GLMM with random-slopes
Hi, it seems that `issingular()` failed for GLMM with random-slopes, but not for GLMM with random-intercept only model (at least when `Bernoulli()` is used). 

For example:
```julia
using MixedModels # version v4.2.0

# intercept-only model
gm_intercept = fit(MixedModel,
	@formula(use ~ 1+urban+(1|dist)),
	MixedModels.dataset(:contra),
	Bernoulli();   
	contrasts=Dict(:dist => Grouping(), :urban => HelmertCoding()),
);

issingular(gm_intercept) # it works well

# add random-slope
gm_slope = fit(MixedModel,
	@formula(use ~ 1+urban+(1+urban|dist)),
	MixedModels.dataset(:contra),
	Bernoulli();   
	contrasts=Dict(:dist => Grouping(), :urban => HelmertCoding()),
);

issingular(gm_slope) # raises the error (see below)
```

`issingular(gm_intercept)` works well but `issingular(gm_slope)` raises the error:
`DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 5 and 3")`

Many thanks for your dedication on this great package!

Pull request description:
Closes #563

Patch:
diff --git a/Project.toml b/Project.toml
index 4fe4abb4b..f24ccbbb9 100644
--- a/Project.toml
+++ b/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.2.0"
+version = "4.2.1"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
diff --git a/src/mixedmodel.jl b/src/mixedmodel.jl
index 8a037df7b..a40d6c4b8 100644
--- a/src/mixedmodel.jl
+++ b/src/mixedmodel.jl
@@ -47,8 +47,13 @@ end
 Test whether the model `m` is singular if the parameter vector is `θ`.
 
 Equality comparisons are used b/c small non-negative θ values are replaced by 0 in `fit!`.
+
+!!! note
+    For `GeneralizedLinearMixedModel`, the entire parameter vector (including
+    β in the case `fast=false`) must be specified if the default is not used.
 """
 issingular(m::MixedModel, θ=m.θ) = any(lowerbd(m) .== θ)
+issingular(m::GeneralizedLinearMixedModel, θ=m.optsum.final) = any(lowerbd(m) .== θ)
 
 # FIXME: better to base this on m.optsum.returnvalue
 StatsBase.isfitted(m::MixedModel) = m.optsum.feval > 0


Test patch:
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
index 64aa067c6..95a326aa3 100644
--- a/test/bootstrap.jl
+++ b/test/bootstrap.jl
@@ -114,7 +114,7 @@ end
         contra = dataset(:contra)
         # need a model with fast=false to test that we only
         # copy the optimizer constraints for θ and not β
-        gm0 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(), fast=false, progress=false)
+        gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(), fast=false, progress=false)
         bs = parametricbootstrap(StableRNG(42), 100, gm0)
         # make sure we're not copying
         @test length(bs.lowerbd) == length(gm0.θ)
diff --git a/test/modelcache.jl b/test/modelcache.jl
index eaaeadee0..f482b4eb2 100644
--- a/test/modelcache.jl
+++ b/test/modelcache.jl
@@ -3,7 +3,8 @@ using MixedModels: dataset
 
 @isdefined(gfms) || const global gfms = Dict(
     :cbpp => [@formula((incid/hsz) ~ 1 + period + (1|herd))],
-    :contra => [@formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist))],
+    :contra => [@formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist)),
+                @formula(use ~ 1+urban+(1+urban|dist))], # see #563
     :grouseticks => [@formula(ticks ~ 1+year+ch+ (1|index) + (1|brood) + (1|location))],
     :verbagg => [@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item))],
 )
diff --git a/test/pirls.jl b/test/pirls.jl
index 0cd7c9fc9..47d4f41a5 100644
--- a/test/pirls.jl
+++ b/test/pirls.jl
@@ -14,7 +14,7 @@ include("modelcache.jl")
 @testset "contra" begin
     contra = dataset(:contra)
     thin = 5
-    gm0 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(); fast=true, progress=false, thin)
+    gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true, progress=false, thin)
     fitlog = gm0.optsum.fitlog
     @test length(fitlog) == (div(gm0.optsum.feval, thin) + 1) # for the initial value
     @test first(fitlog) == (gm0.optsum.initial, gm0.optsum.finitial)
@@ -41,7 +41,7 @@ include("modelcache.jl")
     @test Distribution(gm0) == Distribution(gm0.resp)
     @test Link(gm0) == Link(gm0.resp)
 
-    gm1 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(); progress=false);
+    gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); progress=false);
     @test isapprox(gm1.θ, [0.573054], atol=0.005)
     @test lowerbd(gm1) == vcat(fill(-Inf, 7), 0.)
     @test isapprox(deviance(gm1), 2361.54575, rtol=0.00001)
@@ -51,7 +51,7 @@ include("modelcache.jl")
     @test nobs(gm0) == 1934
     refit!(gm0; fast=true, nAGQ=7, progress=false)
     @test isapprox(deviance(gm0), 2360.9838, atol=0.001)
-    gm1 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(); nAGQ=7, progress=false)
+    gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); nAGQ=7, progress=false)
     @test isapprox(deviance(gm1), 2360.8760, atol=0.001)
     @test gm1.β == gm1.beta
     @test gm1.θ == gm1.theta
@@ -74,6 +74,11 @@ include("modelcache.jl")
     #@test isapprox(sum(gm1.resp.devresid), 2237.349, atol=0.1)
     show(IOBuffer(), gm1)
     show(IOBuffer(), BlockDescription(gm0))
+
+    gm_slope = fit(MixedModel, gfms[:contra][2], contra, Bernoulli(); progress=false);
+    @test !issingular(gm_slope)
+    @test issingular(gm_slope, zeros(5))
+
 end
 
 @testset "cbpp" begin
diff --git a/test/predict.jl b/test/predict.jl
index ba2ac478d..e4c0f1317 100644
--- a/test/predict.jl
+++ b/test/predict.jl
@@ -34,7 +34,7 @@ include("modelcache.jl")
 
     @testset "GLMM" begin
         contra = DataFrame(dataset(:contra))
-        m = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(); fast=true,
+        m = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true,
                 contrasts=Dict(:urban => EffectsCoding()), progress=false)
         mc = deepcopy(m)
         fit!(simulate!(StableRNG(42), mc); progress=false)
@@ -79,7 +79,7 @@ end
     @testset "GLMM" begin
         contra = dataset(:contra)
         for fast in [true, false]
-            gm0 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(); fast, progress=false)
+            gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast, progress=false)
 
             @test_throws ArgumentError predict(gm0, contra; type=:doh)