Missing Data Imputation with Adversarially-trained Graph Convolutional Networks

Introduction

While machine learning and deep learning have achieved tremendous results over the last years (Goodfellow et al. 2016), with new breaktroughs arising constantly (e.g., in drug discovery (Chen et al. 2018)), the vast majority of supervised learning methods still require datasets with complete information. At the same time, many real-world problems require dealing with incomplete data, such as in the biomedical or insurance sectors (Van Buuren 2018). For this reason, flexible missing data imputation (MDI) methods are a fundamental component for widespread adoption of machine learning. An MDI algorithm takes a dataset with missing values in some of its input vectors, and replaces these values with some appropriately predicted ones in order to obtain a full dataset. In the literature, this problem is also called missing value imputation (MVI) (Lin and Tsai 2019). In this paper, we use the terms data and value interchangeably based on context. In particular, there is the need for powerful multivariate imputation methods able to work in a variety of data generation regimes (Yoon et al. 2018).

It has been recognized for a while that data imputation can be framed under a predictive framework, and classical machine learning methods (e.g., for regression and classification) might be adapted for this task (Bertsimas et al. 2017). However, some care must be taken when adapting them, for two fundamental reasons. Firstly, different inputs in general have different missing components, while most machine learning models assume full input vectors. Secondly, MDI might be a simple preprocessing step for downstream learning tasks, and in this case performance in terms of reconstruction might not be a perfect proxy for classification/regression accuracy later on.

In general, the resulting predictive approaches to MDI can be classified depending on whether they try to build a global model for data imputation, or whether they use similar data points to infer the missing components. Algorithms in the latter class include using simple statistics computed from the entire dataset (e.g., medians), or more advanced k-NN strategies (Lakshminarayan et al. 1996). In the former case, instead, we have simple linear models (Lakshminarayan et al. 1996), support vector machines (Wang et al. 2006) or, more recently, deep neural architectures (Yoon et al. 2018). These are surveyed more in-depth in Section 2.1.

We argue that a more powerful technique for MDI should exploit both ideas, i.e., use similar data points for each imputation and global models built from the overall dataset. In fact, recently a large class of neural network techniques have emerged that are able to model and exploit this kind of structured information (in the form of relationships between examples), by working in the domain of graphs (Bronstein et al. 2017; Battaglia et al. 2018). These models have been applied successfully to a wide range of problems, among which recommender systems (Ying et al. 2018), quantum chemistry (Gilmer et al. 2017), entity extraction from relational data (Schlichtkrull et al. 2018), semi-supervised learning (Kipf and Welling 2017), and many others. However, to the best of our knowledge these techniques have never been applied to MDI. To overcome this, in this paper we define an architecture for MDI based on a specific class of graph neural networks, namely, graph convolutional networks (GCN), and empirically evaluate it on a large set of benchmark datasets (Kipf and Welling 2017).

Our generic framework for MDI is shown later on in Figure 1. We frame the overall problem in terms of a GCN autoencoder, We use the term autoencoder to refer to any architecture that learns to map an input (or a corrupted version in the case of denoising autoencoders) to itself. that learns to reconstruct the overall dataset conditioned on some artificial noise added during the training phase (similar to a classical denoising autoencoder (Vincent et al. 2008)). To build a graph of similarities between points we leverage prior literature on manifold regularization (Belkin et al. 2006), and we describe a simple technique that was found to work well in most situations.

After describing the basic architecture, we also detail three extensions to it that are able to improve either the accuracy or the speed of convergence:

Firstly, we train the autoencoder with a mixture of standard loss functions and an adversarial loss, which was shown to provide significant improvements for denoising autoencoders in the non-graph case (Yoon et al. 2018).

Secondly, we motivate another extension with the inclusion of residual connections from the input to the output layer, similar to residual networks (He et al. 2016).

Finally, we also describe how to include global information on the dataset (e.g., means and medians for all feature columns) using a generic context vector in input to the GCN layers.

We test our overall architecture on a large benchmark of datasets with varying levels of artificially-added noise and three real-world datasets with pre-existing missing values (two biomedical datasets and one time-series dataset). For the former, we show that our proposed GINN method is on par or outperforms several existing state-of-the-art approaches, especially when we consider high levels of injected artificial noise, e.g., up to $50\%$ of missing values in the original dataset. For the latter, we show that our method is robust to the selection of a downstream classifier, with an accuracy comparable to any other combination of an imputation method and a classifier.

