Advancing Post-Traumatic Seizure Classification and Biomarker Identification: Information Decomposition Based Multimodal Fusion and Explainable Machine Learning with Missing Neuroimaging Data

Multimodal Fusion

For classification tasks utilizing multiple modality information, three major approaches are found in the literature. They are, namely: late, intermediate, and early fusion. Late fusion involves combining predictions from multiple different base learners (Base), or heterogeneous ensemble classification, and it has been shown to be an efficient way of improving predictive accuracy in machine learning (ML) (Wolpert, 1992). As a form of late fusion, stacking involves learning how to best combine the classifiers, in addition to combination using simple statistics, such as voting or averaging (Džeroski and Ženko, 2004). Džeroski et al. claimed that heterogeneous base learners, stacked with a final meta classifier, can perform comparably, if not better, to the best of the individual classifiers (Džeroski and Ženko, 2004). In the second major approach, intermediate fusion, features or information extracted from the individual modalities are aggregated, and subsequently fed to a single classifier. Since fusion in this method occurs in the feature space, it is often more interpretable and also allows the application of multivariate information theoretic fusion (Brown, 2009). The final approach, early fusion, involves representing the raw individual modality information collectively in a common space. From this aggregated information, an ML pipeline can then extract relevant features and perform the desired classification. Early fusion can be challenging to implement, since meaningful aggregation of different types of modalities (e.g. spatial and temporal) is not always straightforward (Zhang et al., 2020). Past studies involving magnetic resonance imaging (MRI) and electroencephalography (EEG) have found intermediate and late fusion to be less affected by noise, misregistration, and to have lower complexity; compared to early fusion (Lewis et al., 2007; Zhang et al., 2020).

Missing Data

While acquiring data from all modalities for each subject is ideally desirable, it is often not feasible. This creates a challenge of handling such missing modality information. Missing data has historically been dealt with imputation: using the existing data to make estimates for the missing values. Imputation techniques generally take either the form of a univariate approach that considers each missing feature/ modality in isolation across all subjects (such as mean imputation), or a multivariate approach that considers the other available features/ modalities, which could be within a subject or across all subjects (such as k-nearest neighbor, kNN imputation (Scheffer, 2002)). Following imputation, standard fusion and classification pipelines can be implemented. Aside from conventional imputation, contemporary approaches include jointly performing data imputation and self-representation learning (Liu et al., 2021), or graphical model estimators that marginalize missing attributes (Murphy, 2012).

Seizure Classification

TBI causes widespread damage to both functional and structural features. Since each imaging modality tends to pick up only certain items from this entire set of features, it is essential to consider multiple modalities for such classification (Sharp et al., 2011). Whereas many studies have investigated a unimodal strategy for seizure classification, few have involved a multimodal approach. Two such studies explored non-TBI seizure classification from EEG and structural MRI (Memarian et al., 2015), and from EEG and electrocardiography (ECG) (Yang et al., 2022). Akrami et al. examined PTE classification, utilizing functional MRI (fMRI) and T2-weighted fluid-attenuated inversion recovery (T2-FLAIR) for lesion information (Akrami et al., 2021). However, to our knowledge, no earlier work investigated the classification of LPTS from an incomplete multimodal dataset, proposed suitable novel fusion algorithms, and searched for modalityspecific potential biomarkers. Aside, the current work builds on our earlier preliminary works in preprocessing and analyzing the data (Akbar et al., 2020, 2021; Faghihpirayesh et al., 2021) from the ongoing multi-center EpiBioS4Rx cohort.

2.1.1 Acquisition

