dbcsp : User-friendly R package for Distance-Based Common Spacial Patterns

Common Spacial Patterns (CSP) is a widely used method to analyse electroencephalography (EEG) data, concerning the supervised classiﬁcation of brain’s activity. More generally, it can be useful to distinguish between multivariate signals recorded during a time span for two diﬀerent classes. CSP is based on the simultaneous diago-nalization of the average covariance matrices of signals from both classes and it allows to project the data into a low-dimensional subspace. Once data are represented in a low-dimensional subspace, a classiﬁcation step must be carried out. The original CSP method is based on the Euclidean distance between signals and here, we extend it so that it can be applied on any appropriate distance for data at hand. Both, the classical CSP and the new Distance-Based CSP (DB-CSP) are implemented in an R package, called dbcsp .


Background
Eigenvalue and generalized eigenvalue problems are very relevant techniques in data analysis.The well-known Principal Component Analysis with the eigenvalue problem in its roots was already established by late seventies [Mardia et al., 1979].In mathematical terms, Common Spatial Patterns (CSP) is based on the generalized eigenvalue decomposition or the simultaneous diagonalization of two matrices to find projections in a low dimensional space.Although in algebraic terms PCA and CSP share several similarities, their main aims are different: PCA follows a non-supervised approach but CSP is a two-class supervised technique.Besides, PCA is suitable for standard quantitative data arranged in a 'individuals × variables' tables, while CSP is designed to handle multivariate signals time series.That means that while for PCA each individual or unit is represented by a classical numerical vector, for CSP each individual is represented by several signals recorded during a time span, i.e., by a 'number of signals × time span' matrix.CSP allows to represent the individuals in a dimension reduced space, a crucial step given the high dimensional nature of the original data.CSP computes the average covariance matrices of signals from the two classes to yield features whose variances are optimal to discriminate the classes of measurements.
Once data is projected into a low dimensional space, a classification step is carried out.The CSP technique was first proposed under the name Fukunaga-Koontz Transform in Fukunaga and Koontz [1970] as an extension of PCA and Müller-Gerking et al. [1999] used it to discriminate electroencephalography data (EEG) in a movement task.Since then, it is a widely used technique to analyze EEG data and develop Brain Computer Interfaces (BCI), with different variations and extensions [Blankertz et al., 2007a,b, Grosse-Wentrup and Buss, 2008, Lotte and Guan, 2011, Wang et al., 2012, Astigarraga et al., 2016, Darvish Ghanbar et al., 2021].Samek et al. [2014] offer a divergence-based framework including several extentions of CSP.As a general term, CSP filter maximizes the variance of the filtered or projected EEG signals of one class of movements while minimizing it for the signals of the other class.Similarly, it can be used to detect epileptic activities Khalid et al. [2016] or other brain activities.Together with, BCI systems can be of great help to people who suffer from some disorders of cerebral palsy, or who suffer from other diseases or disabilities that prevent the normal use of their motor skills.These systems can considerably improve the quality of life of these people, for which small advances and changes imply big improvements.BCI systems can also contribute to the human vigilance detection, connected with occupations involving sustained attention tasks.Among others, CSP and variations of it have been applied to the vigilance estimation task [Yu et al., 2019].
The original CSP method is based on the Euclidean distance between signals.However, as far as we know, it was not introduced a generalization allowing the use of any appropriate distance.The aim of the present work is to introduce a novel Distance-Based generalization of it (DB-CSP).This generalization is of great interest, since these techniques can also offer good solutions in other fields where multivariate time series data arise beyond pure electroencephalography data [Poppe, 2010, Rodríguez-Moreno et al., 2020].
Although CSP in its classical version is a very well-known technique in the field of BCI, it is not implemented in R. In addition, being DB-CSP a new extension of it, it is worth building an R package that includes both, CSP and DB-CSP techniques.The package offers functions in a user friendly way for the less familiar users of R but it also offers complete information in its objects so that reproducible analysis can be done, as well as more advanced and customised analysis can be performed taking advantage of already well-known packages of R.
The paper is organized as follows.First, we review the mathematical formulation of the Common Spatial Patterns method.Below, the core of our contribution, we describe both the novel CSP' extension based on distances and the package dbcsp.Then, the main functions in dbcsp are introduced along with reproducible examples of their use.Finally, some conclusions are drawn.

