Application of Multiple Imputation for Missing Values in Three-Way Three-Mode Multi-Environment Trial Data

It is a common occurrence in plant breeding programs to observe missing values in three-way three-mode multi-environment trial (MET) data. We proposed modifications of models for estimating missing observations for these data arrays, and developed a novel approach in terms of hierarchical clustering. Multiple imputation (MI) was used in four ways, multiple agglomerative hierarchical clustering, normal distribution model, normal regression model, and predictive mean match. The later three models used both Bayesian analysis and non-Bayesian analysis, while the first approach used a clustering procedure with randomly selected attributes and assigned real values from the nearest neighbour to the one with missing observations. Different proportions of data entries in six complete datasets were randomly selected to be missing and the MI methods were compared based on the efficiency and accuracy of estimating those values. The results indicated that the models using Bayesian analysis had slightly higher accuracy of estimation performance than those using non-Bayesian analysis but they were more time-consuming. However, the novel approach of multiple agglomerative hierarchical clustering demonstrated the overall best performances.


Introduction
Multi-way data analysis has become common in many areas of research involving multivariate data. Three-way three-mode pattern analysis refers to the combined use of such clustering and ordination procedures. Its application to multivariate multi-environment trial (MET) data has provided a comprehensive summary of the patterns of variation and the interactions among the three modes, genotypes, environments and attributes, for plant breeders and other scientists interested in plant improvement [1,2]. However, many multivariate MET datasets are incomplete and the presence of missing values cause complications because most analytical methods developed for multivariate data assume complete data arrays [3,4]. This is the case for (iterative) clustering and ordination procedures where the inability to routinely apply them to incomplete datasets has been an obstacle to their wider usage (as a full data array is needed to provide starting values for any necessary iteration). Thus, it is important to obtain the best possible estimates of missing values to form a complete multi-way MET data array which can then be subjected to multi-way pattern analysis.
There are some statistical methods and mathematical algorithms specifically designed to handle incomplete two-way two-mode data matrices. In one of them, multiple imputation (MI) [5,6] is used to generate different imputed values for each missing value to form different complete datasets. Then the different complete two-way datasets were analysed in order to obtain estimates of the parameters of the corresponding models because these parameters were the main interest for some authors [7]. These different complete datasets were defined as the "estimated data arrays" as they were the complete data arrays containing the estimated missing values using MI approaches.
While we wanted to use multiple imputation to generate different imputed values for each missing cell (and eventually obtain one estimated data array for each incomplete multivariate MET dataset), the estimation of the (different) parameters in the various models used in the imputation process were not of concern to us. Thus, we focused on using different MI approaches to obtain "good" estimates of the missing values to form a complete "estimated data array" which could then be analysed by three-way three-mode pattern analysis, rather than for parameter estimation. The MI methods mentioned above (for two-way two-mode data matrices) were modified to take into account the three-way structure of multivariate MET data. We also introduced one novel MI approach which does not have an underlying model that can be written in a similar format to the others.
To demonstrate the use of MI for estimating missing values in multivariate MET data, two real complete MET datasets and four simulated complete MET datasets were considered. Missing values were generated by randomly deleting values in the full datasets. The methods were assessed by comparing the original complete data arrays with the "estimated data arrays", i.e., the complete data arrays containing estimated missing values. This enabled us to compare all of our methods for imputing missing values. Again, we stress that this was more important to us than the relative efficiency of the various estimators for the parameters in the models used in some of the imputation methods.
Some brief notation about the three-way three-mode data structure is described in the Materials and Methods. The basic algorithms for various MI approaches and corresponding modification in terms of multivariate MET datasets are also described. We then present the six multivariate MET datasets and the random generation of missing values, followed by the results of comparing the original complete data arrays with the complete data arrays containing estimated missing values. We end by discussing the implications of our findings.

Materials and Methods
Three-way three-mode MET data A MET data array generally consists of I genotypes, J environments and K attributes. It can be written as a collection of frontal slices X k (I×J matrix, k = 1,. . .,K), where rows are genotypes and columns are environments, and each cell x ijk is the value measured on the i th genotype in the j th environment for the k th attribute [8,9] (Fig 1).
Each frontal slice X k corresponds to I genotype responses across J environments for a particular attribute, such as yield, moisture or test weight. Each attribute has its own measurement units. There are two types of vectors for each frontal slice, a row vector x 0 ik (i.e. measurements on the i th genotype for the k th attribute for all J environments) and a column vector x jk (i.e. all I genotype responses measured in the j th environment for the k th attribute). We conducted column standardization in order to remove the environmental main effects, but retain the correlation among attributes (over genotypes) for each environment [1,10,11], i.e. we standardized products in development, marketed products or any restrictions whatsoever on the research work reported in this manuscript. Unrestricted funding for a scholarship was the requirement for independent consultancy work by Professor Kaye Basford for Monsanto Company. She has no patents, products in development, or marketed products associated with that company. There are no competing interests in relation to the research being reported in this manuscript. This competing interests statement does not alter or impact on the authors' adherence to PLOS ONE policies on sharing data and materials. each column vector x jk prior to the analysis. It can be defined as: where the environment main effect is removed as X ix ijk ¼ 0 and the correlation among attributes for each environment is retained as: Kroonenberg [9] discussed the various types of missing values in three-way three-mode data, and they are described here in terms of multivariate MET data: • Single observations missing, e.g., individual genotypes in particular environments for a specific attribute are missing. These are missing cells in the three-way three-mode array (Fig 2A).
• Column missing, e.g., a particular attribute is not measured on any genotype in a particular environment. This would correspond to a missing column (x jk ) in our three-way array ( Fig  2B) and is quite common. A missing row (x 0 ik ) where a particular attribute is not measured in any environment for a particular genotype is extremely rare in practice and will not be considered here.
Initially, we rearranged the standardized multivariate MET data array by writing it as a twoway wide matrixX IÂJ K , where the I rows are the I genotypes and JK columns are the J environments nested within each of the K attributes. It can be viewed as the following matrix including missing values: where NA indicates that the observation is not available (missing). Note that there are only missing cells and a missing column (corresponding to a particular attribute not being measured on any genotype in that particular environment), where we only consider one missing column as an example of missing columns.

