Transformers and capsule networks vs classical ML on clinical data for alzheimer classification

Data processing and feature selection

The dataset was acquired in August 2023 and contains 13,205 patient records. It encompasses clinical variables, neuropsychological biomarkers derived from imaging, and APOE-4 status. In total, 113 features were available for modelling.

Following the feature-selection scheme proposed by Yi et al. (2023), we first discarded variables with >50% missing values. We then evaluated the surviving candidates using three complementary selectors—Information Gain (filter), Boruta (wrapper), and Elastic Net (embedded)—thereby covering all three feature-selection paradigms. (i) Information Gain served as a fast entropy-based filter, removing variables with negligible relevance. (ii) Boruta, a random-forest wrapper, retained all relevant features and captured nonlinear interactions overlooked by univariate filters. (iii) An Elastic Net model ( $\alpha$ = 0.5) introduced joint L1/L2 regularisation, shrinking coefficients, and eliminating redundant, highly correlated predictors. Features endorsed by at least two of these methods were kept, yielding 27 clinically meaningful attributes (Table 1). The preprocessing step excluded patients missing selected features, resulting in a final dataset of 1,846 patients: 583 with CN, 122 with AD, 696 with EMCI, and 445 with LMCI. Figure 1 illustrates the class distribution. Stratified sampling preserved this distribution when the data were split into training (80%) and test (20%) sets.

Synthetic minority over sampling technique

Synthetic minority oversampling technique generates synthetic samples for minority classes by interpolating between existing data points. It identifies the k-nearest neighbors of a sample within the minority class and creates new instances along the line segments connecting the original sample to its neighbors. This approach not only balances class distribution but also introduces variability that enhances model generalization (Chawla et al., 2002). As a preprocessing step, SMOTE significantly improved classification performance, particularly for minority classes, by reducing overfitting to the majority classes and enabling models to more effectively capture patterns associated with the early and late stages of cognitive decline.

SMOTE-TOMEK

The SMOTE-Tomek method integrates two complementary techniques, SMOTE and Tomek Links elimination, to address class imbalance in datasets. This hybrid method both balances class distribution through synthetic sample generation and improves separability by removing ambiguous instances. It has proven particularly effective in classification tasks where decision boundaries are unclear and prone to overfitting. In the first stage, SMOTE generates synthetic samples for the minority class by interpolating feature values between selected minority instances and their nearest neighbors in the feature space. While SMOTE significantly improves class balance, placing synthetic points in regions where class overlap occurs may introduce noise near decision boundaries, potentially compromising model generalization. To address this issue, the second stage applies the Tomek Links technique. Tomek Links identifies pairs of instances from opposite classes that are each other’s nearest neighbors, representing regions of class ambiguity. Removing these pairs refines the dataset by eliminating overlapping and conflicting examples, thereby reducing noise and enhancing the clarity of boundaries. As a result, the SMOTE-Tomek method improves data quality by mitigating noise introduced during oversampling and eliminating inherently problematic instances from the original dataset (Batista, Bazzan & Monard, 2003).

Adaptive synthetic sampling

The adaptive synthetic sampling (ADASYN) algorithm addresses the challenge of class imbalance in ML by generating synthetic data points for the minority class, particularly targeting instances that are harder to classify (He et al., 2008). In imbalanced datasets, the dominance of the majority class often leads to biased models and poor predictive performance for the minority class. ADASYN mitigates this issue by adaptively focusing on regions of the feature space where classification is most difficult for minority class samples. The algorithm operates through the following key steps:

• ADASYN identifies minority class instances that are more challenging to classify based on their proximity to majority class samples. Instances in more ambiguous regions receive more synthetic samples, while those in the well-separated areas receive fewer or none.

• It generates synthetic samples by interpolating between selected examples of the minority class and their nearest neighbors. The quantity of synthetic data produced for each instance is proportional to its classification difficulty.

