An 'average information' restricted maximum likelihood algorithm for estimating reduced rank genetic covariance matrices or covariance functions for animal models with equal design matrices

On decrit un algorithme de maximum de vraisemblance restreint de type quasiNewton, qui approche la matrice Hessienne par la moyenne de l'information observee et attendue, dans le but d'estimer les composantes de covariance ou les fonctions de covariance dans un modele lineaire mixte. La strategie de calcul envisagee repose sur les outils de traitement des matrices creuses et sur la differentiation automatique d'une matrice, sans necessiter d'inversions de grandes matrices creuses. Dans le cas particulier d'un modele avec un seul facteur aleatoire et une matrice d'incidence identique pour tous les caracteres, les calculs de la vraisemblance, de ses derivees premieres et secondes «moyennes» peuvent etre effectues caractere par caractere, ce qui ramene les besoins de calcul lies a une analyse multivariate au niveau de ceux d'une serie d'analyses univariates. Ceci est rendu possible par la decomposition canonique des matrices de covariance a partir de la transformation des donnees en caracteres nouveaux, non correles entre eux. Le rang de la matrice de covariance genetique estimee est determine par le nombre de valeurs propres non nulles de la decomposition canonique, et donc peut etre reduit quand on fixe a zero certaines valeurs propres. Le nombre d'analyses univariates est ainsi egal au rang. Ceci est particulierement utile pour l'estimation de la fonction de covariance, qui decrit les covariances entre un tres grand nombre de caracteres tres correles par l'intermediaire d'un polynome d'ordre inferieur.

