Do We Really Need Imputation in AutoML Predictive Modeling?

The concept of a missing mechanism [62] formalizes the generation process of missing data. In this respect, the BI are modeled as random variables and assigned a distribution. There are three types of underlying mechanisms that generate missing data, namely, missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR). For formal definitions of these mechanisms, readers are referred to Reference [40]. Intuitively, MCAR implies that the probability of a value missing is independent of the actual value, the other observed quantities, and any latent variables. MAR implies that the missingness only depends on the observed data (so it can be predicted). MNAR refers to the case that the missing values are related to both the observed and unobserved variables, including the missing value itself. When missingness is MNAR, it is in principle and in general not possible to impute the missing values in a way that follows the unknown underlying data distribution.²

An illustrative example is given in Figure 1, which is adapted from Reference [46]. The missingness mechanisms can be described using a causal graph. Let us assume A and B are observed random variables and O a latent variable. Each variable is depicted as a node of the graph. Assume that A and B have a direct connection to O, which is the variable (node) of interest. The \(R_o\) node is a mask variable that denotes the missingness inserted into O, which causes \(O^*\). \(O^*\) is a surrogate of O but with missing values inserted in the positions specified by \(R_o\). As seen in Figure 1(a), MCAR missing values do not depend on any of the variables A, B, 0. In contrast, missingness depends on B for the MAR mechanism, and in O itself for the MNAR data, as seen in Figure 1(b) and (c).

There are numerous imputation algorithms and approaches in the literature, and we do not attempt a full review. Readers are encouraged to explore comprehensive surveys available in the field for a more in-depth understanding [2, 4, 12, 13, 24]. Imputation approaches can be partitioned into various distinct families/groups of methods. A taxonomy is attempted in Figure 2. First, imputed values can be decided based on only the feature with the missing value (Univariate imputation) or several features (multivariate imputation). The former methods include Mean/Median/Max imputation for continuous data and Mode imputation for categorical data. Multivariate methods can be partitioned into Iterative and Distance-Based, also known as Hot-Deck methods [6]. Distance-based methods employ a distance or a similarity metric for samples to find neighbors or cluster them. A commonly used algorithm in this category is the K-nearest neighbors imputation (KNNi) [71], which imputes values based on the neighbors of the sample with missing values. K-means-based methods cluster the samples before imputation [37].

Iterative methods, start with a simple initial guess (e.g., using MM imputation) and, in each iteration, try to improve the imputed values. We further split iterative methods into Discriminative and Generative. Discriminative methods, build a predictive model per feature with missing values, given the other features in the dataset. This model is used to predict the missing values of the corresponding feature, in each iteration. The Discriminative family can either utilize a (generalized) linear model or a non-linear model. Linear discriminative methods include Multivariate Imputation by Chained Equations (MICE) [76]. Non-linear discriminative methods include the MissForest algorithm [66] employing Random Forests, and Datawig [9] that can impute continuous, categorical, and text data by employing different loss functions according to the missing features’ datatype.

Generative methods try to model the joint distribution of the data and use the generative model to impute values. They can be split into two categories, methods that employ matrix factorization and methods that use neural networks. The matrix-factorization family includes low-rank matrix decomposition methods: First, missing values are imputed with an initial guess, and the matrix is decomposed (factorized) and used to predict the missing values. Imputation is improved in each cycle via expectation-maximization steps. Examples of this family include the PPCA [70], SVDImpute [71], bPCA [10], and SoftImpute [44]. Such algorithms scale better w.r.t. to the number of features than MICE or MissForest that train a different model for each feature with missing values in each iteration. Recently, neural networks have also been tried as generative models. These algorithms are essentially non-linear alternatives to matrix factorization. These methods start with an initial guess and then train a neural network that learns the joint distribution. This family includes methods based on AutoEncoders (AE), such as DAE) [21, 35, 41] and Variational Autoencoder (VAE) [17, 45, 61]. Also, it includes generative adversarial networks (GAIN) [64, 83]). Finally, HoloClean, a data cleaning tool, implements an attention-based neural network for imputation, named Aimnet [82]. A detailed comparison of imputation methods is detailed in Section 8. In the next subsection, we will explain the rationale for our choice to include in our empirical study a subset of the aforementioned imputation methods.