• By adaptively generating samples, ADASYN improves classifier performance, reduces bias toward the majority class, and enhances generalization, particularly for underrepresented categories.

This strategy balances class distributions by generating synthetic samples tailored to reinforce decision boundaries, thereby enhancing the model’s robustness and predictive accuracy.

NearMiss

The NearMiss method is an under-sampling technique designed to address class imbalance by selecting a subset of majority class (negative) examples that are closest to the minority class (positive) instances, based on a predefined distance metric (Mani & Zhang, 2003). By focusing on these strategically chosen samples, NearMiss reduces class-imbalance skew and enhances the model’s ability to learn the distinguishing features of the minority class. The method comprises three primary variations:

• NearMiss-1: Selects majority class samples with the smallest average distances to the three nearest minority class instances. This variation emphasizes boundary regions, capturing the most relevant negative examples for positive class decision margins.

• NearMiss-2: Selects majority class samples with the smallest average distances to the three farthest minority class instances. This approach ensures that the selected negatives are generally close to the overall distribution of the minority class.

• NearMiss-3: For each minority class instance, the algorithm selects a fixed number of the nearest majority class samples. This selection strategy ensures a balanced local distribution by surrounding each positive instance with an equal number of nearby negative examples.

Traditional ML models

• Extra-Trees Classifier (ETC)

  The ETC, or Extremely Randomized Trees, is a variant of the Random Forest algorithm that introduces a higher degree of randomness in the construction of individual decision trees. ETC’s additional randomness can enhance computational efficiency and increase robustness against overfitting in specific scenarios (Géron, 2022). The following section outlines the primary steps in building an ETC:

  • 1.

    Feature selection: At each node, the algorithm randomly selects a subset of features $\mathcal{F}\subset \mathcal{X}$, where $\mathcal{X}$ denotes the full set of $d$ features in the dataset. While this step is similar to that used in Random Forests, ETC increases stochasticity by enforcing a fixed random subset.

  • 2.

    Random split point: For each feature $f\in \mathcal{F}$, a split point $t$ is randomly selected from the range of values observed for $f$ in the current node. ETC differs from Random Forests, which determine the split point by optimizing a criterion to minimize impurity.

  • 3.

    Split criterion: Each candidate pair $(f,t)$ is evaluated using a cost function, such as entropy or the Gini index. Selecting the combination that minimizes node impurity yields the optimal split.

  Formally, the optimal splitting point at a node N is given by: $$\mathrm{O}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{l}=\arg \underset{f\in \mathcal{F},t\in \mathcal{T}}{min}\mathrm{I}\mathrm{m}\mathrm{p}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{y}(N,f,t)$$where $\mathcal{T}$ is the set of possible split points for feature $f$, and the impurity $\mathrm{I}\mathrm{m}\mathrm{p}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{y}(N,f,t)$ is commonly calculated using the Gini index: $$\mathrm{G}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}=1-\sum _{k=1}^K{p}_k^2.$$Here, K represents the number of classes, and $p_k$ denotes the proportion of samples belonging to class $k$ at node N.

• Random Forest (RF)

  The RF is an ensemble learning method that constructs multiple decision trees, each trained on a different subset of the data and features (Breiman, 2001):

  • 1.

    Bootstrap aggregating (Bagging): Multiple bootstrap samples ${\mathcal{B}}_{\mathcal{i}}$ are drawn with replacement from the training set $\mathcal{D}$, and each sample is used to train an individual decision tree $T_i$.

  • 2.

    Random feature selection: At each node of a decision tree, a random subset of features $\mathcal{F}\subset \mathcal{X}$ is selected from the total set of $d$ available features.

  • 3.

    Node splitting: From the subset $\mathcal{F}$, the algorithm selects the optimal feature $f$ and split point $t$ by minimizing an impurity measure, such as entropy or the Gini index.

  The final prediction of the Random Forest is obtained by aggregating the predictions of all individual trees using majority voting for classification tasks or averaging for regression: $$\hat{y}=\frac{1}{n_{\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}}}\sum _{i=1}^{n_{\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}}}{T}_i(\mathbf{x})$$where $T_i(\mathbf{x})$ denotes the prediction of the $i$-th tree for the input instance $\mathbf{x}$.