Organization of the paper

The rest of the paper is organized as follows. In Section 2 we describe the relation of this paper with state-of-the-art methods for MDI (Section 2.1) and graph neural networks (Section 2.2). The GCN, which is the building block of our method, is described in Section 3. Then, our graph imputation neural network (GINN) framework and all its extensions are described in Section 4. After a large experimental evaluation in Section 5, we provide some concluding remarks in Section 6.

Related work

Algorithms for MDI can be categorized depending on whether they perform univariate or multivariate imputation, and on whether they provide one or multiple imputations for each missing datum (Van Buuren 2018). In addition, different algorithms can make different theoretical assumptions on whether the data is missing completely at random (MCAR) or not. In this paper we consider multivariate imputation, which is standard in the neural network’s literature. In the following we briefly review state-of-the-approaches on this topic, including several algorithms that we will compare to, and discuss their relation with our proposal.

A popular technique for MDI is multiple imputation using chained equations (MICE) (Azur et al. 2011; White et al. 2011; Van Buuren 2018). MICE iteratively imputes each variable in the dataset by keeping the other variables fixed, repeating this for multiple cycles, each time drawing one or more observations from some predictive distribution on that variable. Although MICE has shown very good performance in some settings, especially in the bio-medical sector, the assumptions beyond MICE (especially the MCAR assumption) might result in biased predictions and subsequently lower accuracy (Azur et al. 2011).

In the machine learning community, it was recognized very soon that MDI can be framed as a predictive task, on which variants of standard supervised algorithms can be applied, including k-nearest neighbors (k-NN) (Acuna and Rodriguez 2004), decision trees (Lakshminarayan et al. 1996), support vector techniques (Wang et al. 2006), and several others. However, these techniques have always achieved mixed performance in practice compared to simpler strategies such as mean imputation (Bertsimas et al. 2017). k-NN is limited in making weighted averages of similar feature vectors, while other algorithms are required to build a global model of the dataset to be used for imputation. In this paper we also frame the MDI problem in a predictive context, but our proposed model can leverage both global aspects of the dataset and local similarities between different points.

More recently, there has been a surge of interest in applying deep learning techniques to the problem of MDI. These include multiple imputation with deep denoising autoencoders (MIDA) (Gondara and Wang 2018), combinations of deep networks with probabilistic mixture models (Śmieja et al. 2018), recurrent neural networks (Bengio and Gingras 1996; Che et al. 2018), or generative models including generative adversarial networks (Yoon et al. 2018) and variational autoencoders (Nazabal et al. 2018). Generally speaking, these methods are better at capturing complex correlations in the data (and in the missing data process), thanks to their multiple layers of nonlinear computations, but they still require to build a global model from the dataset, while ignoring potentially important contributions from similar points. The method we propose can be seen as an extension both of the MIDA algorithm and of Yoon et al. 2018, but we focus on a more recent class of NNs, graph NNs, to capture local dependencies. We briefly survey the literature on this topic next.

2 Graph neural networks

Some of the earliest works on extending NNs to the domain of generic graphs were presented in Gori et al. 2005; Scarselli et al. 2009, and later reformulated in Li et al. 2015 in a more recent context. These works were mainly motivated by the analogies between unrolled recurrent neural networks and the diffusion of information across a graph.

Another line of work, upon which we build our proposal, considers instead the extension of convolutional neural networks to graph domains under the general term of geometric deep learning (Defferrard et al. 2016; Kipf and Welling 2017; Bronstein et al. 2017) (and (Micheli 2009) for earlier works on a similar context). This is done by exploiting recent ideas in the field of graph signal processing (Sandryhaila and Moura 2013; Sardellitti et al. 2017) to define a more general convolution operator able to work on irregular data structures. In particular, in this work we use the GCN of Kipf and Welling 2017, that for every layer includes a linear diffusion process across neighbors. Additional interesting lines of research in building GNNs that we briefly mention include earlier works on graph autoencoders (Sperduti 1994), graph attention networks (Veličković et al. 2018), non-local NNs (Wang et al. 2018), graph embeddings (Zhang et al. 2018), and tree/graph echo state networks (Gallicchio and Micheli 2010; Gallicchio and Micheli 2013). An overview of many of these ideas is provided in Battaglia et al. 2018. We explore some of the ideas from Battaglia et al. 2018 in our framework by discussing how to include global information about the dataset in the reconstruction process in Section 4.5.

