Boosting the performance of anomalous diffusion classifiers with the proper choice of features

The most popular method of deducing the particles' type of motion from their trajectories is based on the MSD [3, 5, 6, 66], defined as

$$\begin{equation}\text{MSD}(t)=E\left({\Vert}{X}_{t+{t}_{0}}-{X}_{{t}_{0}}{{\Vert}}^{2}\right),\end{equation} \tag{ 1 }$$

where ${({X}_{t})}_{t > 0}$ is a particle's trajectory, ||⋅|| stands for the Euclidean distance and E is the ensemble average. Since in many experiments only a limited number of trajectories is observed, the time-averaged MSD (TAMSD) calculated from a single trajectory is usually used as the estimator of MSD,

$$\begin{equation}\hat{\text{MSD}}(n{\Delta}t)=\frac{1}{N-n+1}\sum\limits _{i=0}^{N-n}{\Vert}{X}_{{t}_{i+n}}-{X}_{{t}_{i}}{{\Vert}}^{2}.\end{equation} \tag{ 2 }$$

The trajectory consists of N consecutive two dimensional positions X_i = (x_i , y_i ) (i = 0, ..., N) recorded with a constant time interval Δt. The quantity n is the time lag between the initial and the final positions of the particle. If the underlying process is ergodic and has stationary increments, TAMSD converges to the theoretical MSD [67].

AnDi refers to a broad class of processes that deviate from the standard Brownian motion. Those deviations display an asymptotic power-law dependence of the MSD on the time lag,

$$\begin{equation}\hat{\text{MSD}}(n{\Delta}t)\sim {K}_{\alpha }{t}^{\alpha },\end{equation} \tag{ 3 }$$

where K_α is the generalized diffusion coefficient and α is the anomalous exponent. As already mentioned in the introduction, the value of the latter may be used to discriminate the trajectories into normal diffusion (α = 1), subdiffusion (α < 1), and superdiffusion (α > 1).

Particles' movements that exhibit deviations from the linear behavior of MSD may be described by multiple models, depending on some specific properties of the corresponding trajectories. Five of such processes have been identified within the AnDi challenge as crucial for the interpretation of the experimental data [58, 59]: CTRW, FBM, LW, ATTM and SBM (see figure 1 for example trajectories). In this section, the main characteristics of those models are briefly summarized.

Zoom In Zoom Out Reset image size

2.2.1. Continuous time random walk

2.2.2. Fractional Brownian motion

FBM [59] is the solution of the stochastic differential equation

$$\begin{equation}\mathrm{d}{X}_{t}^{i}=\sigma \mathrm{d}{B}_{t}^{H,i},\quad i=1,2.\end{equation} \tag{ 4 }$$

Here, σ > 0 is the scale coefficient, which relates to the diffusion coefficient D via $\sigma =\sqrt{2D}$. H ∈ (0, 1) is the Hurst parameter and ${B}_{t}^{H}$ is a continuous-time, zero-mean Gaussian process starting at zero, with the following covariance function

$$\begin{equation}E\left({B}_{t}^{H}{B}_{s}^{H}\right)=\frac{1}{2}\left(\vert t{\vert }^{2H}+\vert s{\vert }^{2H}-\vert t-s{\vert }^{2H}\right).\end{equation} \tag{ 5 }$$

The value of H determines the type of diffusion in the process. For $H< \frac{1}{2}$, FBM produces subdiffusion (regime with a negatively correlated noise). For $H > \frac{1}{2}$, FBM generates superdiffusive motion (positively correlated noise). It reduces to the normal diffusion at $H=\frac{1}{2}$.

2.2.3. Lévy walk

2.2.4. Annealed transient time motion

2.2.5. Scaled Brownian motion

The SBM is a process described by the Langevin equation with a time-dependent diffusivity K(t),

$$\begin{equation}\frac{\mathrm{d}X(t)}{\mathrm{d}t}=\sqrt{2K(t)}\xi (t),\end{equation} \tag{ 6 }$$

where ξ(t) is white Gaussian noise. If the diffusivity has the form K(t) = αK_α t^α−1, the MSD follows the power-law MSD(t) ∝ K_α t^α .

