Participants

Participants were recruited from the Department of Psychiatry at Soonchunhyang University Cheonan Hospital in Korea between December 2018 and August 2021. The study included 3 groups: patients who maintained an MDD diagnosis throughout an approximately 3-year follow-up period (stable MDD group, n = 50), those who transitioned from MDD to BD-II within the same period (converted BD group, n = 16), and those with an established BD diagnosis at enrollment (diagnosed BD group, n = 26). The diagnosed BD group consisted of 7 participants with bipolar I disorder (BD-I) and 19 with BD-II; all patients in this group experienced a depressive episode at the time of EEG recording. The diagnostic conversion from MDD to BD-II in the converted BD group was confirmed based on the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition (DSM-5) criteria for hypomanic episodes.²¹ The DSM-5 stipulates that “A hypomanic episode that emerges during antidepressant treatment and persists at a fully syndromal level beyond the physiological effect of that treatment is sufficient evidence for a hypomanic episode and a diagnosis of bipolar II disorder.”^21, Reference Terao and Tanaka²² All BD and MDD diagnoses were established by a board-certified psychiatrist in accordance with the DSM-5.²¹ Patients with a history of neurological disorders or other severe medical conditions were excluded during the initial screening interviews. The Institutional Review Board (IRB) and Ethics Committee of Soonchunhyang University Cheonan Hospital granted approval for this study (IRB no. 2018-10-032). All procedures were conducted in accordance with the approved guidelines, and all participants provided informed consent.²³

Psychological measures

Anxiety, depression, and impulsivity were assessed using 3 self-report instruments: the State–Trait Anxiety Inventory (STAI),Reference Spielberger, Gorsuch and Lushene²⁴ the Beck Depression Inventory (BDI),Reference Beck, Ward, Mendelson, Mock and Erbaugh²⁵ and the Barratt Impulsiveness Scale, version 11 (BIS-11),Reference Patton, Stanford and Barratt²⁶ respectively. The STAI quantifies the severity of state and trait anxiety, while the BDI measures the severity of depression. The BIS-11 assesses impulsivity as a personality trait, yielding an overall score (BIS-11 total score) and 3 subscale scores: motor impulsiveness (BIS-motor), non-planning impulsiveness (BIS-nonplan), and attentional impulsiveness (BIS-attentional).

EEG acquisition

Resting-state EEG data were collected during the initial hospital visit of each participant. Participants were seated in a sound-attenuated room and asked to keep their eyes open for 5 min. EEG recordings were performed using a NeuroScan SynAmps2 amplifier (Compumedics USA, Charlotte, NC) and 62 Ag/AgCl electrodes arranged according to the extended 10–20 system and mounted on a Quik-Cap. Electrooculogram electrodes were positioned at the outer corners of each eye and above and below the left eye to capture horizontal and vertical movements, respectively. The reference electrodes were placed on both mastoids, and the ground electrode was positioned on the forehead. The EEG signals were sampled at 1000 Hz, band-pass filtered between 0.1 and 100 Hz, and notch filtered at 60 Hz to eliminate power line noise. Electrode impedance was maintained below 5 kΩ throughout the recording.

EEG preprocessing

First, baseline correction was performed by subtracting the mean amplitude of each channel from the corresponding EEG signal to remove the direct current offset.Reference Kim, Song and Kim²⁷ Next, the artifact subspace reconstruction (ASR) algorithm was applied to attenuate high-amplitude transient artifacts. ASR operates by first identifying artifact-free calibration segments within the recording and estimating their statistical distribution. Following this calibration, continuous EEG data were segmented into 500-ms sliding windows and decomposed into subspaces using principal component analysis. Within each window, any principal component whose variance exceeded 20 standard deviations relative to the calibration distribution was identified as artifact-contaminated and reconstructed using the spatial structure of the clean reference data. Corrupted segments were thereby corrected in place rather than rejected, effectively attenuating transient artifacts while preserving the underlying neural activity.Reference Chang, Hsu, Pion-Tonachini and Jung²⁸

