Bugfixes and last set of backports for 2.x (#399)

* warning for GLMM with dispersion parameter (#373)* warning for GLMM with dispersion parameter(cherry picked from commit 66b8863)

* correct coeftable for GLMMs, add tests (#308)

* Move methods for abstract MixedModel struct to separate file.

(cherry picked from commit 5fb6f65)

* chang CoefTable column names to match GLM.jl

* ranef method for GLMM (#336)

(cherry picked from commit 753aa6f)

* make (1|x) + (y|x) equivalent to (1|x) + (0+y|x) (#338)

* fix error with (1|x) + (y|x)

* fix typo

(cherry picked from commit 20c83ef)

* update ranef docstring

Co-authored-by: Douglas Bates <[email protected]>

Author: palday

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

@@ -68,7 +68,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 +105,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)
 
@@ -233,9 +242,9 @@ GeneralizedLinearMixedModel(f::FormulaTerm, tbl,
 GeneralizedLinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable,
                             d::Normal,
                             l::IdentityLink;
-                            wts = [],
-                            offset = [],
-                            contrasts = Dict{Symbol,Any}()) =
+                            wts = wts,
+                            offset = offset,
+                            contrasts = contrasts) =
     throw(ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink"))
 function GeneralizedLinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable,
                                      d::Distribution,
@@ -247,7 +256,14 @@ function GeneralizedLinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable,
         d = Bernoulli()
     end
     (isa(d, Normal) && isa(l, IdentityLink)) &&
-        throw(ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink"))
+    throw(ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink"))
+
+    if !any(isa(d, dist) for dist in (Bernoulli, Binomial, Poisson))
+        @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."""
+    end
+
     LMM = LinearMixedModel(f, tbl, contrasts = contrasts; wts = wts)
     y = copy(LMM.y)
         # the sqrtwts field must be the correct length and type but we don't know those
@@ -288,7 +304,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
@@ -298,18 +318,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))
-        end
-    else
-        for (μ, y) in zip(m.resp.mu, m.resp.y)
-            accum += logpdf(D(μ), y)
+    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
-    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.η)
@@ -460,13 +479,14 @@ varest(m::GeneralizedLinearMixedModel{T}) where {T} = one(T)
 for f in (
     :describeblocks,
     :feL,
+    :fetrm,
     :(LinearAlgebra.logdet),
     :lowerbd,
+    :ranef,
+    :rePCA,
+    :(StatsBase.coefnames),
     :(StatsModels.modelmatrix),
-    :(StatsBase.vcov),
-    :σs,
-    :σρs,
-    )
+)
     @eval begin
         $f(m::GeneralizedLinearMixedModel) = $f(m.LMM)
     end

src/linearmixedmodel.jl

@@ -122,17 +122,15 @@ function StatsBase.coeftable(m::MixedModel)
     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|)"],
-        first(m.feterms).cnames, 4)
+    CoefTable(
+        hcat(co, se, z, pvalue),
+        ["Coef.", "Std. Error", "z", "Pr(>|z|)"],
+        coefnames(m),
+        4, # pvalcol
+        3, # teststatcol
+    )
 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 +229,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])
@@ -506,7 +521,7 @@ end
 ranef!(v::Vector, m::LinearMixedModel, uscale::Bool) = ranef!(v, m, fixef(m), uscale)
 
 """
-    ranef(m::LinearMixedModel; uscale=false) #, named=true)
+    ranef(m::MixedModel; uscale=false, named=false)
 
 Return, as a `Vector{Vector{T}}` (`Vector{NamedVector{T}}` if `named=true`),
 the conditional modes of the random effects in model `m`.
@@ -636,16 +651,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 +672,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 +767,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/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NamedArrays = "86f7a689-2022-50b4-a561-43c23ac3c673"
 RData = "df47a6cb-8c03-5eed-afd8-b6050d6c41da"

test/pirls.jl

@@ -1,5 +1,7 @@
 using DataFrames, LinearAlgebra, MixedModels, RData, Test
 
+using GLM: SqrtLink
+
 if !@isdefined(dat) || !isa(dat, Dict{Symbol, DataFrame})
     const dat = Dict(Symbol(k) => v for (k, v) in
         load(joinpath(dirname(pathof(MixedModels)), "..", "test", "dat.rda")))
@@ -12,12 +14,19 @@ end
     gm0 = fit(MixedModel, contraform, contra, Bernoulli(), fast=true);
     @test gm0.lowerbd == zeros(1)
     @test isapprox(gm0.θ, [0.5720734451352923], atol=0.001)
-    @test isapprox(deviance(gm0,true), 2361.657188518064, atol=0.001)
+    @test isapprox(deviance(gm0), 2361.657188518064, atol=0.001)
+    # the first 9 BLUPs -- I don't think there's much point in testing all 102
+    blups = [-0.9546542382159384, -0.9243280740976597, -0.5994415212823175, 
+             -0.5853637765024964, -0.5582547143159903, -0.5230879394455542, 
+             -0.5180165911218682, -0.5175868781594128, -0.4955043890201569]
+    @test sort(vec(first(ranef(gm0))))[1:9] ≈ blups atol=1e-4
+
     gm1 = fit(MixedModel, contraform, contra, Bernoulli());
     @test isapprox(gm1.θ, [0.573054], atol=0.005)
     @test lowerbd(gm1) == vcat(fill(-Inf, 7), 0.)
     @test isapprox(deviance(gm1,true), 2361.54575, rtol=0.00001)
     @test isapprox(loglikelihood(gm1), -1180.77288, rtol=0.00001)
+
     @test dof(gm0) == length(gm0.β) + length(gm0.θ)
     @test nobs(gm0) == 1934
     fit!(gm0, fast=true, nAGQ=7)
@@ -29,13 +38,11 @@ end
     @test isnan(gm0.σ)
     @test length(gm0.y) == size(gm0.X, 1)
     @test :θ in propertynames(gm0)
-    gm0.β = gm0.beta
-    @test gm0.β == gm0.beta
-    gm0.θ = gm0.theta
-    @test gm0.θ == gm0.theta
-    gm0.βθ = vcat(gm0.β, gm0.theta)
-    @test gm0.β == gm0.beta
-    @test gm0.θ == gm0.theta
+
+    @testset "GLMM rePCA" begin
+        @test length(MixedModels.rePCA(gm0)) == 1
+        @test length(gm0.rePCA) == 1
+    end
     # the next three values are not well defined in the optimization
     #@test isapprox(logdet(gm1), 75.7217, atol=0.1)
     #@test isapprox(sum(abs2, gm1.u[1]), 48.4747, atol=0.1)
@@ -49,7 +56,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))
@@ -76,3 +83,13 @@ end
     #@test isapprox(sum(x -> sum(abs2, x), gm4.u), 196.8695297987013, atol=0.1)
     #@test isapprox(sum(gm4.resp.devresid), 220.92685781326136, atol=0.1)
 end
+
+@testset "dispersion" begin
+
+    form = @formula(Y ~ 1 + U + (1 + U | G))
+
+    @test_logs (:warn, r"dispersion parameter") GeneralizedLinearMixedModel(form, dat[:sleepstudy], Gamma())
+    @test_logs (:warn, r"dispersion parameter") GeneralizedLinearMixedModel(form, dat[:sleepstudy], InverseGaussian())
+    @test_logs (:warn, r"dispersion parameter") GeneralizedLinearMixedModel(form, dat[:sleepstudy], Normal(), SqrtLink())
+
+end