Recognition of Electroencephalography-Related Features of Neuronal Network Organization in Patients With Schizophrenia Using the Generalized Choquet Integrals

In this study, we focused on the verification of suitable aggregation operators enabling accurate differentiation of selected neurophysiological features extracted from resting-state electroencephalographic recordings of patients who were diagnosed with schizophrenia (SZ) or healthy controls (HC). We built the Choquet integral-based operators using traditional classification results as an input to the procedure of establishing the fuzzy measure densities. The dataset applied in the study was a collection of variables characterizing the organization of the neural networks computed using the minimum spanning tree (MST) algorithms obtained from signal-spaced functional connectivity indicators and calculated separately for predefined frequency bands using classical linear Granger causality (GC) measure. In the series of numerical experiments, we reported the results of classification obtained using numerous generalizations of the Choquet integral and other aggregation functions, which were tested to find the most appropriate ones. The obtained results demonstrate that the classification accuracy can be increased by 1.81% using the extended versions of the Choquet integral called in the literature, namely, generalized Choquet integral or pre-aggregation operators.


INTRODUCTION
Mental illnesses are usually long-lasting conditions associated with great psychological suffering, the substantially limited possibility of independent functioning, and social development. Among them, schizophrenia (SZ) is one of the most severe forms of mental health disorder with the complex and multidimensional clinical picture. The onset of SZ occurs most often in adolescence or early adulthood commonly has a slow and hidden course consisting of gradual augmenting of the so-called negative syndromes, i.e., loss of interests, affective blunting, reduced initiative, and social isolation, and more or less delayed phase of active psychotic exacerbation characterized by the presence of delusions, i.e., incorrect judgments of reality and the behavior of other people, as well as hallucinations, i.e., incorrect sensory impressions, most often in the auditory form (Rosen et al., 1984;Heilbronner et al., 2016). In addition to the negative and positive symptoms, there are also various cognitive disorders including disturbances in the course of thinking and deficits in specific cognitive domains, such as attention, memory, cognitive speed, language, and communication, and difficulties with adapting to new circumstances and problem-solving (Szöke et al., 2008;Krukow et al., 2017;Green et al., 2019). It should be noted that long-term pharmacological treatment of the disease is the main form of therapeutic intervention focused mainly on psychotic syndromes. However, even when modern methods of treatment are applied, distortions of cognitive processes improve to a much lesser extent, often causing lifelong constraints in achieving full independence (Keefe, 2019). Personal, social, and economic burdens associated with severe mental illness prompt researchers to search for new therapies and also to develop accurate methods of differential diagnosis, which should be ultimately based on objective, biological markers (Pantelis et al., 2009). The development of new neuroimaging techniques enables researchers to identify neural circuits that underline the human brain integration system. Various neuropsychiatric conditions are correlated with changes in brain communication patterns and pointed as potentially useful biomarkers for clinical applications (Sporns et al., 2005). In accordance with the results of earlier studies focused on brain synchronization, SZ is seen unequivocally as a disconnectivity disorder characterized by abnormal functional and structural connectivity of the brain (Friston and Frith, 1995). Application of diffusion tensor imaging (DTI) methods, such as magnetic resonance imaging (MRI) technique, into SZ research, showed disconnection and multiple microstructural aberrances of brain white matter fibers (Zalesky et al., 2011;Klauser et al., 2017). Studies based on electroencephalography (EEG) and functional MRI (fMRI) also revealed abnormalities in the functional connectivity of the brain, which were also correlated with the clinical picture of the SZ (Skudlarski et al., 2010;Uhlhaas, 2013;Krukow et al., 2018). Nevertheless, to understand the systemic level of the brain organization and to explain neurophysiological processes such as disconnectivity syndrome in the SZ, researchers started to analyze the brain as a complex network (van den Heuvel and Sporns, 2013). The neural network is understood as a system of spatial (anatomical) and temporal (synchronous firing of neuronal assemblies) dimensions, involving different brain regions interconnected with each other (Zalesky et al., 2010). However, to analyze the state of the functional and structural connections from the viewpoint of the entire brain, an infinite number of potential anatomical and functional interactions between a given set of neural regions makes such an analysis a challenge almost impossible to obtain. Therefore, a graph theorem has been introduced to solve this problem and to test the complex whole-brain networks in their global dimension (Van Den Heuvel and Fornito, 2014). Previous studies investigating the neural brain networks in SZ showed significantly changed network organization as indicated by graph-analytical measures of global, short communication paths (Yan et al., 2015), local organization (Alexander-Bloch et al., 2010), and small-worldness (balance between local segregation and global integration) (Shim et al., 2014). Aberrant functional networks in the SZ were also linked with cognitive impairments (Sheffield et al., 2015;Krukow et al., 2020) and the duration of the illness (Jonak et al., 2019).
Previous studies considered the problem of automated classification of altered brain activity in SZ based on the EEG or fMRI data. Among traditional classifiers, methods such as support vector machine (SVM; Shim et al., 2016;Liu et al., 2017;Huang et al., 2018), adaptive boosting (Sabeti et al., 2011), kernel discriminant analysis (KDA; Zhu et al., 2018), or nearest neighbor algorithm (Parvinnia et al., 2014) are used. Some of these studies (Sabeti et al., 2011;Parvinnia et al., 2014) applied time-frequency features obtained from single EEG channels, which is a limited capacity approach as it does not consider interactions between channels understood as a network. Other authors applied a convolutional neural network (Phang et al., 2019) and deep neural networks (DNNs; Plis et al., 2014;Guo et al., 2017). In addition, manifold learning for aggregation was considered in works by Shen et al. (2010); Anderson and Cohen (2013), and Gallos et al. (2021a,b). The idea of applying fuzzy classification into SZ-based data is a relatively new concept, as there are only a few papers on this subject (Sabeti et al., 2007;Silvana et al., 2018). One of the answers to the problems related to the application of single classifiers in the processes of automated disease diagnosis may be using various aggregation models. Aggregation can be carried out at the stage of data analysis in the form of information fusion and the stage of analysis of classification results. Despite some shortcomings such as extending the duration of the diagnosis process or the need to implement additional algorithms, the undoubted advantage of this approach is the increase in the effectiveness of classification, which, combined with the field of application critical to human health, is of key importance. Common examples of aggregation operators are voting, maximum, minimum, and median functions. The methods based on triangular norms (Klement et al., 2000) or ordered weighted averaging operators (OWA; Yager and Kacprzyk, 2012) are somewhat more complex. Various general approaches to the aggregation of classifiers were already presented (e.g., in publications of Alsina et al., 2006;Beliakov et al., 2007;Grabisch et al., 2009;Calvo et al., 2012;Gągolewski, 2015;Dolecki et al., 2016;Baczyński et al., 2017). Recently, one of the dominant techniques is using the Choquet integral or its generalizations or extensions Pedrycz, 2004, 2005;Karczmarek et al., 2014Karczmarek et al., , 2017aKarczmarek et al., ,b, 2018Karczmarek et al., , 2019bAnderson et al., 2018;Rutkowska et al., 2020). In particular, recent studies on the so-called pre-aggregation functions offer hope for the development of this approach (Lucca et al., 2015Bustince et al., 2016;Dimuro et al., 2017;Dias et al., 2018). They are particularly used in computer image analysis and its subdiscipline of facial recognition (Karczmarek, 2018;Karczmarek et al., 2019a). Detailed theoretical and practical analyses of the approach based on pre-aggregation functions, i.e., slightly weakening the classical aggregation (Beliakov et al., 2007) conditions, are still ongoing. Nevertheless, the weakening of these conditions does not have a negative impact on the classification results, which is confirmed by the experimental outcomes from the above-mentioned studies. A survey of the generalizations of Choquet integral can be found in Dimuro et al. (2020).
The problem undertaken in this study was related to the effective automatic distinction between patients diagnosed with SZ and healthy subjects based on EEG-based features of neuronal network organization. The main goal of this study was to find the appropriate operator aggregating the neurophysiological outcomes and categorizing them as patients diagnosed with SZ or healthy controls (HC), i.e., increasing the effectiveness of the classification. For this purpose, various generalizations of the Choquet integral were tested and a set of over a thousand aggregating functions not related to the Choquet integral was verified. In the section on numerical experiments, we indicate the classes of functions and their detailed parameters that work best in terms of the identification of SZ. The dataset applied in this study included data gathered from 40 subjects, i.e., 20 schizophrenic patients and 20 HCs. The Granger causality (GC; Granger, 1969) concept had been applied to particular EEG bands to achieve functional brain connectivity measures. The collected measurements were analyzed using a minimum spanning tree (MST; Stam et al., 2014). The global MST parameters obtained in the analysis were chosen as features in the classified dataset. Applying the MST algorithms enabled grasping the backbone structure of the brain network with only the strongest connections included (González et al., 2016;Van Dellen et al., 2016). Using MST ensured that the link between nodes was not based on an arbitrarily set connectivity strength threshold, which allowed avoiding the bias in network density computations (Tewarie et al., 2015). In general, the MST parameters were chosen because this method lacks some theoretical and mathematical problems incorporated to more typical network organization indicators based on the smallworldness approach and, above all, because the authors wanted to refer to the concept of SZ as a disconnection disease in which pathology of neuronal integration is not isolated only to selected regions or type of synchronizations but has a global dimension. The MST enables the characterization of the global, whole-brain network.