In machine learning, all classification algorithms may be divided into two classes. The traditional machine learning is a set of statistical learning methods, which do not operate on raw data. Instead, each data sample is characterized by a vector of human-engineered features or attributes. Those vectors are then used as input for the classifier. The second class consists of deep learning methods that identify and extract features on their own. The representation of data is constructed automatically and there is no need for its complex preprocessing as in the case of the feature-based methods.

Deep learning is nowadays the state-of-the-art technology for data classification in many domains and overshadows a little bit the qualities of the classical machine learning. However, the choice of a suitable classification method is usually more subtle than simply looking at its performance. A lot of deep learning methods do indeed an excellent job with respect to the predictions but are extremely complex to interpret. To give an example, the ResNet18 architecture we considered in one of our previous attempts to trajectory classification [57] originally had 11, 220, 420 parameters. We were able to reduce their number down to 399, 556 with a positive impact on accuracy. Although it was an impressive achievement, interpretation of all of those remaining parameters is almost impossible. Thus, if one wants to comprehend the decisions made by an ML algorithm, it could be tempted to go for feature-based methods and gain on interpretability at the expense of accuracy.

Being aware of this trade-off, we deliberately decided to follow the interpretability path in our contribution to the AnDi challenge and chose the extreme gradient boosting (XGB) algorithm with decision trees as base learners [68] as the classification method. We already used similar approach [49, 51, 52], however for different sets of data.

XGB is an ensemble method, which combines multiple decision trees to obtain better predictive performance. A single decision tree [69] is fairly simple to build. The original dataset is split into smaller subsets based on the values of a given feature. The process is recursively repeated until the resulting subsets are homogeneous (all samples from the same class) or further splitting does not improve the classification performance. A splitting feature for each step is chosen according to similarity score measure [70].

Single decision trees are famous for their ease of interpretation. However, they are sensitive to even small variations of data and prone to overfitting. That is why one rather uses forests (ensembles) of trees as reliable classifiers. In the case of the XGB algorithm, the trees are built in a stage-wise fashion. We start with a single tree trained on random subsets of data and features and then at every step, a new tree is added to the ensemble. This new tree learns from the mistakes committed by the previous trees.

The main challenge of feature-based classification methods is the proper choice of features that can differentiate between samples belonging to different classes. The impact of their choice on the final performance of the classifier was already discussed in reference [52]. In this work, we continue the search for a robust set of features for the classification of different diffusion modes.

In this section, we briefly introduce features that have been used for the AnDi challenge. Since we were not satisfied with the performance of the resulting classifier in task 2 [59], after the end of the challenge we further searched for features that would improve the classification results. Those additional features will be described in section 2.5. Both the original feature set for the challenge and its extension are listed in table 1.

2.4.1. Anomalous exponent

2.4.2. Diffusion coefficient

2.4.3. Asymmetry

The asymmetry of a trajectory can detect directed motion. According to Saxton [32], it can be derived from the gyration tensor, which describes the second moments of positions of a particle. For a 2D random walk of N steps, the tensor is given by

$$\begin{equation}\mathbf{T}=\left(\begin{matrix}\hfill \frac{1}{N}{\sum }_{j=1}^{N}{({x}_{j}-\langle x\rangle )}^{2}\hfill & \hfill \frac{1}{N}{\sum }_{j=1}^{N}({x}_{j}-\langle x\rangle )({y}_{j}-\langle y\rangle )\hfill \\ \hfill \frac{1}{N}{\sum }_{j=1}^{N}({x}_{j}-\langle x\rangle )({y}_{j}-\langle y\rangle )\hfill & \hfill \frac{1}{N}{\sum }_{j=1}^{N}{({y}_{j}-\langle y\rangle )}^{2}\hfill \end{matrix}\right),\end{equation} \tag{ 10 }$$

where $\langle x\rangle =(1/N){\sum }_{j=1}^{N}{x}_{j}$ is the average of x coordinates over all steps in the random walk. We define the asymmetry as [72]

$$\begin{equation}A=-\mathrm{log}\left(1-\frac{{({\lambda }_{1}-{\lambda }_{2})}^{2}}{2({\lambda }_{1}+{\lambda }_{2})}\right),\end{equation} \tag{ 11 }$$

where λ₁ and λ₂ are the principle radii of gyration, i.e. the eigenvalues of the tensor T.