(Received 21 May 1996; accepted 17 January 1997) Summary -A quasi-Newton restricted maximum likelihood algorithm that approximates the Hessian matrix with the average of observed and expected information is described for the estimation of covariance components or covariance functions under a linear mixed model. The computing strategy outlined relies on sparse matrix tools and automatic differentiation of a matrix, and does not require inversion of large, sparse matrices. For the special case of a model with only one random factor and equal design matrices for all traits, calculations to evaluate the likelihood, first and 'average' second derivatives can be carried out trait by trait, collapsing computational requirements of a multivariate analysis to those of a series of univariate analyses. This is facilitated by a canonical decomposition of the covariance matrices and corresponding transformation of the data to new, uncorrelated traits. The rank of the estimated genetic covariance is determined by the number of nonzero eigenvalues of the canonical decomposition, and thus can be reduced by fixing a number of eigenvalues at zero. This limits the number of univariate analyses needed to the required rank. It is particularly useful for the estimation of covariance function when a potentially large number of highly correlated traits can be described by a low order polynomial. REML / average information / covariance components / reduced rank / covariance function / equal design matrices INTRODUCTION Estimation of (co)variance components by restricted maximum likelihood (REML) fitting an animal model to date is mainly carried out using a derivative-free (DF) algorithm as initially proposed by Graser et al (1987). While this has been found to be slow to converge, especially for multi-trait and multi-parameter analyses, it does not require the inverse of a large matrix and can be implemented efficiently using sparse matrix storage and factorisation techniques, making it computationally feasible for models involving tens of thousands of animals.
Recently there has been renewed interest in algorithms utilising derivatives of the likelihood function to locate its maximum. This has been furthered by technical advances, making computations faster and allowing larger and larger matrices to be stored. Moreover, the rediscovery of Takahashi et al's (1973) algorithm to invert large sparse matrices has removed most of the constraints on algorithms imposed previously by the need to invert large matrices.
In particular, 'average information' (AI) REML, a quasi-Newton algorithm, which requires first derivatives of the likelihood but replaces second derivatives with the average of the observed and expected information, described by Johnson and Thompson (1995) has been found to be computationally highly advantageous over DF procedures.
It is well recognised that for several correlated traits, most information available is contained in a subset of the traits or linear combinations thereof. This subset is the smaller the higher the correlations between traits. More technically, several eigenvalues of the corresponding covariance matrix between traits are very small or zero. If a modified covariance matrix were obtained by setting all small eigenvalues to zero and backtransforming to the original scale (using the eigenvectors corresponding to non-zero eigenvalues), it would have reduced rank.
There has been interest in reduced rank covariance matrices in several areas.
Wiggans et al (1995; unpublished) collapsed the multivariate genetic evaluation for 30 traits (ten test day records each for milk, fat and protein yield in dairy cows) to the equivalent of five univariate analyses by reducing the rank of the genetic covariance matrix and exploiting a transformation to canonical scale. Kirkpatrick and Heckman (1989) introduced the concept of 'covariance functions', expressing the covariance between traits as a higher order polynomial function. Polynomials can be fitted to full or reduced order. In the latter case, the resulting covariance matrix has reduced rank, ie, a number of zero eigenvalues (Kirkpatrick et al, 1990). The covariance function (CF) model was developed with the analysis of 'traits' with potentially infinitely many repeated, or almost repeated records in mind, where the phenotype or genotype of individuals is described by a function rather than a finite number of measurements (Kirkpatrick and Heckman, 1989 While it is possible to modify an estimated covariance matrix to reduce its rank (as done by Kirkpatrick et al, 1990Kirkpatrick et al, , 1994, it would be preferable to impose restrictions on the rank of covariance matrices 'directly' during (REML) estimation. Ideally, this could be achieved by increasing the order of fit (ie, rank allowed) sequentially until an additional non-zero eigenvalue does not significantly increase the likelihood.
Conceptually, this could be implemented simply by reparameterising, to the eigenvalues and corresponding eigenvectors of a covariance matrix, and fixing the required number of eigenvalues at zero. Practical applications of such reparameterisations, however, have been restricted to simple animal models with equal design matrices for all traits; see Jensen and Mao (1988) for a review. For these, a canonical decomposition of the genetic and residual covariance matrix together yields a transformation to uncorrelated variables with unit residual variance, leaving the number of parameters to be estimated unchanged (for full rank). Meyer and Hill (1997) described how REML estimates of CFs or, more precisely, their coefficients could be obtained using a DF algorithm through a simple reparameterisation of the variance component model. However, they found it slow to converge for orders of fit greater than three or four. Moreover, for simulated data sets the DF algorithm failed to locate the maximum of the likelihood accurately in several instances, especially if CFs were fitted to a higher order than simulated. This paper reviews an AI-REML algorithm for the general, multivariate case, presenting a computing strategy that does not require sparse matrix inversion. Subsequently, simplifications for the special case of a simple animal model with equal design matrices for all traits are considered. Additional reductions in computational requirements are shown for the estimation of reduced rank genetic covariance matrices or reduced order CFs.

Model of analysis
Consider the multivariate linear mixed model for t traits with y, 13, u and e denoting the vector of observations, fixed effects, random effects and residual errors, respectively, and X and Z are the incidence matrices pertaining to (3 and u. Let V(u) = G, V(e) = R and Cov(u,e') = 0, so that V(y) = V = ZGZ' + R. For an animal model, u always includes the vector of animals' additive genetic effects (a). In addition, it may contain other random effects, such as animals' maternal genetic effects, permanent environmental effects due to the animal or its dam, or common environmental effects such as litter effects.
Let E A = {a Ai j} ' denote the t x t matrix of additive genetic covariances. For u = a this gives G = E A &reg; A where A is the numerator relationship matrix and 0 denotes the direct matrix product. If other random effects are fitted, G is expanded correspondingly; see Meyer (1991) for a more detailed description. Assuming y is ordered according to traits within animals where N is the number of animals that have records, and 2:: + denotes the direct matrix sum (Searle, 1982). Let E E = {!E!! be the matrix of residual covariances between traits. For t traits, there are a total of W = 2t -1 possible combinations of traits recorded (assuming single records per trait), eg, W = 3 for t = 2. For animal i with combination of traits w, R i is equal to !Ew' the submatrix of E E obtained by deleting rows and columns pertaining to missing records.

Average information REML
Assuming a multivariate normal distribution, ie, y N N(Xb, V), the log of the REML likelihood (G) is (eg, Harville, 1977) where X * denotes a full-rank submatrix of X, and Let e denote the vector of parameters to be estimated with elements O i for i = 1, ... , p. Derivatives of log G are then (Harville, 1977) The latter is commonly called the observed information. It has expectation For V linear in 9, a2V/8Bi8B! = 0, and the average of observed [5] and expected [6] information is (Johnson and Thompson, 1995) The right hand side of [7] is (except for a scale factor) equal to the second derivative of y'Py with respect to B i and 0j , ie, the average information is equal to the data part of the observed information.
REML estimates of e can then be obtained by substituting the average information matrix for the Hessian matrix in a suitable optimisation scheme which uses information from second derivatives of the function to be maximised; see Meyer and Smith (1996) for a detailed discussion of Newton-Raphson-type algorithms in this context.

Calculation of the log likelihood
Calculation of log G pertaining to [1] has been described in detail by Meyer (1991). It relies on rewriting [3] as (Graser et al, 1987;Meyer, 1989) where C is the coefficient matrix in the mixed model equations (MME) for [1] (or a full rank submatrix thereof).
The first two components of log L can usually be evaluated indirectly, requiring only the log determinants of matrices of size equal to the maximum number of records or effects fitted per animal. For u = a where N A denotes the number of animals in the analysis (including parents without records). log IAI is a constant and can be omitted for the purpose of maximising log ,C. Similarly, with N w denoting the number of animals having records for combination of traits w The other two terms in [8], log ICI and y'Py, can be determined in a general way for all models of form !1!. Let M (of size M x M) denote the mixed model matrix (MMM), ie, the coefficient matrix in the MME augmented by the vector of right hand sides (r) and a quadratic in the data vector A Cholesky decomposition of M gives M = LL', with L a lower triangular matrix with elements l ij (l ij = 0 for j > i), and Factorisation of M for large scale animal model analyses is computationally feasible through the use of sparse matrix techniques; see, for instance, George and Liu (1981 (Meyer and Smith, 1996). Let D! = å'5:. A/ å() i be a matrix whose elements are 1 if B i is equal to the klth element of E A and zero otherwise. Further, let 6 ki denote Kronecker's Delta, ie, 6 kl = 1 for k = and zero otherwise, and QA l denotes the klth element of EA 1 . For oj = a A kl Similarly, with D! = BE Ew/ 8Bi and a'lL the klth element of Y. l ¿ while all other first derivatives of log !G! and log !R! are zero. Smith (1995) describes a procedure for automatic differentiation of the Cholesky decomposition. In essence, it is an extension of the Cholesky factorisation which gives not only the Cholesky factor of a matrix but also its derivatives, provided the corresponding derivatives of the original matrix can be specified. In particular, Smith (1995) outlines a 'backwards differentiation' scheme that is applicable when we want to evaluate a scalar function of L, f (L).
It involves computation of a lower triangular matrix F. This is initialised to 10f(L)Ial ijl . On completion of the backwards differentiation, F contains the derivatives of f(L) with respect to the elements of M. Smith (1995) states that the calculation of F (not including the work needed to compute L) requires about twice as much work as one likelihood evaluation. Once F has been determined first derivatives of f (L) can be obtained one at a time as tr(F<9M/<9!t), ie, only one matrix F is required.
Meyer and Smith (1996) describe a REML algorithm utilising this technique to determine first and (observed) second derivatives of log G for the case considered here. For f (L) = log I C + y'Py, the scalar is a function of the diagonal elements of L (see [12] and [13]). Hence, {8f(L)/8l ij } is a diagonal matrix with elements n¡¡ 1 for i = 1, ... , M -1 and 21,!,1,! in row M.
The non-zero derivatives of M have the same structure as the corresponding part (data versus pedigree) part of M As outlined above, R is blockdiagonal for animals. Hence, matrices aR -1 / alJ Ekl have submatrices -E-'Do!E-1 ie, derivatives of M with respect to residual (co)variances can be set up in the same way as the 'data part' of M.
The strategy outlined for the calculation of first derivatives of log L does not require the inverse of the coefficient matrix C. In contrast, Johnson and Thompson (1994,1995) Meyer (1985) described a method of scoring algorithm for this case, exploiting a canonical transformation to reduce a t-variate analysis to t corresponding univariate analyses.
For E E positive definite and E A positive semi-definite, there exists a matrix Q such that A is a diagonal matrix with elements Àii ! 0 which are the eigenvalues of !E/!A, and S2 = It (eg, Graybill, 1969) Transforming the data to then yields t new, 'canonical' traits which are uncorrelated and have unit residual variance. This makes the corresponding coefficient matrix in the MME blockdiagonal for traits, ie, Meyer (1991) described how the log likelihood (on the original scale) in this case can be computed trait by trait as the sum of univariate likelihoods on the canonical scale plus an adjustment for the transformation (last term in !29!) with y* the subvector of y * for trait i and P* the ith diagonal block of the projection matrix on the canonical scale P * which, like C * , is blockdiagonal for traits.
Terms required in [29] can be calculated by setting up and factoring, as described above, univariate MMM (on the canonical scale), Mi , of size M o = (M -1)/t + 1 each Moreover, all first derivatives of log G as well the average information matrix, both on the canonical scale, can be determined trait by trait.
First derivatives on the canonical scale Consider the parameterisation of Meyer (1985) where 0 * , the vector of parameters on the canonical scale, has elements Ag and wi! for i < j = 1, ... , t, ie, parameters are the (co)variances on the canonical scale.
The log likelihood on the canonical scale can be accumulated trait by trait, because Cholesky decompositions of individual MMM, M!, yield the submatrices and subvectors for trait i which are obtained when decomposing M * = L * L * ', ie, On the original scale, L and F have the same sparsity structure (Smith, 1995). However, while Ay = Wij = 0 for given E A and E E , the corresponding derivatives and estimates are not, unless the maximum of the likelihood has been attained.
Hence, while the off-diagonal blocks of LC are zero, the corresponding blocks of F * = 8f (L * )/c7M * are not.
It can be shown that both the diagonal blocks of F * corresponding to L* F* and the row vectors corresponding to 1*', fi are identical to those obtained by backwards differentiation of L!. In other words, first derivatives with respect to the variance components on the canonical scale (!ii and Wii ) can be obtained trait by trait from univariate analyses. Calculation of derivatives with respect to !2! and Wij , however, requires the off-diagonal blocks of F * corresponding to traits i and j, FC_! . Fortunately, as outlined in the A PP endix, matrices FC.! can be determined indirectly from terms arising using the Cholesky decomposition and backwards differentiation for individual traits on the canonical scale. Other terms required to determine the first derivatives on the canonical scale are where G * and R * are the canonical scale equivalents to G and R, respectively.
Average information on the canonical scale For G * = A &reg; A, R * = I tN , and thus V * = Var(y * ) blockdiagonal for traits, [20] and [21] simplify to ie, vectors bi are zero except for subvectors for traits k and l, bi k = si and bi l = sk , with Sj standing in turn for A!Zoa! and 8) .
With P * blockdiagonal for traits, 0g = ).,kl or u ki (k x l) and 0) = A mn or um n (m < n), this gives Hence, calculation of the average information on the canonical scale requires all terms s'P*s j for i ! j = 1,..., 2t and k = 1, t. These can be obtained trait by trait, analogously to the procedure described above for the general case. After the Cholesky factorisation of Mk has been carried out, solution a* for animals' additive genetic effects and residuals ek for trait k are obtained, storing the Cholesky factor L*. Define a matrix S of size N x 2t with columns equal to vectors s i . S is the canonical scale equivalent to B above. Once all columns of S have been evaluated, set up a matrix vi 1for each trait. This is the MMM for trait k on the canonical scale with yk replaced by S. The matrix S'S has elements s's j which are the sum of squares and crossproducts of the vectors of (weighted) solutions and residuals. Absorbing Ck into S'S (using the stored matrix L*) then overwrites this matrix with S'PkS which has elements S!PkSj' After all t traits have been processed bi!P*b!, (twice) the average information, can be 'assembled' according to !38!.
Derivatives on the original scale Let H * , with elements 2 bi ! P* b! , and g * , with elements 810g £* /80g, denote the average information matrix and vector of gradients on the canonical scale, respectively. Corresponding terms on the original scale are then where J with elements <9!/c)! is the Jacobian matrix of e with respect to e * . From [25] and [26], J has non-zero elements A numerical example illustrating calculations is given in the Appendix.

Alternative parameterisation
The above parameterisation requires switching between the original and canonical scale (or accumulating the transformations). Alternatively, as performed by Meyer (1991) for a derivative-free algorithm, 0 * can be defined to have t elements A! (i = 1, ... , t) and t 2 elements q2! (i, j = 1, ... , t), ie, the 'genetic' variances on the canonical scale and the elements of the transformation matrix Q. This allows estimation to be carried out on the canonical scale. Moreover, as outlined by Meyer (1991) for a derivative-free algorithm, evaluation of log G for given Aii requires scalar calculations involving the elements of Q only, ie, maximising the conditional log,C (for given Àii ) with respect to q ij is computationally inexpensive. This allowed the maximum of log G to be located in a two-step procedure with computational requirements equivalent to those of a t parameter (rather than t(t + 1) ) analysis. While derivatives with respect to q2! (not shown) require the off-diagonal blocks of F, derivatives with respect to A, do not (see !32!). Hence a nested maximisation, using an average information step to estimate parameters A zz and a derivative-free procedure to maximize G with respect to q2! for given A zz appears computationally advantageous in this case.

REDUCED RANK COVARIANCE MATRICES
Forcing the estimate of the genetic covariance matrix E A to have reduced rank (k A < t), results in a number of zero eigenvalues on the canonical scale. In this case, several terms contributing to log G or its derivatives are constant or depend on the canonical transformation (Q) only. Thus they need to be evaluated only once per analysis, effectively reducing the computational requirements per round of iteration to those for univariate analyses for the k A non-zero eigenvalues.
Likelihood: for Aii ! 0, contributions to log G [29] become While the former is a constant, determined by the data structure only, the latter depends on the canonical transformation. However, as [45] shows, it can be calculated for any Q from the corresponding residual sums of squares (SS) and crossproducts (CP) on the original scale. Let Y denote a matrix of size N x t with columns equal to vectors of observations y 2 . Both log IX'X O and quadratics ym (I -Xo(XoXo)-Xo) y n can then be determined by factoring the matrix First derivatives: to estimate genetic and residual covariance components, only the derivatives on the canonical scale (with respect to A ij or wi! ) are required for which the eigenvalues for both canonical traits (A ii and A jj ) are non-zero. This reduces the number of Cholesky decompositions (of matrices MZ ) and corresponding backwards differentiations (to obtain matrices F* and vectors f2 ) to be carried out to the number of non-zero eigenvalues, k A . Moreover, only k A (k A -1)/2 instead of t(t-1)/2 off-diagonal blocks Fc * .. 'J need to be evaluated. This can reduce computational requirements per round of iteration dramatically.
Average information: as shown above, the average information matrix can be constructed from residual SS and CP in the vectors of random effects solutions and residuals. For A kk = 0, ak = 0. Residuals on the canonical scale are then a linear combination of residuals on the original scale, ek = !m-1 q km8 m with 8 m = y&dquo;, -Xo(3&dquo;!. Again the latter need to be evaluated only once. Terms required for the AI can then be obtained as before, replacing Ms k in [39] with if A kk = 0. For multiple eigenvalues equal to zero, [47] only needs to be evaluated once (per iterate).

COVARIANCE FUNCTIONS
The 'infinitely-dimensional' model Consider 'repeated' measurements taken at a number of ages (or equivalent), with potentially infinitely many records. The covariance between records taken at ages l and m can be expressed as where 7 is the covariance function and K with elements K ij is the pertaining matrix of coefficients, and a&dquo;, is the mth age, standardised to the interval for which the polynomials are defined. 4! is a matrix of orthogonal polynomial functions evaluated at the given ages with elements !2! _ <!(0t), the jth polynomial evaluated for age i, and k denotes the order of polynomial fit. Conceptually k = oo, but in practice k < t for observations at t ages. Kirkpatrick et al (1990Kirkpatrick et al ( , 1994 suggested the use of the so-called Legendre polynomials (see Abramowitz and Stegun, 1965) which span the interval from -1 to 1.
Note that polynomials include a scalar term, ie, for t records a full order fit involves terms to the power 0, ... , t -1. From [48], a covariance matrix can be rewritten as For a reduced order fit, k < t, 4) is rectangular with k columns and only k(k + 1)/2 coefficients K2! need to be estimated. Let E A = 4l A K A4 l Q for an order of fit k A , and K A the matrix of coefficients for the corresponding covariance function A. Further, partition the error covariance matrix into components owing to permanent and temporary environmental effects, E E = E R + E,. Under the 'finite model', these usually cannot be disentangled unless there are repeated records for the same age. Assume the latter represent independent 'measurement errors', ie, E, = Diag {&OElig;;i } ' and that the former can be described by a CF, R, ie, E R = 4l R K R 4l £ with order of fit k R . Fitting measurement errors separately together with R to the order t -1 yields an equivalent model to a full order fit for E E . Hence the maximum for k R is t -1 rather than t.

General case
Estimates of the elements of the coefficient matrices of the covariance functions (and the measurement errors) can be estimated by REML using algorithms for the multivariate estimation of covariance components together with a simple reparameterisation.
Likelihood: as outlined by Meyer and Hill (1997), log L [9] can be rewritten as Under the CF model, the vector of parameters to be estimated is r l with elements K A ,, for i # j = 1, ... k A , KR!! for i ! j = 1, ... k R and & O E l i g ; ; i for i = 1,..., t. For r CFs to be estimated, it has minimum length r + t (fitting all CFs to order 1 and assuming all & O E l i g ; ; i to be distinct) and maximum length rt(t+1)/2 (fitting all CFs to full order).
for 77 i = <7! and 0 j = aEmm' Note that a reduced order fit for ,A (k A < t) implies a genetic covariance matrix E A of reduced rank. This may lead to computational problems when applying 'standard' methodology to factor the MMM (such as the Cholesky decomposition and its automatic differentiation) as described above. For practical applications, this can be overcome by setting the t -k A zero eigenvalues to an operational zero, ie, a small non-zero value such as 10-5 or 10-6 .
As for the estimation of covariance components, extensions to other models including additional random effects are straightforward.

Equal design matrices
Assuming measurement errors are greater than zero, E E = E R + E E is positive definite, regardless of the order of fit for R. Hence, the transformation to canonical scale exists and the log likelihood can be determined as for the 'finite model', summing terms from univariate analyses on the canonical scale as described above.
For the estimation of covariance components, the canonical transformation doubled as a tool to reduce computational requirements (allowing calculations to be carried out for one trait at a time) and as a reparameterisation, ie, the likelihood was maximised with respect to the eigenvalues and elements of eigenvectors of the canonical decomposition of E A and E E . For the CF model there are three matrices to be considered, K A and K R , or E A and E R derived from them, and E E . Decompositions diagonalising several matrices exist. Lin and Smith (1990), for instance, used the common principal components algorithm of Flury and Constantine (1985) as an equivalent to the canonical decomposition for a multivariate mixed model with several random effects. However, this required the matrices to be transformed to be positive definite, and is thus not suitable for this application. Hence only the first parameterisation described above, namely 0 * having elements A ij and Wij for i < j = 1, ... , t is suitable for the estimation of CFs. First derivatives and average information on the canonical scale can then be determined and transformed back to the original scale as described above. The Jacobian J in this case has non-zero elements where 7 / Jm n are the elements of O = Q8. DISCUSSION A strategy has been outlined for computing the log likelihood in a multivariate mixed model, together with its first and 'average' second derivatives with respect to covariance components or the coefficients of covariance functions. These can be used in a (modified) Newton-type estimation procedure (Marquardt, 1963), together with an additional transformation to ensure estimates to be within the parameter space; see Meyer and Smith (1996) for details.
For a model with only one random factor and equal design matrices for all traits, calculations can be carried out trait by trait, reducing computational requirements to those of a series of univariate analyses. This is feasible through a canonical decomposition of the genetic and residual covariance matrices and a corresponding, linear transformation of the data to new, uncorrelated variables.
In that case, fitting a covariance function to less than full order or forcing estimated genetic covariance matrices to be of reduced rank is equivalent to fixing eigenvalues on the canonical scale at zero. Most terms pertaining to zero eigenvalues then only need to be calculated once per analysis, resulting in considerable reductions of computational requirements. In essence, it reduces work required in each round of an iterative solution scheme to that equivalent to k A univariate analyses. This is particularly useful for a comparatively large number of highly correlated measurements where a covariance function of low order suffices to describe the data adequately. In contrast, for analyses where a transformation to canonical scale is not feasible, the complete t-variate MMM needs to be set up, factored and differentiated in each round of iteration, even if CFs are only fitted to order k A .
Making computational demands proportional to the order of fit for the genetic CF (or matrix) encourages an 'upwards' strategy: increasing the order of fit one by one allows a likelihood ratio test to be performed at each step. The minimum number of parameters describing the data is found when the likelihood ceases to increase significantly when the order of fit is increased. It is envisioned that estimation of reduced order covariance matrices or functions will become the standard procedure for high-dimensional multivariate analyses.