This work was approved by the UCLA Institutional Review Board (IRB#16 − 001576) and the local review boards at each EpiBioS4Rx Study Group institution. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

According to the EpiBioS4Rx protocol (Vespa et al., 2019), moderate-to-severe TBI subjects with frontal and/or temporal hemorrhagic contusion and Glasgow Coma Scale (GCS, severity of injury, lower is severer) score between 3 − 13 were eligible for enrollment. A total of 48 subjects (12 female, 36 male; age = 42.1 ± 19.3 years; GCS = 7.7 ± 4.4), enrolled in 8 different sites, were chosen for this work. Of these subjects, 17(35%) experienced at least one LPTS, whereas the remaining 31 did not experience any for the entire two-year follow-up period, which is the minimum duration needed to identify about 80 − 90% of the subjects who will eventually develop PTE (Temkin, 2009). Even though the modalities of interest for this work are dMRI, EEG, and fMRI, the acquired data from the subjects included other MRI sequences such as T1-weighted, T2-weighted, T2-FLAIR, etc. The following three sections will summarize the preprocessing steps involved in preparing our three data modalities of interest for the extraction of relevant features.

2.1.2 dMRI Preprocessing

dMRI scans corresponding to multiple diffusion gradient values and directions were collected. Acquired dMRI data was processed in the Oxford FMRIB Software Library (FSL) (Jenkinson et al., 2012), to estimate diffusion tensor imaging (DTI) parameters, such as fractional anisotropy (FA) maps, in subjectspecific spaces. Since FA features have been reported earlier to be promising in characterizing PTS (Gupta et al., 2005), we focused on extracting and utilizing FA features in this work.

The collected individual subject FA images were then transformed to a common Montreal Neurological Institute (MNI) space by registering to a standard HCP 1065 DTI FA template, following Pipeline 1 in (Akbar et al., 2021). Finally, tract-based spatial statistic (Smith et al., 2006) of each registered image was carried out using the mask and distance map obtained from the standard template, to extract mean FA values along 63 white matter (WM) tracts and bundles obtained from the JHU-DTI atlas (Zhang et al., 2010). These mean FA values are recorded as the dMRI features (x_d).

2.1.3 EEG Preprocessing

Epileptiform abnormalities (EAs) such as seizures, periodic discharges (PDs), and abnormal rhythmic delta activity (ARDA), are potential biomarkers of epileptogenesis (Vespa et al., 2016; Faghihpirayesh et al., 2021). Kim et al. (Kim et al., 2018) showed that the presence of EAs in the EEG signal during the acute period following TBI independently predicted PTE in the first year post-injury. Clinicians annotated the obtained EEG to denote the presence of three such EAs (Seizure, PD, and ARDA), which formed our EEG features (x_e).

Exclusion Criteria: Of the total 48 subjects in our study, 10 had valid EEG recordings (following review by EEG experts) for analysis. All subjects met the two-year follow-up requirement stipulated by (Temkin, 2009). Within this group, 7 subjects (representing 70% of the sample) developed LPTS, while the remaining 3 subjects (representing 30% of the sample) did not.

2.1.4 fMRI Preprocessing

To incorporate lesion information in the rs-fMRI preprocessing, we manually segmented damaged brain tissue into parenchymal contusions and brain edema with ITK-SNAP (Yushkevich et al., 2006) from the acquired 3D T2-FLAIR scans. Brain contusion was defined as a lesion with abnormal signal intensity and hemorrhagic volume > 1 ml, and brain edema was defined as a region surrounding or in the proximity of the contusion with hyperintense signal compared with the WM signal on T2-FLAIR images (Chang et al., 2006). Each segmentation was carried out by a student research assistant and reviewed for clinical accuracy by one of two medical physicians with research experience in neuroradiology (Bennett et al., 2023). For each TBI subject, we obtained a 3D lesion mask including contusion and edema. For this work, the lesion incorporation served to obtain a more accurate preprocessing pipeline, in which brain alterations related to TBI were not mislabeled as brain tissues such as WM, cerebrospinal fluid (CSF), and grey matter (GM). Therefore, we did not consider the edema and contusion separately within each 3D mask. Then, to use the lesion masks in the rs-fMRI preprocessing pipelines, we performed affine registration on each mask using the MNI 152 template with the Linear Registration Tool (FLIRT) of FSL (Jenkinson et al., 2012). Affine transformations were used because they allowed us to maintain the morphological characteristics of the lesions better while at the same time obtaining robust registrations.

The preprocessing used a modification of the unified segmentation normalization of SPM12 (Ashburner and Friston, 2005) as accessed through the CONN toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012). The literature on normalization suggests that SPM12’s unified segmentation normalization outperforms other pipelines when applied to simulated lesions (Crinion et al., 2007). However, further work (Andersen et al., 2010) showed that to achieve the best performance, information about the lesion should be included, motivating the inclusion of the lesion mask in preprocessing.