2.4.4. Efficiency

Efficiency E measures the linearity of a trajectory. It relates the net squared displacement of a particle to the sum of squared step lengths,

$$\begin{equation}E=\frac{\vert {X}_{N-1}-{X}_{0}{\vert }^{2}}{(N-1){\sum }_{i=1}^{N-1}\vert {X}_{i}-{X}_{i-1}{\vert }^{2}}.\end{equation} \tag{ 12 }$$

It may help to detect directed motion (superdiffusion).

2.4.5. Empirical velocity autocorrelation function

Empirical velocity autocorrelation function [73] for lag 1 and point n is defined as:

$$\begin{equation}{\chi }_{n}=\frac{1}{N-1}\sum\limits _{i=0}^{N-2}({X}_{i+1+n}-{X}_{i+n})({X}_{i+1}-{X}_{i}).\end{equation} \tag{ 13 }$$

It can be used to distinguish subdiffusion processes. In our model, we used χ_n for points n = 1 and n = 2.

2.4.6. Fractal dimension

Fractal dimension is a measure of the space-filling capacity of a pattern (a trajectory in our case). According to Katz and George [74], the fractal dimension of a planar curve may be calculated as

$$\begin{equation}{D}_{f}=\frac{\mathrm{ln}\,N}{\mathrm{ln}(Nd{L}^{-1})},\end{equation} \tag{ 14 }$$

where L is the total length of the trajectory, N is the number of steps, and d is the largest distance between any two positions. It takes values around 1 for straight trajectories (i.e. directed motion), around 2 for random ones (normal diffusion), and around 3 for constrained trajectories (subdiffusion).

2.4.7. Maximal excursion

Maximal excursion of the particle is given by the formula:

$$\begin{equation}\text{ME}=\frac{\mathrm{m}\mathrm{a}\mathrm{x}({X}_{i+1}-{X}_{i})}{{X}_{N-1}-{X}_{0}}.\end{equation} \tag{ 15 }$$

It should detect relatively long jumps (in comparison to the overall displacement).

2.4.8. Mean maximal excursion

According to reference [37], the mean maximal excursion is usually a better observable than MSD to determine the anomalous exponent α. Given the largest distance traveled by a particle,

$$\begin{equation}{D}_{N}=\mathrm{max}\vert {X}_{i}-{X}_{0}\vert ,\end{equation} \tag{ 16 }$$

the mean maximal excursion is defined as its standardized value, i.e.:

$$\begin{equation}{T}_{n}=\frac{\mathrm{max}(\vert {X}_{i}-{X}_{0}\vert )}{\sqrt{{\hat{\sigma }}_{N}^{2}({t}_{N}-{t}_{0})}}.\end{equation} \tag{ 17 }$$

Here, ${\hat{\sigma }}_{N}$ is a consistent estimator of the standard deviation of D_N ,

$$\begin{equation}{\hat{\sigma }}_{N}^{2}=\frac{1}{2N\delta t}\sum\limits _{j=1}^{N}{\Vert}{X}_{j}-{X}_{j-1}{{\Vert}}_{2}^{2}.\end{equation} \tag{ 18 }$$

2.4.9. Mean gaussianity

Gaussianity g(n) [75] checks the Gaussian statistics of increments of a trajectory and is defined as

$$\begin{equation}g(n)=\frac{2\langle {r}_{n}^{4}\rangle }{3{\langle {r}_{n}^{2}\rangle }^{2}},\end{equation} \tag{ 19 }$$

where $\langle {r}_{n}^{k}\rangle$ denotes the kth moment of the trajectory at time lag n:

$$\begin{equation}\langle {r}_{n}^{k}\rangle =\frac{1}{N-n}\sum\limits _{i=1}^{N-n}\vert {X}_{i+n}-{X}_{i}{\vert }^{k}.\end{equation} \tag{ 20 }$$

Gaussianity for normal diffusion is equal to 0. The same result should be obtained for FBM, since its increments follow Gaussian distribution. Other types of motion should show deviations from zero.

Instead of looking at gaussianities at single time lags, we will include the mean over all lags as one of the features:

$$\begin{equation}\langle g\rangle =\frac{1}{N}\sum\limits _{i=1}^{N}g(n).\end{equation} \tag{ 21 }$$