Finally, our work is related to the field of manifold regularization (Belkin et al. 2005; Belkin et al. 2006), a semi-supervised class of methods that exploits similarity information among patterns to enforce a regularization term on the optimization process. We build upon them for the construction of our similarity graph, a necessary step for exploiting the power of GNNs.

Graph convolutional networks

Because the GCN layer is a fundamental building block of our method, we briefly describe it here before moving on to the proposed framework for MDI. Consider a set of $n$ vertices of a directed graph, whose connectivity is described by a (weighted) adjacency matrix $\mathbf{A}\in\R^{n\times n}$, where $A_{ij}$ is different from $0$ if and only if nodes $i$ and $j$ are connected. Each node $i$ has an associated vector of features $\mathbf{x}_{i}\in\R^{d}$, that we collect row-wise in the matrix $\mathbf{X}\in\R^{n\times d}$. We would like to have a generic neural network component able to process simultaneously the features at every node, but also take into consideration their relations, expressed through the adjacency matrix.

One way to extend the idea of convolutional networks to this domain is the so-called graph Fourier transform (Bruna et al. 2013; Sandryhaila and Moura 2013). Define the Laplacian matrix of the graph as $\mathbf{L}=\mathbf{D}-\mathbf{A}$, where $\mathbf{D}$ is the diagonal degree matrix with $D_{ii}=\sum_{j=1}^{n}A_{ij}$. We can perform the eigendecomposition of this matrix as $\mathbf{L}=\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{T}$, where $\mathbf{U}$ is a matrix collecting column-wise the eigenvectors of $\mathbf{L}$, and $\mathbf{\Lambda}$ is a diagonal matrix with the associated eigenvalues. The equivalent of a classical Fourier transform on a signal can be defined in the graph domain as (Sandryhaila and Moura 2013):

and the inverse transform as $\mathbf{X}=\mathbf{U}\hat{\mathbf{X}}$. Using this, a straightforward way to define a convolutional layer on graphs (Bruna et al. 2013) is to first apply the graph Fourier transform, apply a trainable transformation on the frequency components (associated to the eigenvalues of the Laplacian), and then back-transform using the inverse Fourier transform. While viable, this approach is however costly and impractical in most cases.

Later authors (Defferrard et al. 2016; Kipf and Welling 2017) have noted that by applying a restricted class of filters to the frequency components (polynomials), it is possible to work directly in the graph domain using polynomials of the Laplacian itself. In particular, Kipf and Welling 2017 proposed the GCN with the use of linear filters, resulting in the following canonical layer:

where $\mathbf{\Theta}_{1}$ is a matrix of adaptable coefficients, and $g(\cdot)$ a generic element-wise activation function, such as the ReLU $g(s)=\max\left(0,s\right)$. To avoid some numerical instabilities, it is possible to renormalize the Laplacian to properly bound its eigenvalues, ad done in Kipf and Welling 2017. More in general, one can substitute the Laplacian with any valid graph shift operator (Gama et al. 2019). Note that the right-multiplication by $\mathbf{\Theta}_{1}$ is akin to a classical feedforward layer, while the left-multiplication by $\mathbf{L}$ allows to propagate the information across the immediate neighbors of each node. Multiple layers of this form can then be stacked to obtain a complete graph NN. Importantly, for a generic network with $L$ layers of the form (2), the output of node $i$ will depend on the outputs of its neighbors up to degree $L$.

Proposed framework for missing data imputation

In MDI, we are also given a data matrix $\mathbf{X}$, which has the same size and semantic as in the previous section, but in general no graph information associated with it. Some of the values of $\mathbf{X}$, denoted by a binary mask $\mathbf{M}\in\left\{0,1\right\}^{n\times d}$, are missing and need to be imputed for downstream processing or classification/regression. We assume to have either numerical features, which are properly normalized, or categorical features that are represented with one-hot encoding. For other types of features, e.g., text, a previous embedding step is needed (Pennington et al. 2014). When facing a supervised learning problem, for which there is an additional label (e.g., class) associated to each input $\mathbf{x}_{i}$, we can easily include the training labels in the imputation process by concatenating them to the input vector.