The fMRI preprocessing pipeline incorporates information about the lesion by performing SPM12’s unified segmentation normalization using a lesion-modified tissue probability map (TPM). As discussed in (Andersen et al., 2010), modifying the TPM performs implicit cost function normalization because it ultimately causes SPM to ignore the lesion areas during normalization. The TPM normally provides a prior probability of a given voxel belonging to one of 6 tissue classes; GM, WM, CSF, skull, soft tissue, and other (Ashburner, 2009). SPM12’s default TPM was augmented for each subject with a 7th tissue class corresponding to the MNI lesion mask for that subject, thereby setting the prior probability of GM, WM, or CSF in the lesion areas to 0. The modified TPM was created by a new function in the CONN toolbox, conn createtpm. Finally, using the AAL3 atlas (Rolls et al., 2020), the means of the positive (x_{f −p}) and negative (x_{f −n}) parcel-to-parcel correlations (Faria et al., 2012), and the mean of the lesion volume overlap with parcellations (x_{f −o}) were recorded as the fMRI features.

Exclusion Criteria: Of the total 48 subjects in our study, 36 had valid fMRI scans for analysis. From these 36 subjects, we further excluded 4 as they did not meet the two-year followup requirement stipulated by (Temkin, 2009). This left us with a final sample of 32 subjects. Within this group, 12 subjects (representing 38% of the sample) developed LPTS, while the remaining 20 subjects (representing 62% of the sample) did not.

2.2.1 Data Partitions

While the goal of the study was to acquire all three modalities of interest (x_d, x_e, x _f) for each subject, all subjects did not have all modalities. Table 1 shows the number of subjects for each combination of modalities, as well as the relative number of subjects within each modality configuration. Out of the included total N = 48 subjects, only seven had information available from all three modalities. Of the other 41 subjects, 21 had a single missing modality, while the remaining 20 had two missing modalities.

From the total N subjects, we then created train and test partitions, following a five-fold cross-validation (CV). Any missing modalities were handled using methods described in Section 2.2.3. We concatenated all D features (166 × 3 for the three types of fMRI, 63 for dMRI, 3 for EEG) along the columns of , for I ∈ n samples in each CV train fold. Each , representing an individual subject, is then stacked vertically as rows to give the training partition (X_train) for each CV round Each test partition (X_test) is constructed in a similar manner with m = N − n samples. The complete feature set, including the train and the test sets, is denoted by [X = X_train, X_test]^T, where X ∈ ℝ^N×D. Similarly, the train-test partitions for the complete set of LPTS labels y ∈ ℝ^N×1, corresponding to any CV round, could be expressed as y = [y_train, y_test]^T.

2.2.2 Classification Framework

The comprehensive modeling framework employed for the LPTS classification task is illustrated in Fig. 1. To begin, highdimensional features drawn from diverse modalities, symbolized as x_d, x_e, x _f, are concatenated. These combined features then pass through an imputation stage, which estimates missing data points using the available data. The data, now filled in, proceed to the dimension reduction stage, which projects the incoming high-dimensional data into a more tractable lowerdimensional space. This reduced dataset is subsequently fit to a classification stage that predicts the binary seizure label, denoted as ŷ. An all-encompassing fusion stage determines the fusion technique to be utilized at any pertinent stage of the pipeline. The (Yes/ No) choices indicate whether a stage is active or not: if ‘yes’, the stage employs one of several techniques discussed in subsequent sections; if ‘no’, the stage becomes a passive conduit, with its output remaining identical to its input.

• Download figure
• Open in new tab

Fig. 1.

