Concept Drift Detection: Dealing with Missing Values via Fuzzy Distance Estimations

In an evolving data stream, the distribution of available training samples may vary over time [1]. Consider a topological space feature space denoted as $\mathcal{X}\subseteq\mathbb{R}^{n}$, where $n$ is the dimensionality of the feature space. A tuple $(X,y)$ denotes a data instance, where $X\in\mathcal{X}$ is the feature vector, $y\in\{y_{1},\ldots,y_{c}\}$ is the class label, and $c$ is the number of classes. A data stream, denoted as $\mathcal{D}$, can then be represented as a sequence of data instances. Time windows strategy chunks the sequence in a time interval as a batch. These batches are denoted as $D_{T_{i}}\in\mathcal{D}$, where $T_{i}$ is a given time interval that defines the time window. A concept drift has occurred between two time windows $T_{i}$ and $T_{i+1}$ if the joint probability of $X$ and $y$ is different, that is, $p_{T_{i}}(X,y)\neq p_{T_{i+1}}(X,y)$ [5, 6].

There are two approaches to estimating the density discrepancy: parametric and nonparametric [32, 31]. The parametric approach assumes the data is drawn from a known distribution, with the major challenge being how to estimate the parameters to reach the maximum likelihood [31]. Nonparametric approaches assume that it is too difficult to devise an analytic expression of the target distribution’s probability density function [33]. Therefore, the density is estimated empirically rather than by matching the data with any given parametric family. Given that concept drift detection inherently concerns evolving data distributions, it is difficult to rely on one distribution family for an accurate description. Therefore, most drift detection methods are nonparametric [1].

To quantify the distribution discrepancy between two multivariate sample sets, some researchers have proposed new test statistics, while others have proposed novel mapping methodologies that convert multivariate observations into univariate observations. For example, NN-DVI [34] and MMD [35] construct new statistical variables for comparing distributions to determine the differences. The performance of these methods is excellent, but that performance usually comes with a high computational cost. Alternatively, if the distribution of the proposed statistics is unclear, a Monto Carlo method can be applied to approximate the significance level. Examples include kdqTree combines with Kullback-Leibler divergence [36], competence model combined with total variation [7], QuantTree combined with total variation, or QuantTree combined with Chi-square statistics [37]. On yet another track, the solution involves building a model to map the discrepancy in distributions with powerful statistical tools, such as Pearson’s Chi-square [31]. The general premise of these approaches is to convert multivariable data samples into a univariable two-sample test problem, such as an incremental Kolmogorov–Smirnov test [38] or a multi-Wald-Wolfowitz test [39]. However, no research has considered the missing values for concept drift detection.

Missing values are merely observations that were not recorded, perhaps due to human error or a sensor malfunction, which is very common in real-world data [40, 16, 15]. There are three main categories of missing values based on the correlation between the observations and the values that are missing. These are: missing completely at random (MCAR), missing at random (MAR) and missing not at random (MNAR) [27, 41]. A recent survey by Santos et al. [15] provides a set of notations and definitions that demonstrate the difference between each category quite well, as follows:

Assume we have a data set with $m$ number of instances and $n$ number of features. The values of the data set are represented as an $m\times n$ matrix as $X=\{x_{i,j}\}_{i,j}^{m,n}$. Plus, a missing data indicator matrix is defined as an $m\times n$ zero-one matrix, denoted as $M$. $M$ is represented as $M=\{I_{i,j}^{miss}\}_{i,j}^{m,n}$, where $I_{i,j}^{miss}$ is the missing value indicator and equals 1 if $x_{i,j}$ is missing, and 0 otherwise. The relationship describing the missing data is defined as a conditional probability function $p(M|X,\xi)$, where $\xi$ are the parameters of the missing values [26]. In practice, $\xi$ is not important, only the relation between $M$ and the components of $X$ are important, i.e., $X=(X^{complete},X^{miss})$, where $X^{complete}$ is the observations without missing values, and $X^{miss}$ is the the observations with missing values.

MCAR means that a missing value has nothing to do with its hypothetical value or with the values of other variables [15, 16]. The probability of missingness depends only on the parameters $\xi$, i.e., $p(M=1|X,\xi)=p(M=1|\xi)$ [15].

MAR means the missing values are not related to themselves, but they do depend on other features somehow. In other words, the conditional probability of missing values is $p(M=1|X,\xi)=p(M=1|X^{complete},\xi)$ [5].