• Support vector classifier

  Support vector classifier (SVC) is a widely used class of ML algorithms designed to solve both classification and regression tasks involving linear (flat) and nonlinear (curved) decision boundaries. The primary objective of SVM is to maximize the margin between classes, optimizing the separation as a broad “highway” rather than a narrow boundary (Géron, 2022).

  The decision function linearly combines the features of a data point using weights $\omega$ and a bias ${\omega }_0$. The result is passed through a sign function to determine the predicted class: $$f(\omega ;x)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}({\omega }^{T}x+{\omega }_0).$$A positive result assigns the instance to class $+1$; a negative result assigns it to class $-1$.

  SVM employs the hinge loss function to penalize misclassified points and those within the margin boundaries. Hinge loss is given by: $$max(0,1-y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}({\omega }^{T}x+{\omega }_0))$$where:

  • –

    If $y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}=+1$ and ${\omega }^{T}x+{\omega }_0<0$, misclassifying the sample incurs a high penalty.

  • –

    If $y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}=+1$ and $0<{\omega }^{T}x+{\omega }_0<1$, the classification is correct but within the margin, leading to a lower penalty.

  • –

    If $y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}=+1$ and ${\omega }^{T}x+{\omega }_0>1$, no penalty is applied.

  The same logic applies for $y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}=-1$.

  A regularization term is included in the objective function to encourage a wider margin. This regularization term governs the model’s complexity and is weighted by a hyperparameter $\lambda$. The following expression defines the complete loss function: $$L(\omega )=\frac{1}{m}\sum _{i=1}^{m}max(0,1-y_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}^i({\omega }^T{x}^i+{\omega }_0))+\lambda ||\omega |{|}_2^2$$where $m$ is the number of training examples and $\lambda$ regulates the trade-off between margin maximization and classification error.

  The optimization objective involves determining the values of $\omega$ that minimize the loss function. This minimization results in a convex quadratic programming problem, thereby ensuring the existence of a unique global minimum. The problem can be reformulated in its dual form, particularly beneficial in high-dimensional feature spaces. The dual optimization problem is given by: $$\underset{\alpha }{max}\sum _{j=1}^m{\alpha }_j-\frac{1}{2}\sum _{j=1}^m\sum _{k=1}^m{\alpha }_j{\alpha }_k{y}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}^j{y}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}^k(x^j{)}^T{x}^k$$subject to the constraints ${\alpha }_j\ge 0$ and ${\sum }_{j=1}^m{\alpha }_j{y}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}^j=0$.

  Once the optimal values ${\alpha }_j$ are obtained, the weight vector $\omega$ can be recovered using: $$\omega =\sum _{j=1}^m{\alpha }_j{y}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}^j{x}^j.$$With $\omega$ computed, the trained decision function classifies new data points: $$f(x_{\mathrm{n}\mathrm{e}\mathrm{w}})=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}({\omega }^T{x}_{\mathrm{n}\mathrm{e}\mathrm{w}}+{\omega }_0)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\sum _{j=1}^m{\alpha }_j{y}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}^j{(x^j)}^T{x}_{\mathrm{n}\mathrm{e}\mathrm{w}}+{\omega }_0\right).$$

