Local Linear Discriminant Analysis (LLDA) for group and region of interest (ROI)-based fMRI analysis
Introduction
Group analysis in fMRI is typically done in several consecutive steps. First, fMRI data are corrected for motion, despite the fact that most methods cannot easily distinguish changes in fMRI signal from that induced by motion (Liao et al., 2005, Liao et al., 2006). Data are then spatially transformed to a common space such as the atlas by Talaraich (Talairach and Tournoux, 1988) or the probabilistic space suggested by the Montreal Neurological Institute (Collins et al., 1998) to minimize intersubject differences. However, because of the variability in human brain anatomy, the inter-subject registration is typically imperfect, so spatial low-pass filtering (“smoothing”) is performed to de-emphasize anatomical differences (Friston, 1996). Once data have been motion corrected, warped to a common space, and spatially smoothed, the task-related activation of a voxel of a subject k is estimated with linear regression techniques:where Yk is the Tk × 1 time course of the voxel, Xk is the Tk × D design matrix containing the hypothesized activation (often incorporating estimates of the hemodynamic response function) as well as other covariates, εk is the Tk × 1 vector of residuals, σk2 is the homogeneous variance of the residuals, and Vk is the correlation matrix. The subscript k indicates that all the variables are related to subject k.
As fMRI data are typically not temporally white, data are often pre-whitened using a whitening matrix Wk such that:(for an excellent summary the reader is referred to: Mumford and Nichols, 2006). If each term in Eq. (1) is pre-multiplied by Wk, we have:where the superscript *denotes the whitened quantities. The whitening matrix Wk is estimated by the residuals εk and Vk as:The regression estimates of Eq. (3) can then be estimated by Ordinary Least Squares (OLS) to give the Generalized Least Squares estimate of Eq. (1):
Contrasts between conditions are of most interest in an experiment, e.g., contrasting BOLD signal during performance of a given task compared to rest. In the study on a single subject, the null hypothesis is that the contrast between the least-squared estimates is zero:where c is the contrast row vector. For example, if we are interested in the comparison between task 1 and task 2, c is [1, − 1].
Group analyses are usually done using a Summary Statistics method, which is a two-staged approach; first individual models are fit to each subject as described above, and then a second level is applied to make group inferences on the cβk (Mumford and Nichols, 2006). In the usual situation where one is contrasting activation across two groups, the second level is a multivariate regression equation with the design matrix encoded with group inclusion indicators (Fig. 1(a)):where is a binary K × 2 matrix coded to show group inclusion (K is the number of subjects from the two groups), βcont is composed of the contrasts cβk for each individual as defined in the first stage, βg = [βg1, βg2]T is mean activation of the two groups and εg ∼ N(0,δg2Vg) is the residual with the variance δg2 and the correlation matrix Vg being a diagonal matrix, typically just I. Here the null hypothesis is that the group activations for a given voxel in the common spatially transformed space are not significantly different:
A number of different implementations have been proposed to implement the above analysis in a practical way. The fMRIStat method uses Restricted Maximum Likelihood (ReML) to estimate σg2 (Worsley et al., 2002), then smoothes the data to increase its degree of freedom and accuracy and finally tests the hypothesis with t-statistics. The SPM2 package (Friston et al., 2002a, Friston et al., 2002b) estimates the δ2V + δg2Vg term with ReML under a simplifying assumption that all the subjects share a common covariance matrix δk2Vk = δ2V, and then tests the hypotheses with F statistics. The FMRIB software library estimates σg2 with the maximum a posteriori (MAP) criteria, then screens obviously insignificant voxels with Z-statistics and finally performs a Bayesian inference on the significance of the remaining voxels with a slower but more accurate Markov Chain Monte Carlo (MCMC) simulation (Beckmann et al., 1998).
Nevertheless, there are a number of shortcomings with the previously described methods. The above methods work on the voxel level—this assumes that after suitable spatial transformation, there is a perfect correspondence between the same voxel across subjects. While this may be mitigated somewhat by spatial smoothing, such low-pass filtering degrades the spatial resolution of the data. Activation estimates in small, subcortical structures such as the basal ganglia or thalami, which abut functionally different tissues (e.g., the internal capsule), may be particularly affected by mis-registration errors.
One way to partially circumvent the difficulties associated with spatially transforming functional maps to a common space is to manually draw anatomical regions of interest (ROIs) for each subject, and performing analyses at the ROI level—as opposed to the individual voxel level. Using standard atlases, a particular brain region (e.g., the lateral cerebellar hemisphere) is manually circumscribed on the high-resolution structural MRI scans that have been co-registered with the functional data, and the voxels within this region are analyzed. The benefit of this method is that it does not require rigid spatial transformation, preventing possible gross distortion of a particular brain area, as may occur if the anatomy of a given individual differs significantly in size and shape to the homologous area in the exemplar brain. However, drawing ROIs is labor-intensive, subject to human error, and requires the assumption that a functionally active region (the SMA for example) of a given brain will be within an anatomically standardized index (i.e., Broadman's Area 9) which is used to draw the ROI.
In addition to the possibilities of mis-registration, the previously described voxel-based methods do not explicitly model interactions between brain regions. Covarying regions are often of interest, but are not included in the group methods described above. Conceptually, group methods are done in two stages: in the first stage, individually specific regression models are fit to the data; and then the results of these models are used in a group-level analysis. Because the goal of these methods is to test a specific hypothesis, these methods may be conducted sequentially. In contrast, if the goal is to find which combination of brain regions is maximally different between tasks, it is desirable to jointly optimize the individual statistical model and the overall models simultaneously.
In individual-subject fMRI analysis, in addition to hypothesis driven methods, there is a role for data driven methods, such as Independent Component Analysis (ICA) (McKeown et al., 1998, Calhoun et al., 2003), which do not need rigorous a priori specification of activation patterns. In an analogous manner, there may be particular interest in discovering the combinations of brain regions (specified by ROIs) that are maximally contrasted during performance of certain tasks (Fig. 1(b)). There is therefore a need for a multivariate, discriminant analysis approach that works at the region of interest (ROI) level as opposed to the individual voxel level.
Previous work has taken individual activations (or the t-statistics associated with them) and used a multivariate discriminant approach (McKeown and Hanlon, 2004). In order to apply a discriminant approach, we first assume that some statistical analysis has been performed to assign a t-statistic, related to activation, to each voxel, such as simple t-test or regression (e.g., t-statistic of a voxel comparing mean BOLD signal for a given condition to BOLD signal for that voxel at rest). Since most multivariate discriminant approaches assume Gaussian distributions, we extract features from the t-statistics with bootstrap methods to ensure the Gaussian assumption. Consider a subject k, we randomly select a voxel from within each of its N ROIs, and assemble the result into a column t-statistic vector:where t(k, r, vr) is the t-statistic of the vrth voxel in the rth ROI of subject k and the voxel index V = [v1, v2,…,vN]T. The random selection is repeated a number of times, say Bt times, and the ith draw is notated as t(k, Vi). After a reasonable number of draws (e.g., Bt = 30), the normalized mean of the t-statistic vectors is taken as a feature to ensure the data can be modeled as multivariate Gaussian:The above process is repeated Bf times (Bf is several hundreds or thousands) and all feature vectors are then collected:where fk(i) is the ith random sample of a feature vector. This process is then repeated for all S subjects in the groups and all the Fk's are concatenated into a big matrixThe benefit of formulating the problem in this way is two fold; the Fk's will be asymptotically normally distributed (by the multivariate central limit theorem; Flury, 1997) even if the voxel-based statistics are not, and that p-variate linear analyses can now be performed on X (Flury, 1997).
The variability of the amplitude of activation across subjects is well known and typically this is dealt with by using a random-effects analysis (Mumford and Nichols, 2006). In random-effects models, the data are assumed to be derived from a hierarchy of different populations whose differences are constrained by the hierarchy. Each subject from a different group (e.g., group of normal subjects) is considered to be representative of the entire population of normal subjects.
In some cases the magnitude of inter-subject differences in fMRI activation can exceed the task-specific differences within individuals. In order to deal with this situation, yet still maintaining the benefits of linear discriminant analysis, we propose using a recently developed solution, Local Linear Discriminant Analysis (LLDA), that was initially designed to solve a somewhat different, but still related, problem (Kim and Kittler, 2005). In discriminant analysis trying to classify images of faces, often the difference in discriminant features of different poses of the same face can greatly exceed the difference in discriminant features between faces. Finding a classifier that is sensitive to images from different subjects, yet insensitive to different poses from the same subject, has been problematic. Kim and Kittler (2005) suggested LLDA as a solution to this problem. We therefore propose using LLDA to sensitively discriminate between task-dependent ROI-based patterns of activity, while being relatively robust to the differences between subjects (Fig. 1(b)).
We apply this method to data derived from a motor paradigm that would be expected to activate cortical and subcortical structures. We show that the proposed method, consistent with prior neuroscience knowledge derived from animal models, detects significant group activation in subcortical structures that was not present when the same group of data were analyzed using standard methods utilizing spatial normalization.
Section snippets
Experimental methods
To demonstrate the proposed method, we utilized fMRI data that would be expected, based on prior knowledge, to activate subcortical structures. We enrolled 10 healthy volunteers, 5 males and 5 females (range 27–45 years). All subjects were right hand dominant, and had normal neurological examinations, had no history of neurological disease, and were not currently using any psychoactive prescription medications. Handedness was determined according to Edinburgh Handedness Inventory. The paradigm
Results
Using SPM 5, we found activation in the left primary motor cortex during right hand movement, and similarly left primary motor cortex activation during right hand movement. Activation in the left cerebellar hemisphere, right somatosensory cortex, and right thalamus was detected using left hand movement only (Fig. 3(a)).
Similarly, during right hand movement, significant activation was detected with LLDA in the left primary cortex. With left hand movement LLDA detected significant activity in the
Discussion
We have utilized a modification of the local linear discriminant (LLDA) algorithm to perform group-wise analysis on activation maps that have been manually labeled with individually drawn Regions of Interest (ROIs). In addition to finding the contralateral primary motor cortical activation and ipsilateral cerebellar activation similar to the SPM approach, we found a number of regions that were significantly active, including bilateral dorsolateral prefrontal cortex during externally guided
Acknowledgments
This work was supported by a National Parkinson's Foundation Center of Excellence grant (MJM), NSERC Grant CHRPJ 323602-06 (MJM), CIHR grant CPN-80080 (MJM), NSERC grant RGPIN 312490-05 (ZJW), a Parkinson's fellowship grant from His Excellency Sheik Abdulmohsin Al-Sheikh (MML) and NIH grants AG21491 (XH) and RR00046 (GCRC).
References (25)
- et al.
Classical and Bayesian inference in neuroimaging: applications
NeuroImage
(2002) - et al.
Classical and Bayesian inference in neuroimaging: theory
NeuroImage
(2002) - et al.
Group analysis in functional neuroimaging: selecting subjects using similarity measures
NeuroImage
(2003) - et al.
An automated method for neuroanatomic and cytoarchitectonic atlas-based interrogation of fMRI data sets
NeuroImage
(2003) - et al.
A post-processing/region of interest (ROI) method for discriminating patterns of activity in statistical maps of fMRI data
J. Neurosci. Methods
(2004) - et al.
Region of interest based analysis of functional imaging data
NeuroImage
(2003) Prefrontal cortex function, quasi-physiological stimuli, and synaptic plasticity
J. Physiol. (Paris)
(2003)- et al.
A general statistical analysis for fMRI data
NeuroImage
(2002) - et al.
General multilevel linear modeling for group analysis in FMRI
NeuroImage
(1998) - et al.
ICA of Functional MRI Data: An Overview
Design and construction of a realistic digital brain phantom
IEEE Trans. Med. Imag.
Cited by (14)
Separating 4D multi-task fMRI data of multiple subjects by independent component analysis with projection
2013, Magnetic Resonance ImagingCitation Excerpt :Complementary to univariate approaches, multivariate statistical methods consider brain voxels or regions of interest as an integrated network and jointly analyze them [3,4]. Among various multivariate approaches, methods that have been successfully applied to the group analysis of fMRI data are independent component analysis (ICA) [5–7], partial least squares (PLS) [8,9], local linear discriminant analysis [10], independent vector analysis [11] and support vector machines [12]. With the exception of ICA and PLS, the application of most of these methods is limited to single-task fMRI group analyses.
Joint amplitude and connectivity compensatory mechanisms in Parkinson's disease
2010, NeuroscienceCitation Excerpt :Subsequent analysis was restricted to voxels that had a t-value >1.96 in at least one of the task frequencies. The LLDA method is fully described elsewhere (McKeown et al., 2007). In brief, the method involves repetitively sampling t-statistics from each ROI to create feature vectors, equal in length to the number of ROIs, to create a matrix.
Levodopa-sensitive, dynamic changes in effective connectivity during simultaneous movements in Parkinson's disease
2009, NeuroscienceCitation Excerpt :A full description of the procedures for calculating SEM and MAR networks can be found in the technical appendix. A description of the procedure for LLDA has been published elsewhere (McKeown et al., 2007). While it would also be of interest to calculate the second-order MAR, that is, to predict current fMRI signals from not just the prior time point, but the previous two time points, this would require estimation of twice the number of parameters.
Analysis of Alzheimer’s Disease Based on the Random Neural Network Cluster in fMRI
2018, Frontiers in Neuroinformatics