Multiple imputation approaches
We are interested in two patterns of missing values in the two-way wide matrixX IÂJ K of the multivariate MET data, i.e. missing cells (e.g. a data value missing for a particular attribute for a particular genotype in a particular environment) and missing columns (e.g. data values missing for a particular attribute for all genotypes in a particular environment). We consider the different MI approaches in terms of how they take into account these missing patterns when: (1) the estimation task is the same for missing cells and columns; (2) the estimation task is different for missing cells and columns. Under (1), we study imputation approaches based on multiple agglomerative hierarchical clustering (MAHC) and the normal distribution model (NORM) [6]. Under (2), we study imputation approaches based on the normal regression model (NRM) [3] and predictive mean matching (PMM) [12][13][14]. The latter three common MI approaches (NORM, NRM and PMM) were modified in terms of the multivariate MET data structure. In each of these, Bayesian analysis under Gibbs sampling [15] and non-Bayesian analysis were used in implementing the estimation task.

Multiple Agglomerative Hierarchical Clustering (MAHC)
When investigating the genotype response pattern in multivariate MET data, it is most common to cluster the genotypes in terms of the measurements on the attributes in each of the environments. If there were only missing cells, then we could use an agglomerative hierarchical clustering procedure on the two-way wide standardized matrixX IÂJ K to provide an estimate of a missing value for a particular genotype (by using the attribute value within the same environment of the genotype (or the genotype group) with which the genotype with the missing value first merged in the agglomerative process). However, if an attribute is not measured on any genotype in a particular environment, we cannot do this. We were not able to devise any procedure for estimating missing values by clustering genotypes when a whole column of theX IÂJ K matrix was missing.
It was decided to rearrange theX IÂJ K matrix as anX J ÂIK matrix and cluster the environments. Then even when an attribute was not measured on any genotype within a particular environment, there would be no missing columns in the matrix. Hence we could use agglomerative hierarchical clustering of the environments to replace missing values for a particular attribute for all genotypes in a particular environment by the values of the same attribute for each genotype within the environment (or environment group) with which the environment with missing values first merged. This process of clustering environments to estimate missing values can also be used to estimate missing cells, i.e., when there was a missing value for a particular attribute for a particular genotype for a particular environment.
The process of estimating missing values using agglomerative hierarchical clustering of the environments, as described in the previous paragraph, only provided a single estimate for each missing value. We wanted a multiple imputation approach. In achieve this, we chose different subsets (or combinations) of attributes for each imputation of hierarchical agglomerative clustering of the environments, as that would result in (potentially) different hierarchies of the environments and subsequent estimates of missing values.
We initially tried to implement this MI procedure by choosing K h attributes (by a uniform sample from the K measured attributes) for imputation h, h = 1,. . .,H, (H being the total number of imputations) with the value of K h being any value from 1 to K. We conducted an agglomerative hierarchical clustering of the environments using the measurements of these K h attributes on the I genotypes in the J environments expressed in the form of a J×IK h two-way wide array. The environment response for a missing genotype-attribute value was estimated by the non-missing genotype-attribute value in the environment (or environment group) to which the environment with the missing value first joined. However, there was the possibility of not being able to estimate particular missing values because there were missing values on an attribute which was not chosen in that specific K h attributes.
Hence we modified our MI procedure to ensure that all attributes were chosen in each imputation h. Thus for each h, we decided to choose all attributes (as in agglomerative hierarchical clustering) as well as K h attributes with the value of K h being any of the values from 1 to K-1. We conducted one agglomerative hierarchical clustering of the environments using all of the K attributes (with data in the form of a J×IK two-way wide matrix) and another agglomerative hierarchical clustering of the environments using the measurements of these K h attributes (with data in the form of a J×IK h two-way wide array). In each case, the environment response for a missing genotype-attribute value was estimated by the non-missing genotype-attribute value in the environment (or environment group) to which the environment with the missing value first joined. If two estimates were obtained for a particular missing value, they were averaged. This process was repeated H times, and the final estimates of missing values were the averages of the corresponding estimated values from each of the H imputations.
This multiple agglomerative hierarchical clustering used a dissimilarity measure between environments (and environment groups) based on squared Euclidean distance and a grouping strategy which minimized incremental sum of squares [16]. It takes into account different combinations of attributes at each imputation, leading to average values (over the imputations) of the estimated missing values for each cell.

