Research Article  Open Access
Multiple Subject Barycentric Discriminant Analysis (MUSUBADA): How to Assign Scans to Categories without Using Spatial Normalization
Abstract
We present a new discriminant analysis (DA) method called Multiple Subject Barycentric Discriminant Analysis (MUSUBADA) suited for analyzing fMRI data because it handles datasets with multiple participants that each provides different number of variables (i.e., voxels) that are themselves grouped into regions of interest (ROIs). Like DA, MUSUBADA (1) assigns observations to predefined categories, (2) gives factorial maps displaying observations and categories, and (3) optimally assigns observations to categories. MUSUBADA handles cases with more variables than observations and can project portions of the data table (e.g., subtables, which can represent participants or ROIs) on the factorial maps. Therefore MUSUBADA can analyze datasets with different voxel numbers per participant and, so does not require spatial normalization. MUSUBADA statistical inferences are implemented with crossvalidation techniques (e.g., jackknife and bootstrap), its performance is evaluated with confusion matrices (for fixed and random models) and represented with prediction, tolerance, and confidence intervals. We present an example where we predict the image categories (houses, shoes, chairs, and human, monkey, dog, faces,) of images watched by participants whose brains were scanned. This example corresponds to a DA question in which the data table is made of subtables (one per subject) and with more variables than observations.
1. Introduction
A standard problem in neuroimaging is to predict category membership from a scan. Called “brain reading” by Cox and Savoy [1], and more generally multivoxel pattern analysis (MVPA, for a comprehensive review see, e.g., [2]), this approach is used when we want to “guess” the type of category of stimuli processed when a participant was scanned and when we want to find the similarity structure of these stimulus categories (for a review, see, e.g., [3]). For datasets with the appropriate structure, this type of problem is addressed in multivariate analysis with discriminant analysis (DA). However, the structure of neuroimaging data precludes, in general, the use of DA. First, neuroimaging data often comprise more variables (e.g., voxels) than observations (e.g., scans). In addition, in the MVPA framework (see, e.g., the collection of articles reported in [2]), MRI data are collected as multiple scans per category of stimuli and the goal is to assign a particular scan to its category. These MRI data do not easily fit into the standard framework of DA, because DA assumes that one row is one observation (e.g., a scan or a participant) and one column is a variable (e.g., voxel). This corresponds to designs in which one participant is scanned multiple times or multiple participants are scanned once (assuming that the data are spatially normalized—e.g., put in Talairach space). These designs can fit PET or SPECT experiments but do not fit standard MRI experiments where, typically, multiple participants are scanned multiple times. In particular, DA cannot accommodate datasets with different numbers of variables per participant (a case that occurs when we do not use spatial normalization). Finally, statistical inference procedures of DA are limited by unrealistic assumptions, such as normality and homogeneity of variance and covariance.
In this paper, we present a new discriminant analysis method called Multiple Subject Barycentric Discriminant Analysis (MUSUBADA) which implements a DAlike approach suitable for neuroimaging data. Like standard discriminant analysis, MUSUBADA is used to assign observations to predefined categories and gives factorial maps in which observations and categories are represented as points, with observations being assigned to the closest category. But, unlike DA, MUSUBADA can handle datasets with multiple participants when each participant provides a different number of variables (e.g., voxels). Each participant is considered as a subtable of the whole data table and the data of one participant can also be further subdivided into more subtables which can constitute, for example, regions of interest (ROIs). MUSUBADA processes these subtables by projecting portions of the subtables on the factorial maps. Consequently, MUSUBADA does not require spatial normalization in order to handle “group analysis.” In addition, MUSUBADA can handle datasets with a number of variables (i.e., voxels) larger than the number of observations.
We illustrate MUSUBADA with an example in which we predict the type of images that participants’ were watching when they were scanned. For each participant, one anatomical ROI was used in the analysis. Because the ROIs were anatomically defined, the brain scans were not spatially normalized and the corresponding number of voxels, as well as their locations, were different for each participant.
We decided to use handdrawn ROIs, we could also have use functional localizers as these are still widely used (and constitute a valid approach as long as the localizer is not confounded with the experimental tasks). Hand drawing these ROIs is obviously very time consuming and could restrict the use of our technique to only small studies (or to very dedicated teams) and therefore could also make studies harder to replicate. Fortunately, there are ways of obtaining ROIs without drawing them by hand. Specifically, as an alternative to manual tracing and functional localizers, recent methods have been developed to choose ROIs for each subject which are both automatic and a priori. These methods take labels from a standard anatomical atlas such as AAL [4], Tailarach [5], or Brodman [6] and warp these labels to coordinates within each subject’s anatomical space. The anatomical coordinates are then downsampled to the subject’s functional space. These steps can either be completely automated via extensions to standard software, such as IBASPM [7] or by reusing builtin tools, such as the linear FLIRT [8] and nonlinear FNIRT [9] of FSL. Because standard singlesubject atlases do not account for betweensubject variation [10], it may be preferable to use probabilistic atlases determined on multiple subjects (e.g., [11]). As an alternative to anatomical atlases entirely, stereotactic coordinates can also be taken from a metaanalysis and warped into coordinates within the subject’s functional space. Although metaanalyses are generally performed for the task at hand, methods exist for automating even metaanalyses using keywords in published articles (see, e.g., NeuroSynth: [12]).
1.1. Overview of the Method
MUSUBADA comprises two steps: barycentric discriminant analysis (BADA) analyzes a data table in which observations (i.e., scans) are rows and in which variables (i.e., voxels) are columns and where each participant is represented by a subset of the voxels (i.e., one participant is a “subtable” of the whole data table), and projection of the subtables representing the participants on the solution computed by BADA (this is the “MUSU” step in MUSUBADA). In addition, the subtable representing one participant could also be further subdivided into subtables representing, for example, the participant’s ROIs (the ROIs can differ with the participants).
BADA generalizes discriminant analysis and, like DA, it is performed when measurements made on some observations are combined to assign observations to apriori defined categories. BADA is, actually, a class of methods which all rely on the same principle: each category of interest is represented by the barycenter of its observations (i.e., the weighted average; the barycenter is also called the center of gravity of the observations of a given category), and, then, a generalized principal component analysis (GPCA) is performed on the category by variable matrix. This analysis gives a set of discriminant factor scores for the categories and another set of factor scores for the variables. The original observations are then projected onto the category factor space, providing a set of factor scores for the observations. The distance of each observation to the set of categories is computed from the factor scores and each observation is assigned to the closest category.
The comparison between the apriori and aposteriori category assignments is used to assess the quality of the discriminant procedure. When the quality of the model is evaluated for the observations used to build the model, we have a fixed effect model. When we want to estimate the performance of the model for new or future observations, we have a random effect model. In order to estimate the quality of the random effect model, the analysis is performed on a subset of the observations called the training set and the predictive performance is evaluated with a different set of observations called the testing set. A specific case of this approach is the “leaveoneout” technique (also called jackknife) in which each observation is used, in turn, as the testing set whereas the rest of the observations play the role of the training set. This scheme has the advantage of providing an approximately unbiased estimate of the generalization performance of the model [13]. The quality of the discrimination can also be evaluated with an type statistic which expresses the proportion of the data variance explained by the model. Its significance can be evaluated with standard permutation tests.
The stability of the discriminant model can be assessed by a resampling strategy such as the bootstrap (see [13, 14]). In this procedure, multiple sets of observations are generated by sampling with replacement from the original set of observations, and by computing new category barycenters, which are then projected onto the original discriminant factor scores. For convenience, the confidence intervals of the barycenters can be represented graphically as a confidence ellipsoid that encompasses a given proportion (say 95%) of the barycenters. When two category ellipsoids do not intersect, these groups are significantly different.
The problem of multiple tables corresponds to MUSUBADA per se and it is implemented after the BADA step. In the MUSUBADA step, each subtable is projected onto the factor scores computed for the whole data table. These projections are also barycentric as their average gives the factor scores of the whole table. This last step integrates other multitable techniques such as multiple factor analysis or STATIS [15–19] which have also been used in brain imaging (see, e.g., for recent examples [20–22]). In addition to providing subtable factor scores, MUSUBADA evaluates and represents graphically the importance (often called the contribution) of each subtable to the overall discrimination. A sketch of the main steps of MUSUBADA is shown in Figure 1.
MUSUBADA incorporates BADA which, itself, is a GPCA performed on the category barycenters and as such MUSUBADA implements a discriminant analysis version of different multivariate techniques such as correspondence analysis, biplot analysis, Hellinger distance analysis, and canonical variate analysis (see, e.g., [23–26]). In fact, for each specific type of GPCA, there is a corresponding version of BADA. For example, when the GPCA is correspondence analysis, this gives the most wellknown version of BADA: discriminant correspondence analysis (DICA, sometimes also called correspondence discriminant analysis; see [23, 27–31]).
2. Notations
Matrices are denoted with bold uppercase letters (i.e., ) with generic element denoted with the corresponding lowercase italic letter (i.e., ). The identity matrix is denoted . Vectors are denoted with bold lowercase letter (i.e., ) with generic element denoted with the corresponding lower case italic (i.e., ).
The original data matrix is an observation by variables matrix denoted . Prior to the analysis, the matrix can be preprocessed by centering (i.e., subtracting the column mean from each column), by transforming each column into a score, or by normalizing each row so that the sum of its elements or the sum of its squared elements is equal to one (the rationale behind these different types of normalization is discussed later on). The observations in are partitioned into apriori categories of interest with being the number of observations of the th category (and so . The columns of matrix can be arranged in a priori subtables. The numbers of columns of the th subtable are denoted (and so . So, the matrix can be decomposed into by blocks as
2.1. Notations for the Categories (Rows)
We denote by the by design (aka dummy) matrix for the categories describing the rows of : when row belongs to category and , , otherwise. We denote by the by 1 vector of masses for the rows of and by the by diagonal matrix whose diagonal elements are the elements of (i.e., using the diag operator which transforms a vector into a diagonal matrix, we have ). Masses are positive numbers and it is convenient (but not necessary) to have the sum of the masses equal to one. The default value for the mass of each observation is often . We denote by the by 1 vector of masses for the categories describing the rows of and by the by diagonal matrix whose diagonal elements are the elements of .
2.2. Notations for the Subtables (Columns)
We denote by the by design matrix for the subtables from the columns of : if column belongs to subtable , , otherwise. We denote by the by 1 vector of weights for the columns of and by the by diagonal matrix whose diagonal elements are the elements of . We denote by the by 1 vector of weights for the subtables of and by the by diagonal matrix whose diagonal elements are the elements of . The default value for the weight of each variable is , a more general case requires only to be positive definite and this includes nondiagonal weight matrices.
3. Barycentric Discriminant Analysis (BADA)
The first step of BADA is to compute the barycenter of each of the categories describing the rows. The barycenter of a category is the weighted average of the rows in which the weights are the masses rescaled such that the sum of the weights for each category is equal to one. Specifically, the by matrix of barycenters, denoted , is computed as where is an by 1 vector of 1 s and the diagonal matrix rescales the masses of the rows such that their sum is equal to one for each category.
3.1. Masses and Weights
The type of preprocessing and the choice of the matrix of masses for the categories () and the matrix of weights for the variables () is crucial because these choices determine the type of GPCA used.
For example, discriminant correspondence analysis is used when the data are counts. In this case, the preprocessing is obtained by transforming the rows of into relative frequencies, and by using the relative frequencies of the barycenters as the masses of the rows and the inverse of the column frequencies as the weights of the variables. Another example of GPCA, standard discriminant analysis, is obtained when is equal to the inverse of the within group variancecovariance matrix (which can be computed only when this matrix is full rank). Hellinger distance analysis (also called “spherical factorial analysis”; [32–35]) is obtained by taking the square root of the relative frequencies for the rows of and by using equal weights and masses for the matrices and . Interestingly, the choice of weight matrix is equivalent to defining a generalized Euclidean distance between dimensional vectors [32]. Specifically, if and are two dimensional vectors, the generalized Euclidean squared distance between these two vectors is
3.2. GPCA of the Barycenter Matrix
Essentially, BADA boils down to a GPCA of the barycenter matrix under the constraints provided by the matrices (for the categories) and (for the columns). Specifically, the GPCA is implemented by performing a generalized singular value decomposition of matrix [23, 26, 36, 37], which is expressed as where is the by diagonal matrix of the singular values (with being the number of nonzero singular values), and (resp., ) being the by (resp., by ) matrix of the left (resp., right) generalized singular vectors of .
3.3. Factor Scores
The by matrix of factor scores for the categories is obtained as These factor scores are the projections of the categories on the GPCA space and they provide the best separation between the categories because they have the largest possible variance. In order to show this property, recall that the variance of the columns of is given by the square of the corresponding singular values (i.e., the “eigenvalue” denoted ) and are stored in the diagonal matrix ). This can be shown by combining (4) and (5) to give Because the singular vectors of the SVD are ordered by size, the first factor extracts the largest possible variance, the second factor extracts the largest variance left after the first factor has been extracted, and so forth.
3.3.1. Supplementary Elements
The rows of matrix can be projected (as “supplementary” or “illustrative” elements) onto the space defined by the factor scores of the barycenters. Note that the matrix from (5) is a projection matrix. Therefore, the by matrix of the factor scores for the rows of can be computed as These projections are barycentric, which means that the weighted average of the factor scores of the rows of a category gives the factors scores of this category. This property is shown by first computing the barycenters of the row factor scores as (cf. (2)) as then plugging in (7) and developing. Taking this into account, (5) gives
3.4. Loadings
The loadings describe the variables of the barycentric data matrix and are used to identify the variables important for the separation between the categories. As for standard PCA, there are several ways of defining the loadings. The loadings can be defined as the correlation between the columns of matrix and the factor scores. They can also be defined as the matrix or as variable “factor scores” which are computed as (Note that and differ only by a normalizing factor).
4. Quality of the Prediction
The performance, or quality of the prediction of a discriminant analysis, is assessed by predicting the category membership of the observations and by comparing the predicted with the actual category membership. The pattern of correct and incorrect classifications can be stored in a confusion matrix in which the columns represent the actual categories and in which the rows represent the predicted categories. At the intersection of a row and a column is the number of observations from the column category assigned to the row category.
The performance of the model can be assessed for the observations (e.g., scans or participants) actually used to compute the categories (the set of observations used to generate the model is sometimes called the training set). In this case, the performance of the model corresponds to a fixed effect model because this assumes that a replication of the experiment would use the same observations (i.e., the same participants and the same stimuli). In order to assess the quality of the model for new observations, its performance, however, needs to be evaluated using observations that were not used to generate the model (the set of “new observations” used to evaluate the model is sometimes called the testing set). In this case, the performance of the model corresponds to a random effect model because this assumes that a replication of the experiment would use the different observations (i.e., different participants and stimuli).
4.1. Fixed Effect Model
The observations in the fixed effect model are used to compute the barycenters of the categories. In order to assign an observation to a category, the first step is to compute the distance between this observation and all categories. Then, the observation is assigned to the closest category. Several possible distances can be chosen, but a natural choice is the Euclidean distance computed in the factor space. If we denote by the vector of factor scores for the th observation, and by the vector of factor scores for the th category, then the squared Euclidean distance between the th observation and the th category is computed as Obviously, other distances are possible (e.g., Mahalanobis distance), but the Euclidean distance has the advantage of being “directly read” on the map.
4.1.1. Tolerance Intervals
The quality of the category assignment of the actual observations can be displayed using tolerance intervals. A tolerance interval encompasses a given proportion of a sample or a population. When displayed in two dimensions, these intervals have the shape of an ellipse and are called tolerance ellipsoids. For BADA, a category tolerance ellipsoid is plotted on the category factor score map. This ellipsoid is obtained by fitting an ellipse which includes a given percentage (e.g., 95%) of the observations. Tolerance ellipsoids are centered on their categories. The overlap of the tolerance ellipsoids of two categories reflects the proportion of misclassifications between these two categories for the fixed effect model.
4.2. Random Effect Model
The observations of the random effect model are not used to compute the barycenters but are used only to evaluate the quality of the assignment of new observations to categories. A convenient variation of this approach is “leaveoneout” (aka jackknife) approach: Each observation is taken out from the dataset and, in turn, is then projected onto the factor space of the remaining observations in order to predict its category membership. For the estimation to be unbiased, the leftout observation should not be used in any way in the analysis. In particular, if the data matrix is preprocessed, the leftout observation should not be used in the preprocessing. So, for example, if the columns of the data matrix are transformed into scores, the leftout observation should not be used to compute the means and standard deviations of the columns of the matrix to be analyzed, but these means and standard deviations will be used to compute the score for the leftout observation.
The assignment of a new observation to a category follows the same procedure as for an observation from the fixed effect model. The observation is projected onto the original category factor scores and is assigned to the closest category. Specifically, we denote by the data matrix without the th observation, and by the 1 by row vector representing the th observation. If is preprocessed (e.g., centered and normalized), the preprocessing parameters will be estimated without (e.g., the mean and standard deviation of are computed without ) and will be preprocessed with the parameters estimated for (e.g., will be centered and normalized using the means and standard deviations of the columns of ). Then, the matrix of barycenters is computed and its generalized eigendecomposition is obtained as (cf. (4)) (with and being the mass and weight matrices for ). The matrix of factor scores denoted is obtained as (cf. (5)) The projection of the th observation, considered as the “testing” or “new observation,” is denoted and it is obtained as (cf. (7)) Distances between this th observation and the categories can be computed from the factor scores (cf. (11)). The observation is then assigned to the closest category. In addition, the jackknife approach can provide an (unbiased) estimate of the position barycenters as well as their standard error (see, e.g., [38], for this approach).
Often in MRI experiments, observations are structured in blocks in which observations are not independent of each others (this is the case in most “block designs”). In such cases, a standard leaveoneout approach will overestimate the quality of prediction and should be replaced by a “leave one block out” procedure.
4.2.1. Prediction Intervals
In order to display the quality of the prediction for new observations, we use prediction intervals. Recall that a “leave one out” or jackknife (or “leave one block out”) procedure is used to predict each observation from the other observations. In order to compute prediction intervals, the first step is to project the leftout observations onto the original complete factor space. There are several ways to project a leftout observation onto the factor score space. Here, we propose a twostep procedure. First, the observation is projected onto the factor space of the remaining observations. This provides factor scores for the leftout observation and these factor scores are used to reconstructed the observation from its projection (in general, the leftout observation is imperfectly reconstructed and the difference between the observation and its reconstruction reflects the lack of fit of the model). Then, the reconstructed observation, denoted , is projected onto the full factor score solution. Specifically, a leftout observation is reconstructed from its factor scores as (cf. (4) and (14)) The projection of the leftout observation is denoted and is obtained by projecting as a supplementary element in the original solution. Specifically, is computed as Prediction ellipsoids are not necessarily centered on their categories (the distance between the center of the ellipse and the category represents the estimation bias). Overlap of two predictions intervals directly reflects the proportion of misclassifications for the new observations.
5. Quality of the Category Separation
5.1. Explained Inertia () and Permutation Test
In order to evaluate the quality of the discriminant model, we use a coefficient inspired by the coefficient of correlation. Because BADA is a barycentric technique, the total inertia (i.e., the “variance”) of the observations to the grand barycenter (i.e., the barycenter of all categories) can be decomposed into two additive quantities: the inertia of the observations relative to the barycenter of their own category, and the inertia of the category barycenters to the grand barycenter. Specifically, if we denote by the vector of the coordinates of the grand barycenter (i.e., each component of this vector is the average of the corresponding components of the barycenters), the total inertia, denoted , is computed as the sum of the squared distances of the observations to the grand barycenter (cf. (11)): The inertia of the observations relative to the barycenter of their own category is abbreviated as the “inertia within.” It is denoted and computed as The inertia of the barycenters to the grand barycenter is abbreviated as the “inertia between.” It is denoted and computed as So the additive decomposition of the inertia can be expressed as This decomposition is similar to the familiar decomposition of the sum of squares in the analysis of variance. This suggests that the intensity of the discriminant model can be tested by the ratio of between inertia by the total inertia, as is done in analysis of variance and regression. This ratio is denoted and it is computed as The ratio takes values between 0 and 1, the closer to one, the better the model. The significance of can be assessed by permutation tests and confidence intervals can be computed using crossvalidation techniques such as the jackknife (see [19]).
5.2. Confidence Intervals
The stability of the position of the categories can be displayed using confidence intervals. A confidence interval reflects the variability of a population parameter or its estimate. In two dimensions, this interval becomes a confidence ellipsoid. The problem of estimating the variability of the position of the categories cannot, in general, be solved analytically and crossvalidation techniques need to be used. Specifically, the variability of the position of the categories is estimated by generating bootstrapped samples from the sample of observations. A bootstrapped sample is obtained by sampling with replacement from the observations. The “bootstrapped barycenters” obtained from these samples are then projected onto the original discriminant factor space and, finally, an ellipse is plotted such that it comprises a given percentage (e.g., 95%) of these bootstrapped barycenters. When the confidence intervals of two categories do not overlap, these two categories are “significantly different” at the corresponding alpha level (e.g., ).
It is important to take into account the structure of the design when implementing a bootstrap scheme because the bootstrap provides an unbiased estimate only when the bootstrapped observations are independent [39]. Therefore, when the observations are structured into subtables, these subtables should be bootstrapped in addition to the observations. Conversely, the fixed part of a design should be kept fixed. So, for example, when scans are organized into blocks, the block structure should be considered as fixed [40].
5.2.1. Interpreting Overlapping Multidimensional Confidence Intervals
When a confidence interval involves only one dimension (e.g., when using a confidence interval to compare the mean of two categories), the relationship between hypothesis testing and confidence intervals is straightforward. If two confidence intervals do not overlap, then the null hypothesis is rejected. Conversely, if two confidence intervals overlap, then the null hypothesis cannot be rejected. The same simple relationship holds with a 2dimensional display as long as all the variance of the data can be described with only two dimensions. If two confidence ellipses (i.e., the 2dimensional expression of an interval) do not overlap, then the null hypothesis is rejected, whereas if the ellipses overlap, then the null hypothesis cannot be rejected.
In most multivariate analyses (such as BADA), however, the 2dimensional maps used to display the data represent only part of the variance of the data and these displays can potentially be misleading because the position of a confidence ellipse in a 2dimensional map gives only an approximation of the real position of the intervals in the complete space. Now, if two confidence ellipses do not overlap in at least one display, then the two corresponding categories do not overlap in the whole space (and the null hypothesis can be rejected). However, when two confidence ellipses overlap in a given display, then the two categories may or may not overlap in the whole space (because the overlap may be due to a projection artifact). In this case, the null hypothesis may or may not be rejected depending upon the relative position of the ellipses in the other dimensions of the space. Strictly speaking, when analyzing data laying in a multidimensional space, the interpretations of confidence intervals are correct only when performed in the whole space and the 2dimensional representations give only an approximation. As a palliative, it is important to examine additional graphs obtained from other dimensions than the first two.
5.2.2. Confidence Intervals: The Multiple Comparison Problem
As mentioned earlier, a confidence interval generalizes a null hypothesis test. And, just like a standard test, the level chosen is correct only when there are only two categories to be compared because the problem of the inflation of Type I error occurs when there are more than two categories. A possible solution to this problem is to use a Bonferonni or a Šidàk correction (see, [41] for more details). Specifically, for categories, a Bonferonnicorrected confidence interval at the overall level is obtained as Along the same lines, a Šidàkcorrected confidence level for all pairs of comparisons is expressed as
6. Multiple Subject Barycentric Discriminant Analysis (MUSUBADA)
In a multitable analysis, the subtables can be analyzed by projecting the categories and the observations for each subtable. As was the case for the categories, these projections are barycentric because the barycenter of the all the subtables gives the coordinates of the whole table.
6.1. Partial Projection
Each subtable can be projected in the common solution. The procedure starts by rewriting (4) as where is the th subtable (comprising the columns of corresponding to the columns of the th block). Then, to get the projection for the th subtable, (5) is rewritten as where is the weight matrix for the columns of the th block.
Equation (25) can also be used to project supplementary rows corresponding to a specific subtable. Specifically, if is a row vector of a supplementary element to be projected according to the th subtable, its factor scores are computed as ( is supposed to have been preprocessed as )
6.2. Inertia of a Subtable
Recall from (6) that, for a given dimension, the variance of the factor scores of all the columns of matrix is equal to the eigenvalue of this dimension. Each subtable comprises a set of columns, and the contribution of a subtable to a dimension is defined as the sum of this dimension squared factor scores of the columns comprising this subtable. Precisely, the inertia for the th table and the th dimension is computed as Note that the sum of the inertia of the blocks gives back the total inertia:
6.2.1. Finding the Important Subtables
The subtables that contribute to the discrimination between the classes are identified from their partial contributions to the inertia (see (27)).
6.3. Subtable Normalization and Relationship with Other Approaches
In the current example, the subtables are simply considered as sets of variables. Therefore, the influence of a given subtable reflects not only its number of variables, but also its factorial structure because, everything else being equal, a subtable with a large first eigenvalue will have a large influence on the first factor of the whole table. By contrast, a subtable with a small first eigenvalue will have a small effect on the first factor of the whole table. In order to eliminate such discrepancies, the multiple factor analysis approach [15, 17, 36] normalizes each subtable by dividing each element of a subtable by its first singular value.
An alternative normalization can be derived from the STATIS approach [16, 18, 19, 42, 43]. In this framework, each subtable is normalized from an analysis of the by matrix of the betweensubtable matrix (recall that the coefficient plays a role analogous to the coefficient of correlation for semipositive definite matrices, see [41, 44, 45]). In this framework, subtables are normalized in order to reflect how much information they share with the other subtables.
In the context of MUSUBADA, these normalization schemes can be performed on the subtables of the matrix , but their effects are easier to analyze if it is performed on the barycentric matrix . These two normalizing schemes can also be combined with the subtables being first normalized according the multiple factor analysis approach and then using the STATIS approach.
6.4. Programs
MATLAB and R programs are available from the home page of the first author (at http://www.utdallas.edu/~herve/).
6.5. MUSUBADA and Other Classifiers
MUSUBADA is a linear classifier that can handle multiple participants when these participants are represented by different numbers of voxels or ROIs. As such, its performance is likely to be closely related to the performance of other linear classifiers (see [46] for a general presentation in the context of MVPA models) such as (linear) support vector machines (SVM, see, e.g., [47] for an application of SVM to MRI data) or linear discriminant analysis (in general performed on the components of a PCA, see [48]). When the structure of the data (i.e., spatial normalization) allows the use of different linear techniques, these techniques will likely perform similarly (numerical simulations performed with some examples confirm this conjecture). However, as indicated previously, most of these techniques will not be able to handle the integration of participants with different number of voxels or ROIs. Canonical STATIS (CANOSTATIS [49]) can integrate discriminant analysis problems obtained with a different number of variables, but in its current implementation, it requires that the data of each participant be full rank, a restriction that makes difficult for this technique to be used with large brain imaging datasets. It has been suggested [19] to use approaches such as ridge or regularization techniques (see, e.g., [50], for a review) to make CANOSTATIS able to handle multicolinearity, and such a development could be of interest for brain imaging. However, CANOSTATIS works on Mahanalobis distance matrices between categories and does not have—in its current implementation—ways of predicting membership of scans to categories and cannot identify the voxels important for the discrimination between categories.
The closest statistical technique to MUSUBADA is probably mean centered partial leastsquare correlation (MCPLSC, see [51–53], for presentations and reviews), which is widely used in brain imaging. In fact, when the data are spatially normalized, MUSUBADA and MCPLSC will analyze (i.e., compute the singular value decomposition of) the same matrix of barycenters. MUSUBADA adds to MCPLSC the possibility of handling different numbers of voxels per participant, ROI analysis, and also an explicit prediction of group categories which is lacking from the standard implementation of MCPLSC. So, MUSUBADA can be seen as a generalization of MCPLSC that can handle multiple ROIs, different numbers of voxels per participant, and can predict category membership.
7. An Example
As an illustration, we are using the dataset originally collected by Connolly et al. [54]. In this experiment, 10 participants performed a oneback recognition memory task on visual stimuli from seven different categories: female faces, male faces, monkey faces, dog faces, houses, chairs, and shoes. The stimuli were presented in 8 runs of 7 blocks—one block per category per run. A block consisted of 16 image presentations from a given category with each image presented for a duration of 500 ms separated by a 2 s interstimulus interval (2.5 s stimulus onset asynchrony) for a total block duration of 40 s. Full brain scans were acquired every 2 s (TR = 2000 ms). Blocks were separated by 12 s intervals of fixation. A different random ordering of category blocks was assigned to each run, and each subject saw the same set of runs (i.e., the same set of random category block orderings), however, each subject saw the runs in a different randomly assigned run order. For analysis, the runs for each subject were reordered to a canonical run order so that a single set of timepoint labels could be used for all subjects. Timepoint labels were coded such that the first four TRs of each block were discarded so that only the maximal evoked BOLD response was coded for each category block. The timecourse for each run was thus coded with 16 consecutive fullbrain scans (covering 32 sec starting 8 secs after the onset of the block) assigned to each stimulus block. No regression analyses were performed to estimate the average BOLD response for the stimulus categories; rather, each brain image (16 per stimulus block) contributed to a unique row in the data matrix.
7.1. Data Acquisition
BOLD MRI images were acquired with gradient echoplanar imaging using a Siemens Allegra head only 3T scanner. The MRI images consisted of thirtytwo 3 mm thick axial images (TR = 2000, TE = 30 ms, flip angle = 90, 64 × 64 matrix, FOV = 192 × 192 mm) and included all of the occipital and temporal lobes and the dorsal parts of the frontal and parietal lobes. Highresolution T1weighted structural images were also acquired.
7.2. Imaging Preprocessing
Preprocessing of the MRI data included slice timing correction, volume registration to correct for minor head movements, correction of extreme values, and mean correction for each run. No spatial normalization or coregistration was performed.
7.3. Region of Interest
For each participant, a mask was hand drawn based on anatomical landmarks obtained from the structural scan. The mask was used to identify one ROI the Ventrotemporal (VT) area ROI which included the inferior temporal, fusiform, and parahippocampal gyri. Because the ROI was anatomically defined, its size (in number of voxels) differed from participant to participant (from 2791 to 4815). The total number of voxels (for all 10 participants) was equal to 39,163.
7.4. Statistical Preprocessing
Prior to the analysis, the data were centered by removing the average scan. In addition, a subtable normalization was performed on each participant. For each participant, a singular value decomposition (SVD) was performed and all the voxels of a given participant were divided by the first singular value obtained from this SVD. In addition, each row of the data table was normalized such that the sum of the squares of its elements was equal to one.
7.5. Final Data Matrices
The data matrix, , is a 896 × 39,163 matrix. The rows of correspond to the 896 scans which were organized in 8 blocks of 7 stimulus categories comprising 16 brain images each. Therefore, each of the seven stimulus groups was represented by 16 × 8 = 128 brain images. The 39,163 columns of represented the voxel activations of all the participants and are organized in 10 subtables (one per participant; the data could be fitted in such a simple format because both the scan block orders were constant for all participants). The design matrix was a 896 × 7 matrix coding for the categories. A schematic of the data matrix is shown in Figure 2.
7.6. Goal: Finding the Categories
The goal of the analysis was to find out if it was possible to discriminate between scans coming from the different categories of images. To do this, each category in was summed to create the 7 × 39,163 barycentric matrix (see (2)). We then computed the BADA on the matrix which provided the map shown in Figure 3. The first dimension of this map shows a clear opposition of faces (i.e., female, male, and dog faces) and objects (i.e., houses, chairs, and shoes). This dimension can be interpreted as a semantic dimension: faces versus nonfaces. The second dimension, in contrast, separates the smaller from the larger objects (i.e., chairs and shoes versus houses) and also reflects variation among the faces, as it separates monkey from human faces (with dog faces inbetween). It is also possible that dimension 2 reflects a combination of features (e.g., low level visual features). In order to support this interpretation, we ran a BADA on the stimuli (see, e.g., [55] for a similar idea) and obtained a set of factor scores for the picture groups. We then projected the picture categories, as supplementary elements, onto the factor space obtained from the MRI data. We used the procedure described in Abdi [56] after rescaling the picture factor scores so that the factor scores for the first dimension have the same variance (i.e., eigenvalue) as the first factor scores of the analysis performed on the participants. The results are shown in Figure 4. In agreement with our interpretation, we found that the categories obtained from the pictures lined up with dimension 2 and had very little variance on dimension 1.
(a)
(b)
7.7. Stability of the Discriminant Model
7.7.1. Fixed Effect Model
To examine the reliability of the fixed effect model, we looked at the fixed effect confusion matrix and the tolerance intervals. The fixed effect confusion matrix is shown in Table 1. We see that, for the fixed effect model, classification was very good: the model correctly classified 126/128 female faces, 127/128 male faces, all 128 monkey faces, 127/128 houses, all 128 chairs, 124/128 shoes, and all 128 dog faces.

To answer the question of separability, we looked at the tolerance intervals (shown in Figure 5(a)). The tolerance intervals include 95% of the projections of the observations centered around their respective categories. Recall that tolerance intervals reflect the accuracy of the assignment of scans to categories. Therefore, when two categories do not overlap for at least one dimension, they can be considered as separable (note, that, again all dimensions need to be considered and that a twodimensional display is accurate only when all the variance of the projections is two dimensional). We show here the tolerance intervals only for dimensions 1 and 2, but we based our interpretation on the whole space.
(a)
(b)
Large differences were found between houses and female faces, between houses and male faces, between houses and dog faces, between chairs and female faces, between chairs and male faces, between chairs and monkey faces, and between chairs and dog faces. The variance of the projection of the individual scans (as represented by the area of the ellipsoids) varies with the categories with monkey faces and shoes displaying the largest variance and with human faces showing the smallest variance.
Interestingly, the female and male faces are close to each other and overlap on most dimensions, this pattern suggests that the scans from these categories were very similar.
7.7.2. Random Effect Model
The random effect model evaluates the performance of the model for new observations. The current experiment used a block design, and the scans within a block are generally not independent because of the large time correlation typical of MRI experiments, whereas scans between blocks can be considered independent (because of the resting interval between blocks). Therefore, in order to compute the random effect confusion matrix, we used a “leaveoneblockout” procedure in which we left the whole block out in order to avoid confounding prediction with time correlation (see [38]). This resulted in the confusion matrix shown in Table 2. As expected, classification for the random effect model was not as good as that for the fixed effect model, with 47/128 female faces, 34/128 male faces, 53/128 monkey faces, 81/128 houses, 96/128 chairs, 45/128 shoes, and 53/128 dog faces correctly classified. So, overall the categories are well separated in the random effect model. However, the human females and male faces were not separated and seem to constitute a single category (with female faces showing more variability the male faces).

Like the fixed effect model, we can examine the stability of the random effect model. The prediction intervals include 95% of the bootstrap sample projections. Note that prediction intervals may not be centered around the category mean, because new observations are unlikely to have the exact same mean as the old observations. Recall that prediction intervals reflect the accuracy of assignment of new scans to categories and represent where observations would fall in the population. The prediction intervals for the example are shown in Figure 5(b). Note that the large overlap of the categories ellipsoids suggests that some scans of a given categories can be misclassified for scans from any other category.
7.8. Quality of Category Separation
7.8.1. and the Permutation Test
The quality of the model was also evaluated by computing the percentage of variance explained by the model. This gave a value of with a probability of (by permutation test) which confirms the quality of the model.
7.8.2. Confidence Intervals
To test whether category discrimination was significant, we used 95% confidence intervals. Figure 6shows the 95% confidence ellipses for each category on the maps made from Dimensions 12 and 34. The configuration of the ellipses indicates that all categories are reliably separated, but the separation of some categories is stronger than some others. Specifically, the male and female faces are separated only on the 34 map and their distance is relatively small, and this indicates that their separation—even if significant—is small.
(a)
(b)
7.9. Subtable Projections
7.9.1. Partial Inertias of the Participants
The respective importance of the participants is expressed by the partial inertia (cf. (27)) of their groups of voxels (i.e., the subtable associated to each subject). The partial inertias are plotted in Figure 7. We can see that there is some difference between the participants in the way they express the effect. For example, participant 7 shows the least overall separation between categories, participant 9 shows the clearest overall separation between categories.
7.9.2. Projecting the Subtables as Supplementary Elements
Interestingly, the overall factorial configuration is very similar for all participants, but some participants exhibit clearer configurations. In order to illustrate this point, in Figure 8, we show the map of participant 7 who shows the least overall separation between categories, participant 9 who shows the clearest overall separation between categories, and participant 6 who shows the clearest separation between categories on dimension 1 (see Figure 7).
(a)
(b)
(c)
8. Conclusion
Multiple subject barycentric discriminant analysis is particularly well suited for the analysis of neuroimaging data because it does not require brains to be spatially normalized. In addition, MUSUBADA can handle very large data sets (with more variables than observations). Even though we did not illustrate this property here, with MUSUBADA, we can also plot (i.e., with “glass brains”) the voxels important for discriminating between categories. Therefore, MUSUBADA can be used to study the similarity structure of the categories as well as that of the voxels. In addition, MUSUBADA could accommodate more sophisticated designs than the one illustrated here. For example, MUSUBADA could analyze design for which several ROIs are defined per participant (see, e.g., [57]), or different categories of participants (e.g., old versus young participants as in [58]). Finally, because MUSUBADA incorporates inferential components, it complements other popular approaches such as partial least squares [53, 59] which are widely used for brain imaging data.
References
 D. D. Cox and R. L. Savoy, “Functional magnetic resonance imaging (fMRI) “brain reading”: detecting and classifying distributed patterns of fMRI activity in human visual cortex,” NeuroImage, vol. 19, no. 2, pp. 261–270, 2003. View at: Publisher Site  Google Scholar
 N. Kriegestkorte and G. Kreiman, Visual Population Codes: Towards a Common Multivariate Framework for Cell Recording and Functional Imaging, MIT Press, Cambridge, Mass, USA, 2012.
 A. J. O'Toole, F. Jiang, H. Abdi, N. Pénard, J. P. Dunlop, and M. A. Parent, “Theoretical, statistical, and practical perspectives on patternbased classification approaches to the analysis of functional neuroimaging data,” Journal of Cognitive Neuroscience, vol. 19, no. 11, pp. 1735–1752, 2007. View at: Publisher Site  Google Scholar
 N. TzourioMazoyer, B. Landeau, D. Papathanassiou et al., “Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI singlesubject brain,” NeuroImage, vol. 15, no. 1, pp. 273–289, 2002. View at: Publisher Site  Google Scholar
 J. L. Lancaster, D. TordesillasGutiérrez, M. Martinez et al., “Bias between MNI and talairach coordinates analyzed using the ICBM152 brain template,” Human Brain Mapping, vol. 28, no. 11, pp. 1194–1205, 2007. View at: Publisher Site  Google Scholar
 S. B. Eickhoff, S. Heim, K. Zilles, and K. Amunts, “Testing anatomically specified hypotheses in functional imaging using cytoarchitectonic maps,” NeuroImage, vol. 32, no. 2, pp. 570–582, 2006. View at: Publisher Site  Google Scholar
 Y. AlemánGómez, L. MelieGarcná, and P. ValdésHernandez, “IBASPmml: toolbox for automatic parcellation of brain structures,” in Proceedings of the 12th Annual Meeting of the Organization for Human Brain Mapping, vol. 27, NeuroImage, Florence, Italy, 2006. View at: Google Scholar
 M. Jenkinson and S. Smith, “A global optimisation method for robust affine registration of brain images,” Medical Image Analysis, vol. 5, no. 2, pp. 143–156, 2001. View at: Publisher Site  Google Scholar
 A. Klein, J. Andersson, B. A. Ardekani et al., “Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration,” NeuroImage, vol. 46, no. 3, pp. 786–802, 2009. View at: Publisher Site  Google Scholar
 R. A. Poldrack, “Region of interest analysis for fMRI,” Social Cognitive and Affective Neuroscience, vol. 2, no. 1, pp. 67–70, 2007. View at: Publisher Site  Google Scholar
 A. Hammers, R. Allom, M. J. Koepp et al., “Threedimensional maximum probability atlas of the human brain, with particular reference to the temporal lobe,” Human Brain Mapping, vol. 19, no. 4, pp. 224–247, 2003. View at: Publisher Site  Google Scholar
 T. Yarkoni, R. A. Poldrack, T. E. Nichols, D. C. Van Essen, and T. D. Wager, “Largescale automated synthesis of human functional neuroimaging data,” Nature Methods, vol. 8, no. 8, pp. 665–670, 2011. View at: Publisher Site  Google Scholar
 T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Springer, New York, NY, USA, 2nd edition, 2008.
 B. Efron and R. Tibshirani, An Introduction to the Bootstrap, Champman & Hall, New York, NY, USA, 1993.
 B. Escofier and J. Pagès, “Multiple factor analysis (AFMULT package),” Computational Statistics and Data Analysis, vol. 18, no. 1, pp. 121–140, 1994. View at: Google Scholar
 C. Lavit, Y. Escoufier, R. Sabatier, and P. Traissac, “The ACT (STATIS method),” Computational Statistics and Data Analysis, vol. 18, no. 1, pp. 97–119, 1994. View at: Google Scholar
 H. Abdi and D. Valentin, “Multiple factor analysis,” in Encyclopedia of Measurement and Statistics, N. Salkind., Ed., pp. 657–663, Sage, Thousand Oaks, Calif, USA, 2007. View at: Google Scholar
 H. Abdi and D. Valentin, “STATIS,” in Encyclopedia of Measurement and Statistics, N. Salkind., Ed., pp. 284–290, Sage, Thousand Oaks, Calif, USA, 2007. View at: Google Scholar
 H. Abdi, L. J. Williams, D. Valentin, and M. BennaniDosse, “STATIS and DISTATIS: optimum multitable principal component analysis and three way metric multidimensional scaling,” Wiley Interdisciplinary Reviews: Computational Statistics, vol. 4, no. 2, pp. 124–167, 2012. View at: Google Scholar
 S. V. Shinkareva, H. C. Ombao, B. P. Sutton, A. Mohanty, and G. A. Miller, “Classification of functional brain images with a spatiotemporal dissimilarity map,” NeuroImage, vol. 33, no. 1, pp. 63–71, 2006. View at: Publisher Site  Google Scholar
 S. V. Shinkareva, R. A. Mason, V. L. Malave, W. Wang, T. M. Mitchell, and M. A. Just, “Using fMRI brain activation to identify cognitive states associated with perception of tools and dwellings,” PLoS ONE, vol. 3, no. 1, Article ID e1394, 2008. View at: Publisher Site  Google Scholar
 S. V. Shinkareva, V. L. Malave, M. A. Just, and T. M. Mitchell, “Exploring commonalities across participants in the neural representation of objects,” Human Brain Mapping. In press. View at: Google Scholar
 H. Abdi, “Singular value decomposition (SVD) and generalized singular value decomposition (GSVD),” in Encyclopedia of Measurement and Statistics, N. Salkind, Ed., pp. 907–912, Sage, Thousand Oaks, Calif, USA, 2007. View at: Google Scholar
 H. Abdi, “Discriminant correspondence analysis (DCA),” in Encyclopedia of Measurement and Statistics, N. Salkind, Ed., pp. 280–284, Sage, Thousand Oaks, Calif, USA, 2007. View at: Google Scholar
 R. Gittins, Canonical Analysis: A Review with Applications in Ecology, Springer, New York, NY, USA, 1980.
 M. Greenacre, Theory and Applications of Correspondence Analysis, Academic Press, London, UK, 1984.
 P. Celeux and J. Nakache, Analyse Discriminante sur Variables Qualitatives, Polytechnica, Paris, France, 1994.
 A. Leclerc, “Une étude de la relation entre une variable qualitative et un groupe de variables qualitatives,” International Statistical Review, vol. 44, pp. 241–248, 1976. View at: Google Scholar
 G. Perrière, J. R. Lobry, and J. Thioulouse, “Correspondence discriminant analysis: a multivariate method for comparing classes of protein and nucleic acid sequences,” Computer Applications in the Biosciences, vol. 12, no. 6, pp. 519–524, 1996. View at: Google Scholar
 G. Perrière and J. Thioulouse, “Use of correspondence discriminant analysis to predict the subcellular location of bacterial proteins,” Computer Methods and Programs in Biomedicine, vol. 70, no. 2, pp. 99–105, 2003. View at: Publisher Site  Google Scholar
 G. Saporta and N. Niang, “Correspondence analysis and classification,” in Multiple Correspondence Analysis and Related Methods, M. Greenacre and J. Blasius, Eds., pp. 371–392, Chapman & Hall/CRC, Boca Raton, Fla, USA, 2006. View at: Google Scholar
 H. Abdi, “Distance,” in Encyclopedia of Measurement and Statistics, N. Salkind, Ed., pp. 270–275, Sage, Thousand Oaks, Calif, USA, 2007. View at: Google Scholar
 D. Domenges and M. Volle, “Analyse factorielle sphérique: une exploration,” Annales de l’INSEE, vol. 35, pp. 3–83, 1979. View at: Google Scholar
 B. Escofier, “Analyse factorielle et distances répondant au principe d’équivalence distributionnelle,” Revue de Statistiques Appliquées, vol. 26, pp. 29–37, 1978. View at: Google Scholar
 C. Rao, “Use of hellinger distance in graphical displays,” in Multivariate Statistics and Matrices in Statistics, E.M. Tiit, T. Kollo, and H. Niemi, Eds., pp. 143–161, Brill Academic Publisher, Leiden, The Netherlands, 1995. View at: Google Scholar
 H. Abdi and L. J. Williams, “Principal component analysis,” Wiley Interdisciplinary Reviews: Computational Statistics, vol. 2, pp. 433–459, 2010. View at: Google Scholar
 Y. Takane, “Relationships among various kinds of eigenvalue and singular value decompositions,” in New Developments in Psychometrics, H. Yanai, A. Okada, K. Shigemasu, Y. Kano, and Meulman, Eds., pp. 45–56, Springer, Tokyo, Japan, 2002. View at: Google Scholar
 H. Abdi and L. J. Williams, “Jackknife,” in Encyclopedia of Research Design, N. Salkind and B. Frey, Eds., Sage, Thousand Oaks, Calif, USA, 2010. View at: Google Scholar
 M. R. Chernick, Bootstrap Methods: A Guide for Practitioners and Researchers, Wiley, New York, NY, USA, 2008.
 D. Freedman, R. Pisani, and R. Purves, Statistics, Norton, New York, NY, USA, 2007.
 H. Abdi, J. P. Dunlop, and L. J. Williams, “How to compute reliability estimates and display confidence and tolerance intervals for pattern classifiers using the Bootstrap and 3way multidimensional scaling (DISTATIS),” NeuroImage, vol. 45, no. 1, pp. 89–95, 2009. View at: Publisher Site  Google Scholar
 Y. Escoufier, “L’analyse conjointe de plusieurs matrices de données,” in Biométrie et Temps, M. Jolivet, Ed., pp. 59–76, Société Francaise de Biométrie, Paris, France, 1980. View at: Google Scholar
 Y. Escoufier, “Operator related to a data matrix: a survey,” in Proceedings of the 17th Symposium on Computational Statistics (COMPSTAT '07), pp. 285–297, Physica, New York, NY, USA, 2007. View at: Google Scholar
 P. Robert and Y. Escoufier, “A unifying tool for linear multivariate statistical methods: the ${r}_{\upsilon}$coefficient,” Applied Statistics, vol. 25, pp. 257–265, 1979. View at: Google Scholar
 H. Abdi, “Congruence: congruence coefficient, RV coefficient, and Mantel coefficient,” in Encyclopedia of Research Design, N. Salkind and B. Frey, Eds., Sage, Thousand Oaks, Calif, USA, 2010. View at: Google Scholar
 M. Mur and N. Kriegeskorte, “Tutorial on pattern classification,” in Visual Population Codes: Towards a Common Multivariate Framework for Cell Recording and Functional Imaging, N. Kriegeskorte and G. Kreiman, Eds., pp. 539–565, MIT Press, Cambridge, Mass, USA, 2012. View at: Google Scholar
 S. LaConte, S. Strother, V. Cherkassky, J. Anderson, and X. Hu, “Support vector machines for temporal classification of block design fMRI data,” NeuroImage, vol. 26, no. 2, pp. 317–329, 2005. View at: Publisher Site  Google Scholar
 S. C. Strother, J. Anderson, L. K. Hansen et al., “The quantitative evaluation of functional neuroimaging experiments: the NPAIRS data analysis framework,” NeuroImage, vol. 15, no. 4, pp. 747–771, 2002. View at: Publisher Site  Google Scholar
 A. VallejoArboleda, J. L. VicenteVillardón, and M. P. GalindoVillardón, “Canonical STATIS: biplot analysis of multitable group structured data based on STATISACT methodology,” Computational Statistics and Data Analysis, vol. 51, no. 9, pp. 4193–4205, 2007. View at: Publisher Site  Google Scholar
 S. Dudoit, J. Fridlyand, and T. P. Speed, “Comparison of discrimination methods for the classification of tumors using gene expression data,” Journal of the American Statistical Association, vol. 97, no. 457, pp. 77–86, 2002. View at: Publisher Site  Google Scholar
 A. R. McIntosh, F. L. Bookstein, J. V. Haxby, and C. L. Grady, “Spatial pattern analysis of functional brain images using partial least squares,” NeuroImage, vol. 3, no. 3 I, pp. 143–157, 1996. View at: Publisher Site  Google Scholar
 A. R. McIntosh and N. J. Lobaugh, “Partial least squares analysis of neuroimaging data: applications and advances,” NeuroImage, vol. 23, no. 1, pp. S250–S263, 2004. View at: Publisher Site  Google Scholar
 A. Krishnan, L. J. Williams, A. R. McIntosh, and H. Abdi, “Partial Least Squares (PLS) methods for neuroimaging: a tutorial and review,” NeuroImage, vol. 56, pp. 455–475, 2011. View at: Publisher Site  Google Scholar
 A. Connolly, I. Gobbini, and J. Haxby, “Three virtues of similiaritybased multi voxel pattern analysis,” in Visual Population Codes: Towards a Common Multivariate Framework for Cell Recording and Functional Imaging, N. Kriegeskorte and G. Kreiman, Eds., pp. 335–355, MIT Press, Cambridge, Mass, USA, 2012. View at: Google Scholar
 A. J. O'Toole, F. Jiang, H. Abdi, and J. V. Haxby, “Partially distributed representations of objects and faces in ventral temporal cortex,” Journal of Cognitive Neuroscience, vol. 17, no. 4, pp. 580–590, 2005. View at: Publisher Site  Google Scholar
 H. Abdi, “Multidimensional scaling,” in Encyclopedia of Measurement and Statistics, N. Salkind, Ed., pp. 598–605, Sage, Thousand Oaks, Calif, USA, 2007. View at: Google Scholar
 A. Connolly, J. S. Guntupalli, G. Gors et al., “Representation of biological classes in the human brain,” Journal of Neuroscience, vol. 32, no. 8, pp. 2606–2618, 2012. View at: Google Scholar
 M. StLaurent, H. Abdi, H. Burianová, and C. L. Grady, “Influence of aging on the neural correlates of autobiographical, episodic, and semantic memory retrieval,” Journal of Cognitive Neuroscience, vol. 23, no. 12, pp. 4150–4163, 2011. View at: Publisher Site  Google Scholar
 H. Abdi, “Partial least square regression, projection on latent structure regression, PLSregression,” Wiley Interdisciplinary Reviews: Computational Statistics, vol. 2, pp. 97–106, 2010. View at: Google Scholar
Copyright
Copyright © 2012 Hervé Abdi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.