new rePCA: now on correlations and with optional extra info (#264)

* Added PCA property with correlation or covariance components and loadings.

* fix rowsnormalization for rows with all zeros

* PCA methods for ReMat, show method for PCA

* properties for PCA structure

* new rePCA

* remove redundant, misplaced definitions

* PCA methods for model and set appropriate signatures

* GLMM rePCA

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

Author: palday

src/generalizedlinearmixedmodel.jl

@@ -373,7 +373,7 @@ 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 ∈ (:A, :L, :λ, :lowerbd, :PCA, :rePCA, :optsum, :X, :reterms, :feterms, :formula, :σs, :σρs)
         getproperty(m.LMM, s)
     elseif s == :y
         m.resp.y
@@ -566,6 +566,8 @@ for f in (
     :feL,
     :(LinearAlgebra.logdet),
     :lowerbd,
+    :PCA,
+    :rePCA,
     :(StatsModels.modelmatrix),
     :(StatsBase.vcov),
     :σs,

src/linearmixedmodel.jl

@@ -450,6 +450,8 @@ function Base.getproperty(m::LinearMixedModel, s::Symbol)
         filter(Base.Fix2(isa, FeMat), getfield(m, :allterms))
     elseif s == :objective
         objective(m)
+    elseif s == :PCA
+        NamedTuple{fnames(m)}(PCA.(m.reterms))
     elseif s == :pvalues
         ccdf.(Chisq(1), abs2.(coef(m) ./ stderror(m)))
     elseif s == :reterms
@@ -535,6 +537,7 @@ Base.propertynames(m::LinearMixedModel, private = false) = (
     :lowerbd,
     :X,
     :y,
+    :PCA,
     :rePCA,
     :reterms,
     :feterms,
@@ -616,11 +619,31 @@ end
 
 LinearAlgebra.rank(m::LinearMixedModel) = first(m.feterms).rank
 
-function rePCA(m::LinearMixedModel{T}) where {T}
-    re = m.reterms
-    nt = length(re)
-    tup = ntuple(i -> normalized_variance_cumsum(re[i].λ), nt)
-    NamedTuple{ntuple(i -> re[i].trm.sym, nt),typeof(tup)}(tup)
+"""
+    rePCA(m::LinearMixedModel; corr::Bool=true)
+
+Return a named tuple of the normalized cumulative variance of a principal components
+analysis of the random effects covariance matrices or correlation
+matrices when `corr` is `true`.
+
+The normalized cumulative variance is the proportion of the variance for the first
+principal component, the first two principal components, etc.  The last element is
+always 1.0 representing the complete proportion of the variance.
+"""
+function rePCA(m::LinearMixedModel; corr::Bool=true)
+    pca = PCA.(m.reterms, corr=corr)
+    NamedTuple{fnames(m)}(getproperty.(pca,:cumvar))
+end
+
+"""
+    PCA(m::LinearMixedModel; corr::Bool=true)
+
+Return a named tuple of the principal components analysis of the random effects
+covariance matrices or correlation matrices when `corr` is `true`.
+"""
+
+function PCA(m::LinearMixedModel; corr::Bool=true)
+    NamedTuple{fnames(m)}(PCA.(m.reterms, corr=corr))
 end
 
 """

src/remat.jl

@@ -374,6 +374,11 @@ function *(adjA::Adjoint{T,<:ReMat{T,S}}, B::ReMat{T,P}) where {T,S,P}
     BlockedSparse{T,S,P}(cscmat, reshape(cscmat.nzval, S, :), cscmat.colptr[1:P:(cscmat.n + 1)])
 end
 
+PCA(A::ReMat{T,1}; corr::Bool=true) where {T} =
+    PCA(corr ? reshape(abs.(vec(A.λ)),1,1) : reshape(ones(T,1),1,1), corr=corr)
+
+PCA(A::ReMat{T,S}; corr::Bool=true) where {T,S} = PCA(A.λ, corr=corr)
+
 function reweight!(A::ReMat, sqrtwts::Vector)
     if length(sqrtwts) > 0
         if A.z === A.wtz

src/utilities.jl

@@ -59,14 +59,14 @@ function Base.sum!(s::AbstractVector{T}, a::RaggedArray{T}) where {T}
     s
 end
 
-"""
-    normalized_variance_cumsum(A::AbstractMatrix)
-
-Return the cumulative sum of the squared singular values of `A` normalized to sum to 1
-"""
-function normalized_variance_cumsum(A::AbstractMatrix)
-    vars = cumsum(abs2.(svdvals(A)))
-    vars ./ last(vars)
+function rownormalize!(A::AbstractMatrix)
+    for r in eachrow(A)
+        # all zeros arise in zerocorr situations
+        if !iszero(r)
+            normalize!(r)
+        end
+    end
+    A
 end
 
 """
@@ -137,3 +137,106 @@ dataset(nm::Symbol) = dataset(string(nm))
 Return a vector of names of the available test data sets
 """
 datasets() = first.(Base.Filesystem.splitext.(filter(Base.Fix2(endswith, ".feather"), readdir(TestData))))
+
+
+"""
+    PCA{T<:AbstractFloat}
+
+Principal Components Analysis
+
+## Fields
+
+* `covcorr` covariance or correlation matrix
+* `sv` singular value decomposition
+* `corr` is this a correlation matrix
+"""
+struct PCA{T<:AbstractFloat}
+    covcor::Symmetric{T,Matrix{T}}
+    sv::SVD{T,T,Matrix{T}}
+    corr::Bool
+end
+
+"""
+    PCA(::AbstractMatrix; corr::Bool=true)
+    PCA(::ReMat; corr::Bool=true)
+    PCA(::LinearMixedModel; corr::Bool=true)
+
+Constructs a [`MixedModels.PCA`](@ref]) object from a covariance matrix.
+
+For `LinearMixedModel`, a named tuple of PCA on each of the random-effects terms is returned.
+
+If `corr=true`, then the covariance is first standardized to the correlation scale.
+"""
+
+function PCA(covfac::AbstractMatrix; corr::Bool=true)
+    covf = corr ? rownormalize!(copy(covfac)) : copy(covfac)
+    PCA(Symmetric(covf*covf', :L), svd(covf), corr)
+end
+
+function Base.getproperty(pca::PCA, s::Symbol)
+    if s == :cumvar
+        cumvv = cumsum(abs2.(pca.sv.S))
+        cumvv ./ last(cumvv)
+    elseif s == :loadings
+        pca.sv.U
+    else
+        getfield(pca, s)
+    end
+end
+
+Base.propertynames(pca::PCA, private = false) = (
+    :covcor,
+    :sv,
+    :corr,
+    :cumvar,
+    :loadings,
+#    :rotation,
+)
+
+Base.show(pca::PCA;
+        ndigitsmat=2, ndigitsvec=2, ndigitscum=4,
+        covcor=true, loadings=true, variances=false, stddevs=false) =
+        Base.show(Base.stdout, pca,
+                    ndigitsmat=ndigitsmat,
+                    ndigitsvec=ndigitsvec,
+                    ndigitscum=ndigitscum,
+                    covcor=covcor,
+                    loadings=loadings,
+                    variances=variances,
+                    stddevs=stddevs)
+
+function Base.show(io::IO, pca::PCA;
+        ndigitsmat=2, ndigitsvec=2, ndigitscum=4,
+        covcor=true, loadings=true, variances=false, stddevs=false)
+    println(io)
+    if covcor
+        println(io,
+                "Principal components based on ",
+                pca.corr ? "correlation" : "(relative) covariance",
+                " matrix")
+        # only display the lower triangle of symmetric matrix
+        Base.print_matrix(io, round.(LowerTriangular(pca.covcor), digits=ndigitsmat))
+        println(io)
+    end
+    if stddevs
+        println(io, "Standard deviations:")
+        sv = pca.sv
+        show(io, round.(sv.S, digits=ndigitsvec))
+        println(io)
+    end
+    if variances
+        println(io, "Variances:")
+        vv = abs2.(sv.S)
+        show(io, round.(vv, digits=ndigitsvec))
+        println(io)
+    end
+    println(io, "Normalized cumulative variances:")
+    show(io, round.(pca.cumvar, digits=ndigitscum))
+    println(io)
+    if loadings
+        println(io, "Component loadings")
+        Base.print_matrix(io, round.(pca.loadings, digits=ndigitsmat))
+    end
+
+    nothing
+end

test/pirls.jl

@@ -23,6 +23,11 @@ using DataFrames, Feather, LinearAlgebra, MixedModels, Test
     @test isnan(gm1.σ)
     @test length(gm1.y) == size(gm1.X, 1)
     @test :θ in propertynames(gm0)
+    @testset "GLMM rePCA" begin
+        @test length(MixedModels.PCA(gm0)) == 1
+        @test length(MixedModels.rePCA(gm0)) == 1
+        @test length(gm0.rePCA) == 1
+    end
     # gm0.βθ = vcat(gm0.β, gm0.theta)
     # the next three values are not well defined in the optimization
     #@test isapprox(logdet(gm1), 75.7217, atol=0.1)

test/pls.jl

@@ -159,7 +159,11 @@ end
     resid1 = residuals(fm1);
     @test size(resid1) == (73421, )
     @test first(resid1) ≈ 1.82124 atol=0.00001
-    @test length(fm1.rePCA) == 3
+
+    @testset "PCA" begin
+        @test length(fm1.rePCA) == 3
+        @test length(MixedModels.PCA(fm1)) == 3
+    end
 
     fm2 = fit!(LinearMixedModel(@formula(y ~ 1 + service*dept + (1|s) + (1|d)), insteval))
     @test objective(fm2) ≈ 237585.5534151694 atol=0.001
@@ -228,6 +232,13 @@ end
     @test fmnc.θ == ones(2)
     @test lowerbd(fmnc) == zeros(2)
 
+    @testset "zerocorr PCA" begin
+        @test length(fmnc.rePCA) == 1
+        @test fmnc.rePCA.subj == [0.5, 1.0]
+        @test MixedModels.PCA(fmnc).subj.loadings == I(2)
+        @test show(IOBuffer(), MixedModels.PCA(fmnc)) === nothing
+    end
+
     fit!(fmnc)
 
     @test deviance(fmnc) ≈ 1752.0032551398835 atol=0.001