• Gradient boosting

  GB is an ensemble technique that constructs a strong predictive model by sequentially combining multiple weak learners (Friedman, 2002). The core idea is to build models iteratively, correcting the residual errors of the previous models at each step. The model achieves this by minimizing a specified loss function using a gradient descent approach (Géron, 2022). The mathematical formulation is:

  • 1.

    Initial model: The process begins with a constant base model $F_0(x)$, typically defined as the value that minimizes the loss function over the training data—often the mean of the target variable. The following presents the mathematical formulation:

    $$F_0(x)=\arg \underset{\gamma }{min}\sum _{i=1}^{n}L(y_i,\gamma )$$where $L(y,F(x))$ is the loss function measuring the discrepancy between the true value $y$ and the model prediction $F(x)$.

  • 2.

    Residual computation: At each iteration $m$, the algorithm computes the pseudo-residuals, which correspond to the negative gradients of the loss function concerning the current predictions $F_{m-1}(x)$:

    $$r_{im}=-{\left[\frac{\mathrm{\partial }L(y_i,F(x_i))}{\mathrm{\partial }F(x_i)}\right]}_{F(x)=F_{m-1}(x)}.$$These residuals indicate how the model should adjust to minimize the prediction error.

  • 3.

    Fitting the new model: A new weak learner $h_m(x)$ is fitted to the residuals using a least squares approach:

    $$h_m(x)=\arg \underset{h}{min}\sum _{i=1}^n{\left[r_{im}-h(x_i)\right]}^2.$$Aligning the model $h_m(x)$ with the structure of the current residual errors

  • 4.

    Model update: The ensemble model is updated by adding the newly fitted model $h_m(x)$, scaled by a learning rate $\nu$, to the previous model:

    $$F_m(x)=F_{m-1}(x)+\nu h_m(x).$$The learning rate $\nu \in (0,1]$ regulates the contribution of each new model, providing a trade-off between training speed and performance.

  • 5.

    Final prediction: After M iterations, the final prediction function is defined as:

    $$F_M(x)=F_0(x)+\sum _{m=1}^M\nu h_m(x).$$This cumulative model integrates all weak learners constructed during training to form a strong predictor.

DL models

• CNNs

  CNNs extract meaningful features from input data, forming the foundation of their functionality. Researchers have developed numerous CNN architectures in response to advances in computational power and DL; however, most architectures share a standard structure consisting of three main components: a feature extraction layer, a fully connected layer, and an output layer.

  The feature extraction layer consists of multiple convolutional kernels that map various features from the input data. Kernels learn and capture low-level to high-level patterns, followed by the application of an activation function. The fully connected layer connects each neuron from the previous layer to every neuron in the current layer, transforming spatial features into semantic representations. The output layer, commonly used for classification tasks, typically includes an optimizer and a softmax activation function to produce probabilistic outputs (Gu et al., 2018). (1) $$M^l=\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}\left(f\left(\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(\sum _{i=1}^n\left(M_i^{l-1}\ast K_i^l\right)+b^l\right)\right)\right).$$In this equation, $M^l$ denotes the output feature map at layer $l$, $M_i^{l-1}$ is the input feature map from the previous layer, $K_i^l$ represents the convolutional kernel, $b^l$ is the bias term, and * indicates the convolution operation.

  Figure 2 illustrates the feature-extraction block. The network first applies three convolutional layers with 8, 16, and 32 filters, each employing the “selu” activation function (scaled exponential linear unit, SELU). It then applies a max pooling layer and a dropout layer with a rate of 0.5 to mitigate overfitting. Next, the network applies additional convolutional layers with 64, 128, and 256 filters, each followed by a max pooling layer and dropout to further refine the extracted features.

  Figure 3 presents the fully connected architecture. The network includes Fully connected layers with 1,024 to 16 neurons in descending order, respectively. The final output layer adapts according to the number of target classes. Each layer utilizes the “selu” activation function and incorporates batch normalization to stabilize training and accelerate convergence.

