Abstract
Studies of functional MRI data are increasingly concerned with the estimation of differences in spatiotemporal networks across groups of subjects or experimental conditions. Unsupervised clustering and independent component analysis (ICA) have been used to identify such spatiotemporal networks. While these approaches have been useful for estimating these networks at the subjectlevel, comparisons over groups or experimental conditions require further methodological development. In this paper, we tackle this problem by showing how selforganizing maps (SOMs) can be compared within a Frechean inferential framework. Here, we summarize the mean SOM in each group as a Frechet mean with respect to a metric on the space of SOMs. The advantage of this approach is twofold. Firstly, it allows the visualization of the mean SOM in each experimental condition. Secondly, this Frechean approach permits one to draw inference on group differences, using permutation of the group labels. We consider the use of different distance functions, and introduce two extensions of the classical sum of minimum distance (SMD) between two SOMs, which take into account the spatiotemporal pattern of the fMRI data. The validity of these methods is illustrated on synthetic data. Through these simulations, we show that the three distance functions of interest behave as expected, in the sense that the ones capturing temporal, spatial and spatiotemporal aspects of the SOMs are more likely to reach significance under simulated scenarios characterized by temporal, spatial and spatiotemporal differences, respectively. In addition, a reanalysis of a classical experiment on visuallytriggered emotions demonstrates the usefulness of this methodology. In this study, the multivariate functional patterns typical of the subjects exposed to pleasant and unpleasant stimuli are found to be more similar than the ones of the subjects exposed to emotionally neutral stimuli. In this reanalysis, the grouplevel SOM output units with the smallest sample Jaccard indices were compared with standard GLM groupspecific score maps, and provided considerable levels of agreement. Taken together, these results indicate that our proposed methods can cast new light on existing data by adopting a global analytical perspective on functional MRI paradigms.
18 \MakePerPage[2]footnote \cbcolorblue \VerbatimFootnotes
Running Head: GROUP ANALYSIS OF SELFORGANIZING MAPS
Group Analysis of Selforganizing Maps based on Functional MRI using Restricted Frechet Means
Arnaud P. Fournel, Emanuelle Reynaud, Michael J. Brammer,
Andrew Simmons, and Cedric E. Ginestet.
Department of Neuroimaging, Institute of Psychiatry, King’s College London, UK, Laboratoire d’Etude des Mécanismes Cognitifs (EMC), EA 3082, Université Lumière Lyon II, France, National Institute of Health Research (NIHR) Biomedical Research Centre for Mental Health.
Acknowledgments
This work was supported by a fellowship and core funds from the UK National Institute for Health Research (NIHR) Biomedical Research Centre for Mental Health (BRCMH) at the South London and Maudsley NHS Foundation Trust and King’s College London. This work has also been funded by the Guy’s and St Thomas’ Charitable Foundation as well as the South London and Maudsley Trustees. APF has received financial support from the Region RhôneAlpes and the Université Lumière Lyon 2 through an ExploraâDoc grant. We would also like to thank three anonymous reviewers for their valuable inputs.
Correspondence
Correspondence concerning this article should be sent to Cedric Ginestet at the Centre for Neuroimaging Sciences, NIHR Biomedical Research Centre, Institute of Psychiatry, Box P089, King’s College London, De Crespigny Park, London, SE5 8AF, UK. Email may be sent to cedric.ginestet@kcl.ac.uk
KEYWORDS: Barycentre, Frechet Mean, fMRI, Group Comparison, Karcher mean, Multivariate analysis, SelfOrganizing Maps, Unsupervised Learning.
Introduction
Selforganizing Maps (SOMs) were originally introduced by Teuvo Kohonen et al. (2000). A SOM is an unsupervised artificial neural network that describes a training data set as a (typically planar) layer of neurons or output units. Each neuron learns to become a prototype for a number of input units, until convergence of the algorithm. The resulting SOM therefore represents a projection of the inputs into a twodimensional grid. In this sense, SOMs can be regarded as a dimensionreduction clustering algorithm. One of the main advantages of this unsupervised method is that the relative position of the neurons on the grid can be directly interpreted, in the sense that proximity of two units on the map indicates similarity of the prototypal profiles of these units.
SOMs have proved to be useful for datadriven analysis and have become popular tools in the machine learning community (Tarca et al., 2007). In neuroimaging, these methods have been successfully applied to the detection of fMRI response patterns related to different cognitive tasks (Liao et al., 2008, Ngan et al., 2002, Ngan and Hu, 1999, Wismüller et al., 2004). Since SOMs are nonparametric unsupervised neural networks, they do not require the specification of temporal signal profiles, such as haemodynamic response function or anatomical regions of interest in order to generate meaningful summaries of spatiotemporal patterns of brain activity. As a result, these methods have also been used to identify variations in lowfrequency functional connectivity (Peltier et al., 2003). Statistically, however, observe that the absence of a probabilistic model can also be a limitation, as this does not allow for a formal evaluation of the goodnessoffit of the method.
When used for clustering, the SOMs have the following main advantages. Firstly, starting with a sufficient number of neurons, the SOM procedure is able to identify features in the data even when these features are only typical of a small number of input vectors. Secondly, the resulting layer of neurons is arranged according to the similarity of these prototypes in the original data space. This ‘topologypreserving’ property is generally not available in other datareduction techniques, such as independent component analysis (ICA) or means clustering. This is one desirable property that SOMs share with multidimensional scaling. This topological structure facilitates the merging of nodes in order to form ‘superclusters’, which provide a way to visualize and compare highdimensional fMRI data sets. Fischer and Hennig (1999), for instance, have demonstrated the specific relevance of these advantages to the analysis of experimental fMRI data.
One of the outstanding questions in the application of SOMs to fMRI data is whether one can summarize several subjectspecific SOMs into a ‘mean map’, which would pool information over several subjects. In addition, it may be of interest to draw inference over group differences, by comparing the mean maps of several groups of subjects. Here, the term group is used interchangeably with the concept of experimental condition. Hence, two distinct groups need not be composed of different subjects, but may only represent different sets of measurements on the same individuals. One may, for instance, be interested in extracting the SOM that summarizes functional brain activity during a particular cognitive task; or in comparing the restingstate SOM signature of schizophrenic patients with that of normal subjects. Few studies, however, have tackled the problem of formally comparing two or more families of SOMs. Although several authors have proposed distance functions on spaces of SOMs (Kaski and Lagus, 1996, Deng, 2007, Kirt et al., 2007), to the best of our knowledge, none of these researchers have attempted to draw statistical inference on the basis of such comparisons. This lack of methodology highlights a pressing need for developing new strategies that would permit the extension of such multivariate methods from singlesubject analysis to multiple group comparison.
SOM group analysis can naturally be articulated within a Frechean statistical framework. In 1948, Frechet introduced the concept of Frechet mean, sometimes referred to as barycentre or Karcher mean in the context of Euclidean and Riemannian geometry, respectively (Karcher, 1977). The Frechet mean extends this concept to any metric space. This quantity is a generalization of the traditional arithmetic mean, applied to abstractvalued random variables, defined over a metric space. The definition of a generalized notion of the arithmetic mean therefore solely relies on the specification of a metric on the data space of interest. Once such a metric has been specified, the Frechet mean is simply the element that minimizes a convex combination of the squared distances from all the elements in the space of interest. Hence, we can construct a metric space of SOMs by choosing a metric on that space, which permits the comparison of any two given SOMs in that space. Note that, in that context, the chosen pairwise distance function should be a proper metric in the sense that it should satisfy the four metric axioms: (i) nonnegativity, (ii) coincidence, (iii) symmetry and (iv) the triangle inequality (see Searcóid, 2007, for an introduction to metric spaces). In the sequel, we consider the use of different distance functions on spaces of SOMs, which do not satisfy the triangle inequality. Nonetheless, we will show that such distance functions can easily be transformed into proper metrics, using a straightforward manipulation (Mannila and Eiter, 1997). See appendix A. The concept of the Frechet mean has proved to be useful in several domains of applications, including image analysis (Thorstensen et al., 2009, Bigot and Charlier, 2011), statistical shape analysis (Dryden and Mardia, 1998), and in the study of phylogenetic trees (Balding et al., 2009).
In this paper, our purpose is to use the concept of the Frechet mean for drawing statistical inference over several families of subjectspecific SOMs. We thus construct Frechean independent and pairedsample statistics, by analogy with the classical treatment of realvalued random variables. Statistical inference for these different tests are then drawn using permutation of the group labels. In the paper at hand, these statistics will be constructed using the restricted Frechet mean, which has been shown to have desirable asymptotic properties (SverdrupThygeson, 1981), but is more convenient to use from a computational perspective. The restricted Frechet mean is defined as the element in the sample space, which minimizes the squared distances from all the elements in the sample. This formal approach to group inference on families of subjectspecific SOMs has the advantage of allowing a direct representation of the mean SOM in each group, thereby pooling together subjectspecific information. In addition, the proposed methods also allow to formally draw inference at the grouplevel in terms of the chosen distance function.
The paper is organized as follows. In the next section, we give a general introduction to SOMs, and how they are computed highlighting the specific algorithm, which will be used throughout the rest of the paper. In a third section, we describe our proposed Frechean framework for drawing inference on several groups of SOMs. This strategy is entirely reliant on the choice of metric for comparing two given SOMs, and we therefore dedicate a fourth section to the description of several distance functions on spaces of SOMs, which appear particularly wellsuited for the analysis of fMRI data. These methods are tested on synthetic data, under a range of different conditions in section five, and on a real data set in section six. We close the paper by discussing the potential usefulness of this statistical strategy with an emphasis on the critical importance of the choice of the distance function.
SelfOrganizing Maps (SOMs)
We assume here that an fMRI data set is available, which consists of several spatiotemporal volumes , with , for subjects in the experimental group. Each is a matrix, with voxels and time points. In the sequel, it will be of interest to compare several families of such volumes, such that , for experimental conditions. When describing the SOM inference algorithms, however, we will focus on a single subjectspecific data set, X.
Sequential algorithm
A SOM, denoted M, consists of output units or neurons arranged in a twodimensional rectangular grid of size where . For convenience, we here assume that the grids of interest are square grids, such that . Thus, the units of a SOM will be indexed by , where ‘reads’ the units in SOM from left to right and top to bottom. Each entry in M is hence denoted by , and corresponds to the coordinates of that unit in M. That is, is a twodimensional vector representing the position of in M, such that, for instance, , and , and so forth. Each output unit has an associated weight vector , which is, in our case, a time series over data points.
The sequential SOM algorithm takes a set of input units, ’s, corresponding to the rows of the input data X. The steps of the procedure will be indexed by , which denote the iterations of the algorithm, and is the final step at which a stopping condition is satisfied. In our case, the stopping rule is simply the number of iterations, but more sophisticated convergencebased criteria can be used. We firstly initialize the output units in M as random draws from a uniform distribution on . Secondly, an input vector, denoted , is randomly chosen amongst the time series. All voxels are selected at each step of the algorithm, and these input vectors are therefore dependent on . We will thus denote this dependence on the iterations by . For each input vector presented to M, we identify the unit in M, which is the ‘closest’ to the input . Here, closeness is generally measured in terms of Euclidean distance with respect to the values taken by the input vectors. The unit in M, which is the closest to is referred to as the Best Matching Unit (BMU). The index of that BMU, for a given input vector at iteration , is defined as follows,
(1) 
with denoting the Euclidean norm on . Here, and are dimensional time series. Thirdly, we update the BMU and its neighbors. The new values of these units are defined as a linear relationship of the input vector . For a given , the updating rule for the BMU and its neighbors is the following,
(2) 
for every . After updating their weights, the BMU and its neighbours are closer to in the sense that they constitute a better representation of that input vector. These steps are repeated for a fixed number of iterations, . The updating rule in equation (2) contains two key parameters: (i) the learning rate, denoted and (ii) the kernel function represented by , which grows smaller as we consider units in M, which are further away from the BMU in the space of coordinates of M. We describe these two quantities in turn.
The learning rate, in equation (2), is a decreasing function of the number of iterations, , which controls the amount of learning accomplished by the algorithm –that is, the dependence of the values of the units in M on the inputs. By convention, we have for every . Three common choices for are a linear function, a function inversely proportional to the number of iterations and a power function, such as the following recursive definition,
(3) 
for every . A popular initialization for the learning rate is (Peltier et al., 2003, González and Dasgupta, 2003). Clearly, as the algorithm progresses towards , the value of decreases towards . Note, however, that, in the paper at hand, we use the batch version of this algorithm, which does not require the specification of a learning rate.
In equation (2), we have also made use of the neighborhood kernel, . As for the learning rate, the value taken by this kernel decreases with the number of iterations, and is therefore dependent on . This dependence on has been emphasized through a subscript on . For a given output unit in the map M, the neighborhood kernel, , quantifies how ‘close’ is to the BMU, which has index . Observe that this closeness is expressed in terms of Euclidean distances on the grid coordinates. As commonly done in this field, we here choose a standard Gaussian kernel to formalize the dependence of each unit on the values of its neighbors, such that
(4) 
where represents the twodimensional dot product. Here, is a linear function of the number of iterations, which controls the size of the neighborhood around the BMU. This function is defined recursively as , where is a parameter value that represents the initial neighborhood radius. This parameter is commonly initialized with respect to the size of the twodimensional grid, M, such that , which is the ‘height’ of the output SOM.
Batch algorithm
A popular alternative to the sequential SOM algorithm described in the previous section is the batch SOM algorithm, which has the advantage of being more computationally efficient than its sequential counterpart (Vesanto and Alhoniemi, 2000). It has been successfully used in the context of fMRI analysis (Ngan et al., 2002), in natural language processing (Kohonen et al., 2000), and in the face recognition literature (Tan et al., 2005). The main difference between these two approaches is that the entire training set is considered at once in the batch SOM algorithm, which permits the updating of the target SOM with the net effect of all the inputs.
This ‘global’ updating is performed by replacing the input vector, denoted in the previous section, with a weighted average of the input vectors, where the relative weight of each input vector is proportional to the neighborhood kernel values. At the step of the algorithm, we are therefore conducting the following global updating,
(5) 
for every . It can easily be seen from equation (5) that is a convex linear combination of the input vectors, ’s, where each of the inputs is weighted by , and the sum of these weights is equal to 1. Another nonnegligible advantage of the batch SOM algorithm is that it removes the dependence of the outputs on the learning rate parameter, denoted , as stated in the previous section.
Throughout the rest of the paper, we will make use of the batch algorithm, with , and , thereby producing SOMs of dimensions . Output units in all SOMs are initialized randomly. Other groups of researchers in neuroimaging have used square SOMs (Peltier et al., 2003, Liao et al., 2008). However, we have also investigated using simulated data, whether the specification of rectangular maps had a significant impact on our proposed inferential methods (see appendix B).
SOM Group Frechean Inference
The question of inferring the statistical significance of the difference between two families of SOMs can be addressed through the use of abstractvalued random variables as advocated by Fréchet (1948). In this approach, random variables are solely defined with respect to a probability measure on a metric space, (see Parthasarathy, 1967, chap. 2). Hence, it suffices to define a metric on the space of interest, in order to obtain a valid statistical framework. Once such a metric has been chosen, one can construct the mean element in that space, which is commonly referred to as the Frechet mean.
In the paper at hand, we are considering a space of SOMs, which we may denote by , where is a metric on that space. A range of different distance functions for such spaces of SOMs will be described in the next section. As in most standard fMRI designs, we assume that we have experimental conditions, with subjects in each condition, thereby allowing for a different number of subjects in each experimental condition. A full data set will be summarized as an array of SOMs, , with and , such that corresponds to the SOM of the subject in the condition. Given such a sample of SOMs, we can then define the Frechet mean for the condition as follows,
(6) 
where we have used the Bessel’s correction (i.e. ) by analogy with the realvalued setting. Given the complexity of the underlying space of SOMs, such a minimization may be unwieldy. As a result, it is computationally more practical to consider the restricted Frechet mean, as introduced by SverdrupThygeson (1981). The classical Frechet mean in equation (6) is obtained by identifying the element in the population of SOMs, which minimizes the average squared distances from all the elements in the sample. The restricted Frechet mean, by contrast, is obtained by identifying the element in the sample, which has this property. Hence, the restricted Frechet mean is computed as follows,
(7) 
where denotes the sampled SOMs in the condition, such that . The restricted Frechet mean has been shown to be consistent, through a generalization of the strong law of large numbers due to SverdrupThygeson (1981). Asymptotically, converges almost surely to a subset of the theoretical restricted mean, which takes values in the support of the target population distribution. In the sequel, the theoretical restricted Frechet mean for the condition will be denoted by , following standard convention.
Similarly, one can define the conditionspecific sample Frechet variances. These quantities are simply the values taken by the criteria, which are minimized in equation (7), such that the (restricted) Frechet variance for the condition is defined with respect to the restricted Frechet mean in the following manner, for every ,
(8) 
Using the restricted Frechet mean and variance, it is now possible to construct a nonparametric test on the metric space of SOMs. Here, we therefore assume that we solely have two experimental conditions, such that . The null hypothesis stating that the (restricted) Frechet means of these two distributions are separated, can be formally expressed as follows, . Naturally, our interest will especially lie in testing the null hypothesis stating that there is no difference between the theoretical restricted Frechet means, which corresponds to . This can be tested using the following Frechet statistic,
(9) 
where the denominator, , is the classical pooled sample variance, which is defined by analogy with the realvalued setting as
In addition, if one is considering two samples of equal sizes and assuming equal Frechet variances, then the aforementioned statistic for such a mean difference can be defined as follows
(10) 
where, in this case, the pooled variance is simply the sum of the variances of the two samples, such that , and with . Statistical inference is then conducted using permutation on the group labels. Although our proposed statistic is a realvalued random variable, its asymptotic distribution is unknown. Indeed, the behavior of this statistic depends on a large number of other random variables, which are combined using the nonlinear procedure for obtaining grouplevel SOMs. As a result, the permutationbased distribution of under the null hypothesis is not expected to follow a standard distribution. In particular, the null distribution of need not be symmetric. Since we are here solely considering a generalization of the test, but will be applying this statistic to more than two experimental conditions, we will also make use of the standard Bonferroni correction for multiple testing.
Choice of Distance Functions
In our proposed approach to group comparison, the choice of the metric on the space of SOMs is paramount. Different distance functions capture different aspects of the SOMs under scrutiny. It is therefore of interest to evaluate group differences with respect to several choices of distance functions. We here review the main distance functions, which have been previously proposed in the literature for comparing two given SOMs. In addition, we introduce a spatiotemporal sum of minimum distances, which is especially relevant for the study of fMRIbased SOMs.
Quantization Error and Other Measures
This measure is not a metric, but a popular tool for evaluating the accuracy of the SOM generated from a given data set. The socalled quantization error measures the average quantization error of the target SOM (Kohonen, 2001). It is defined as the sum of the Euclidean distances between each input unit, , and its best matching prototype on M –that is, the BMU of . The quantization error, denoted , is thus formally defined as follows,
where, as before, denotes the dimensional Euclidean norm and is the index of the BMU in M with respect to as described in equation (1). This measure is a good indicator of the convergence of a SOM, and is often used when assessing the behavior of the algorithms described in the previous section. In this paper, we will use a variant of the quantization error in order to identify the output units, which explain the largest amount of betweensubject ‘variance’ in the data. However, the quantization error does not allow the computation of the distance between two given SOMs.
Kaski and Lagus (1996) have proposed a measure of dissimilarity between two SOMs. They proceeded by comparing the shortest path on each SOM after matching a given pair of input vectors. This dissimilarity measure is computed by comparing the distances between all pairs of data samples on the feature maps. This method, however, is not computationally efficient, and would be especially challenging when considering fMRI data sets, where neuroscientists are commonly handling about 100,000 input vectors –that is, the voxelspecific time series– for every subject.
Sum of Minimum Distances (TSMD)
The Sum of Minimum Distances (SMD) was originally introduced by Mannila and Eiter (1997) and has been widely used in image recognition and retrieval (Kriegel, 2004, Takala et al., 2005, Tungaraza et al., 2009). Moreover, the SMD function and some of its variants have already been used in order to tackle the problem of comparing several SOMs (Deng, 2007). Given two SOMs, denoted and for input data sets X and Y, respectively, the SMD can be computed as follows. For every unit, in , we calculate the Euclidean distance between and every unit in in order to retain the unit in that minimizes this distance. These minimal distances are summed and then normalized by the total number of input vectors, denoted , in our case. This gives an to score. The same procedure is performed in the opposite direction in order to produce an to score. The average of the to and the to scores is then defined as the overall SMD between and . Therefore, this distance function compares SOMs on the basis of the dissimilarity of the time series underlying each output unit. It follows that this procedure mainly emphasizes temporal differences between the fMRI volumes of interest. Thus, we will label this classical SMD as temporal SMD, and denote it by TSMD. It is formally defined as
(11) 
where is the Euclidean distance between and on , where and represent dimensional prototypal time series for maps and , respectively.
It is important to note that the SMD function can be rewritten by treating a map, , as a set of weight vectors, . In this case, we consider the metric space of all weight vectors, . This metric space is . By a slight abuse of notation, the SOM, , will be used to denote the set of all output vectors, associated with the units in . Therefore, we have . As a result, we can apply the classical definition of the distance between the subset of a metric space and an element of that space, , such that . Using these conventions, it becomes possible to reformulate the SMD function in equation (11) in the following manner as stated by Mannila and Eiter (1997),
In addition, observe that the SMD function is not in general a proper metric, in the sense that the triangle inequality may fail to be satisfied (see Mannila and Eiter, 1997, for a counterexample). However, one can easily produce a proper metric through the identification of the shortest paths between any two elements in the space of interest, and then define a new metric with respect to these shortest paths (see appendix A). It can easily be shown that such a transformation necessarily produces proper metrics, when considering metrics based on the SMD function (Mannila and Eiter, 1997). This particular procedure can easily be implemented in our case, because we have focused our attention on the restricted Frechet mean, where the minimization required to identify the mean element is solely conducted over the space of the sampled elements. As a result, there exists a small number of possible shortest paths between every pair of elements in the sample, which greatly facilitates the required transformation for producing a proper metric. This procedure was systematically conducted in the sequel, and therefore all the variants of the SMD function utilized in this paper are indeed proper metrics. We will thus assume throughout this paper that all distance functions have been adequately transformed. We now introduce two novel variants of the SMD function, which take into account the spatial and spatiotemporal properties of the fMRI data.
Spatial SMD
One may also be interested in quantifying the amount of ‘spatial overlap’ between two given SOMs. This question is especially pertinent when analyzing SOMs based on fMRI data sets. Here, we therefore wish to evaluate whether the units in two different maps correspond to similar subsets of voxels in the original images. Such a distance can be quantified through a slight modification of the aforementioned SMD metric, where the Hamming distance is used in the place of the Euclidean distance.
(12)  
with denoting the number of voxels in the fMRI volumes of interest, and where denotes the binarized index vector of the voxels, which have been assigned to unit in . That is, if the voxel has been assigned to , then , otherwise, we have . In addition, we have here made use of the celebrated Hamming distance, which takes the following form (Hamming, 1950),
(13) 
where stands for the indicator function taking a value of if is true, and otherwise. Here, the term spatial refers to the spatial distribution of the voxels allocated to a particular output unit. Hence, SSMD does not emphasize the spatial location of the output units, as these allocations are arbitrary, but rather the spatial distribution of the voxels allocated to that output unit. Observe that we have here solely considered differences in the spatial distributions of the best matched pair of SOM units, where that matching is done through the minimization reported in equation (12). This approach, however, omits to take into account the similarity of these SOM units as prototypal time series. Both the spatial and the temporal aspects of these maps can nonetheless be combined, as described in our proposed spatiotemporal SMD.
SpatioTemporal SMD
In this novel variant of the classical SMD, we are quantifying the amount of spatial overlap between any pair of output units in two distinct maps. In contrast to the SSMD described in the previous section, however, we are here comparing the spatial distributions (i.e. the sets of voxel indexes assigned to a particular unit) of the units that are the closest in terms of time series profiles, thereby combining the temporal and spatial properties of the data. This spatiotemporal version of the SMD function is defined as follows,
(14) 
where
where, again, is the Euclidean distance on , and is the index set of the voxels in X, whose time series are best represented by . Albeit the formulae in equation (14) is somewhat convoluted, it corresponds to, perhaps, the most intuitive perspective on the problem of comparing SOMs, when using fMRI data. Indeed, we are quantifying the amount of spatial overlap between units, which are as similar as possible in terms of their temporal profiles.
We summarize this section with a concise description of these three different types of SOM distance functions:

Temporal SMD (TSMD) is based on the sum of the minimum Euclidean distances between the time series of the output units.

Spatial SMD (SSMD) is based the sum of the minimum Hamming distances between the sets of voxels allocated to the output units.

Spatiotemporal SMD (STSMD) is based on the sum of Hamming distances between the sets of voxels allocated to the output units, which are the most similar in terms of their time series.
Synthetic Data Simulations
The proposed methods were tested on three different simulated scenarios, with varying degrees of difficulty. In particular, by isolating different types of differences between the two groups of interest, each of these scenarios emphasizes the need for the use of specific metrics, capturing different aspects of the spatiotemporal process under investigation. In this sense, our simulations strive to produce realistic differences between the two families of fMRI volumes, where group differences may be either spatial, temporal or both spatial and temporal.
Simulation Scenarios
For each simulated data set we have constructed two groups of 20 subjects, where each subjectspecific data set is composed of twodimensional images with voxels, over 50 time points, as represented in panels (ac) of figure 1. The time series at each voxel can be of three different types, as illustrated in panels (df) of figure 1, composed of two different signals and one background time series. The two signals represented in panels (d) and (e) are sinusoids oscillating over , with a frequency of Hz and Hz, respectively. We have then added a vector of Gaussian random noise, z, to these two types of time series, such that , for every , and for different choices of . The background noise time series, in panel (f), is solely composed of the random noise for a given .
The three scenarios in panels (ac) of figure 1 are ordered in terms of the degree of ‘separability’ of the two groups, where the easiest scenario is on the left and the most difficult one is on the right. The first data set (SC1) was built with three different time series at three different locations, corresponding to the purple, blue and gray colors. In this scenario, the groups differ both in terms of the temporal profiles of some of their voxels and in terms of the spatial locations of these voxels. The second scenario (SC2) was constructed with two different types of time series. Here, the two groups solely differ in terms of the temporal profiles of some of their voxels. Finally, in the third scenario (SC3), the two groups only differ in terms of the spatial locations of the voxels, which have been assigned the second temporal profile.
In addition, we also studied the effect of the signaltonoise ratio (SNR) on the performance of our inferential methods. In particular, we varied our choice of , when generating the different time series displayed in panels (df) of figure 1, in order to produce different SNRs. In these simulations, the ‘signal’ of interest was defined as the amplitude of the original sinusoids, which oscillated between 1 and 1, thus giving an amplitude of . Since the noise affecting this signal was specified to be Gaussian, the SNR was defined, in our case, as . Thus, by setting to either , or , we produced three different SNRs of , and , respectively.
These synthetic data sets were analyzed using our proposed inferential framework. For each simulated subjectspecific volume, a SOM was computed, and the restricted Frechet mean was identified for each group. In all scenarios, SOMs were produced by using the batch SOM algorithm. The output grid was of size with ; the number of iterations was set to steps; and we used a decreasing neighborhood kernel of size , as commonly done in this field (Kohonen et al., 2000). For computational convenience, statistical inference was drawn in each scenario after 100 permutations. Each simulated scenario was reproduced 100 times, and constructed for the three different levels of SNR, thereby totalling distinct simulations.
Scenarios and Factors  TSMD  SSMD  STSMD  