Normal Distribution Model (NORM)
Letx ik be the J×1 vector corresponding to the (standardized) responses for the i th genotype grown in all J environments for attribute k. The values are taken from the two-way wide standardized matrixX IÂJ K and would have been part of the i th row in that matrix. Then it is assumed that this vector,x ik is multivariate normally distributed with mean vector μ ik (a J×1 vector) and covariance matrix P ik ¼s 2 ik I (a J×J matrix), with I denoting an Identity matrix of size J×J. The corresponding mean vector μ ik and covariance matrix S ik are unknown. This could be written as: We obtain the estimated values of unknown parameters μ ik and S ik using Bayesian analysis and non-Bayesian analysis as follows.
To conduct the Bayesian analysis, we needed to construct the prior distributions of parameters μ ik and S ik . Here, we assumed that there was no strong prior information, so the Jeffrey's prior distribution [17] would be as follows: x ik jμ ik ; s 2 ik $ Nðμ ik ; s 2 ik IÞ; i ¼ 1; . . . ; I: As a constant value (here 1 is not a very realistic probability density function, it is an improper prior distribution for f ðμ ik ; s 2 ik Þ. However, when applying Bayes' rule, this constant value prior distribution leads to a proper posterior probability density function, introducing some information about μ ik and s 2 ik [17]. Then the posterior distribution was derived by Then the full conditionals of f ðs 2 ik jx ik ; μ ik Þ and f ðμ ik jx ik ; s 2 ik Þ were: is the sample covariance matrix for known μ ik , and The posterior distributions of parameters μ ik and s 2 ik were the normal distribution and inverse-gamma distribution, respectively. Using Gibbs sampling [15], the estimates of the parameters μ ik and s 2 ik were obtained. Gibbs sampling is Markov Chain Monte Carlo (MCMC) methodology. It was used to generate Z = (Z 1 ,Z 2 ,. . .,Z n ) from a target probability density function (pdf) f(z), given the conditional pdf f(z i | z 1 ,z 2 ,. . .z i−1 ,z i+1 ,. . .,z n ). During the process of Gibbs sampling, the Markov Chain was generated from a sequence of conditional distributions.
From the above derivation, the conditional posterior distributions of parameters μ ik and s 2 ik were the normal and inverse-gamma distribution, respectively. After the full distributions of parameters μ ik and s 2 ik were obtained, the estimates of μ ik and s 2 ik were obtained by generating a large number of extracted samples (here 5500) from those distributions and setting the estimate equal to their expectation values, i.e.μ ik ¼ Eðμ ik Þ; To conduct the non-Bayesian analysis, we also found the log-likelihood function to obtain maximum likelihood estimators (MLEs) for parameters μ ik and s 2 ik . The log-likelihood function forx ik was: Thus, the MLEs for parameters μ ik and s 2 ik were: Therefore, the missing valuesx m ik for the i th genotype measured on J environments for each attribute k could be drawn from the following equation: where z is random standard Normal vector. As a result, H imputed values for each missing values could be obtained by implementing the above process H times.