• DigitCapsule-Net

  The Capsule Network (Capsule-Net) architecture significantly advances modeling spatial and hierarchical relationships within data. Unlike traditional CNNs, which aggregate features via convolution and pooling operations, Capsule-Net employs convolutional blocks to learn features and then routes these features through capsules. This approach enhances classification accuracy by emphasizing critical activations and maintaining a richer input representation. Capsule-Net demonstrates robustness to minor variations in input data, leading to more accurate and reliable image interpretations.

  The DigitCapsule-Net layer is a core component of the Capsule-Net architecture and follows the primary (or central) capsules. This layer contains multiple capsules, each functioning as a small neural network that generates an activation vector. Each capsule actively represents a specific class, such as digits in image recognition tasks, and produces a vector that captures detailed information about the input. This vector-based representation enables the model to express and interpret the data more comprehensively than traditional scalar outputs (Sabour, Frosst & Hinton, 2017).

  A key innovation of the DigitCapsule-Net is its use of dynamic routing, a process that determines the strength of the connections between the outputs of primary capsules and the higher-level capsules in the DigitCapsule-Net layer. During dynamic routing, the predictions from the lower-level capsules are iteratively refined based on their agreement with the outputs of the higher-level capsules. This mechanism allows the network to preserve more information than max-pooling, resulting in more precise data representations.

  DigitCapsule-Net derives class probabilities from the length of each output vector, which reflects the likelihood of a particular class being present. Additionally, the orientation of these vectors encodes detailed properties of the represented objects, such as their pose, size, and deformation. This dual encoding enhances the interpretability and accuracy of classification tasks, offering a more nuanced understanding of the input data. The DigitCapsule-Net exhibits strong robustness to variations in object orientation, illumination, and position, making it particularly suitable for complex recognition scenarios (Sabour, Frosst & Hinton, 2017; Dombetzki, 2018; Holguin-Garcia et al., 2024).

  In the proposed model, the DigitCapsule-Net layer is integrated following the convolutional and pooling layers, as shown in Fig. 4. The dimensions are adjusted to facilitate the efficient processing of capsule vectors. The model first normalizes and compresses the output values using an activation function to ensure they remain within an appropriate range. The model then passes this compressed output to DigitCapsule-Net, which applies linear transformations and produces a multidimensional tensor representing the learned class-specific features.

• Transformer encoder

  The Transformer encoder captures and models contextual information by establishing relationships among input vectors. Unlike recurrent language models such as long short-term memory (LSTM) networks, it employs self-attention mechanisms that effectively capture the contextual representation of features within a sequence of vectors. Each encoder layer includes multiple attention heads, feed-forward layers, and positional encoding, which retains information about the order of the input elements (Yan et al., 2019).

  Multiple attention heads apply self-attention in parallel; each head partitions the input into a query (Q), key (K), and value (V) vectors, allowing the model to process different representation subspaces simultaneously. The outputs of all attention heads are concatenated and passed through a learned linear transformation to form the final output, a process known as multi-head attention (Vaswani et al., 2017; Bravo-Ortiz et al., 2024c, 2024a).

  Figure 5 illustrates the Transformer encoder architecture used in this study. The following equation defines the self-attention mechanism: (2) $$\mathrm{A}\mathrm{t}\mathrm{t}\mathrm{n}(\mathbf{Q},\mathbf{K},\mathbf{V})=\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{x}\left(\frac{\mathbf{Q}{\mathbf{K}}^T}{\sqrt{d_k}}\right)\mathbf{V}.$$In this equation, Q represents the query vector corresponding to the current token, K denotes the key vector, and V is the value vector containing relevant contextual information. The term $\sqrt{d_k}$ serves as a normalization factor for scaling the dot product, where $d_k$ is the dimensionality of the key vectors.

Convolutional neural network–fully connected

The CNN with a fully connected layer (Fig. 6) represents the most basic configuration and comprises two main components: a feature extraction block and a classification block implemented via a fully connected layer at the output.

The input to the model is a vector with dimensions (None, 27, 1). The convolutional filters process the data into a shape of (None, 256), which is then reduced to (None, 16) before it enters the classification layer. The term “None” denotes the variable batch size.

