Formation of a feature space for the anomaly detection problems in the behavior of objects with the use of data streams

The current systems to collect and process information in the economic, technical and social spheres generate large amounts of data, recording the state of observed objects at discrete points in time. A typical example is cash register equipment (CRE), transmitting fiscal data to tax authorities on a daily basis – information on revenue, number of transactions, share of cash in circulation, as well as on the range of goods and services sold [1–4]. Such data form multidimensional time series with low intensity – one measurement per day or several hours [5; 6].

One of the applied problems solved within the framework of tax administration is the identification of objects whose behavior deviates from the norm. The anomalies are specifically complex if they appear not in the form of isolated outliers, which are easily detected by statistical methods, but in the form of long-term structural changes in the dynamics of indicators [7; 8]. An example of such an anomaly is the use of cash register machines in schemes aimed to conceal real financial turnover (the so-called “rental” schemes, when the taxpayer provides their cash register machine to another person) [9; 10].

The theme of our research is the behavior of a technical object controlled by a person (subject). In our practical example, the cash register is controlled by both a salesperson, who sets the operating hours and product range at the outlet, and a customer, who makes purchasing decisions in accordance with their preferences. The process of interaction between the seller and buyers is a non-stationary temporal process, the parameters of which, under normal conditions, depend on the season, the current market situation, the effectiveness of assortment management, pricing, and logistics processes on the seller's side, that is, the pattern of their behavior, which adapts to the pattern of buyers' behavior. We assume for an anomalous object the seller's behavior pattern changes, and a specific behavioral anomaly, such as the "rental" pattern of using a cash register, is reflected in the set of values of the parameters (statistical features) of the process. In some occasions, the buyer's behavior may simply be simulated by the seller (fictitious sales for small amounts to simulate normal operations). In other cases, the location and product range of the retail outlet change, meaning both the customer base and the actual seller alter, which is also reflected in the process parameter values. Moreover, the seller can use techniques aimed at smoothing out the changes that arise, which complicates the detection of anomalous behavior using simple statistical methods. In addition to our practical example with cash registers, similar complex controlled objects generating streaming data include, for example, traffic in communication networks – mobile or fixed. In this case, the technical object is the communication line, as in our example, only records are the fact of interaction between two subjects – the information provider and its consumer (subscriber). Moreover, anomalous behavior also indicates the substitution of one of the interacting entities, for example, as a result of hacking the systems of one of them.

Traditional statistical methods to detect anomalies are Z-score, interquartile ranges, distributionbased methods [11] – have low sensitivity to such patterns, especially in highly unbalanced samples (the proportion of anomalous objects is less than 1 %) [12; 13]. Machine learning methods [14–16], on the contrary, demonstrate high efficiency, but require significant volumes of labeled data, which is a serious limitation in practice.

While solving a number of problems, for example, in problems of forecasting the development of ecosystems [17], time series with an extremely low sampling frequency (one measurement per month) and pronounced seasonality are used, due to which successive measurements of the same indicator cannot be considered as a single time series. In this occasion, each measurement of every indicator is studied as a separate informative feature, which is only acceptable at extremely low data flow rates.

This paper proposes an approach based on explicit hypothetical modeling of target anomalies. The advanced hypotheses are based on a specific practical example – data from a cash register, but without being tied to specific values of the indicators of a given object, they can be applied to other objects and tasks, such as identifying anomalies in the behavior of subscribers of communication networks [18] or some tasks of medical diagnostics [19; 20]. Each assumption concerning the behavior of an anomalous object is then formalized in the form of one or more quantitative features extracted from the original time series. The proposed approach ensures interpretability: each feature has a substantive justification, feasibility means features are not introduced “at random”, but correspond to expert logic, and a reduction in dependence on the size of the training sample – features purposefully highlight the most significant differences.

Formal problem statement

We specify a set of objects Ω = {ω 1 , …, ω N }. For each object ω i a time series is defined.

Y i ( t ) = ( y i ⁽¹⁾ ( t ), …, y i ^{( m )} ( t )) ^T ∈ R^m , t = 1, …, T ,