2.4.10. Mean-squared displacement ratio

MSD ratio gives information about the shape of the corresponding MSD curve. We will define it as

$$\begin{equation}\text{MSDR}({n}_{1},{n}_{2})=\frac{\langle {r}_{{n}_{1}}^{2}\rangle }{\langle {r}_{{n}_{2}}^{2}\rangle }-\frac{{n}_{1}}{{n}_{2}},\end{equation} \tag{ 22 }$$

where n₁ < n₂. MSDR = 0 is zero for normal diffusion (α = 1). We should get MSDR ⩽ 0 for sub- and MSDR ⩾ 0 for superdiffusion. We simply took n₂ = n₁ + Δt and calculate an averaged ratio for every trajectory.

2.4.11. Kurtosis

Kurtosis gives insight into the asymmetry and peakedness of the distribution of points within a trajectory [72]. To calculate it, the position vectors X_i are projected onto the dominant eigenvector $\vec{r}$ of the gyration tensor (10),

$$\begin{equation}{x}_{i}^{p}={X}_{i}\cdot \vec{r}.\end{equation} \tag{ 23 }$$

Kurtosis is then defined as the fourth moment of ${x}_{i}^{p}$,

$$\begin{equation}K=\frac{1}{N}\sum\limits _{i=1}^{N}\frac{{({x}_{i}^{p}-{\bar{x}}^{p})}^{4}}{{\sigma }_{{x}^{p}}^{4}},\end{equation} \tag{ 24 }$$

with ${\bar{x}}^{p}$ being the mean projected position and ${\sigma }_{{x}^{p}}$—the standard deviation of x^p .

2.4.12. Statistics based on p-variation

The empirical p-variation is given by the formula [76, 77]:

$$\begin{equation}{V}_{m}^{(p)}=\sum\limits _{i=0}^{\frac{N}{m}-1}\vert {X}_{(i+1)m}-{X}_{im}{\vert }^{p}.\end{equation} \tag{ 25 }$$

These statistics can be used to detect the fractional Lévy stable motion (including FBM). We defined 6 features based on ${V}_{m}^{(p)}$:

• power γ^p fitted to p-variation for lags 1 to 5 [76],
• statistics P used in reference [52], based on the monotonicity changes of ${V}_{m}^{(p)}$ as a function of m:
  $$\begin{equation}P=\begin{cases}0\hfill & \text{if}\ {V}_{m}^{(p)}\ \text{does}\;\text{not}\;\text{change}\;\text{the}\;\text{monotonicity},\hfill \\ 1\hfill & \text{if}\;\ {V}_{m}^{(p)}\ \text{is}\;\text{convex}\;\text{for}\;\text{the}\;\text{highest}\ p\ \text{for}\;\text{which}\;\text{it}\;\text{is}\;\text{not}\;\text{monotonous},\hfill \\ -1\hfill & \text{if}\;\ {V}_{m}^{(p)}\ \text{is}\;\text{concave}\;\text{for}\;\text{the}\;\text{highest}\ p\ \text{for}\;\text{which}\;\text{it}\;\text{is}\;\text{not}\;\text{monotonous}.\hfill \end{cases}\end{equation} \tag{ 26 }$$

2.4.13. Straightness

Straightness S measures the average direction change between subsequent steps. It relates the net displacement of a particle to the sum of all step lengths,

$$\begin{equation}S=\frac{\vert {X}_{N-1}-{X}_{0}\vert }{{\sum }_{i=1}^{N-1}\vert {X}_{i}-{X}_{i-1}\vert }.\end{equation} \tag{ 27 }$$

2.4.14. Trappedness

With trappedness we will refer to the probability that a diffusing particle is trapped in a bounded region with radius r₀. A comparison of analytical and Monte Carlo results for confined diffusion allowed Saxton [32] to estimate this probability with

$$\begin{equation}P(D,t,{r}_{0})=1-\mathrm{exp}\left(0.2048-0.251\,17\left(\frac{Dt}{{r}_{0}^{2}}\right)\right).\end{equation} \tag{ 28 }$$

Here, r₀ is approximated by half of the maximum distance between any two positions along a given trajectory. D in equation (28) is estimated by fitting the first two points of the MSD curve (i.e. it is the so called short-time diffusion coefficient).