Convolutional neural network–DigitCapsule-Net

Convolutional neural network–DigitCapsule-Net (CNN+DigitCapsule-Net), illustrated in Fig. 7, receives input vectors with dimensions (None, 27, 1). These inputs are first processed through a feature extraction layer, resulting in a tensor of shape (None, 1, 256). The model passes the output to the DigitCapsule-Net, which produces a final representation of shape (None, 3, 16) for classification. As before, “None” indicates the variable input batch size.

Convolutional neural network–Transformer encoder

The convolutional neural network–transformer encoder (CNN+TF) shown in Fig. 8, takes input vectors of shape (None, 27, 1). These vectors are passed through feature extraction layers, producing an output of shape (None, 1, 256). Within the Transformer block, the data undergo processing through multiple attention heads while maintaining the same dimensionality. A multilayer perceptron performs the classification in the final stage.

Convolutional neural network–Transformer encoder + DigitCapsule-Net (CNN+TF+DigitCapsule-Net)

The hybrid model illustrated in Fig. 9 integrates several advanced components for data processing. Input vectors of shape (None, 27, 1) are first passed through a shared feature extraction block, resulting in an intermediate representation of shape (None, 1, 256). This output passes through a Transformer encoder, which captures complex feature relationships via attention mechanisms without altering dimensionality. Finally, the DigitCapsule-Net receives the output and produces a final representation of shape (None, 3, 16) for classification. This architecture combines the efficiency of CNNs, the contextual modeling capabilities of Transformers, and the hierarchical feature representation of capsule networks, thereby leveraging the strengths of each.

Metrics

By emphasizing the differences between false positives (FP), false negatives (FN), true positives (TP), and true negatives (TN), Tabares-Soto et al. (2021) illustrate the importance of metrics in the evaluation of a model. The following are the most crucial metrics:

• Accuracy

  Accuracy evaluates a classification model’s predictive performance. To compute it, one divides the total number of correct predictions by the total number of predictions. Accuracy is a value ranging from 0 to 1, representing the percentage of accurate predictions. (3) $$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}.$$

• Precision

  The precision metric gauges the ratio of correctly identified positive cases by a model to the total cases identified as positive, encompassing both true and false positives. (4) $$Precision=\frac{TP}{TP+FP}.$$

• Recall

  Also known as sensitivity, it shows the ability of the classifier to display correct predictions. (5) $$Recall=\frac{TP}{TP+FN}.$$

• F1

  F1 is a metric used to assess a model’s ability to precisely identify positive and negative cases, particularly in scenarios where the data is imbalanced, and the positive class is infrequent. It is determined as the harmonic mean of precision and recall, proving especially valuable in situations with uneven class distribution. (6) $$F1=2\ast \frac{Precision\ast Recall}{Precision+Recall}.$$

• Confusion matrix

  The confusion matrix is a tabular representation that encapsulates the correspondence between a model’s predictions and the actual labels of the data. It incorporates true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN). This matrix proves valuable for visually summarizing a model’s performance, highlighting both its successes and errors. Each row in the matrix corresponds to the actual class, and each column indicates the number of predictions made for each class. Additionally, it helps identify instances where the model misclassifies one class as another.

Cross-validation

Cross-validation (CV) is a resampling technique used to evaluate a model’s performance by partitioning the dataset into $k$ mutually exclusive subsets (folds) of approximately equal size. At each iteration, one fold serves as the validation set, while the remaining $k-1$ folds constitute the training set. The process continues until all folds have individually served as the validation set. Averaging the $k$ resulting performance scores provides a more reliable estimate of the model’s generalization ability. To preserve the original class distribution within every fold and avoid biased estimates in the presence of class imbalance, we adopt stratified 10-fold cross-validation, as illustrated in Eq. (7).