Following ASR, a wavelet-based artifact removal approach was applied.Reference Lopez, Monachino, Morales, Leach, Bowers and Gabard-Durnam²⁹ Each channel was decomposed using a 10-level discrete wavelet transform with a biorthogonal 4.4 wavelet basis.Reference Hwang, Choi, Kim and Lee¹¹ Level-dependent noise estimates were obtained through Bayesian estimation, and a soft thresholding rule was applied to separate artifact-related wavelet coefficients from neural components.Reference Johnstone and Silverman^30, Reference Donoho³¹ The isolated artifact-related coefficients were reconstructed into a time-domain artifact signal via the inverse wavelet transform.Reference Johnstone and Silverman^30, Reference Donoho³¹ This artifact signal was subsequently subtracted from the recording to eliminate artifact components while retaining neural activity.

After wavelet denoising, independent component analysis (ICA) was performed using the extended Infomax algorithm to identify and remove residual artifacts.Reference Lee, Girolami and Sejnowski³² Component classification was conducted using an automated classifier trained on a large independent component database, which assigned probability estimates for each independent component across predefined categories (brain, muscle, eye, heart, line noise, channel noise, and other).Reference Pion-Tonachini, Kreutz-Delgado and Makeig³³ Components were rejected if the estimated probability of muscle artifact exceeded 0.50 or if the estimated probability of any other non-brain category exceeded 0.90.Reference Kim, Song and Kim²⁷ To preserve neural information, the maximum number of rejected components was constrained to 20% of the total independent components per participant (i.e., a maximum of 12 out of 62 components).Reference Choi, Hwang, Yang, Jung, Im and Lee³⁴ Once artifact components had been removed, a common average reference was applied across all channels to standardize the voltage reference, completing the artifact removal phase.Reference Bigdely-Shamlo, Mullen, Kothe, Su and Robbins³⁵ Importantly, the entire preprocessing pipeline described above—including ASR, wavelet denoising, ICA decomposition, and component classification—was executed as a fully automated procedure, thereby ensuring that all artifact rejection decisions were made without knowledge of diagnostic group membership.

For each participant, preprocessed continuous EEG recordings were divided into non-overlapping 30-s candidate segments. A segment was initially deemed artifact-free if the maximum absolute amplitude across all channels did not exceed 90 μV.Reference Hwang, Kim, Shim and Kim³⁶ Each candidate segment was further evaluated using a multi-metric quality control procedure that assessed peak-to-peak amplitude, channel-wise variance homogeneity, kurtosis, spectral slope, and temporal gradient stability; only segments that passed all criteria were retained for subsequent analysis.

The retained segments were then processed differently depending on the target entropy measure. For band-scale entropy analysis, each retained 30-s segment was subjected to MSE computation without further segmentation, and the subject-level band-scale entropy value for each feature was determined as the median across all retained segments. No additional frequency-domain filtering was applied prior to this computation, as the coarse-graining process inherent in MSE analysis functions as a low-pass filter.Reference Bosl, Loddenkemper and Vieluf³⁷ In contrast, for CSE analysis, each 30-s segment was divided into 15 non-overlapping epochs of 2 s each; the CSE was computed for each epoch, and the segment-level CSE value was determined as the median across the 15 epochs. Subject-level CSE values were likewise obtained by computing the median across all retained segments.

To assess robustness to segment selection, the sensitivity of entropy values to varying numbers of these segments was evaluated through repeated random sampling from the pool of eligible data. Specifically, varying subsamples were randomly drawn across multiple iterations, and the resulting entropy estimates were compared to quantify their stability as a function of the number of segments included. All preprocessing, segmentation, and entropy-robustness procedures described in this section were implemented in MATLAB R2021b (MathWorks, Natick, MA, USA); classifier robustness analyses are described separately in section “Feature selection and classification”.

Band-scale entropy analysis

