Improved handling of empty and rank-deficient fixed effects (#253)

* rankUpdate! for diagonal A

* add warnings for rank-deficient FE

* fix ordering of FE names in coeftable (undo pivot)

* full FeMat support for empty FE

* fix coefnames and add fixefnames for accessing cnames and respecting pivots

* coefnames and fixefnames for GLMM

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

Author: palday

src/MixedModels.jl

@@ -75,6 +75,7 @@ export @formula,
        fit!,
        fitted,
        fixef,
+       fixefnames,
        fulldummy,
        fnames,
        GHnorm,

src/femat.jl

@@ -20,10 +20,24 @@ end
 
 function FeMat(X::AbstractMatrix, cnms)
     T = eltype(X)
-    ch = statscholesky(Symmetric(X'X))
-    pivot = ch.piv
-    Xp = all(pivot .== 1:size(X, 2)) ? X : X[:, ch.piv]
-    FeMat{T,typeof(X)}(Xp, Xp, pivot, ch.rank, cnms[pivot])
+    if size(X,2) > 0
+        ch = statscholesky(Symmetric(X'X))
+        pivot = ch.piv
+        rank = ch.rank
+        Xp = all(pivot .== 1:size(X, 2)) ? X : X[:, ch.piv]
+        if rank < length(pivot)
+            @warn "fixed-effects matrix is rank deficient"
+        end
+    else
+        # although it doesn't take long for an empty matrix,
+        # we can still skip the Cholesky step, which gets the rank
+        # wrong anyway
+        pivot = Int[]
+        Xp = X
+        rank = 0
+    end
+    # ch.rank is wrong for empty FE
+    FeMat{T,typeof(X)}(Xp, Xp, pivot, rank, cnms[pivot])
 end
 
 function reweight!(A::FeMat{T}, sqrtwts::Vector{T}) where {T}

src/generalizedlinearmixedmodel.jl

@@ -53,6 +53,9 @@ struct GeneralizedLinearMixedModel{T<:AbstractFloat} <: MixedModel{T}
     mult::Vector{T}
 end
 
+StatsBase.coefnames(m::GeneralizedLinearMixedModel) = fixefnames(m, false)
+
+
 """
     deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1)::T where {T}
 
@@ -271,6 +274,10 @@ function fixef(m::GeneralizedLinearMixedModel{T}, permuted = true) where {T}
     invpermute!(v, piv)
 end
 
+function fixefnames(m::GeneralizedLinearMixedModel{T}, permuted = true) where {T}
+    fixefnames(m.LMM, permuted)
+end
+
 GeneralizedLinearMixedModel(
     f::FormulaTerm,
     tbl,

src/linalg/rankUpdate.jl

@@ -122,3 +122,20 @@ function rankUpdate!(
     C.diag .= β .* C.diag .+ α .* abs2.(A.diag)
     C
 end
+
+
+function rankUpdate!(
+    C::HermOrSym{T,Matrix{T}},
+    A::Diagonal{T},
+    α::Number = true,
+    β::Number = true,
+) where {T}
+    m, n = size(A)
+    m == size(C, 2) || throw(DimensionMismatch(""))
+    adiag = A.diag
+    cdata = C.data
+    for (i, a) in zip(diagind(cdata), adiag)
+        cdata[i] = β * cdata[i] + α * abs2(a)
+    end
+    C
+end

src/linearmixedmodel.jl

@@ -161,19 +161,21 @@ fit(
 
 StatsBase.coef(m::MixedModel) = fixef(m, false)
 
-βs(m::LinearMixedModel) = NamedTuple{(Symbol.(first(m.feterms).cnames)...,)}((fixef(m)...,))
+βs(m::LinearMixedModel) = NamedTuple{(Symbol.(fixefnames(m))...,)}((fixef(m)...,))
 
-StatsBase.coefnames(m::LinearMixedModel) = first(m.feterms).cnames
+StatsBase.coefnames(m::LinearMixedModel) = fixefnames(m, false)
 
 function StatsBase.coeftable(m::MixedModel)
     co = coef(m)
     se = stderror(m)
     z = co ./ se
     pvalue = ccdf.(Chisq(1), abs2.(z))
+    names = coefnames(m)
+
     CoefTable(
         hcat(co, se, z, pvalue),
         ["Estimate", "Std.Error", "z value", "P(>|z|)"],
-        first(m.feterms).cnames,
+        names,
         4,
     )
 end
@@ -380,6 +382,28 @@ function fixef(m::LinearMixedModel{T}, permuted = true) where {T}
     val
 end
 
+"""
+    fixefnames(m::MixedModel, permuted=true)
+
+Return the associated term names for fixed-effects parameter vector estimate of `m`.
+
+If `permuted` is `true` the vector elements are permuted according to
+`first(m.feterms).piv` and truncated to the rank of that term.
+"""
+function fixefnames(m::LinearMixedModel{T}, permuted = true) where {T}
+    Xtrm = first(m.feterms)
+    val = copy(Xtrm.cnames)
+    if !permuted
+        piv = Xtrm.piv
+        p = length(piv)
+        invpermute!(val, piv)
+    else
+        val = val[1:Xtrm.rank]
+    end
+
+    val
+end
+
 """
     fnames(m::MixedModel)
 

test/linalg.jl

@@ -46,6 +46,9 @@ end
     d2 = Diagonal(fill(2., 2))
     @test Diagonal(fill(5.,2)) == MixedModels.rankUpdate!(Diagonal(ones(2)), d2, 1.)
     @test Diagonal(fill(-3.,2)) == MixedModels.rankUpdate!(Diagonal(ones(2)), d2, -1.)
+
+    # when converting straight from diagonal to symmetric, the type is different
+    @test Diagonal(fill(5.,2)) == MixedModels.rankUpdate!(Symmetric(Matrix(1. * I(2)), :L), d2)
 end
 
 @testset "lmulλ!" begin

test/matrixterm.jl

@@ -15,4 +15,12 @@ using LinearAlgebra, MixedModels, Random, Test
     wts = rand(MersenneTwister(123454321), 30)
     MixedModels.reweight!(trm, wts)
     @test mul!(prd, trm', trm)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)
+
+    # empty fixed effects
+    trm = MixedModels.FeMat(ones(10,0), String[])
+    @test size(trm) == (10, 0)
+    @test length(trm) == 0
+    @test size(trm') == (0, 10)
+    @test eltype(trm) == Float64
+    @test trm.rank == 0
 end

test/pls.jl

@@ -340,11 +340,27 @@ end
 end
 
 @testset "Rank deficient" begin
-    Random.seed!(0)
-    x = rand(100)
-    data = columntable((x = x, x2 = 1.5 .* x, y = rand(100), z = repeat('A':'T', 5)))
+    rng = MersenneTwister(0)
+    x = rand(rng, 100)
+    data = columntable((x = x, x2 = 1.5 .* x, y = rand(rng, 100), z = repeat('A':'T', 5)))
     model = fit!(LinearMixedModel(@formula(y ~ x + x2 + (1|z)), data))
     @test length(fixef(model)) == 2
     @test rank(model) == 2
     @test length(coef(model)) == 3
+    ct = coeftable(model)
+    @test ct.rownms ==  ["(Intercept)", "x", "x2"]
+    @test fixefnames(model, true) == ["(Intercept)", "x"]
+    @test fixefnames(model, false) == ["(Intercept)", "x", "x2"]
+    @test last(coef(model)) == -0.0
+
+    # check preserving of name ordering in coeftable and placement of
+    # pivoted-out element
+    fill!(data.x, 0)
+    model2 = fit!(LinearMixedModel(@formula(y ~ x + x2 + (1|z)), data))
+    ct = coeftable(model2)
+    @test ct.rownms ==  ["(Intercept)", "x", "x2"]
+    @test fixefnames(model2, true) == ["(Intercept)", "x2"]
+    @test fixefnames(model2, false) == ["(Intercept)", "x", "x2"]
+    @test coef(model2)[2] == -0.0
+    @test last(fixef(model)) ≈ (last(fixef(model2)) * 1.5)
 end