With MNAR, the missing values depend on both the observed and the unobserved information – $X^{complete}$ and $X^{miss}$ – and the conditional probability of the missing values cannot be simplified, i.e., $p(M=1|X,\xi)=p(M=1|X^{complete},X^{miss},\xi)$ [15]. Conventional methods of handling missing values include complete case analysis, available case analysis, dummy variable adjustment, and imputation [26].

Imputation is currently the most widely used method for handling missing values [16]. Like dummy variable adjustment, missing values are replaced but this time with plausible values, not constants or a mean. The model is then trained just as though no values were missing [26]. However, there are two serious problems with most imputation methods. First, variances tend to be underestimated, which can lead to bias in the learning parameters [26]. Second, the standard error calculations presume that all data are real, which means the inherent uncertainty and sampling variability in the imputed values is not taken into account [26]. Of course, this results in confidence intervals that are too narrow.

Pearson’s Chi-square Test has been widely used to determine whether there are significant differences between expected frequencies and observed frequencies [30]. When there is no significant difference (i.e., when the null hypothesis is true), the statistic follows a chi-square distribution. Thus, the premise of the test is to assume the null hypothesis is true and then evaluate how likely a specific observation would be.

The standard process of the chi-square test is to use sample data to find: the degrees of freedom, the expected frequencies, the test statistic, and the p-value associated with the test statistic [30]. Given a contingency table with $i$ rows and $j$ columns of categorical variables, the degrees of freedom are equal to

The expected frequency counts are computed separately for each level of one categorical variable at each level of the other categorical variable. The $i$-th and $j$-th expected frequencies of the contingency table are calculated with the equation

where $n_{i}$ is the sum of the frequencies for all columns in row $i$, $n_{j}$ is the sum of the frequencies for all rows of columns $j$, and $n$ is the sum of all rows and columns. The test statistic is a chi-square random variable $\chi^{2}$ defined by

where $O_{i,j}$ is the observed frequency count at row $i$ and column $j$, and $E_{i,j}$ is the expected frequency count at row $i$ and column $j$. The p-value is the probability of observing a sample statistic as extreme as the test statistic. Since the p-value is a $\chi^{2}$ test statistic, it can be computed with the chi-square probability distribution function.

Pearson’s chi-square test should be used with a sufficiently large sample set. According to the central limit theorem, an $\chi^{2}$ distribution is the sum of $O_{i,j}$ independent random variables with a finite mean and variance that converges to a normal distribution for large $O_{i,j}$. For practical purposes, Box et al. [30] claim that, for $O_{i,j}>50$ and $E_{i,j}>5$, the distribution of the estimated test statistics is sufficiently close to a normal distribution to ignore the difference.

Problem Statement. Recall the drift detection problem outlined in Section II-A, i.e., Let $\mathcal{A}$ and $\mathcal{B}$ be random variables defined on a topological space $\mathcal{X}\subseteq\mathbb{R}^{n}$ at different periods in a data stream with respect to $p_{\mathcal{A}},p_{\mathcal{B}}\in\mathcal{P}(\mathcal{X})$, where $\mathcal{P}(\mathcal{X})$ consists of all Borel probability measures on $\mathcal{X}$. Given two observation sets drawn from $\mathcal{A}$ and $\mathcal{B}$, $A=\{X_{i}^{\mathcal{A}}\}_{i=1}^{m_{1}}$ and $B=\{X_{i}^{\mathcal{B}}\}_{i=1}^{m_{2}}$, the problem is: How confident are we to claim that $\mathcal{A}=\mathcal{B}$ based on $A$ and $B$?

To simplify this problem, it is commonly assumed that the observations $A$, $B$ are i.i.d., which makes the research objective equivalent to a two-sample test problem. The general procedure is to measure the difference between the observations $A$ and $B$ (sample sets), then infer the significance of the difference between $\mathcal{A}$ and $\mathcal{B}$ (population).

Remark: Recall the definition of concept drift, $p_{T_{i}}(X,y)\neq p_{T_{i+1}}(X,y)$, a data instance is defined as a tuple of ($X$, $y$). For supervised learning setting, the target variable $y$ can be considered as the $(n+1)$-th feature of the observations while for unsupervised setting, the drift detection will be equivalent to covariate shift detection [42].

