Plant breeding selection tools built on factor analytic mixed models for multi-environment trial data

An advanced and widely used method of analysis for multi-environment trial data involves a linear mixed model with factor analytic (FA) variance structures for the variety by environment effects. This model can accommodate unbalanced data, that is, not all varieties in all environments, it allows the use of pedigree information and appropriate accommodation of individual trial experimental designs, and most importantly the FA structure for the variety by environment effects is parsimonious and regularly results in a good fit to the data. The model provides accurate predictions of the variety effects for every environment in the data-set but this constitutes a large and unwieldly amount of information to process for the purpose of variety selection. We address this issue in the current paper by proposing factor analytic selection tools to summarise the predictions in a concise yet informative manner. The tools, which are natural derivatives of the FA structure, result in measures of overall performance and stability across the environments in the data-set. All measures are expressed on the same scale as the trait under consideration and can easily be combined to form an index for selection.


Introduction
Plant breeders use information from the analysis of multi-environment trial (MET) data to select superior varieties. A commonly occurring impediment to selection is the presence of variety by environment interaction (VEI), which is characterised by the differential response of varieties to a change in environment. A statistical analysis of MET data must therefore aim to accurately encapsulate VEI. To highlight some of the difficulties in analysing MET data, and the short-comings of current approaches we consider an example from a tree breeding programme.
The Radiata Pine Breeding Company (RPBC) conducts genetic trials in which progeny trees obtained from selected parental trees are grown in test plantations. The MET data-set under study in this paper comprises a series of 92 genetic trials planted in different years and in various geographic locations in New Zealand and Australia. This is an extension of the data-set used in Cullis et al. (2014) and includes 15 more recently planted trials. The trait of interest is tree stem diameter (cm) measured at breast height (DBH). Each trial comprised a number of plots (areas of land measuring approximately 3 m 9 3 m) to which progeny trees were allocated. In the majority of trials there was no true replication of progeny trees so that each plot was allocated a genetically distinct tree. In eight trials, henceforth called clonal trials, the progeny trees were cloned so that there were multiple plots (true replicates) of each progeny tree. The experimental designs varied between trials and the possible types of designs are described in Cullis et al. (2014). The full MET data-set comprised 348,806 plots (and thence data records) corresponding to 336,528 trees, with single data records for each of the 334,977 nonclonal progeny trees and multiple data records (replicates) for each of the 1551 clonal progeny trees. Pedigree information was available on 339,589 trees which comprised the 336,528 progeny trees that were grown in the trials and 3061 parental trees. Note that none of the parental trees were grown in the trials.
The aim of the analysis was to obtain predicted additive genetic effects (also known as estimated breeding values, EBVs) for use not only within the breeding programme but also the New Zealand and Australian forest industry as a whole. The breeding programme uses EBVs to select parents for crossing, with the aim of producing genetically improved germplasm. At an industry level, forest owners purchase seed from nurseries and pay a premium for seed from superior parents. The quality of the parents is reflected in a rating (GF Plus 2006) that is derived from their EBVs. In the data-set under study, EBVs were required for the majority (3057) of the existing parents (so-called backward selections) and also the 1551 progeny trees in the clonal trials since they are potential new parents (forward selections). Here-after, these 4608 trees will be called ''varieties'' in order to align with the standard terminology of METs, in particular ''variety by environment interaction''. We stress that this is for pedagogical reasons only since Pinus radiata is an outcrossing species so these trees are not varieties in the sense of inbred crops. In terms of the trait of DBH, both the breeding programme and industry are primarily interested in varieties that will produce progeny that are likely to grow well across a wide range of environments as represented by the geographic locations and planting dates of the trials in the data-set.
This example illustrates some important requirements and difficulties that arise in practice. First is that pedigree information must be included in the analysis in order to investigate additive VEI and obtain EBVs.
Second is the potential computational burden associated with modelling VEI given the large numbers of trees and environments. Cullis et al. (2014) overcame this by proposing an approximate reduced animal (ARA) model which meant that additive VEI was modelled using only the 4608 varieties rather than the full set of 339,589 trees in the pedigree. This reduced the size of the problem but we note that there remained a substantial degree of imbalance in the data since not all of these varieties were represented in all trials. In fact, only 4% of all possible variety by environment combinations were represented. The inclusion of both clonal and non-clonal trials and the use of a range of experimental designs added to the complexity. The analysis must therefore be sophisticated enough to allow all of these issues to be adequately accommodated. Finally, it is noted that the breeding programme and industry require informative summaries of EBVs across environments in order to make selections and produce ratings.
It is instructive to consider any statistical analysis as comprising several distinct, but linked, components (Nelder 1994) which can be condensed into (1) model fitting and checking and (2) inference and prediction of effects of interest. In the context of MET data, component (1) requires the use of a model with appropriate genetic and non-genetic effects. The former must encapsulate VEI and the latter must reflect sources of variation and correlation associated with individual trial experimental designs. Component (2) requires summaries that provide breeders with concise and accurate information on which to base selection. In the presence of VEI this typically includes some measure of overall performance and ''stability'' across environments for each variety.
Many current approaches for the analysis of MET data focus on component (2) at the expense of component (1). So they aim to provide simple summaries by fitting simple models. We provide a brief history of these models here but the reader is referred to Smith et al. (2005) for a comprehensive review. Early methods of analysis for MET data involved an Analysis of Variance (ANOVA) of the two-way table of variety by environment means. The total variation in the data was partitioned into sources due to varieties, environments (trials) and residual variation which is a composite of VEI and within-trial error. Overall performance for a variety was obtained as the estimate of the variety main effect. Various 123 stability measures have been derived for this model (see Lin et al. 1986, for a review) and are typically some function of the residuals. A commonly used measure is the stability variance of Shukla (1972) which is the sample variance of the residuals for individual varieties.
A greater emphasis on interpreting VEI lead to the use of more complex models than ANOVA, in particular models involving regressions onto an independent environmental variable or onto the marginal (environment) means of the two-way table (Yates and Cochran 1938;Finlay and Wilkinson 1963). In these models the stability measure is the slope of the response for each variety. The motivation to examine more general patterns in VEI lead to the use of principal component analysis (PCA). Kempton (1984) used PCA on the residual effects from the two-way ANOVA model and displayed the results using biplots. This method of analysis was subsequently badged as AMMI (Additive Main effects and Multiplicative Interaction, Gauch 1992) and is still one of the most widely used methods for MET data. Note that the aim of this method is to be able to visualise relationships between varieties and environments. It does not provide simple numerical summaries that could be used for variety selection.
The models discussed thus far have numerous deficiencies. Some key issues are that they involve piecemeal approaches (typically first requiring analyses of individual trials to obtain variety by evironment means for use as data in a subsequent analysis) so are inherently inefficient (Welham et al. 2010;Gogel et al. 2018); most require balanced data, that is, all varieties in all environments; they assume variety effects to be fixed rather than random (see Smith et al. 2005, for a discussion), which has particular limitations for our example since it is not possible to include pedigree information and finally, they rarely provide a good fit to the data.
These methods therefore do not satisfactorily address the model fitting component of the MET analysis. In contrast, we consider the linear mixed model approach of Smith et al. (2001) and its extensions, in particular Oakey et al. (2007) and Beeck et al. (2010) who include pedigree information in order to partition genetic effects into additive and nonadditive effects and Cullis et al. (2014) who use a modification for estimating additive genetic effects in out-crossing plant species. Underpinning these approaches is a one-stage analysis of individual plot data combined across trials, and factor analytic (FA) structures for the variety effects in individual environments. Additionally separate non-genetic models are used for individual trials. The approach has been used in Australia for over 15 years and is now the preferred method of analysis in all major plant breeding programmes. The FA structure has been found to perform extremely well in terms of providing a good fit to the data and a parsimonious model for VEI (Kelly et al. 2007). The success of the approach for plant breeding programmes has also led to its adoption within the Australian National Variety Trials (NVT) system (Smith et al. 2015;Gogel et al. 2018).
In terms of component (2) of the analysis, various summaries from the FA model have been used to aid in examining VEI. These include heatmaps for visualising estimated genetic correlations between environments, and so-called latent regression plots which provide visual representations of the multiple regression implicit in the FA model (see Cullis et al. 2010Cullis et al. , 2014, for example). The FA model can also be regarded as a random effects analogue of the AMMI model so that bi-plots as per Kempton (1984) or Gauch (1992) can be constructed. Although these graphics are informative they do not directly address the fundamental issue of variety selection. The FA mixed model provides predictions of genetic effects for each variety in each environment in the data-set. These predictions provide a complete inventory of the twoway table of variety by environment effects which need to be summarised to facilitate selection. Until now there has been no statistically or biologically satisfactory way to achieve this without sacrificing information on VEI.
To summarise, the Smith et al. (2001) approach out-performs others in terms of the model fitting component of a MET analysis but it has failed to deliver on the prediction component, in the sense of providing concise information to aid with variety selection. We rectify this in the current paper by presenting factor analytic selection tools (FAST) which exploit the underlying form of the FA structure, in particular its analogy with multiple linear regression. FAST include natural measures of overall performance, stability and sensitivity for each variety which can then be used to form a relatively simple index for selection. We note that the measures accommodate selection for both broad adaptation and also for specific adaptation as guided by patterns of VEI revealed in the analysis itself. FAST can be routinely implemented in a plant breeding programme to identify both superior varieties and superior parents.
The paper is arranged as follows. In ''Statistical models'' the general linear mixed model for MET data with the inclusion of pedigree information is described, with particular attention given to the FA structures for variety by environment effects. The specific model for the motivating example is given. Methods for the investigation of VEI following the fitting of an FA structure are given in ''Post-processing to investigate variety by environment interaction''. This includes a detailed development of the new tools (FAST). The application of FAST to the motivating example and some remarks on the general applicability of the tools are given in ''Results and discussion''.