where T is the observation horizon (the number of discrete time intervals studied); m is the number of primary indicators recorded. It is assumed that this number is relatively small. In our practical example, the number of primary indicators is extremely limited, m = 4, namely: y_i ⁽¹⁾ ( t ) is total revenue per day; y_i ⁽²⁾ ( t ) is number of transactions (receipts) per day; y_i ⁽³⁾ ( t ) is the share of cash in revenue; y_i ⁽⁴⁾ ( t ) is average price of a good (service) per day. The object is represented by a vector of initial characteristics.

It is required to construct a decision rule F:Ω → {0,1}, where F(ωi) = 1 is interpreted as the presence of characteristics of the target anomaly in the object (in our practical example, participation in a rental scheme).

Let a limited-volume training sample L ⊂ Ω be given with true labels y_i ∈ {0,1}. The volume | L | does not exceed several dozen objects. The proportion of the positive class in the general population does not exceed 0.01. Time series y_i ⁽¹⁾ ( t )– y_i ⁽⁴⁾ ( t ) are characterized by non-stationarity (the mean and variance can change over time), the presence of gaps (lack of data for individual days), and the presence of noise (random outliers not associated with an anomaly) [21].

The decision rule F is the object of this study, and the subject of the study is the methods for forming the space of values of informative features f i,j = f j ( Y i ( t )) of objects, where i is the object number, j is the feature number. Our goal is to construct a decision rule that operates on the basis of a generated set of informative features, ensuring an increase in the efficiency of anomaly detection in accordance with known efficiency criteria.

In our occasion, the decision rule is applied not directly to the initial data Y i ( t ), but to the values of some statistical features (statistics) f i,j calculated on their basis. Each of the features is a functional (statistic) defined on a sample (time series) corresponding to the objects under consideration. The researched problem is reduced to the formation of such a set of features (statistics) f 1 ( Y i ( t )) … f m ( Y i ( t )), which would ensure the best (in practice, it is sufficient) efficiency to solve the problem of classifying objects.

Therefore, our problem can be formalized as follows:

L(F*(f1(Y1(t)), … fm(Y1(t)), …, F*(f1(YN(t)), … fm(YN(t))) → maxf1,...,fm, where F* is a decision rule operating on the basis of informative features f1(Yi(t)), … fm(Yi(t)), L(.) is the selected measure of efficiency. Practically, we have a multi-criteria problem where, by increasing the selected measure of efficiency (for example, accuracy), some acceptable values of other measures of efficiency must be ensured.

The decision rule F ^* is based on a certain method of classifying objects (an algorithm or algorithmic combination) that processes secondary data – the calculated values of the features (statistics) of a set of objects f_i, j . We assume that such a method is known; our task focuses only on the selection of a set of secondary informative features (statistics).

The current study proposes a two-stage procedure for forming a feature space. The first stage advances some rational hypotheses concerning the possible behavior of anomalous objects, on the basis of which a set (probably redundant) of features is formed, based on their values each of the hypotheses could be manifested. Therefore, we form a complete set of features f^* 1 ( Y i ( t )), … f^* M ( Y i ( t )). The second stage proposes to apply a set of methods for selecting (filters) features that form a subset { f 1 ( Y i ( t )), … f m ( Y i ( t ))} ⊂ { f^* 1 ( Y i ( t )), … f^* M ( Y i ( t ))}.

Assumptions about the manifestations of anomalous behavior of objects

Unlike automatic feature generation methods (e.g. based on deep convolutional or recurrent architectures [22]), the proposed approach is based on explicitly formulated and testable assumptions about how exactly the target anomaly manifests itself in the observed data.

Below there are five main hypotheses H1–H5, each generates a family of quantitative time series features for each of the primary indicators.

H1: Hypothesis of a change in the stationary regime (structural shift).

Contensive formulation. An anomalous object changes its average level of functioning at some point in time. If, up to a certain point in time τ, the values of the indicator y(t) fluctuated around one average value μ₁, then after τ they remain stable near another average value μ₂. This behavior may correspond to the start of using the rented device at a new point of sale (increase in revenue) or the end of use at the previous point (decrease in revenue).

Mathematical formalization. For each object and each primary indicator, the duration of the maximum series of consecutive reference periods with a value above the overall average is calculated:

L above = max{ k : ∀ j ∈ {1, …, k }, y ( t j ) >  y },

where y = -1 E T _i y (¹) •

In addition, the duration of the maximum series of days with a value below the overall average is determined:

L below = max{ k : ∀ j ∈ {1, …, k }, y ( t j ) <  y }.                          (3)

Also similar indicators are introduced, but relative to the moving average, yiilding (t )=1E T=i y (t) ,                                         (4)

they are designed to reduce the number of false alarms at the initial stage of the series.

H2: Hypothesis on the concentration of the total volume (the "burst" nature of the object's behavior).

Contensive formulation. Anomalous complex objects are characterized by behavior in which the main share of the total value of an indicator (for example, revenue) falls on a small number of time samples. This distinguishes them from facilities with stable, uniform activity. Regarding our practical example, rental devices are often used intensively for short periods (for example, during the operation of a temporary retail outlet), followed by extended periods of downtime.

Mathematical formalization. Let y (1) ≥ y (2) ≥ … ≥ y ( T ) be the values of the indicator in descending order. Then the feature “the proportion of k time samples with maximum values” is calculated as

E ti y W E ‘=1 y (j)'

In our practical example, we used values of k = 5, 10, 15, 30, 60 days, which allows to estimate the concentration on different scales.

H3: Hypothesis of a binary (two-mode) distribution structure.

Contensive formulation. The time series of an anomalous object can be approximated by a mixture of two distributions with different mathematical expectations: low (idle or inactive use mode) and high (active operation mode within the anomalous scheme). Normal objects operating in a stable mode reveal either a unimodal distribution or two closely spaced centers.

Mathematical formalization. For each object, the indicator values are clustered into two clusters using the k -means method [23–25] with k = 2. We denote:

μ1, μ2 – average values in clusters (μ1 < μ2);

n ₁, n ₂ – number of objects in clusters;

σ₁ ² , σ₂ ² – intracluster variances of the indicator.

The following signs are formed:

– absolute difference of means: Δμ = μ2 - μ1;

– relative difference of means: Δμ rel =

H 2 - H i ;

/ °2 + °2

– small cluster share: fsmall = min (n1, n 2 ) ni - n 2    ’

max (n,, n2)

– large cluster share: f large =         ^{1 2} .

n i - n 2

High values of Δμ and Δμrel at moderate cluster fractions indicate the presence of two well-separated operating modes, one of which (or both) may correspond to the target anomaly.

H 4 : Hypothesis of the presence of long periods of inactivity.

Contensive formulation. Anomalous objects often exhibit long intervals during which the indicator value is zero (or close to zero). This reflects the downtime of the facility between periods of one-off abnormal use (in our practical example – between tenant changes), periods of transportation or reconfiguration.

For normal facilities operating on a daily basis, such long pauses are not typical.

Mathematical formalization. The following is calculated to each object:

– total number of zero periods Z _total = E ■ 1 y (¹ ) = 0 ) ^, where ¹⁽ -) – an identity function equal to 1 if the condition is met and 0 if it is not met;

– maximum length of a continuous sequence of zero periods of the count Z max streak = max{ k : ∀ j ∈ {1, …, k }, y ( t j ) = 0};

– proportion of zero reference periods: f zeros = Z total / T .

H 5 : Hypothesis of increased variability (instability).

Contensive formulation. Anomalous complex objects generally reveal greater variability in their indicator values compared to typical objects operating in a predictable, rhythmic mode. This is due to the irregular nature of abnormal usage: periods of high activity are followed by complete inactivity, which generates high variance.

Mathematical formalization. Classical variation indices are used:

• -    standard deviation: a = ^ T E T 1 ( У ( t ) — У ) 2 ;

• –    variation coefficient CV = σ/ y (if y > 0);

• –    range R = max_t y ( t ) - min _t y ( t );

• –    interquartile range IQR = Q₃ - Q₁.

Additional hypotheses (auxiliary features).

In addition to the five main hypotheses, the feature space includes additional features of time series reflecting the distribution asymmetry (skewness) to identify a shift towards high or low values, excess (kurtosis) to identify the “heavy tailed” distribution, autocorrelation on lags (in the practical example, the lags are 1, 7, 14 days) to identify periodic structures, and the ratio of the maximum value to the mean to identify abnormal bursts.

Complete system of features

Based on the formulated hypotheses, a multidimensional feature space was formed (the result of the first stage of the feature space formation procedure). All features are calculated for each of the primary indicators m , that is the number of features is always a multiple of the number of primary indicators, as well as for some of their combinations. The features are divided into six groups according to the method used to determine their values.

The group of features based on means and deviations includes the most commonly used statistics for each primary feature:

• –    the arithmetic mean of the primary indicator;

• –    median of the primary indicator;

• –    standard deviation of the primary indicator;

• –    coefficient of variation of the primary indicator;

• –    asymmetry (skewness) of the primary indicator;

• –    excess (kurtosis) of the primary indicator;

• –    minimum value of the primary indicator;

• –    the maximum value of the primary indicator;

• –    the range of the primary indicator;

• –    interquartile range of the primary indicator;

• –    the ratio of the maximum of the primary indicator to the average;

• –    the ratio of the maximum of the primary indicator to the median;

• –    the proportion of time readings with a primary indicator value above average;

• –    the proportion of time readings with a primary indicator value above the median;

• –    the average of the logarithmic values of the primary indicator (for highly skewed series);

• –    median of logarithmic values of the primary indicator.

The feature group based on series of values (streaks) forms eight features for each primary feature, according to two types of comparison (with a general mean and with a running mean):

• –    the maximum length of a series of values of the primary indicator above the average;

• –    the maximum length of a series of values of the primary indicator below the average;

• –    the number of series of values of the primary indicator above the average;

• –    the number of series of the primary indicator is below average.

A group of features based on concentration (proportions of time readings with maximum values) includes features for each indicator and a certain number k of time readings. In our practical example with CRE k ∈ {5, 10, 15, 30, 60, 90}. The group includes the following features:

• –    the share of the sum of k time readings with maximum values of the indicator in the total sum;

• –    the proportion of the sum of k time readings with zero or minimum values of the indicator (to identify long downtimes).

A feature group based on a bi-cluster structure is formed according to the results of biclustering (separately for each object) of time intervals for each primary feature. The group includes the following indicators:

• –    the absolute difference between the mean values of the primary feature for two clusters;

• –    the relative difference between the average clusters of values of the primary characteristic (normalized by the variance);

• –    share of a small cluster;

• –    share of a large cluster.

The group of features based on zero values includes for each primary feature:

• –    the total number of time readings with a zero (or minimum) value of the primary indicator;

• –    the maximum length of a series of zero time readings with a zero (or minimum) value of the primary feature;

• – the proportion of time readings with zero (or minimal) value of the primary feature.

The group of features based on autocorrelation includes features calculated for each primary feature and certain values of the length of time lags, in our example –τ ∈ {1, 7, 14}:

– autocorrelation coefficient ρ(τ).

Combined method of feature selection

The original feature space has redundancy: some features may be related weakly to the target variable, while others are highly correlated with each other (multicollinearity). To improve the efficiency of subsequent classification models and reduce the risk of overfitting, feature selection methods are used.

We selected and studied experimentally widespread methods of reducing the feature space, presented in Table 1, where the rationale for using the feature and choosing its threshold values is also provided.

Preliminary experiments with the listed basic methods show their ability to reduce significantly the number of features for our practical example. Thus, based on the results of the correlation analysis, 40–45 features remain out of 96 (depending on the sample). After removing multicollinear features, the dimensionality is reduced to 25–35. PCA allows to reduce the dimensionality from 96 to 20–30 components, however, the obtained components do not have a direct interpretation, which is a disadvantage for problems requiring explanation of solutions.

Table 1

Method	Operating principle	Justification for the choice of threshold values
Correlation analysis with the target variable	For each feature fj , the Pearson linear correlation coefficient r ( fj , y ) with the binary target variable y is calculated. Features with | r | below a given threshold are considered weakly related to the anomaly and are excluded	In problems with highly unbalanced classes, even a moderate correlation (|r| ≈ 0.15–0.20) can be statistically significant. The threshold θ = 0.2 was chosen empirically
Selection based on the Kolmogorov– Smirnov (K-S) test	The K-S test is a nonparametric method to test the hypothesis that two samples (in this case, feature values for anomalous and normal objects) are drawn from the same distribution. The statistic is used Dn ,m = supx|Fn(x) - Gm(x)|, where Fn и Gm – empirical distribution functions	Signs with p-value ≥ 0.05 (the null hypothesis of equality of distributions is not rejected) are excluded
ANOVA	The ANOVA method (selection is based on singlefactor analysis of variance) tests the hypothesis of equality of the mean values of a feature in two groups (anomalous and normal objects). The F-statistic is introduced: F = S between / S within , where S between – inter-group dispersion; S within – intragroup variance of a feature. It is efficient for features with a nearly normal distribution. For highly skewed features, a logarithmic transformation is used	Signs with p-value ≥ 0.01 are excluded (it is a more stringent threshold compared to the K-S test to compensate for multiple comparisons)
Removal of multicollinear features	The matrix of Pearson pairwise correlations between all remaining features is calculated. From each pair of features with a correlation coefficient r > θcorr, the one with a lower correlation with the target variable is removed. Multicollinearity worsens the stability of models (especially linear ones) and complicates interpretation	The threshold θcorr = 0.8 is chosen as standard in feature selection problems
Principal component analysis (PCA)	To compare, an alternative approach is used: dimensionality reduction without removing individual features. Standardization of features is performed (normalization to zero mean and unit variance), then PCA is applied	No threshold is used. The number of components explaining at least 99% of the total variance is selected

Basic methods of feature selection

As shown above, each of the considered feature selection methods has both advantages and disadvantages. Correlation analysis effectively filters out features that do not relate linearly to the target variable, but it is insensitive to nonlinear dependencies. The Kolmogorov–Smirnov test captures any differences in distributions, but it can produce false positives with small samples. ANOVA is sensitive to differences in means but it requires normal distributions [26]. Removing multicollinear features reduces redundancy but it does not consider the relationship with the target variable.

The proposed combined method consistently applies all four methods, which allows to use the advantages of each of them and obtain a minimum set of features with maximum predictive power, while ensuring robustness to outliers and anomalies in the training set.

Step 1 forms the initial (complete) feature space. For each object ω_i ∈ Ω, all features within the complete system of features are calculated. The X ∈ R^{N × m} matrix and a vector of target labels y ∈ {0,1} ^N are formed for the marked part of the sample.

Step 2 filters the features based on their correlation with the target variable. For each feature f j , j = 1, …, m , the Pearson correlation coefficient is calculated:

r j =

Xio (xy - xj)(y - y) io J X L1( y - y )2 ’

where L is a set of marked objects. Features with | r j | ≤ 0,2 are excluded. The threshold of 0.2 was chosen empirically as ensuring the removal of features with a statistically insignificant linear relationship while simultaneously preserving potentially useful features with a nonlinear relationship.

Step 3 performs filtering using the Kolmogorov – Smirnov test. A two-sample K-S test is performed to each remaining feature f j .

The empirical distribution functions are calculated for a group of anomalous F j ^anom ( x ) and normal objects F j ^norm ( x ):

F j “°"( x ) =    ¹ X     1( X j <  x ), F j — ( x ) = X     1( X j <  x ),           (7)

i : ^y =1                                        i—H : y i =1

anom                                    norm then the statistics are evaluated

D j = x sup | F j ^anom ( x ) - F j ^norm ( x )|.

The p -value (level of significance) is calculated. Features with p ≥ 0.05 are excluded, as the distributions of the feature in the two groups are statistically indistinguishable.

Step 4 performs filtering based on ANOVA. For each remaining feature, single-factor analysis of variance is performed. The between-group sum of squares is calculated SSB = n _anom( x _anom - x ) ² + n norm ⁽ X „orm ^- X ^{) 2} , within-group sum Of squares SSW = X i _{: y/ = 1} ( ^xi ^{- x}_anOm ) ² ■ X i _{: y}i = 1 ( ^xi - ^Xn_Orm ) ² and F-statistics:

