In dmbates/RePsychLing: Data sets from Psychology and Linguistics experiments

library(lme4)
library(RePsychLing)
library(knitr)
opts_chunk$set(cache=FALSE,comment=NA)
options(show.signif.stars=FALSE,width=92,digits = 5)

In a linear mixed model (LMM) incorporating vector-valued random effects, say by-subject random effects for intercept and for slope, the variance component parameters determine a variance-covariance matrix for these random effects.

As described in @Bates:Maechler:Bolker:Walker:2015 the parameters used in fitting the model are the entries in the Cholesky factor, $\Lambda$, of the relative variance-covariance matrix of the unconditional distribution of the random effects.

Consider a maximal model fit to the KWDYZ data from @Kliegl:Wei:Dambacher:Yan:Zhou:2011, available in this package,

summary(m0 <- lmer(rt ~ 1+c1+c2+c3 + (1+c1+c2+c3|subj), KWDYZ, REML=FALSE))

The parameter vector, $\theta$, for this model

zapsmall(getME(m0,"theta"))

has a value that is effectively zero in the last position.

(The zapsmall function, used here and in what follows, sets the format for printing a numeric vector or matrix so that very small elements do not cause a shift to scientific notation. In scientific notation it is more difficult to see at a glance which numbers are large and which are small.)

These 10 parameter values are the values on and below the diagonal of a lower triangular Cholesky factor

zapsmall(chf <- getME(m0,"Tlist")[[1]])

The $\theta$ vector elements fill the lower triangular matrix in column major order. That is, the first four elements of the vector are the first column, the next three are the elements on and below the diagonal in the second column, and so on.

The relative covariance matrix for the random effects is $\Lambda\Lambda'$, which can be evaluated as

tcrossprod(chf)

(The expression crossprod(X) forms $X'X$ and tcrossprod(X) forms $XX'$.)

To reproduce the covariance matrix tcrossprod(chf) must be scaled by $s^2$

tcrossprod(getME(m0,"sigma")*chf)

(Compare the diagonal entries of this result with the variance components of the random effects listed in the model summary.)

The (unconditional) correlation matrix of the random effects is obtained by scaling each row of $\Lambda$ to have unit length, then applying tcrossprod.

(rowlengths <- sqrt(rowSums(chf*chf)))

tcrossprod(chf/rowlengths)

(Compare the off-diagonal elements of this matrix with the correlations in the model summary.)

Singularity in the Cholesky factor

Some of the material in this section gets a bit technical. Don't be too concerned if parts seem unintelligible.

Because the last column of chf is the zero vector, chf is rank-deficient. That is, although chf is a 4x4 matrix, the linear subspace formed by all possible linear combinations of the columns is 3-dimensional. The random-effects vectors that can be generated from this fitted model must lie in this 3-dimensional subspace. Thus there will be no variability in one direction of the space of random effects.

The singular value decomposition of chf

(svd0 <- svd(chf,nv=0))

returns the singular values, d, and the matrix of left singular vectors, u. In the language of principal components analysis (PCA), the columns of u are the component loadings. These columns have unit length and are mutually orthogonal.

zapsmall(crossprod(svd0$u))

The singular values, svd0$d, are on the standard deviation scale. To convert them to standard deviations on the scale of the response, multiply by the estimate of $\sigma$.

zapsmall(getME(m0,"sigma")*svd0$d)

A more meaningful scaling, as used in PCA, is to consider the proportion of the overall variance accounted for by each component.

vc <- svd0$d^2   # variances of principal components
zapsmall(vc/sum(vc))

The principal components are ordered so that the first component accounts for the largest proportion of the variance, the second component accounts for the second largest proportion, and so on. This ordering is enforced by the singular values which are defined to be non-negative values in decreasing order.

The cumulative proportion of the variance shows how what proportion is in the subspace spanned by the first principal component, the first two components, the first three components and so on, indicating how many components we should retain.

cumsum(vc/sum(vc))

We see that 100% of the variance is accounted for by the first three principal components, which is another way of saying that the covariance matrix is of rank 3. Furthermore, 96.7% of the variance in the random effects is in the first two prinipal components, indicating that it should be possible to reduce the model from three to two dimensions without too much of a drop in the quality of the fit.

Using the `rePCA` function

The rePCA (random-effects Principal Components Analysis) function takes a object of class lmerMod (i.e. a model fit by lmer) and performs the steps decribed above to produce a list of prcomp objects.

prc <- rePCA(m0)
class(prc)
length(prc)
names(prc)
class(prc$subj)
prc$subj
summary(prc$subj)

The prcomplist class has its own summary method, which simply applies summary to each element of the list.

Multiple random-effects terms for the same grouping factor

When several random-effects terms have the same grouping, the Cholesky factors from the terms are accumulated in a block-diagonal matrix. For example, in the zero correlation parameter model

VarCorr(m1 <- lmer(rt ~ 1+c1+c2+c3 + (1+c1+c2+c3||subj), KWDYZ, REML=FALSE))

the "Tlist" consists of four 1x1 matrices

getME(m1,"Tlist")

all named "subj". The block diagonal matrix, which in this case is simply a diagonal matrix, is created by

bdiag(getME(m1,"Tlist"))

The singular value decomposition of a diagonal matrix is trivial; d is the diagonal of the matrix and u is the identity.

svd(bdiag(getME(m1,"Tlist")),nv=0)

summary(rePCA(m1))

The block diagonal nature of the Cholesky factor is better illustrated with the results of the model

VarCorr(m2 <- lmer(rt ~ 1+c1+c2+c3 + (1+c1+c2|subj) + (0+c3|subj), KWDYZ, REML=FALSE))

(chf <- bdiag(getME(m2,"Tlist")))
svd(chf,nu=0,nv=0)
summary(rePCA(m2))

Two or more grouping factors

When a model incorporates random effects with repect to two or more grouping factors, such as this model fit to the kb07 data from @Kronmuller:Barr:2007 and also discussed in @Barr:Levy:Scheepers:Tily:13

m3 <- lmer(RTtrunc ~ 1+S+P+C+SP+SC+PC+SPC + (1|subj) + (1+P|item), kb07, REML=FALSE)
print(summary(m3),corr=FALSE)

the rePCA function produces a list of prcomp objects, one for each grouping factor.

summary(rePCA(m3))

Summary

Principal components analysis (PCA) of the estimated covariance matrices for the random effects in a linear mixed model allows for simple assessment of the dimensionality of the random effects distribution. As shown in other vignettes in this RePsychLing package, the maximal model in many analyses of data from Psychology and Linguistics experiments, is almost always shown by this analysis to be degenerate.