correlation of fixed effects (#274)

* correlation of fixed effects

* corr and vcov properties, add names to GLMM getproperties

Author: palday

src/MixedModels.jl

@@ -64,8 +64,8 @@ export @formula,
        coefnames,
        coeftable,
        cond,
-       describeblocks,
        condVar,
+       describeblocks,
        deviance,
        dispersion,
        dispersion_parameter,

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, :PCA, :rePCA, :optsum, :X, :reterms, :feterms, :formula, :σs, :σρs)
+    elseif s ∈ (:A, :L, :λ, :lowerbd, :corr, :vcov, :PCA, :rePCA, :optsum, :X, :reterms, :feterms, :formula, :σs, :σρs)
         getproperty(m.LMM, s)
     elseif s == :y
         m.resp.y
@@ -416,6 +416,10 @@ Base.propertynames(m::GeneralizedLinearMixedModel, private = false) = (
     :lowerbd,
     :σρs,
     :σs,
+    :corr,
+    :vcov,
+    :PCA,
+    :rePCA,
     fieldnames(typeof(m))...,
 )
 

src/linearmixedmodel.jl

@@ -450,6 +450,10 @@ function Base.getproperty(m::LinearMixedModel, s::Symbol)
         filter(Base.Fix2(isa, FeMat), getfield(m, :allterms))
     elseif s == :objective
         objective(m)
+    elseif s == :corr
+        vcov(m, corr=true)
+    elseif s == :vcov
+        vcov(m, corr=false)
     elseif s == :PCA
         NamedTuple{fnames(m)}(PCA.(m.reterms))
     elseif s == :pvalues
@@ -537,6 +541,8 @@ Base.propertynames(m::LinearMixedModel, private = false) = (
     :lowerbd,
     :X,
     :y,
+    :corr,
+    :vcov,
     :PCA,
     :rePCA,
     :reterms,
@@ -851,23 +857,30 @@ end
 Returns the estimate of σ², the variance of the conditional distribution of Y given B.
 """
 varest(m::LinearMixedModel) = pwrss(m) / dof_residual(m)
+"""
+    vcov(m::LinearMixedModel)
 
-function StatsBase.vcov(m::LinearMixedModel{T}) where {T}
+Returns the variance-covariance matrix of the fixed effects.
+If `corr=true`, then correlation of fixed effects is returned instead.
+"""
+function StatsBase.vcov(m::LinearMixedModel{T}; corr=false) 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]
+        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
-        covmat[iperm, iperm]
+        vv = covmat[iperm, iperm]
     end
+
+    corr ?  StatsBase.cov2cor!(vv,stderror(m)) : vv
 end
 
 """

test/pls.jl

@@ -207,6 +207,9 @@ end
     @test last(values(first(σs))) ≈ 5.716827808126002 atol=0.0001
     show(IOBuffer(), fm)
 
+    @test fm.vcov ≈ [43.986843599069985 -1.370391879315539; -1.370391879315539 2.2567113314443645] atol=0.0001
+    @test fm.corr ≈ [1.0 -0.1375451787621904; -0.1375451787621904 1.0] atol=0.0001
+
     u3 = ranef(fm, uscale=true)
     @test length(u3) == 1
     @test size(first(u3)) == (2, 18)