F = SSW / ( N - k )

SSB / (k -1) , where k = 2 is the number of groups. Features with p ≥ 0.01 are excluded. A more stringent threshold (0.01 instead of 0.05) was chosen to compensate of multiple comparisons.

Finally, multicollinear features are removed, for which the Pearson pairwise correlation matrix R is calculated for the remaining features. For each pair of features ( f_a , f_b ) with | r_ab | > 0,8, the feature that has a smaller absolute value of correlation with the target variable (see step 2) or, if it is equal, a smaller variance is removed. The process is repeated iteratively, since removing one feature can change the correlations between the remaining ones.

The result of the method is a set of selected features F _selected ⊂ { f ₁, …, f_m }, | F _selected| ≪ m and the corresponding matrix X selected ∈ R^N ×| F selected |.

We can estimate the computational complexity of the proposed method. Let N labeled be the number of labeled objects, of which n anom are anomalous objects and n norm are objects without signs of anomalousness, M = 96 is the number of primary features. We are able to estimate the computational complexity of the main steps.

Step 2 (correlation): O ( M ⋅ N _labeled) – linear complexity in the number of features and objects.

Step 3 (K-S test): For each feature, sorting the values into two groups is required, which results in O ( M ⋅ ( n _anom log n _anom + n _norm log n _norm)). In the worst case, n _anom + n _norm = N l _a b _e l _e d .