Normal Regression Model (NRM)
Let I m jk < I be the number of missing genotypes measurements in the j th environment for the k th attribute. Under the NRM, the particular standardised column vector (an I×1 vector)x jk containing I m jk missing values can be expressed as a partition into two vectors, one containing non-missing valuesx obs jk (an ðI À I m jk Þ Â 1 vector), and the other containing missing valuesx m jk (an I m jk Â 1 vector). The response (or measurement) for one attribute k measured on a particular genotype i in a particular environment j is not independent of the responses for the other attributes measured on that genotype in that environment. Hence, the column vectors (genotype responses) for attribute k in environment j are not independent of the column vectors (genotype responses) for the other attributes k 0 (k 0 = 1,. . .,K, k 0 6 ¼k) in that environment.
We assumed a standardized column vectorx jk satisfied the regression model is the usual regression coefficient, I is an I×I identity matrix, and X Ã is the I×N design matrix containing elements from the I genotype responses for the N columns in the two-way wide matrix.
In order to modify the models used with two-way data arrays to take account of the multivariate nature of the measurements on each genotype in each environment, we employed the design matrix X Ã which could be expressed in four ways (i.e. there were four options for the N columns). In these options, the "zero" elements in matrix X Ã (or X Ã m defined shortly) substituted for missing values. Thus, the missing values do not contribute to the normal regression model. However, there were a few "true" zero values in the data array because of the column standardization conducted prior to the analysis. These zero values do not contribute to the model either, as their original values corresponded to the average environment effect.
The four ways to express the design matrix X Ã were as follows: 1. The first option contained elements from the I genotype responses over the (JK-1) columns, i.e. I×(JK-1) matrix.
3. The third option was a combination of independent elements in I×K(J-1) columns and adjusted dependent elements in (K−1) × r kk 0 columns.
x I11 Á Á Áx Ij1 r 1kðjÞ Á Á Áx IJ 1 ::::x I1k Á Á Á Ã Á Á Á 0 ::::x I1K Á Á Áx IjK r kKðjÞ Á Á Áx IJ K where the column corresponding to the particular environment j for particular attribute k was discarded, but the other columns which were discarded in the second option were replaced by those columns multiplied by their corresponding correlation coefficient between attribute k and k 0 (k 0 = 1,. . .,K, k 0 6 ¼k) for the same environment j, i.e. there were (KJ-1) columns in the design matrix. The correlation coefficients r kk 0 (j) are described in the correlation matrix R (shown later).
4. The fourth option contained elements from the I genotype responses for the same environment over the (K-1) attributes, i.e. an I×(K-1) matrix.
The above design matrices could also be expressed as a partition into two matrices, one containing (I À I m jk ) rows in X Ã obs , and the other containing I m jk rows in X Ã m . Estimation of missing cells. When the parameters β jk and s 2 jk were determined (explained below), the missing cellsx m jk in one environment j for one attribute k could be drawn from Eq where z is random standard Normal vector, and X Ã m is the corresponding design matrix described above. For the third way of determining the design matrix X Cor j among attributes for any environment j, j = 1,. . .,J is: : We took the distinct upper off-diagonal elements and wrote them as a K(K-1)/2 row vector r 0 j ¼ ðr 12ðjÞ ; Á Á Á r k 0 ðk 0 þ1ÞðjÞ . . . ; r ðKÀ1ÞKðjÞ Þ, where k 0 = 1,. . .,K,. Then we combined all such vectors for each of the J environments to obtain a J×K(K-1)/2 full correlation matrix R, shown as follows: ; where NA indicates that the correlation coefficients between attribute k' and the other attributes for environment j 0 are not available (i.e. the genotype responses for environment j 0 for attribute k' are missing). Estimation of missing columns. When there is a missing column (x m j 0 k 0 , an I×1 vector), there are no observations for that particular attribute in that particular environment. It is impossible to estimate parameters β j 0 k 0 and s 2 j 0 k 0 from the non-missing values in this column (as there are no non-missing values). Therefore, we propose using the correlation coefficients to estimate the missing column.
However, as shown in the correlation matrix R, the correlation coefficients between attribute k 0 and the other attributes for environment j 0 were not available. Thus, the replacement of these correlation coefficients was computed by the average correlation ( r kk 0 ¼ , j6 ¼j 0 ) between attribute k 0 and the other attributes over the (J-1) environments. Then, the average correlation between attribute k 0 and the other attributes acted as the linear regression coefficient for this particular environment j 0 . Therefore, the particular missing column was estimated as follows:x where s is the average value of the σ jk from the other non-missing environments for each attribute (i.e. K(J-1)).
Alternatively, missing correlation coefficients in the matrix R could be considered as some of the elements in the response vector r 0 j 0 (j 0 = 1,. . .,J, j 0 6 ¼j) which is assumed to be linearly related to the respective correlation coefficients for the different environments j (j = 1,. . .,J, j6 ¼j 0 ) in R. That is equivalent to where β jj 0 could be obtained by the MLE, and then the (K-1) missing correlation coefficients (r k 0 1(j 0 ) ,. . .,r k 0 K(j 0 ) , k 0 = 1,. . .,K, k 0 6 ¼k) were estimated by linear regression. Therefore, the particular missing column was estimated by following equation: Eqs (3) and (4) were used to estimate the values in the other missing columns. We wanted to obtain H imputed values for the missing values, so the process of using Eqs (2), (3) and (4) was repeated H independent times.
Using the above process, the estimation of values in the missing columns (corresponding to a particular environment-attribute combination) is based on measurements of the observed attributes in the same environment and they are not independent of those other attributes measured on the genotypes in that environment. On the other hand, the estimation of values for single missing cells in particular columns is based on one of four ways of combining other observations, and the optimum combination is determined by the accuracy of estimation performance.
To obtain estimates of parameters β jk and s 2 jk in the above, both Bayesian analysis and non-Bayesian analysis were used.
To conduct Bayesian analysis, we assumed that the prior distribution of β jk was normally distributed with β jk * N(β 0 , V β ) and s 2 jk was inverse-gamma distributed with s 2 jk $ IGðn 0 =2; S 0 =2Þ. We followed Gelman [18] in assuming that V β = τ 2 I, and the hyperparameters β 0 , ν 0 , S 0 and τ 2 were fixed and known. As the sample variance of each observed column vectorx obs jk is 1, the shape parameter ν 0 /2 and scale parameter S 0 /2 have ν 0 and S 0 set to 4 and 2, respectively, Thus the mean of s 2 jk with inverse-gamma prior distribution is 1 ). In addition, the sample mean of each observed column vectorx obs jk is zero, hence, the mean of β jk with normal prior distribution is 0 (= β 0 ).
Thus, the full conditional distributions f ðs 2 jk jx jk ; β jk Þ and f ðβ jk jx jk ; s 2 jk Þ, are derived based on the Bayes' rule as follows: Here, both f ðs 2 jk jx jk ; β jk Þ and f ðβ jk jx jk ; s 2 jk Þ are proportional to the same product of distributions. Thus, Thus, the full distribution f ðβ jk ; s 2 jk Þ can be obtained using Gibbs sampling. From the above equations, when the prior distributions of parameters β jk and s 2 jk were set up as the normal and inverse-gamma distribution, respectively, their corresponding conditional posterior distributions were also normal and inverse-gamma distribution, respectively. After the full distributions of parameters β jk and s 2 jk were obtained, the estimates of β jk and s 2 jk were obtained by generating a large number of extracted samples (here 5500) from those distributions and setting the estimate equal to their mean values, i.e.: Then the estimation of single missing values could be obtained by Eq (2) using the above estimates of parameters β jk and s 2 jk , and the estimation of a missing column could be obtained by Eqs (3) and (4).
Alternatively, to conduct the non-Bayesian analysis, we needed to obtain the MLEs for parameters β jk and s 2 jk . The log-likelihood function forx jk is: Thus, the MLEs for parameters β jk and s 2 jk are: Therefore, the estimation of single missing values could be obtained by Eq (2) using MLEs for parameters β jk and s 2 jk , and the estimation of a missing column could be obtained using Eqs (3) and (4).

