BACKGROUND
due to the increase in available data information  <cit> , integrating large amounts of heterogeneous data is currently one of the major challenges in systems biology.
biological data integration provides scientists with a deeper insight into complex biological processes.
however, when dealing with multiple data tables, the presence of missing values is a common situation for a variety of reasons.
in omics data integration studies, it is common for some individuals to not be present in all data tables, resulting in a specific missing data pattern for multiple tables, as shown in fig.
<dig>  for instance, in clinical studies, this can occur when a patient forgets to fill out a form.
it also can be attributable to the study design if individual data are expensive or difficult to measure.
fig.
<dig> pattern of missingness.
pattern of missing values specific to multiple data tables with missing rows for some tables of variables.
the colored parts represent strata of the observed data, whereas the empty parts indicate missing data
missing row values for a table of variables are challenging to handle because most statistical methods cannot be directly applied to incomplete datasets.
in the multiple multivariate framework, several approaches have already been proposed to deal with missing row values  <cit> .
the only methods widely available for analyzing incomplete data focus on removal of the missing rows, either by ignoring subjects with incomplete information or by replacing the missing items with plausible values .
in multivariate statistical analysis, case deletion procedures can be very inefficient, discarding an unacceptably high proportion of subjects because even if the per-table rates of missing rows are low, only a few subjects may have complete data for all tables.
in addition, case-deletion procedures may bias the results if the remaining subjects providing the complete data are unrepresentative of the entire sample  <cit> .
on the other hand, simple mean substitution can seriously distort the marginal and joint distribution of the variables  <cit>  and be an issue because many statistical methods rely on estimation of the variance-covariance matrix.
recently, two approaches have been proposed to deal with missing row values in multiple multivariate analysis.
the first method, introduced by van de velden and bijmolt   <cit> , was developed in the context of generalized canonical correlation analysis.
its application in the omics framework is often limited by the size, noise and multicollinearity of the data  <cit> .
the second method, described in husson and josse   <cit> , was developed in the context of multiple factor analysis .
this method, designated regularized iterative mfa , was derived from a method available in principal component analysis  and consists of alternating the estimation of axes and components, and the estimation of missing values  <cit> .
here we consider an alternative method, involving a multiple imputation approach adapted to the mfa framework, and called mi-mfa.
multiple imputation  was proposed by rubin   <cit>  in order to estimate both the parameters of interest and their variability in a data missingness framework.
it relies on the principle that a single value cannot reflect the uncertainty of the estimation of a missing value.
first, mi is used to generate plausible synthetic data values, called imputations, for missing values in the data.
this step results in a number  of imputed datasets in which the missing data are replaced by random draws of plausible values according to a specific statistical model.
the second step consists of analyzing each imputed dataset using a statistical method that estimates the parameters of interest.
this step results in m analyses  which differ only because the imputations differ.
finally, mi combines all the results together to obtain a single consensus estimate, thereby combining variation within and across the m imputed datasets.
under fairly liberal conditions, this last step results in statistically valid estimates that properly reflect sampling variability.
the major challenge in mi involves generating possible values for each missing observation.
statistically advanced imputation procedures can therefore be used for this.
two general approaches are often used for imputing multivariate data: joint modeling   <cit>  and fully conditional specification , also known as multivariate imputation by chained equations   <cit> .
jm involves specifying a multivariate distribution for the missing data and drawing imputations from their conditional distributions by markov chain monte carlo  techniques.
fcs specifies the multivariate imputation model on a variable-by-variable basis using a set of conditional models, one for each incomplete variable.
the key issue in jm is appropriate specification of the multivariate distribution.
a multivariate normal model has often been used as it is computationally tractable .
this model has even been used when some of the variables are not gaussian.
however, the main weakness of jm is that it can only be applied when the imputation involves a small number of variables.
this is not very common in omics datasets that are often composed of tens of thousands of variables, or more.
fcs allows greater flexibility than jm in creating multivariate models.
indeed, fcs can use specialized imputation models by separately defining the conditional densities for each variable, even if this can require a considerable amount of work.
when the number of variables is large, it is often impractical and too computer-intensive to test and develop the best models for each variable.
as an alternative to the jm and fcs approaches, we propose using the hot-deck imputation approach  <cit> .
this approach is a nonparametric imputation method that resolves the most important limitation of the jm and fcs approaches as it can be applied to data tables containing more than just a few variables.
when using the mi method, special attention must be given to the process that gave rise to the missing data, referred to as the missing data mechanism.
most methods for generating multiple imputations, fully-, semi- and non-parametric methods, assume that the mechanism responsible for missing data is ignorable  <cit> .
briefly, if the missing data mechanism is ignorable, then the analysis can focus on the observed values rather that also having to model the process that resulted in certain values being observed and certain values being missing.
if the assumption of an ignorable missing-data mechanism is valid, then statistical methods that rely on that assumption can be expected to produce results with minimal bias.
one way that the missing-data mechanism can be viewed as ignorable is if the missing data are missing completely at random .
for data to be mcar, there must not be any systematic differences between the cases that have missing items and the cases that are fully observed.
in microarray experiments, technical failure, low signal-to-noise ratio and measurement errors can for instance be considered as sources of mcar patterns.
the missing-data mechanism can also be viewed as ignorable under the less restrictive missing at random  scenario, which allows missingness to depend on observed variables but not on unobserved variables.
late-stage cancer patients, as compared to early-stage cancer patients, unfortunately have more chance of dropping out of follow-up studies, which may result in a mar pattern in a clinical data table.
when the ignorability assumption does not hold, the imputation needs to be drawn from the posterior distribution of the missing data given the complete data and the missingness mechanism.
non-ignorable missing data occurs frequently in mass-spectrometry-based experiments.
measures too close to the limit of detection of the instrument are censored, resulting in a higher rate of missing values.
the probability of being missing is, in this particular case, directly dependent on the intensity value.
in this paper, we decided to focus on models with ignorable missing-data mechanisms.
methods
mathematical basis of multiple factor analysis
mfa  <cit>  is devoted to the simultaneous exploration of multiple data tables where the same individuals are described by several tables of variables.
in mfa, the number and the type of variables  may vary from one table to another, but within each table the nature of the variables is the same.
here we focus on quantitative variables.
the aims of mfa are similar to those of pca, namely to study the similarities between individuals from a multidimensional point of view, to analyze the relationships between variables and characterize individuals based on these relationships.
however, beyond these conventional uses, mfa can also be used to study the links between tables of variables and to compare the information contributed by each table.
mfa analyzes a set of j data tables k <dig> …,kj, where each kj corresponds to a table of quantitative variables measured on the same i individuals .
the core of mfa is a pca in which weights are assigned to variables.
more formally, the matrix of variance-covariance associated with each data table kj is decomposed by pca and its largest eigenvalue λ1j is derived.
then, each variable belonging to kj is weighted by 1/λ1j.
finally, a global pca is performed on the merged and weighted data table k= to obtain the configuration f .
the main reason for the weighting step is to remove from each table all information related to its own dimensionality or variance.
therefore, no single table can dominate the first dimension of the global analysis.
mfa provides the same graphical representations as pca  but also, due to the table structure, specific representations such as the table representation and the superimposed representation are available  <cit> .
the multiple imputation multiple factor analysis approach
to deal with multiple tables with missing rows, we propose the mi-mfa approach, a multiple imputation  adapted to the framework of mfa.
the aim of our method is not to get the best possible estimations of the missing values, but to replace them with plausible values in order to provide estimates of the mfa configurations.
according to mi methodology, the mi-mfa approach is carried out by performing the following three steps:
<dig>  imputation: generate m different imputed datasets k,…,k,…,k of k.
<dig>  mfa analysis: perform an mfa on each k imputed dataset leading to m different configurations f <dig> …,fm,…,fm.
<dig>  combination: find a consensus configuration between all f <dig> …,fm configurations.
these steps are outlined in fig.
<dig> and described in detail in the following sections.
fig.
<dig> overview of our mi-mfa approach to handling missing rows in multi-omics data integration.
the top part of the graphic indicates that analysis starts with observed, incomplete data tables k. in a second step, multiple imputation is performed using the hot-deck imputation approach: m imputed versions k
,…,k
of k are obtained by replacing the missing values by plausible data values.
these plausible values are drawn from donor pools.
the imputed sets are identical for the non-missing data entries, but differ in the imputed values.
the third step is to estimate the configuration matrix f
m for each imputed dataset k
using mfa.
the estimated configurations differ from each other because their input data differ.
the last step is to combine the m estimated configurations f
<dig> …,f
m into a compromise configuration f
c using the statis method
generating imputed data: multiple hot-deck imputation
hot-deck imputation involves replacing missing values of one or more variables with available values from a similar unit  <cit> .
the observation from which these available values are taken for imputation is called the donor and the observation with the missing value, which receives the donor’s value, is the recipient.
the donor can be randomly selected from a set of potential donors, called the donor pool.
selection of a suitable donor pool is not an easy task and is beyond the scope of this article  <cit> .
the general principle is to choose donor units that are as close as possible to the recipient with respect to some affinity score.
affinity is defined in terms of the degree to which each potential donor matches the recipient’s values across all variables other than the one being imputed.
intuitively, in the framework of stratified multiple omics tables, the donor pool can be formed of available individuals belonging to the same stratum  and the same omics table as the recipient.
multiple hot decking differs from other forms of hot decking by using several donors for a single recipient  <cit> .
multiple hot-deck imputation proceeds as follows.
let k= be the merged data table containing missing rows with strata s= <dig> …,s, then carry out the following steps :
step  <dig>  create donor pools by taking donors belonging to the same stratum s and the same table kj as the recipient.
recipients within the same stratum have the same donor pool.
suppose always that there is a large enough number of donors for recipients in each stratum.
step  <dig>  for each recipient in kj, impute the missing individual by drawing randomly with replacement a donor from the corresponding donor pool.
repeat this procedure until all missing individuals in the j tables have been imputed.
step  <dig>  repeat step  <dig> until m different imputed datasets k,…,k of k are obtained.
by conducting the imputations in this way, it is reasonable to assume that the within-unit between-variables multivariate relationships are preserved.
the combination procedure: the statis method
the question that arises after using mi in an mfa framework is how should all the configurations resulting from the analyses be combined to obtain a single unique estimate of the consensus configuration?
while averaging is an appropriate combination procedure in many other statistical techniques, it is not recommended for mfa due to possible reflection, dilation or rotation of the different configurations with respect to each other  <cit> .
here we consider an alternative approach by implementing the statis method which provides a compromise configuration balancing all configurations.
the statis method  <cit>   is a generalization of pca used to simultaneously study several tables of variables collected on the same individuals.
the goal of this method is to analyze the structure of the individual tables  and to derive from this structure an optimal set of weights for computing a common configuration of the observations.
the solution obtained, called the compromise, is the configuration agrees the most with all other configurations.
an overview of the statis method is presented in additional file 1: figure s <dig>
statis analyzes a set of n tables x <dig> …,xn, where each xn is a table of quantitative variables measured on the same individuals.
the first stage of statis consists of calculating a matrix of cross-products between individuals for each table wn=xnxnt  reflecting the similarities between individuals within this table.
the use of matrices wn instead of xn simplifies the computation because it obviates the determination of rotations when matching the xn.
the basic idea in statis is then to find a compromise space wc=∑n=1nαnwn that globally balances these cross-product matrices by choosing a suitable optimal set of weights α <dig> …,αn.
these weights are obtained from the pca of the matrix r whose generic term rjk gives the cosine between tables  defined as:
rjk=tracetrace·trace,  where the trace is the sum of the main diagonal elements of a square matrix.
the first eigenvector obtained from the pca of r represents the “agreement between tables”.
its elements are normalized in such a way that their sum is equal to  <dig> and used as weights αn in order to define wc.
tables with larger values of αn are more similar to the other tables and therefore will have a larger weight, while the weight of the “outlier tables” will be closer to zero with respect to the other weights.
the principal components from the pca of wc then gives the coordinates of the individuals in the compromise space, called the compromise configuration.
implementation of mi-mfa
the mi-mfa algorithm can be summarized as follows :
step  <dig>  start with an observed, incomplete dataset k. define the number of imputations m and the dimensionality d of the compromise configuration.
step  <dig>  perform multiple hot-deck imputation.
for m= <dig> …,m:
obtain an imputed version k of k, such that, k≠k for m≠m′.
the imputed datasets are identical for the non-missing data entries, but differ for the imputed values.
the imputed version of the data is obtained by using the hot-deck imputation approach.
perform an mfa using d components on the imputed dataset k to obtain the configuration fm.
step  <dig>  perform a statis on the set of configurations f <dig> …,fm to obtain fc, the compromise configuration.
note that the number of dimensions d used in the algorithm has to be chosen a priori.
however, the number of dimensions does not affect the estimation of the imputed values and the estimation of the compromise configuration.
moreover, for given k,…,k imputed datasets, solutions provided by the algorithms are nested .
since the core of mi-mfa is a weighted pca, the strategies suggested to choose the number of components in pca can be adapted to mi-mfa, but work needs to be done to validate the quality of these extensions.
how many imputations?
when using mi, one of the uncertainties concerns the number m of imputed datasets needed to obtain satisfactory results.
the number of imputed datasets in mi depends to a large extent on the proportion of missing data.
the greater the missingness, the larger the number of imputations needed to obtain stable results.
however, in multiple hot-deck imputation, the number of imputed datasets is limited by the size of the donor pools.
in any case, the total number of possible imputations mtotal can be calculated before applying the imputation approach .
if mtotal is small , then m=mtotal can be used in mi-mfa.
the appropriate number of imputations can be informally determined by carrying out mi-mfa on n replicate sets of ml imputations for l= <dig> , <dig> …, with m0<m1<m2<⋯<mtotal, until the estimate compromise configurations are stabilized.
more precisely, this approach can be carried out by applying the following steps:
step  <dig>  start with an observed, incomplete dataset k. define the number of imputations ml with m0<m1<m2<⋯<mtotal and the number n of replicate sets of ml imputations.
step  <dig>  create collections inml, n= <dig> …,n, each one containing ml different imputed datasets of k, such that, inml≠in′ml, for n≠n′ and inml−1⊂inml for m0<m1<m2<⋯<mtotal.
step  <dig>  for n= <dig> …,n, perform an mi-mfa using inm <dig>  to obtain n different compromise configurations fc1m <dig> …,fcnm <dig>
step  <dig>  let l= <dig>
for n= <dig> …,n,
perform an mi-mfa using the collection inml, to obtain a compromise configuration fcnml;
calculate rnl=r, a measure of the distance or correlation between configurations .
step  <dig>  calculate r¯l=1n∑n=1nrnl  the standard error of r¯l).
step  <dig>  repeat steps  <dig> to  <dig> for l= <dig> ,… until the differences between two subsequent r¯l ) become smaller than a certain convergence criterion.
uncertainty of mi-mfa solutions
in an mi-mfa framework, after estimating the configurations from the imputed datasets, a new source of variability due to missing values can be taken into account.
here we describe two approaches to visualize the uncertainty of the estimated mfa configurations attributable to missing row values.
first, an individual plot for all estimated mfa configurations is constructed.
the individual plot is obtained by projecting each estimated mfa configuration onto the compromise configuration .
each individual is represented by m points, each corresponding to one of the m mfa configurations.
confidence ellipses and convex hulls can then be constructed for the m configurations for each individual.
the computed convex hull results in a polygon containing all m solutions.
all individuals have confidence areas, even those without missing values.
indeed, even if only the estimation of missing values is the only change, this will have a possible impact on all mfa parameters.
therefore, the area of such an ellipse  provides an insight into the uncertainty of the estimated configuration.
the larger the area of an ellipse , the more uncertain the exact location of the individual.
thus, when the area of an ellipse is large, the scientist should remain really careful regarding its interpretation.
performance of the method
we conducted two case studies to assess the performance of our method.
instead of using theoretical distributions to generate simulated data, our studies were based on two real datasets, denoted as the original datasets.
subsequently, specific patterns of missingness were created in these datasets as illustrated in fig.
<dig>  resulting in what we called the incomplete datasets.
this approach was used in order to more closely mirror situations that may occur in the omics context.
next, missing row values were estimated and the resulting complete datasets were referred to the imputed datasets.
we then compared our mi-mfa method to the ri-mfa approach  <cit>  and the mean variable imputation mfa  method, in which the missing values are simply replaced by the mean of each variable after which an mfa is performed on the imputed dataset.
this latter approach was considered as the common base for comparing the mi-mfa and ri-mfa methods.
for each configuration obtained using mi-, ri- and mvi-mfa, the similarity between the configuration solution and the true configuration  was assessed from the rv coefficient  <cit> .
the rv coefficient, which ranges from zero to one, can be interpreted as a correlation coefficient between two matrices, which allows the relative positions of objects to be compared from one configuration to another.
we also provide comprehensive graphical comparisons of the true vs. the mi-, ri- or mvi- mfa configurations.
the individuals from both configurations are drawn in a same plot and connected by an arrow, the length of which indicates the divergence between the two configurations.
implementation of the analyses
all analyses were performed using the r computing environment  <cit> .
mfa was performed using the mfa function of the factominer r package  <cit> .
the statis function of the ade <dig> r package  <cit>  was used to determine the compromise configuration.
the ri-mfa method is implemented in the imputemfa function available in the missmda r package  <cit> .
note that the number of components ncp used to predict the missing entries in the imputemfa function has to be chosen a priori.
this choice is crucial and difficult  <cit> .
as the true configuration was known in our case, the number of components ncp was chosen to minimize the rv coefficient between the true and the imputemfa configurations.
the appropriateness of the results from mi-, ri- and mvi- mfa was then determined by comparing the configurations resulting from these three strategies with the true mfa configuration.
due to a possible lack of alignment  between two configurations , it was necessary to align them before being compared.
ordinary procrustes analysis  <cit>  was used to align these configurations prior to their comparison.
datasets
liver toxicity
the datasets originated from a liver toxicity study  <cit>  in which  <dig> male rats of the inbred fisher f344/n strain were exposed to toxic doses of acetaminophen  in a controlled experiment.
necropsies were performed  <dig>   <dig>   <dig> and  <dig> h after exposure and mrna was extracted.
the data consisted of the expression of  <dig>  genes and  <dig> clinical variables considered to be markers of liver injury.
the  <dig> subjects  were cross-classified in eight strata  according to two factors:
exposure time:  <dig>   <dig>   <dig> and  <dig> h;
toxic doses of acetaminophen: high  or low .
eight subjects per stratum were included.
these datasets were downloaded from the mixomics r package  <cit> .
nci- <dig> data
the nci- <dig> dataset contained transcriptomic  <cit>  and proteomic  <cit>  tables for a collection of  <dig> cell lines from the national cancer institute .
the nci- <dig> panel included cell lines derived from various cancer types: colon , renal , ovarian , breast , prostate , lung  and central nervous system , as well as leukemia  and melanoma .
the gene expression profiles used here were generated using an agilent platform  <cit>  and downloaded from cellminer  <cit> .
data were log2-transformed.
to facilitate data interpretation and computation, the transcriptomic data were filtered to exclude probes that did not map to an official hugo gene symbol and to retain only the probe with the highest average value when multiple probes mapped to the same gene, as previously described in  <cit> .
gene invariants across all  <dig> cell lines, corresponding to genes without any effect between cancer types, were also removed.
filtering produced a dataset of  <dig>  genes.
the nci- <dig> proteome table was also downloaded from cellminer  <cit> .
proteomic data were obtained using high-density reverse-phase lysate microarrays  <cit> .
data were log2-transformed and protein abundance levels were available for  <dig> proteins  <cit> .
RESULTS
liver toxicity data analysis
a specific pattern of missing values was created as illustrated in fig.
<dig>  to obtain an incomplete dataset, three individuals per stratum were randomly removed from the transcriptomic table.
for this specific pattern, there were 3×8= <dig> missing individuals.
mi-mfa was then performed on the incomplete dataset using m= <dig> imputed datasets.
ri-mfa, and mvi-mfa were also performed.
figure  <dig> shows the divergence of the mi-, ri- and mvi-mfa configurations from the true configuration.
as can be seen, the configuration obtained with mi-mfa was very close to the true configuration .
this result was confirmed by the high rv coefficient .
the configurations obtained with ri-mfa and mvi-mfa were more distorted and less close to the true configuration with rv coefficients of  <dig>  and  <dig>  respectively .
fig.
<dig> comparison of mfa configurations for the liver toxicity data analysis.
representation of the individuals on the two-dimensional configuration obtained by performing mfa on the complete dataset  and, mi-mfa , ri-mfa  and mvi-mfa  versus the true configuration.
dots and arrows represent the projected coordinates from the true configuration and the imputation mfa method, respectively.
the length of the line joining dots and arrows is proportional to the divergence between the projected coordinates.
the color of each individual  reflects the treatment .
imputed individuals from the transcriptomic dataset are represented by empty circles
the number m of imputed datasets in mi-mfa for the above incomplete liver toxicity data was determined as described in the how many imputations section.
collections of size n= <dig> were generated for each of the following numbers of imputations: ml= <dig> l, for l= <dig> …, <dig>  the stability of the estimated mi-mfa configurations was then determined by calculating the rv coefficient between the configurations obtained using ml and ml+ <dig> imputations .
as the true configuration was known, we also described the stability of the estimated mi-mfa configurations by calculating the rv coefficients between the true configuration and those obtained using ml imputations .
although the missing information is substantial, fig.
<dig> shows that only a slight increase in precision was obtained by using more than  <dig> imputations.
fig.
<dig> stability of the estimated mi-mfa configurations using different numbers of imputations.
mi-mfa configurations were obtained for the following numbers of imputations: m
l= <dig> l, for l= <dig> …, <dig>  left: rv coefficient between configurations obtained by mi-mfa with m
l and m
l+ <dig> imputations.
right: rv coefficient between the configuration obtained by mfa on the complete dataset and mi-mfa on the incomplete dataset with m
l imputations.
the values shown are the mean rv coefficients for the n= <dig> two-dimensional configurations as a function of the number of imputations.
error bars represent the standard deviation of the rv coefficients
mi-mfa was also applied to different scenarios of missingness.
first, mfa was performed on the original dataset to obtain the true configuration.
secondly, individuals were randomly removed from each stratum in the original transcriptomic table.
for this, three scenarios were considered for the number of missing rows , in which there were respectively one, two and three missing rows per stratum.
the total number of missing rows per scenario was therefore  <dig>   <dig> and  <dig> respectively.
thirdly, for each missingness scenario,  <dig> incomplete datasets were randomly chosen for analysis.
as previously, mi-mfa was performed on each incomplete dataset using m= <dig> imputed datasets.
ri-mfa, as well as mvi-mfa, were also computed.
the rv coefficients between the true configuration and the configurations obtained using each method  were then calculated.
figure  <dig> shows the mean of the rv coefficients for the  <dig> two-dimensional configurations as a function of the missingness scenario for each method.
note that the average results using mi-mfa were always better than with ri-mfa or mvi-mfa, whatever the scenario.
the rv coefficients between the true configuration and the mi-mfa configuration were close to one and remained satisfactory even when the number of missing row values was high, and the results obtained with the ri-mfa and mvi-mfa decreased significantly.
fig.
<dig> performance on liver toxicity data according to the number of missing rows.
rv coefficients between the configuration obtained by mfa on the complete dataset and either mi-mfa , ri-mfa  or mvi-mfa  on the incomplete dataset.
the values shown are the mean rv coefficients for the  <dig> two-dimensional configurations for each missingness scenario: low, medium and high .
error bars represent the standard deviation of the rv coefficients
the performance of the mi-mfa procedure was then further investigated for more complex scenarios of missingness.
more precisely, missing row values were inserted into each stratum  of the original dataset  according to the scenarios illustrated in table  <dig>
number of missing rows inserted in each stratum of the original dataset, including both transcriptomic and clinical data, for incomplete data creation.
the # cases indicate the number of possibilities of incomplete cases per stratum
twenty incomplete datasets were then selected at random from each stratum  and each scenario.
all analyses  and the calculation of the rv coefficient were performed in the same way as previously described.
figure  <dig>  shows the means of rv coefficients for the two-dimensional configurations as a function of the scenarios for each method.
fig.
<dig> performance on liver toxicity data with missing individuals in the high-dose treatment.
average rv coefficient between the configuration obtained by mfa on the complete dataset and either mi-mfa , ri-mfa  or mvi-mfa  on the incomplete dataset.
results are given for all of the first three two-dimensional possible configurations as a function of the scenarios presented in table  <dig>  the discrete rv values are joined by lines for ease of understanding.
the performance of the different methods on the liver toxicity data with missing individuals in the low-dose treatment are presented in additional file 1: figure s3
for almost all the scenarios, the average results obtained with mi-mfa were better than with the other methods.
as the number of missing row values increased , the results obtained with the ri-mfa and mvi-mfa algorithms worsened rapidly, especially when the dimension increased, whereas the results with the mi-mfa approach were still satisfactory.
nci- <dig> data analysis
to confirm the performance of our method, mi-mfa was also performed on the nci- <dig> dataset.
a pattern of missing values was created as illustrated in table  <dig>  one or two individuals were removed per table for all types of cancer cell lines, except for the prostate cancer line .
this specific pattern would reflect a study in which a lot of rows were missing.
number of missing rows inserted in each stratum  of the original dataset, including both transcriptomic and proteomic data, for incomplete data creation.
the # cases indicate the number of possibilities of incomplete cases per stratum
to compare the mfa configurations, one incomplete dataset was chosen from a large range of possibilities of incomplete datasets  according to the scenario of missingness illustrated in table  <dig>  we then computed our mi-mfa method on this incomplete dataset by using m= <dig> imputed datasets.
as with the liver toxicity data, ri-mfa and mvi-mfa were also performed.
we chose m= <dig> in order to achieve stable results.
figure  <dig> shows the divergence of the mi-, ri- and mvi-mfa configurations from the true configuration.
for this specific example, the mi-mfa configuration was closest to the true configuration  with a rv coefficient of  <dig> , whereas the configurations obtained with ri-mfa and mvi-mfa were more distorted with rv coefficients of  <dig>  and  <dig>  respectively .
fig.
<dig> comparison of mfa configurations for nci- <dig> data analysis.
representation of the individuals on the two-dimensional configuration obtained by mfa on the complete dataset  and, mi-mfa , ri-mfa  and mvi-mfa  versus the true configuration.
dots and arrows represent the projected coordinates from the true configuration and the imputation mfa method respectively.
the length of the line joining dots and arrows is proportional to the divergence between the projected coordinates.
each individual  is colored according to the cancer type .
imputed individuals from transcriptomic and proteomic data are represented by empty circles and diamonds, respectively
to broaden the scope of assessment to more than a single case,  <dig> possible cases of incomplete datasets were chosen at random among the 6× <dig> possibilities according to the specific missingness scenario shown in table  <dig>  for each incomplete dataset, mi-mfa was performed using m= <dig> imputed datasets.
ri-mfa and mvi-mfa were also computed.
figure  <dig> shows the rv coefficient for the two-dimensional configurations as a function of each case.
the results obtained with mi-mfa and ri-mfa were similar, with rv coefficients of approx.
<dig> , whereas the results obtained with mvi-mfa were much further from the true configuration.
thus, even with complex patterns of missingness, the mi-mfa approach still provided satisfactory results, as did ri-mfa in this case.
however, unlike ri-mfa, mi-mfa took into account the variability of missing row values, as demonstrated in the following section.
fig.
<dig> performance on nci- <dig> data.
rv coefficient between the configuration obtained by mfa on the complete dataset and either mi-mfa , ri-mfa  or mvi-mfa  on the incomplete dataset.
the results shown are the mean rv coefficients for the  <dig> two-dimensional configurations as a function of the cases according to the scenario represented in table  <dig>  the discrete rv values are joined by lines for ease of understanding
why is it essential to evaluate uncertainty?
this question was addressed through an example using the nci- <dig> dataset.
a specific pattern of missing values was created.
the missing rows were randomly introduced for four melanoma and two leukemia cancer lines in the transcriptomic table.
the six inserted missing rows represented  <dig> % of the total number of individuals.
the missing rows were inserted for specific groups of individuals that contributed substantially to the construction of the first two dimensions of the mfa on the original dataset .
mi-mfa was performed on the incomplete dataset using m= <dig> imputed datasets.
confidence ellipses and convex hulls were then computed from the  <dig> configurations projected on the compromise configuration.
figure  <dig>  shows the uncertainty due to missing rows around individuals on the compromise configuration.
the use of different imputed individuals in each dataset implied slightly different configurations.
consequently, since the configurations changed, the positions of all the individuals also changed and thus all the individuals had confidence areas, even individuals not being imputed .
however, the greatest uncertainty occurred around the imputed individuals .
fig.
<dig> visualization of the uncertainty induced by the missing individuals.
projection of the  <dig> obtained configurations on the compromise configuration from mi-mfa on the incomplete data with  <dig> % missing rows  and  <dig> % missing rows .
the  <dig> % confidence ellipses  and convex hulls  show the uncertainty for each individual .
empty circles represent imputed individuals from the transcriptomic data.
for ease of understanding, not all individuals for the  <dig> configurations obtained were plotted
in order to highlight the importance of the uncertainty of mi-mfa configurations induced by missing rows, additional rows were removed from the transcriptomic table resulting in  <dig> % missing rows.
we then carried out mi-mfa on the incomplete dataset using m= <dig> imputed datasets.
figure  <dig>  shows the impact of the missing rows around individuals on the compromise configuration.
as expected, the size of the ellipses  of the additional missing individuals was increased.
however, the size of the ellipses and convex hulls was not excessive even when  <dig> % of the rows were missing from the transcriptomic table.
discussion
in the present paper, we propose a mi approach to handle missing row values in an mfa framework and therefore resolve one of the major issues associated with multiple omics data tables.
the proposed method, which we called mi-mfa, provides point estimates of mfa configurations with a notion of uncertainty due to the missing values.
the aim of our method was not to obtain the best possible estimates for the missing values, but to replace them so as to be able to estimate mfa configurations.
the mi-mfa method generates m imputed datasets from an mfa model, where multiple hot-deck imputation is used to fill in the missing values.
the hot-deck approach resolves the most important limitation of other model-based techniques  in that it can be applied to large datasets.
furthermore, other major advantages of this method are that:  it is not necessary to define an explicit model for the distribution of the missing values,  imputations tend to be realistic since they are based on observed values, and  it is flexible in the sense that it can preserve complex within-unit and between-variable associations.
however, a weakness is that it requires good donor-recipient matches that reflect the available covariate information.
finding such good matches is not an easy task and is beyond the scope of this article  <cit> .
in an ideal framework of stratified multiple omics data tables, donor pools would consist of the available individuals belonging to the same stratum and the same omics table as the recipient.
a potential shortcoming of the method is that the donor pools might contain too few donor observations and thus introduce a risk of bias on the mfa results.
likewise, if the overall sample size is very small, then typically there are also too few potential donors.
however, all imputation techniques are challenged by small sample sizes since these reduce the availability of information required to create suitable conditional statements  <cit> .
an important aspect of our strategy is the choice of the number m of imputed datasets.
as shown in the two case studies, we think that this number should be a good compromise between the need to obtain stable estimates and to avoid computation bottlenecks.
the statis method was proposed to combine the results of mi-mfa.
as the core of mfa is a pca, combining the results from mfa is the same as combining the results from pca.
several procedures have previously been proposed to combine results from pca such as the mean varimax method  or mean correlation matrix  approaches .
however, van ginkel and kroonenberg  <cit>  demonstrated that generalized procrustes analysis  was more suited to this purpose.
gpa fits the pca configurations obtained from the imputed datasets to a single fixed reference configuration to produce the final solution, the centroid configuration .
one advantage of statis, as compared to mvm and mcm, is that it automatically corrects for possible reflection, dilation or rotation of the different configurations.
additionally, and contrary to the gpa procedure, the statis algorithm is very computer-time efficient since it is a non-iterative process.
another appealing feature of statis is its robust properties.
as this algorithm includes weights proportional to the agreement between configurations, the results do not seem to be affected by the presence of large outliers.
two approaches have been proposed to visualize the uncertainty of the estimated mfa configurations due to missing row values: confidence ellipses and convex hulls.
these graphical representations provide scientists with considerable guidance when interpreting the significance of mfa results in a missing data framework.
indeed, ellipses and convex hull areas offer great assistance by either supporting the mfa results if they are small or suggesting that caution be exercised otherwise.
it should be noted that ri-mfa  also provides a configuration of individuals whatever the missingness pattern; however, there is no way of knowing if the results obtained are plausible and if the user can interpret the results without making any mistakes.
we have illustrated our approach by applying it to two real case studies using the liver toxicity and nci- <dig> datasets.
incomplete artificial datasets with different patterns of missingness were created within these datasets.
the configurations resulting from mi-, ri- and mvi-mfa were compared with the mfa configurations of the corresponding original population .
performance of the methods was assessed by considering the rv coefficient with respect to the true configuration.
in the liver toxicity study, the performances of the methods were compared in two different missingness settings.
first the number of missing rows in each stratum of the transcriptomic table was chosen to be low , medium  or high .
secondly, seven scenarios were created by inserting missing rows in the original dataset which included both transcriptomic and clinical tables.
the overall results showed that mi-mfa clearly outperformed the ri-mfa and mvi-mfa approaches in nearly all settings.
in the nci- <dig> study, we illustrated the performance of our method on complex patterns of missingness where substantial numbers of rows were missing from both tables of the nci- <dig> dataset.
as previously, this study showed that mi-mfa clearly performed better than mvi-mfa.
we also demonstrated that ri-mfa performed better than mvi-mfa.
the differences between mi-mfa and ri-mfa were small, but on average slightly in favor of our method.
as the purpose of this study was also to illustrate the uncertainty of mi-mfa configurations induced by missing rows, we demonstrated how the areas of the confidence ellipses and convex hulls got larger as the number of missing rows increased.
CONCLUSIONS
we propose here a new method, mi-mfa, an extension of mfa, designed to deal with multiple tables with missing row values.
mi-mfa is a useful and attractive method to estimate the coordinates of individuals for mfa configurations despite the missing rows.
the study cases showed that the other proposed methods either encountered serious problems or were unable to adequately assess the accuracy due to missing data.
the configurations obtained with our method were closer to the true configuration even when a significant number of individuals were missing, and thus provided better results.
moreover, the uncertainty due to the missing rows could be visualized on the compromise configuration.
the software for our mi-mfa method is available in an easy-to-use code for the r statistical environment.
additional files
additional file  <dig> supplementary figures.
figures s1–s <dig>
additional file  <dig> calculation of the total number of possible imputations in mi-mfa.
additional file  <dig> r code implementing the mi-mfa method.
r  is required.
abbreviations
fcsfully conditional specification
jmjoint modeling
marmissing at random
mcarmissing completely at random
mfamultiple factor analysis
mimultiple imputation
micemultivariate imputation by chained equations
mi-mfamultiple imputation multiple factor analysis
mvi-mfamean variable imputation multiple factor analysis
pcaprincipal component analysis
ri-mfaregularized iterative multiple factor analysis
statisstructuration des tableaux à trois indices de la statistique

 due to the increase in available data information, integrating large amounts of heterogeneous data is currently one of the major challenges in systems biology. biological data integration provides scientists with a deeper insight into complex biological processes. however, when dealing with multiple data tables, the presence of missing values is a common situation for a variety of reasons. in omics data integration studies, it is common for some individuals to not be present in all data tables, resulting in a specific missing data pattern for multiple tables, as shown in fig. for instance, in clinical studies, this can occur when a patient forgets to fill out a form. it also can be attributable to the study design if individual data are expensive or difficult to measure. fig. pattern of missingness. pattern of missing values specific to multiple data tables with missing rows for some tables of variables. the colored parts represent strata of the observed data, whereas the empty parts indicate missing data missing row values for a table of variables are challenging to handle because most statistical methods cannot be directly applied to incomplete datasets. in the multiple multivariate framework, several approaches have already been proposed to deal with missing row values. the only methods widely available for analyzing incomplete data focus on removal of the missing rows, either by ignoring subjects with incomplete information or by replacing the missing items with plausible values. in multivariate statistical analysis, case deletion procedures can be very inefficient, discarding an unacceptably high proportion of subjects because even if the per-table rates of missing rows are low, only a few subjects may have complete data for all tables. in addition, case-deletion procedures may bias the results if the remaining subjects providing the complete data are unrepresentative of the entire sample. on the other hand, simple mean substitution can seriously distort the marginal and joint distribution of the variables and be an issue because many statistical methods rely on estimation of the variance-covariance matrix. recently, two approaches have been proposed to deal with missing row values in multiple multivariate analysis. the first method, introduced by van de velden and bijmolt, was developed in the context of generalized canonical correlation analysis. its application in the omics framework is often limited by the size, noise and multicollinearity of the data. the second method, described in husson and josse, was developed in the context of multiple factor analysis. this method, designated regularized iterative mfa, was derived from a method available in principal component analysis and consists of alternating the estimation of axes and components, and the estimation of missing values. here we consider an alternative method, involving a multiple imputation approach adapted to the mfa framework, and called mi-mfa. multiple imputation was proposed by rubin in order to estimate both the parameters of interest and their variability in a data missingness framework. it relies on the principle that a single value cannot reflect the uncertainty of the estimation of a missing value. first, mi is used to generate plausible synthetic data values, called imputations, for missing values in the data. this step results in a number of imputed datasets in which the missing data are replaced by random draws of plausible values according to a specific statistical model. the second step consists of analyzing each imputed dataset using a statistical method that estimates the parameters of interest. this step results in m analyses which differ only because the imputations differ. finally, mi combines all the results together to obtain a single consensus estimate, thereby combining variation within and across the m imputed datasets. under fairly liberal conditions, this last step results in statistically valid estimates that properly reflect sampling variability. the major challenge in mi involves generating possible values for each missing observation. statistically advanced imputation procedures can therefore be used for this. two general approaches are often used for imputing multivariate data: joint modeling and fully conditional specification, also known as multivariate imputation by chained equations. jm involves specifying a multivariate distribution for the missing data and drawing imputations from their conditional distributions by markov chain monte carlo techniques. fcs specifies the multivariate imputation model on a variable-by-variable basis using a set of conditional models, one for each incomplete variable. the key issue in jm is appropriate specification of the multivariate distribution. a multivariate normal model has often been used as it is computationally tractable. this model has even been used when some of the variables are not gaussian. however, the main weakness of jm is that it can only be applied when the imputation involves a small number of variables. this is not very common in omics datasets that are often composed of tens of thousands of variables, or more. fcs allows greater flexibility than jm in creating multivariate models. indeed, fcs can use specialized imputation models by separately defining the conditional densities for each variable, even if this can require a considerable amount of work. when the number of variables is large, it is often impractical and too computer-intensive to test and develop the best models for each variable. as an alternative to the jm and fcs approaches, we propose using the hot-deck imputation approach. this approach is a nonparametric imputation method that resolves the most important limitation of the jm and fcs approaches as it can be applied to data tables containing more than just a few variables. when using the mi method, special attention must be given to the process that gave rise to the missing data, referred to as the missing data mechanism. most methods for generating multiple imputations, fully-, semi- and non-parametric methods, assume that the mechanism responsible for missing data is ignorable. briefly, if the missing data mechanism is ignorable, then the analysis can focus on the observed values rather that also having to model the process that resulted in certain values being observed and certain values being missing. if the assumption of an ignorable missing-data mechanism is valid, then statistical methods that rely on that assumption can be expected to produce results with minimal bias. one way that the missing-data mechanism can be viewed as ignorable is if the missing data are missing completely at random. for data to be mcar, there must not be any systematic differences between the cases that have missing items and the cases that are fully observed. in microarray experiments, technical failure, low signal-to-noise ratio and measurement errors can for instance be considered as sources of mcar patterns. the missing-data mechanism can also be viewed as ignorable under the less restrictive missing at random scenario, which allows missingness to depend on observed variables but not on unobserved variables. late-stage cancer patients, as compared to early-stage cancer patients, unfortunately have more chance of dropping out of follow-up studies, which may result in a mar pattern in a clinical data table. when the ignorability assumption does not hold, the imputation needs to be drawn from the posterior distribution of the missing data given the complete data and the missingness mechanism. non-ignorable missing data occurs frequently in mass-spectrometry-based experiments. measures too close to the limit of detection of the instrument are censored, resulting in a higher rate of missing values. the probability of being missing is, in this particular case, directly dependent on the intensity value. in this paper, we decided to focus on models with ignorable missing-data mechanisms. methods mathematical basis of multiple factor analysis mfa is devoted to the simultaneous exploration of multiple data tables where the same individuals are described by several tables of variables. in mfa, the number and the type of variables may vary from one table to another, but within each table the nature of the variables is the same. here we focus on quantitative variables. the aims of mfa are similar to those of pca, namely to study the similarities between individuals from a multidimensional point of view, to analyze the relationships between variables and characterize individuals based on these relationships. however, beyond these conventional uses, mfa can also be used to study the links between tables of variables and to compare the information contributed by each table. mfa analyzes a set of j data tables k …,kj, where each kj corresponds to a table of quantitative variables measured on the same i individuals. the core of mfa is a pca in which weights are assigned to variables. more formally, the matrix of variance-covariance associated with each data table kj is decomposed by pca and its largest eigenvalue λ1j is derived. then, each variable belonging to kj is weighted by 1/λ1j. finally, a global pca is performed on the merged and weighted data table k= to obtain the configuration f. the main reason for the weighting step is to remove from each table all information related to its own dimensionality or variance. therefore, no single table can dominate the first dimension of the global analysis. mfa provides the same graphical representations as pca but also, due to the table structure, specific representations such as the table representation and the superimposed representation are available. the multiple imputation multiple factor analysis approach to deal with multiple tables with missing rows, we propose the mi-mfa approach, a multiple imputation adapted to the framework of mfa. the aim of our method is not to get the best possible estimations of the missing values, but to replace them with plausible values in order to provide estimates of the mfa configurations. according to mi methodology, the mi-mfa approach is carried out by performing the following three steps: imputation: generate m different imputed datasets k,…,k,…,k of k. mfa analysis: perform an mfa on each k imputed dataset leading to m different configurations f …,fm,…,fm. combination: find a consensus configuration between all f …,fm configurations. these steps are outlined in fig. and described in detail in the following sections. fig. overview of our mi-mfa approach to handling missing rows in multi-omics data integration. the top part of the graphic indicates that analysis starts with observed, incomplete data tables k. in a second step, multiple imputation is performed using the hot-deck imputation approach: m imputed versions k,…,k of k are obtained by replacing the missing values by plausible data values. these plausible values are drawn from donor pools. the imputed sets are identical for the non-missing data entries, but differ in the imputed values. the third step is to estimate the configuration matrix f m for each imputed dataset k using mfa. the estimated configurations differ from each other because their input data differ. the last step is to combine the m estimated configurations f …,f m into a compromise configuration f c using the statis method generating imputed data: multiple hot-deck imputation hot-deck imputation involves replacing missing values of one or more variables with available values from a similar unit. the observation from which these available values are taken for imputation is called the donor and the observation with the missing value, which receives the donor’s value, is the recipient. the donor can be randomly selected from a set of potential donors, called the donor pool. selection of a suitable donor pool is not an easy task and is beyond the scope of this article. the general principle is to choose donor units that are as close as possible to the recipient with respect to some affinity score. affinity is defined in terms of the degree to which each potential donor matches the recipient’s values across all variables other than the one being imputed. intuitively, in the framework of stratified multiple omics tables, the donor pool can be formed of available individuals belonging to the same stratum and the same omics table as the recipient. multiple hot decking differs from other forms of hot decking by using several donors for a single recipient. multiple hot-deck imputation proceeds as follows. let k= be the merged data table containing missing rows with strata s= …,s, then carry out the following steps: step create donor pools by taking donors belonging to the same stratum s and the same table kj as the recipient. recipients within the same stratum have the same donor pool. suppose always that there is a large enough number of donors for recipients in each stratum. step for each recipient in kj, impute the missing individual by drawing randomly with replacement a donor from the corresponding donor pool. repeat this procedure until all missing individuals in the j tables have been imputed. step repeat step until m different imputed datasets k,…,k of k are obtained. by conducting the imputations in this way, it is reasonable to assume that the within-unit between-variables multivariate relationships are preserved. the combination procedure: the statis method the question that arises after using mi in an mfa framework is how should all the configurations resulting from the analyses be combined to obtain a single unique estimate of the consensus configuration? while averaging is an appropriate combination procedure in many other statistical techniques, it is not recommended for mfa due to possible reflection, dilation or rotation of the different configurations with respect to each other. here we consider an alternative approach by implementing the statis method which provides a compromise configuration balancing all configurations. the statis method is a generalization of pca used to simultaneously study several tables of variables collected on the same individuals. the goal of this method is to analyze the structure of the individual tables and to derive from this structure an optimal set of weights for computing a common configuration of the observations. the solution obtained, called the compromise, is the configuration agrees the most with all other configurations. an overview of the statis method is presented in additional file 1: figure s statis analyzes a set of n tables x …,xn, where each xn is a table of quantitative variables measured on the same individuals. the first stage of statis consists of calculating a matrix of cross-products between individuals for each table wn=xnxnt reflecting the similarities between individuals within this table. the use of matrices wn instead of xn simplifies the computation because it obviates the determination of rotations when matching the xn. the basic idea in statis is then to find a compromise space wc=∑n=1nαnwn that globally balances these cross-product matrices by choosing a suitable optimal set of weights α …,αn. these weights are obtained from the pca of the matrix r whose generic term rjk gives the cosine between tables defined as: rjk=tracetrace·trace, where the trace is the sum of the main diagonal elements of a square matrix. the first eigenvector obtained from the pca of r represents the “agreement between tables”. its elements are normalized in such a way that their sum is equal to and used as weights αn in order to define wc. tables with larger values of αn are more similar to the other tables and therefore will have a larger weight, while the weight of the “outlier tables” will be closer to zero with respect to the other weights. the principal components from the pca of wc then gives the coordinates of the individuals in the compromise space, called the compromise configuration. implementation of mi-mfa the mi-mfa algorithm can be summarized as follows: step start with an observed, incomplete dataset k. define the number of imputations m and the dimensionality d of the compromise configuration. step perform multiple hot-deck imputation. for m= …,m: obtain an imputed version k of k, such that, k≠k for m≠m′. the imputed datasets are identical for the non-missing data entries, but differ for the imputed values. the imputed version of the data is obtained by using the hot-deck imputation approach. perform an mfa using d components on the imputed dataset k to obtain the configuration fm. step perform a statis on the set of configurations f …,fm to obtain fc, the compromise configuration. note that the number of dimensions d used in the algorithm has to be chosen a priori. however, the number of dimensions does not affect the estimation of the imputed values and the estimation of the compromise configuration. moreover, for given k,…,k imputed datasets, solutions provided by the algorithms are nested. since the core of mi-mfa is a weighted pca, the strategies suggested to choose the number of components in pca can be adapted to mi-mfa, but work needs to be done to validate the quality of these extensions. how many imputations? when using mi, one of the uncertainties concerns the number m of imputed datasets needed to obtain satisfactory results. the number of imputed datasets in mi depends to a large extent on the proportion of missing data. the greater the missingness, the larger the number of imputations needed to obtain stable results. however, in multiple hot-deck imputation, the number of imputed datasets is limited by the size of the donor pools. in any case, the total number of possible imputations mtotal can be calculated before applying the imputation approach. if mtotal is small, then m=mtotal can be used in mi-mfa. the appropriate number of imputations can be informally determined by carrying out mi-mfa on n replicate sets of ml imputations for l=, …, with m0<m1<m2<⋯<mtotal, until the estimate compromise configurations are stabilized. more precisely, this approach can be carried out by applying the following steps: step start with an observed, incomplete dataset k. define the number of imputations ml with m0<m1<m2<⋯<mtotal and the number n of replicate sets of ml imputations. step create collections inml, n= …,n, each one containing ml different imputed datasets of k, such that, inml≠in′ml, for n≠n′ and inml−1⊂inml for m0<m1<m2<⋯<mtotal. step for n= …,n, perform an mi-mfa using inm to obtain n different compromise configurations fc1m …,fcnm step let l= for n= …,n, perform an mi-mfa using the collection inml, to obtain a compromise configuration fcnml; calculate rnl=r, a measure of the distance or correlation between configurations. step calculate r¯l=1n∑n=1nrnl the standard error of r¯l). step repeat steps to for l=,… until the differences between two subsequent r¯l ) become smaller than a certain convergence criterion. uncertainty of mi-mfa solutions in an mi-mfa framework, after estimating the configurations from the imputed datasets, a new source of variability due to missing values can be taken into account. here we describe two approaches to visualize the uncertainty of the estimated mfa configurations attributable to missing row values. first, an individual plot for all estimated mfa configurations is constructed. the individual plot is obtained by projecting each estimated mfa configuration onto the compromise configuration. each individual is represented by m points, each corresponding to one of the m mfa configurations. confidence ellipses and convex hulls can then be constructed for the m configurations for each individual. the computed convex hull results in a polygon containing all m solutions. all individuals have confidence areas, even those without missing values. indeed, even if only the estimation of missing values is the only change, this will have a possible impact on all mfa parameters. therefore, the area of such an ellipse provides an insight into the uncertainty of the estimated configuration. the larger the area of an ellipse, the more uncertain the exact location of the individual. thus, when the area of an ellipse is large, the scientist should remain really careful regarding its interpretation. performance of the method we conducted two case studies to assess the performance of our method. instead of using theoretical distributions to generate simulated data, our studies were based on two real datasets, denoted as the original datasets. subsequently, specific patterns of missingness were created in these datasets as illustrated in fig. resulting in what we called the incomplete datasets. this approach was used in order to more closely mirror situations that may occur in the omics context. next, missing row values were estimated and the resulting complete datasets were referred to the imputed datasets. we then compared our mi-mfa method to the ri-mfa approach and the mean variable imputation mfa method, in which the missing values are simply replaced by the mean of each variable after which an mfa is performed on the imputed dataset. this latter approach was considered as the common base for comparing the mi-mfa and ri-mfa methods. for each configuration obtained using mi-, ri- and mvi-mfa, the similarity between the configuration solution and the true configuration was assessed from the rv coefficient. the rv coefficient, which ranges from zero to one, can be interpreted as a correlation coefficient between two matrices, which allows the relative positions of objects to be compared from one configuration to another. we also provide comprehensive graphical comparisons of the true vs. the mi-, ri- or mvi- mfa configurations. the individuals from both configurations are drawn in a same plot and connected by an arrow, the length of which indicates the divergence between the two configurations. implementation of the analyses all analyses were performed using the r computing environment. mfa was performed using the mfa function of the factominer r package. the statis function of the ade r package was used to determine the compromise configuration. the ri-mfa method is implemented in the imputemfa function available in the missmda r package. note that the number of components ncp used to predict the missing entries in the imputemfa function has to be chosen a priori. this choice is crucial and difficult. as the true configuration was known in our case, the number of components ncp was chosen to minimize the rv coefficient between the true and the imputemfa configurations. the appropriateness of the results from mi-, ri- and mvi- mfa was then determined by comparing the configurations resulting from these three strategies with the true mfa configuration. due to a possible lack of alignment between two configurations, it was necessary to align them before being compared. ordinary procrustes analysis was used to align these configurations prior to their comparison. datasets liver toxicity the datasets originated from a liver toxicity study in which male rats of the inbred fisher f344/n strain were exposed to toxic doses of acetaminophen in a controlled experiment. necropsies were performed and h after exposure and mrna was extracted. the data consisted of the expression of genes and clinical variables considered to be markers of liver injury. the subjects were cross-classified in eight strata according to two factors: exposure time: and h; toxic doses of acetaminophen: high or low. eight subjects per stratum were included. these datasets were downloaded from the mixomics r package. nci- data the nci- dataset contained transcriptomic and proteomic tables for a collection of cell lines from the national cancer institute. the nci- panel included cell lines derived from various cancer types: colon, renal, ovarian, breast, prostate, lung and central nervous system, as well as leukemia and melanoma. the gene expression profiles used here were generated using an agilent platform and downloaded from cellminer. data were log2-transformed. to facilitate data interpretation and computation, the transcriptomic data were filtered to exclude probes that did not map to an official hugo gene symbol and to retain only the probe with the highest average value when multiple probes mapped to the same gene, as previously described in. gene invariants across all cell lines, corresponding to genes without any effect between cancer types, were also removed. filtering produced a dataset of genes. the nci- proteome table was also downloaded from cellminer. proteomic data were obtained using high-density reverse-phase lysate microarrays. data were log2-transformed and protein abundance levels were available for proteins.  liver toxicity data analysis a specific pattern of missing values was created as illustrated in fig. to obtain an incomplete dataset, three individuals per stratum were randomly removed from the transcriptomic table. for this specific pattern, there were 3×8= missing individuals. mi-mfa was then performed on the incomplete dataset using m= imputed datasets. ri-mfa, and mvi-mfa were also performed. figure shows the divergence of the mi-, ri- and mvi-mfa configurations from the true configuration. as can be seen, the configuration obtained with mi-mfa was very close to the true configuration. this result was confirmed by the high rv coefficient. the configurations obtained with ri-mfa and mvi-mfa were more distorted and less close to the true configuration with rv coefficients of and respectively. fig. comparison of mfa configurations for the liver toxicity data analysis. representation of the individuals on the two-dimensional configuration obtained by performing mfa on the complete dataset and, mi-mfa, ri-mfa and mvi-mfa versus the true configuration. dots and arrows represent the projected coordinates from the true configuration and the imputation mfa method, respectively. the length of the line joining dots and arrows is proportional to the divergence between the projected coordinates. the color of each individual reflects the treatment. imputed individuals from the transcriptomic dataset are represented by empty circles the number m of imputed datasets in mi-mfa for the above incomplete liver toxicity data was determined as described in the how many imputations section. collections of size n= were generated for each of the following numbers of imputations: ml= l, for l= …, the stability of the estimated mi-mfa configurations was then determined by calculating the rv coefficient between the configurations obtained using ml and ml+ imputations. as the true configuration was known, we also described the stability of the estimated mi-mfa configurations by calculating the rv coefficients between the true configuration and those obtained using ml imputations. although the missing information is substantial, fig. shows that only a slight increase in precision was obtained by using more than imputations. fig. stability of the estimated mi-mfa configurations using different numbers of imputations. mi-mfa configurations were obtained for the following numbers of imputations: m l= l, for l= …, left: rv coefficient between configurations obtained by mi-mfa with m l and m l+ imputations. right: rv coefficient between the configuration obtained by mfa on the complete dataset and mi-mfa on the incomplete dataset with m l imputations. the values shown are the mean rv coefficients for the n= two-dimensional configurations as a function of the number of imputations. error bars represent the standard deviation of the rv coefficients mi-mfa was also applied to different scenarios of missingness. first, mfa was performed on the original dataset to obtain the true configuration. secondly, individuals were randomly removed from each stratum in the original transcriptomic table. for this, three scenarios were considered for the number of missing rows, in which there were respectively one, two and three missing rows per stratum. the total number of missing rows per scenario was therefore and respectively. thirdly, for each missingness scenario, incomplete datasets were randomly chosen for analysis. as previously, mi-mfa was performed on each incomplete dataset using m= imputed datasets. ri-mfa, as well as mvi-mfa, were also computed. the rv coefficients between the true configuration and the configurations obtained using each method were then calculated. figure shows the mean of the rv coefficients for the two-dimensional configurations as a function of the missingness scenario for each method. note that the average results using mi-mfa were always better than with ri-mfa or mvi-mfa, whatever the scenario. the rv coefficients between the true configuration and the mi-mfa configuration were close to one and remained satisfactory even when the number of missing row values was high, and the results obtained with the ri-mfa and mvi-mfa decreased significantly. fig. performance on liver toxicity data according to the number of missing rows. rv coefficients between the configuration obtained by mfa on the complete dataset and either mi-mfa, ri-mfa or mvi-mfa on the incomplete dataset. the values shown are the mean rv coefficients for the two-dimensional configurations for each missingness scenario: low, medium and high. error bars represent the standard deviation of the rv coefficients the performance of the mi-mfa procedure was then further investigated for more complex scenarios of missingness. more precisely, missing row values were inserted into each stratum of the original dataset according to the scenarios illustrated in table number of missing rows inserted in each stratum of the original dataset, including both transcriptomic and clinical data, for incomplete data creation. the # cases indicate the number of possibilities of incomplete cases per stratum twenty incomplete datasets were then selected at random from each stratum and each scenario. all analyses and the calculation of the rv coefficient were performed in the same way as previously described. figure shows the means of rv coefficients for the two-dimensional configurations as a function of the scenarios for each method. fig. performance on liver toxicity data with missing individuals in the high-dose treatment. average rv coefficient between the configuration obtained by mfa on the complete dataset and either mi-mfa, ri-mfa or mvi-mfa on the incomplete dataset. results are given for all of the first three two-dimensional possible configurations as a function of the scenarios presented in table the discrete rv values are joined by lines for ease of understanding. the performance of the different methods on the liver toxicity data with missing individuals in the low-dose treatment are presented in additional file 1: figure s3 for almost all the scenarios, the average results obtained with mi-mfa were better than with the other methods. as the number of missing row values increased, the results obtained with the ri-mfa and mvi-mfa algorithms worsened rapidly, especially when the dimension increased, whereas the results with the mi-mfa approach were still satisfactory. nci- data analysis to confirm the performance of our method, mi-mfa was also performed on the nci- dataset. a pattern of missing values was created as illustrated in table one or two individuals were removed per table for all types of cancer cell lines, except for the prostate cancer line. this specific pattern would reflect a study in which a lot of rows were missing. number of missing rows inserted in each stratum of the original dataset, including both transcriptomic and proteomic data, for incomplete data creation. the # cases indicate the number of possibilities of incomplete cases per stratum to compare the mfa configurations, one incomplete dataset was chosen from a large range of possibilities of incomplete datasets according to the scenario of missingness illustrated in table we then computed our mi-mfa method on this incomplete dataset by using m= imputed datasets. as with the liver toxicity data, ri-mfa and mvi-mfa were also performed. we chose m= in order to achieve stable results. figure shows the divergence of the mi-, ri- and mvi-mfa configurations from the true configuration. for this specific example, the mi-mfa configuration was closest to the true configuration with a rv coefficient of, whereas the configurations obtained with ri-mfa and mvi-mfa were more distorted with rv coefficients of and respectively. fig. comparison of mfa configurations for nci- data analysis. representation of the individuals on the two-dimensional configuration obtained by mfa on the complete dataset and, mi-mfa, ri-mfa and mvi-mfa versus the true configuration. dots and arrows represent the projected coordinates from the true configuration and the imputation mfa method respectively. the length of the line joining dots and arrows is proportional to the divergence between the projected coordinates. each individual is colored according to the cancer type. imputed individuals from transcriptomic and proteomic data are represented by empty circles and diamonds, respectively to broaden the scope of assessment to more than a single case, possible cases of incomplete datasets were chosen at random among the 6× possibilities according to the specific missingness scenario shown in table for each incomplete dataset, mi-mfa was performed using m= imputed datasets. ri-mfa and mvi-mfa were also computed. figure shows the rv coefficient for the two-dimensional configurations as a function of each case. the results obtained with mi-mfa and ri-mfa were similar, with rv coefficients of approx., whereas the results obtained with mvi-mfa were much further from the true configuration. thus, even with complex patterns of missingness, the mi-mfa approach still provided satisfactory results, as did ri-mfa in this case. however, unlike ri-mfa, mi-mfa took into account the variability of missing row values, as demonstrated in the following section. fig. performance on nci- data. rv coefficient between the configuration obtained by mfa on the complete dataset and either mi-mfa, ri-mfa or mvi-mfa on the incomplete dataset. the results shown are the mean rv coefficients for the two-dimensional configurations as a function of the cases according to the scenario represented in table the discrete rv values are joined by lines for ease of understanding why is it essential to evaluate uncertainty? this question was addressed through an example using the nci- dataset. a specific pattern of missing values was created. the missing rows were randomly introduced for four melanoma and two leukemia cancer lines in the transcriptomic table. the six inserted missing rows represented % of the total number of individuals. the missing rows were inserted for specific groups of individuals that contributed substantially to the construction of the first two dimensions of the mfa on the original dataset. mi-mfa was performed on the incomplete dataset using m= imputed datasets. confidence ellipses and convex hulls were then computed from the configurations projected on the compromise configuration. figure shows the uncertainty due to missing rows around individuals on the compromise configuration. the use of different imputed individuals in each dataset implied slightly different configurations. consequently, since the configurations changed, the positions of all the individuals also changed and thus all the individuals had confidence areas, even individuals not being imputed. however, the greatest uncertainty occurred around the imputed individuals. fig. visualization of the uncertainty induced by the missing individuals. projection of the obtained configurations on the compromise configuration from mi-mfa on the incomplete data with % missing rows and % missing rows. the % confidence ellipses and convex hulls show the uncertainty for each individual. empty circles represent imputed individuals from the transcriptomic data. for ease of understanding, not all individuals for the configurations obtained were plotted in order to highlight the importance of the uncertainty of mi-mfa configurations induced by missing rows, additional rows were removed from the transcriptomic table resulting in % missing rows. we then carried out mi-mfa on the incomplete dataset using m= imputed datasets. figure shows the impact of the missing rows around individuals on the compromise configuration. as expected, the size of the ellipses of the additional missing individuals was increased. however, the size of the ellipses and convex hulls was not excessive even when % of the rows were missing from the transcriptomic table. discussion in the present paper, we propose a mi approach to handle missing row values in an mfa framework and therefore resolve one of the major issues associated with multiple omics data tables. the proposed method, which we called mi-mfa, provides point estimates of mfa configurations with a notion of uncertainty due to the missing values. the aim of our method was not to obtain the best possible estimates for the missing values, but to replace them so as to be able to estimate mfa configurations. the mi-mfa method generates m imputed datasets from an mfa model, where multiple hot-deck imputation is used to fill in the missing values. the hot-deck approach resolves the most important limitation of other model-based techniques in that it can be applied to large datasets. furthermore, other major advantages of this method are that: it is not necessary to define an explicit model for the distribution of the missing values, imputations tend to be realistic since they are based on observed values, and it is flexible in the sense that it can preserve complex within-unit and between-variable associations. however, a weakness is that it requires good donor-recipient matches that reflect the available covariate information. finding such good matches is not an easy task and is beyond the scope of this article. in an ideal framework of stratified multiple omics data tables, donor pools would consist of the available individuals belonging to the same stratum and the same omics table as the recipient. a potential shortcoming of the method is that the donor pools might contain too few donor observations and thus introduce a risk of bias on the mfa results. likewise, if the overall sample size is very small, then typically there are also too few potential donors. however, all imputation techniques are challenged by small sample sizes since these reduce the availability of information required to create suitable conditional statements. an important aspect of our strategy is the choice of the number m of imputed datasets. as shown in the two case studies, we think that this number should be a good compromise between the need to obtain stable estimates and to avoid computation bottlenecks. the statis method was proposed to combine the results of mi-mfa. as the core of mfa is a pca, combining the results from mfa is the same as combining the results from pca. several procedures have previously been proposed to combine results from pca such as the mean varimax method or mean correlation matrix approaches. however, van ginkel and kroonenberg demonstrated that generalized procrustes analysis was more suited to this purpose. gpa fits the pca configurations obtained from the imputed datasets to a single fixed reference configuration to produce the final solution, the centroid configuration. one advantage of statis, as compared to mvm and mcm, is that it automatically corrects for possible reflection, dilation or rotation of the different configurations. additionally, and contrary to the gpa procedure, the statis algorithm is very computer-time efficient since it is a non-iterative process. another appealing feature of statis is its robust properties. as this algorithm includes weights proportional to the agreement between configurations, the results do not seem to be affected by the presence of large outliers. two approaches have been proposed to visualize the uncertainty of the estimated mfa configurations due to missing row values: confidence ellipses and convex hulls. these graphical representations provide scientists with considerable guidance when interpreting the significance of mfa results in a missing data framework. indeed, ellipses and convex hull areas offer great assistance by either supporting the mfa results if they are small or suggesting that caution be exercised otherwise. it should be noted that ri-mfa also provides a configuration of individuals whatever the missingness pattern; however, there is no way of knowing if the results obtained are plausible and if the user can interpret the results without making any mistakes. we have illustrated our approach by applying it to two real case studies using the liver toxicity and nci- datasets. incomplete artificial datasets with different patterns of missingness were created within these datasets. the configurations resulting from mi-, ri- and mvi-mfa were compared with the mfa configurations of the corresponding original population. performance of the methods was assessed by considering the rv coefficient with respect to the true configuration. in the liver toxicity study, the performances of the methods were compared in two different missingness settings. first the number of missing rows in each stratum of the transcriptomic table was chosen to be low, medium or high. secondly, seven scenarios were created by inserting missing rows in the original dataset which included both transcriptomic and clinical tables. the overall results showed that mi-mfa clearly outperformed the ri-mfa and mvi-mfa approaches in nearly all settings. in the nci- study, we illustrated the performance of our method on complex patterns of missingness where substantial numbers of rows were missing from both tables of the nci- dataset. as previously, this study showed that mi-mfa clearly performed better than mvi-mfa. we also demonstrated that ri-mfa performed better than mvi-mfa. the differences between mi-mfa and ri-mfa were small, but on average slightly in favor of our method. as the purpose of this study was also to illustrate the uncertainty of mi-mfa configurations induced by missing rows, we demonstrated how the areas of the confidence ellipses and convex hulls got larger as the number of missing rows increased.  we propose here a new method, mi-mfa, an extension of mfa, designed to deal with multiple tables with missing row values. mi-mfa is a useful and attractive method to estimate the coordinates of individuals for mfa configurations despite the missing rows. the study cases showed that the other proposed methods either encountered serious problems or were unable to adequately assess the accuracy due to missing data. the configurations obtained with our method were closer to the true configuration even when a significant number of individuals were missing, and thus provided better results. moreover, the uncertainty due to the missing rows could be visualized on the compromise configuration. the software for our mi-mfa method is available in an easy-to-use code for the r statistical environment. additional files additional file supplementary figures. figures s1–s additional file calculation of the total number of possible imputations in mi-mfa. additional file r code implementing the mi-mfa method. r is required. abbreviations fcsfully conditional specification jmjoint modeling marmissing at random mcarmissing completely at random mfamultiple factor analysis mimultiple imputation micemultivariate imputation by chained equations mi-mfamultiple imputation multiple factor analysis mvi-mfamean variable imputation multiple factor analysis pcaprincipal component analysis ri-mfaregularized iterative multiple factor analysis statisstructuration des tableaux à trois indices de la statistique 
