High-dimensionality time-structured data are extremely common in fields such as finance, biophysics and biomedical data. Very often, experimental limitations lead to uneven sampling (i.e. a different number of time points in terms of frequency or duration) [1] and this poses problems for many types of analysis (e.g. sample classification). As a consequence, clipping or padding techniques are applied, altering the underlying temporal structure. In recent years, studies on such data have seen an increasing popularity in a wide range of fields, from functional magnetic resonance imaging (fMRI) [2, 3, 4] to time series exploration for critical transition prediction in clinical scenarios [5, 6]. The common goal of this type of research is to develop models and algorithms capable of reaching the highest possible classification and prediction performances, for diagnostic and real time applications, while unveiling underlying information about a system. Reproducibility and generalization issues of commonly applied methods are in part caused by ad-hoc preprocessing of data, due to the lack of simple null models, often substituted by reshuffling-based null models. We introduce a method based on the statistical distribution of symmetric positive-definite matrices (i.e covariance and correlation matrices) extracted from data, using the Wishart distributon as a null model, as a possible way to overcome some of the aforementioned issues. Properties of distribution of random symmetric positive-definite matrices have proven to be useful in fields such as condensed matter, especially in the study of disordered systems [7, 8]. The WISDoM method exploits the properties of the Wishart distribution in order to compute limit distributions for the classes of samples in a classification problem, and a log-likelihood based score is defined for the single variables to quantify their relevance in the classification task.
2.1 The Wishart Distribution
The is a probability distribution of random nonnegative-definite
matrices that is used to model random covariance matrices.
The parameter n is the number of degrees of freedom (e.g. the number of points in the time series), and is a nonnegative-definite symmetric
matrix (with p the number of variables, or features, of the time series) called the scale matrix.
distribuited vectors, forming a data matrix
. The distribution of a
random matrix is a Wishart distribution. [9]
We have then by definition:
so that is the distribution of a sum of n rank-one matrices defined by independent normal
In particular, it holds for the present case:
2.2 PDF Computation for Invertible Σ
In general, any can be represented as
so that
where Assuming that
is invertible, the density of the random
matrix M in 1 can be written
so that symmetric and positive-definite. [10] Note that in 6 we define
generalized gamma function:
2.3 Estimation of the Wishart Parameters from Empirical Covariance
We justify the use of the Wishart distribution under the assumption of Multivariate Gaussian distributed data scenarios. This kind of assumption is indeed generally good for a wide range of problems. Furthermore, the use of the average covariance matrix (obtained from all the elements of one class) to compute the scale matrix for the class estimated distribution will be proven to be a good approximation of a complete Bayesian model.
This is done by considering that the Wishart Distribution is the conjugate prior of a multivariate Gaussian distribution, such as the Gamma distribution for the univariate Gaussian case. By considering a Gaussian model with known mean , so that the free parameter is the variance
, as in [11], the likelihood function is defined as follows:
The conjugate prior is an inverse Gamma distribution. Recall that inverse Gamma distribution with parameters
The density then takes the form
Using this prior, the posterior distribution of is given by
In the multidimensional setting, the inverse Wishart takes the place of the inverse Gamma. It has already been stated that the Wishart distribution is a distribution over more compact form of the density is given by
where the parameters are the degrees of freedom and the positive-definite
If we can then state that W has an Inverse Wishart Distribution, whose density has the form
Let distributed observed data. Then an inverse Wishart prior multiplying the likelihood
where S is the empirical covariance :
Thus, an a posteriori distribution with the form
is obtained. Similarly, it can be stated that for the inverse covariance (precision) matrix the conjugate prior is a Wishart distribution.
2.4 Class-Wise Estimated Distribution
The core idea of the WISDoM method is to represent each element undergoing classification as a covariance matrix of its features. Nominally, each element can be characterized by the covariance matrix extracted by the repeated observations of the vector of its features, for example derived by a time series. The aim is to use the free parameters of the Wishart distribution (the scale matrix and the number n of the degrees of freedom, as shown in 6) to compute an estimation of the distribution for a certain class of elements, and then assign a single element to a given class by computing a log-likelihood between the element being analyzed and each class. Furhermore, a score can be assigned to each feature by estimating the variation in terms of log-likelihood, due to its removal from the feature set. If the removal of a feature causes significant increase (or decrease) in the log-likelihood, it can be stated that such feature is highly representative of the system analyzed. Thus, the WISDoM approach allows not only to assign a given element to a class, but also to identify the features with the highest relevance in the classification process.
Covariance matrices are a good choice for a distance metrics in a classi-fication task, both for the way they represent a system and for the property that the average of a set of covariance matrices is a covariance matrix itself. If each element of a given class C is represented by a covariance matrix its features, this property allows us to estimate a distribution for the class by choosing
The other necessary parameter for the estimation is the number of degrees of freedom n. Assume that an features is associated to each element i of a given class, while having n observations for this vector. The covariance matrix
computed over the n observations will represent the interactions between the features of element i. The number of degrees of freedom n of the Wishart distribution is then given by the number of times
is observed.
Let us give an example tied to functional MR brain imaging. An image of patient i’s brain is acquired; as usual these images are divided in a certain number p of zones (voxel, pixel etc.), each zone being sampled n times over a given time interval in order to observe a certain type of brain activity and functionality. In this example, the features contained in vector associated to patient i are indeed the zones chosen to divide the brain image, each zone having been sampled n times during an acquisition interval. The
correlation matrix
is then representative of the functional correlation between the p brain areas. Repeating this procedure for the N patients of a known class C (i.e. a diagnostic group) and computing the
scale matrix for the class, will allow us to estimate a Wishart distribution for that class and draw samples from it.
2.5 Log-Likelihood Ratio Score
After defining how to represent classes distribution, WISDoM allows to compute the log-likelihood of each element to belong to one of the classes. Moreover, WISDoM allows to compute the variation of log-likelihood ratio scores due to the removal of features, singularly or in groups, thus estimating how much the classification performance changes. Uninformative (or less informative) features can thus be pruned, allowing for a dimensionality reduction of the initial feature set. The whole process can be seen as a feature transformation, mapping the covariance matrix of subject i to a score vector formed by the change in log-likelihood for each feature.
Complete Matrix Score
The WISDoM Classifier relies upon computing the log-likelihood of a matrix with respect to the Wishart distribution estimated for a class
as the scale marix. If a problem concerning two given classes
taken into account, the score assigned to each
can be defined as follows:
tively, and is the logarithm of the probability of
belonging to the Wishart distribution estimated for one of the two classes A, B.
Single Feature Score
WISDoM allows to obtain information about the features used for classifica-tion by reducing the matrix A to its principal submatrices (see Appendix). An important property for the principal submatrices of a symmetric positive definite matrix is that positive definite.
By removing one feature from the dataset, calculating the WISDoM scores, and iterating this process over all the features (i.e. analyzing all the principal submatrices of
) the method can assign a score to each feature, representing its relevance in the decision for
assigned to one class or another. Let
be a principal submatrix of order
, of the matrix
computed on the observation of
subject
. Similarly,
let be a principal submatrix of order
, of the matrix
for the class C. The score assigned to each feature of
given by eq.(23).
A generalization to dimensionality reduction can be found in the supplemental materials.
3.1 Eye state detection via EEG
The dataset used was downloaded from the UCI Machine Learning Repository (http://archive.ics.uci.edu/ml/). This dataset has been chosen for many reasons: it’s openly accessible, contains records from 14 electrodes with standard headset placement (fig.1), thus making the features of our problem directly linked to brain topology and a published classification performance benchmark on the dataset exists [12]. The data consisted in a series of 14980 time points, sampled for each one of the 14 electrodes and labelled with a 1 or a 0 to mark wether the eyes of the subject are open or closed at that time point. The time series has been split into batches of different length according to eye state changes. In this way, a correlation matrix can be extracted for each batch (the "elements" for this classification problem), while the length of each batch is used for computing the degrees of freedom of each class Wishart distribution during training. A total of 140 batches with various lenghts, 70 with eye state 1 and 70 with eye state 0, were obtained.
The representative matrix for each class is computed as the average (weighted on the length of each batch) of matrices of the elements belonging to eye state 0 or eye state 1, excluding the element to be predicted in a Leave One Out fashion in order to avoid overfitting. By doing this, we verify that the method is independent from the sampling window chosen when applied to time series data, with the only constraint that the length of such window cannot be less than the number of the features of the system.
After undergoing feature score computation, a stochastic grid search on a set of classifiers has been performed in order to obtain the best prediction performance with the transformed features. All the classification tasks are validated through a 10-fold crossvalidation. Versions and references for all Python packages used can be found in supplemental and at [13, 14, 15, 16, 17, 18, 19, 20].
We first tried to assess eye state using complete matrix score, as in eq. (21). Classifiers reported in fig.(2) were trained and tuned, with the aim of obtaining the best performance possible. However, in this scenario the resulting classification performances were poor, reaching an accuracy of in the best cases. We then proceeded to compute single feature scoring, as in eq.(23), obtaining a feature transformation. As in fig.(2), different classi-fiers belonging to two main categories (decision trees and linear classifiers) have been trained on the transformed features. The best performance has been achieved with a C-support Vector Machine (Python 3.6 SciKitLearn implementation) resulting in a 0.85% ROC AUC score and an accuracy score of 84.3%, comparable with the benchmark of 83.5% accuracy set by Rajesh Kanna et al.[12].
To assess which features contain the largest amount of useful information for prediction, a set of single feature C-SVM classifications has been performed (fig.3): a performance of 75% accuracy is obtained by using only the top three ranking electrodes (fig. 4). Training the classifier with the top three electrodes yields a local maxima in the landascape performance, highlighting the importance of the information recorded by these three electrodes about the state of the whole system.
3.2 Autism classification via fMRI
The Autism Brain Imaging Data Exchange [21] is a consortium effort, aggregating fMRI datasets from individuals with autism spectrum disorders (ASD) and age-matched typical controls (TC). Data from 17 acquisition sites were merged, using different preprocessing tools and pipelines [22]. Complete information about the dataset is found at http://preprocessed-connectomes-project.org/abide/. For our classification task , we focused on male subjects of the "Autism" diagnostic group (AUT): we analyzed a total of 369 TC and 220 AUT subjects, with 200 time points each (number of degrees of freedom of the Wishart distribution). The chosen preprocessing pipeline for the extraction of the average time series of the ROIs is the Connectome Computation System (CCS), with a global signal correction and the application of a bandpass filter (0,01-0,1 Hz). The 116 ROIs (features), of which covariance matrices are extracted, are labelled according to the Automated Anatomical Labeling of the Brain (AAL2) [23].
The representative matrix for each class is computed as the average of matrices of the elements belonging to class AUT or TC, excluding the sets of element to be predicted. For this task, a shuffled 10-fold splitting of the dataset for feature score computation has been used to avoid overfitting. After undergoing feature score computation, the parameters of a C-Support Vector Machine Classifier have been fine tuned in order to obtain the best prediction performance with the transformed features. All the classification tasks are validated through a stratified 10-fold crossvalidation, in order to minimize the effects of class imbalance in train and test sets. Versions and references for all Python packages used can be found in supplemental and at [13, 14, 15, 16, 17, 18, 19, 20].
The C-SVM classifier trained on the transformed features resulted in an accuracy score of 72.1% and a ROC AUC score of 0.76. Furhtermore, we obtained a ROC AUC score of 0.79 and an accuracy of 73.5% with a fine tuned Random Forest Classifier. We also compared the classifiers in fig.2 training them with WISDoM transformed features and non-transformed features (nominally, the elements of the lower triangle of each covariance matrix). Results in fig.5 show an overall improvement of classification performances when using transformed features. As a comparison, the state of art of clas-sification on the entire ABIDE dataset is set at 70% accuracy obtained with a deep learning architecture built by Heinsfeld et al. [24]. This result on the whole spectrum of autism required the use of various stacked denoising autoencoders and hidden layers, resulting in a large time-consuming training routine (hours), while WISDoM obtained satisfying classification performance in much smaller time (
minutes, including feature transformation which is the most time-consuming step of the pipeline).
The WISDoM framework is introduced: a method for modelling symmetric positive definite matrices, such as covariance and correlation matrices, used in a wide array of problems. It can provide a null model for classification purposes in which each sample is represented as a covariance/correlation matrix, even if the number of observations (e.g. the length of the time series) is different from sample to sample. This property makes the WISDoM method suitable for problems with non-homogeneous data size, for example time series with uneven lengths, missing points or irregularly sampled data. Moreover, we show that a feature transformation based on WISDoM scores can be used for dimensionality reduction, providing a ranking for the most important variables in the dataset. While showing good generalization capabilities with time-series data and non-homogeneous sampling related issues, the method is not suitable when the number of features exceeds the sampling (p > n). This is a theoretical limit tied to the invertibility of the scale matrix required to compute the Wishart probability density function. At present, WISDoM cannot thus be applied to problem involving the so called "long data", such as gene expression tables, unless considering corrections such as matrix regularization methods and hierchical methods such as power priors [25]
The method has been tested on the EEG eye state prediction dataset of the open UCI Machine Learning Repository, slightly improving the previous classification benchmark with little to no preprocessing, and giving useful insights on the minimum number and location of electrodes needed to record sufficient information for the task. Moreover, the method has been applied to the classification of a subset of the ABIDE dataset, using brain functional connectivity data. We obtained satisfying classification scores, comparable with the state of art classification results on the dataset, with very simple classifiers and without the use of additional time-consuming processing routines. Furthermore, the Bayesian-like framework of scores-computation through log-likelihood, could allow for a sort of inline learning by continuously updating the estimation of each class Wishart distribution. This property makes the WISDoM method also suitable for real-time learning during
data acquisition.
EG, CM, GC and DR designed the research. CM and EG analyzed the data and implemented the method. CM, DR, GC and EG wrote the paper.
[1] A. Bazzani, G. Castellani, and L. Cooper, “Eigenvalue distributions for a class of covariance matrices with application to bienenstock-cooper-munro neurons under noisy conditions,” Physical Review E, vol. 81, 2010.
[2] M. Van Den Heuvel and H. E. H. Pol, “Exploring the brain network: a review on resting-state fmri functional connectivity,” European Neuropsychopharmacology, vol. 20, pp. 519–534, 2010.
[3] W. Yan, V. Calhoun, M. Song, Y. Cui, H. Yan, S. Liu, L. Fan, N. Zuo, Z. Yang, K. Xu, J. Yan, L. Lv, J. Chen, Y. Chen, H. Guo, L. Li, Pand Lu, P. Wan, H. Wang, H. Wang, Y. Yang, H. Zhang, D. Zhang, T. Jiang, and J. Sui, “Discriminating schizophrenia using recurrent neural network applied on time courses of multi-site fmri data,” EBioMedicine, vol. 47, pp. 543–552, 2019.
[4] F. Azarmi, S. Miri Ashtiani, A. Shalbaf, H. Behnam, and M. Daliri, “Granger causality analysis in combination with directed network measures for classification of ms patients and healthy controls using taskrelated fmri,” Comput Biol Med, vol. 115, 2019.
[5] D. Cuesta-Frau, P. Miró-Martínez, S. Oltra-Crespo, A. Molina-Picó, P. Dakappa, C. Mahabala, B. Vargas, and P. González, “Classification of fever patterns using a single extracted entropy feature: A feasibility study based on sample entropy,” Math Biosci Eng., vol. 17, pp. 235–249, 2019.
[6] P. Ghalati, S. Samal, J. Bhat, R. Deisz, G. Marx, and A. Schuppert, “Critical transitions in intensive care units: A sepsis case study,” Sci Rep., vol. 9, 2019.
[7] H. Zhu, S. Zou, and T. Hastie, “Multi-class adaboost,” Statistics and Its Interface, vol. 2, p. 349–360, 2009.
[8] A. Crisanti, G. Paladin, M. Serva, and A. Vulpiani, “Products of random matrices for disordered systems,” Physical Review E, vol. 49, p. R953–R955, 1994.
[9] Hardle, Wolfgang, and L. Simar, Applied Multivariate Statistical Analysis. Springer Berlin, 2003.
[10] Anderson, T. W.Hardle, Wolfgang, and L. Simar, An Introduction to Multivariate Statistical Analysis. John Wiley and Sons, 2012.
[11] H. Liu and L. Wasserman, Statistical Machine Learning. CMU University, 2014.
[12] K. Rajesh, V. Sabarinathan, V. Sarath Kumar, and V. Sugumaran, “Eye state prediction using eeg signal and c4.5 decision tree algorithm,” International Journal of Applied Engineering Research, vol. 10, 2015.
[13] E. Jones, T. Oliphant, P. Peterson, et al., “SciPy: Open source scientific tools for Python,” 2001–.
[14] W. McKinney, “Data structures for statistical computing in python,” in Proceedings of the 9th Python in Science Conference (S. van der Walt and J. Millman, eds.), pp. 51 – 56, 2010.
[15] J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing in Science & Engineering, vol. 9, no. 3, pp. 90–95, 2007.
[16] O. Travis E, A guide to NumPy. Trelgol Publishing, 2006.
[17] A. Meurer, C. P. Smith, M. Paprocki, O. Certík, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh, T. Rathnayake, S. Vig, B. E. Granger, R. P. Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson, F. Pedregosa, M. J. Curry, A. R. Terrel, S. Roucka, A. Saboo, I. Fernando, S. Kulal, R. Cimrman, and A. Scopatz, “Sympy: symbolic computing in python,” PeerJ Computer Science, vol. 3, p. e103, Jan. 2017.
[18] T. Kluyver, B. Ragan-Kelley, F. Pérez, B. Granger, M. Bussonnier, J. Frederic, K. Kelley, J. Hamrick, J. Grout, S. Corlay, P. Ivanov, D. Avila, S. Abdalla, and C. Willing, “Jupyter notebooks – a publishing format for reproducible computational workflows,” in Positioning and Power in Academic Publishing: Players, Agents and Agendas (F. Loizides and B. Schmidt, eds.), pp. 87 – 90, IOS Press, 2016.
[19] F. Pérez and B. E. Granger, “Ipython: A system for interactive scientific computing,” Computing in Science & Engineering, vol. 9, no. 3, pp. 21– 29, 2007.
[20] M. Waskom, O. Botvinnik, D. O’Kane, P. Hobson, S. Lukauskas, D. C. Gemperline, T. Augspurger, Y. Halchenko, J. B. Cole, J. Warmenhoven, J. de Ruiter, C. Pye, S. Hoyer, J. Vanderplas, S. Villalba, G. Kunter, E. Quintero, P. Bachant, M. Martin, K. Meyer, A. Miles, Y. Ram, T. Yarkoni, M. L. Williams, C. Evans, C. Fitzgerald, Brian, C. Fonnesbeck, A. Lee, and A. Qalieh, “mwaskom/seaborn: v0.8.1 (september 2017),” Sept. 2017.
[21] A. Di Martino, C. Yan, Q. Li, E. Denio, F. Castellanos, K. Alaerts, J. Anderson, M. Assaf, S. Bookheimer, M. Dapretto, B. Deen, S. Delmonte, I. Dinstein, B. Ertl-Wagner, D. Fair, L. Gallagher, D. Kennedy, C. Keown, C. Keysers, J. Lainhart, C. Lord, B. Luna, V. Menon, N. Minshew, C. Monk, S. Mueller, R. Muller, M. Nebel, J. Nigg, K. O’Hearn, K. Pelphrey, S. Peltier, J. Rudie, S. Sunaert, M. Thioux, J. Tyszka, L. Uddin, J. Verhoeven, N. Wenderoth, J. Wiggins, S. Mostofsky, and M. Milham, “The autism brain imaging data exchange: towards a largescale evaluation of the intrinsic brain architecture in autism,” Molecular Psychiatry, vol. 19, pp. 659–667, 2014.
[22] C. Craddock, Y. Benhajali, C. Chu, F. Chouinard, A. Evans, A. Jakab, B. S. Khundrakpam, J. D. Lewis, Q. Li, M. Milham, C. Yan, and P. Bellec, “The neuro bureau preprocessing initiative: open sharing of preprocessed neuroimaging data and derivatives,” Frontiers in Neuroinformatics, no. 41.
[23] N. Tzourio-Mazoyer, B. Landeau, D. Papathanassiou, F. Crivello, O. Etard, N. Delcroix, B. Mazoyer, and M. Joliot, “Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni mri single-subject brain,” NeuroImage, vol. 15, no. 1, pp. 273 – 289, 2002.
[24] A. Heinsfeld, A. Franco, R. Craddock, A. Buchweitz, and F. Meneguzzi, “Identification of autism spectrum disorder using deep learning and the abide dataset,” Neuroimage Clin., vol. 17, pp. 16–23, 2017.
[25] J. Ibrahim, M. Chen, Y. Gwon, and F. Chen, The power prior: theory and applications. Wiley, 2015.
[26] L. Breiman, “Random forests,” Machine Learning, vol. 45, pp. 5–32, 2001.
[27] L. Breiman, J. Friedman, R. Olshen, and C. Stone, Classification and Regression Trees. Wadsworth, 1984.
[28] C. Barros Rodrigo, M. P. Basgalupp, A. C. P. L. F. Carvalho, and A. Freitas Alex, “A survey of evolutionary algorithms for decision-tree induction,” IEEE Transactions on Systems, Man and Cybernetics. Part C: Applications and Reviews, vol. 42, p. 291–312, 2012.
[29] Y. Freund and R. Schapire, “A decision-theoretic generalization of on- line learning and an application to boosting,” Journal of Computer and System Sciences, vol. 55, pp. 119–139, 1997.
[30] M. G. J., Discriminant Analysis and Statistical Pattern Recognition. Wiley Interscience, 2004.
[31] Y. Hsiang-Fu, H. Fang-Lan, and L. Chih-Jen, “Dual coordinate descent methods for logistic regression and maximum entropy models,” Machine Learning, vol. 85, pp. 41–75, 2011.
[32] F. Y. and S. R. E., “Large margin classification using the perceptron algorithm,” Machine Learning, vol. 37, p. 277–296, 1999.
[33] C.-C. Chang and C.-J. Lin, “Libsvm: A library for support vector ma- chines,” 2013–. Department of Computer Science National Taiwan University.
[34] J. Platt, “Probabilistic outputs for support vector machines and com- parison to regularized likelihood methods,” in Advances in Large Margin Classifiers, 1999.
Figure captions
Figure 1: Electrodes position in the headset used for EEG dataset acquisition.
Figure 2: Performance comparison of different classifiers on the WISDoM EEG transformed features. The classifiers are reported as follows: RFC: Random Forest Classifier [26], DTC: Decision Tree Classifier[27, 28], ADA: ADA Boosting Tree Classifier [29, 7], LDA: Linear Discriminant Analysis Classification [30], LogReg: Logistic Regression Classifier [31], Perc: Perceptron Classifier [32], SVM: C-Support Vector Machine[33, 34]. All classifiers are SciKitLearn implementations.
Figure 3: EEG feature ranking based classification performance (ROC AUC score). Temporal (T8) and outer frontal (F8, F7) electrodes seems to convey the most important signals for eye state prediction.
Figure 4: EEG features performance landscape at an increasing number of ranked features used for classification. Labels on the X-axis point out which feature is being added to the previous ones.
Figure 5: Performance comparison of different classifiers on the WISDoM transformed features and non-transformed feature in the ABIDE dataset. The classifiers are reported as follows: RFC: Random Forest Classifier [26], DTC: Decision Tree Classifier[27, 28], ADA: ADA Boosting Tree Classifier [29, 7], LDA: Linear Discriminant Analysis Classification [30], LogReg: Logistic Regression Classifier [31], Perc: Perceptron Classifier [32], SVM: C-Support Vector Machine[33, 34]. All classifiers are SciKitLearn implementations.