Band-scale entropy was employed to assess the complexity of the EEG data.Reference Hwang, Choi, Im, Yang, Kim and Lee¹² This measure is derived from MSE analysis, which computes sample entropy (SampEn) across multiple temporal scales.Reference Costa, Goldberger and Peng¹⁹ The first step of MSE analysis involves generating a coarse-grained time series from the original EEG signal.Reference Costa, Goldberger and Peng¹⁹ For an original EEG series x of length N, a coarse-grained time series $ {y}^{\tau } $ is generated for each scale factor τ, with each element $ {y}_i^{\tau } $ calculated as follows:Reference Costa, Goldberger and Peng¹⁹

(1) $$ {y}_i^{\tau }=\frac{1}{\tau}\sum_{n=\tau \left(i-1\right)+1}^{\tau i}{x}_n,1\le i\le \left\lfloor \frac{N}{\tau}\right\rfloor =L $$

In this equation, L denotes the length of the coarse-grained series, and i represents the index of its elements. For example, when τ = 1, the time series $ {y}^1 $ is identical to the original series x. The present analysis employed a total of 70 scale factors, resulting in the shortest coarse-grained time series (corresponding to the largest τ value) containing 428 data points. This sample size is considered sufficient for reliable SampEn estimation, given the general recommendation for a time series length of $ {10}^m $ to $ {20}^m $ points, where m is the embedding dimension (set to 2 in this study; see below).Reference Awan, Aziz and Shah³⁸

Subsequently, the SampEn values for each coarse-grained time series $ {y}^{\tau } $ were computed, denoted as $ MSE(\tau )=SampEn({y}^{\tau },\hskip2pt m,\hskip2pt r) $ (2).Reference Costa, Goldberger and Peng¹⁹ SampEn is defined as the negative natural logarithm of the conditional probability that 2 sequences of m consecutive data points that are similar within tolerance r also remain similar when extended by one data point (to length m + 1).Reference Costa, Goldberger and Peng¹⁹ Sequence similarity was assessed using the Chebyshev distance, with parameters set to the commonly used values of m = 2 and r = 0.15 × the standard deviation of each time series. ^{Reference Hwang, Choi, Im, Yang, Kim and Lee12,Reference Costa, Goldberger and Peng19} Once SampEn was obtained at each scale factor, region-level measures were derived by averaging the MSE(τ) values across electrodes within each of 6 electrode clusters, with each cluster representing a distinct brain region.Reference Hwang, Choi, Im, Yang, Kim and Lee¹² These clusters included left frontal (AF3, F1, F3, and F5), right frontal (AF4, F2, F4, and F6), left central (FC1, FC3, FC5, C1, C3, C5, CP1, CP3, and CP5), right central (FC2, FC4, FC6, C2, C4, C6, CP2, CP4, and CP6), left parietal (P1, P3, P5, PO3, and PO5), and right parietal (P2, P4, P6, PO4, and PO6).Reference Hwang, Choi, Kim and Lee¹¹ This process yielded an MSE curve for each region.

The definition of frequency-band-specific scale factor ranges was grounded in previous research demonstrating that lower scale factors reflect the complexity of higher-frequency components, whereas higher scale factors reflect that of lower-frequency components.Reference Hwang, Choi, Kim and Lee11, Reference Hwang, Choi, Im, Yang, Kim and Lee12 Specifically, prior studies found that the MSE curves for simulated sinusoidal rhythms of 50, 30, 12, and 8 Hz peaked at scale factors 7, 13, 30, and 41, respectively.Reference Hwang, Choi, Im, Yang, Kim and Lee¹² Based on these findings, scale factor ranges were defined for the following frequency bands: gamma (scale factors 7–12), beta (scale factors 13–29), alpha (scale factors 30–40), and theta (scale factors 41–70).Reference Hwang, Choi, Im, Yang, Kim and Lee¹² The band-scale entropy for each brain region and frequency band was then calculated by averaging the MSE(τ) values within the corresponding range of scale factors.Reference Hwang, Choi, Im, Yang, Kim and Lee¹²

Intra- and inter-regional cross-sample entropy analysis