In this section, we present the main characteristics of the imputation methods given in Table 1 that we included in our testbed. In the analysis of their computational complexity, n denotes the number of samples, m the number of features, \(\#comp\) the number of principal components, \(\#sing\) the number of singular values, and \(\#trees\) for the number of trees.

MM imputation is selected as a baseline and one of the most commonly used methods. MissForest is selected as a representative of a multivariate iterative imputation method over MICE, based on the results on imputation accuracy presented in Reference [79]. Encoding missingness information using BI is also experimentally evaluated, as it performed better than other methods in Reference [43]. Distance-based methods are excluded for various reasons. First, they need to memorize the full dataset to produce imputation as they do not learn a model. k-means and KNN imputation were not included in our testbed as according to other empirical studies are outperformed by MissForest [30, 42, 66].

As representatives of matrix completion-based methods, we chose PPCA that have shown the best performance according to previous empirical studies [25, 31, 58]. SoftImpute was also selected based on the experimental results presented in Reference [85]. GAIN and DAE were also included as representatives of neural network–based methods as they excel in several studies [11, 54]. The former is based on using a generative adversarial network to learn the probability distribution to impute and the latter on autoencoder. VAE and Aimnet were not included based on the inferior or comparable results to GAIN and MM, respectively [11, 38]. Finally, Datawig [9] is not included as, according to the results reported in Reference [9], it is outperformed by MissForest for both continuous and categorical data while it comes with a high cost to fit one neural network per feature with missing values.

AutoML platforms employ imputation methods, as well as modeling algorithms that directly treat missing values as a separate category. The current versions of JADBio (version v1.4) and AutoSklearn [16] employ MM by default, while DataRobot³ may also include BI variables. AutoSklearn allows the user to specify additional imputation methods to optimize over as part of the pipeline. TPOT employs median imputation for all missing features [36]. Auger.AI⁴ does iterative regression or mean imputation for numerical features depending on the dataset size and creates a new category for the categorical features. BigML⁵ by default does not impute the missing value; the missing values are handled internally by their predictive models, which are based on trees only. Autoprognosis optimizes the ML pipeline over a variety of missing data imputation algorithms. Specifically, it employs MICE, MissForest, Bootstrapped Expectation-Maximization imputation, Soft-Impute, and MM [3]. DriverlessAI by H2O creates a new value to express missingness when the XGBoost, LightGBM, and RuleFit algorithms are used. For generalized linear models, it performs MM imputation, while for tensorflow models missing values are treated as outliers [22]. GAMA [20] does not impute missing values by default. Autogluon [14] uses median imputation for continuous features and introduces a new “Unknown” category for categorical features.

Incomplete Real Datasets. There are currently 364 datasets with missing values in the OpenML repository [78], we restricted our selection to binary classification datasets. We selected 25 binary classification datasets in an effort to cover a range of various dataset characteristics. The datasets contain both continuous and discrete features. The number of features ranges from 7 to 69, the sample size ranges from 155 to 31,406, the prevalence of the minority class ranges from 0.06 to 0.48, the number of features with at least 1% missing values ranges from 1 to 32, and, finally, the percentage of missing values ranges from 1.11% to 71.64%. Table 6 presents the characteristics of the datasets, along with their OpenML id.

Complete Datasets: We selected 10 complete datasets from OpenML, where we introduce and simulate missingness. The number of features ranges from 9 to 135, the sample size ranges from 101 to 5,473, and the prevalence of the minority class ranges from 10% to 49%. Table 7 contains these values for each dataset, along with their OpenML id.

We note that the evaluation concerns only binary classification. The main metric of predictive performance is the Area Under the ROC curve (AUC). To save space and make interpretation easier, we report classification accuracy and F1-score results in the Appendices C. The datasets are split to 50% training and 50% hold-out test set used only for performance evaluation. Our experiments were conducted only once, due to the computational complexity of the experimental procedure (see Section 3.6). We applied statistical tests to compensate for the lack of repeated experiments. This allows reliable conclusions to be drawn from the experimental results.

To experiment with different imputation algorithms when CASH optimization is taking place, we employed the JADBio AutoML platform [73]. JADBio is a commercial product (a version of JADBio with basic functionality is freely available) but was offered to us for research purposes. JADBio includes feature selection as part of the ML pipeline and, thus, it can be used to study the effect of feature selection on imputation.

A quick description of JADBio’s architecture now follows. For each dataset to analyze, an internal knowledge base system, called Algorithm and Hyper-Parameter Space selection (AHPS) in Reference [73], selects the feature construction, preprocessing, feature selection, and modeling algorithms to try, along with a set of values for their hyper-parameters. The AHPS also selects the configuration evaluation protocol, e.g., 10-fold cross-validation, repeated cross-validation, or hold-out to estimate the performance of each configuration and select the winning one. The knowledge in AHPS is engineered by experienced analysts but also induced by meta-level learning algorithms.

The choices of the AHPS are based on the meta-features of the dataset (e.g., sample size, number of features), as well as the user preferences. For example, an algorithm that does not scale to the number of samples in the current dataset, will not be selected by AHPS. The choice of the evaluation protocol also depends on the meta-features: For a typical-sized dataset, JADBio may run a 10-fold cross-validation, for a large balanced dataset a hold-out, while for a small sample or an imbalanced dataset, it may run a repeated cross-validation protocol.

Subsequently, JADBio executes all configurations effectively performing a grid search for CASH optimization. However, JADBio includes pruning heuristics that may drop a configuration in the early folds of cross-validation if it is not deemed promising, departing from a pure grid search strategy [74]. Once configurations execute, the final model is built on all available data using the winning configuration.

The final performance of the model producing with the winning configuration is the cross-validated AUC adjusted for the bias incurred due to multiple tries (called “winner’s curse” in statistics). This adjustment is conceptually equivalent to adjusting p-values in multiple hypotheses testing. JADBio uses the BBC-CV algorithm for the performance estimate adjustment [74]. In Reference [73], experiments on 360 omics datasets of small sample size show that this estimation protocol returns slightly conservative out-of-sample AUC performances of the returned model. Nevertheless, for the purposes of this article, JADBio’s performance estimation was not used; instead, the performances on the 50% held-out set are reported.

Regarding the settings of JADBio employed in this set of experiments, we note the following. One of the user preferences indirectly controls the execution time and the number of configurations to try and has the settings Preliminary, Typical, Extensive, with Extensive trying more configurations and performing a more thorough optimization. All subsequent experiments were run using the Preliminary setting to make the computational requirements manageable. The number of configurations may vary between datasets depending on their meta-features, but in our experiments, it ranges from 900 to over 1,000. The training protocol of JADBio depends on the sample size, the class imbalance, and other factors. For typical-size datasets, JADBio uses a repeated 10-fold cross-validation with #repeats from 1 to 20. A heuristic procedure stops repetitions of cross-validation if no progress is detected. Overall, JADBio uses estimation protocols that execute each configuration between 10 to 200 times per dataset to choose the winning configuration and produce a model.

JADBio optimizes over the following set of algorithms. For feature selection, JADBio uses the Lasso [69] and a variant of the SES algorithm [75] with an upper bound on the number of conditional independence tests to perform. For classification, it optimizes over Ridge Logistic Regression, Decision Tree, Random Forests, and Support Vector Machines with polynomial, linear, and radial basis kernels.

To evaluate imputation algorithms, we embedded them into the JADBio configurations as the second step, after the standardization of continuous features and before feature selection, using the API provided. It is important to note that configurations are cross-validated as an atom, and hence, learning to impute is based only on the training data. This is necessary to avoid overestimating the performances of configurations and correspondingly, the imputation methods. Each imputation method returns an imputation model that is used to impute the test data before modeling is applied. It is worth noting that even if the feature selection step selects a small subset S of features when some values of S are missing in the test set, the imputation model may impute them based on other features. Hence, even if the predictive model requires just the features in S, the predictive pipeline may require more features. Specifically, all multivariate algorithms selected in the article require all features to impute. Hence, the predictive pipeline always requires all features when these algorithms are employed, even with feature selection.

We used the JadBio version 1.4.0 for our experiments. MM and BI methods were already implemented by the developing team of the tool used. For PPCA and SoftImpute, we relied on third-party implementations in R from the PCA methods 1.64.0 [65] and ‘softImpute’ package version 1.4.1 respectively. We implemented MissForest in python 3.8.4 using the iterativeImputer and RandomForest models from sci-kit learn 1.0.1 [53]. Pytorch 1.7.1 version [52] was utilized for the implementation of GAIN and DAE. We adapted the DAE implementation found at https://github.com/Harry24k/MIDA-pytorch to closely follow the description of DAE by the original authors in Reference [21]. We employed the GAIN from https://github.com/dhanajitb/GAIN-Pytorch.

The predictive performance experiments of the article were conducted on a fedora-powered VM using 8-core AMD Threadripper 3970x at 3.7 GHz with 12 GB RAM. The neural networks were trained using CPUs. The execution time results reported were measured on an eight-core AMD Ryzen-3600x at 4.6 GHz with 16 GB RAM and Windows 11 OS.

During the experiments, more than 41 days of CPU time have been spent training more than 80,000 configurations to conduct the experiments mentioned in the article.

The code is available on the Github repository: https://github.com/mensxmachina/Imputation_in_AutoML. The code in the repository consists of scripts for the plots, the datasets, the meta-level analysis as well as the basic implementation of each imputation algorithm.

In the experiments, 24 hyper-parameter (hp) value sets were tried for the imputation algorithms: MM (1 hp set), MissForest (2 hp sets), SoftImpute (3 hp sets), PPCA (3 hp sets), DAE (9 hp sets), and GAIN (6 hp sets). The values tried for each hyper-parameter are shown in Table 2. These choices were based on the algorithm’s authors’ defaults and suggestions. These 24 hp sets were coupled with all other choices of JADBio multiplying by 24 the number of configurations normally tried. In subsequent experiments, each of these 24 hp sets is run on the original dataset, as well as the dataset with the inclusion of the BI features, leading to 48 different combinations. MM has no parameters and therefore does not need tuning. For MissForest, we train RF models with 250 trees, which offers higher imputation accuracy according to Reference [66]. However, we restrict the maximum depth of the tree and maximum leaf nodes, because the trained model had storing memory issues (see Section 2.3.2). SoftImpute and PPCA require selecting the number of principal components to use as a hyper-parameter. The majority of papers in the literature fails to report the tuning of the aforementioned methods that led us to develop the following heuristic: We select as many components required to explain \(x\%\) of the data variance. The values of x are shown in Table 2 as the values of “variance-explained.” The default hyper-parameters are used for DAE with the exception of the dropout layer and the hidden layers’ dimensions. The range of the dropout layer is based on Reference [63], while the theta value is tuned within a neighborhood of the author’s suggested default value. In the current implementation, we have three hidden layers for the encoder and the decoder. For each successive layer in the encoder, \(\theta\) hidden layer nodes are added and hyperbolic tangent is used as the activation function, as it produces better results for small and medium-sized datasets [21]. The model is trained using Stochastic Gradient Descent with an adaptive learning rate with a time decay factor of 0.99 and Nesterov’s accelerated gradient. GAIN architecture consists of three hidden layers for the discriminator and the generator while using Rectified Linear Unit as the activation function. For GAIN, we tune two hyper-parameters; alpha and hint rate. These hyper-parameters are considered the most important for GAIN. Alpha balances the loss between the discriminator and the generator, while the hint rate is responsible for the training of the discriminator. Both DAE and GAIN are trained at the specific epochs as suggested by authors and use the sigmoid activation function for the output layer.

Under MCAR, missing values are missing with a given probability (percentage) independently of any other factors such as the value itself or the values of other features. To simulate missing values at a realistic missingness percentage we sampled 64 real-world datasets from the OpenML repository with varying characteristics (see Section B.1). We then computed the 25%, 50%, and 75% quantiles of missingness percentages. Features with less than 1% of missing values were excluded from the calculation, as they probably point to features that missing values from typos or non-systematic sources. The quantile values turn out to be about 10%, 25%, and 50% of missingness. The quantile values are then used to vary the missingness percentages in both MCAR and MAR simulation experiments. We then introduced missing values with the given percentages at the 10 complete datasets described in Section 3.1. Even though it is trivial to introduce MCAR missing values, for consistency reasons, we employed the code available in Reference [48], which is also used for the MAR simulations below. To simulate the MCAR mechanism the software discards values uniformly at random from the dataset at the specified missingness percentage.

Under MAR, missing values are missing with a probability (percentage) that depends (is conditional) on other observed features, i.e., \(P(I_j = 1|F_{k_1}, \ldots , F_{k_m})\). To realistically simulate data under MAR, one needs to decide (a) the number of features upon which the probability depends, (b) the functional form of the conditional probability function, and (c) the set \(\lbrace F_{k_1}, \ldots , F_{k_m}\rbrace\). To answer (a) we needed a realistic estimate of the number m of features in the conditional probability. To that end, in the corpus of the 10 real-world binary datasets of Table 7, we randomly selected one feature with missing values as the target feature and then performed predictive modeling using JADBio including feature selection.⁶ These experimental results suggest that, on average, a feature with missing values is dependent on 12 features, so we set \(m=12\). Subsequently, for each \(F_j\) we randomly selected a set of m other features with uniform probability. Finally, for the functional form of P, we used a logistic regression model: \(P(I_j = 1|F_{k_1}=f_1, \ldots , F_{k_m}=f_m) = \frac{1}{1+e^{\langle -w, f\rangle }}\), where w is a set of randomly chosen coefficients from a normal Gaussian distribution, f is the vector of values of the features \(F_{k_l}\), and \(\langle \cdot , \cdot \rangle\) denotes the inner product. For the simulation, the software [48] was also used. The software allows the simulation of MAR missing data, as described above, with prespecified missingess percentages. The same percentages as in the MCAR case were used.

First, we partition results achieved when optimizing over any single imputation algorithm. Specifically, for each imputation algorithm, on a given dataset, the best AUC was selected over all configurations that include the specific algorithm. We will refer to this best AUC simply as the AUC of a given imputation algorithm, in all subsequently reported results. Figure 3(a) shows the difference in AUC performance when Binary Indicators are used versus when excluded. As we can see, MM and GAIN have the largest average increase by 0.0074 AUC and 0.0056 AUC, respectively. MF when extended with BI shows an average AUC increase of 0.0046, while PPCA shows an increase of 0.0038 AUC. The lowest average improvement is achieved by SOFT, which improves by 0.00022 AUC. Contrary to the above observations, DAE is the only method that does not benefit from the addition of BIs with a negligible decrease of 0.0004 AUC when BI’s are included. Figure 3(b) offers a complementary view. It illustrates the count of datasets per imputation method where the inclusion of BI is beneficial to the downstream performance. We observe that for every imputation method, including BI is beneficial in most instances. SOFT exhibits improvement across 19 of 25 datasets. GAIN, MM, and PPCA in 17 datasets. Finally, including BI, leads to enhancements in DAE and MF across 16 and 15 datasets, respectively.

To determine the statistical significance of the results, we performed a paired matched t-test for each algorithm with the null hypothesis H0 being that the BI+base has worse performance than the base method. The resulting p-values were converted to q-values with the Benjamini/Hochberg [7] method, to control for multiple testing. Table 3 shows the results. Using a q-value threshold of 0.25 there are four statistically significant results, resulting in accepting the alternative hypotheses that MF, GAIN, PPCA, and MM improve their performance when BIs are present. At the level of \(q=0.25=\frac{1}{4}\) this implies that, in the worse case, we expect one of these four discoveries to be false. While the inclusion of BIs may, in the worse case, double the dimensionality of the dataset, based on the above results, we would recommend their inclusion when the above imputation algorithms are employed.

Figure 4 shows the average ranking achieved by each algorithm using the Autorank tool [26] (lower ranking is better). To avoid clutter, and based on the results of Section 5.1, we only show results when BIs are included. The horizontal black bars in the graph connect tools with non-statistically different ranks, according to a non-parametric Friedman test and post hoc Nemenyi test.

Results show that BI+DAE is the highest ranking algorithm with an average rank of 2.84, followed by BI+MM with a 2.94 average ranking, BI+MF with 3.5, and BI+GAIN with 3.56, although their rank difference is not statistically significant at the 0.05 level. The two lowest-ranked methods are BI+SOFT and BI+PPCA, with 3.74 and 4.42 average rankings, respectively. BI+DAE’s rank is statistically significantly lower compared to BI+PPCA.

We now study the performance effectiveness vs. the computational efficiency tradeoff of the algorithms. In Figure 5(a) we use MM as the baseline. A point (execution run) corresponds to AutoML predictive modeling on a dataset with a given imputation algorithm. This results in \(5 \times 25 = 125\) points. The x-axis shows the effectiveness ratio defined as the ratio of the AUC corresponding to the point divided by the corresponding performance of BI+MM. Similarly, the y-axis shows the efficiency ratio defined as the training time of the point divided by the corresponding time of BI+MM. Hence, points in the first/fourth quadrant (top-left/bottom-right) correspond to runs where BI+MM dominates/is-dominated by other algorithms on the same datasets in both time and AUC. Notice that the scale of the y-axis is logarithmic. Larger points correspond to the mean value of an imputation method over all datasets.

In total BI+MM is inferior in terms of predictive performance in 42 cases (16 of the 25 datasets) and, unsurprisingly, never gets dominated in terms of AUC performance and efficiency at the same time. The computational time of the other algorithms is orders of magnitude slower than BI+MM. However, in 83 of 125 points, BI+MM is both more efficient and effective than the compared method. Only BI+DAE scores on average higher than the AUC score. All the other imputation methods are on average slower to train and worse in terms of predictive performance.

Figure 5(b) shows the same exact results with BI+DAE as the baseline. In contrast to BI+MM above, BI+DAE dominates the other imputation methods on predictive performance and training time, in only 35 of 125 combinations. In 41 of 125 cases, it provides better predictive performance but at a higher computation cost. In 15 points, BI+DAE is faster but has lower predictive performance than the compared imputation method. Finally, 34 times it is dominated in both metrics. In conclusion, when if a single imputation algorithm is to be used, BI+MM arguably provides the best tradeoff between computational time and predictive performance.

In this section, we examine the results from a different perspective, trying to answer the question: What is the minimal-size subset of algorithms to try to achieve close-to-maximum AUC performance? To answer this question, we have implemented a simple greedy algorithm, where we assume the analyst starts with the subset \(\lbrace\)BI+MM\(\rbrace\) as an efficient baseline and adds algorithms to consider. In each iteration, the algorithm that leads to the largest AUC improvement of the subset when added is selected for inclusion. The maximum AUC performance is the sum of the maximum AUC for each dataset when including all imputation methods in the optimization pipeline, averaged across all datasets.

The results are shown in Figure 6 and quantitatively in Table 11. The x-axis shows the imputation algorithms in order of addition to the subset. For each algorithm, several hyper-parameter combinations are tried and combined with all other feature selection and modeling choices by AutoML. Hence the total number of configurations tried is multiplied by this factor. At each tick, the multiplication factor for the whole set is depicted in the parenthesis next to the name of the algorithm added to the set in that step. For example, BI+MM has no hyper-parameters (\(1\times\)), while BI+DAE has 9, so the multiplicative factor of the set \(\lbrace BI+MM, BI+DAE\rbrace\) is \(10\times\). The y-axis is the average (over all datasets) relative AUC achieved when performance is optimized over all algorithms and their hyper-parameters in the corresponding subset.

BI+MM, by itself, accounts for 98.69% of the maximum AUC. When BI+DAE is added to the mix, relative performance reaches 99.69%. The next best algorithm to add is BI+MF; 100% of AUC is reached when invoking all imputation algorithms. In summary, the addition of BI+MF, BI+PPCA, BI+GAIN, and BI+SOFT provide only marginal gains to the set \(\lbrace\)BI+MM, BI+DAE\(\rbrace\).

Feature selection algorithms try to reduce the number of features that enter the model without sacrificing predictive performance. Feature selection is often the primary task in analysis, while the predictive model may be just a side-benefit. For example, a medical doctor may be more interested in the quantities that determine the risk of disease and may reveal new medical knowledge, rather than the risk model itself. Feature selection leads to more interpretable models that provide intuition into the domain. In fact, the solution to the feature selection problem is directly linked to the causal model that underlies data generation [72]. In other circumstances, it is important to reduce the cost of measuring the features to provide predictions. The cost may be measured in monetary units, the computational cost to compute the features or risk to a patient from medical procedures that measure these features.

Figure 7 shows the impact of feature selection for each imputation algorithm on the real dataset. The drop in AUC performance when feature selection is enforced vs. not enforced (i.e., optimizing over all configurations) in the final configuration is shown. For each algorithm is about two to three AUC points. In other extensive experiments with hundreds of complete (no missing values), small-sample, high-dimensional omics datasets, JADBio has been shown to reduce the number of features by a factor of 4,000 without a noticeable drop in AUC performance [73]. The results provide evidence that feature selection may be more challenging in the presence of missing values.

In any case, the problem of including both imputation and a feature selection step in the ML pipeline is that imputation invalidates feature selection, in some sense. Let us explain this statement with an example. Let us assume the pipeline that produces the final model consists of MF imputation, Lasso feature selection, and RF predictive modeling. Let us assume that Lasso selects the features \(\lbrace A, B, C\rbrace\). If any of these values (say the value of A) is missing on a new sample, then the MF imputation model will impute them using a Random Forest for A using some other subsets of features. If any of those are also missing, then MF will invoke its Random Forests for each value that is missing, and so on, recursively. Hence, if there are missing values on the test samples, one may need to measure an arbitrarily large feature subset, not just the selected features. The storage required to apply the ML pipeline includes both the RF as well as the MF model, which in turn includes a Random Forest for every feature that may need imputation.

Figure 8 presents the AUC performance results (see Figure 16(b) for MCAR results). The AUC performance denoted is the absolute difference in performance at the specified missingness percentage minus the AUC performance of the complete dataset. First, we note that the figure illustrates that as the missingness percentage increases the average predictive performance for every imputation method used decreases, as expected. As we can see, increasing the missingness from 25% to 50% leads to a sizable performance drop for all imputation methods. Specifically, methods based on linear dimensionality reduction, namely BI+PPCA and BI+SOFT are the most affected by this increase in missing values.

Figures 8 and 16(b) illustrate that in both MAR and MCAR data, respectively, MissForest combined with Binary Indicators is, on average, the best-performing method. Additionally, we note that PPCA and SOFT are the two worst imputation methods, especially as the missingness percentage increases. Table 12(a) in Appendix C.2 contains the quantitative results in detail and a detailed discussion on the ranking of the algorithms.

We now identify the minimal-size algorithm subset with close-to-optimal performance for simulated missing data. We use again the simple greedy algorithm introduced in Section 5.4 and apply it to the MCAR and MAR simulated data results. The results for MAR are in Figure 9, which is similar to Figure 6. The quantitative results are shown in Table 13(b). As shown in the figure, the \(\lbrace BI+MM, BI+MF\rbrace\) subset can score over 99% of the total max AUC for MAR data and would be the suggested set of algorithms to run in such problems. The results for MCAR are in Appendix C.2.6. They are qualitatively similar. The results for simulated missing values are somewhat different than the ones in the real datasets, namely BI+MF scores better than BI+DAE, which is now placed in third place. Possible reasons why are discussed in Section 9.

Figures 10 and 16(a), show the effectiveness vs. efficiency tradeoff of the algorithms. The aformentioned figures are similar to Figure 5(a) above for the real datasets. We repeat the explanation of the figure: The reference (baseline) algorithm is BI+MM. The x-axis shows the effectiveness ratio defined as the ratio of the AUC corresponding to the point divided by the corresponding performance of BI+MM. Similarly, the y-axis shows the efficiency ratio defined as the training time of the point divided by the corresponding time of BI+MM. Hence, points in the first/fourth quadrant (top-left/bottom-right) correspond to runs where BI+MM dominates/is-dominated-by other algorithms on the same datasets in both time and AUC. Notice that the scale of the y-axis is logarithmic. Larger points correspond to the mean value of an imputation method over all datasets. There are five imputation methods to compare against MM for 10 datasets over 3 percentages of missing values. This will naturally result in 150 points. However, MissForest did not run in three of the datasets (image dataset variations) due to its dimensionality; see Section 2.3.2 for details. The resulting plot will consist of 147 points.

For MAR data, BI+MM is never dominated in both metrics, as it is by far the most efficient method. In 104 of 147 cases, it dominates the opposing imputation methods in terms of both effectiveness and efficiency. However, 43 times it is dominated in efficiency. Only BI+MF has on average better predictive performance than BI+MM. However, it is 23,000 times slower to train on average. All the other imputation methods are worse than BI+MM on average while also taking more time to train. The results for MCAR data are qualitatively similar (see Figure 16(a) and discussion in Appendix C.2).

In total, BI+MM is again found to provide the best, arguably, tradeoff between efficiency and effectiveness. The results in the simulated data are further validating the results in the real-world data, verifying that BI+MM is indeed a decent imputation method all around. BI+MM, on average, is on par with more sophisticated methods such as BI+MF, BI+GAIN, and BI+DAE, while being thousands of times faster to train.

Compared to the related work, we contribute in various ways. Our work can be directly compared to two other works that are conducted in an AutoML tool [19, 49]. Compared to the mentioned works, we include more datasets, more missingness mechanisms, neural network methods, and binary indicators in the experimental setup. For the first time, deep learning methods are compared to simple imputation methods in an AutoML predictive setting. One of the deep learning methods (BI+DAE) has the best average performance on real-world data with native missing values. Additionally, for the first time, the effect of the imputation methods on predictive performance is measured on datasets with generated missing values and native missing values. Until now, comparisons were conducted on only one of the two settings, specifically half of the papers use real-world datasets with missing values in them, while the other half use complete datasets with generated missing values. Contrary to the majority of the literature, we tune both imputation and predictive methods to fairly evaluate them. Only two of the eight related mention tuning both imputation and predictive modeling methods [38, 49]. We also conducted experiments on more datasets compared to the majority of the literature, while unlike [30], our simulation setting is applied to all features in the datasets and not only one. Also, only one of the eight aforementioned works includes feature selection as part of the ML pipeline. Finally, meta-learning, for the first time, is used to identify useful data characteristics that could give insights into the choice of a simple vs. a sophisticated imputation method. In general, as shown in Table 5, our testbed is the most complete overall in terms of dataset selection, missingness selection, imputation method selection, and pipeline steps. This allows us to fairly evaluate imputation methods in a state-of-the-art AutoML environment.