Consider a multivariate linear mixed model for *q* traits:

$$\mathbf{y}= \mathbf{X} {\boldsymbol{\upbeta}} + \mathbf{Zu} + \mathbf{e} ,$$

(1)

with **y**, \({\varvec{\upbeta}}\), **u** and **e** denoting the vectors of observations, fixed effects, animals’ additive genetic effects and residuals and **X** and **Z** the pertaining design matrices, where **y** is assumed to have a multivariate normal distribution and **X** is assumed to have full rank. Let both **y** and **u** be ordered by traits within individuals. This gives \({\text{Var}}(\mathbf{e}) = \mathbf{R}\) which is block-diagonal for animals with blocks equal to submatrices of \({\varvec{\Sigma}_{e}}\) corresponding to the (subset of) traits recorded, with \({\varvec{\Sigma}_{e}}\) the \(q \times q\) matrix of residual covariances between traits. Furthermore, \({\text{Var}}\left(\mathbf{u} \right) = \mathbf{H} \otimes {\varvec{\Sigma}_{u}}\), with **H** the relationship matrix between animals, \({\varvec{\Sigma}_{u}}\) the \(q \times q\) genetic covariance matrix between traits and ‘\(\otimes\)’ denoting the Kronecker product. Among other options, **H** can represent the pedigree-based numerator relationship matrix, the genomic relationship matrix or the joint relationship matrix between genotyped and non-genotyped individuals [1].

The pertaining mixed model equations (MME) equations are then:

$$\left[ \begin{array}{cc} \mathbf{X}^{\prime} \mathbf{R}^{-1}\mathbf{X} & \mathbf{X}^{\prime} \mathbf{R}^{-1}\mathbf{Z}\\ \mathbf{Z}^{\prime} \mathbf{R}^{-1}\mathbf{X} & \mathbf{Z}^{\prime} \mathbf{R}^{-1}\mathbf{Z} + \mathbf{H}^{-1}\otimes {\varvec{\Sigma}}_{u}^{-1} \end{array} \right] \left[\begin{array}{ll} \hat{{\varvec{\upbeta}}} \\ \hat{\mathbf{u}} \end{array} \right] = \left[ \begin{array} {ll} \mathbf{X}^{\prime} \mathbf{R}^{-1}\mathbf{y} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1}\mathbf{y} \end{array} \right].$$

(2)

Let **C**, partitioned as in Eq. (2), denote the coefficient matrix in the MME:

$$\mathbf{C} = \left[\begin{array} {cc} \mathbf{X}^{\prime} \mathbf{R}^{-1}\mathbf{X} & \mathbf{X}^{\prime} \mathbf{R}^{-1}\mathbf{Z}\\ \mathbf{Z}^{\prime} \mathbf{R}^{-1}\mathbf{X} & \mathbf{Z}^{\prime} \mathbf{R}^{-1}\mathbf{Z} + \mathbf{H}^{-1}\otimes {\varvec{\Sigma}}_{u}^{-1} \end{array} \right]= \left[\begin{array}{cc} \mathbf{C}_{\beta \beta} & \mathbf{C}_{\beta u} \\ \mathbf{C}_{\beta u}^{\prime} & \mathbf{C}_{uu}\end{array} \right].$$

(3)

Each of the four submatrices has non-zero elements contributed by the data, i.e. arising from (weighted) products of the design matrices **X** and **Z**. For \(\mathbf{C}_{uu}\), \(\mathbf{Z}^{\prime} \mathbf{R}^{-1}\mathbf{Z}\) is sparse, consisting of diagonal blocks of size \(q \times q\) for *N* individuals. For traits that are not recorded, corresponding rows and columns (of the diagonal blocks) have zero elements. The second part of \(\mathbf{C}_{uu}\), \(\mathbf{H}^{-1} \otimes {\varvec{\Sigma}}_{u}^{-1}\), contributes a \(q \times q\) block for each non-zero element of \(\mathbf{H}^{-1}\).

An alternative formulation is obtained by rewriting Eq. (1) as [10, 15]:

$$\mathbf{y}= \mathbf{X} {\varvec{\upbeta}} + \mathbf{Z} \left ( \mathbf{I}_N \otimes \mathbf{Q} \right ) \left ( \mathbf{I}_N \otimes \mathbf{Q}^{-1} \right ) \mathbf{u} + \mathbf{e}= \mathbf{X} {\varvec{\upbeta}} + \mathbf{Z}^{\star}\mathbf{u}^{\star}+ \mathbf{e} ,$$

(4)