Predictive models for MDI described previously in Section 2.1 build a function $f(\mathbf{x}_{i})$ for imputing missing values of a single example $\mathbf{x}_{i}$, but in general, do not exploit directly the potentially important information contained in points that might be similar to it. Here, we propose to model this constraint explicitly by building $f$ using GCN blocks, as shown schematically in Figure 1. To do this, we first need to build a graph of inter-patterns similarities from $\mathbf{X}$, as described in the next section.

The first fundamental step of our method is the discovery of the graph structure underneath the tabular representation of the data. In the resulting graph, each feature vector in the dataset is encoded as a node of the graph, while the adjacency matrix $\mathbf{A}$ is derived from a similarity matrix $\mathbf{S}$ of the features vectors. As we stated before, constructing a similarity graph from the data is a known problem in the literature, and here we leverage some work from the field of manifold regularization (Belkin et al. 2005; Belkin et al. 2006), adapting it for handling the presence of missing data. A small overview on possible alternatives is provided at the end of this section.

The similarity matrix is computed pairwise for all features vectors using the Euclidean distance, but each time only the non-missing elements of both vectors are used for the computation (Eirola et al. 2013):

where $\odot$ stands for the Hadamard product between vectors, $\mathbf{M}_{i}$ is the $i$th column of the binary matrix $\mathbf{M}$, and $d$ is the Euclidean distance.

In practice, it is common to have sparse similarity matrices (Belkin et al. 2006), most notably for efficiency and computational reasons (see in particular the discussion on computational cost later on), allowing most operations to scale at-most linearly in the number of neighbors. In order to have a sparse graph, with meaningful connections between similar nodes, we apply a pruning step over the similarity matrix $\mathbf{S}$ inspired from the large-scale manifold learning algorithm in Talwalkar et al. 2013. A threshold is applied independently over every row of $\mathbf{S}$, by computing a percentile for each row that will act as a threshold, such that only the connections above this threshold inside every row are kept. To make this process more robust, we iterate it twice. The result is used as adjacency matrix for the GCN in the next section.

We found that the $97.72$nd percentile provides good results over all the datasets used in the benchmark of Section 5, allowing us to discard at each step around $95\%$ of all the possible connections (in accordance with Talwalkar et al. 2013). Figure 2 shows a subset of the graph produced with this procedure starting from the Iris dataset; here we display the relationship between the number of missing features in a node and its degree. It can be seen from the figure how nodes with very few non-zero elements are correlated with a higher number of connections. This is due to the absence of a re-normalization step in the computation of the similarity matrix (Eirola et al. 2013). In practice, we have found this setup to work better than renormalizing all distance measures, possibly because of the increased degree of elements with multiple missing values.

On the construction of the similarity graph: the method described in this section, which is the one we follow in our experiments and in our open-source implementation, was found to provide good empirical performance. Nonetheless, we underline that in the proposed GINN framework, the similarity graph can be built according to any guideline or method available in the literature. For example, classical alternative choices in the semi-supervised literature include binary weights on the edges, heat kernel similarity (Belkin et al. 2006), or selecting a fixed number of neighbors instead of a percentile (Geng et al. 2012). If the inputs contain text, images, or similar data, cosine similarity on custom pre-trained embeddings are also a popular choice (Bui et al. 2018). We leave an analysis of these different alternatives to future work.

2 Autoencoder architecture

Autoencoders are composed by an encoder which maps the inputs to an intermediate representation in a different dimensional space $\mathbf{h}=\text{encode}(\mathbf{x})$, and a decoder that maps $\mathbf{h}\in\R^{m}$ to the original dimensional space $\hat{\mathbf{x}}=\text{decode}(\mathbf{h})$. We use $m>d$ for an overcomplete representation, thus mapping the input into a higher dimensional space with the aim of helping data recovery. Our graph imputer neural network (GINN) will thus be defined as follows:

where $\mathbf{L}$ has been defined in Section 3 (the extension to networks with multiple hidden layers being straightforward).

Note that we cannot trivially train the autoencoder on the missing values, because they are not known in the training stage. To solve this, we adopt a denoising version of the autoencoder (Vincent et al. 2008), in which for every optimization step we add additional masking noise on the input, by the means of an inverted dropout layer applied directly on the input of the network. In particular, at each optimization step, we randomly remove $50\%$ of the remaining inputs. The only exception: whenever training labels are used as inputs for the imputation process, we do not apply dropout on them. In this way, the autoencoder learns to reconstruct any part of the input matrix, similarly to (Gondara and Wang 2018).