Predictive Mean Matching (PMM)
Rubin [12] proposed a statistical matching method for univariate nonresponse data while Little [13] developed and modified the method for multivariate nonresponse data, calling it predictive mean matching, where the respondent genotype vector (an I 0 ×1 vector, I 0 <I) satisfied a regression model. The parameters in the model (i.e. the regression coefficients β and residual variance σ 2 ) were determined by non-missing values. The predicted values of respondent genotype vectors, including non-missing and missing values, were obtained from the regression model. All predicted values for non-missing values were compared with the predicted values for the missing values by a distance function [13]. Then the C (C<<I 0 , e.g. 1, 2, and 3. . ., say 3) closest predicted values to a predicted missing observation implied that these particular C actual non-missing values could be used to estimate the missing value. One of these C values was chosen at random for each such missing value.
Estimation of missing cells. The regression model we used was again the normal regression model (NRM) asx jk $ N ðX Ã β jk ; s 2 jk IÞ; 8j; k. Here, the design matrix employed one of the four options described above, i.e. the one determined to have the most accurate estimation performance. We drew a bootstrap sample of I 0 observations (I 0 (I-I m )) from the non-missing (I-I m ) values in the vectorx obs jk (i.e. an (I-I m )×1 vector) for each imputation h, and put these values into a new vector _ x jk , which also contained the missing values, so it was of size (I 0 +I m ). It contained two components, _ x obs jk of size I 0 ×1, andx m jk of size I m ×1. The estimators of β jk and s 2 jk were obtained using the same procedure as we described in the NRM imputation above but with a design matrix _ X Ã which has (I 0 +I m ) rows.
Each of the predicted values within the_ x obs jk vector were compared with each of the predicted values within thex m jk vector using the following function [13]: Estimation of missing columns. For a missing column (environment j 0 in which attribute k 0 was not measured), the predict function for this particular column isx m j 0 k 0 ¼ , hence the predicted values of the observations in this missing column are related to the observations within the same environment j 0 measured for the other attributes k (k = 1,. . .,K, k6 ¼k 0 ). Then the predict function for these observations (x j 0 k ) iŝ where N k is a randomly selected uniform sample from the K attributes, and β Ã j 0 k is a N k (J-1)×1 vector. The computation of estimators of parameters β Ã j 0 k and s 2 j 0 k were the same as for the NRM imputation above using Gibbs sampling and MLEs.
As predicted values of observations within such a whole missing column werex m j 0 k 0 , the distance d 2 ij 0 ðkk 0 Þ was calculated between the elements ofx m j 0 k 0 andx j 0 k . Thus the estimate of each of the missing valuesx m j 0 k 0 in the missing column was randomly drawn from one of the C closest corresponding actual values inx obs j 0 k (k⊰N k , k6 ¼k 0 ), where these C estimates were determined from the C smallest values of d 2 ij 0 ðkk 0 Þ .