The overall modeling framework for the seizure classification task. High-dimensional features from the three different input modalities (x_d, x_e, x _f) are first concatenated, and then passed on to the first stage. The data is processed in each stage, and proceeds along in the direction of the arrows. The final output on the right is the predicted binary seizure label (ŷ). The (Yes/ No) choices dictate whether a certain stage is active or not: if yes, it uses one of the several techniques discussed later; if no, it becomes transparent where the output is the same as the input.

2.2.3 Imputation Techniques

Most classification techniques documented in the literature are not equipped to handle missing data. To establish baselines for individual modalities and apply several relevant classification pipelines, we explored four standard techniques for imputing missing modality features: mean, median, k-nearest neighbors (kNN), and iterative imputation.

The mean and median imputation strategies substitute each missing feature in with the mean and median, respectively, of the non-missing values in the corresponding feature column of X_train. For kNN imputation, each missing feature in is replaced with the mean value of the feature found in a specified number of nearest neighbors s_ne (based on Euclidean distance among the non-missing features) of within X_train (Troyanskaya et al., 2001). Iterative imputation, on the other hand, fills the missing values in by modeling the d features with missing values as a function of the other D−d features in a roundrobin fashion, for a predefined number of iterations. For all the imputation techniques utilized, we extracted information solely from X_train (in an unsupervised learning environment), and not from y_train.

2.2.4 Dimensionality Reduction and Classifiers

When considering X_train as input, the feature dimension D is considerably larger than n. Aside, by the inclusion of all features, we risk hurting the final classification performance by adding noise and overfitting by increasing complexity (Kohavi and John, 1997). Consequently, it is advisable to reduce the dimensionality of such high-dimensional input. As such, we cascaded a standard dimension reduction technique at the start of each classifier for all fusion approaches, except for the intermediate fusion techniques (Other) outlined in Sections 2.6.4-2.6.7, which inherently serve the purpose of dimension reduction. The three standard reduction procedures adopted in this work are: principal component analysis (PCA) (Hotelling, 1992), and Kbest features chosen by maximizing either the χ²-test score or the Fischer score (f_score), between X_train and y_train (Pe-dregosa et al., 2011).

For all classification tasks, either as an interim individual learner or as the final estimator, we evaluated both SVM and AdaBoost classifiers. The choice of these two classifiers aligns with previous studies, which have found both SVM (Jie et al., 2015) and AdaBoost (Akbar et al., 2021) to be suitable for classification tasks involving neuroimaging features. For SVM, we explored both linear and radial basis function (RBF) kernel variants.

The top K or k features/components to be chosen, for either one of the standard techniques or the intermediate fusion techniques, were varied in the search space as shown in Fig. 2. The implementation difference lies in the selection of different values for K per CV fold for the standard reduction techniques, following a nested 5-fold CV within the training set, whereas the same value for k was utilized (for computational considerations) for all outer CV folds involving the intermediate reduction techniques. For the classifiers, the search space for SVM RBF’s gamma (γ), L-2 regularization (𝓁₂), and the number of trees (tree) in AdaBoost can also be seen in Fig. 2. All models were trained and tested in Python 3.7.4, with support from scikit-learn (Pedregosa et al., 2011) version 1.0.2, CCA-Zoo(Chapman and Wang, 2021), and MMIDimReduction (Özdenizci and Erdoğmuş, 2021). The source code for this work is open-source and available to the public. ¹

• Download figure
• Open in new tab

Fig. 2.

Model pipeline showing the different dimension reduction techniques and the classifiers, along with their parameter search spaces. The Other block within dimension reduction contains our extension of an existing technique (Section 2.3), as well as our two proposed algorithms (Sections 2.6.6 and 2.6.7).

2.6.1 Feature Union

In this approach, all D features from the three modalities were concatenated to give X_train ∈ ℝ^n×D. Subsequently, X_train was passed through a standard reduction technique to extract K features, and then fitted to a classifier.

2.6.2 Sequential Feature Selection (SFS)