(7) $$\hat{\mathcal{M}}=\frac{1}{k}\sum _{i=1}^k{\mathcal{M}}_i,$$where ${\mathcal{M}}_{\mathcal{i}}$ denotes the performance metric computed on the $i$-th validation fold.

Statistical comparison of model performance

This study rigorously evaluated model performance differences using several non-parametric statistical tests, as recommended by Rainio, Teuho & Klén (2024). These methods are particularly appropriate in scenarios where the assumptions underlying parametric tests—such as normality and homoscedasticity—are violated, which is often the case with cross-validated performance metrics.

• Friedman test

  The Friedman test serves as a non-parametric alternative to repeated-measures ANOVA, designed to detect differences in treatments across multiple testing conditions. Within the context of ML, it assesses whether statistically significant performance differences exist among various models evaluated on the same datasets. The test ranks model performance scores per dataset and analyzes these ranks to determine the significance of observed differences. Its robustness to violations of normality makes it particularly suitable for comparing more than two models simultaneously (Rainio, Teuho & Klén, 2024).

• Wilcoxon signed-rank test

  This study applied the Wilcoxon signed-rank test for pairwise model comparisons. This non-parametric method evaluates whether a significant difference exists between two related samples, such as the performance scores of two models across the same cross-validation folds. By accounting for both the magnitude and direction of the differences, the test provides a more comprehensive assessment than the sign test. It is especially well-suited for situations in which the assumption of normality is not upheld (Rainio, Teuho & Klén, 2024).

• Levene’s test for equality of variances

  Analyzing performance variability is essential, as models with similar mean accuracies may differ substantially in consistency. Levene’s test was employed to assess the homogeneity of variances across model performance metrics. Unlike other variance equality tests, Levene’s method is less sensitive to deviations from normality, making it an appropriate choice for this analysis (Rainio, Teuho & Klén, 2024).

ML hyperparameter tuning with grid search and pipelines

This study employed a combination of grid search and pipeline integration to enhance the performance of traditional ML models. Grid search is a widely used methodology that systematically explores a predefined subset of the hyperparameter space to identify the optimal configuration. Although effective, grid search can be computationally expensive due to the exponential growth in combinations as the number of hyperparameters increases. Consequently, practitioners often restrict the search to carefully selected regions of the parameter space, typically assuming mutual independence among the hyperparameters (Liashchynskyi & Liashchynskyi, 2019).

We ensured reproducibility and efficiency by embedding all hyperparameter tuning procedures within a pipeline architecture, consistent with the implementation described in Pedregosa et al. (2011). Pipelines facilitate the seamless integration of preprocessing steps (feature selection, scaling), model training, and evaluation. This approach minimizes the risk of data leakage and promotes methodological consistency across experiments.

Each model underwent independent optimization for both balanced and imbalanced datasets. The optimized hyperparameters for the three-class classification task are presented in Table 2.

DL model tuning

The hyperparameters of DL models were carefully tuned to ensure optimal configuration during the training process. The categorical cross-entropy loss function was employed in conjunction with the Adam optimizer, configured with a learning rate of 0.001. All models were trained with a batch size of 256 for a total of 500 epochs, incorporating early stopping to prevent overfitting.

The following section provides a detailed description of the hyperparameters used.

• Batch normalization (BN)

  The features in each data map are normalized using the BN to have a mean of 0 and a variance of 1, enabling rescaling and retranslating the distribution. A faster rate of learning is made possible by this training method (Ioffe & Szegedy, 2015). The following equation denotes it: (8) $$BN(x,\gamma ,\beta )=\beta +\gamma \frac{x-\text{E[X]}}{\sqrt{Var[X]+ϵ}}$$where:

  $\gamma$ = The re-scaling scalar.

  $\beta$ = re-translation scalar.

  E[X] = expectation.

  $Var[X]$ = variance.