CSE is a refinement of cross-approximate entropy that exhibits improved statistical consistency across varying data lengths and parameter settings.Reference Richman and Moorman¹³ This measure quantifies the degree of asynchrony between 2 distinct yet related time series by assessing the regularity of their cross-recurrence patterns.Reference Richman and Moorman¹³ Given 2 sequences of vectors, $ {X}_i^{(m)}=\left\{{x}_i,{x}_{i+1},\dots, {x}_{i+m-1}\right\}\hskip0.3em \mathrm{and}\hskip0.3em {Y}_j^{(m)}=\left\{{y}_j,{y}_{j+1},\dots, {y}_{j+m-1}\right\} $ , where $ 1\le i,\hskip0.45em j\le \hskip0.15em N-m $ , the CSE can be calculated as follows:Reference Richman and Moorman¹³

(3) $$ CSE(m,\hskip2pt r,\hskip2pt N)=-\mathrm{ln}\frac{\sum \limits_{i=1}^{N-m}{n}_i^{(m+1)}}{\sum \limits_{i=1}^{N-m}{n}_i^{(m)}} $$

In this formula, $ {n}_i^{(m)} $ represents the number of template matches obtained by comparing vector $ {X}_i^{(m)} $ from time series X with all vectors $ {Y}_j^{(m)} $ from time series Y. Reference Richman and Moorman¹³ A template match occurs between $ {X}_i^{(m)} $ and $ {Y}_j^{(m)} $ when the Chebyshev distance between them is less than or equal to the threshold r. Reference Richman and Moorman¹³ Correspondingly, $ {n}_i^{\left(m+1\right)} $ represents the number of template matches found when comparing vectors of length m + 1 using the same criterion.Reference Richman and Moorman¹³

The natural logarithm transforms the ratio of template matches into an entropy measure. Lower CSE values indicate greater temporal synchrony between time series, reflecting more predictable cross-signal interactions.Reference Richman and Moorman¹³ Based on these properties, the CSE analysis in the present study was performed on the same 6 brain regions defined in the MSE analysis, with parameters m = 2 and r = 0.15, within commonly recommended ranges for (cross-) SampEn-type estimators.Reference Yentes, Hunt, Schmid, Kaipust, McGrath and Stergiou³⁹ Before the CSE was calculated, each time series was normalized using a z-score transformation to achieve a mean of zero and a standard deviation of one. On the basis of the normalized signals, 2 distinct CSE measures were computed to assess neural synchrony at multiple spatial scales: intra-regional CSE and inter-regional CSE. Intra-regional CSE quantified the synchrony among electrodes within the same brain region and was calculated as the mean CSE across all unique electrode pairs within each region, excluding self-pairs. As a complementary measure, inter-regional CSE quantified synchrony between different brain regions and was calculated as the mean CSE across all unique electrode pairs connecting 2 distinct regions.

Statistical analysis

Group differences in EEG entropy measures were evaluated across the 3 diagnostic groups (stable MDD, converted BD, and diagnosed BD) using analysis of covariance (ANCOVA), with BDI score included as a continuous covariate to control for depression severity. BDI scores were mean-centered to aid interpretation of covariate-adjusted group means.Reference Schielzeth⁴⁰ The feature set comprised 45 EEG-derived measures, consisting of 24 band-scale entropy features and 21 CSE features. For each feature, an ANCOVA model was fitted with diagnostic group as a 3-level categorical factor and mean-centered BDI as a continuous covariate. The omnibus group effect was tested using an F statistic comparing a full model (intercept, 2 group indicators with stable MDD as the reference category, and mean-centered BDI) with a reduced model containing only the intercept and mean-centered BDI (numerator degrees of freedom = 2; denominator degrees of freedom = N − 4).

To account for unequal sample sizes, potential heteroscedasticity, and residual non-normality without relying on strict parametric assumptions, inferential p-values were obtained using a restricted-residual wild bootstrap approach.Reference Davidson and Flachaire⁴¹ Under the null hypothesis of no group effect, residuals from the reduced model were multiplied by independently drawn wild weights from the Mammen 2-point distribution and added to the reduced-model fitted values to generate bootstrap response vectors.Reference Davidson and Flachaire^41, Reference Mammen⁴² For each replicate, both reduced and full models were refitted and feature-wise F statistics were recomputed.Reference Davidson and Flachaire^41, Reference Mammen⁴² A total of 20,000 replicates (B) were generated with a fixed random seed. Upper-tail p-values were computed with a + 1 correction, (1 + count)/ (1 + B), to avoid zero p-values.Reference Phipson and Smyth⁴³

