Repository: JuliaStats/MixedModels.jl (julia)

Issue:
Make Distribution/Family Part of the Type Parameter for GLMM
This will allow dispatching on particular families.

Pull request:
Closes #463

Patch:
diff --git a/NEWS.md b/NEWS.md
index 5aa63403d..075f9fb94 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -18,6 +18,9 @@ MixedModels v4.0.0 Release Notes
   or `Diagonal`. This representation is both more memory efficient and
   enables additional computational optimizations for particular covariance
   structures.[#489]
+* `GeneralizedLinearMixedModel` now includes the response distribution as one
+  of its type parameters. This will allow dispatching on the model family and may allow
+  additional specialization in the future.[#490]
 
 Run-time formula syntax
 -----------------------
@@ -198,3 +201,4 @@ Package dependencies
 [#484]: https://github.com/JuliaStats/MixedModels.jl/issues/484
 [#486]: https://github.com/JuliaStats/MixedModels.jl/issues/486
 [#489]: https://github.com/JuliaStats/MixedModels.jl/issues/489
+[#490]: https://github.com/JuliaStats/MixedModels.jl/issues/490
diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl
index 7cd82362a..ceae64122 100644
--- a/src/generalizedlinearmixedmodel.jl
+++ b/src/generalizedlinearmixedmodel.jl
@@ -36,7 +36,7 @@ In addition to the fieldnames, the following names are also accessible through t
 - `y`: response vector
 
 """
-struct GeneralizedLinearMixedModel{T<:AbstractFloat} <: MixedModel{T}
+struct GeneralizedLinearMixedModel{T<:AbstractFloat, D<:Distribution} <: MixedModel{T}
     LMM::LinearMixedModel{T}
     β::Vector{T}
     β₀::Vector{T}
@@ -132,7 +132,7 @@ end
 
 """
     deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)
-    
+
 Update `m.η`, `m.μ`, etc., install the working response and working weights in
 `m.LMM`, update `m.LMM.A` and `m.LMM.R`, then evaluate the [`deviance`](@ref).
 """
@@ -157,6 +157,8 @@ end
 
 GLM.dispersion_parameter(m::GeneralizedLinearMixedModel) = dispersion_parameter(m.resp.d)
 
+Distributions.Distribution(m::GeneralizedLinearMixedModel{T,D}) where {T,D} = D
+
 function StatsBase.dof(m::GeneralizedLinearMixedModel)::Int
     length(m.β) + length(m.θ) + GLM.dispersion_parameter(m.resp.d)
 end
@@ -391,7 +393,7 @@ function GeneralizedLinearMixedModel(
         # it is empty unless there is a single random-effects term
     vv = length(u) == 1 ? vec(first(u)) : similar(y, 0)
 
-    res = GeneralizedLinearMixedModel(
+    res = GeneralizedLinearMixedModel{T, typeof(d)}(
         LMM,
         β,
         copy(β),
@@ -443,6 +445,8 @@ end
 getθ(m::GeneralizedLinearMixedModel)  = copy(m.θ)
 getθ!(v::AbstractVector{T}, m::GeneralizedLinearMixedModel{T}) where {T} = copyto!(v, m.θ)
 
+GLM.Link(m::GeneralizedLinearMixedModel) = GLM.Link(m.resp)
+
 function StatsBase.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T}
     accum = zero(T)
     # adapted from GLM.jl
@@ -647,7 +651,7 @@ For Gaussian models, this parameter is often called σ.
 """
 sdest(m::GeneralizedLinearMixedModel{T}) where {T} =  dispersion_parameter(m) ? dispersion(m, false) : missing
 
-function Base.show(io::IO, ::MIME"text/plain", m::GeneralizedLinearMixedModel)
+function Base.show(io::IO, ::MIME"text/plain", m::GeneralizedLinearMixedModel{T,D}) where {T,D}
     if m.optsum.feval < 0
         @warn("Model has not been fit")
         return nothing
@@ -655,8 +659,8 @@ function Base.show(io::IO, ::MIME"text/plain", m::GeneralizedLinearMixedModel)
     nAGQ = m.LMM.optsum.nAGQ
     println(io, "Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = $nAGQ)")
     println(io, "  ", m.LMM.formula)
-    println(io, "  Distribution: ", Distribution(m.resp))
-    println(io, "  Link: ", GLM.Link(m.resp), "\n")
+    println(io, "  Distribution: ", D)
+    println(io, "  Link: ", Link(m), "\n")
     println(io)
     nums = Ryu.writefixed.([loglikelihood(m), deviance(m), aic(m), aicc(m), bic(m)], 4)
     fieldwd = max(maximum(textwidth.(nums)) + 1, 11)


Test patch:
diff --git a/src/likelihoodratiotest.jl b/src/likelihoodratiotest.jl
index 8a11bcc9d..507435009 100644
--- a/src/likelihoodratiotest.jl
+++ b/src/likelihoodratiotest.jl
@@ -168,12 +168,10 @@ end
 
 function _iscomparable(m::GeneralizedLinearMixedModel...)
         # TODO: test that all models are fit with same fast/nAGQ option?
-        glms = getproperty.(m,:resp);
-
-        allequal(Distribution.(glms)) ||
+        allequal(Distribution.(m)) ||
             throw(ArgumentError("Models must be fit to the same distribution"))
 
-        allequal(string.(Link.(glms))) ||
+        allequal(string.(Link.(m))) ||
             throw(ArgumentError("Models must have the same link function"))
 
         allequal(nobs.(m)) ||
diff --git a/test/Project.toml b/test/Project.toml
index e7b63dd80..281ed3891 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/test/likelihoodratiotest.jl b/test/likelihoodratiotest.jl
index 17c626970..80ec31096 100644
--- a/test/likelihoodratiotest.jl
+++ b/test/likelihoodratiotest.jl
@@ -10,10 +10,10 @@ include("modelcache.jl")
 @testset "isnested" begin
     slp = dataset(:sleepstudy)
 
-    # these tests don't actually depend on the models being fit, 
+    # these tests don't actually depend on the models being fit,
     # so we just construct them
 
-    # mismatched RE terms 
+    # mismatched RE terms
     m1 = LinearMixedModel(@formula(reaction ~ 1 + days + (1+days|subj)), slp)
     m2 = LinearMixedModel(@formula(reaction ~ 1 + days + (0+days|subj)), slp)
     @test !isnested(m1, m2)
@@ -29,7 +29,7 @@ include("modelcache.jl")
     m2 = LinearMixedModel(@formula(rt_trunc ~ 1 + (1|item)), kb07)
     @test !isnested(m1, m2)
 
-    # fixed-effects specification in REML and 
+    # fixed-effects specification in REML and
     # conversion of internal ArgumentError into @error for StatsModels.isnested
     kb07  = dataset(:kb07)
     m1 = fit(MixedModel, @formula(rt_trunc ~ 1 + prec + (1|subj)), kb07, REML=true)
@@ -42,8 +42,6 @@ end
 
 @testset "likelihoodratio test" begin
     slp = dataset(:sleepstudy);
-    
-
 
     fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp);
     fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp);
@@ -67,7 +65,7 @@ end
 
     # differing FE with REML
     fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp, REML=true);
- 
+
     @test_throws ArgumentError likelihoodratiotest(fm0,fm1)
 
     contra = MixedModels.dataset(:contra);
diff --git a/test/pirls.jl b/test/pirls.jl
index 9d979dc4e..809fd295e 100644
--- a/test/pirls.jl
+++ b/test/pirls.jl
@@ -1,10 +1,12 @@
 using DataFrames
+using Distributions
 using MixedModels
 using PooledArrays
 using StableRNGs
 using Tables
 using Test
 
+using GLM: Link
 using MixedModels: dataset
 
 include("modelcache.jl")
@@ -32,6 +34,8 @@ include("modelcache.jl")
     @test sdest(gm0) === missing
     @test varest(gm0) === missing
     @test gm0.σ === missing
+    @test Distribution(gm0) == Distribution(gm0.resp)
+    @test Link(gm0) == Link(gm0.resp)
 
     gm1 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli());
     @test isapprox(gm1.θ, [0.573054], atol=0.005)