Participants
Twenty patients, who met the DSM-5 (Structured Clinical Interview for DSM-5) criteria for SZ, were involved from the Department of Psychiatry, Psychotherapy and Early Intervention, Medical University of Lublin. Additionally, the other criteria were as follows: age over 18 years; minimum 10 years of regular education; not more than 5 psychiatric hospitalizations associated with exacerbation of psychosis; no markers of structural brain abnormalities visualized on MRI, indicative of surviving craniocerebral trauma or neurovascular episodes; and lack of serious somatic diseases needing intense pharmacotherapy that would impact the EEG recordings. During testing, all patients were on stable doses of atypical antipsychotics. Using anticholinergic agents, benzodiazepines, and mood stabilizers up to 3 months before the assessment was an exclusionary factor for all participants. The patients participated in the study during the last week of psychiatric hospitalization, after obtaining a significant clinical and functional improvement, being fully able to give consent and undergo EEG examination. The control group consisted of HC, demographically matched to the clinical group, who were chosen from the local community. Additionally, HC had no history of psychiatric diagnoses, per the Structured Clinical Interview for DSM-5, brain disease, or neurological injury as well as no family history of psychosis. All patients consented to the study in accordance with the protocol approved by the Bioethical Commission of the Medical University of Lublin. The Commission also validated the methods used in the study. Demographical and clinical data are presented in Table 1. The groups did not differ significantly in terms of age (SZ = 34.41, SD = 8.41; HC = 31.63, SD = 6.42), number of the years of education (SZ = 12.43, SD = 2.94; HC = 14.87, SD = 1.68), and gender (SZ = 50% of men; HC = 50% of men). In the SZ group, the duration of illness lasted for about 12 years.