The primary analysis was a global set-level test across all 45 entropy features. The global test statistic was defined as the sum of the feature-wise F statistics, and the global p-value was estimated from the bootstrap distribution of this sum statistic.Reference Davidson and Flachaire⁴¹ Secondary (exploratory) analyses were conducted at the feature level by reporting wild bootstrap p-values for the omnibus group effect, partial eta squared (η ²ₚ) as an effect-size index, and covariate-adjusted group means evaluated at the sample mean BDI level. ^{Reference Richardson44,Reference Searle, Speed and Milliken45}

For features showing a significant omnibus group effect at the feature level (p < 0.05), exploratory post hoc comparisons were conducted using 3 prespecified linear contrasts (converted BD versus stable MDD, diagnosed BD versus stable MDD, and diagnosed BD versus converted BD). Contrast-specific p-values were obtained from the wild bootstrap distribution of the corresponding t statistics using 2-sided testing based on absolute values. ^{Reference Davidson and Flachaire41,Reference Mammen42} The Holm–Bonferroni method was applied to correct for multiple comparisons across the 3 contrasts within each feature.Reference Holm⁴⁶

To examine the clinical relevance of group-discriminating entropy features, within-group Spearman rank correlation analyses were conducted for features showing significant post hoc group differences after Holm–Bonferroni correction.Reference Spearman⁴⁷ For each such feature, the analysis was performed separately within each group involved in the significant contrast. Each selected entropy feature was correlated with 7 psychological measures: state anxiety and trait anxiety from the STAI, depression severity from the BDI, and BIS-attentional, BIS-motor, BIS-nonplan, and BIS-11 total score from the BIS-11. The Spearman rank correlation coefficient (ρ) was used as the effect-size measure.Reference Spearman⁴⁷

Inferential p-values for each feature–clinical-measure pair were obtained using a restricted-residual wild bootstrap procedure applied to rank-transformed variables.Reference Davidson and Flachaire^41, Reference Conover and Iman⁴⁸ Both the entropy feature and the clinical variable were rank-transformed, and a simple linear regression of the ranked feature on the ranked clinical variable was fitted.Reference Conover and Iman⁴⁸ Under the null hypothesis of no association, residuals from an intercept-only reduced model were multiplied by wild weights drawn from the Mammen 2-point distribution and added to the reduced-model fitted values; the full model was then refitted and the slope t statistic was recomputed for each replicate.Reference Davidson and Flachaire^41, Reference Mammen⁴² A total of 20,000 replicates were generated with a fixed random seed, and 2-sided p-values were computed with the +1 correction.Reference Phipson and Smyth⁴³ A global test across all feature–clinical-measure pairs within each group was conducted using the sum of squared t statistics as the global test statistic, with the global p-value estimated from the bootstrap distribution of this sum.Reference Davidson and Flachaire⁴¹ Within each feature, Holm–Bonferroni correction was applied across the 7 clinical scales.Reference Holm⁴⁶

As a supplementary sensitivity analysis, the stability of feature-level group separation between stable MDD and diagnosed BD under varying segment counts was assessed using unadjusted Mann–Whitney U tests. For each segment count k (ranging from 1 to 8), k segments were randomly drawn without replacement from the pool of retained artifact-free 30-s segments for each participant across 100 iterations; the subject-level entropy value for each feature at each iteration was defined as the median across the sampled segments. At each k and iteration, 2-sided Mann–Whitney U tests were performed separately for each of the entropy features that showed significant covariate-adjusted post hoc group differences in the primary ANCOVA. These tests were conducted without covariate adjustment to provide a conservative, model-free assessment of whether raw group separation persisted across segment-selection conditions. Stability was quantified as the proportion of iterations in which unadjusted significance (p < 0.05) was observed at each k. All statistical analyses described in this section were performed using MATLAB R2021b.