Step 4 (ANOVA): O ( M ⋅ N labeled ) – calculation of means and sums of squares.

Step 5 (removing multicollinearity): in the worst case, this may require O ( k ³ ) operations, where k is the number of features remaining after steps 2–4. In practice, the method converges in 2–3 iterations.

The overall complexity is O ( M ⋅ N _labeled), which is acceptable for N _labeled ≈ 2000 and M = 96, as in our practical example. The execution time using a typical modern personal computer requires a few seconds.

Experimental results

The current research uses the following machine learning models to classify objects (detect anomalies): Logistic Regression [27], Decision Tree [28], naive Bayes classifier (GaussianNB) [29] – a probabilistic model that assumes normal distributions of features, k-nearest neighbors (KNN) [30–32], support vector machine (SVC) [33], Random Forest [34], XGBoost [35] – gradient boosting optimized for speed and accuracy, LightGBM [36] – gradient boosting with leaf growth of trees.

To evaluate the comparative effectiveness of the methods used with different variants of the feature space, standard metrics to assess the quality of binary classification are applied:

– Accuracy: ( TP + TN ) / ( TP + TN + FP + FN ), where TP is a number of correctly assigned positive classifications, TN is a number of correctly assigned negative classifications, FP is number of false positive classifications, FN is a number of false negative classifications;