with \(\mathbf{I}_N\) denoting an identity matrix of size *N*, \(\mathbf{Z}^{\star}= \mathbf{Z} \left (\mathbf{I}_N \otimes \mathbf{Q}\right )\) and \(\mathbf{u}^{\star}= \left ( \mathbf{I}_N \otimes \mathbf{Q}^{-1} \right ) \mathbf{u}\). This gives \({\text{Var}}\left( \mathbf{u}^{\star}\right) = \mathbf{H} \otimes \mathbf{Q}^{-1}{\varvec{\Sigma}_{u}}\mathbf{Q}^{-T}\) and it follows that \(\mathbf{Z} {\text{Var}}\left( \mathbf{u} \right) \mathbf{Z}^{\prime} = \mathbf{Z}^{\star}{\text{Var}}\left( \mathbf{u}^{\star}\right) \mathbf{Z}^{\star}{^{\prime}}\). For computational efficiency, we can choose **Q** so that \(\mathbf{Q}^{-1}{\varvec{\Sigma}_{u}}\mathbf{Q}^{-T} = \mathbf{I}_q\). This is an equivalent model to Eq. (1) with the coefficient matrix in the pertaining MME:

$$\mathbf{C}^{\star}= \left[\begin{array}{cc} \mathbf{X}^{\prime} \mathbf{R}^{-1}\mathbf{X} & \mathbf{X}^{\prime} \mathbf{R}^{-1}\mathbf{Z}^{\star}\\ \mathbf{Z}^{\star}{^{\prime}} \mathbf{R}^{-1}\mathbf{X} & \mathbf{Z}^{\star}{^{\prime}} \mathbf{R}^{-1}\mathbf{Z}^{\star}+ \mathbf{H}^{-1}\otimes \mathbf{I}_q \end{array}\right] = \left[\begin{array}{ll} \mathbf{C}_{\beta \beta} & \mathbf{C}_{\beta u}^{\star} \\ {\mathbf{C}_{\beta u}^{\star}}^{\prime} & \mathbf{C}_{uu}^{\star}\end{array}\right] .$$

(5)

The transformation can be thought of as transferring genetic links between traits from the part due to the covariance of the random effect to the ‘data part’ of the coefficient matrix. For single records per trait, each row of **Z** typically has a single non-zero element of unity. In contrast, \(\mathbf{Z}^{\star}\) contains up to *q* non-zero coefficients per row, positioned in the columns representing the respective animals’ genetic effects for the *q* traits. Thus, the first part of \(\mathbf{C}_{uu}^{\star}\) is again block-diagonal for individuals, though — compared to \(\mathbf{C}_{uu}\) – some additional non-zero elements arise for missing records. Similarly, some additional non-zero elements can be generated in \(\mathbf{C}_{\beta u}^{\star}\) compared to \(\mathbf{C}_{\beta u}\). The second component (\(\mathbf{H}^{-1}\otimes \mathbf{I}_q\)) of \(\mathbf{C}_{uu}^{\star}\), however, is substantially sparser than its counterpart in \(\mathbf{C}_{uu}\), the more so the higher *q* and the denser \(\mathbf{H}^{-1}\). Here each non-zero off-diagonal element of \(\mathbf{H}^{-1}\) contributes only *q* additional non-zero coefficients in Eq. (5) compared to \(q^2\) in the above Eq. (3).

Suitable choices for **Q** arise from the eigen-decomposition of the genetic covariance matrix, \({\varvec{\Sigma}_{u}} = \mathbf{E} {\boldsymbol{\Lambda}} \mathbf{E}^{\prime}\). For **Q** equal to the matrix of eigenvectors, **E**, \(\mathbf{u}^{\star}\) would be equal to the vector of genetic principal component scores, which can be scaled to variances of unity by setting \(\mathbf{Q} = \mathbf{E} {\boldsymbol{\Lambda}}^{1/2}\) (with \({\boldsymbol{\Lambda}}\) the diagonal matrix of eigenvalues). Thus, we refer to Eq. (4) as the ‘principal components’ model. In general, this form of **Q** has \(q^2\) non-zero elements. This can be reduced to \(q(q+1)/2\) by an orthogonal transformation to lower triangular form or, more generally, by setting \(\mathbf{Q}\) to be a Cholesky factor of \({\varvec{\Sigma}_{u}}\). Furthermore, the latter integrates well with REML implementations which often parameterise to estimate the elements of the Cholesky factor of covariance matrices rather than the covariance components directly.

REML estimation based on the transformation to principal components is implemented in our mixed model package \(\mathbb {WOMBAT}\) [16, 17] for selected models of analysis.