We train the whole model end-to-end minimizing the reconstruction error over the non-missing elements of the dataset. The loss function is thus defined as the combination of a mean squared error (MSE) for the numerical variables and the cross-entropy (CE) for the categorical variables:

where MSE always returns $0$ for categorical values and vice versa for CE. $\alpha$ is an additional hyper-parameter that we initialize as the ratio between the number of numerical columns of the dataset and the total number of columns. Alternatively, it can be tuned like the other hyper-parameters of the network, although we have not found any definite improvement in doing so.

The computational cost of our approach is related mostly to (a) the one-time cost of constructing the similarity graph, and (b) the use of a GCN layer instead of a standard feedforward layer as in alternative autoencoder architectures (Gondara and Wang 2018). The cost of point (a) is well-studied in the literature on large-scale similarity search, e.g., (Dong et al. 2011), and many techniques and implementations exist for making it efficient. The cost of the GCN layer is discussed in Kipf and Welling 2017. In particular, using a sparse representation for the adjacency matrix reduces the memory requirement to $\mathcal{O}(|E|)$, where $E$ is the number of edges in the graph, and the cost of computing (2) to $\mathcal{O}(|E|CF)$, where $C$ and $F$ are the number of input and output features respectively. In practice, when using an early-stopping strategy for training, we have found the training time for our GINN algorithm (including point (a)) to be significantly faster than alternative neural approaches and on-par with non-neural competitors such as MICE, e.g., see Fig. 4 in the additional materials for a comparison.

3 Adversarial training of the autoencoder

In order to speed up training, we use an additional adversarial training strategy where a critic, a feedforward network in our case, learns to distinguish between imputed and real data. This is inspired by generative adversarial networks (Goodfellow et al. 2014) and was found to have significant effects in several reconstructions tasks, particularly in the medical domain (Ker et al. 2018). In particular, having an adversarial loss during reconstruction forces the reconstructed vector to lie close to the natural distribution of the original patterns (Shen et al. 2019).

To train jointly autoencoder and critic and have a stable training we used the Wasserstein distance introduced in Arjovsky et al. 2017, which is informally defined as the minimum cost of transporting mass in order to transform a distribution $q$ into a distribution $p$. Using the Kantorovich-Rubinstein duality (Villani 2008) the objective function is obtained as follows:

The original loss in Arjovsky et al. 2017 used weight clipping to force the Lipschitz property. A further step towards training stability, introduced in Gulrajani et al. 2017, is to use a gradient penalty instead of the weight clipping, obtaining the final loss:

since it must fool the critic and minimize the reconstruction error at the same time.

4 Including skip connections in the model

The autoencoder itself generates an approximate reconstruction of the dataset, while the critic loss guides the autoencoder in this process. However, our main scope is the imputation of values not present in the data. For this task we want a greater contribution coming from the most similar nodes. For this reason, we introduce an additional skip layer which consists always in a graph convolution operation but propagating the information across the immediate neighbors of each node without the node itself. This prevents the autoencoder from learning the identity function.

5 Including global statistics from the dataset

Another extension we explore is the possibility of including global information on the dataset during the computation of the autoencoder. The inclusion of a global set of attributes, in the context of graph neural networks, was described in-depth by Battaglia et al. 2018.

As a proof of concept, in our case we set as global attribute vector $\mathbf{g}$ for the graph some statistical information of the dataset, including mean or mode of every attribute. The global component can be taken into account in the last layer, weighting their contribution for the update of each node:

In addition, if the computation of the global information is differentiable, we can compute a loss term with respect to the global attributes of the original dataset:

where $\gamma$ is some additional weighting term.

Experimental evaluation

We divide our experimental evaluation in five subsections. Firstly, following common literature, we evaluate the proposed GINN framework on 20 real-world datasets from the UCI Machine Learning Repository (Dua and Graff 2017) to which we artificially add some desired level of missing values, to evaluate the imputation performances. The characteristics of these 20 datasets are summarized in Table 1. This selection contains categorical, numerical, and mixed datasets, ranging from 150 observation to 30000 and from only 4 attributes to almost 40. Every dataset is divided into training 70% and test 30% sets. Missingness is introduced completely at random on the training set with 4 different levels of noise: 10%, 20%, 30%, and 50%. Our evaluation will focus first on imputation performance as described in Section 5.1, then on the accuracy of post-imputation prediction in Section 5.2. We then perform a comprehensive ablation study of the architecture in Section 5.3, and an evaluation of the performance of the algorithm on new data in Section 5.4.