– Precision (prediction accuracy of the positive class): TP / ( TP + FP );

– Recall (completeness): ( TP + FN ) / TP ;

– F1-score: (Precision ⋅ Recall) / (Precision+Recall);

– ROC-AUC – area under the ROC curve (True Positive Rate vs False Positive Rate).

Under the circumstances of strong class imbalance, the focus is on F1-score and ROC-AUC, since Accuracy can be high when trivially predicting all objects as normal.

The experiments were conducted using fiscal data provided by the Federal Tax Service of Russia. The sample includes data from 2,200 cash registers, within a time horizon of 365 days. Expert marking: 40 devices were identified as participants in rental schemes (anomalous, that is positive class). Additionally, there are 9 devices with documented anomalous behavior used for validation.

Table 2 shows the classification results (Random Forest, 100 Trees) for different feature selection options.

Table 2

Comparing the feature selection methods

Selection method	Number of features	Accuracy	Precision	Recall	F1-score	ROC-AUC
Without selection (all 96)	96	0.82	0.76	0.68	0.72	0.85
Correlation with target (θ = 0.2)	42	0.87	0.82	0.74	0.78	0.90
K-S test ( p < 0.05)	38	0.86	0.81	0.73	0.77	0.89
ANOVA ( p < 0.01)	35	0.86	0.80	0.72	0.76	0.88
Removing multicollinearity (8θ = 0.8)	58	0.84	0.78	0.70	0.74	0.87
Combined (proposed)	28	0.91	0.87	0.82	0.84	0.94
PCA (99 % of variance)	24	0.89	0.84	0.79	0.81	0.92