SC1 (Spatiotemporal)  
SC2 (Temporal)  
SC3 (Spatial)  
Simulation Results
The summary results of the analysis of these synthetic data sets are reported in table 1 and figure 2. Overall, the different metrics of interest were found to successfully capture the aspects of the simulated SOMs that they were expected to identify. That is, in the first column of table 1, one can see that TSMD, which solely takes into account the differences in voxelspecific temporal profiles, attains its most significant values under the temporal scenario, SC2. Similarly, in the second column of table 1, the spatial version of the SMD metric, denoted SSMD exhibits its best performance under the first and third scenarios, denoted SC1 and SC3, respectively. Indeed, these two scenarios are the only ones, where the two groups can be discriminated in terms of the spatial locations of the different types of time series. Finally, in the third column of table 1, the spatiotemporal metric, STSMD, appears to be optimal under the first scenario, where group differences can be characterized through the spatiotemporal properties of the simulated images.
In addition, we have also evaluated the effect of sample size on the capacity of the metrics to detect group differences. These results are not reported in this paper, but we have observed, as expected, that the statistical power of all the studied metrics improved as the number of subjects in each group increases. In particular, it was noted that for the STSMD, we solely needed in each group to identify significant differences under SNR=1, and group sizes of under SNR=0.5; for all scenarios. Exemplary null distributions for statistics in the three different scenarios and the three different levels of SNR with are reported in figure 2.
In sum, one can note that these three distance functions exhibit different levels of robustness. In particular, SSMD appeared to be especially sensitive to noise. Although SSMD succeeded to capture the spatial differences simulated in SC3, it only outperformed STSMD for high SNRs. Also, in the spatiotemporal scenario (SC1), TSMD behaved as well or better than STSMD. This suggests that the TSMD function is more ‘powerful’ than the STSMD, even when the group differences are characterized by both spatial and temporal properties. However, the use of STSMD remains justified, because it also succeeds to capture spatial differences, whereas TSMD fails to do so. In general, we therefore recommend the joint use of TSMD and STSMD: If only TSMD indicates the presence of group differences, then one can conclude that such differences are mainly of a temporal nature; whereas when the use of STSMD indicates greater group differences, this suggests that such differences also have a spatial character. Overall, these simulations highlight the importance of using several types of distance functions, as there may not exist a single type of metric, which would capture all of the aspects of the data of interest.
Experimental Data
We also evaluated our methods with the reanalysis of a classical data set, originally published by MouraoMiranda et al. (2006). Since this first publication, this data set has been reanalyzed several times with different machine learning algorithms, as conducted by MouraoMiranda et al. (2007) using spatiotemporal support vector machine (SVM), and Hardoon et al. (2007) with unsupervised methods.
Subjects and Experimental Design
This data set consists of fMRI data from 16 righthanded males with a mean age of 23 years. All participants had normal eyesight and no history of neurological or psychiatric disorders, and gave written informed consent to participate in the study, in accordance with the local ethics committee of the University of North Carolina (see MouraoMiranda et al., 2006). Data were acquired using an experimental block design, composed of three different conditions: (i) exposure to unpleasant visual stimuli (i.e. photos of dermatological diseases), (ii) exposure to neutral visual stimuli (i.e. photos of neutral daytoday scenes including human actors) and (iii) exposure to malerelevant pleasant visual stimuli (i.e. scantly dressed women or women in swimsuits). The entire experimental design consisted of six blocks, where each block contained seven images, which were each presented to the subjects for three seconds. Each block was followed by a resting block period where subjects were solely exposed to a fixation cross. All blocks were of 21 in length.
Metrics  Tests  values 