incomplete artificial datasets with different patterns of missingness were created from these data. we believe that mi-mfa provides a useful and attractive method for estimating the coordinates of individuals on the first mfa components despite missing rows. the mi-mfa results were compared with two other approaches i.e., regularized iterative mfa and mean variable imputation. this method takes into account the uncertainty of mi-mfa configurations induced by the missing rows, thereby allowing the reliability of the results to be evaluated. for each configuration resulting from these three strategies, the suitability of the solution was determined against the true mfa configuration obtained from the original data and a comprehensive graphical comparison showing how the mi-, ri- or mvi-mfa configurations diverge from the true configuration was produced. missing row values are challenging to deal with because most statistical methods cannot be directly applied to incomplete datasets. mfa is then applied to each completed dataset to produce m different configurations. finally, the m configurations are combined to yield a single consensus solution. a free and easy-to-use code was proposed to implement the mi-mfa method in the r statistical environment. in omics data integration studies, it is common, for a variety of reasons, for some individuals to not be present in all data tables. mi-mfa configurations were close to the true configuration even when many individuals were missing in several data tables. in this study, we focus on multiple factor analysis as a tool to compare and integrate multiple layers of information. mi involves filling the missing rows with plausible values, resulting in m completed datasets. two approaches i.e., confidence ellipses and convex hulls, to visualize and assess the uncertainty due to missing values were also described. we showed how the areas of ellipses and convex hulls increased with the number of missing individuals.
 in omics data integration studies, it is common, for a variety of reasons, for some individuals to not be present in all data tables. missing row values are challenging to deal with because most statistical methods cannot be directly applied to incomplete datasets. to overcome this issue, we propose a multiple imputation approach in a multivariate framework. in this study, we focus on multiple factor analysis as a tool to compare and integrate multiple layers of information. mi involves filling the missing rows with plausible values, resulting in m completed datasets. mfa is then applied to each completed dataset to produce m different configurations. finally, the m configurations are combined to yield a single consensus solution.  we assessed the performance of our method, named mi-mfa, on two real omics datasets. incomplete artificial datasets with different patterns of missingness were created from these data. the mi-mfa results were compared with two other approaches i.e., regularized iterative mfa and mean variable imputation. for each configuration resulting from these three strategies, the suitability of the solution was determined against the true mfa configuration obtained from the original data and a comprehensive graphical comparison showing how the mi-, ri- or mvi-mfa configurations diverge from the true configuration was produced. two approaches i.e., confidence ellipses and convex hulls, to visualize and assess the uncertainty due to missing values were also described. we showed how the areas of ellipses and convex hulls increased with the number of missing individuals. a free and easy-to-use code was proposed to implement the mi-mfa method in the r statistical environment.  we believe that mi-mfa provides a useful and attractive method for estimating the coordinates of individuals on the first mfa components despite missing rows. mi-mfa configurations were close to the true configuration even when many individuals were missing in several data tables. this method takes into account the uncertainty of mi-mfa configurations induced by the missing rows, thereby allowing the reliability of the results to be evaluated. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords multiple omics data integrationmultivariate factor analysismissing individualsmultiple imputationhot-deck imputationissue-copyright-statement© the author 2016 