use slow loglikelihood for Gaussian GLMM

Author: palday

src/MixedModels.jl

@@ -47,6 +47,7 @@ export @formula,
        IdentityLink,
        InverseGaussian,
        InverseLink,
+       InverseSquareLink,
        LinearMixedModel,
        LogitLink,
        LogLink,

src/generalizedlinearmixedmodel.jl

@@ -83,7 +83,8 @@ is the squared length of the conditional modes, `u`, plus the determinant
 of `Λ'Z'WZΛ + I`, plus the sum of the squared deviance residuals.
 """
 function StatsBase.deviance(m::GeneralizedLinearMixedModel{T}, nAGQ = 1) where {T}
-    nAGQ == 1 && return T(sum(m.resp.devresid) + logdet(m) + sum(u -> sum(abs2, u), m.u))
+    m.resp.d isa Normal && return -2 * loglikelihood(m)
+    nAGQ == 1 && return T(sum(x -> x^2, m.resp.devresid) + logdet(m) + sum(u -> sum(abs2, u), m.u))
     u = vec(first(m.u))
     u₀ = vec(first(m.u₀))
     copyto!(u₀, u)
@@ -341,10 +342,16 @@ function GeneralizedLinearMixedModel(
     (isa(d, Normal) && isa(l, IdentityLink)) &&
     throw(ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink"))
 
-    if !any(isa(d, dist) for dist in (Bernoulli, Binomial, Poisson))
+    if any(isa(d, dist) for dist in (InverseGaussian, Gamma))
         @warn """Results for families with a dispersion parameter are not reliable.
                  It is best to avoid trying to fit such models in MixedModels until
                  the authors get a better understanding of those cases."""
+    elseif isa(d, Normal)
+        @warn """Results for Normal distribution with non-identity link
+                 use a slow loglikelihood computation to approximate the deviance.
+                 As such fitting such models in MixedModels will be substantially
+                 slower than other families until the authors get a better
+                 understanding of those cases."""
     end
 
     LMM = LinearMixedModel(f, tbl, contrasts = contrasts; wts = wts)

test/pirls.jl

@@ -1,7 +1,6 @@
 using DataFrames
 using Distributions
 using MixedModels
-using Random
 using Tables
 using Test
 
@@ -139,25 +138,27 @@ end
     @test only(m11.θ) ≈ 1.838245201739852 atol=1.e-5
 end
 
+
+# NB: "deviance" is complicated in lme4
+# there are several "deviances" defined:
+# https://github.com/lme4/lme4/issues/375#issuecomment-214494445
+# this is the way deviance(glmm) is computed in lme4
+# also called devianceCondRel in
+# https://github.com/lme4/lme4/blob/master/misc/logLikGLMM/logLikGLMM.R
+lme4deviance(x) = sum(x.resp.devresid)
+
+
 @testset "dispersion parameter" begin
-    @testset "gaussian with non identity link" begin
+    @testset "Gaussian with non identity link" begin
         dyestuff = MixedModels.dataset(:dyestuff)
-
-        # gauss = fit(MixedModel, only(gfms[:dyestuff]), dyestuff, Normal(), LogLink(), fast=true)
-        # @test only(gauss.θ) ≈atol=0.001 # value from lme4gaus
-        # # this corresponds to the deviance(gauss)  from lme4
-        # # NB: "deviance" is complicated in lme4
-        # # this is the way deviance(glmm) is computed in lme4
-        # # but there are several "deviances" defined:
-        # # https://github.com/lme4/lme4/issues/375#issuecomment-214494445
-        # @test sum(gauss.resp.devresid) ≈ . atol=0.001 # value from lme4
-
-        # # bobyqa fails here
-        # neldermead = GeneralizedLinearMixedModel(only(gfms[:dyestuff]), dyestuff, Normal(), LogLink());
-        # neldermead.optsum.optimizer = :LN_NELDERMEAD;
-        # fit!(neldermead)
-        # @test only(neldermead.θ) ≈  atol=0.001 # value from lme4
-        # @test sum(neldermead.resp.devresid) ≈ . atol=0.001 # value from lme4
+        gauss = fit(MixedModel, only(gfms[:dyestuff]), dyestuff, Normal(), LogLink(), fast=true)
+
+        # reference values from lme4
+        # note that the LL is shifted by exactly (?!) one in Julia
+        # I don't know why
+        @test loglikelihood(gauss) ≈ -167.36494298739052 + 1
+        @test lme4deviance(gauss) ≈ 115187.5
+        @test only(gauss.θ) ≈ 0.0
+        @test only(gauss.beta) ≈ 7.331387691046023
     end
-
 end