The most intuitive method for capturing this difference is a histogram. A histogram is a set of intervals, i.e., bins, and density is estimated by counting the number of observations in each bin. A general framework for using a histogram and Pearson’s chi-square test to detect concept drift is shown in Fig. 2.

The next question is how to partition the feature space as bins to build a histogram. In the literature, the partitions of most multivariate histograms are based on a grid or a tree structure [37, 31]. We have used clusters because, given our problem, this method makes it easy to determine which bin an observation belongs to. More specifically, with grid- and tree-based histograms, observations are allocated to a bin according to the value of each feature, whereas, with cluster-based histograms, the bin is determined by the distance between the observations and the centroids, as demonstrated in Fig. 3.

In practice, there are three advantages to using clusters as histogram bins.

1) The shape of the bins is flexible. Both grid and tree-based histograms have hyper-rectangular-shaped bins. However, the shape of the bins in a cluster-based histogram dynamically adapts to the data distribution, which is beneficial when trying to maintain a similar number of observations per bin and not leaving any empty bins.

2) Missing values in an observation do not need to be estimated one-by-one with cluster-based histograms. Rather, only the observation-to-centroid distances need to be modeled to determine the appropriate observation-to-cluster membership. In addition, the error of the distance estimation can be characterized by fuzzifying the distance then using that fuzzy distance to calculate the observation’s degree of membership to a cluster.

3) Since our algorithm is task-oriented imputation rather than value-oriented, our drift detection algorithm has no assumption on the data type, i.e., both numerical and categorical data can be addressed.

This section sets out our formulation for a distance learning method to estimate observation-to-observation distances with missing values. The intuition behind the idea is to retrieve the observations without missing values, denoted as $X^{complete}$, and artificially introduce null values as a mask to simulate the missing values, that is,

where $f_{mask}(X)$ denotes the missing value masking function. Since our target variable $d$ is known, i.e., the actual paired distances between the retrieved observations, we can use the masked samples $X^{mask}$ and their actual distances $d$ as a training set to build a regression model, denoted as $L_{mask}$. The standard deviation $\sigma_{mask}$ of $L_{mask}$ can be used as a coefficient for fuzzifying observation-to-observation distances.

The workflow of this masked distance learning method comprises three major steps, as shown in Fig. 4. Put simply, masked distance learning asks and answers three questions: 1) How must the function $f_{mask}(X)$ be defined so that the observations with missing values have the same distribution as the original sample set? 2) Which method of building the learning model $L_{mask}$ will minimize the prediction errors? 3) How should the error of the $L_{mask}$ be estimated? Step 1-3 are proposed to address these questions.

Step 1 is a sample selection process. Retrieving observations with no missing values can be complete by a for loop over the sample set. The runtime complexity is $\mathcal{O}(mn)$, where $m$ is the number of observations, and $n$ is the dimension of the feature space. The number of observations with missing values is then counted and divided by the total number of observations to produce a missing rate vector, denoted as $V_{r}=\{r_{i}\}_{i=1}^{n}$, where $r_{i}$ is the missing rate of the $i$-th feature. At this stage, we assume the missing values are MCAR. Based on this assumption, the probability mass function of the masking variable for each feature is

and we have an masked observations as $X_{1}^{mask}=\{x_{i}^{mask}\}_{i=1}^{n}$. The actual pairwise distances of $X^{complete}$ for Step 2 is also calculated at this time.

Step 2 is the prediction model training. In this step, we want to build a learning model that maps the masked observations to the actual pairwise distance while minimizing the error between the predictions and the actual values, i.e.,

where $\mathcal{H}$ is the hypothesis space, and $h\in\mathcal{H}$ is the hypothesis of mapping observations to a real value $h:X\rightarrow d$. Few existing data imputation methods are applied to the masked observations, and the masked distances are calculated as new training features. Then these new training features, along with the actual distances prior to masking, are used to train a masked distance learning model. In our solution, we use the gradient boosted decision tree model to perform the learning task, which is a non-linear regression model. We consider that non-linear regression models often outperform linear models in complicated cases.

The basic assumption is that the true distance is close to the imputed distances but which imputation method provides the best results depends on the data. For example, mean replacement is the best option for or Gaussian-distributed features, whereas most frequent fills could be the best for Poisson-distributed features. So, to make the optimal choice, the imputed distances are passed to a regression model and learned. The masked training set contains $\frac{|X^{mask}|^{2}}{2}$ instances for training, where $|\cdot|$ denotes the cardinality of the set. On the MCAR assumption, the masked distance learning model should perform similarly across the entire sample set.

