{'fraction_non_alphanumeric': 0.026910543466834858, 'fraction_numerical': 0.0124391183993833, 'mean_word_length': 4.8043522473052676, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MissForest - nonparametric missing value imputation for mixed-type data

Modern data acquisition based on high-throughput technology is often facing the problem of missing data. Algorithms commonly used in the analysis of such large-scale data often depend on a complete set. Missing value imputation offers a solution to this problem. However, the majority of available imputation methods are restricted to one type of variable only: continuous or categorical. For mixed-type data the different types are usually handled separately. Therefore, these methods ignore possible relations between variable types. We propose a nonparametric method which can cope with different types of variables simultaneously. We compare several state of the art methods for the imputation of missing values. We propose and evaluate an iterative imputation method (missForest) based on a random forest. By averaging over many unpruned classification or regression trees random forest intrinsically constitutes a multiple imputation scheme. Using the built-in out-of-bag error estimates of random forest we are able to estimate the imputation error without the need of a test set. Evaluation is performed on multiple data sets coming from a diverse selection of biological fields with artificially introduced missing values ranging from 10% to 30%. We show that missForest can successfully handle missing values, particularly in data sets including different types of variables. In our comparative study missForest outperforms other methods of imputation especially in data settings where complex interactions and nonlinear relations are suspected. The out-of-bag imputation error estimates of missForest prove to be adequate in all settings. Additionally, missForest exhibits attractive computational efficiency and can cope with high-dimensional data.


Introduction
Imputation of missing values is often a crucial step in data analysis. Many established methods of analysis require fully observed data sets without any missing values. However, this is seldom the case in medical and biological research today. The ongoing development of new and enhanced measurement techniques in these fields provides data analysts with challenges prompted not only by high-dimensional multivariate data where the number of variables may greatly exceed the number of observations but also by mixed data types where continuous and categorical variables are present. In our context categorical variables can arise as any kind ranging from technical settings in a mass spectrometer to a diagnostic expert opinion on a disease state. Additionally, such data sets often contain complex interactions and nonlinear relation structures which are notoriously hard to capture with parametric procedures.
Most prevalent imputation methods, like k nearest neighbours (KNNimpute, Troyanskaya et al. (2001)) for continuous data, saturated multinomial model (Schafer (1997)) for categorical data and multivariate imputation by chained equations (MICE, Van Buuren and Oudshoorn (1999)) for mixed data types depend on tuning parameters or specification of a parametric model. The choice of such tuning parameters or models without prior knowledge is difficult and might have a dramatic effect on a method's performance. Excluding MICE the above methods and the majority of other imputation methods are restricted to one type of variable. Furthermore, all these methods make assumptions about the distribution of the data or subsets of the variables, leading to questionable situations, e.g. assuming normal distributions.
The literature on mixed-type data imputation is rather scarce. Its first appearance was in the developing field of multiple imputation brought up by Rubin (1978). Little and Schluchter (1985) presented an approach based on maximum likelihood estimation combining the multivariate normal model for continuous and the Poisson/multinomial model for categorical data. This idea was later on extended in the book of Little and Rubin (1987). See also Li (1988), Rubin and Schafer (1990) and Schafer (1997). A more refined method to combine different regression models for mixed-type data was proposed by Van Buuren and Oudshoorn (1999) using chained equations. The conditional model in MICE can be specified for the missing data in each incomplete variable. Therefore no multivariate model covering the entire data set has to be specified. However, it is assumed that such a full multivariate distribution exists and missing values are sampled from conditional distributions based on this full distribution (for more details see Section 3). Another similar method using variable-wise conditional distributions was proposed by Raghunathan et al. (2001) called sequential regression multivariate imputation. Unlike in MICE the predictors must not be incomplete. The method is focussed on survey data and therefore includes strategies to incorporate restrictions on subsamples of individuals and logical bounds based on domain knowledge about the variables, e.g., only women can have a number of pregnancies recorded.
Our motivation is to introduce a method of imputation which can handle any type of input data and makes as few as possible assumptions about structural aspects of the data. Random forest (RF, Breiman (2001)) is able to deal with mixed-type data and as a nonparametric method it allows for interactive and nonlinear (regression) effects. We address the missing data problem using an iterative imputation scheme by training a RF on observed values in a first step, followed by predicting the missing values and then proceeding iteratively. Mazumder et al. (2010) use a similar approach for the matrix completion problem using a soft-thresholded SVD iteratively replacing the missing values. We choose RF because it can handle mixed-type data and is known to perform very well under barren conditions like high dimensions, complex interactions and nonlinear data structures. Due to its accuracy and robustness RF is well suited for the use in applied research often harbouring such conditions. Furthermore, the RF algorithm allows for estimating out-of-bag (OOB) error rates without the need for a test set. For further details see Breiman (2001).
Here we compare our method with k-nearest neighbour imputation (KNNimpute, Troyanskaya et al. (2001)) and the Missingness Pattern Alternating Lasso (MissPALasso) algorithm by Städler and Bühlmann (2010) on data sets having continuous variables only. For the cases of categorical and mixed type of variables we compare our method with the MICE algorithm by Van Buuren and Oudshoorn (1999) and a dummy variable encoded KNNimpute. Comparisons are performed on several data sets coming from different fields of life sciences and using different proportions of missing values.
We show that our approach is competitive to or outperforms the compared methods on the used data sets irrespectively of the variable type composition, the data dimensionality, the source of the data or the amount of missing values. In some cases the decrease of imputation error is up to 50%. This performance is typically reached within only a few iterations which makes our method also computationally attractive. The OOB imputation error estimates give a very good approximation of the true imputation error having on average a proportional deviation of no more than 10 -15%. Furthermore, our approach needs no tuning parameter, hence, is easy to use and needs no prior knowledge about the data.

Approach
We assume X = (X 1 , X 2 , . . . , X p ) to be a n × p-dimensional data matrix. We propose using a random forest to impute the missing values due to its earlier mentioned advantages as a regression method. The random forest algorithm has a built-in routine to handle missing values by weighting the frequency of the observed values in a variable with the random forest proximities after being trained on the initially mean imputed data set (Breiman (2001)). However, this approach requires a complete response variable for training the forest.
Instead, we directly predict the missing values using a random forest trained on the observed parts of the data set. For an arbitrary variable X s including missing values at entries i Note that x To begin, make an initial guess for the missing values in X using mean imputation or another imputation method. Then, sort the variables X s , s = 1, . . . , p according to the amount of missing values starting with the lowest amount. For each variable X s the missing values are imputed by first fitting a random forest with response y   The stopping criterion γ is met as soon as the difference between the newly imputed data matrix and the previous one increases for the first time with respect to both variable types, if present. Here, the difference for the set of continuous variables N is defined as and for the set of categorical variables F as where #NA is the number of missing values in the categorical variables. After imputing the missing values the performance is assessed using the normalised root mean squared error (NRMSE, Oba et al. (2003)) for the continuous variables which is defined by where X true is the complete data matrix and X imp the imputed data matrix. We use mean and var as short notation for empirical mean and variance computed over the continuous missing values only. For categorical variables we use the proportion of falsely classified entries (PFC) over the categorical missing values, ∆ F . In both cases good performance leads to a value close to 0 and bad performance to a value around 1. When a RF is fit to the observed part of a variable we also get an OOB error estimate for that variable. After the stopping criterion γ was met we average over the set of variables of the same type to approximate the true imputation errors. We assess the performance of this estimation by comparing the absolute difference between true imputation error and OOB imputation error estimate in all simulation runs.

Methods
We compare missForest with four methods on ten different data sets where we distinguish between situations with continuous variables only, categorical variables only and mixed variable types.
The most well-known method for imputation of continuous data sets especially in the field of gene expression analysis is the KNNimpute algorithm by Troyanskaya et al. (2001). A missing value variable X j is imputed by finding its k nearest observed variables and taking a weighted mean of these k variables for imputation. Thereby, the weights depend on the distance of the variable X j . The distance itself is usually chosen to be the Euclidean distance.
When using KNNimpute the choice of the tuning parameter k can have a large effect on the performance of the imputation. However, this parameter is not known beforehand. Since our method includes no such parameter we implement a cross-validation (see Algorithm 2) to obtain a suitable k.
Require: X an n × p matrix, number of validation sets l, range of suitable number of nearest neighbours K 1: X CV ← initial imputation using mean imputation; 2: for t in 1, . . . , l do 3: X CV mis,t ← artificially introduce missing values to X CV ; 4: for k in K do 5: X CV KN N,t ← KNN imputation of X CV mis,t using k nearest neighbours; 6: ε k,t ← error of KNN imputation for k and t; 7: end for 8: end for 9: k best ← argmin k 1 l l t=1 ε k,t ; 10: X imp ← KNN imputation of X using k best nearest neighbours.
In the original paper of Troyanskaya et al. (2001) the data was not standardized before applying the KNNimpute algorithm. This constitutes no issue in the case of gene expression data because such data generally consists of variables on similar scales. However, we are applying the KNNimpute algorithm to data sets with varying scales in the variables. To avoid variance based weighting of the variables we scale them to a unit standard deviation. We also center the variables at zero. After imputation the data is retransformed such that the error is computed on the original scales. This last step is performed because missForest does not need any transformation of the data and we want to compare the performance of the methods on the original scales of the data.
Another approach for continuous data, especially in the case of high-dimensional normal data matrices, is presented by Städler and Bühlmann (2010) using an EM-type algorithm. In their Missingness Pattern Alternating Imputation and l 1 -penalty (MissPALasso) algorithm the missing variables are regressed on the observed ones using the lasso penalty by Tibshirani (1996). In the following E step the obtained regression coefficients are used to partially update the latent distribution. The MissPALasso has also a tuning parameter λ for the penalty. As with KNNimpute we use cross-validation to tune λ (cf. Algorithm 2). When applying MissPALasso the data is standardized as regularization with a single λ requires the different regressions to be on the same scale.
In the comparative experiments with categorical or mixed-type variables we use the MICE algorithm by Van Buuren and Oudshoorn (1999) based on the multivariate multiple imputation scheme of Schafer (1997). In contrast to the latter the conditional distribution for the missing data in each incomplete variable is specified in MICE, a feature called fully conditional specification by Van Buuren (2007). However, the existence of a multivariate distribution from which the conditional distribution can be easily derived is assumed. Furthermore, iterative Gibbs sampling from the conditional distributions can generate draws from the multivariate distribution. We want to point out that MICE in its default setup is not mainly intended for simple missing value imputation. Using the multiple imputation scheme MICE allows for assessing the uncertainty of the imputed values. It includes features to pool multiple imputations, choose individual sampling procedures and allows for passive imputation controlling the sync of transformed variables. In our experiments we used MICE with either linear regression with normal errors or mean imputation for continuous variables, logistic regression for binary variables and polytomous logistic regression for categorical variables with more than two categories.
For comparison across different types of variables we apply the KNNimpute algorithm with dummy coding for the categorical variables. This is done by coding a categorical variable X j into m dichotomous variablesX j,m ∈ {−1, 1}. Application of the KNNimpute algorithm for categorical data can be summarized as: 1. Code all categorical variables into {−1, 1}-dummy variables; 2. standardize all variables to mean 0 and standard deviation 1; 3. apply the cross-validated KNNimpute method from Algorithm 2; 4. retransform the imputed data matrix to the original scales; 5. code the dummy variables back to categorical variables; 6. computed the imputation error.
For each experiment we perform 50 independent simulations where 10%, 20% or 30% of the values are removed completely at random. Each method is then applied and the NRMSE, the PFC or both are computed (see Section 2). We perform a paired Wilcoxon test of the error rates of the compared methods versus the error rates of missForest. In addition, the OOB error estimates of missForest is recorded in each simulation.

Continuous variables only
First, we focus on continuous data. We investigate the following four publicly available data sets: • Isoprenoid gene network in A. thaliana: This gene network includes p = 39 genes each with n = 118 gene expression profiles corresponding to different experimental conditions. For more details on this data set see Wille et al. (2004).
• Voice measures in Parkinson's patients: The data described by Little et al. (2008) contains a range of biomedical voice measurements from 31 individuals, 23 with Parkinson's disease (PD). There are p = 22 particular voice measurements and n = 195 voice recordings from these individuals. The data set also contains a response variable giving the health status. Dealing only with continuous variables the response was removed from the data. We will return to this later on.
• Shapes of musk molecules: This data set describes 92 molecules of which 47 are musks and 45 are non-musks. For each molecule p = 166 features describe its conformation, but since a molecule can have many conformations due to rotating bonds, there are n = 476 different lowenergy conformations in the set. The classification into musk and non-musk molecules is removed.  Figure 1: Continuous data. Average NRMSE for KNNimpute (grey), MissPALasso (white) and miss-Forest (black) on four different data sets and three different amounts of missing values, i.e., 10%, 20% and 30%. Standard errors are in the order of magnitude of 10 −4 . Significance levels for the paired Wilcoxon tests in favour of missForest are encoded as "*" <0.05, "**" <0.01 and "***" <0.001. If the average error of the compared method is smaller than that of missForest the significance level is encoded by a hash (#) instead of an asterisk. In the lowermost data set results for MissPALasso are missing due to the methods limited capability with regard to high dimensions.
• Insulin gene expression: This high-dimensional data set originates from an analysis by Wu et al. (2007) of vastus lateralis muscle biopsies from three different types of patients following insulin treatment. The three types are insulin-sensitive, insulin-resistant and diabetic patients. The analysis involves p = 12 626 genes whose expression levels were measured from n = 110 muscle biopsies. Due to computation time we only perform 10 simulations instead of 50.
Results are given in Figure 1. We can see that missForest performs well, sometimes reducing the average NRMSE by up to 25% with respect to KNNimpute. In case of the musk molecules data the reduction is even above 50%. The MissPALasso performs slightly better than missForest on the gene expression data. However, there are no results for the MissPALasso in case of the Insulin data set because the high dimension makes computation not feasible.
For continuous data the missForest algorithm typically reaches the stopping criterion quite fast needing about 5 iterations. The imputation takes about 10 times as long as performing the cross-validated KNNimpute where {1, . . . , 15} is the set of possible numbers of neighbours. For the Insulin data set an imputation takes on average 2 hours on a customary available desktop computer.

Categorical variables only
We also consider data sets with only categorical variables. Here, we use the MICE algorithm described in Section 3 instead of the MissPALasso. We use a dummy implementation of the KNNimpute algorithm to deal with categorical variables (see Section 3). We apply the methods to the following data sets: • Cardiac single photon emission computed tomography (SPECT) images: Kurgan et al. (2001) discuss this processed data set summarizing over 3000 2D SPECT images from n = 267 patients in p = 22 binary feature patterns. • Promoter gene sequences in E. coli: The data set contains sequences found by Harley and Reynolds (1987) for promoters and sequences found by Towell et al. (1990) for non-promoters totalling n = 106. For each candidate a sequence of 57 base pairs was recorded. Each variable can take one of four DNA nucleotides, i.e., adenine, thymine, guanine or cytosine. Another variable distinguishes between promoter and non-promoter instances.
• Lymphography domain data: The observations were obtained from patients suffering from cancer in the lymphatic of the immune system. For each of the n = 148 lymphoma p = 19 different properties were recorded mainly in a nominal fashion. There are nine binary variables. The rest of the variables have three or more levels.
In Figure 2 we can see that missForest is always imputing the missing values better than the compared methods. In some cases -namely for the SPECT data -the decrease of PFC compared to MICE is up to 60%. However, for the other data sets the decrease is less pronounced ranging around 10 -20% -but there still is a decrease. The amount of missing values on the other hand seems to have only a minor influence on the performance of all methods. Except for MICE on the SPECT data, error rates remain almost constant increasing only by 1 -2%. We pointed out earlier that MICE is not primarily tailored for imputation performance but offers additional possibilities of assessing uncertainty of the imputed values due to the multiple imputation scheme. Anyhow, the results using the cross-validated KNNimpute (see Algorithm 2) on the dummy-coded categorical variables is surprising. The imputation for missForest needs on average 5 times as long as a cross-validated imputation using KNNimpute.

Mixed-type variables
In the following we investigate four data sets where the first one has already been introduced, i.e. musk molecules data including the categorical response yielding the classification. The other data sets are: • Proteomics biomarkers for Gaucher's disease: Gaucher's disease is a rare inherited enzyme deficiency. In this data set Smit et al. (2007) present protein arrays for biomarkers (p = 590) from blood serum samples (n = 40). The binary response distinguishes between disease status.
• Gene finding over prediction (GFOP) peptide search: This data set comprises mass-spectrometric measurements of n = 595 peptides from two shotgun proteomics experiments on the nematode Caenorhabditis elegans. The collection of p = 18 biological, technical and analytical variables had the aim of novel peptide detection in a search on an extended database using established gene prediction methods.
• Children's Hospital data: This data set is the product of a systematic long-term review of children with congenital heart defects after open-heart surgery. Next to defect and surgery related variables also long-term psychological adjustment and health-related quality of life was assessed. After removing observations with missing values the data set consists of n = 55 patients and p = 124 variables of which 48 are continuous and 76 are categorical. For further details see Latal et al. (2009).
The results of this comparison are given in Figure 3. We can see that missForest performs better than the other two methods, again reducing imputation error in many cases by more than 50%. For the GFOP data, KNNimpute has a slightly smaller NRMSE than missForest but makes twice as much error on the categorical variables. Generally, with respect to the amount of missing values the NRMSE tends to have a greater variability than the PFC which remains largely the same.
The imputation results for MICE on the Children's Hospital data have to be treated cautiously. Since this data set contains ill-distributed and nearly dependent variables, e.g., binary variables with very few observations in one category, the missingness pattern has a direct influence on the operability of the MICE implementation in the statistical software R. The imputation error illustrated in Figure 3 was computed from 50 successful simulations by randomly generating missingness patterns, which did not include only complete cases or no complete cases at all within the categories of the variables. Therefore, the actual numbers of simulations were larger than 50 for all three missing value amounts. Furthermore, nearly dependent variables were removed after each introduction of missing values. This leads to an average of 7 removed variables in each simulation. Due to this ad-hoc manipulation for making the MICE implementation work, we do not report significance statements for the imputation error.

Estimating imputation error
In each experiment we get for each simulation run an OOB estimate for the imputation error. In Figure  4 the differences of true imputation error, err true , and OOB error estimates, err OOB , are illustrated for the continuous and the categorical data sets. Also, the mean of the true imputation error and the OOB error estimate over all simulations is depicted.
We can see that for the Isoprenoid and Musk data sets the OOB estimates are very accurate only differing from the true imputation error by a few percents. In case of the Parkinson's data set the OOB estimates exhibit a lot more variability than in all other data sets. However, on average the estimation is comparably good. For the categorical data sets the estimation accuracy behaves similarly over all scenarios. The OOB estimates tend to underestimate the imputation error with increasing amount of missing values. Apparently, the absolute size of the imputation error seems to play a minor role in the accuracy of the OOB estimates which can be seen nicely when comparing the SPECT and the Promoter data.

Computational efficiency
We assess the computational cost of missForest by comparing the runtimes of imputation on the previous data sets. Table 1 shows the runtimes in seconds of all methods on the analyzed data sets. We can see that KNNimpute is by far the fastest method. However, missForest runs considerably faster than MICE and the MissPALasso. In addition, applying missForest did not require antecedent standardization of the data, laborious dummy coding of categorical variables nor implementation of CV choices for tuning parameters. Significance levels for the paired Wilcoxon tests in favour of missForest are encoded as "*" <0.05, "**" <0.01 and "***" <0.001. If the average error of the compared method is smaller than that of missForest the significance level is encoded by a hash (#) instead of an asterisk. Note that, due to ill-distribution and near dependence in the Child hospital data, the results for MICE have to be treated with caution (see Section 4.3).
There are two possible ways to speed up computation. The first one is to reduce the number of trees grown in each forest. In all comparative studies the number of trees was set to 100 which offers high precision but increased runtime. In Table 2 we can see that changing the number of trees in the forest has a stagnating influence on imputation error, but a strong influence on computation time which is approximately linear in the number of trees.
The second one is to reduce the number of variables randomly selected at each node (m try ) to set up the split. Table 2 shows that increasing m try has limited effect on imputation error, but computation time is strongly increased. Note that for m try = 1 we do not longer have a random forest, since there is no more choice between variables to split on. This leads to a much higher imputation error, especially for the cases with low numbers of bootstrapped trees. We use for all experiments √ p as default value, e.g., in the GFOP data this equals 4.

Conclusion
Our new algorithm, missForest, allows for missing value imputation on basically any kind of data. In particular, it can handle multivariate data consisting of continuous and categorical variables simultaneously. MissForest has no need for tuning parameters nor does it require assumptions about distributional aspects of the data. We show on several real data sets coming from different biological and medical fields that missForest outperforms established imputation methods like k-nearest neighbours imputation or multivariate imputation using chained equations. Using our OOB imputation error estimates miss-Forest offers a way to assess the quality of an imputation without the need of setting aside test data nor performing laborious cross-validations. For subsequent analysis these error estimates represent a mean Lymphography SPECT Promoter q q q q q q q q q q mean(err OOB) mean(errtrue) Figure 4: Difference of true imputation error err true and OOB imputation error estimate err OOB for the continuous data sets (left) and the categorical data sets (right) and three different amounts of missing values, i.e., 0.1, 0.2 and 0.3. In each case the average err true (circle) and the average err OOB (plus) over all simulations is given.  Table 2: Average imputation error (NRMSE/PFC in percent) and runtime (in seconds) with different numbers of trees (n tree ) grown in each forest and variables tried (m try ) at each node of the trees. Here, we consider the GFOP data set with artificially introduced 10% of missing values. For each comparison 50 simulation runs were performed using always the same missing value matrix for all numbers of trees/randomly selected variables for a single simulation.