• Scaled exponential linear unit

  This nonlinear activation function, proposed by Klambauer et al. (2017), behaves linearly for positive input values and exhibits exponential behavior for negative ones. The function incorporates two constants: $\lambda$, approximately equal to 1.0507, and $\alpha$, the negative slope coefficient, with an approximate value of 1.67326. These parameters enable effective scaling and propagation of signals across multiple neural network layers. Furthermore, researchers consider the function self-normalizing because it maintains a constant mean and variance during forward propagation, thereby enhancing model stability (Rasamoelina, Adjailia & Sinčák, 2020).

  The following equation gives it: (9) $$\mathrm{S}\mathrm{E}\mathrm{L}\mathrm{U}(x)=\lambda \{\begin{array}{cc}x\hfill & \mathrm{i}\mathrm{f}x>0\hfill \\ \alpha \cdot (\exp (x)-1)\hfill & \mathrm{i}\mathrm{f}x\le 0\hfill \end{array}.$$

• Dropout

  The regularization technique known as Dropout is employed to prevent overfitting in neural networks. This method randomly drops out selected nodes and their connections during training, reducing the risk of co-adaptation among neurons and preventing the network from relying excessively on specific pathways. As a result, the model is less likely to memorize the training data and more capable of generalizing to unseen samples.

  Beyond mitigating overfitting, Dropout contributes to the development of more robust networks by forcing them to function under varying structural configurations during training. This variability promotes model resilience and improves generalization. Additionally, Dropout facilitates the averaging of predictions and helps reduce variance at test time (Srivastava et al., 2014). Dropout is incorporated to prevent overfitting in this study, as illustrated in Fig. 2.

Feature importance

The importance of features is a fundamental technique in interpreting ML models. It helps understand the model’s behavior and identify biases and critical features. This method is crucial as artificial intelligence models have become increasingly complex and challenging to interpret due to scientific advancements (Adler & Painsky, 2022).

Gradient-weighted class activation mapping

Gradient-weighted class activation mapping (Grad-CAM) is a widely adopted technique for interpreting CNNs by visualizing input regions that most influence the model’s predictions for making predictions based on the information provided by the gradient. It computes the gradient of the predicted class score concerning the feature maps of the last convolutional layer, averages these gradients to derive weights, and multiplies them by the corresponding feature maps to produce a saliency map that highlights the essential input regions (Selvaraju et al., 2016).

Shapley additive explanations values

Shapley additive explanations (SHAP) values constitute a widely adopted technique for interpreting the output of ML models. They provide insight into how individual input features contribute to a specific prediction. SHAP values is particularly valued for its consistency and fairness, as it is grounded in principles from cooperative game theory (Marcílio & Eler, 2020; Meng et al., 2020).

The theoretical foundation of SHAP lies in Shapley values, a concept from cooperative game theory that defines a fair method for distributing the payoff of a task among its participants. In the context of ML, the “participants” are the input features, and the “payoff” corresponds to the model’s predicted output. SHAP values quantify each feature’s contribution to the prediction by considering all possible combinations of features (Marcílio & Eler, 2020; Meng et al., 2020).

Formally, SHAP values assigns a value to feature $x_i$ by averaging its marginal contributions across all possible feature subsets.

Given a set of input features $\{x_1,x_2,\dots ,x_n\}$, the SHAP value ${\varphi }_i$ is computed as:

$${\varphi }_i=\sum _{S\subseteq N\setminus \{i\}}\frac{|S|!(|N|-|S|-1)!}{|N|!}\left(f(S\cup \{i\})-f(S)\right)$$where:

• ${\varphi }_i$ is the SHAP value for feature $x_i$.

• N denotes the full set of input features.

• S is any subset of features that does not include $x_i$.

• $f(S)$ is the model’s prediction using only the features in subset S.

• $f(S\cup \{i\})$ is the model’s prediction when feature $x_i$ is added to subset S.

• $\frac{|S|!(|N|-|S|-1)!}{|N|!}$ is a weighting factor that ensures a fair distribution based on all possible feature coalitions involving $x_i$.