Feature selection and classification

Machine-learning analyses were performed using the full cohort (total n = 92) to evaluate the diagnostic classification utility of entropy-derived EEG biomarkers under class imbalance. Three binary classification tasks were evaluated with task-specific sample sizes, as described below. First, stable MDD was compared with the BD spectrum (n = 92), where the converted BD and diagnosed BD groups were combined into a single positive class. Second, stable MDD was compared with converted BD (n = 66), with the diagnosed BD group (n = 26) excluded. Third, stable MDD was compared with diagnosed BD (n = 76), with the converted BD group (n = 16) excluded.

The predictor set comprised 45 EEG variables derived from band-scale entropy (6 regions × 4 frequency bands) and intra- and inter-regional CSE. To examine the contribution of depressive symptom severity, 3 feature configurations were defined for the first and third tasks (stable MDD versus BD spectrum; stable MDD versus diagnosed BD): EEG-only, EEG + BDI, and BDI-only. Because BDI scores in these tasks spanned clinically distinct ranges across diagnostic groups, the inclusion of BDI was motivated by its potential role as a complementary clinical predictor. For the second task (stable MDD versus converted BD), only the EEG-only configuration was evaluated; BDI was omitted a priori on the basis that both groups presented with comparable depressive symptom severity at the time of EEG recording, limiting its expected discriminative value. Feature configuration comparisons in the first and third tasks were conducted entirely within the nested cross-validation framework described below.

Model development and evaluation were conducted using nested cross-validation designed to prevent data leakage. ^{Reference Varma and Simon49,Reference Varoquaux, Raamana, Engemann, Hoyos-Idrobo, Schwartz and Thirion50} Within each cross-validation split, features were standardized to zero mean and unit variance; feature selection, class-imbalance handling, decision-threshold determination, and classifier fitting were then performed exclusively on training data and applied to the corresponding test data. ^{Reference Ambroise and McLachlan51,Reference Santos, Soares, Abreu, Araujo and Santos52} Hyperparameter optimization was conducted in the inner loop using stratified 5-fold cross-validation, and generalization performance was estimated in the outer loop using repeated stratified k-fold procedures (5 folds × 20 repeats).Reference Varoquaux, Raamana, Engemann, Hoyos-Idrobo, Schwartz and Thirion⁵⁰ Hyperparameters were optimized by randomized search (60 iterations) using the area under the receiver operating characteristic curve (ROC-AUC) as the primary optimization criterion.Reference Bergstra and Bengio⁵³ To ensure reproducibility, random seeds were fixed across all stochastic operations.

Class imbalance was addressed within the training data via class-balanced learning (class-weighted losses for classifiers supporting this mechanism) and resampling-based balancing, ^{Reference He and Garcia54,Reference Haixiang, Yijing, Shang, Mingyun, Yuanyue and Bing55} including synthetic minority over-sampling technique (SMOTE) and borderline-SMOTE. ^{Reference Chawla, Bowyer, Hall and Kegelmeyer56,Reference Han, Wang and Mao57} The option of retaining the original class distribution without resampling was also included as a candidate within the hyperparameter search. Resampling was confined to training partitions both during inner-loop tuning and during refitting on the outer training data, to prevent data leakage.Reference Santos, Soares, Abreu, Araujo and Santos⁵²