TSMD  Pleasant vs. Neutral  0.0 
Unpleasant vs. Neutral  0.0  
Pleasant vs. Unpleasant  0.258  
SSMD  Pleasant vs. Neutral  0.412 
Unpleasant vs. Neutral  0.518  
Pleasant vs. Unpleasant  0.423  
STSMD  Pleasant vs. Neutral  0.021 
Unpleasant vs. Neutral  0.013  
Pleasant vs. Unpleasant  0.128 
Data Acquisition and PreProcessing
Bloodoxygenationleveldependent (BOLD) signal was measured using a 3Tesla Allegra headonly MRI System at the Magnetic Resonance Imaging Research Center in the University of North Carolina. The scanning parameters were specified as follows. Voxel size was , TR was , TE was , FA was , FOV was and each MRI volume had dimensions . In each experiment, a total of functional volumes were acquired for each participant.
Data were preprocessed using the FSL Software suite (Smith et al., 2004a); through the use of the Nipype Python Library (Gorgolewski et al., 2011). All fMRI volumes were first motioncorrected and the skulls were removed, after tissue segmentation. The voxel time series were detrended and filtered in time and space: that is, lowfrequency (drift) fluctuations were reduced using a bandpass temporal filter comprised between 0.008Hz and 0.1Hz. The use of a such a bandpass filter is typical of restingstate connectivity analysis (Cordes et al., 2001). In addition, spatial smoothing was performed using an 8mm fullwidth at halfmaximum Gaussian kernel.
The first two volumes of each block were discarded, to allow for the betweenblock lag in haemodynamic response. The remaining volumes were concatenated to form three distinct time series representing the three different conditions. Time series concatenation in the context of functional connectivity has been introduced by Fair et al. (2007) and has been implemented by several authors for the study of functional MRI networks (Ginestet and Simmons, 2011, Ginestet et al., 2011).
Sample Jaccard Index  
Ranked SOM Units  Neutral  Pleasant  Unpleasant 
Units up to 1  0.68  0.32  0.55 
Units up to 2  0.11  0.28  0.16 
Units up to 3  0.08  0.04  0.03 
Units up to 4  0.01  0.07  0.07 
Units up to 5  0.03  0.08  0.03 
Units up to 6  0.01  0.05  0.04 
Units up to 7  0.05  0.08  0.05 
Units up to 8  0.01  0.02  0.05 
Units up to 9  0.02  0.06  0.02 
Results of Realdata Analysis
The significance results of all the pairwise comparisons are reported in table 2, after applying the Bonferroni correction for multiple testing. Our reanalysis of this data set has highlighted a substantial degree of difference between the neutral and pleasant conditions, on one hand, and between the neutral and unpleasant conditions, on the other hand. These pairwise differences were found to be highly significant under the TSMD distance function. The mean SOM in the unpleasant condition was also found to be significantly different from the mean SOM in the neutral condition with respect to the STSMD metric (), albeit to a lesser extent than under the TSMD function. The difference between the mean SOMs in the pleasant condition and the neutral one was found to be just above significance level under the STSMD metric (). By contrast, none of these differences reached significance under the SSMD, indicating that differences in spatial allocation of the different output units of the SOMs in these experimental conditions were not sufficient to distinguish between the grouplevel SOMs. Importantly, the mean SOMs under the pleasant and unpleasant conditions were not found to be significantly different under all three metrics.
As noted in the analysis of the synthetic data, the fact that the SOMs are computed on the basis of the similarities between the voxelspecific time series is likely to be responsible for the important differences that we are reporting between the metrics capturing the temporal aspects of the process and SSMD, which does not emphasize differences in temporal profiles. Intriguingly, it should also be noted that the differences between the SOMs in the pleasant and unpleasant conditions under STSMD is lower than under the TSMD, thereby perhaps suggesting that the differences between these two conditions is characterized by an ‘interaction’ between the temporal and spatial properties of these functional patterns.
Visualization of Grouplevel SOMs
In each of the three conditions, we have represented the subjectspecific restricted Frechet mean in order to produce robust visual summaries of the different output units in each condition. This was conducted by identifying the output units that explained the largest amount of ‘variance’ in the data, in terms of sample Jaccard index. This measure quantifies how representative is a mean output unit in terms of the overlap of this unit with the output units of all the subjects in that group.
For each experimental condition, we identified the output unit in the subjectspecific SOMs that explained the largest amount of Jaccard index. For each experimental condition, we have plotted in figures 3, 4 and 5, the three output units that are associated with the least amount of Jaccard index over all subjects. That is, for the condition, the sample Jaccard index of the output unit in the restricted Frechet mean for that condition was defined as follows,
Here, denotes the global Jaccard distance between a mean SOM and a subjectspecific SOM. Moreover, the classical vectorspecific Jaccard distance is
where is the number of elements satisfying and , for running from to the length of these vectors; and and are defined similarly. The three output units minimizing the sample Jaccard index in each of the three experimental conditions have been plotted in figures 3, 4 and 5, for the neutral, pleasant, and unpleasant conditions, respectively. Allocation of a voxel to the particular output unit of a SOM is ‘hard’, in the sense, that either a voxel is included into an output unit or it is included in a different one. Thus, in figures 3, 4 and 5, we have provided the spatial location of the voxels that have been assigned to particular output units in each condition.
From these three figures, one can observe that the neutral condition is characterized by a distinct network of brain regions identified as the third output unit in the neutral condition, as can be seen from panel (c) of figure 3. The set of regions associated with this particular output unit may be interpreted as a visual network, since it contains a considerable number of regions located in the occipital lobe. This particular output unit was not found to be present in the three units with the least Jaccard index in either the pleasant or unpleasant condition, as can be noted from figures 4 and 5, respectively.
Comparison with Standard GLM Maps
The grouplevel SOM output units selected using the sample Jaccard index were compared with standard general linear model (GLM) score maps. Separate GLM analyses were conducted for each experimental condition, using the FEAT function provided in the FSL software suite (Smith et al., 2004b). The score maps were thresholded at , for comparison purposes. These binarized score maps were compared with the maps of the three ‘best’ grouplevel output units, selected on the basis of their sample Jaccard indices, as described in the previous section. These thresholded GLM maps have been overlaid in blue over the selected output units in figures 3, 4 and 5.
In each experimental condition, we computed the percentage overlap between the maps obtained using these two different techniques. That is, in each condition, we evaluated how many voxels were present in both the thresholded score maps and each output unit, normalized by the number of voxels included in the score maps. These numerical comparisons are reported in table 3, where we have described the individual percentage overlap of the nine different output units, ranked with respect to their sample Jaccard indices. The combined percentage overlap of the three output units exhibiting the lowest sample Jaccard indices with the thresholded GLM score maps was , , and for the neutral, pleasant and unpleasant conditions, respectively. Although our proposed SOMbased methods summarize such fMRI volumes in a nonlinear manner, these numerical comparisons show that the resulting output units exhibit a considerable degree of agreement with standard GLM analysis.
Discussion
Advantages of Proposed Methods
The main contribution of this paper is the construction of an inferential framework for the comparison of grouplevel SOMs. Although some previous researchers have considered various ways of comparing SOMs (Kaski and Lagus, 1996, Deng, 2007, Kirt et al., 2007), to the best of our knowledge, no authors have yet treated the problem of evaluating the statistical significance of such group differences. In addition, observe that our Frechean approach to statistical inference can be conducted for any choice of metrics. Indeed, the idea of defining a group distance statistic, such as the Frechean test in equation (9) and then evaluating its significance using permutation can be implemented for any choice of distance functions. Therefore, this allows the specification of a rich array of different distance functions capturing different aspects of the SOMs under scrutiny. As illustrated in the main body of the paper, we have shown that classical distance functions such as the SMD can be modified in order to emphasize spatial, temporal or spatiotemporal differences between the groups of interest. Here, the choice of hypothesis is therefore superseded by the choice of metrics over the space of SOMs. In particular, this inferential framework allows to test hitherto untestable grouplevel hypotheses.
Another substantial advantage of combining a Frechean approach with the computation of subjectspecific SOMs is that this bypasses the problem of multiple testing correction. In standard massunivariate analyses of MRI volumes, one needs to control for the inflation of the number of false positives introduced by performing a large number of nonindependent statistical tests. By contrast, we are here conducting a single test, which identifies whether the volumes of interest are different at a multivariate level, through the comparison of two nonparametric unsupervised representations of the original data.
Finally, one should also note that the use of the restricted Frechet mean in our proposed framework is advantageous for several reasons. On the one hand, the restricted Frechet mean greatly reduces the computational cost of our overall analytical procedure. This is especially true, because inference was drawn using permutations of the group labels, and it is not clear whether such a large number of permutations would have been possible, lest for the use of the restricted Frechet mean. On the other hand, the restricted Frechet mean has also the advantage of quasiautomatically transforming any distance function into a proper metric that satisfies the four metric axioms. This results in a nonnegligible simplification of the probabilistic theory needed to justify our inferential approach. Indeed, most of the asymptotic results, which have been previously established relies on the postulate that the distance function of interest is a proper metric (Ziezold, 1977, SverdrupThygeson, 1981).
Limitations of the Frechean Framework
One can identify three substantial limitations to our proposed Frechean inferential framework for SOMs, which are (i) a lack of contrast maps, (ii) a reliance on permutation for statistical inference, and (iii) the use of the restricted Frechet mean in the place of the unrestricted mean element in the space of all SOMs. We address these three limitations in turn.
Firstly, one of the important limitations of our current method is that it does not directly permit the production of a ‘groupdifference SOM’ representing the difference between two groupspecific Frechet means. In particular, this implies that we cannot represent such differences by plotting a differential pattern of activation or connectivity, as is commonly done using standard massunivariate approaches. See the statistical parametric network (SPN) approach, advocated by Ginestet and Simmons (2011) for example, when conducting functional network analyses. From a neuroscientific perspective, this is a considerable limitation, as it diminishes the interpretability of the results. We will consider different ways of tackling this issue and producing groupdifference SOMs in future work.
Secondly, our inferential framework has relied on permutation for evaluating the statistical significance of the test statistics under scrutiny. This was made computationally feasible, because this choice of inferential method was used in conjunction with the restricted Frechet mean. That is, for each permutation of the group labels, the computation of the groupspecific Frechet means was straightforward, because the identification of the restricted Frechet mean can be conducted by using the margins of the dissimilarity matrix of the sample points –that is, the dissimilarity matrix of all the subjectspecific SOMs. Hence, the cost of calculating the group means at each permutation was small, and the full inference could be drawn within a couple of hours on a standard desktop computer. Such a level of computational efficiency may not have been achieved if we had attempted to derive the unrestricted Frechet mean, as described in equation (6), which would have necessitated to perform a minimization over the space of all possible SOMs.
However, although the use of the restricted Frechet mean was advantageous from a computational perspective, this particular methodological choice has also its limitations. Indeed, the use of the restricted Frechet mean in the place of the unrestricted mean results in a loss of the classical benefits usually associated with computing an average of real numbers. In decisiontheoretic parlance and when considering realvalued random variables, the arithmetic mean is the quantity that minimizes the squared error loss (SEL) (see Berger, 1980, for an introduction to decision theory). The restricted version of the arithmetic mean for realvalued random variables would also minimize the restricted SEL. However, the restricted arithmetic mean would necessarily achieve a sample variance greater or equal to the one of the unrestricted frechet mean. Note, however, that the problems associated with the utilization of the restricted Frechet mean are also mitigated by the fact that computing this quantity quasiautomatically makes the space of interest a metric space, regardless of the particular choice of distance function.
We have here used the sample Jaccard index for selecting the most ‘relevant’ output units in the grouplevel SOMs. It certainly does not follow from such a selection criterion that these output units are of greater physiological relevance. This criterion is entirely statistical, and the resulting interpretation of these output units should remain statistical. In practice, it is advisable to visualize the entire set of output units obtained after this type of SOM analysis, in order to identify relevant physiological differences based on prior neuroanatomical knowledge.
Possible Extensions of these Methods
Our proposed Frechean inferential framework could be extended in a range of different directions. One of the most natural extensions of this method would be to devise an test, which would generalize the aforementioned twosample statistic. A Frechet statistic may take the following form. Let a data set of the form , where labels the objects in the group with . By analogy with the classical realvalued setting, the statistic can be defined as the ratio of the betweengroup to withingroup variances, , where these quantities are here defined with respect to the Frechet moments, such that , and , using standard notation for the Frechet sample group means, , and grand mean, . One can then test for the null hypothesis that , where and are the theoretical equivalents of and , respectively. Statistical inference can, again, be conducted using permutation of the group labels.
In addition, the analytical strategy that we have here described could also be improved through the use of different types of SOM algorithm. In the present paper, we have made use of the batch SOM algorithm. However, several other alternatives to the traditional sequential SOM algorithm have been proposed in the literature. In particular, Vesanto and Alhoniemi (2000) have showed that the SOMs obtained when using the batch SOM algorithm with an initialization of the maps based on the eigenvectors of the input data can produce more robust results. Since every SOM is computed independently for each subject, such an improvement of the existing batch SOM algorithm could easily be incorporated in our proposed inferential framework.
One of the outstanding questions that is implicitly raised in this paper is the possibility of separately weighting the individual contributions of the spatial and temporal properties contributing to the overall SOM difference. Such a question is likely to be arduous to answer, however, since the temporal and spatial properties of the fMRI volumes of interest necessarily live in distinct abstract spaces. On one hand, the temporal differences in TSMD were quantified using a Euclidean distance in a dimensional vector space; whereas, on the other hand, the spatial differences in SSMD were quantified using the Hamming distance on binary vectors of varying sizes. It is unclear whether the magnitude of the distances in these different metric spaces could be normalized in order to ensure a modicum of comparability.
Conclusions
In this paper, we have described a formal framework for drawing grouplevel inference between unsupervised multivariate summaries of fMRI data. Our proposed approach proceeds by computing subjectspecific SOMs, and computing the sample Frechet mean in each group of subjects. Despite the unwieldy nature of the space of all possible SOMs, this can be done efficiently by identifying the restricted Frechet mean. Statistical inference on the difference between the group restricted Frechet means can be conducted using permutation on the group labels. This framework can be implemented for any choice of metrics quantifying the difference between pairs of SOMs. As such, the specification of a particular distance function is equivalent to the choice of a particular hypothesis test. Different researchers may therefore be interested in evaluating different metrics, which capture different aspects of the SOMs.
We have hence described and evaluated several types of distance functions for SOMs based on fMRI data. In particular, we have considered variants of the classic SMD function, which has previously been used to compare pairs of SOMs. Our proposed variants distinguish between the temporal, spatial and spatiotemporal properties of the data under scrutiny. Our inferential framework and these metrics were tested on both synthetic and real data, Our analysis of the simulated data showed that the distance functions of interest were indeed capturing the aspects of the data that they were purported to measure. In addition, the findings of the reanalysis of an fMRI experiment has demonstrated the capacity of these methods to extract new information from existing data sets. In this paradigm, the differences of the restricted mean SOMs in the pleasant and unpleasant conditions were found to be smaller than the differences between the mean SOMs in any of these two conditions with respect to the one in the neutral condition.
Taken together, the analyses of these synthetic and real data sets have underlined the robustness and potential usefulness of these methods. It is hoped that this type of global inferential perspective on neuroimaging data will inspire other neuroscientists to follow this research avenue. One could imagine a range of other subjectspecific abstractvalued random variables that could be suitably analyzed using this type of Frechet inferential framework. In fact, the very use of massunivariate approaches in the context of neuroimaging could be superseded by a more global perspective, where a single statistical test is conducted, thereby bypassing the need for exacting multiple testing penalties.
Appendices
A. From Distance Functions to Metrics
Let be a distance function on a finite space of SOMs, , which satisfies the positivity, coincidence and symmetry axioms. In order to transform the distance function into a proper metric satisfying the triangle inequality, we need to construct a saturated graph representing the topology of . The vertex set of is defined as . Its edge set is composed of all the possible links between the elements of . That is, is a saturated graph, in the sense that it contains the maximal number of edges. Each of these edges is denoted by , for any .
A path in is a nonempty subgraph of the form and , where the ’s are all distinct. Following an idea proposed by Mannila and Eiter (1997), it is now possible to construct a new distance function, denoted , defined as the set of shortest paths in , such that for any , we have
where is the set of all paths in between M and . By construction, it immediately follows that satisfies the triangle inequality. Therefore, forms a proper metric space.
B. Choice of SOM Dimensions
A supplemental set of simulations was conducted in order to investigate the effect of the choice of SOM dimensions on grouplevel statistical inference. We assessed the effect of rectangular SOMs, as well as the effect of increasing the dimensions of these maps. The synthetic data used for these simulations followed the design described in the section entitled Synthetic Data Simulations, based on the three different scenarios, and using the three types of SMD functions described in this paper, and setting . These results are reported in table 4.
The results of these simulations were consistent with the ones described in our first analysis of these synthetic data. In particular, for any choice of SOM dimensions, we obtained strong corroborations of the previous findings. Under both SC1 and SC2, the TSMD tended to outperform its counterparts for any choice of SOM dimensions. As before, SSMD performed poorly throughout these simulations, irrespective of the choice of SOM dimensions. Finally, the STSMD function exhibited good performance on all scenarios, and outperformed the TSMD under SC3, although STSMD did not reach significance level for this particular scenario.
Scenarios  SOM Dimensions  TSMD  SSMD  STSMD 