Step 3 is to estimate the performance of the trained distance learning model $L_{mask}$. According to the assumption of linear regression, i.e., that $d_{i}=a+\theta X_{i}+\epsilon_{i}$ where the error item $\epsilon\sim\mathcal{N}(0,\sigma^{2})$, we can use the error item $\epsilon$ (residual) to fuzzify the predicted distance. To estimate the variance, the residual standard error calculated by the equation below and then stored,

For each predicted distance $\hat{d}_{i}$, we have high confidence that the true distance is located in the interval $[\hat{d}_{i}-3\sigma_{mask},\hat{d}_{i}+3\sigma_{mask}]$ according to the three-sigma rule. To avoid sampling bias, we apply cross-validation to the $X^{mask}$, which splits $X^{mask}$ into kFold paired subsets $\{X_{train_{i}}^{mask},X_{test_{i}}^{mask}\}_{i=1}^{n_{split}}$, after which the $\sigma_{mask}$ is computed, where the $n_{split}$ is the number of split,.

Alg. 1. presents the pseudocode of masked distance learning. The inputs are the observations, a default regression model, a parameter for cross-validation, and a set of data imputation methods. The outputs are a trained $L_{mask}$ and its residual standard deviation $\sigma_{mask}$. The overall runtime complexity is $\mathcal{O}(n_{split}m^{2}n)$ broken down as follows. The complexity for Lines 1-3 to initialize the data set statistics is $\mathcal{O}(mn)$. The pairwise distance has a complexity of $\mathcal{O}(m^{2}n)$ without any optimization. The worst case for the observation masking, i.e., when every value is missing, is $\mathcal{O}(mn)$. The worst case for the data imputation is multiple iterative imputations, which would be $\mathcal{O}(m^{2}n)$. The kFold training-testing complexity is $\mathcal{O}(n_{split}m^{2}n)$. Calculating standard deviations is $\mathcal{O}(m)$. Therefore, the runtime complexity for Alg. 1. is feasible within the complexity of $\mathcal{O}(n_{split}m^{2}n)$

After building the missing value distance prediction model $L_{mask}$ and obtaining the standard deviation $\sigma_{mask}$ of the errors, we can now start estimating the frequency of observations in histogram bins. For any two observations with missing values, we use the same masking function to convert them then pass them to $L_{mask}$. The returned result $\hat{d}$ with the standard deviation $\sigma_{mask}$ will be used to construct the fuzzy distance to present its membership to a bin. In fact, there are two options with the predicted distance $\hat{d}$. One is simply to use the predictions $\hat{d}$ as crisp distances without considering the errors of $L_{mask}$. The observations would then be assigned to histogram bins in a many-to-one manner, which means each observation would only be counted once by a bin. Another option is to fuzzify the distance according to the $\sigma_{mask}$, that is, to consider the estimated distance as a fuzzy number, as shown in Fig. 5.

The workflow of fuzzy-weighted frequency is as follows. Step 1) model the estimated distance between $X_{i}$ and $C_{z}$ as fuzzy distance, $\tilde{d}_{i,z}$; Step 2) calculate the degree of membership $\mu_{i,k}$ of every observation to each bin based on the fuzzy distances. Step 3) normalize the degree of membership as observation’s weight. Step 4) Accumulate the weights in each bin as the frequency. The novelty of the FWF compare to fuzzy clustering is that the degree of membership is calculated based on fuzzy distance rather than crisp distance.

In this case, the sample set $X$ is the universe of discourse. The value $\mu(X_{i})$ for each $X_{i}\in X$ is the grade of membership of $X_{i}$ in $(X,\mu)$. The function $\mu(X_{i})$ is the membership function of a fuzzy set $C=(X,\mu)$, where a fuzzy set is a cluster denoted as $C_{k}$. When measuring the frequency of observations, the fuzzy distance will be used to determine the degree of membership of an observation to each clusters.

We propose two types of fuzzy numbers to present the estimated distance: the Triangular fuzzy numbers (TFN) and Gaussian fuzzy number (GFN). We chose the TFN as a baseline for comparison because it is the most widely used [12], and we chose the GFN as it is the most reasonable since the predictions and errors are assumed to follow a Gaussian distribution. Once the fuzzy distance formulation is determined, we can leverage them to calculate the degree of membership. The contingency table is then built with the degree of membership.