CSP and DB-CSP
Let us consider we have n statistical individuals or units classified in two classes C 1 and C 2 , with #C 1 = n 1 and #C 2 = n 2 .For each unit i in class C k data from c sources or signals are collected during T time units and therefore unit i is represented in matrix X ik (i = 1, . . ., n k ; k = 1, 2).For instance, for electroencephalograms, data are recorded by a c-sensor cap each t time units (t = 1, . . ., T ).As usual, we consider that each X ik is already scaled or with the appropriate pre-processing in the context of application; for instance, if working with EGG data, each signal should be band-pass filtered before its use.
The goal is to classify a new unit X in C 1 or C 2 .To this end, first a projection into a low-dimensional subspace is carried out.Then, classically the Linear Discriminant classifier (LDA) is applied taking as input data for the classifier the log-variance of the projections obtained in the first step.It is obvious that the importance of the technique lies mainly in the first step, and once it is done, LDA or any other classifiers could be applied.Based on that, we focus on how this projection into a low-dimensional space is done, from the classical CSP point of view as well as its novel extension DB-CSP (see Figure 1).

Classical CSP
The main idea is to use a linear transform to project or filter data into low-dimensional subspace with a projection matrix, in such a way that each row consists of weights for signals.This transformation maximizes the variance of two-class signal matrices.The method performs a simultaneous diagonalization of the covariance matrices of both classes.Given data X 11 , . . ., X n 1 1 (matrices c × T ) from class C 1 and X 12 , . . ., X n 2 2 (also matrices c × T ) from class C 2 , the following steps are needed: • All matrices are standardized so that traces of X ik X ′ ik are the same.• Compute average covariance matrices: • Look for directions W = (w 1 , . . ., w c ) ∈ R c×c according to the criterion: The solution is given by the generalized spectral decomposition B 1 w = λB 2 w choosing the first and the last q eigenvectors: W CSP = (w 1 , . . ., w q , w c−q+1 , . . ., w c ).
Vectors w j offer weights so that new signals X ′ i1 w j and X ′ i2 w j have big and low variability for the first q vectors (j = 1, . . ., q) respectively, and vice-versa for the last q vectors (j = c − q + 1, . . ., c).To clarify the notation and interpretation, let us denote a j = w j the first q vectors and b j = w c+1−j the last q.That way, and broadly speaking, variability of elements in C 1 is big when projecting on vectors a j and low on vectors b j , and vice-versa, for elements in class C 2 .
Finally, the log-variability of these new and few 2q signals are considered as input for the classification, which classically is the Linear Discriminant Analysis.Obviously, any other classification technique can be used.

Distance-based CSP
Following the commented ideas, the Distance-Based CSP (DB-CSP) is an extension of the classical CSP method.In the same way as the classical CSP, DB-CSP gives some weights to the original sources or signals and obtains new and few 2q signals useful for the discrimination between the two classes.Nevertheless, the considered distance between the signals sources can be any other than the Euclidean.The steps are the following: • Compute an appropriate distance measure between sources and the double-centered inner product: where H stands for the centering matrix and the superindex in brackets (2) for squared elements in the matrix.Again, all matrices are standardized so that all traces of X ik X ′ ik are the same.• Compute average distance-based covariance matrices: Once we have the covariance matrices related to the chosen distance matrix, the directions are found as in classical CSP and new signals X ′ ik a j , X ′ ik b j are built (j = 1, . . ., q).Again, for individuals in class C 1 the projections on vectors a and b are big and low respectively; for individuals in class C 2 the other way round.
It is important to note that if the chosen distance does not produce a positive definite covariance matrix, it must be replaced by a similar one that is positive definite.
When the selected distance is the Euclidean, then, DB-CSP reduces to classical CSP.
Once the q directions a j and b j are calculated, new 2q signals are built.Many interesting characteristics of the new signals could be extracted, although the most important in the procedure is the variance.Those characteristics of the new signals are the input data for the classification step.