For the selection of features from such high D dimensional data, standard univariate feature selection techniques ignore the mutual information among features. Multivariate feature selection techniques, however, consider a group of features in its entirety. Unfortunately, searching for the globally optimal subset with exhaustive search is O(2ⁿ), and can be computationally intractable. As a result, it is common practice to resort to algorithms striving to obtain a locally optimal, but perhaps globally sub-optimal, feature set with a lower complexity (Dash and Liu, 1997). The wrapper method, a multivariate feature selection technique, attempts to find a set of highly important features by fitting a particular classifier on all features in a nested CV within the train set (Erdogmus et al., 2007). Wrapper methods have been demonstrated to provide superior performance (compared to filter methods in feature selection) in classification tasks involving high-dimensional neuroimaging data with limited samples (Jie et al., 2015).

In SFS, a multivariate wrapper technique, we greedily choose features for classification, either by forward selection or by backward elimination. For the first round in forward selection, based on the score of the wrapper classifier, the most informative feature is added. In the case of backward selection, the least informative feature is removed. Forward selection thus involves a bottom-up search strategy, which begins with an empty set, and during each iteration, a new feature is added to the current set so the loss function is reduced (Kohavi and John, 1997). On the contrary, backward elimination, follows a top-down approach, starting with the complete set D, features are removed one at a time such that the reduction in performance is kept at a minimum (Kohavi and John, 1997). Choosing either forward or backward selection, we retain k_SFS features. Selected features are collected as , and subsequently fitted to a final classifier. SFS has an average complexity of O(n²).

2.6.3 Stochastic Mutual Information Gradient (SMIG)

In this method, we aim to learn a feature transformation network such that the high n × D-dimensional input feature space is mapped to a lower n × k_SMIG-dimensional transformed feature space (Özdenizci and Erdoğmuş, 2021). This mapping is done while maximizing the mutual information between the transformed set and the train labels y_train as where I(·) is the information/ mutual information function, Z_SMIG contains the transformed training samples, and Ω denotes the function space for possible feature mappings ϕ. The parameters of this ϕ^⋆(·) network is updated iteratively by a technique known as SMIG (Özdenizci and Erdoğmuş, 2021).

2.6.4 Canonical Correlation Analysis (CCA)

We hypothesized that by looking in the directions of the maximal correlations of the input modalities, we would be able to extract information about their shared space, which might help the classification task, since the y was assumed to be the underlying common latent state in Fig. 3. In CCA, two sets of variables are linearly combined to obtain canonical variates such that the correlation between the linear combinations is maximized (Hotelling, 1992). This technique was thus utilized to combine our input modalities in a three-fold combination: (x_d, x_e); (x_d, x _f); and (x_e, x _f).

For mathematical formulation, let us assume the combination (x_d, x _f). If we can express two linear functions u = a^Tx _f and v = b^Tx_d, and if the covariance (Hotelling, 1992) it can be shown that their correlation is maximized when where the eigenvalues λ ∈ { ρ₁≥ … ≥ ρ_s+t } are in descending order. The maximum correlation was obtained with λ = ρ₁, and it provided our first pair of variates and . From the obtained D such pair of variates, we then selected k_CCA pairs to form , and fitted Z_CCA to a final classifier.

2.6.5 Generalized Canonical Correlation Analysis (GCCA)

CCA can only linearly combine two sets of variables. GCCA, however, can combine more than two variables (Tenenhaus and Tenenhaus, 2011). In our case, GCCA combined (x_d, x_e, x _f) directly. GCCA attempts to solve the following optimization problem where corr(·) is the Pearson correlation function, g(·) is the identity function, c _jq = 1 represents the connections between modality feature matrices { x_d, X_e, X _f }, and J=3.

Similar to CCA, we then selected k_GCCA pairs of variates from the total D pairs, to give , and fitted Z_GCCA to a final classifier.

Algorithm 1:

RECC for a given selection algorithm

• Download figure
• Open in new tab

2.6.6 Recursive Elimination of Correlated Components (RECC)

2.6.7 Information Decomposition and Selective Fusion (IDSF)

Background

Our earlier hypothesis in Section 2.6.4 was that by utilizing information from the shared space among the

Algorithm 2:

IDSF for one decomposition algorithm and one selection algorithm prepared by RECC

• Download figure
• Open in new tab

modalities, we would hope to observe a better classification performance, since the LPTS label y was also assumed to be a shared latent state. However, the possibility exists that a considerable fraction of the information, useful for classification, is present uniquely in a space not shared by the individual modalities (Williams and Beer, 2010). This is illustrated in Fig.4. The total information provided by two sample modalities x₁ and x₂ about y may be denoted by the mutual information measure I(y; x₁, x₂). However, it is also important to understand how exactly x₁ and x₂ contribute to this aggregate information. Three primary scenarios may arise (Williams and Beer, 2010) when examining the individual contributions of x₁ and x₂. Firstly, x₁ and x₂ may offer common shared information as

where ∩ is the intersection operator. Secondly, x₁ might offer unique information, Unique(), not provided by x₂, or vice versa. An example of the former case is Finally, the combination of x₁ and x₂ may provide synergistic information that is not available from either alone, but only obtained when both are considered together (e.g. an XOR gate for binary inputs)

• Download figure
• Open in new tab

Fig. 4.

Structure of multivariate information that two sample modalities x₁ and x₂ would likely provide about the LPTS label y.

Feature Contribution

For interpreting the top models, the SHapley Additive exPlanations (SHAP) framework is used. Introduced by (Lundberg and Lee, 2017), it builds upon the concept of Shapley values from cooperative game theory. By using SHAP values, one can fairly distribute the contribution of each feature to a prediction for a specific instance. These values enable a better understanding of the inner workings of complex ML models, such as decision trees, SVMs, and even deep learning networks. Mathematically, the SHAP value ϕ _j of feature j is computed as the average marginal contribution of the feature across all possible combinations of features: where F is the set of all features, S indicates all feature subsets obtained from F that do not include feature j, and f is the ML model. We use the kernel method of SHAP calculation (Lundberg and Lee, 2017).

Imputation

Based on our preliminary classification results with the four different imputation techniques, kNN imputation outperformed the other three. For this reason, all subsequent results in this section, wherever data was imputed, was done using the kNN imputer, with s_ne = 1.

NB-extend

The comparison of the mean AUC and the mean weighted F1 performances, on the test sets from all five CV folds, for the different fusion techniques (except IDSF) are summarized in Table 2. The first row shows the classification performance of our extended NB-extend estimator with late fusion, where the original input features (no imputation) were used to fit individual modality base learners. NB-extend (AUC=0.710) outperformed the individual modality base learners (AUC= {0.664,0.670,0.619}), using imputed features without any fusion, and which form the next three rows of Table 2.

Late Fusion

Of the techniques using imputed features and fusion, we first list the performance of the soft late fusion estimator. This probabilistic estimator also performed better (AUC=0.756) than the individual modality base learners. It even performed better than both the late fusion NB-extend estimator (no imputation), and the intermediate fusion (union of all modality features) classifier (AUC=0.736), in the subsequent row of Table 2.

Intermediate Fusion

The last six rows of Table 2 correspond to the dimension reduction techniques involving intermediate fusion, as discussed earlier in Sections 2.6.4-2.6.6. Among them, the first four rows presented existing techniques in literature, and the best performance (AUC=0.784) was achieved with CCA. Even though GCCA took into account all three modalities (x_d, x_e, x _{f −n}), compared to CCA’s two (x_d, x _{f −p}), we noticed it lagged in performance (AUC=0.654) substantially compared to that of CCA. Of the two SFS implementations discussed (forward and backward), forward selection performed better (AUC=0.676), and was thus recorded. Similarly, between the two implementations of SMIG (linear and non-linear), linear performed better (AUC=0.699) and was recorded. Finally, the last two rows in Table 2 demonstrated the improvement in performances with our proposed RECC (AUC= { 0.710, 0.753}), when applied on the original SFS and SMIG implementations, respectively. The best F1 score (0.753) was obtained with RECC-SMIG, indicating a potential choice for the best model, if F1 score is more important in the model selection.

Proposed RECC