In general, the membership degree of an observation with missing values to a cluster is determined by its estimated distances to the cluster centroids and the standard deviation of the errors. This is formulated by the membership function $\mu\left(X_{i},C_{z}\right)$, namely $\mu_{i,z}$. Without any prior knowledge of the sample set, we assume the observations all have equal weights of $1$. Therefore, the fuzzy weighting function can be formulated as a normalized degree of membership, denoted as

Fig. 6 illustrates how to calculate the weights based on the Triangular fuzzy distance. In this example, the estimated distance between an observation $X_{i}$ and the cluster centroids $\{C_{k}\}_{k=1}^{4}$ are represented with a Triangular fuzzy number. As shown, $X_{i}$ is closest to the centroid of $C_{1}$ and furthest from the centroid of $C_{4}$. Even considering a possible error by $L_{mask}$, it is still very unlikely that $X_{i}$ should fall into cluster $C_{4}$, which can be formulated by the condition that $a_{4}\geq c_{1}$. To calculate how likely it is that observation should fall into $C_{1}$, $C_{2}$, $C_{3}$ and $C_{4}$, we use the overlapping areas as a membership function. The most likely cluster for $X_{i}$ is $C_{1}$ with the overlapping area $\mu_{i,1}$, as the areas for $\mu_{i,2}$ and $\mu_{i,3}$ are smaller, and $\mu_{i,4}=0$. The sum of the total overlapped area is $\mu_{i,1}+\mu_{i,2}+\mu_{i,3}+\mu_{i,4}$. Therefore, we can normalize the weight of $X_{i}$ to each cluster as $w_{C_{1}}(X_{i})=\frac{\mu_{i,1}}{\mu_{i,1}+\mu_{i,2}+\mu_{i,3}+\mu_{i,4}}$, $w_{C_{2}}(X_{i})=\frac{\mu_{i,2}}{\mu_{i,1}+\mu_{i,2}+\mu_{i,3}+\mu_{i,4}}$, $w_{C_{3}}(X_{i})=\frac{\mu_{i,3}}{\mu_{i,1}+\mu_{i,2}+\mu_{i,3}+\mu_{i,4}}$, and $w_{C_{4}}(X_{i})=0$.

The degree of membership of $\mu_{i,k}$ is calculated based on the fuzzy distance. We define the triangular fuzzy distance as a $\mathrm{TFN}(a,b,c)$, that is,

where $\hat{d}_{i,z}$ is the predicted distance between observation $X_{i}$ and the cluster $C_{z}$. Since all triangles are congruent triangles, the triangles of $\mu_{i,z}$ and $\mu_{i,1}$ must be similar. Hence, for simplicity in explaining how the degrees of membership are normalized, we have assumed that the degree of membership of an observation to the closet cluster is equal to 1, i.e., $\mu_{i,1}=1$. Based on the area theorem of similar triangles, if two triangles are similar, then the ratio of the area of both triangles is proportional to the square of the ratio of their corresponding sides. Therefore, the degree of membership of observation $i$ to its $z$-th closest cluster is

where $\left(\frac{c_{1}-a_{z}}{c_{1}-a_{1}}\right)^{2}$ is the square of the ratio of their corresponding sides. Substitute to Eq. (4) we have the triangular fuzzy-weighted frequency represented as:

Because of their smoothness and concise notation, Gaussian and bell membership functions are popular methods for specifying fuzzy sets. With these functions, the distance between the cluster centroids and the observations is a Gaussian fuzzy number, $\mathrm{GFN}(a,b)$, and the prediction value $\hat{d}$ plus standard deviation $\sigma_{mask}$ are used to model the distance, that is,

The degree of membership of $X_{i}$ to $C_{z}$ is the area under the curve and is formulated as

And, similar to the triangular fuzzy number, the Gaussian fuzzy-weighted frequency is

In contrast to distance fuzzification, the crisp frequency is measured based on $\hat{d}$, and the frequency measuring is proceeded by finding the closest cluster.

To simplify the computation and to avoid the standard deviation $\sigma_{mask}$ become too large to provide reasonable fuzzy distances. We propose a Top-Q cluster constrain for normalization the weight, that is, only the closest $Q$ clusters will be considered for frequency measuring, that is, the Eq. (4) is updated as follow.