Feature selection was embedded within the cross-validation pipeline as a single step positioned before resampling. Three candidate strategies were compared as mutually exclusive alternatives within the hyperparameter search: no feature selection (i.e., retaining the complete feature set); a univariate filter-based approach using analysis of variance (ANOVA) F-statistics to rank features by linear association with the target class (with the top-k features retained, where k was optimized within the inner cross-validation);Reference Guyon and Elisseeff⁵⁸ and wrapper-based selection using recursive feature elimination (RFE) with a logistic regression estimator.Reference Guyon, Weston, Barnhill and Vapnik⁵⁹ This wrapper method iteratively pruned 10% of features per step based on the absolute magnitude of model coefficients until reaching a target number of features (also optimized within the inner cross-validation).Reference Han, Huang and Zhou⁶⁰ Across all strategies, the selection method and its associated hyperparameters were jointly optimized within the inner cross-validation to prevent selection bias. ^{Reference Varma and Simon49,Reference Ambroise and McLachlan51} Following feature selection and resampling, 7 classifier families were compared: support vector machine with a radial basis function (RBF) kernel (RBF-SVM), linear SVM,Reference Cortes and Vapnik⁶¹ regularized logistic regression (L1/L2), ^{Reference Krishnapuram, Carin, Figueiredo and Hartemink62,Reference le Cessie and Houwelingen63} ridge classifier,Reference Grüning, Kropf, Spiliopoulou, Kruse, Borgelt, Nürnberger and Gaul⁶⁴ linear discriminant analysis,Reference Fisher⁶⁵ random forest,Reference Breiman⁶⁶ and multilayer perceptron (MLP).Reference Rumelhart, Hinton and Williams⁶⁷

Performance metrics (accuracy, balanced accuracy, sensitivity, specificity, and F1-score) were computed for each outer fold and summarized as means and standard deviations to characterize fold-level variability.Reference Sokolova and Lapalme⁶⁸ In addition, pooled subject-level out-of-fold predictions were aggregated across all outer folds and repeats to derive pooled performance estimates, which served as the primary measure of overall generalization.Reference Varoquaux, Raamana, Engemann, Hoyos-Idrobo, Schwartz and Thirion⁵⁰ Receiver operating characteristic (ROC) curves and the area under the ROC curve (ROC-AUC) were calculated using these pooled out-of-fold decision scores.Reference Fawcett⁶⁹ Bootstrap 95% confidence intervals (2000 iterations) were estimated for all pooled performance metrics to quantify uncertainty.Reference Efron, Kotz and Johnson⁷⁰ Separately, classification thresholds were determined on training data by maximizing the Youden J statistic within the nested procedure to avoid reliance on a fixed default threshold.Reference Youden⁷¹ The derived thresholds were then applied to the corresponding outer test folds.

Permutation testing was performed for the model family selected by composite mean-rank aggregation across 6 pooled performance metrics (balanced accuracy, ROC-AUC, accuracy, sensitivity, specificity, and F1-score). Balanced accuracy was used as the primary metric for permutation inference, and the complete nested cross-validation procedure—including hyperparameter optimization in the inner loop—was re-executed for each of 200 random label permutations to generate its null distribution.Reference Ojala and Garriga⁷² For each permuted evaluation, decision thresholds were re-determined from the corresponding outer training partition by maximizing the Youden J statistic. A one-sided permutation p-value was computed as the proportion of permutation scores greater than or equal to the observed mean score (with a plus-one adjustment).Reference Phipson and Smyth⁴³

Ablation analyses were performed under the same nested evaluation framework to characterize the contribution of entropy feature families. Band-scale entropy, CSE, and the combined feature set were evaluated, with BDI score integrated into each feature configuration for the first and third tasks, in which BDI was included as a complementary clinical predictor in the main analysis. In addition to feature-family ablation, robustness checks were conducted to assess the stability of the results. Sensitivity to segment choice was quantified by re-evaluating the best-performing model under repeated random sampling from the pool of preprocessed, artifact-free, non-overlapping 30-s segments. This procedure was repeated for 100 iterations, and the resulting performance distribution was summarized. Finally, robustness to randomization was assessed by repeating the full analysis across 10 distinct random seeds (0–9) and aggregating performance across seeds.

Model interpretability was evaluated using Shapley Additive Explanations (SHAP) values computed for the best-performing tuned model to characterize global feature contributions.Reference Lundberg and Lee⁷³For reproducibility, all machine learning analyses described in this section—from nested cross-validation through interpretability assessment—were implemented in Python (version 3.11.14) using scikit-learn (version 1.8.0), imbalanced-learn (version 0.14.1), and SHAP (version 0.50.0). ^{Reference Lundberg and Lee73–Reference Lemaître, Nogueira and Aridas75}