Elsevier

NeuroImage

Volume 37, Issue 3, 1 September 2007, Pages 855-865
NeuroImage

Local Linear Discriminant Analysis (LLDA) for group and region of interest (ROI)-based fMRI analysis

https://doi.org/10.1016/j.neuroimage.2007.04.072Get rights and content

Abstract

A post-processing method for group discriminant analysis of fMRI is proposed. It assumes that the fMRI data have been pre-processed and analyzed so that each voxel is given a statistic specifying task-related activation(s), and that individually specific regions of interest (ROIs) have been drawn for each subject. The method then utilizes Local Linear Discriminant Analysis (LLDA) to jointly optimize the individually-specific and group linear combinations of ROIs that maximally discriminates between groups (or between tasks, if using the same subjects). LLDA tries to linearly transform each subject's voxel-based activation statistics within ROIs to a common vector space of ROI combinations, enabling the relative similarity of different subjects' activation to be assessed. We applied the method to data recorded from 10 normal subjects during a motor task expected to activate both cortical and subcortical structures. The proposed method detected activation in multiple cortical and subcortical structures that were not present when the data were analyzed by warping the data to a common space. We suggest that the method be applied to group fMRI data when warping to a common space may be ill-advised, such as examining activation in small subcortical structures susceptible to mis-registration, or examining older or neurological patient populations.

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:Yk=Xkβk+εk,andCov(εk)=σk2Vkwhere 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:WkVkWkT=I(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:Yk*=Xk*βk+εk*where the superscript *denotes the whitened quantities. The whitening matrix Wk is estimated by the residuals εk and Vk as:Wk=Vk12.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):βˆkGLS=(Xk*TXk*)1Xk*TYk*Cov(βˆkGLS)=σk2(Xk*TXk*)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:H0:cβk=0where 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)):βcont=Xgβg+εgwhere Xg=(10100101) 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:H0:βg1βg2=0.

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:t(k,V)=[t(k,1,v1),t(k,2,v2),,t(k,N,vN)]T,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:fk=1Bti=1Btt(K,Vi).The above process is repeated Bf times (Bf is several hundreds or thousands) and all feature vectors are then collected:Fk=[fk(1),fk(2),,fk(Bf)],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 matrixX=[F1,F2,FS].The 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)

  • D.L. Collins et al.

    Design and construction of a realistic digital brain phantom

    IEEE Trans. Med. Imag.

    (1998)
  • Damasio, H., 2005. Human Brain Anatomy in Computerized...
  • Cited by (14)

    • Separating 4D multi-task fMRI data of multiple subjects by independent component analysis with projection

      2013, Magnetic Resonance Imaging
      Citation 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, Neuroscience
      Citation 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, Neuroscience
      Citation 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.

    View all citing articles on Scopus
    View full text