where, $Q\leq K$, the $I_{z\leq Q}$ is the indicator function that $I_{z\leq Q}=1$ if $z\leq Q$, otherwise $0$. When $Q=1$, the measured frequency will be equal to the crisp distance frequency. The Alg 2 presents the pseudocode of fuzzy-weighted frequency measuring. With the frequency, we can apply the chi-square test to accept or reject the null hypothesis.

Alg. 2 presents the pseudocode of the fuzzy-weighted frequency. The runtime complexity is determined by $k$-means clustering (Line 1-3), imputation (Line 4-8) and measuring the fuzzy-weighted frequency (Line 9-15). For the clustering, according to Lloyd’s $k$-means algorithm, the runtime complexity of the clustering procedure in our algorithm is $\mathcal{O}(m^{2}n)$. The complexity of imputing the data is based on the worst case of the imputation algorithm in $\Xi$. Given our default settings, the worst case is an iterative imputation with $\mathcal{O}(m^{2}n)$ on a $m\times n$ data set. The fuzzy weighting process has linear runtime complexity of $\mathcal{O}(mn)$. Therefore, the runtime complexity of Alg. 2 is $\mathcal{O}(m^{2}n)$. Combining Algs. 1 and 2, the overall runtime complexity of MDL-FWF is $\mathcal{O}(n_{split}m^{2}n+m^{2}n)=\mathcal{O}(n_{split}m^{2}n)$, which is bounded by the MDL algorithm.

This experiment evaluate the extent to which masked distance learning can reduce estimation errors.

For each distribution, we generated 500 observations and randomly removed some values in different configurations, as summarized in Table I. For the uniform, Gaussian, and exponential distributions, we generated $n=10$ dimensional data. The first $5$ dimensions had a missing value rate of 20%. For the Poisson, and uniform categorical distributions, we generated $n=5$ dimensional data where the first $3$ dimensions had a missing value rate of 20%.

Missing values were replaced with five kinds of imputed values: zeros, means, medians, most frequent, and iterative. We generated five sample sets with different distributions and stored the true pairwise distance matrix for each. Then we randomly masked some values as null and applied different distance estimation methods to estimate the pairwise distance matrix. The error between the estimated distance and the true distance was used to evaluate performance.

We used both root mean squared error (RMSE) and mean absolute error (MAE) as the metric for distance estimation evaluation. RMSE gives a relatively high weight to large errors. This means the RMSE should be more useful when large errors are undesirable. MSE does not penalize huge errors as badly as RMSE, but it is a good metric to reflect the overall prediction accuracy.

The distance estimation results are shown in Table III, and the feature importance is plotted in Fig. 7. In this experiment, we find that combining multiple imputation methods and training a regression model improves the overall estimations. It is easy to understand that regression considers the different imputation method will reach the best result. Compared to calculating distance from imputed values, using masked distance learning to predict the distance had lower RMSE and MAE. This inspired us to think that combining multiple instances of uncertainty into one problem and searching for the best optimization solution may produce a better result than solving each uncertainty issue individually.

Regarding the feature importance returned by the gradient boosting decision tree, we noticed that different imputation methods resulted in different importance ranks for different distributions. For example, mean imputation contributed the most for the distance etismtaion with missing values on uniform and Poisson distributions while, for exponential distributions, zero imputation was more important. The overall estimation results show that integrating diverse imputation method for distance estimation with missing values is effectual.

In this experiment, we analyzed the contribution of each component’s performance of MDL-FWF.

First, we needed to simulate drift, so we generated data sets using three different distributions and five different types of drift. We adopt the commonly used drift synthetic methods for drift detection [7, 36, 37, 39].

Config. 1 – No missing values (Complete). The first configuration is the MDL-FWF baseline, i.e., drift detection with no missing values in the data set. In this configuration, we created k-Means space partitions and counted the observations in each cluster as the frequency. The presence of a drift was determined from a Pearson’s chi-square test.

Config. 2 – Impute missing values (MV Impute). In this setting, we duplicated the data set generated in Config. 1 and introduced MCAR missing values by randomly replacing some values with null. Those missing values were then imputed with iterative imputation. We chose iterative imputation because after comparing the five imputation methods, the iterative approach gave the best performance except MDL. k-Means and Pearson’s chi-square tests were used to detect drifts.

Config. 3 – Masked distance learning only (MDL). In this configuration, we used the masked distance learning algorithm to estimate the distance between the cluster centroids and the observations with missing values. Those estimated distances were then used to compute the contingency table. Again, Pearson’s chi-square test was used to detect drifts.