SC1 (Spatiotemporal)  
SC2 (Temporal)  
SC3 (Spatial)  
References
 Balding et al. (2009) Balding, D., Ferrari, P., Fraiman, R., and Sued, M. (2009). Limit theorems for sequences of random trees. TEST, 18(2), 302–315.
 Berger (1980) Berger, J. (1980). A Robust Generalized Bayes Estimator and Confidence Region for a Multivariate Normal Mean. The Annals of Statistics, 8(4), 716–761.
 Bigot and Charlier (2011) Bigot, J. and Charlier, B. (2011). On the consistency of Frechet means in deformable models for curve and image analysis. Electronic Journal of Statistics, 5, 1054–1089.
 Cordes et al. (2001) Cordes, D., Haughton, V., Arfanakis, K., Carew, J., Turski, P., Moritz, C., Quigley, M., and Meyerand, M. (2001). Frequencies contributing to functional connectivity in the cerebral cortex in ârestingstateâ data. American Journal of Neuroradiology, 22(7), 1326–1333.
 Deng (2007) Deng, D. (2007). Contentbased image collection summarization and comparison using selforganizing maps. Pattern Recognition, 40(2), 718–727.
 Dryden and Mardia (1998) Dryden, I. and Mardia, K. (1998). Statistical Analysis of Shape. Wiley, London.
 Fair et al. (2007) Fair, D., Schlaggar, B., Cohen, A., Miezin, F., Dosenbach, N., Wenger, K., Fox, M., Snyder, A., Raichle, M., and Petersen, S. (2007). A method for using blocked and eventrelated fMRI data to study “resting state” functional connectivity. NeuroImage, 35(1), 396–405.
 Fischer and Hennig (1999) Fischer, H. and Hennig, J. (1999). Neural networkbased analysis of MR time series. Magnetic resonance in medicine, 41(1), 124–131.
 Fréchet (1948) Fréchet, M. (1948). Les éléments aléatoires de nature quelconque dans un espace distancié. Annales de L’Institut Henri Poincaré, 10(4), 215–310.
 Ginestet et al. (2011) Ginestet, C., Nichols, T., Bullmore, E., and Simmons, A. (2011). Brain network analysis: Separating differences in cost from differences in topology. PLoS ONE, 6(7), e21570.
 Ginestet and Simmons (2011) Ginestet, C. and Simmons, A. (2011). Statistical parametric network analysis of functional connectivity dynamics during a working memory task. NeuroImage, 5(2), 688–704.
 González and Dasgupta (2003) González, F. and Dasgupta, D. (2003). Anomaly detection using realvalued negative selection. Genetic Programming and Evolvable Machines, 4(4), 383–403.
 Gorgolewski et al. (2011) Gorgolewski, K., Burns, C.D., Madison, C., Clark, D., Halchenko, Y.O., Waskom, M.L., and Ghosh, S.S. (2011). Nipype: a flexible, lightweight and extensible neuroimaging data processing framework in python. Front Neuroinform, 5.
 Hamming (1950) Hamming, R. (1950). Error detecting and error correcting codes. Bell System Technical Journal, 29, 147–160.
 Hardoon et al. (2007) Hardoon, D., Mourão Miranda, J., Brammer, M., and ShaweTaylor, J. (2007). Unsupervised analysis of fmri data using kernel canonical correlation. NeuroImage, 37(4), 1250–1259.
 Karcher (1977) Karcher, H. (1977). Riemannian center of mass and mollifier smoothing. Communication in Pure and Applied Mathematics, 30, 509–541.
 Kaski and Lagus (1996) Kaski, S. and Lagus, K. (1996). Comparing selforganizing maps. Artificial Neural Networks—ICANN 96, 809–814.
 Kirt et al. (2007) Kirt, T., Vainik, E., and Võhandu, L. (2007). A method for comparing selforganizing maps: case studies of banking and linguistic data. In Eleventh EastEuropean Conference on Advances in Databases and Information Systems ADBIS, Varna, Bulgaria, Technical University of Varna, 107–115.
 Kohonen (2001) Kohonen, T. (2001). SelfOrganizing Maps. Springer Series in Information Sciences. Springer.
 Kohonen et al. (2000) Kohonen, T., Kaski, S., Lagus, K., Salojarvi, J., Honkela, J., Paatero, V., and Saarela, A. (2000). Self organization of a massive document collection. IEEE Transactions on Neural Networks, 11(3), 574–585.
 Kriegel (2004) Kriegel, H. (2004). Classification of websites as sets of feature vectors. Proc. IASTED DBA, 127–132.
 Liao et al. (2008) Liao, W., Chen, H., Yang, Q., and Lei, X. (2008). Analysis of fMRI data using improved selforganizing mapping and spatiotemporal metric hierarchical clustering. IEEE Transactions on Medical Imaging, 27(10), 1472–1483.
 Mannila and Eiter (1997) Mannila, H. and Eiter, T. (1997). Distance Measures for Point Sets and Their Computation. Acta Informatica, 34, 109–133.
 MouraoMiranda et al. (2007) MouraoMiranda, J., Friston, K., and Brammer, M. (2007). Dynamic discrimination analysis: a spatiotemporal SVM. NeuroImage, 36(1), 88–99.
 MouraoMiranda et al. (2006) MouraoMiranda, J., Reynaud, E., McGlone, F., Calvert, G., and Brammer, M. (2006). The impact of temporal compression and space selection on svm analysis of singlesubject and multisubject fmri data. NeuroImage, 33(4), 1055–1065.
 Ngan and Hu (1999) Ngan, S.C. and Hu, X. (1999). Analysis of functional magnetic resonance imaging data using selforganizing mapping with spatial connectivity. Magnetic resonance in medicine : official journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine, 41(5), 939–946.
 Ngan et al. (2002) Ngan, S.C., Yacoub, E.S., Auffermann, W.F., and Hu, X. (2002). Node merging in Kohonen’s selforganizing mapping of fMRI data. Artificial intelligence in medicine, 25(1), 19–33.
 Parthasarathy (1967) Parthasarathy, K. (1967). Probability Measures on Metric Spaces. American Mathematical Society, London.
 Peltier et al. (2003) Peltier, S.J., Polk, T.A., and Noll, D.C. (2003). Detecting lowfrequency functional connectivity in fMRI using a selforganizing map (SOM) algorithm. Human Brain Mapping, 20(4), 220–226.
 Searcóid (2007) Searcóid, M. (2007). Metric spaces. Springer, London.
 Smith et al. (2004a) Smith, S.M., Jenkinson, M., Woolrich, M.W., Beckmann, C.F., Behrens, T.E.J., JohansenBerg, H., Bannister, P.R., De Luca, M., Drobnjak, I., Flitney, D.E., Niazy, R.K., Saunders, J., Vickers, J., Zhang, Y., De Stefano, N., Brady, J.M., and Matthews, P.M. (2004a). Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage, 23 Suppl 1, S208–19.
 Smith et al. (2004b) Smith, S.M., Jenkinson, M., Woolrich, M.W., Beckmann, C.F., Behrens, T.E., JohansenBerg, H., Bannister, P.R., De Luca, M., Drobnjak, I., Flitney, D.E., Niazy, R.K., Saunders, J., Vickers, J., Zhang, Y., De Stefano, N., Brady, J.M., and Matthews, P.M. (2004b). Advances in functional and structural mr image analysis and implementation as fsl. NeuroImage, 23(Supplement 1), S208–S219.
 SverdrupThygeson (1981) SverdrupThygeson, H. (1981). Strong law of large numbers for measures of central tendency and dispersion of random variables in compact metric spaces. The Annals of Statistics, 9(1), 141–145.
 Takala et al. (2005) Takala, V., Ahonen, T., and Pietikainen, M. (2005). Blockbased methods for image retrieval using local binary patterns. In Proc. 14th Scandinavian Conference on Image Analysis (SCIA, 882–891.
 Tan et al. (2005) Tan, X., Chen, S., Zhou, Z., and Zhang, F. (2005). Recognizing partially occluded, expression variant faces from single training image per person with som and soft knn ensemble. Neural Networks, IEEE Transactions on, 16(4), 875–886.
 Tarca et al. (2007) Tarca, A.L., Carey, V.J., Chen, X.w., Romero, R., and DrÄghici, S. (2007). Machine learning and its applications to biology. PLoS Comput Biol, 3(6), e116.
 Thorstensen et al. (2009) Thorstensen, N., Segonne, F., and Keriven, R. (2009). Scale Space and Variational Methods in Computer Vision, vol. 5567, chap. Preimage as Karcher Mean Using Diffusion Maps: Application to Shape and Image Denoising, 721–732. Springer.
 Tungaraza et al. (2009) Tungaraza, R., Guan, J., Rolfe, S., Atmosukarto, I., Poliakov, A., Kleinhans, N., Aylward, E., Ojemann, J., Brinkley, J., and Shapiro, L. (2009). A similarity retrieval method for functional magnetic resonance imaging (fmri) statistical maps. In Proceedings of SPIE, vol. 7259, 72590E.
 Vesanto and Alhoniemi (2000) Vesanto, J. and Alhoniemi, E. (2000). Clustering of the selforganizing map. Neural Networks, IEEE Transactions on, 11(3), 586–600.
 Wismüller et al. (2004) Wismüller, A., MeyerBäse, A., Lange, O., Auer, D., Reiser, M., and Sumners, D. (2004). Modelfree functional mri analysis based on unsupervised clustering. Journal of Biomedical Informatics, 37(1), 10–18.
 Ziezold (1977) Ziezold, H. (1977). On expected figures and a strong law of large numbers for random elements in quasimetric spaces. Transactions of the Seventh Prague Conference on Information Theory, Statistical Decision Functions, Random Processes and of the 1974 European Meeting of Statisticians.