Comparison of methods
For each MI method discussed above, there were H imputed complete datasets, called "estimated data arrays", containing the observed values in the incomplete data array and the imputed estimates of the missing values in that array. The overall or final "estimated data array" was obtained by averaging the cells in the H imputed data arrays. It could also be determined by averaging the H estimates of the individual missing values and putting them into the incomplete data array to form a complete array. Comparisons between the original data arrays and the "estimated data arrays" were conducted using the normal root mean square error (NRMSE) [19] N RMSE ¼ wherex ijk is a standardized element in the original MET data array, andx ijk is a standardized element in the "estimated data array", as the standardized non-missing elements in the "estimated data array" are different from the standardized corresponding non-missing elements in the original MET data array. By considering each element of the data array in the NRMSE computation, we investigated both the influence of column standardisation prior to the analysis and the estimation performance. Also, the missing values could be estimated using the EM algorithm [9,20,21] in the three-mode ordination, referring to as the Tucker3 model [9,22]. This method of estimating the missing values is defined as single imputation [9]. We therefore needed to compare the techniques we are proposing for estimating missing values (estimates for missing cells and estimates for a missing column) with those generated by the Tucker3 model. The NRMSE values for each MI method and the EM algorithm were compared to determine which method was more accurate for estimating missing values.
In the above, we considered both missing cells and missing columns simultaneously. However, we could consider the two patterns of missing values separately, i.e. missing cells alone and a missing column alone. We did that by comparing the estimation performance in each case using the following criterion to test the efficiency of multiple imputation.
There were H different imputed datasets for missing values in the full "estimated data arrays". The corresponding imputed data values for the missing values,Q h , and variance of all missing values for each imputation h, U h , were obtained. For MI analysis, there are two types of variance [5,6,9,23]. One is called the within-imputation variance and defined by The other is referred to as the between-imputation variance and defined by h . Then the total variance associated with the overall estimate is T ¼ U þ ð1 þ 1 H ÞB [5]. Because we assessed these MI methods using complete data arrays from which we discarded values (those designated as "missing"), we know the original "true" or "actual" values of the missing values. Each element of the difference between the "true" values of missing values Q h and overall estimate Q divided by its overall standard deviation (ðQ h À QÞ= ffiffiffiffi T p ) has an approximate t distribution with degrees of freedom n H ¼ ðH À 1Þ½1 þ U ð1þH À1 ÞB 2 [5,24]. For such a t distribution, we could calculate the 95% confidence interval (CI) for the "actual" values, and then determine the percentage of these CIs which contained their corresponding "actual" values (of the total number of estimated values). This is referred to as Coverage CI [25]. The higher the value of Coverage CI, the more efficient the MI method is judged to be. It was useful to repeat the MI process, as a different selection of cells would be designated as missing each time. We arbitrarily chose 10 repetitions and subsequently calculated a mean Coverage CI and standard error of that mean for each MI method.
MI approaches were considered in two ways, i.e. the same estimation procedure for missing cells and a missing column, and different estimation procedures for missing cells and a missing column. Thus, we divided our investigation of missing values for each repetition into missing cells alone, missing column alone, and combined (equivalent to all missing values). For missing cells alone and all missing values, we considered 5%, 10%, 15%, 20% and 25% missing values. This enabled us to evaluate the estimation performance of the MI approaches for each of these situations. For a missing column alone, we either deleted each environment for each attribute or deleted 10 distinct random environments for each attribute.
Firstly, the estimation procedure for missing cells alone was conducted for each MI approach, to give 100 imputed values for each missing cell. These 100 values were randomly allocated into 5 sets of 20 values, and the average over each of the 20 imputed values for each percentage of missing cells for each repetition gave the final 5 (H) imputed values for each missing cell. Based on those 5 imputed values for each missing cell, we calculated the 95% confidence interval for the "true" value. Then we obtained the percentage of all estimated missing cells where the 95% CI contained the "actual" value of the missing cell (Coverage CI). As the MI process was repeated 10 times, we calculated the mean Coverage CI and its standard error. Note that the percentage of missing cells for each repetition was slightly less than the percentage of missing values being quoted as there was a designated missing column (in the incomplete two-way wide array) whose cells were not included here.
Secondly, a missing column (corresponding to an attribute not measured on any genotype in a particular environment) was included in each percentage of missing values for each repetition. As it was randomly specified, it could be different across repetitions. We decided to investigate the variability of missing columns in the estimation procedure by carefully considering the choice of columns across replicates. For the small datasets (i.e. when the number of environments was 10 or less) each column was sequentially selected as the missing column (environment) for each attribute for each imputation in turn. The 100 imputed values of each missing data value in the missing column were obtained for each replicate. Again, 5 average imputed values (obtained by dividing the 100 values randomly into 5 sets of 20) of each missing value in the missing column were used to calculate the 95% Coverage CI. The repetitions gave J Coverage CIs (for each attribute) from which a mean and standard error were obtained. For the larger datasets (i.e. when the number of environments was greater than 10), we did not sequentially select each column (environment) to be missing in the replicates. Instead, we randomly selected 10 environments without replacement to be the missing column (environment) for each attribute for each imputation in turn. Then the same calculation was applied to obtain the 95% Coverage CI for each replicate. This gave 10 Coverage CIs (for each attribute) from which a mean and standard error were obtained.
Finally, we considered all missing values (missing cells and a missing column simultaneously). We obtained 100 imputed values of each missing value (whether it corresponded to a missing cell or to a cell in a missing column) for each percentage of missing values for each repetition of the MI process. To compare with the EM algorithm for estimating missing values, we took the average over the 100 imputed values of each missing value to form the final "estimated data array" for each percentage of missing values for each of 10 repetitions. We already had these final "estimated data arrays" for each of 10 replicates for the EM algorithm. Then the estimation performances of the different MI procedures and the EM algorithm were compared using the normal root mean squared error (NRMSE) criterion for each of the 10 replicates. These were presented as box plots for each method for each percentage of missing values.