To obtain a better understanding of how the choice of ρ_thresh affected the AUC performances of the RECC combinations corresponding to the last two rows in Table 2, we plotted AUC vs. ρ_thresh in Fig. 5. The standard implementations of SFS and SMIG corresponded to ρ_thresh = 1. As we lowered ρ_thresh, we noticed the AUC performances go up until they peaked at ρ_thresh = 0.5 for both RECC implementations, and then started to go down. At ρ_thresh = 0.15, the performance of RECC-SMIG almost returned to its SMIG baseline, whereas that of RECC-SFS remained higher than its SFS baseline.

• Download figure
• Open in new tab

Fig. 5.

AUC performance of RECC-SFS and RECC-SMIG combinations from the last two rows in Table 2, as a variation of the Pearson coefficient threshold ρ_thresh, with SVM classifier and modalities x_d, x_e, x _{f −o}.

Proposed IDSF

In Table 3, we recorded the best AUC performances for each combination of the two decomposition algorithms (CCA and GCCA) with the two selection algorithms (SFS and SMIG) prepared by RECC, using intermediate fusion with IDSF. We found the best IDSF performance (AUC=0.792), also the best among all tested models, was obtained with CCA and RECC-SFS. Just as seen earlier with CCA and GCCA, IDSF with CCA was superior to IDSF with GCCA. Interestingly, while individually both SMIG and RECC-SMIG performed considerably better than SFS and RECC-SFS earlier, their IDSF performances were closely matched. Aside, while IDSF improved on the performances of the decomposition algorithms on all occasions, it did not necessarily do so over that of the selection algorithms (e.g. IDSF with GCCA and RECC-SMIG). Finally, the best F1 score (0.677) for IDSF was obtained with GCCA and RECC-SFS.

ROC Analysis

To examine the effect of the proposed algorithms RECC and IDSF in terms of operating points, mean ROC curves (from the 5-fold CV) are plotted in Fig. 6 for the best model from Table 3, as well as for the individual techniques that went into the best model from Table 2. The variation of the best ROC curve (IDSF) is estimated by a binomial distribution approximation, for a 95% level of confidence (Sourati et al., 2015). We noticed that just like RECC-SFS improved on standard SFS, IDSF (AUC=0.79) improved on the performance over its two input techniques: CCA (AUC=0.78) and RECC-SFS (AUC=0.71). It also obtained the best operating point (0.01, 0.67: shown in a red circle), according to the highest Youden index (0.68) and the lowest d_UL (0.23).

SHAP Analysis

To delve deeper into the role of individual features in LPTS prediction, we employed the SHAP toolbox (Lundberg and Lee, 2017), fitted with the Xtrain from each fold to interpret the Xtest for the corresponding fold. The leading models without and with imputation are NB-extend and IDSF (CCA-RECC-SFS), as earlier seen in Table 4. It is worth noting that SHAP values cannot be computed for missing entries in NB-extend, and IDSF only offers access to a limited number of features per fold. To establish a middle ground for easy comparison between these models, SHAP values were computed on the imputed data (used by IDSF), utilizing the AdaBoost classifier (same configuration as NB-extend). The resultant SHAP values for all the five folds (the entire dataset with x_d, x_e, x _{f −o}) are depicted in Fig.7. The mean SHAP values illustrated in Fig.7(a) represent the overall influence of individ-ual features, whereas the swarm plot in Fig.7(b) demonstrates the effects of the value of individual feature instances on predicting the binary LPTS labels. From Fig.7(a), it is observed that ALIC-L (anterior limb of internal capsule left) from x_d had the most profound impact and was also previously identified as a potential biomarker in Table 4. Fig.7(b) further reveals that higher values of ALIC-L tended to predict no LPTS, while lower values predicted LPTS. Cluster 90 from x _{f −o} had the second most significant effect on the seizure outcome prediction task, as it previously possessed the smallest unadjusted p-value among all tested features in Table4. Furthermore, Fig. 7(b) indicates that higher values of cluster 90 suggest potential LPTS, while lower values suggest otherwise.