Configs. 4.1 & 4.2 – MDL-FWF. Similar to Config. 3, we predicted the distance between the cluster centroids and the observations with missing values via MDL. However, as opposed to using the estimated distances to build the contingency table, we used the fuzzy-weighted frequency measurements via a Gaussian fuzzy membership function (Config. 4.1) and a triangular fuzzy membership function (Config. 4.2). Pearson’s chi-square tests determined drift.

Configs. 5.1 & 5.2 – MDL-FWF with the top-Q clusters (MDL-FWF-TopQ). Based on Configs. 4.1 and 4.2, we measured the fuzzy-weighted frequency of the top-Q clusters, instead of considering the degree of the membership to all clusters. The same membership functions as in Configs. 4.1 and 4.2 apply to Configs. 5.1 and Config. 5.2.

The conventional drift detection evaluation metrics are Type-I and Type-II errors. Type-I errors are rejections of a true null hypothesis (also known as “false positives”), and Type II errors are the false null hypothesis rates (i.e., “false negatives”). For concept drift detection, the null hypothesis is that, there is no concept drift between the populations. Commonly, the Type-I and Type-II errors are shown separately and it is difficult to visualise them in one plot. Therefore, to better visualize the results of each drift detection method, we used the drift detection ratio as the metric [34], which is the number of rejected hypotheses divided by the total number of tests.

To quantify the performance of each configuration, we used Pearson’s correlation coefficient as an indicator of the similarity between different drift detection configurations. Since the base model is kMeans, the drift detection ratios of Config. 1 (No MV) forms the baseline and Configs. 2 to 5 show the impact of each component. A higher Pearson’s correlation coefficient to Config. 1 means that the component is less affected on the missing values.

Fig. 8 summarizes the results of Experiment 2, and Table IV records the average drift detection ratio in terms of the drift severity level. Further, the Pearson’s correlation coefficients for each configuration follows. Bold indicates the best configuration.

• •

  Config. 1 vs. 2: 0.798

• •

  Config. 1 vs. 3: 0.922

• •

  Config. 1 vs. 4.1: 0.651

• •

  Config. 1 vs. 4.2: 0.694

• •

  Config. 1 vs. 5.1: 0.924

• •

  Config. 1 vs. 5.2: 0.944

These coefficients reveal the impact of each configuration on the drift detection results. Noting that the closer the coefficient to 1, the better the missing value handling method is.

According to the coefficient, it is clear that simply applying data imputation then concept drift detection (Config. 2, MV Impute) did not produce very accurate results. Table IV shows that imputation might increase sensitivity for drifts of low severity but, with very severity drifts, it could have a opsite impact. Config. 3, where the distance between observations was estimated with masked distance learning, had a coefficient (0.922) very close to the baseline. However, the average drift detection ratios shown in Table IV tell us that MDL increased the drift sensitivity in most cases. This opens the possibility that MDL may have a high false alarm when there is no drift and, therefore, will not reach the required confidence level. By contrast, Config. 4.1 and 4.2 (MDL-FWF with Gaussian and triangular functions) performed very poorly. These configurations treat distance estimation errors as an uncertainty issue and apply a fuzzy membership function to model the weight of observations. However, in some cases, the distance estimation errors can be very large, which results in an overestimation of the observations’ weights, in turn reducing the sensitivity of drift detection. Our solution was to choose the top-Q clusters only to handle these overestimated weights. As Configs. 5.1 and 5.2 show, performance with this approach was almost as good as the baseline without increasing the false alarm rate, where the baseline was the drift detection on the complete data sets.

These experiments were designed to evaluate MDL-FWF’s performance in comparison to other concept drift detection algorithms. We chose four state-of-the-arts. The evaluation data sets and metric are the same as Experiment 2.

Multivariate Wald-Wolfowitz test (MWW test). This test was developed by Friedman and Rafsky [39] as a multivariate generalization of the Wald-Wolfowitz test. This is the baseline for this experiment. The implementation is publicly available online¹¹1https://gist.github.com/vmonaco/e9ff0ac61fcb3b1b60ba. No specific parameters are required to run this test.

QuantTree with Pearson statistics is a histogram-based change detection method for multivariate data streams [37]. Pearson’s statistic is selected because it has the best reported performance. The implementation is available online²²2http://home.deib.polimi.it/boracchi/Projects/projects.html. The partition size was set to the sample size divided by $50$, which is the same as the default setting for MDL-FWF.