MET data arrays
We employed two real complete multivariate MET datasets and four simulated multivariate MET datasets. The first real multivariate MET dataset, described by Basford and Tukey [26], consisted of 58 soybean lines evaluated in 4 sites in Queensland Australia in each of 2 years (denoted as 8 environments) for 6 attributes, i.e. a 58×8×6 data array. The second was from the 2014 CIMMYT wheat breeding program. It contained 50 wheat lines grown in 31 environments with measurements on 4 attributes (i.e. a 50×31×4 data array). These were denoted as Datasets 1 and 2, respectively. The other four datasets were simulated to be of various sizes. They were based on other trials in the CIMMYT wheat breeding program, the first of size (60×10×6) and the second of size (80×15×6), denoted as Datasets 3 and 4, respectively, and maize trials from a commercial company, the third of size (100×20×5) and the fourth of size (120×60×4), denoted as Datasets 5 and 6, respectively.
For simulated MET datasets, we computed the variance components for each random effect (ŝ 2 G ;ŝ 2 E ;ŝ 2 GE ;ŝ 2 ε ) within each attribute from an analysis using the mixed linear model for threeway three-mode real MET datasets. The values of the variance components were estimated using restricted maximum likelihood [27,28] as implemented in the ASREML software [29]. As the estimated variance components for genotype, environment, genotype×environment, and residuals were obtained, the data values were randomly multivariate normally distributed with mean zero and corresponding estimated variance-covariance matrix from the real multivariate MET data (with only genotypes having non-zero off-diagonal terms). Then these values were used to form the full simulated datasets.

Results
The results for missing cells, a missing column, and overall missing values will be discussed in turn.
Firstly, Table 1 contains the mean and standard of Coverage CI (calculated over the 10 repetitions) for each MI estimation method for each percentage of randomly generated missing cells. Note that the percentage of missing cells is less than the percentage of missing values being quoted (because only missing cells were considered, not the missing column, i.e. the percentages in the table indicated the number of missing cells including those in the missing column).
As expected, the MAHC imputation was efficient, as the average coverage rates and their standard errors were good (means 74% or higher, standard errors less than 3.7%) for the real datasets (Datasets 1 and 2). For the simulated datasets (Datasets 3 to 6), the average coverage rates of MAHC imputation were similar to or slightly smaller than those for Datasets 1 to 2 (means over 71%), but the corresponding standard errors were somewhat higher (standard errors less than 4%) than those for the real datasets.  For the other imputation methods (NORM, NRM, and PMM) based on both analysis (Bayesian and non-Bayesian), the average coverage rates were lower and their standard errors were higher than those for MAHC at every percentage of missing cells for each dataset. Overall, both forms of NORM imputation had relatively higher coverage rates than those for the other two imputations, especially for larger datasets. The next best was PMM imputation.
On average, the results of Bayesian analysis of NORM imputations had slightly larger coverage rates than those for non-Bayesian analysis, but there was not much difference in coverage rates between Bayesian and non-Bayesian analysis for NRM and PMM imputations. The standard errors of the means for Bayesian and non-Bayesian analysis for the three imputation methods were quite similar at every percentage of missing cells. Also, the values of the standard errors of the coverage rates did not differ across the percentages of missing cells for most datasets.
Secondly, each environment was sequentially designated as the missing column for each attribute for the smaller datasets, while each of 10 random environments (chosen without replacement) was designated as the missing column for each attribute for the larger datasets. The mean and standard error of Coverage CI calculated over each environment missing in turn for each attribute for each estimation method for the smaller datasets and calculated over 10 environments (chosen at random without replacement) for each attribute for each estimation method for the larger datasets were presented in Table 2. The number of repetitions used to calculate average CI coverage for each dataset was 8 (for Dataset 1) and 10 (for Dataset 3), while the number of repetitions for Datasets 2, 4, 5 and 6 was 10. Therefore, Table 2 showed the average estimation performance for a missing column for each attribute for each dataset. The estimation procedures to estimate a missing column (environment) were the same as those to estimate missing cells for MAHC and both analysis of NORM imputation. For NRM and PMM imputations, we applied two relationships, "average" correlation coefficients and "linear" correlation coefficients, to estimate a missing column. The average Coverage CI of a missing column was expected to be similar across attributes for each dataset.
The MAHC imputation was very efficient as the average CI coverage rates were good (above 80%) for each dataset (Table 2). This was especially so for Dataset 1 where the average CI coverage rates were close to 85% for some attributes. In general, the average CI coverage rates for NORM imputation for both analysis were the next largest, followed by the other two imputations (NRM and PMM). The results of Bayesian analysis of NORM imputations had slightly larger average CI coverage values than those for non-Bayesian analysis. The standard errors of average CI coverage rates for MAHC imputation were smallest for a missing column for each attribute compared with other imputation methods. The standard errors for NORM imputation for Bayesian analysis were higher than those for non-Bayesian analysis for some attributes and lower for others. For NRM and NRM imputations, average correlation analysis had larger average CI coverage values and the lower standard errors than those from linear correlation analysis.
In general, the average coverage rates for both analysis of NORM imputation were relatively larger than those for NRM and PMM imputations. Again, the results of Bayesian analysis of NORM imputations had slightly larger average coverage rates than those from non-Bayesian analysis. The results of average correlation analysis of NRM and PMM imputations also had slightly larger average coverage rates and lower standard errors than those from linear correlation analysis.
Using the linear correlation analysis to calculate the missing column had slightly lower coverage rates than using the average correlation analysis (Table 2). Therefore, we subsequently used the average correlation analysis to compute the missing column for NRM and PMM imputation methods, in conjunction with Bayesian and non-Bayesian analysis to estimate missing cells for these two methods.
Consequently for missing values overall, we compared MAHC, NORM imputation with Bayesian and non-Bayesian analysis, NRM and PMM imputation with Bayesian and non-Bayesian analysis where the estimation of a missing column was via the average correlation analysis. For each MI approach and the EM algorithm for each percentage of missing values (including missing cells and a missing column), we obtained final estimated data values (by averaging over the 100 imputed estimates for each missing value) for each of 10 repetitions. The NRMSE values for assessing the accuracy of the MI approaches and the EM method for estimating missing values (by calculating the differences between the original "true" values and the estimated values) for each of the 10 replications for each percentage of missing values for each of Datasets 1 to 6 were displayed using a boxplot with the addition of the mean value . Low values of NRMSE correspond to better performance. For Datasets 1 and 2 (the real datasets), the values of NRMSE for estimating the performance of EM, NRM-NBA, and NRM-BA imputation were comparably larger than those for the other imputations at every percentage of missing values (Fig 3). Overall, MAHC, NORM-BA and NORM-NBA imputations had relatively lower values of NRMSE, so performed better.
For Datasets 3 and 4 (of size 60×10×6 and 80×15×6, respectively), the values of NRMSE for estimating the performance of EM, NRM-NBA, and NRM-BA imputations were relatively higher than those for the other methods (Fig 4). Generally, MAHC imputation had much better performance than the other imputation methods for these two datasets. The variability of NRMSE values for EM imputation for Dataset 3 was much larger than for the other imputation methods. Overall, the variability of NRMSE values for Dataset 4 was smaller than for Dataset 3.
For Datasets 5 and 6 (of size 100×20×5 and 120×60×4, respectively), the values of NRMSE for estimating the performance of EM, NRM-NBA, and NRM-BA imputations were relatively higher than those for the other imputations (Fig 5). MAHC imputation had much better performance than the other imputation methods. The variability of NRMSE for these two datasets was much smaller than for the other datasets.
The overall best performance for estimating randomly generated missing values (including both missing cells and a missing column) for these six datasets was MAHC imputation, followed by NORM imputation for both Bayesian and non-Bayesian analysis. The number of missing values in relation to the size of the data array had an impact on the accuracy of estimating performance, with larger values of NRSME corresponding to larger percentages of missing data. However, the variability of the NRMSE values was smaller when the size of the dataset increased.