Implementation
In this section, the structure of the package and the functions implemented are explained.The dbcsp package was developed for the free statistical R environment (http://www.r-project.org) and it is available from the Comprehensive R Archive Network (CRAN) at https://cran.r-project.org/web/packages/dbcsp/index.html.

Input
The input data are the corresponding n 1 and n 2 matrices X ik of the n units classified in classes C 1 and C 2 , respectively (i = 1, . . ., n k ; k = 1, 2).Let x1 and x2 be two lists of length n 1 and n 2 , respectively, with X ik matrices (c × T ) as elements of the lists.We also consider that new items to be classified are in list xt.The aforementioned first step of the method is carried out by building the object called "dbcsp".

dbcsp object
The dbcsp object is an S4 class created to compute the projection vectors W .The object has the following slots: • Slots X1 = "list", X2 = "list", the lists X1 and X2 (lengths n 1 and n 2 ) containing the matrices X ik for the two classes C 1 and C 2 , respectively (i = 1, . . ., n k ; k = 1, 2).q = "integer", to determine the number of pairs of eigenvectors a j and b j that are kept.By default q =15.labels = "character", vector of two strings indicating labels names, by default names of elements in X1 and X2.type = "character", to set the type of distance to be considered, by default type='EUCL'.The supported distances are these ones: - out= "list", list containing elements of the output.Mainly, matrix W with vectors a j and b j in element vectors, log-variances of filtered signals in proy and partitions considered in the k-fold approach with reproducibility purposes.
• Usage Following the standard procedure in R, an instance of a class dbcsp is created via the new() constructor function: Slots X1 and X2 are compulsory since they contain the original data.When a slot is not specified, the default value is considered.First, the S4 object of class dbcsp must be created.By default, the Euclidean distance is used, nevertheless it can be changed.For instance, "Dynamic Transform Distance" [Giorgino et al., 2009] can be set: or a mixture between this distance and the euclidean can be indicated by: mydbcsp.mix<-new('dbcsp', X1=x1, X2=x2, labels=c("C1", "C2"), mixture=TRUE, w=0.4,type="dtw") The object contains all the information to carry out classification task in a lower dimension space.

Functions plot and boxplot
For exploratory and descriptive purposes, the original signals X ik and the projected ones can be plotted for selected individual i in class k, and selected pair of dimensions a j and b j (i = 1, . . ., n k , k = 1, 2).
• Usage plot(mydbcsp) Besides, the log-variances of the projected signals of both classes can be shown in boxplots.This graphic can help to understand the discriminative power that is in the low-dimension space.It is worthy to take into account that in the aforementioned functions, values in argument vectors must lie between 1 and 2q, being q the number of dimensions used to perform the DB-CSP algorithm when creating the dbcsp object.Therefore, values 1 to q correspond to vectors a 1 to a q and values q + 1 to 2q correspond to vectors b 1 to b q .Then, if pairs=TRUE, it is recommended that values in argument vectors are in {1, . . ., q}, since their pairs are plotted as well.In case that values are above q, it should be noted that they correspond to vectors b 1 to b q .For instance, if q=15 and boxplot(object, vectors=16, pairs=FALSE), vector b 1 (16 − q = 1) is shown.

Function selectQ, Function train and Function predict
The functions in this section help the classification step in the procedure.Function selectQ helps to find an appropriate dimension needed for the classification.Given different values of dimensions, the accuracy related to each dimension is calculated so that the user can assess which dimension of the reduced space can be sufficient.A k-fold cross-validation approach or a holdout approach can be followed.Function train performs the Linear Discriminant classification based on the log-variances of the dimensions built in the dbcsp object.The accuracy of the classifier is computed based on k-fold validation procedure.Finally, function predict performs the classification of new individuals.
• Usage of selectQ selectQ(mydbcsp) • Arguments object, an object of class dbcsp Q, vector of integers which represents the dimensions to use, by default Q=c(1,2,3,5,10,15).train_size, float between 0.0 and 1.0 representing the proportion of the data set to include in the train split, by default train_size=0.75.
CV, logical indicating whether a k-fold cross validation must be performed or a hold-out approach (if TRUE, train_size is not used), by default CV=FALSE.
folds integer, number of folds to use if CV is performed.
seed numeric value, by default seed=NULL.Set a seed in case you want to be able to replicate the results.
This functions returns the accuracies related to each dimension set in Q.If CV=TRUE, the mean accuracy as well as the standard deviation among folds is also returned.
• Usage of train train(mydbcsp) or embedded as a parameter in: new('dbcsp', X1=x1, X2=x2, training=TRUE, type="dtw") • Arguments object, an object of class dbcsp selected_q, integer value indicating the number of vectors to use when training the model.By default all dimensions considered when creating the object dbcsp.
Besides, arguments seed and fold are available.
It is important to note that in this way a classical analysis can be carried out, in the sense of: • LDA is applied based on the log-variances of the dimensions indicated by the user in, select_q; • percentage of good classification is obtained via k-fold cross validation.
However, it is evident that it may be of interest to use other classifiers or other characteristics in addition to or different from log-variances.This more advanced procedure is explained below.See the basic analysis of the "User guide with real example" section in order to visualize and follow the process of a first basic/classic analysis.
• Usage of predict predict(mydbcsp, X_test=xt) • object, an object of class dbcsp X_test, list of matrices for be classified.
true_targets, optional, if available, vector of true labels of the instances.Note that they must match the name of the labels used when training the model.

User guide with a real example
To show an example beyond pure electroencephalography data, an Action Recognition data is considered.Besides, in order to have a reproducible example to show the use of the implemented functions and the results they offer, this Action Recognition data set is included in the package.The data set contains the skeleton data extracted from videos of people performing six different actions, recorded by a semi-humanoid robot.
It consists of a total of 272 videos with 6 action categories.There are around 45 clips in each category, performed by 46 different people.Each instance is composed of 50 signals (xy coordinates for 25 body key points extracted using OpenPose [Cao et al., 2019]), where each signal has 92 values, one per frame.These are the six categories included in the data set: 1. Come: gesture for telling the robot to come to you.There are 46 instances for this class.
2. Five: gesture of 'high five'.There are 45 instances for this class.
3. Handshake: gesture of handshaking with the robot.There are 45 instances for this class.
4. Hello: gesture for telling hello to the robot.There are 44 instances for this class.
5. Ignore: ignore the robot, pass by.There are 46 instances for this class.

Basic/classic analysis
Suppose an analysis using 15-dimensional projections and the Euclidean distance.At a first step the user can obtain vectors W by: x1 <-AR.data$comex2 <-AR.data$fivemydbcsp <-new('dbcsp', X1=x1, X2=x2, q=15, labels=c("C1", "C2")) summary(mydbcsp) Creating the object mydbcsp, the vectors W are calculated.As indicated in parameter q=15, the first and last 15 eigenvectors are retained.With summary, the obtained output is: There are 46 instances of class C1 with [50x92] dimension.There are 45 instances of class C2 with [50x92] dimension.The DB-CSP method has used 15 vectors for the projection.EUCL distance has been used.Training has not been performed yet.Now, if the user knows from the beginning that 3 is an appropriate dimension, the classification step could be done while creating the object.Using classical analysis, with for instance 10-fold, LDA as classifier and log-variances as characteristics, the corresponding input and summary output are: mydbcsp <-new('dbcsp', X1=x1, X2=x2, q=3, labels=c("C1", "C2"), training=TRUE, fold = 10, summary(mydbcsp) There are 46 instances of class C1 with [50x92] dimension.There are 45 instances of class C2 with [50x92] dimension.The DB-CSP method has used 3 vectors for the projection.EUCL distance has been used.An accuracy of 0.9130556 has been obtained with 10 fold cross validation and using 3 vector If a closer view of the accuracies among the folds are needed, the user can get them from the out slot of the object: # Accuracy in each fold mydbcsp@out$folds_acc # Intances belonging to each fold mydbcsp@out$used_folds 4.2 Basic/classic analysis selecting the value of q Furthermore, it is clear that the optimal value of q should be chosen based on the percentages of good classification.It is worth mentioning that the LDA is applied on the 2q projections, as set in the object building step.It is interesting to measure how many dimensions would be enough using selectQ function: mydbcsp <-new('dbcsp', X1=x1, X2=x2, labels=c("C1", "C2")) selectDim <-selectQ(mydbcsp, seed=30, CV=TRUE, fold = 10) selectDim Q acc sd 1 1 0.7663889 0.12607868 2 2 0.9033333 0.09428818 3 3 0.8686111 0.11314534 4 5 0.8750000 0.13289537 5 10 0.8797222 0.09513230 6 15 0.8250000 0.05257433 Since the 10-fold cross-validation approach is chosen, the mean accuracies as well as the corresponding standard deviations are returned.Thus, with Linear Discriminant Analysis (LDA), log-variances as characteristics, it seems that dimensions related to first and last q = 2 eigenvectors (2 × 2 dimensions in total) are enough to obtain a good classification, with an accuracy of 90%.Nevertheless, it is also observed that variation among folds can be relevant.
To visualize what is the representation in the reduced dimension space function plot can be used.For instance, to visualize the first unit of the first class, based on projections along the 2 first and last vectors (a 1 , a 2 and b 1 , b 2 ): plot (mydbcsp, index=1, class=1, vectors=1:2) In the top graphic of Figure 2 it can be seen what is the representation of the first video of class C 1 given by matrix X 11 , where the horizontal axis represents the frames of the video and the lines are the positions of the body key points (50 lines).In the bottom graphic, the same video is represented in a reduced space where the video is represented by the new signals (only 4 lines).
To have a better insight of the discriminating power of the new signals in the reduced dimension space, we can plot the corresponding log-variances of the new signals.
Parameter vectors in function boxplot sets which are the eigenvectors considered to plot.

boxplot(mydbcsp, vectors=1:2)
In Figure 3 it can be seen that variability of projections on the first eigenvector direction (log(V AR(X ′ ik a 1 ))) are big for elements in x1, but small for elements in x2.Analogously, projecting on the last dimension (log(V AR(X ′ ik b 1 ))), low variability is hold in x1 and big variability in x2.The same pattern holds when projecting on vectors a 2 and b 2 .

Basic/classic analysis new unit classification
Once the value of q has been decided and the accuracy of the classification is known, the classifier should be built (through train()) so that the user can proceed to predict the class a new action hold in a video belongs to, using the function predict.For instance, with only illustrative purpose, we can classify the first 5 videos which are stored in x1.In case it is known the label of the testing items, the latter function returns the accuracy.

Extending the example
In the previous section a basic workflow to use functions implemented in dbcsp is presented.Nevertheless, it is straightforward to extend the procedure.Once the interesting directions in W are calculated through dbcsp, other summarizing characteristics beyond the variance could be extracted from the projected signals, as well as other classifiers could be used in the classification step.Here it is shown how once the eigenvectors are extracted from an object dbcsp, several characteristics could be extracted from the signals and a new data.framecan be built so that any other classification technique could be applied.In this example we worked with caret package to apply different classifiers.It is important to pay attention to which the train and test sets are, so that the vectors are computed based only on training set instances.
Then, imagine we want to perform our classification step with the interquartile range information besides the log-variance.

Conclusions
In this work a new Distance-Based Common Spatial Pattern is introduced.It allows to perform the classical Common Spatial Pattern when the Euclidean distance between signals is considered, but it can be extended to the use of any other appropriate distance between signals as well.All of it is included in package dbcsp.The package is easy to use for non specialized users but, in the sake of flexibility, more advanced analysis can be carried out combining the created object and obtained results with already well known R packages, as caret, for instance.

Figure 1 :
Figure 1: The flow-chart shows how the consecutive steps of filtering, feature extraction and the construction of a classification model can be used to classify a new data.

•
Usage boxplot(mydbcsp) • Arguments object, an object of class dbcsp vectors, integer or vector of integers, indicating the index of the projected vectors to plot, by default index=1.pairs logical, if TRUE the pairs a j and b j of the indicated indices are also shown, by default pairs=TRUE.show_loglogical, if TRUE the logarithms of the variances are displayed, else the variances, by default show_log=TRUE.

6.
Look at: stare at the robot in front of it.There are 46 instances for this class.The data set is accessible via AR.data.Each class is a list of matrices of [K × num_f rames] dimensions, where K = 50 signals and num_f rames = 92 values.As mentioned before, the 50 signals represent the xy coordinates of 25 body key points extracted by OpenPose.For example, two different classes can be accessed this way: x1 <-AR.data$comex2 <-AR.data$fivewhere, x1 is a list of 46 instances of [50x92] matrices and x2 is a list of 45 instances of [50x92] matrices.First a basic/classic analysis is performed.

Figure 2 :
Figure 2: First video of class C 1 represented as its original version (top) and as the projections on vectors a 1 and a 2 (continuous lines) and b 1 and b 2 (dotted lines).

Figure 3 :
Figure 3: Log-variabilities of the projected signals on vectors a 1 and a 2 and b 1 and b 2 and separated by classes C 1 and C 2 .
fold = "integer", integer value, by default fold=10.It controls the number of partitions for the k-fold validation procedure, if the classification is done.seed = "numeric", numeric value, by default seed=NULL.Set a seed in case you want to be able to replicate the results.eig.tol= "numeric", numeric value, by default eig.tol=1e-06.If the minimum eigenvalue is below this tolerance, average covariance matrices are replaced by the most similar matrix that is definite positive.It is done via function nearPD in Matrix.