Statistical models
It is assumed that the MET data-set comprises t trials that have been conducted in p environments. Often, trials are synonymous with environments, so that p ¼ t. However there are situations in which there are multiple trials in an environment, so that p\t. Let y j denote the n j -vector of data for the jth trial. We then let y denote the n-vector of data combined across all trials in the MET, so write y ¼ ðy > 1 ; y > 2 ; . . .; y > t Þ > . Note that n ¼ P t j¼1 n j . The linear mixed model for y can be written as where s is a vector of fixed effects with associated design matrix X (assumed to have full column rank); u g is the vector of random genetic effects with associated design matrix Z g ; u p is a vector of random non-genetic (or peripheral) effects with associated design matrix Z p and e ¼ ðe > 1 ; e > 2 ; . . .; e > t Þ > is the combined vector of residuals from all trials. The vector of fixed effects includes mean parameters for individual environments. The vector of random peripheral effects includes effects associated with the experimental designs of individual trials.
The random genetic effects comprise the variety effects nested within environments, and will be referred to as the variety by environment (VE) effects. If we let m denote the total number of unique varieties across all environments, then the vector u g has length mp. Typically, not all varieties are grown in all environments so that the design matrix Z g will contain columns in which all the elements are zero. We assume the VE effects to be ordered as varieties within environments so they can be written as . . .; u g > p Þ > where u g j is the m-vector of VE effects for environment j.
The effects u g ; u p and e are assumed to be mutually independent, and distributed as multivariate Gaussian, with zero means. The variance matrix for u g will be described in detail below. The variance matrix for u p is typically given by is the number of components in u p and q i is the number of effects in (length of) u p i . The variance matrix for the residuals is assumed to be block diagonal, so that R is the variance matrix for the residuals for the jth trial.
In this paper we allow for the inclusion of pedigree information in the analysis, so partition the VE effects into additive and non-additive (residual VE) effects as follows: It is assumed that var u a ð Þ ¼ G a A where A is the numerator relationship matrix (see Oakey et al. 2007;Beeck et al. 2010, for example), and G a is a p Â p symmetric positive (semi)-definite matrix that will be referred to as the between environment additive genetic variance matrix. The diagonal elements are the variances of the additive VE effects for individual environments and the off-diagonal elements are the covariances between additive VE effects in different environments. In terms of the non-additive effects, it is assumed that var u e ð Þ ¼ G e I m where G e is a p Â p symmetric positive (semi)-definite matrix that will be referred to as the between environment non-additive genetic variance matrix. The variance matrix of the total VE effects (that is, additive plus non-additive) is therefore given by Note that if pedigree information is not included in the analysis, the VE effects simplify to u g ¼ u e with var u g À Á ¼ G e I m .