The combined method, including the sequential application of correlation filtering, K-S test, ANOVA and removal of multicollinear features, demonstrated the best results for all metrics. Reducing the dimensionality from 96 to 28 features is accompanied by an increase in ROC-AUC from 0,85 to 0,94.

Table 3 shows a comparison of the proposed method with other combined feature selection approaches.

Table 3

Comparing the combined feature selection methods

Method	Composition	Number of selected features	ROC-AUC	Interpretability
Filter + cover (RFECV)	Correlation + recursive removal	32	0.92	average
Hybrid (mutual info + SFS)	Mutual information + sequential selection	35	0.91	low
LASSO (regularization L1)	L1-regularization	41	0.89	high
Proposed combined	Correlation + K-S + + ANOVA + multicollinearity removal	28	0.94	high

The proposed method demonstrates better classification quality with a smaller number of selected features, while maintaining high interpretability.

To evaluate the robustness of the proposed method to variations in the training sample, an experiment with bootstrapping (100 iterations) was conducted. For each iteration, 80% of the labeled objects were randomly selected with replacement, feature selection was performed, and the classifier (Random Forest) was trained. Results:

– average number of selected features: 28,3±3,1;

– selection stability (percentage of features included in selection in more than 90 % of iterations): 76 %;

– average ROC-AUC: 0,935±0,008.

The high stability of the method confirms its suitability for practical use.

Table 4 presents the results of training eight models on an expanded sample (after applying the label propagation method [37; 38], which expanded the labeling to 2200 objects). Combined feature selection was used (28 features).

Table 4

Results of object classification using different models

Model	Accuracy	Precision	Recall	F1-score	ROC-AUC
Logistic regression	0.89	0.85	0.78	0.81	0.91
Decision tree (max_depth = 5)	0.84	0.79	0.71	0.75	0.82
Naive Bayes (GaussianNB)	0.81	0.74	0.65	0.69	0.78
K-nearest neighbors (k = 7)	0.83	0.77	0.69	0.73	0.80
Support Vector Machine (RBF)	0.87	0.83	0.75	0.79	0.88
Random Forest (100 trees)	0.91	0.87	0.82	0.84	0.94
XGBoost	0.92	0.88	0.83	0.85	0.94
LightGBM	0.91	0.86	0.81	0.83	0.93

Ensemble methods (Random Forest, XGBoost, LightGBM) outperform linear and simple nonlinear models across all metrics due to their ability to model complex nonlinear dependencies between features and robustness to noise.

Table 5 provides a comparison of the proposed approach with some alternative anomaly detection methods on the same data.

The proposed approach significantly outperforms alternatives in terms of the proportion of confirmed anomalies and ROC-AUC, while maintaining high interpretability due to the explicit generation of features based on hypotheses.

Table 5

Comparison with alternative approaches

Method	Requires labeling	Proportion of confirmed anomalies	ROC-AUC	Interpretability (link to hypotheses H1–H5)
Z-score (statistical)	No	< 30 %	0.65	High
IQR (interquartile range)	No	< 30 %	0.62	High
K-means + target cluster	Partially	40 %	0.75	Average
Supervised learning (original labeling)	Yes (40 objects)	–	0.72	High
The proposed hybrid approach	Yes (40 objects + Label Spreading)	84 %	0.94	High