Secondly, in Section 5.5 we evaluate the performance of the algorithm on three real-world datasets with pre-existing (i.e., non artificially induced) missing values. In this part we also evaluate the computational cost of the method when compared to other state-of-the-art approaches. In one case, missing values are also present in the training labels, in which case we provide an application of our method to a semi-supervised scenario (as described later on).

For all the benchmarks, we use an embedding dimension of the hidden layer of $128$, sufficient for an overcomplete representation for all the datasets involved, and we train the model for a maximum of $10000$ iterations with an early stopping strategy for the reconstruction loss over the known elements. The critic used is a simple 3-layer feed-forward network trained 5 times for each optimization step of the autoencoder. We used the Adam optimizer (Kingma and Ba 2014) for both networks with a learning rate of $1\times 10^{-3}$ and $1\times 10^{-5}$ respectively for autoencoder and critic. When label information is available for the datasets, we consider the training labels as an additional feature of each input vector, but we remove this information when processing new (validation or test) data. All experiments are repeated five times and we collect average performance and standard deviation.

All the code for replicating our experiments and using the GINN algorithm is released as an open-source library on the web. https://github.com/spindro/GINN

This evaluation focuses on the comparison of MAE and RMSE between GINN and 6 other state-of-the-art imputation algorithms: MICE (van Buuren and Groothuis-Oudshoorn 2011), MIDA (Gondara and Wang 2018), MissForest (Stekhoven and Buehlmann 2012), mean (Little and Rubin 1986), matrix factorization and k-NN imputation (Botstein et al. 2001). Concerning MICE and MissForest we used the default parameters discussed in the corresponding papers. For the other methods, we fine-tuned the hyper-parameters according to the corresponding literature to provide a fair comparison. In particular, we used $1\times 10^{-3}$ and $1\times 10^{-4}$ as learning rate for matrix factorization and MIDA with the latter being a 2-layer with 128 units per layer network like our embedding dimension. Finally, we let $k=5$ for the k-NN. The imputation accuracy for each dataset is presented in Table 2 for the scenario in which 30% of the entries are missing, while the results for all other levels of missingness (both in terms of RMSE and MAE) are presented in the supplementary material. We can see from Table 2 that the proposed GINN method obtains the best imputation performance in almost half of the datasets, being the second-best in almost all the remaining ones.

To provide a more schematic comparison, in Figure 4 we show the summary of those results for every level of missing data. In these histograms, we provide the number of times that each method achieves the best imputation (on average), i.e. the lowest RMSE in Figure 4(a) and MAE in Figure 4(b). Concerning the lower percentage of missing elements (10%, 20%) our method is almost on par with the best among the algorithms tested, i.e., MissForest. When those percentages increase our method brings a huge improvement over the state-of-the-art with the highest difference at 50% where our method significantly outperforms all other techniques. Aggregating the results obtained at 30% and 50% percentage of missing features, we have that our method is the best in 50% of the cases against the 27.5% of its best competitor MissForest, when looking at the MAE, and 47.5% against 20% for the RMSE. We defer a statistical analysis of these results to the next subsection, where we analyze also the results for a downstream predictive task.

2 Predictive Performance

Now we evaluate the performance of standard machine learning algorithms for classification, both binary and multi-class, trained on the various imputations analyzed previously. We consider 4 different classifiers: a k-NN classifier with $k=5$, regularized logistic regression, C-Support Vector Classification with an RBF kernel and a random forest classifier with 10 estimators and a maximum depth of 5. All hyper-parameters are initialized with their default values in the scikit-learn implementation. https://scikit-learn.org/stable/modules/classes.html

The classification accuracy is presented in Table 3 for the scenario in which 30% of the entries in the data matrix are missing, assuming MCAR, with a random forest classifier. In Figure 5 we show the summary of the results for each noise level and for each classifier. In these histograms we analyze the number of times each imputation technique allows the classifier to achieve the best average accuracy. In this comparison, we consider also the draws.