Data Acquisition
Using a 21-scalp position, electro-cap electroencephalograph (Electro-Cap International Inc., OH, United States) and Ag/AgCl disk electrodes, in 10 min of resting-state, EEG data were recorded for each participant. Electrodes were distributed according to the 10-20 International system (Fp1, Fp2, F3, F4, C3, C4, P3, P4, O1, O2, A1, A2, F7, F8, T3, T4, T5, T6, Fz, Pz, and, Cz). Subjects were seated with eyes closed and restricted head movement. The electrode impedances were kept below 5 k, and the data were filtered from 0.5 to 70 Hz (with active notch filter set at 50 Hz) when the sampling rate was 512 Hz. The data were exported into ASCII format after recording. Post-processing procedures were carried out in the EEGLAB program, which is a MATLAB toolbox. First, the signal was filtered using the bandpass filter at 0.5-45 Hz (second-order Butterworth filter). Second, the reference was changed offline into the averaged. Next, from the processed signal, 25 epochs lasting for 8 s (4,096 samples) without artifacts were extracted for each patient by a clinical neurophysiologist. Last, EEG signals were divided into six

Data Processing
Several steps were involved in the data processing procedure. Feature extraction was associated with the calculation of functional brain connectivity (FC) measures separately for particular electrodes and each EEG frequency band. The FC was calculated to indicate the statistical dependence between the spatially distributed neurophysiological time series such as EEG signals stemming from separate units of a nervous system (Cheng et al., 2015). There were several metrics used in assessing FC strength, such as classical measures (e.g., Pearson's correlation coefficient, cross-correlation function, or coherence), phase synchronization indexes (e.g., phase lag index or phaselocking value), and GC measures. We chose the classical linear GC. The idea of GC (Granger, 1969) was based on the assumption that having two simultaneously determined signals (X and Y), the signal X could be better explained using information from the signal Y than using only information from the signal X.
In such a situation, signal Y could be specified as "causal" to signal X. The GC measure is widely applied as a statistical tool to detect the influence of particular system components (Nolte et al., 2010). For the GC calculation, we used the Matlab MVGC toolkit (Sackler Centre for Consciousness Science, University of Sussex, Brighton, United Kingdom; Barnett and Seth, 2014), which was based on advanced VAR (vector autoregressive) model theory. To optimize auto-covariance delays, the Akaike information criterion was used to estimate the optimal model order. MVGC algorithms were used to convert EEG signals into auto-covariance data. The observed auto-covariance sequences were then subjected to paired spectral GC. In the case study, the calculated GC measures were used in the next feature extraction procedure step consisting of deriving MST. The MST was built based on Kruskal's algorithm. First, the weight of all edges was sorted, and then the stronger edges were connected, i.e., those with the highest connectivity values, eliminating those that form loops. These stages were repeated as many times as necessary until the final tree had 19 nodes and 18 edges, which corresponded to the total number of electrodes used in our EEG recording. The MST metrics were generated separately for each frequency band. Global MST parameters were used as the features for the classification procedure, which included the following: • maximal degree, a maximal degree in the MST tree, • maximal BC, maximal betweenness centrality in the MST tree, • leaf fraction (Lf), the ratio between leaf vertex number (further denoted as L) and the total vertex number, • diameter (d), the longest distance between any two vertices in the MST tree, • hierarchy (Th), the measure describing the optimality of the tree topology.
After the feature extraction procedure, the classification was done using classical classifiers.

A General Processing Scheme
The classification procedure was performed to assign observations into one of two classes: schizophrenic patient and HC. Several classical classification methods were used: cubic SVM, linear SVM, decision tree, logistic regression, multilayer perceptron (MLP), random forest, and k-nearest neighbors (kNN). Data were split into training and testing datasets based on the 80:20 ratio.
The results of this classification process, in the form of the probabilities of belonging to considered classes, were taken as the inputs to the aggregation functions. The next step of the analysis procedure was to generate fuzzy measure density values to apply aggregation operators. The aggregation operators allow combining the predictions of multiple classifiers to further improve the results. The fuzzy measures can be interpreted as the degree of trust (weights or level of importance) to predictions of the individual classifier.
In general, there are several methods of fuzzy measure generation, such as expert assumption, optimization, and the heuristic one. In our study, the cross validation-based heuristic one was applied. N-fold cross-validation was run on the training set to obtain a density measure for a classifier.
It is a well-known fact that the results of different classifications can be aggregated. This situation can be easily illustrated by the example of various sports competitions held in the form of a Grand Prix cycle, where the points are added together to determine the final winner. EEG signal-based features can be similarly aggregated. In general, the results obtained using different kinds of classifiers can be added, averaged, or, generally speaking, transformed with the help of different kinds of aggregation operators. Typical operators are median, minimum, and maximum functions, etc. The general scheme of classification using aggregation methods is presented in Figure 1. It is worth noting that the values that are input to the aggregation operator can change the distances between the training and testing representation of an EEG signal in the case of k-nearest neighborbased classifiers, likelihoods of belonging to a class in the case of neural network-like methods, etc. Here, it is worth stressing that the weights presented in this diagram can be obtained from experts but also based on the quality of classification of individual classifiers (e.g., their accuracy measures).

Aggregation of Classifiers Using the Choquet Integral and Its Extensions
One of the best known and most efficient classifiers is the Choquet integral. Hence, let us recall the main properties of the fuzzy measure, Choquet integral, and its generalizations. Let us denote a set as X. Then P(X) = 2 X is a family of all its subsets. 2 X is a σ-algebra; i.e., the empty set belongs to it, the complement of a set belonging to σ-algebra belongs to it, and the sum of countable many sets from the σ-algebra also belongs to it. Generally speaking, in the context of classification tasks, the elements of the set X are the individual classifiers (methods, parts of the images under consideration, etc.). In the context of this specific study, they are particular classifiers, see the experimental section for details of the methods. These classifiers are denoted as x 1 , x n , n 1. Now, one can define (Sugeno, 1974) a fuzzy measure as a set function g : P(X) → R satisfying the following conditions: g (∅) = 0 and g (X) = 1 (1) Here, {U n }, n = 1, 2, . . . means increasing set sequence. Recall that Sugeno λ-fuzzy measure realizes the above conditions and with λ >− 1. Here, U and W are not overlapping. In addition, we have for U i = {x 1 , ..., x i }, U i + 1 = {x 1 , ..., x i + 1 }. The following notation is used commonly: g i = g({x i }), i = 1, . . ., n. Now, let us introduce a function h(x) and let the series h(x i ), i = 1, . . ., n, be ordered non-increasingly and let h(x n+1 ) = 0. In the context of this study, the function h(·) represents a value of classifier describing the probability of belonging to a specific class. Next, the U i set is, in fact, only an abstract object. The real importance has the value of g (U i ) appearing in (5) which can be easily found recursively starting from the values of g i . The value g i represents a significance (or importance) of a particular classifier x i . Its value can be commonly defined twofold: (1) based on the opinions of experts and (2) based on initial tests. In this study, we applied the second method. Finally, n is a number of classifiers. The last parameter to be found is λ, which can be obtained from the following equation: see Sugeno (1974). For such assumptions, the Choquet integral is defined as From this function, many generalizations and extensions can be delivered as follows: and for any t-norm M (·,·), see Lucca et al. (2014), (Lucca et al., 2014(Lucca et al., , 2015, ) see , C Min (Lucca et al., 2015), where the role of the function M is played by the minimum, or C O [see ] with a so-called overlap function under the integral sign. Newer functions were proposed in Karczmarek (2018) and Karczmarek et al. (2019a). They are and the integrals inspired by some numerical analysis formulae such as and It is worth noting that M(·,·) can be any triangular norm which, as an intersection or conjunction operator in many application areas, is a counterpart to a classic product operator appearing in the original Choquet integral.

Individual Classifiers
In this study, we described particular classifiers that were considered in the series of numerical experiments and determined their accuracy. We applied the following classical machine learning models: SVM with linear and cubic kernels, logistic regression, kNN, decision tree, random forest, and MLP. The classical machine learning models were used due to a low number of observations available for training and testing. In order to obtain the fuzzy density that is necessary for aggregation, the following approach can be adapted. According to the holdout validation procedure, the data were split into training and validation subsets, where 20% of the dataset was used for validation. The fuzzy density was calculated as a mean accuracy measure obtained in the process of a fivefold cross-validation run using the training data. The resulting classification quality of separate models was tested on the validation subset after fitting the models on the complete training set. The classification accuracy values obtained with separate classifiers are presented in Figure 2.

Aggregation of Classifiers
The experimental results of the aggregation scheme used for the classifiers discussed in the previous section, namely decision tree, k-nearest neighbor, quadratic SVM, cubic SVM, linear SVM, logistic regression, random forest, and MLP, are discussed. The accuracies of the individual classifiers obtained in the initial series of experiments are the input to establish the fuzzy measure densities g i . The values of the function h are the results of the classification of testing elements being the probabilities of belonging to the two classes, namely, healthy and SZ patients. The validation procedure described in the previous section was repeated 200 times. After each run, a value of fuzzy density and 8 probability vectors (20% out of 40 observations) were obtained per classification model. The details of aggregation algorithm implementation required a single estimation of classification accuracy per model and aggregation method. Hence, all 1,600 (200 × 8) classification results were analyzed. It is worth noting that the models were fitted and the fuzzy densities were obtained independently in every separate run of the experiment. In the series of experiments, we have evaluated 25 classes (families) of popular and commonly considered in the literature triangular norms (Alsina et al., 2006, page 72). The monograph can be treated as a compendium of the t-norms to be applied in more advanced aggregation operators. In this particular approach, the t-norms serve as the integer functions M with parameters −10, −9.9, . . ., 0, . . ., 9.9, 10, but only if the parameter is in the range allowed for the t-norm. Such a choice of the parameter range seems to be optimal and emphasizes the most important properties of each of the t-norm classes. The maximal accuracy was obtained for the classifier C D2 for triangular norm from the family no. 8, namely   In this case, the best option to choose is α = 0.1. The plot illustrating the recognition rate for the combination of C D2 and the function (18) is given in Figure 3. It is obvious that satisfying accuracy can be obtained only for relatively small values of the parameter α.
It is important to understand the process of calculation of the value of the aggregation function. Let us consider, for example, the process of C D2 finding. Here, the individual classifiers x i , i = 1, ..., n = 5 (five classifiers are discussed in this experiment) should be analyzed, and for their significance measures (simply, weights) g i = g ({x i }), see Figure 2. Next, the parameter λ appearing in (6) is calculated. Based on the value of λ, the values of g (U i ) appearing in Eq. (5) are found recursively. Next, using the values of h (x i−1 ), which are the likelihoods of belongings of a given probe to a specific class, the final sum (16) can be obtained, taking into account that M (·, ·) is any t-norm, in particular a function given by (19).
Very good results were obtained also for the aggregation functions C M and very similar C FM . Maximal yielded values were 98.81% for the function number 12 serving as integer function and α = 0.2. The formula of the t-norm is as follows: It is worth stressing that the best average result among the operators C M , C FM , C CM , C MC , C MMin , C MMin2 , C D1 , C D2 , and C D3 was also obtained for the triangular norm no. 8 given by the formula (18) and its parameter α = 0.1. The plot presenting the values of the combination of all the aggregation functions with this t-norm and parameter is presented in Figure 4. It is obvious that the function (18) works well with almost all of the aggregation operators except C CM and C D3 .
As a supplement to the results, it is worth noting that Table 3 presents the information for which the best results of the t-norms Frontiers in Neuroinformatics | www.frontiersin.org were obtained by aggregation operators. It can help match t-norms and generalizations of the Choquet integral by experts conducting similar research. For instance, function no. 4, see Alsina et al. (2006), works well when it is combined with a few Choquet integral-based operators.
Finally, what should be emphasized, we have conducted the same tests for over 1,000 functions that are not Choquet-like integrals. The functions that were used in this competition were selected based on the studies by Alsina et al. (2006), Beliakov et al. (2007); Grabisch et al. (2009), andCalvo et al. (2012). The best results were obtained for the ordinary weighted averaging operator at the level of 98.81%.
where y j is the j-th largest of the x i , and the weights are ω 1 = 1, ω 2 = 1 -1 n , . . ., ω n = 1 n with or without normalization to their sum. Despite it being a very good result, it is obvious that it is hard to find the function giving the results more satisfying than Choquet integral-based operators. Moreover, to find the proper form of OWA, similar to the Choquet integral case, a proper heuristic can be used.

CONCLUSION AND FUTURE WORK
In this study, we have indicated the most appropriate operator aggregating the results of binary classification of patients to efficiently distinguish individuals with SZ and healthy subjects using a set of neural network organization features extracted from EEG-based functional connectivity measures. A series of both types of functions, generalizations of the Choquet integral and other aggregating functions, have been verified to determine the classes of functions and their parameters, which are the most effective in the classification of SZ. As an input to the main analysis, the results of classification were performed with classical methods such as decision tree, k-nearest neighbor, quadratic SVM, cubic SVM, linear SVM, logistic regression, random forest, and MLP were applied. The original results obtained in the study of classical methods classification reached 97% for logistic regression. Although the initially obtained results were high, we decided to verify if there is the possibility to reach even higher results using the fuzzy-based classifier.
The results prove that applying various classification models in combination with aggregation functions enable further improvement of classification results. This approach allows us to take advantage of the additional knowledge cumulated in the parameters of the trained models.
Detailed results show that several aggregation functions enabled to give promising results [presented in the study as Eqs (9-13) and (15, 16)], which increase the classification result by more than 1%. Among numerous functions evaluated and implemented in the thorough comparison, the best accuracy was reached for the aggregating integral C D2 with triangular norm appearing under the integral sign given by the formula (19).
Very good results (classification accuracy higher than 98.8%) were reached also for aggregation functions CM and CFM. It is worth noting that the obtained results occurred to be better than the original accuracy reached with classical methods by 1.81%. Although the obtained improvement is not very high (less than 2%), the overall increase in classification accuracy from 97% (for the best classical classifier) to as high as almost 99% (for the properly selected pre-aggregation operators) is relevant. Nevertheless, we have not done double cross-validation analysis, so this limitation can influence classification accuracy, hence the classification rate could be slightly overestimated. Results show the usefulness of this method especially if the role of aggregation function is an extended version of the Choquet integral. In contrast, the application of aggregation functions could give a relatively better improvement in case of weaker initial individual classification results. In future, we have planned to extend the analysis to consider more phases and stages of SZ. Moreover, we are interested in the application of other classes of aggregation operators and the determination of their weights (significance in the process of aggregation) based on the opinions of medical experts. More theoretically, it is still an interesting as well as difficult task to find the optimal parameters of the integral operators only based on their results according to various classification tasks with no relation to the accuracies or expert opinions. Finally, an application of aggregation techniques in other medical pattern recognition or classification problems will be worth analyzing.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Bioethics Committee at the Medical University of Lublin. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
MP-W and MK: contributed in methodology, software, validation, writing-original draft preparation, and visualization. PKa: contributed in conceptualization methodology, software, writing-original draft preparation, and visualization. MT: contributed in methodology, software, writing-original draft preparation, and visualization. PKr: contributed in methodology, validation, and writing-original draft preparation. KJ: contributed in methodology and writing-original draft preparation. All authors contributed to the article and approved the submitted version.