Discussion
We investigated one multiple hierarchical clustering method (MAHC imputation) and three multiple imputation approaches (NORM, NRM and PMM) with and without Bayesian analysis for estimating missing values in three-way three-mode MET datasets.
To conduct Bayesian analysis using Gibbs sampling, we needed to assume some prior distributions for the parameters in the normal distribution model (NORM) and normal regression model (NRM). Also, we needed to set up some hyper-parameters in the prior distribution for NRM imputation. Thus, compared with NORM imputation, NRM imputation introduced more uncertainty and produced lower CI coverage rates and higher average values of NRMSE.
For the PMM imputation with and without Bayesian analysis, the estimated missing values were actual non-missing values from another cell plus a random standard normal term multiplied by an appropriate variance, where these actual non-missing values were those for which their predicted values were closest to the predicted missing values. The predicted values of non-missing values were obtained using the NRM model. The closest distance between predicted values of each cell based on the NRM model may differ from the closest distance between actual values of each cell. When that is the case, PMM imputation had lower coverage rates and higher average values of NRMSE than NORM imputation.
During the implementation of MAHC imputation, we employed different combinations of attributes for each imputation. Thus, the final estimated values using MAHC imputation were based on the average of the corresponding values from the same genotype in the environment (or environment group) deemed most like that particular environment when the "similarity" assessment used different combinations of attributes. It took advantage of multivariate measurements by combining the attributes in various ways in the estimation procedure. Another advantage was that there was no other uncertainty in MAHC imputation. The three other imputation approaches required the estimation of parameters in their underlying models. As a result, it is probably not surprising that MAHC imputation performed best.
For the three imputation models (NORM, NRM, and PMM), imputations with Bayesian analysis had slightly higher accuracy than those corresponding imputation models without Bayesian analysis, but they took more computer time to implement. For all comparisons with these imputations, MAHC imputation had the smallest implementation time and the highest accuracy.
For estimating a missing column, we made two different assumptions for the correlation among attributes for each environment. The correlation coefficients for the missing environment for a particular attribute were not available. Therefore, we decided that the missing correlation coefficients could be the average value of non-missing correlation coefficients or a linear combination of non-missing correlation coefficients. However, according to the results of our investigation, the first assumption (based on the average correlation coefficients) had better estimation performance than the second assumption (based on a linear combination of correlation coefficients). The linear combination of correlation coefficients was calculated using all other non-missing correlation coefficients among the attributes, while the average value of correlation coefficients was calculated from other non-missing correlation coefficients with this particular attribute. This is likely to be the reason that the average value introduced less uncertainty and produced higher CI coverage.
Multivariate multi-environment trial datasets are the focus of our research. As many of these are incomplete, it is important to estimate the missing values to form a complete array that can then be analysed by three-way pattern analysis methodology which has proven valuable for that situation. However, the above imputation methods could apply to other types of datasets where various measurements are made on the same entities under different conditions.
Supporting Information S1 Datasets. Datasets used in this study. This compressed file contains six directories for each of six datasets. A description of each text file is included in the corresponding directory. (ZIP) S1 File. R source code for multiple imputation approaches. (ZIP)