In order to improve the performance of the original classifier, we searched for further feature candidates after the AnDi challenge. The additional features extending the basic set are described below.

2.5.1. D'Agostino-Pearson test statistic

D'Agostino–Pearson κ² test statistic [78, 79] is a goodness-of-fit measure that aims to establish whether or not a given sample comes from a normally distributed sample. It is defined as

$$\begin{equation}{\kappa }^{2}={Z}_{1}({g}_{1})+{Z}_{2}(K),\end{equation} \tag{ 29 }$$

where K is the sample kurtosis given by equation (24) and ${g}_{1}={m}_{3}/{m}_{2}^{3/2}$ is the sample skewness with m_j being the jth sample central moment. The transformations Z₁ and Z₂ should bring the distributions of the skewness and kurtosis as close to the standard normal as possible. Their definitions may be found elsewhere [78, 79]. This feature should help to distinguish ATTM and SBM from other motions.

2.5.2. Kolmogorov–Smirnov statistic against χ² distribution

The Kolmogorov–Smirnov statistic quantifies the distance between the empirical distribution function of the sample F_n (X) and the cumulative distribution function G_n (X) of a reference distribution,

$$\begin{equation}{D}_{n}=\underset{X}{\mathrm{sup}}\vert {F}_{n}(X)-{G}_{n}(X)\vert .\end{equation} \tag{ 30 }$$

Here, n is the number of observations (i.e. the length of a trajectory). The value of this statistic for the empirical distribution of squared increments of a trajectory against the sampled χ² distribution has been taken as the next feature. The rationale of such choice is that for a Gaussian trajectory the theoretical distribution of squared increments is the mentioned χ² distribution.

2.5.3. Noah, Moses, and Joseph exponents

2.5.4. Detrending moving average

The detrending moving average (DMA) statistic [81, 82] is given by

$$\begin{equation}\text{DMA}(\tau )=\frac{1}{N-\tau +1}\sum\limits _{i=\tau }^{N}{\left({X}_{i}-{\bar{X}}_{i}^{\tau }\right)}^{2},\end{equation} \tag{ 31 }$$

for τ = 1, 2, ..., where ${\bar{X}}_{i}^{\tau }$ is a moving average of τ observations, i.e. ${\bar{X}}_{i}^{\tau }=\frac{1}{\tau +1}{\sum }_{j=0}^{\tau }{X}_{i-j}$.

As mentioned in reference [82], the DMA-based statistical tests can help in the detection of the SBM. In our model, we used two values of DMA for each trajectory as input features, namely DMA(1) and DMA(2).

2.5.5. Average moving window characteristics

As the moving average methods have already been successfully applied to many problems (see for example reference [83] for change point detection), we decided to accommodate them in our feature set as well. We added eight features based on the formula

$$\begin{equation}\text{MW}=\frac{1}{2(N+1)}\sum\limits _{t=0}^{N}\left\vert \mathrm{sgn}\left({\bar{X}}_{i+1}^{(m)}-{\bar{X}}_{i}^{(m)}\right)-\mathrm{sgn}\left({\bar{X}}_{i+1}^{(m)}-{\bar{X}}_{i}^{(m)}\right)\right\vert ,\end{equation} \tag{ 32 }$$

where ${\bar{X}}^{(m)}$ denotes a statistic of the process calculated within the window of length m and sgn is the signum function. In particular, we used the mean and the standard deviation for $\bar{X}$ and calculated moving window (MW) with windows of lengths m = 10 and m = 20 separately for x and y coordinates.

2.5.6. Maximum standard deviation

The idea of the MW helped us to introduce another two features based on the standard deviation σ_m of the process calculated within a window of length m. They are given by

$$\begin{equation}\text{MXM}=\frac{\mathrm{max}\left({\sigma }_{m}(t)\right)}{\mathrm{min}\left({\sigma }_{m}(t)\right)}\end{equation} \tag{ 33 }$$

and

$$\begin{equation}\text{MXC}=\frac{\mathrm{max}\left\vert {\sigma }_{m}(t+1)-{\sigma }_{m}(t)\right\vert }{\sigma },\end{equation} \tag{ 34 }$$