Our method outperforms competitors with every classification algorithm tested, and it has the highest number of wins, winning in $85.62\%$ of the cases. Our method worked very well when paired with SVC and Random forest, improving the accuracy for every percentage of missing elements. With the logistic regression and k-NN classifiers, at the lowest missing percentage (10%), our method is slightly below the state-of-the-art. As missing percentages increase, we have a huge improvement over the other competitors. This reflects in the findings of the previous Section and confirms the ability of our method of being very successful in the case of moderate to severely damaged datasets.

We corroborate the results of this and the previous section by performing a statistical analysis of the algorithms according to the guidelines in Demšar 2006. A Friedman rank test confirms that there are statistical significant differences both with respect to the RMSE of Figure 4 (p-value of $1.11e^{-11}$), and with respect to the accuracy of the classifiers in Figure 5 (e.g., p-value for random forest is $1.01e^{-10}$). A successive set of Nemenyi post-hoc tests further confirms statistical significant differences between GINN and all other methods for the random forest classification in Figure 5, and between GINN and all other methods except RF for the imputation results in Figure 4. The full set of rankings and of p-values for these tests can be found in the additional material for this paper.

3 Ablation study

To investigate how much each step described in Section 4 improves our method, we supervised the quality of imputation and convergence by starting from a basic autoencoder. As the starting point we used a 2-layer denoising autoencoder (DAE), obtained by setting the adjacency matrix to be the identity, obtaining a method similar to (Gondara and Wang 2018). Then we introduced the graph and the graph convolutional layer in our imputer (GINN) followed by the addition of the critic and the adversarial training (A-GINN), the skip connection (A-GINN skip) and finally the global attributes (A-GINN skip global).

Regarding imputation performances, Table 4 shows how the introduction of the graph and the graph-convolution operation over three randomly selected datasets makes a huge difference against a standard autoencoder. After that, each following step helps to refine even further the imputation accuracy. The reconstruction loss in Eq. (5), more precisely its logarithm, shows a similar behaviour, shown in Figure 6. After each step we have a better convergence. Similar results are obtained for all the other datasets.

The results in Figure 6 are shown with respect to the number of iterations. In the supplementary material, we also provide a similar analysis with respect to a fixed computational budget, while an analysis of the overall computational time when compared to the other algorithms is provided later on in Section 5.5.

4 Imputation over unseen data

In this section we test the ability of the model of imputing a new damaged portion of the dataset that was not available at training time. In order to impute these new values, we have first to inject the new data in the graph, adding nodes and edges. We compute a new similarity matrix for the new features (not considering the labels) and also the similarity between these new features and the older ones. Then we add the new nodes to the graph and the edges resulting from the double threshold procedure described in Section 4.1.

To evaluate this, we introduced missing values and evaluated the imputation performances with and without fine-tuning the model on a second randomly kept portion of the datasets. In Table 5 we show the MAE of the imputation over this new portion of dataset and compare against the other MDI algorithms. We used the same dataset and settings of the ablation study. The fine-tuned version consists of an additional 500 epochs of training over the new graph. It can be seen how our method is able to perform a state-of-the-art imputation on new unseen data without performing additional training and how it improves in case of a small number of additional optimization steps.

5 Evaluation on datasets with pre-existing missing values

To evaluate GINN in a real-world scenario, we compared its performance against the other state-of-the-art algorithms over three datasets with pre-existing missing values. When the data is not missing completely at random, the problem gets more complex, because there may exist relationships between the probability of a variable to be missing and other observed data. We consider two biomedical datasets: mammographic mass introduced in Elter et al. 2007 and cervical cancer by Fernandes et al. 2017 respectively with 4% and 13% of missing elements. We performed the imputation without considering the additional information of the label. Then we solved the binary classification task as done in Section 5.2. The third dataset is a time-series of air quality measurements (De Vito et al. 2008) with 13% of missing elements. Differently from the other two, missing values can also be found in the labels, making this a semi-supervised task. We select the three initial most damaged months as training data and perform imputation on all data, including the target variable. Then, we performed the downstream task of classifying, in the two successive months, the target variable discretized in three bins. A summary of the three datasets is provided in Table 6.

In Table 7 we show the average accuracy (computed over 10 trials) for every combination of imputation method and downstream classifier. For clarity, we highlight in bold the best result, and we underline the second-best one. While our proposed algorithm does not necessarily achieve the best accuracy overall, it can be seen from Table 7 that it is, in average, the most resilient to the choice of an external classifier. Overall, when combined with a logistic regression we obtain the best accuracy for the cervical cancer dataset (on par with several other algorithms), while we obtain the second-best result when combined with a random forest or an SVC in the other two cases. In order to highlight the resilience of the algorithm, in the supplementary material we provide a statistical analysis when the results from the different classifiers are aggregated.

