correct coeftable for GLMMs, add tests (#308)

* correct coeftable for GLMMs, add tests

* Move methods for abstract MixedModel struct to separate file.

(cherry picked from commit 5fb6f65)

Author: dmbates

src/MixedModels.jl

@@ -117,6 +117,7 @@ include("linearmixedmodel.jl")
 include("gausshermite.jl")
 include("generalizedlinearmixedmodel.jl")
 include("mixed.jl")
+include("mixedmodel.jl")
 include("linalg/statschol.jl")
 include("linalg/cholUnblocked.jl")
 include("linalg/rankUpdate.jl")

src/generalizedlinearmixedmodel.jl

@@ -53,6 +53,20 @@ struct GeneralizedLinearMixedModel{T <: AbstractFloat} <: MixedModel{T}
     mult::Vector{T}
 end
 
+function StatsBase.coeftable(m::GeneralizedLinearMixedModel)
+    co = fixef(m)
+    se = stderror(m)
+    z = co ./ se
+    pvalue = ccdf.(Chisq(1), abs2.(z))
+    CoefTable(
+        hcat(co, se, z, pvalue),
+        ["Estimate", "Std.Error", "z value", "P(>|z|)"],
+        coefnames(m),
+        4, # pvalcol
+        3, # teststatcol
+    )
+end
+
 """
     deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1)::T where {T}
 
@@ -68,7 +82,7 @@ function StatsBase.deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1) where {T}
     u = vec(first(m.u))
     u₀ = vec(first(m.u₀))
     copyto!(u₀, u)
-    ra = RaggedArray(m.resp.devresid, first(m.LMM.reterms).refs)
+    ra = RaggedArray(m.resp.devresid, first(m.LMM.allterms).refs)
     devc0 = sum!(map!(abs2, m.devc0, u), ra)  # the deviance components at z = 0
     sd = map!(inv, m.sd, m.LMM.L[Block(1,1)].diag)
     mult = fill!(m.mult, 0)
@@ -105,8 +119,17 @@ function deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)
     deviance(m, nAGQ)
 end
 
-GLM.dispersion(m::GeneralizedLinearMixedModel, sqr::Bool=false) =
-    dispersion(m.resp, dof_residual(m), sqr)
+function GLM.dispersion(m::GeneralizedLinearMixedModel{T}, sqr::Bool = false) where {T}
+# adapted from GLM.dispersion(::AbstractGLM, ::Bool)
+# TODO: PR for a GLM.dispersion(resp::GLM.GlmResp, dof_residual::Int, sqr::Bool)
+    r = m.resp
+    if dispersion_parameter(r.d)
+        s = sum(wt * abs2(re) for (wt, re) in zip(r.wrkwt, r.wrkresid)) / dof_residual(m)
+        sqr ? s : sqrt(s)
+    else
+        one(T)
+    end
+end
 
 GLM.dispersion_parameter(m::GeneralizedLinearMixedModel) = dispersion_parameter(m.resp.d)
 
@@ -295,7 +318,11 @@ function Base.getproperty(m::GeneralizedLinearMixedModel, s::Symbol)
         m.β
     elseif s ∈ (:σ, :sigma)
         sdest(m)
-    elseif s ∈ (:A, :L, :λ, :lowerbd, :optsum, :X, :reterms, :feterms, :formula, :σs, :σρs)
+    elseif s == :σs
+        σs(m)
+    elseif s == :σρs
+        σρs(m)
+    elseif s ∈ (:A, :L, :λ, :lowerbd, :corr, :PCA, :rePCA, :optsum, :X, :reterms, :feterms, :formula)
         getproperty(m.LMM, s)
     elseif s == :y
         m.resp.y
@@ -305,18 +332,17 @@ function Base.getproperty(m::GeneralizedLinearMixedModel, s::Symbol)
 end
 
 function StatsBase.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T}
-    accum = zero(T)
+    r = m.resp
     D = Distribution(m.resp)
-    if D <: Binomial
-        for (μ, y, n) in zip(m.resp.mu, m.resp.y, m.wt)
-            accum += logpdf(D(round(Int, n), μ), round(Int, y * n))
+    accum = (
+        if D <: Binomial
+            sum(logpdf(D(round(Int, n), μ), round(Int, y * n)) 
+                for (μ, y, n) in zip(r.mu, r.y, m.wt))
+        else
+            sum(logpdf(D(μ), y) for (μ, y) in zip(r.mu, r.y))
         end
-    else
-        for (μ, y) in zip(m.resp.mu, m.resp.y)
-            accum += logpdf(D(μ), y)
-        end
-    end
-    accum - (mapreduce(u -> sum(abs2, u), + , m.u) + logdet(m)) / 2
+    )
+    accum - (sum(sum(abs2, u) for u in m.u) + logdet(m)) / 2
 end
 
 StatsBase.nobs(m::GeneralizedLinearMixedModel) = length(m.η)
@@ -467,13 +493,14 @@ varest(m::GeneralizedLinearMixedModel{T}) where {T} = one(T)
 for f in (
     :describeblocks,
     :feL,
+    :fetrm,
     :(LinearAlgebra.logdet),
     :lowerbd,
+    :PCA,
+    :rePCA,
+    :(StatsBase.coefnames),
     :(StatsModels.modelmatrix),
-    :(StatsBase.vcov),
-    :σs,
-    :σρs,
-    )
+)
     @eval begin
         $f(m::GeneralizedLinearMixedModel) = $f(m.LMM)
     end

src/linearmixedmodel.jl

@@ -126,13 +126,6 @@ function StatsBase.coeftable(m::MixedModel)
         first(m.feterms).cnames, 4)
 end
 
-"""
-    cond(m::MixedModel)
-
-Return a vector of condition numbers of the λ matrices for the random-effects terms
-"""
-LinearAlgebra.cond(m::MixedModel) = cond.(m.λ)
-
 """
     condVar(m::LinearMixedModel)
 
@@ -231,12 +224,29 @@ function StatsBase.dof_residual(m::LinearMixedModel)::Int
 end
 
 """
-    feL(m::MixedModel)
+    feind(m::LinearMixedModel)
+
+An internal utility to return the index in `m.allterms` of the fixed-effects term.
+"""
+feind(m::LinearMixedModel) = findfirst(Base.Fix2(isa, FeMat), m.allterms)
+
+"""
+    feL(m::LinearMixedModel)
 
 Return the lower Cholesky factor for the fixed-effects parameters, as an `LowerTriangular`
 `p × p` matrix.
 """
-feL(m::LinearMixedModel) = LowerTriangular(m.L.blocks[end - 1, end - 1])
+function feL(m::LinearMixedModel)
+    k = feind(m)
+    LowerTriangular(m.L.blocks[k, k])
+end
+
+"""
+    fetrm(m::LinearMixedModel)
+
+Return the fixed-effects term from `m.allterms`
+"""
+fetrm(m::LinearMixedModel) = m.allterms[feind(m)]
 
 """
     fit!(m::LinearMixedModel[; verbose::Bool=false, REML::Bool=false])
@@ -636,16 +646,13 @@ function Base.show(io::IO, m::LinearMixedModel)
     show(io,coeftable(m))
 end
 
-function σs(m::LinearMixedModel)
-    σ = sdest(m)
-    NamedTuple{fnames(m)}(((σs(t, σ) for t in m.reterms)...,))
-end
-
-function σρs(m::LinearMixedModel)
-    σ = sdest(m)
-    NamedTuple{fnames(m)}(((σρs(t, σ) for t in m.reterms)...,))
-end
+"""
+    size(m::LinearMixedModel)
 
+Returns the size of a mixed model as a tuple of length four:
+the number of observations, the number of (non-singular) fixed-effects parameters,
+the number of conditional modes (random effects), the number of grouping variables
+"""
 function Base.size(m::LinearMixedModel)
     n, p = size(first(m.feterms))
     n, p, sum(size.(m.reterms, 2)), length(m.reterms)
@@ -660,10 +667,26 @@ This value is the contents of the `1 × 1` bottom right block of `m.L`
 """
 sqrtpwrss(m::LinearMixedModel) = first(m.L.blocks[end, end])
 
+
+"""
+    ssqdenom(m::LinearMixedModel)
+Return the denominator for penalized sums-of-squares.
+For MLE, this value is the number of observations. For REML, this
+value is the number of observations minus the rank of the fixed-effects matrix.
+The difference is analagous to the use of n or n-1 in the denominator when
+calculating the variance.
+"""
+function ssqdenom(m::LinearMixedModel)::Int
+    (n, p, q, k)  = size(m)
+    m.optsum.REML ? n - p : n
+end
+
 """
     std(m::MixedModel)
 
 Return the estimated standard deviations of the random effects as a `Vector{Vector{T}}`.
+
+FIXME: This uses an old convention of isfinite(sdest(m)).  Probably drop in favor of m.σs
 """
 function Statistics.std(m::LinearMixedModel)
     rl = rowlengths.(m.reterms)
@@ -739,25 +762,7 @@ end
 
 Returns the estimate of σ², the variance of the conditional distribution of Y given B.
 """
-varest(m::LinearMixedModel) = pwrss(m) / dof_residual(m)
-
-function StatsBase.vcov(m::LinearMixedModel{T}) where {T}
-    Xtrm = first(m.feterms)
-    iperm = invperm(Xtrm.piv)
-    p = length(iperm)
-    r = Xtrm.rank
-    Linv = inv(feL(m))
-    permvcov = varest(m) * (Linv'Linv)
-    if p == Xtrm.rank
-        permvcov[iperm, iperm]
-    else
-        covmat = fill(zero(T)/zero(T), (p, p))
-        for j in 1:r, i in 1:r
-            covmat[i,j] = permvcov[i, j]
-        end
-        covmat[iperm, iperm]
-    end
-end
+varest(m::LinearMixedModel) = pwrss(m) / ssqdenom(m)
 
 """
     zerocorr!(m::LinearMixedModel[, trmnms::Vector{Symbol}])

src/mixedmodel.jl

@@ -0,0 +1,44 @@
+
+"""
+cond(m::MixedModel)
+
+Return a vector of condition numbers of the λ matrices for the random-effects terms
+"""
+LinearAlgebra.cond(m::MixedModel) = cond.(m.λ)
+
+function σs(m::MixedModel)
+    σ = dispersion(m)
+    NamedTuple{fnames(m)}(((σs(t, σ) for t in m.reterms)...,))
+end
+
+function σρs(m::MixedModel)
+    σ = dispersion(m)
+    NamedTuple{fnames(m)}(((σρs(t, σ) for t in m.reterms)...,))
+end
+
+"""
+    vcov(m::LinearMixedModel)
+
+Returns the variance-covariance matrix of the fixed effects.
+If `corr=true`, then correlation of fixed effects is returned instead.
+"""
+function StatsBase.vcov(m::MixedModel; corr=false)
+    Xtrm = fetrm(m)
+    iperm = invperm(Xtrm.piv)
+    p = length(iperm)
+    r = Xtrm.rank
+    Linv = inv(feL(m))
+    T = eltype(Linv)
+    permvcov = dispersion(m, true) * (Linv'Linv)
+    if p == Xtrm.rank
+        vv = permvcov[iperm, iperm]
+    else
+        covmat = fill(zero(T) / zero(T), (p, p))
+        for j = 1:r, i = 1:r
+            covmat[i, j] = permvcov[i, j]
+        end
+        vv = covmat[iperm, iperm]
+    end
+
+    corr ?  StatsBase.cov2cor!(vv, stderror(m)) : vv
+end

test/pirls.jl

@@ -51,7 +51,7 @@ end
     gm2 = fit(MixedModel, @formula(prop ~ 1 + p + (1|h)), cbpp, Binomial(), wts = cbpp[!,:s])
     @test isapprox(deviance(gm2,true), 100.09585619324639, atol=0.0001)
     @test isapprox(sum(abs2, gm2.u[1]), 9.723175126731014, atol=0.0001)
-    @test isapprox(logdet(gm2), 16.90099, atol=0.0001)
+    @test isapprox(logdet(gm2), 16.90099, atol=0.001)
     @test isapprox(sum(gm2.resp.devresid), 73.47179193718736, atol=0.001)
     @test isapprox(loglikelihood(gm2), -92.02628186555876, atol=0.001)
     @test isnan(sdest(gm2))