where σ denotes the overall standard deviation of the sample. We used the window of length m = 3 and calculated the features for both coordinates separately. They should improve the detection of ATTM type of movements.

The performance of the classifiers in task 2 of the AnDi challenge was evaluated with the F₁ score,

$$\begin{equation}{F}_{1}=\frac{\text{TP}}{\text{TP}+0.5(\text{FP}+\text{FN})},\end{equation} \tag{ 35 }$$

where TP, FP, and FN are the true positive, false positive, and false negative rates, respectively. Since we are dealing here with a multilabel classification problem, a micro-averaged F₁ score was calculated to assess the aggregated contributions of all classes (i.e. the sums of all TP, FP and FN were first determined and then inserted into equation (35). For balanced sets (same number of samples for each class), like those of the AnDi challenge, the micro-averaged F₁ score coincides with the accuracy of the classifier.

Results for the performance of the classifiers are shown in table 4. As we can see, the model with the original set of features achieves an accuracy of 0.73. After adding the features described in section 2.5 we observe a significant improvement of its performance. The accuracy of the extended version (0.83) is comparable with the best methods of the AnDi challenge [59].

Zoom In Zoom Out Reset image size

Inspection of the confusion matrices of the classifiers may give us additional insight into partial contributions of the overall performance. To recall, an element c_ij of the confusion matrix indicates how many observations known to belong to class i (true label) are predicted to be in class j (predicted labels) [68]. Normalized confusion matrices for both models are shown in figure 2. The best performance of the base model is observed for LWs. Among all LW samples in the test data, 92% of them were correctly recognized. Most of the misclassified LW trajectories were assigned to FBM motion (5%) and SBM (2%). The accuracy of the classifier for CTRW (85%) and SBM (81%) is smaller, but still quite good. The incorrectly classified CTRW trajectories have been usually confused with ATTM. In the case of SBM, many trajectories have been assigned to ATTM (11%) and FBM (6%) classes.

We have not expected such a moderate performance of the original classifier for FBM (72%). In our earlier attempts [49, 51, 52], we usually got much higher performance for this type of motion. However, it should be stressed that previously we trained our classifiers on different set of trajectories than used for the AnDi challenge. Moreover, we usually used longer trajectories. Thus, it seems that the features in the original set are lacking the power to differentiate FBM successfully from other types of motion considered in this paper, in particular from SBM (19% of misclassifications) and LW (7%).

Finally, it should be emphasized that the performance of the base classifier for ATTM is unsatisfactory (37%). Most of the ATTM trajectories are confused with SBM (37%) and CTRW (22%). It is clear that the base model cannot separate the different origins of trajectories from models with time-varying diffusion coefficient D.

Adding the features from section 2.5 to the input vectors has significantly improved the overall accuracy of the classifier. As can be seen in the right plot of figure 2, all of the partial contributions to the accuracy improved as well. Although the performance for ATTM remains weak (61% only), it is much higher than in the previous case. As for the other classes, the modified classifier achieves very good accuracies for LW (98%) and CTRW (94%) and good ones for SBM (83%) and FBM (82%). The biggest challenge is still to distinguish ATTM and FBM from SBM trajectories.

In figure 3, the F₁ scores for different ranges of the trajectory lengths in the case of the extended model are shown. Not surprisingly, the classifier struggles a lot with very short trajectories. In this regime, many of the paths are indistinguishable from normal diffusion. Moreover, in the case of ATTM and SBM models, the information contained in the short trajectories may be not enough for detecting the time-varying diffusion coefficient D. Interestingly, the accuracy of the FBM classification increases with the sample length to achieve around 90% for trajectories longer than 500 steps.

Zoom In Zoom Out Reset image size

It is worth to mention that the difficulties with very short trajectories are not specific to the XGB in particular or to feature-based methods in general. All methods presented in the challenge were afflicted by similar problems in this regime [59].

In figure 4, the F₁ scores for different levels of noise (i.e. different values of its standard deviation) are presented in case of our extended model. While the classifier performs quite well for the CTRW and LW even in the case of high noise levels, its accuracy significantly decreases for the remaining diffusion models. It seems that the classifier is not able to capture enough characteristics of ATTM, SBM and FBM types of movements from noise trajectories in order to successfully distinguish them from each other. Interestingly, in a series of computer experiments we observed that a classifier trained on very noisy data performs still reasonably well on a test set with low levels of noise.

Zoom In Zoom Out Reset image size

One of the benefits of feature-based methods is the relative ease of their interpretation. In particular, it is possible to assign a score to each of the features indicating how important these features are to the model's predictions. This importance aids to extract relevant knowledge concerning the relationships contained in the data or learned by the model. But it can also help in the selection of the meaningful features for a simplified version of the classifier.

One of the popular feature assessment techniques is the permutation feature importance [68]. It is defined as the decrease in a classifier's score when a single feature values are randomly shuffled. Since this procedure breaks the relationship between the feature and the target, the drop in the score indicates how much the classifier depends on the feature. In table 5, the 20 most important features in our model according to the permutation method are presented. It is worth to note that the features added in the extended version of the model play the key role among them. In particular, the average MW features have the highest score, followed by the D'Agostino–Pearson test statistic κ².

Since the permutation importance is known to have several drawbacks while working with correlated features [86] and this is the case in our model (e.g. we have several features based on MSD), we decided to check the gain importance as well [87]. In this approach, the importance of a feature represents the relative contribution of the feature to the model, calculated by taking each feature's contribution for each tree in the model. A higher value of this metric when compared to other features implies it is more important for generating predictions. The top 20 features according to gain are listed in table 6. Although the moving average window and D'Agostino–Pearson characteristics are still very important, the anomalous exponent α and the MSD ratio get now the highest scores. But again, many of the most important features were added in the extended version of the model. Thus it seems we made some reasonable choices to improve the classification accuracy.

Last but not least, to further elaborate on that issue, we will employ another well known technique—the SHAP values based on feature attribution, i.e. the impact of having a certain value of a given feature in comparison to the prediction we would make if that feature took some baseline value [88, 89]. Compared the previous methods, SHAP values may be used to explain individual predictions. Moreover, they offer a more consistent approach to assess the overall contribution of features to the final outcome of the model.

In table 7, SHAP values for the most important features in different diffusion models are shown, in the case of the extended model. It should be noted that for every single diffusion model, the additional features turned out to be the most important ones, in agreement with the previous methods. In other words, the additional features prove to be an excellent choice. This also shows how critical the choice of features can be for a classifier (see reference [52] for further details).

The Moses exponent M, which relates to non-stationary increments, turned out to be the essential attribute for distinguishing ATTM and SBM from the rest of models; it seems to be quite important for FBM as well. However, in the latter case, the anomalous exponent is the most meaningful feature, followed by the D'Agostino–Pearson test statistic κ². Characteristics based on the average MW have the highest importance for CTRW. And finally, the maximum standard deviation feature is the top one for LW.

In figures 5 and 6, we may check how each of the 20 most relevant features impacts the decisions of the classifier for each diffusion model. A dot in the plot corresponds to a single instance of the explanation given by the classifier. The characteristics on the Y axis are sorted by their mean importance. The SHAP values on the horizontal axis reflect the positive or negative influence of each factor on the model prediction. A positive SHAP value increases the likelihood that the sample will be classified as a specific model, whereas a negative SHAP value—that the observation belongs to one of the other classes. Larger values of SHAP (positive or negative) indicate a greater impact on the final prediction. Additionally, the colormap represents the feature values—red being high and blue being low (the color values are relative). The dependent variable Y is positively correlated with high SHAP values on the positive side of the X axis. In other words, those values have a beneficial effect on prediction. For instance, in the case of CTRW, the larger the value of mw_y_mean_10, the higher SHAP attribution, thus the likelihood of classification as CTRW increases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As can be seen, the SHAP values for CTRW and LW indicate that our classifier is more likely to distinguish these classes from the rest (some features have a long-tailed distribution of values and the colors indicating positive and negative impact are well separated from each other). In contrast, ATTM, FBM and SBM models reveal a high concentration of samples near 0 for many features. Moreover, the colors are mixed with each other indicating that there is no unique relationship between the SHAP value of a feature and its impact on the prediction. In this case, only a combination of all features can give us a reasonable outcome. In other words, it would be rather impossible to build a statistical test to distinguish these three models that operates with a single feature (among the ones present in our set).