another potential solution for non-id gauss

JuliaStats · Oct 14, 2020 · c1bd575 · c1bd575
1 parent 035b4a1
commit c1bd575
Show file tree

Hide file tree

Showing 2 changed files with 42 additions and 12 deletions.
diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl
@@ -83,8 +83,11 @@ 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}
-    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))
+    dispersion_parameter(m) && nAGQ == 1 &&
+        return T(2 * (-loglikelihood(m) + dof(m) - 1) +  logdet(m) + sum(u -> sum(abs2, u), m.u))
+
+    nAGQ == 1 && return T(sum(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)
@@ -529,7 +532,7 @@ function pirls!(
             end
             varyβ && map!(average, β, β, β₀)
             obj = deviance!(m)
-            verbose && println(lpad(nhalf, 8), ", ", obj)
+            verbose && println(lpad(nhalf, 8), ", ", obj, " ", β)
         end
         if isapprox(obj, obj₀; atol = 0.00001)
             break

diff --git a/test/pirls.jl b/test/pirls.jl
@@ -1,10 +1,13 @@
 using DataFrames
 using Distributions
 using MixedModels
+using PooledArrays
+using StableRNGs
 using Tables
 using Test
 
 using MixedModels: dataset
+using GLM: linkinv
 
 const gfms = Dict(
     :cbpp => [@formula((incid/hsz) ~ 1 + period + (1|herd))],
@@ -138,17 +141,28 @@ end
     @test only(m11.θ) ≈ 1.838245201739852 atol=1.e-5
 end
 
+@testset "dispersion parameter" begin
 
-# 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)
-
+    # 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)
+
+    ## simulate some data ##
+    rng = StableRNG(42)
+    ng = 25
+    ns = 500
+    # random effect
+    σ = 1
+    σre  = 0.5
+    u = repeat(randn(rng, ns) .* σre, ng)
+    id = PooledArray(string.(repeat(1:ns, ng)))
+    # fixed effect
+    x = rand(rng, ns*ng)
 
-@testset "dispersion parameter" begin
     @testset "Gaussian with non identity link" begin
         dyestuff = MixedModels.dataset(:dyestuff)
         gauss = fit(MixedModel, only(gfms[:dyestuff]), dyestuff, Normal(), LogLink(), fast=true)
@@ -160,5 +174,18 @@ lme4deviance(x) = sum(x.resp.devresid)
         @test lme4deviance(gauss) ≈ 115187.5
         @test only(gauss.θ) ≈ 0.0
         @test only(gauss.beta) ≈ 7.331387691046023
+        # this works, but starting from a clean fit doesn't
+        # I'm not sure why the optimizer gets lost
+        refit!(gauss, fast=false)
+
+        rng = StableRNG(42)
+        link = InverseLink()
+        offset = abs(minimum(u)) + 1
+        y =  map(d ->  rand(rng, d), Normal.(linkinv.(link, u .+ offset), σ))
+        dat = (u=u, id=id, x=x, y=y)
+        gaussim = GeneralizedLinearMixedModel(@formula(y ~ 1 + (1|id)), dat,
+                                              Normal(), InverseLink())
+        fit!(gaussim, fast=true)
+        @test only(gaussim.β) ≈ linkfun(link, mean(y))
     end
 end