123
Variance models for genetic effects Following Smith et al. (2001) we propose a factor analytic model for the between environment additive genetic variance matrix. The aim is to account for the covariances of the additive VE effects between environments in terms of a small number, k a , of (unknown) common factors. The number, k a , is called the order of the model and we let FAk a denote the FA model with this order. The model is postulated in terms of the additive VE effects as linear combinations of the common factors, plus an error term. Thus the additive VE effect for variety i and environment j is written as Each of the first k a terms is the product of a variety effect (f a ri ), which is known as a score, and an environment effect (k a rj ), which is known as a loading. The final term, d a ij , represents the error or lack of fit in the model. This model can be written in vector notation as where k a r is the p-vector of environment loadings for the rth common factor and f a r is the associated mvector of variety scores; the mk a -vector of variety scores and d a ¼ . . .; d > a p Þ > is the mp-vector of lack of fit effects where d a j is the m-vector for the jth environment.
It is assumed that f a and d a are independent and are distributed as multivariate Gaussian with zero means and variance matrices given by where W a is a p Â p diagonal matrix with elements w aj , which are the so-called additive specific variances for individual environments. These assumptions lead to a factor analytic form for the between environment additive genetic variance matrix, namely . Note that the variance matrix for the additive VE effects for variety i is given by a ii G a , where a ii is the ith diagonal element of A and is defined to be 1 þ F i , where F i is the inbreeding coefficient for variety i. The FA model in Eq.
(3) can be naturally separated in two parts. We let b a ¼ K a I m ð Þ f a so that The effects in b a will be called the common additive VE effects because they are associated with the common factors that explain the additive genetic covariance between environments. To show this we partition b a conformably with u a so that b a j denotes the m-vector for the jth environment. Table 1 shows that the covariance between u aj and u ah is identical to, and therefore entirely defined by, the covariance between b a j and b a h . The common VE effects may also be thought of as the correlated VE effects since the effects for one environment are correlated with those in at least one other environment. In contrast, the lack of fit effects, d a j are uncorrelated between environments (see Table 1). In other words these effects are specific to the environment so that d a will be called the vector of specific additive VE effects. In a similar manner to the additive VE effects, a factor analytic model may also be assumed for the non-additive effects. The order of the model is denoted by k e and this may differ from k a . The FAk e model for the non-additive effects is then given by where the terms are defined in an analogous manner to those for the additive VE effects. It is assumed that f e and d e are independent and are distributed as multivariate Gaussian with zero means and variance matrices given by Model fitting, estimation and prediction The model fitting process commences with the estimation of the variance parameters using residual maximum likelihood (REML). Then given these estimates, empirical best linear unbiased estimates (EBLUEs) and empirical best linear unbiased predictions (EBLUPs) of the fixed and random effects, respectively, can be obtained. All models in this paper have been fitted using ASReml-R (Butler et al. 2009). In terms of the FA models, a sparse implementation of the average information algorithm (Thompson et al. 2003) is used so that reduced rank (RR) variance models are accommodated. In this paper we have exploited this by splitting the FA models into their component parts and explicitly fitting the common VE effects (with RR variance structure) separately from the specific VE effects. The associated variance parameters are the loadings and specific variances. The REML estimates of these parameters will be denoted byk s rj andŵ s j (for s ¼ a; e). It is important to note that when k s [ 1, the matrix of loadings is not unique so that constraints are required to ensure identifiability. The constraints used in ASReml-R (Butler et al. 2009) are to fix the elements in the upper triangle ofK s to zero.
EBLUEs and EBLUPs of the fixed and random effects are obtained as solutions to the mixed model equations (MME) (Henderson 1950). Given the sparse RR formulation of the FA model (Thompson et al. 2003), this provides EBLUPs of the common and specific VE effects and the variety scores, namelỹ b s ;d s andf s . Note that EBLUPs of the VE effects can then be obtained asũ s ¼b s þd s .
Note also, that if C is used to denote the coefficient matrix of the MME, then C À1 provides prediction error variances of effects. These can then be used to calculate a measure of accuracy for predictions of genetic effects. We letũ c ¼ ðũ > c a ;ũ > c e Þ > be the vector of EBLUPs of all genetic effects in the MME and wherẽ We consider scalar predictions of The accuracy of this prediction is given by where pevp ð Þ is the prediction error variance which is given by d > C cc d where C cc is the partition of C À1 that relates toũ c .

Statistical model for motivating example
The data from the motivating example were analysed using the approximate reduced animal (ARA) model of Cullis et al. (2014). In the current paper the analysis itself is not of primary interest, but rather the postprocessing of the results to facilitate selection. The interested reader is therefore referred to Cullis et al. (2014) for full details of the linear mixed model. The key features of the linear mixed model for the motivating example were the inclusion of a fixed main effect for each trial and random effects associated with experimental design terms. In accordance with the ARA model, additive VE effects were included for parents and clones (across all trials) and non-additive VE effects were included for clones (across clonal trials only). The former were modelled using a factor analytic model of order 3 (FA3). The variance matrix for the latter was assumed to have a diagonal form since the estimated non-additive genetic variances were found to be relatively small. A separate residual variance was fitted for each trial.

Post-processing to investigate variety by environment interaction
Varietal selection in METs is made more complex by the presence of VEI. VEI is characterised by the differential response of varieties to environments and 123 can be broadly categorised as either crossover or noncrossover. The former is often regarded as the most important for varietal selection since it is associated with changes in the rank of varieties between environments. In this section tools for exploring VEI will be developed in the context of the additive genetic effects. The methods are directly applicable to the non-additive effects but the total VE effects require special attention and will be discussed in ''FAST for total VE effects''. It is useful to write the additive VE effects as a twoway structure, namely the m Â p matrix inter-relationships between the columns (environments) and rows (varieties) of this matrix and may be considered from either perspective. In terms of variety selection, it is crucial to focus on VEI from the variety perspective. This can be done by exploring various aspects of the FA model, in particular the analogy with multiple regression. Given the FAk a model, the EBLUPs of the additive VE effects can be written as Equation (6) can then be viewed as a series of multiple regressions in which the independent variables are the estimated environment loadings,k a 1 . . .k a ka . There is a separate regression for each variety and the predicted regression coefficients are given by the predicted variety scores,f a 1 . . .f a ka :

Rotation of REML estimates of loadings
In later sections, individual terms in the regression will be explored. When k a [ 1, this first requires the independent variables, namely the estimated loadings, to be rotated to a meaningful solution. Recall from ''Model fitting, estimation and prediction'' that when k a [ 1, the loadings are estimated subject to constraints that are imposed for computational convenience. In order for the estimated loadings to have a meaningful interpretation we choose to rotate to a principal component solution so that the first rotated estimated loading accounts for the maximum amount of covariance in the VE effects, the second accounts for the next largest amount and is orthogonal to the first, and so on. The properties of orthogonality and decreasing contributions to covariance will be shown to be important for the selection tools developed in ''Factor analytic selection tools (FAST)''. Thus, after model fitting, we obtain the singular value decomposition ofK a : where L a is a diagonal matrix with elements given by the square roots of the eigenvalues ofK aK > a (arranged in decreasing order) and B a and V a are orthogonal matrices with columns given by the eigenvectors of K aK > a andK > aK a , respectively. Then we obtain the rotated estimated loadings aŝ where c is a constant that is either 1 or -1. The sign is chosen to ensure the majority of first rotated loadings are positive rather than negative. This aids with interpretation, particularly in the context of the selection tools developed in ''Factor analytic selection tools (FAST)''. After rotation it is meaningful to consider the percentage of additive genetic variance that is explained by individual factors. This can be computed for factor r and environment j as In this way the mean contribution of factor r to additive genetic variance can be computed across environments as P p j¼1 v a rj =p. Given the rotation of the estimated loadings, the predicted variety scores must also be rotated to: Note that the EBLUPs of both the VE effects and the common VE effects are invariant to the rotation.

Rotation of REML estimates of loadings for motivating example
In the motivating example k a ¼ 3 factors were fitted. The estimated loadings after rotation are given in the ''Appendix''. In summary, the first vector of loadings ranges from 1.0 to 19.1 with a mean of 10.0; the second vector ranges from -18.4 to 13.6 with a mean of 0.9; the third vector ranges from -11.8 to 14.7 with a mean of -0.4. The mean percentage variance accounted for by individual factors was 50.7, 14.8 and 13.7% for factors 1, 2 and 3, respectively, and 79.3% for the FA3 model as a whole (that is, the mean across all factors and environments). Note that the order of k a ¼ 3 was chosen using a pragmatic approach that aimed to balance parsimony and goodness of fit. The latter was assessed in terms of both the overall percentage variance accounted for and the distribution of individual environment values (also see Cullis et al. 2014). With respect to the latter, a large number (53) of the environments had a variance accounted for greater than 80%, whilst only a small number (14) had a variance accounted for less than 50%.

VE effects for individual environments
It would seem intuitive to use the EBLUPs,ũ a , of the additive VE effects to examine variety performance in individual environments. However, the EBLUPs,b a , of the common additive VE effects provide an alternative that have several appealing features. First, by definition, they represent VEI that is driven by influences that are common to several environments and therefore exclude isolated VEI that is specific to a single environment. They also have a natural interpretation as a set of ''smoothed'' VE effects. This is clear using the regression analogy of an FA model. In a similar manner to a standard multiple regression problem, predictions of the dependent variable in Eq. (6) may be obtained as fitted values along the regression surface. These are given by the sum of the first k a terms which is equivalent to the EBLUPs,b a , of the common VE effects. Finally we note that the EBLUPs of the common VE effects provides a set of predictions that are compatible across environments, irrespective of whether there is data on the variety in the environment. This is because they are fitted values on the regression surface. In contrast, the EBLUPs of the VE effects also include the EBLUPs of the specific VE effects which are the residuals in the regression. Predictions of the specific VE effects are intrinsically different for varieties with and without data in an environment. In the case of non-additive effects, the EBLUP of the specific VE effect will be zero when there is no data on the variety in that environment. Cullis et al. (2010) provide a more theoretical discussion of this issue and recommend the use of predictions that are ''marginal'' to the specific VE effects. Thus in the remainder of this paper it will be assumed that the EBLUPs of the common VE effects will be used for selections. Generalisations for the use of VE effects are straight-forward.
Note that all predictions can be accompanied by a measure of accuracy as detailed in ''Model fitting, estimation and prediction'' so that, in particular, the accuracy of VE predictions with and without data can be assessed.

Factor analytic selection tools (FAST)
From a breeding perspective, varietal selection using the full matrix of predicted common VE effects may be a formidable task unless the number of environments is very small. Broadly speaking, there is a need to obtain measures of overall performance and stability across environments for each variety. Cullis et al. (2014) suggested the use of latent regression plots to examine variety stability. These are similar to added variable plots but exploit the orthogonality property of the rotated factors. They comprise a series of k a plots for each variety in which the y -and x-axes for variety i are defined by Plot 1: y j ¼b a ij and x j ¼k Ã a 1j Plot 2: y j ¼b a ij Àk Ã a 1jf Ã a 1i and x j ¼k Ã a 2j . . .
Plot k a : y j ¼b a ij À P k a À1 r¼1k Ã a rjf Ã a ri and x j ¼k Ã a ka j The points on plot r (¼ 1. . .k a ) can be supplemented with a line which has slope given by the EBLUP of the variety score for that factor, that is,f Ã a ri . In the examination of the latent regression plots from their analysis, Cullis et al. (2014) state that ''Since all the estimated loadings for [the first] factor are positive this then means that large positive regression coefficients [scores] for this factor are desirable for DBH.'' This is an oblique reference to overall performance. In terms of stability, Cullis et al.

123
(2014) comment on the ''sensitivity'' of varieties to individual factors as depicted in the plots.
Latent regression plots are very informative, but are not an ideal tool for selection since they only provide an informal examination of overall performance and stability. Additionally, it can be a laborious task to examine and compare the plots for all varieties under consideration. Typically, breeders require one or two relevant measures to incorporate into a selection index. In the following we summarise information in the latent regression plots to obtain formal measures of overall performance and stability that can be used for this purpose.
In order to develop these measures we separate the EBLUPs of the common VE effects in Eq. (6) into the effects associated with the first factor and the remainder. Thus for variety i and environment j we havẽ If represented graphically for a single variety, this corresponds to the first latent regression plot. If all varieties are graphed together, this represents a series of m straight lines, with the slope for variety i being given byf Ã a 1i . The points along the regression lines are the fitted values for the first factor, which are given bŷ k Ã a 1jf Ã a 1i , and the deviations from the regression lines are given by Ã a ij . There are no explicit intercepts included in the model so that all lines pass through the origin, that is the point (0, 0). Thus if all the (rotated) estimated loadings for the first factor are positive, the lines for any pair of varieties do not intersect within the range of the data. This means that the fitted values associated with the first factor represent non-crossover interaction and this characteristic can be exploited to obtain measures of both overall performance and stability. The (relatively infrequent) scenario in which some of the rotated estimated loadings for the first factor are negative will be discussed in ''General applicability of FAST''.
The concepts are illustrated in Fig. 1 for two varieties, labelled V1 and V6, from the example. In terms of the fitted values for the first factor, that is, k Ã a 1jf Ã a 1i , it is clear that the rankings of the varieties do not change between environments because the regression lines do not intersect within the range of the loadings. This is due to the fact that the rotated loadings are all positive (see ''Appendix''). Furthermore, the lines diverge which indicates that the difference between the fitted values for the two varieties increases as the loadings increase. This is the classic representation of non-crossover VEI and is typically linked to changes in scale or the so-called ''discriminating ability'' of environments. It is therefore natural to use the fitted values for the first factor to form a measure of overall performance for each variety. The predicted scoresf Ã a 1i could be used for this purpose, but because they are standardised, it may be preferable to obtain a measure that is on the same scale as the trait being analysed. If the analogy with regression is continued, an obvious choice is the fitted value at the mean value of the ''regressor'' which is also the mean of the fitted values. If we let k 1 denote the mean of the rotated estimated loadings for the first factor, then the overall performance (OP) measure for variety i is computed as k 1f In the example, k 1 ¼ 10:0 and the OP of the two varieties is given by 34.9 cm for V1 and 18.8 cm for V6 (also see Fig. 1).
If the fitted values from the first factor regression represent non-crossover VEI for pairs of varieties, it is also natural to base measures of variety stability on the remaining factors. In this way, changes in variety performance due primarily to changes in scale are eliminated from the examination of stability. A single global measure of stability for each variety can be obtained as the root mean square deviation (RMSD) from the regression line associated with the first factor. This is given for variety i by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 p As with OP this measure is on the scale of the data. This measure of stability is easily visualised using the first latent regression plot. For example, Fig. 1 shows a much larger spread about the line for variety V6 compared with V1 and this is formally quantified by the RMSD of 5.6 cm for V1 and 12.7 cm for V6. In addition to the global measure of stability, responses of varieties to individual factors (excluding the first) may be of interest. These can be assessed visually using the latent regression plots 2. . .k a (see Cullis et al. 2014, for example). The essential information on each plot is the slope of the line which is given byf Ã a ri for variety i and factor r. These slopes could be used to rank varieties in terms of their responsiveness to the factors but they do not give any indication of the magnitude of the response with reference to the data. Given that the (rotated) estimated loadings for each of the factors 2. . .k a typically include both positive and negative values, (see ''Appendix'', for example), they reflect contrasts between environments. In this case, the magnitude of the response of a variety to the factor may be quantified as the average contrast in terms of the fitted values, namely the mean of the fitted values for positive loadings minus the mean of the fitted values for negative loadings. If we let k rþ and k rÀ denote the mean of the positive and negative estimated loadings for factor r, then the responsiveness of variety i to factor r is computed as This is illustrated graphically for varieties V1 and V6 for the second factor in Fig. 2. The means of the positive and negative estimated loadings are given by k rþ ¼ 4:4 and k rÀ ¼ À 5:3. The responsiveness to the second factor for each variety is represented by the vertical distance from the downward pointing open triangle to the upward pointing open triangle. This is given by -1.1-1.4 = -2.5 cm for V1 and 5.8-(-7) = 12.8 cm for V6. Thus, on average, the predicted common VE effects for variety V6 increase by 12.8 cm in response to the covariate implicit in the second factor, whereas they decrease by 2.5 cm for V1.

FAST for total VE effects
Given that factor analytic variance structures may be used for both the additive and non-additive VE effects, it is possible to extend the concepts developed in ''Factor analytic selection tools (FAST)'' for total (additive plus non-additive) VE effects. This is achieved by combining the models for the additive and non-additive VE effects in Eqs. (3) and (4) to obtain a special factor analytic form for the total VE effects. Due to the non-identity variance matrix for the additive effects it is necessary to re-scale the loadings and scores for these effects so they are comparable to the loadings and scores for the non-additive effects.
In order to apply selection tools, we first obtain EBLUPs of the common total VE effects as b g ij ¼b a ij þb e ij . The REML estimate of K g in Eq. (11) is obtained by replacing k a r and k e r with their (unrotated) REML estimates. Similarly the EBLUP of f g is obtained by replacing f a r and f e r with their EBLUPs. The matrix of estimated loadings is then rotated to a principal component solution as described in ''Rotation of REML estimates of loadings''. This providesK Ã g and alsof Ã g . The tools of ''Factor analytic selection tools (FAST)'' are then directly applicable. The implementation of this approach is the subject of future work.

Results and discussion
Application of FAST to motivating example Recall from ''Introduction'' that the aim of the analysis of the motivating example is to obtain EBVs to enable the breeding programme to select varieties to use as parents and to provide ratings for use by industry. In the RPBC breeding programme the emphasis on selection for the trait of DBH is on high EBVs across a wide range of environments. Thus OP and RMSD are two main drivers in the selection process. They are easily jointly assessed using an x À y plot for all the varieties under consideration for selection (see Fig. 3). Note that the selection process can be further enhanced by considering the accuracy of individual OP values. To this end, the points in Fig. 3 have been shaded lighter if the OP accuracy is less than 0.8. An examination of Fig. 3 allows easy and quick identification of varieties of interest. For example, variety V1 has the highest OP (34.9 cm) and an average level of stability (RMSD of 5.6 cm). We note that this is an existing parent and although it has the highest OP for DBH it is not widely deployed as it has an issue with spiral grain. Varieties V2 and V3 are two potential new parents (forward selections). Variety V2 has the second highest OP (33.7 cm) but is relatively unstable (RMSD of 11.7 cm), whereas V3 has a slightly lower OP (31.7 cm) and is very stable (RMSD of 1.4 cm). It is noteworthy that V2 and V3 share a common parent which is V1.
Once varieties of interest have been identified on Fig. 3, we then recommend examining the EBVs (EBLUPs of common VE effects) for these varieties in greater detail using first latent regression plots. These plots show the full set of EBVs for a variety and lend visual support to the OP and RMSD values. Varietal comparisons are aided by drawing pairs of varieties (often a test and control variety) on the same graph. Figure 4 contains the latent regression plots for the ten varieties labelled on Fig. 3, with variety V1 appearing on every panel for comparison. Varieties V2 and V3 were chosen to be graphed as they are potential new parents. The remaining varieties are existing parents and V5, V6, V8 and V10 were chosen as they are the four most widely deployed varieties; the other varieties were chosen because they show extremes in terms of stability and responsiveness. Note that the panels are ordered on OP, with the last panel corresponding to V10 which has the lowest OP of the ten selected varieties. It is clear from Fig. 4 that the OP of varieties V2 and V3 is very similar to that of V1 since the regression lines are almost co-incident (and thence the fitted values at the mean of the loadings are very similar). The regression lines for varieties V7 -V10 reflect fitted values that are consistently lower than for V1, with substantial differences at the mean of the loadings and hence much lower OPs compared with V1. In terms of stability, varieties V9 and V7 have a large scatter of points about their regression lines, and hence have large RMSD, whereas the points for V4 and V3 lie very close to the lines, so they have small RMSD (also see Fig. 3). As with Fig. 3, interpretation should take into account major variations in accuracy. Hence individual EBVs on Fig. 4 have been shaded lighter if their accuracy is less than for DBH for all 4608 varieties under consideration for selection as parents. Varieties V1-V10 labelled. Lighter coloured points correspond to varieties with an accuracy for OP of less than 0.8 123 0.8. We note then, for example, the accuracy of many of the EBVs for V2 and V3 are much lower than those for the existing parents due to the relatively limited testing of V2 and V3. In addition to examining overall variety stability using RMSD, it may be of interest to consider the individual responsiveness measures. Figures 5 and 6 plot OP against responsiveness for the second and third factors. These measures may be particularly useful if the associated factors represent meaningful environmental characteristics that can be exploited for specific adaptation. A method that is often proposed to assess this is to compute correlations between individual rotated vectors of estimated loadings and measured environmental covariates. In our experience this is rarely successful since the factors typically reflect complex combinations of environmental stresses. A potentially more fruitful approach involves a reversal of the focus in the sense of using variety rather than environment information. In this strategy, reference or probe varieties which are known to have differential performance in the presence of certain environmental conditions, are used to characterise the environments in the MET (Mathews et al. 2011). This has been successful for domesticated crops such as wheat, but is currently of limited use for Pinus radiata since it has a more complex genome and has only undergone one or two selection cycles. However, for pedagogical reasons, in the following we discuss the method in the context of the motivating example. Figure 5 shows that V9 has a large positive response to the second factor. Thus V9 had relatively better EBVs for environments with large positive loadings for this factor (eg. E83, E5, E61 as shown in the ''Appendix'') compared with environments with large negative loadings (eg. E65, E64, E17). In contrast, V7 exhibits the opposite behaviour. Additionally, neither V7 nor V9 show much response to the third factor (see Fig. 6). Thus the genetic architecture of V7 compared with V9 could be used to characterise the environments with extreme positive and negative loadings for the second factor. We note that these are very old varieties and have only average OP for DBH so their responsiveness is of interest mainly for characterising environments and possibly for the purpose of maintaining genetic diversity. In contrast we consider variety V2 which has a high OP and also a large negative response to the third factor (see Fig. 6) which means it had relatively better EBVs for environments with large negative loadings for this REML estimates of loadings for first factor EBVs (EBLUPs of common VE effect) Fig. 4 Motivating example: first latent regression plots for DBH for 10 varieties. Panels correspond to different varieties (V2-V10, coloured blue) and variety V1 is shown on every panel (coloured orange). Lighter coloured points have an accuracy of less than 0.8. The vertical dotted line corresponds to the mean value of the estimated loadings for the first factor factor (eg. E17, E83, E79 as shown in the ''Appendix'') compared with environments with large positive loadings (eg. E21, E76, E58). Thus the previously mentioned high RMSD for this variety may represent exploitable variation rather than undesirable instability. It will be particularly interesting to track this variety as more data are collected since it has the potential for both high OP and boosted performance under specific environmental conditions.
As previously mentioned, the interpretation of individual responsiveness measures in this example is problematic so that the focus in terms of stability is RMSD. In the RPBC breeding program it has been suggested that both OP and RMSD represent characteristics with economic importance. The fact that they are on the same scale, namely the scale of the trait being analysed, should aid in assigning economic weights. Given the weights, it will be straight-forward to combine OP and RMSD measures across key traits (including DBH) to form a simple, concise selection index.
In terms of information for industry, a certificate is issued for each seedlot purchased by forest owners and provides a rating for individual traits, including DBH. The rating is based on the genetic quality of the parents and their proportion in the seedlot. Currently the genetic quality of a parent is obtained in a piecemeal manner using the results of the MET analysis. In the first step, the estimated additive genetic correlations between trials are clustered in order to identify trials that lack correlation with others. EBVs for an individual parent are then averaged across all remaining trials to provide a single value which is then converted to a rating. With the advent of FAST as described in this paper, there is potential for a more coherent, objective and informative scheme for industry. The natural choice for a single measure of the genetic quality of a parent is OP. This avoids the need Responsiveness to second factor (cm) OP (cm) to ignore trials, which is a some-what subjective procedure, and also avoids the loss of information incurred by converting to a rating. In addition to OP, it may be important to consider stability, so that RMSD may also be reported to industry. We note that it is straight-forward to compute OP, RMSD and associated measures of accuracy on a seedlot basis given the specific mixture of parents. The challenge remains as to how this information may be disemminated and thence adopted by industry.

General applicability of FAST
The development and application of FAST thus far has assumed that all of the rotated estimated loadings for the first factor are positive. Here we discuss implications and departures from this assumption. First, it is instructive to make the comparison between OP and the more traditional concept of a variety main effect. The latter is typically obtained by fitting a linear mixed model that partitions VE effects into variety main effects and VEI. Thus the additive VE effects are given by where 1 p is the p-vector with all values equal to unity, u v is the m-vector of additive variety main effects (which has associated variance r 2 v A) and u ve is the mpvector of additive VEI effects (which has associated variance r 2 ve I p A). We note that the model in Eq. (12) can be re-written in the form of a factor analytic model, namely Experience has shown that the general FAk a model provides a far superior fit to most MET data-sets and hence our development of FAST in this paper. If all rotated estimated loadings for the first factor are positive, then the first factor represents a generalised version of the main effect part of the model in Eq. (12) in which the loadings are positive but no longer equal. OP for a variety may then be thought of as a ''generalised main effect'' which allows for heterogeneity of scale between environments.
Another key link with the model in Eq. (12) is in terms of the between environment additive genetic covariance structure. The model in Eq. (12) leads to a genetic variance matrix, G a , in which all diagonal elements are given by r 2 v þ r 2 ve and all off-diagonal elements by r 2 v . Thus the variety main effect variance can also be interpreted as the (common) genetic covariance between every pair of environments. In terms of the FAk a model, the use of the principal component rotation (see ''Rotation of REML estimates of loadings'') means that the first rotated factor accounts for the maximum amount of genetic covariance in the data. The estimate of this source of covariance for two environments j and h is given by the product of the corresponding two estimated loadings for the first factor, namelyk Ã a 1jk Ã a 1h . Thus if all the rotated estimated loadings for the first factor are positive, there is a dominant component of positive genetic covariance between all pairs of environments. Once again this is a generalisation of the main effect scenario in the sense that heterogeneity is accommodated.
In our experience in analysing MET data, the vast majority of results are characterised by rotated estimated loadings for the first factor that are either all positive or include a few small negative values. In terms of the analyses of the five key traits for RPBC, four traits, including DBH as presented in this paper, had all positive loadings in the first factor and the remaining trait had a single negative value. In the most recent annual analyses of grain yield data from four Australian pulse breeding programs, three out of eight analyses had all positive loadings in the first factor and the remainder had a few small negative values.
The presence of small negative estimated loadings in the first factor indicates that the first (and therefore dominant) factor contains cross-over VEI, but of such a small magnitude to be of no practical importance. We illustrate this graphically using the results of the eight pulse breeding analyses. Figure 7 contains schematic representations of first latent regression plots for these analyses. In each panel, two extreme varieties (A and B) are graphed to show the maximum cross-over VEI represented in the first factor. Only the fitted values (and not the EBLUPs of the common VE effects) are shown on the panels and they are plotted both as the regression lines and also the underlying points (to show the position of the environment loadings). Due to the presence of negative loadings in analyses 4-8, Fig. 7 shows that this causes the lines to intersect, which translates to cross-over VEI. However, the magnitude of the cross-over is very small, so that the superiority of B over A is probably not biologically or economically important and the dominant feature is the superiority of A over B, both in terms of frequency of environments and magnitude of differences. This is captured in OP which is higher for A than B (see fitted values at mean loadings in Fig. 7). Thus in cases where there are a few small negative loadings in the first rotated factor, OP and RMSD are still meaningful and we would proceed with FAST as described in this paper.
As the proportion and magnitude of the negative loadings in the first factor increases, there reaches a point where OP and RMSD are no longer meaningful. For example, if roughly half of the loadings in the first factor are negative, then the dominant feature in the VE effects is cross-over VEI so that an overall performance measure is inappropriate. The analogy with the main effect model in Eq. (12) is a reduction in the main effect variance estimate to a point where it is zero and thence the variety main effects effectively ''do not exist'', that is, the associated EBLUPs are all zero. This is a very rare scenario and requires a different approach for selection. Within the framework of FAST, one possibility is to compute a responsiveness to the first factor in addition to the remaining factors. Then graphs of all pairs of responsiveness measures, that is, of the form given in Figs. 5 and 6, could be used to aid with selection decisions. The magnitude of VEI dictates that selection must be 123 for specific rather than broad adaptation. It is key to note that FAST not only identifies that this is the case, but also offers a way forward for informed selection decisions, particularly if an interpretation can be ascribed to the loadings.