MDL-FWF-TopQ Our proposed algorithm, again, implemented with two different membership functions, i.e., Gaussian and Triangle. The implementation code is available online³³3https://github.com/Anjin-Liu/TFS-MDL-FWF.

ME Test is a generic normalized mean embedding (ME) test using a specified kernel. This is used in [43, 44]. The implementation of the ME test is available online⁴⁴4https://github.com/wittawatj/interpretable-test. The algorithm parameters were set to the defaults, except for $J$, which was set to $5$ following the authors’ recommendation.

Quadratic MMD Test [45], where the null distribution is computed by permutation. The implementation is available online⁵⁵5https://github.com/wittawatj/interpretable-test. The parameters were set following the authors’ recommendation.

The drift detection results are plotted in Fig. 9, and the average ratios on each of the data sets are shown in Tables V and VI. The Pearson’s correlation coefficients of the 50 drift detection ratios were as follows:

• •

  ME missing vs. complete 0.235

• •

  MWW missing vs. complete 0.667

• •

  QuantTree missing vs. complete 0.785

• •

  MMD missing vs. complete 0.887

• •

  MDL-FWF-Gau-TopQ vs. complete: 0.924

• •

  MDL-FWF-Tri-TopQ vs. complete: 0.944

Missing values impacted the ME test the most, followed by the MWW test, QuantTree, MMD test and the MDL-FWF methods. These results confirm our assumption that missing values introduce uncertainty into distribution-based concept drift detection. Although it is difficult to examine their general performance on all distributions, these results, at least, demonstrate the need to consider missing values when dealing with distribution changes. Another conclusion we draw from these results is that fuzzy theory has good potential for handle missing values in tasks that required accurate and sensitive concept drift detection.

Our final test was to compare MLD-FWF with state-of-the-art algorithms in real-world scenarios.

HIGGS Bosons and Background Data Set[35]. The objective of this data set is to distinguish signatures of the processes that produce Higgs boson particles from background processes that do not. We selected four low-level indicators of azimuthal angular momenta for four particle jets as features as the same as [35]. The jet momenta distributions of the background processes are denoted as Back ($\mathcal{B}$), and the processes that produce Higgs bosons are denoted as Higgs ($\mathcal{A}$). There were four types of data integration: Back-Back (Type-I Real I.a), where both sample sets were drawn from Back; Higgs-Higgs (Type-I Real I.b), where both sample sets were drawn from Higgs; Back-Higgs (Type-II Real I.a); and Higgs-Back (Type-II Real I.b)

Arabic Digit Mixture Data Set [38]. This data set contains audio features of 88 people (44 females and 44 males) pronouncing Arabic digits between $0$ and $9$. We applied the same configuration as Denis et al. [38]. The data set has $26$ attributes which are a replacement mean and standard deviation for $13$ time series. Mixture distribution $\mathcal{A}$ contained randomly selected samples of both males and females, with female labels from $0$ to $4$ and male labels from $5$ to $9$. Mixture distribution $\mathcal{B}$ reversed the labels, i.e., males from $0$ to $4$ and females from $5$ to $9$. The drift detection on $\mathcal{A}-\mathcal{A}$ is Type-I Real II.a; $\mathcal{B}-\mathcal{B}$ is Type-I Real II.b; $\mathcal{A}-\mathcal{B}$ is Type-II Real II.a; and $\mathcal{B}-\mathcal{A}$ is Type-II Real II.b.

Insects Mixture Data Set [38]. This data set contains 49 features from a laser sensor. The task is to distinguish between 5 possible specimens of flying insects that pass through a laser in a controlled environment (Flies, Aedes, Tarsalis, Quinx, and Fruit). A preliminary analysis showed no drift in the feature space. However, the class distributions gradually change over time. To simulate drift in multiple clusters, we selected the samples from different insects and grouped them together at different percentages; thus, the data distribution varies. The populations are generated baBsed on the ratios that $\mathcal{A}$=$\{Flies:0.2,Aedes:0.2,Tarsalis:0.2,Quinx:0.2,Fruit:0.2\}$ and $\mathcal{B}$$=$$\{Flies:0.14,Aedes:0.14,Tarsalis:0.2,Quinx:0.2,Fruit:0.32\}$