Concerning computational performance, we provide execution times for all the algorithms (for simplicity, in the random forest case) in the supplementary material. Briefly, our algorithm is faster than alternative neural approaches (as already described in Section 5.3), but slower than non-neural approaches such as k-NN.

Conclusions and future work

In this paper we introduced a novel technique for missing data imputation, where we used a novel graph convolutional autoencoder to reconstruct the full dataset. We also describe several improvement to our technique, including the use of an adversarial loss, and the inclusion of global information from the dataset in the reconstruction phase. We show through an extensive numerical simulation that our method has good imputation performance, and the results are robust to the selection of an additional classifier later on. In experiments with a large level of artificial noise, our method is also shown to significantly outperform competitors.

Future work can consider the adoption of different graph neural architectures for the autoencoding process (such as those mentioned in Section 2), or the extension to other types of noisy data beyond vector-valued data and different types of similarity measures. In addition, in order to further improve accuracy and training time, we can think of training our imputation module together with a classification step in a end-to-end fashion.

Currently, the major drawbacks of our method are the need for computing the similarity matrix of the data, and the difficulty of performing mini-batching in the presence of graph-based data. Both problems are well-known in the corresponding literature, and in future work we plan on investigating techniques for speeding up similarity search and mini-batching on the graph to improve the computational complexity of the method.

References

References

Additional material

In Table 8 we provide the individual RMSE imputation values, where each row is one of the $20$ datasets and each column is an imputation method (abbreviations are explained in the main text). We separate the results with respect to the level of artificial corruption of the original dataset: the suffix $\_$xx means $xx\%$ of missing values that have been artificially added. These results are aggregated and commented in Fig. 4 and Tab. 3 of the main text. Results for MAE are similar and can be found on our online repository.

Detailed accuracy for regression/classification (Section 5.2)

In Table 9 we report the accuracy of the downstream classification/regression task for each choice of classification/regression technique when using a random forest technique (results are similar for other methods, and for brevity we provide them on our online repository). Like before, we separate the results with respect to the level of artificial corruption of the original dataset: the suffix $\_$xx means $xx\%$ of missing values that have been artificially added. These results are aggregated and commented in Fig. 5(c) and Tab. 4 of the main text.

Detailed results for the statistical tests

In Fig. 7 we provide the average rankings and p-values for the statistical test performed on the results of Section 5.1 of the paper (Nemenyi post-hoc tests on all pairs of algorithms), corresponding to Tab. 2. In Fig. 8 we instead provide average rankings and corresponding p-values for the statistical tests performed on the random forest classifier of Section 5.2, corresponding to Tab. 9. The results are discussed more in-depth in the main paper in Section 5.2.

Ablation study with a computational budget (Section 5.3)

In this section we replicate the ablation study of Section 5.3, but we evaluate the convergence of each variant with respect to a fixed computational budget, as shown in Fig. 9. As can be seen, our baseline method with GCN but no adversarial loss can converge to a significantly better result than a standard DAE in a fraction of the time. Including the adversarial loss slightly improves the final result, requiring however $3-4$ times the budget of the baseline method.

Detailed results for the evaluation on real-world datasets

In Figure 10 we report the scores obtained with a random forest classifier and the time in seconds needed for the imputation, including for GINN, the time needed for the similarity-graph construction. As can be seen, GINN’s imputation allows to achieve the best accuracy and consistency across all three datasets, even in the case where the imputed training labels are used to train the classification model. Regarding execution times, GINN is considerably faster than alternative neural approaches and missForest, but it is slower than the approaches that do not perform a training phase.

We corroborate the results by performing a statistical analysis. To focus on the impact of the imputation techniques, we compute the relative ranking for each combination of imputation method and classification algorithm, and we then aggregate the results with respect to the latter. A Friedman rank test confirms that there are statistical significant differences with respect to the accuracy of the classifiers with p-values of $1.6e^{-7}$, $1.3e^{-7}$, $1.2e^{-9}$ respectively for the Cervical cancer, Mammographic mass and Air quality datasets. A successive set of Nemenyi post-hoc tests, between all pairs of algorithms, further confirms statistical significant differences between GINN and all other methods as reported in Figure 11.