To confirm the practical feasibility of the proposed approach, 785 objects classified as anomalous (Random Forest, combined feature selection) were submitted to the Federal Tax Service of Russia for expert review. The results are in Table 6.

Table 6

Results of validation on independent samples

Category of objects	Number of verified	Confirmed anomalies	Confirmation rate, %
Objects inside known risk zones	915	915 (known previously)	100
Objects outside known risk zones	145	122	84

The rate of confirmed anomalies in previously unknown zones was 84 %, significantly exceeding the results achieved using clustering alone (40 % confirmation). The false-positive rate was 16 %.

The most significant features were identified due to the analysis of feature importance in the Random Forest model (the average decrease in accuracy when a feature is excluded). The results are in Table 7.

Table 7

The most informative features

Rank	Description	Accuracy decrease
1	Small cluster revenue share	0.087
2	Relative difference of average clusters	0.076
3	Maximum series length below average	0.068
4	Share of 15 maximum days in revenue	0.059
5	Number of zero days	0.052
6	Maximum series length above moving average	0.048
7	Absolute difference of average clusters	0.044
8	Coefficient of variation of the number of receits	0.041

These features directly correspond to the hypotheses of binary structure (ranks 1, 2, 7), regime change (rank 3, 6), revenue concentration (rank 4), and long-term downtime (rank 5).

Conclusion

The obtained results confirm the efficiency of the proposed methodology. The combined feature selection not only reduced the dimensionality of the space by 3.4 times, but also improved classification quality (ROC-AUC increased from 0.85 to 0.94), which is explained by the removal of noisy and redundant features that reduced the generalization ability of the models. The high proportion of confirmed anomalies in validation (84%) indicates that the proposed feature system actually captures the essential characteristics of the target anomalous behavior, rather than random correlations in the training set.

Therefore, the research has solved the problem of forming a feature space to analyse the low-intensity streaming data in order to identify anomalies in the behavior of objects using CRE as an example. A system of features from six groups has been developed, based on five clearly formulated hypotheses about the nature of anomalies (change of regime, concentration of indicators, binary (bicluster) structure, long pauses, increased variability). Mathematical formulations of each type of feature and a detailed description of all six groups of features have been provided. A combined method of feature selection has been proposed, including the sequential application of correlation analysis with the target variable (θ = 0.2), the Kolmogorov–Smirnov test (p < 0.05), ANOVA (p < 0.01) and the removal of multicollinear features (r > 0.8). For our practical example, the method allowed to reduce the dimensionality from 96 to 28 while simultaneously increasing the ROC-AUC from 0.85 to 0.94. The most informative features have been determined: the proportion of a small cluster, the relative difference between average clusters, the maximum length of a series below the average value, the proportion of time samples with maximum values of indicators, and the number of time samples with zero indicators.

The proposed approach has got a number of limitations, in particular:

• –    sensitivity to the choice of thresholds – when binarizing features (used in clustering) and selecting features, the results depend on the chosen thresholds;

• –    fixed time horizon – the approach involves analyzing data for the entire period with a fixed number within the time frame.

The proposed methodology can be applied to other problems of detecting anomalies in technical and economic systems with low-intensity data flows (equipment monitoring, financial transaction analysis, telemetry data control, detection of abnormal operating modes of technical systems, medical diagnostics). Moreover, the following areas of further research seem promising:

• –    automatic expansion of the hypothesis system using methods of symbolic regression and genetic programming [39] to identify non-obvious patterns in the data;

• –    considering conceptual drift [40] – the behavior of anomalous objects can change over time; adaptation of models in a streaming mode with periodic retraining is required;

• –    extraction of text features from the product nomenclature using natural language processing methods (word2vec, FastText, BERT) to identify discrepancies with the declared type of economic activity;

• –    graph methods to take into account the relationships between objects (migration of objects between retail outlets, common tenants) can improve the accuracy of detection through the use of structural information;

• –    the use of recurrent neural networks (LSTM, GRU) and transformers [41] for automatic feature extraction from time series, which can complement the proposed system of hypotheses.