A Neighborhood Method for Statistical Analysis of fMRI Data

In an effort to cope with the fact that functional magnetic resonance imaging (fMRI) data are spatiotemporally correlated, we propose a novel statistical method with a view to improve the detection of brain regions with increased neuronal activity in fMRI. In this method, we make use of information from neighboring voxels of a voxel, for estimation at the voxel. We examined performance of the method against the statistical parametric mapping (SPM) method using both simulated and real data. The proposed method is shown to be considerably better than the SPM in the context of receiver operating characteristics (ROC) curves.


Introduction
Functional magnetic resonance imaging data analysis are evolving quickly in a fast growing community, because of the excellent temporal and spatial resolution of these data and the innocuous aspect of their acquisition in human brains.One of the major objectives of fMRI studies is to determine which areas of the brain are activated in response to a stimulus or task.Functional imaging data are based on the principles of magnetic resonance and the fact that increases in neural activity are accompanied by changes in regional cerebral blood flow (rCBF) and blood oxygenation.This blood oxygenation level dependent (BOLD) effect is the basis for most of the fMRI studies for mapping patterns of activation in the working human brain.Lueck et al. [1] and Friston et al. [2] are among the first to develop a statistical parametric mapping package (SPM) by using MATLAB (programming language) to map patterns of activation based on a fMRI data, which are available in the public domain (www.fil.ion.ucl.ac.uk/spm).
fMRI data have a time sampling that allows for a precise estimation of the neuronal response to various experimental conditions.SPM is typically analyzed fMRI data to various conditions with a massively univariate approach, where a univariate general linear model (GLM) is fit independently at every voxel.In SPM, there is a major problem for fitting the univariate GLM, the successive fMRI scans of every voxel of the brain are not independent or spatiotemporally correlated.A solution to this problem was given by Worsley and Friston [3,4] which provided a general framework for the spatiotemporally correlated fMRI scans.In their framework, they estimate the error covariance matrix l l l of size s × s from the fMRI data, by employing the restricted maximum likelihood estimate method (ReMLE) [5].They used where s I is an identity of size s (total number of time points) of fMRI data, 1, , i s   and 1, , j s   In the V matrix, they estimate global hyperparameters l λ ( error variance components ) from the autoregressive model ( AR(1) + wn ) based on pooled data from all the interested voxels.They used this estimated error covariance matrix in SPM package in the form of ˆ W V ˆ to make independent successive fMRI scans of every voxel of the whole brain.Finally, they used this independent fMRI data for fitting the GLM of every voxel of the brain to get SPMs or activity map.
The fMRI data are temporally correlated or autocorrelated and functional imaging data have some spatial correlation.This spatial correlation is further enhanced by some operations with the analysis of SPM such as smoothing and re-slicing fMRI data, and also, fMRI data of low resolution from an individual voxel will contain some signal from the tissue around that voxel.This spatiotemporal autocorrelation Valdes-Sosa [6] has led us to propose a method in which we incorporate estimation results at neighboring voxels of a voxel, v say, in the es-F.AHMAD ET AL. 16 timation at voxel v.This proposed method is shown to improve the detection of brain regions with increased neuronal activity which is statistically more powerful as compared to the ReMLE approach of SPM.

The Proposed Method
We propose regularizing the autocorrelations based on neighboring voxel as described in details in the next Section (2.1).In practice the estimated autocorrelation parameters  vary considerably about their true values.Even when no correlation is present, it is well known that the standard deviation of the sample autocorrelation ϱ is about  2 , ξ ρ σ 1 n by [7].Some method of regularization that reduces this variability seems desirable.Purdon et al. [8] achieve this by spatial smoothing of the likelihood before estimating the autocorrelation parameters.
We analyzed the fMRI data with regard to temporal correlation, the SPM does not consider any relationship among voxels to make independent data for fitting GLM.For the efficient analysis, SPM estimates a global value of ϱ by using model AR( 1 2 can also be clearly observe in the form of bird-view.Due to these variation of the parameters and incorporates the spatial correlation between neighboring voxels, we propose an accurate pre-whitening strategy with estimation made at every voxel in the brain.We will call this strategy a neighborhood method (NH).We compare our approach by modifying the SPM code to functional imaging with the ReMLE method which has been established by Friston et al. [4] and Worsley et al. [5].

Neighborhood Algorithm
In our method, we denote an fMRI data set consisting of s time points or scans at n voxels as the s × n matrix Y.In mass-univariate GLM model [Friston et al. 1995], these data are explained in terms of a s × m design matrix X, containing the values of m regressors at s time points, and a m × n matrix of regression coefficients β, containing m regression coefficients at each of the n voxels.The model is written where is a ε n s  error matrix.1) The GLM of the variable Y at voxel from n voxels is given by where is error covariance matrix.where 3) Obtain by using the following simplest model of autocorrelation AR (1) proposed by Bullmore et al. [9] and [10][11][12][13].
In the model (3) scans are equally spaced in time and the errors from previous time point, , are mixed up with a white noise into the error of the current time point, .  Multiply the general linear model ( 2) by where a) For a neighboring voxel of , set initial values of the estimates as: where is computed based on by using model (3).
convergence takes place (e.g., c) If a voxel has more than one voxels where estimation is already made, then we take averages of the estimates from the neighboring voxels for the initial values of the parameters 1 v  ,  , V, and W. In other words, the right-hand sides of the equations in step 4a are re-placed with the corresponding averages of the estima-tes.
5) Apply step 4 to all the voxels of the brain that are involved in the data.
6) Finally, to detect the effect or activity in terms of optimal * β from step (4b), we test the null hypothesis that the effect is zero for every voxel of the brain.An effect of interest, such as a difference between stimuli, can be specified by * Ĉβ , where C is a row-vector of contrast of m conditions.It is estimated by using estimated error covariance matrix which is an exactly t-distribution because S is a square root of 2 χ distribution with degrees of freedom (p is number of parameters ) whereas Worsley and Friston [3] used Satterthwaite approximation [14] for T-statistic and effective degrees of freedom in SPM because in their ReMLE method denominator of T is not a square root of

Simulated fMRI Data
In order to evaluate our proposed method, we generated a Gaussian artificial signals in a 4-D fMRI human brain data by using the Fayyaz et al. [15] 1.Moreover, for better visualization of four ROIs, we showed three image of slices of a volume of the data in Figure 3.In this figure, the regions inside black circles demonstrated artificial activation areas in the dimension of   2 ; ~,5 , 1, ,150 where   7) and ( 8) volumes, beginning with nonactive volumes.80 acquisitions or volumes were made: 8 blocks and 10 reps in each block for rest and task conditions.In the 10 repetitions block, we assume a task stimulus, and rest conditions for the next 10 repetitions block.Successive blocks alternated between rest and task settled up a blocked paradigm as shown Figure 5, and the design matrix X of size 80 × 2 represented as.

Analysis of Simulated fMRI Data
ed and ana-w times in the first row of X for rest and task conditions and all 1's in second row used only for computation of GLM.
The simulated fMRI data were preprocess lyzed using SPM latest package (SPM8b; Welcome Department of Cognitive Neurology London) which is based on MATLAB.To avoid the errors, preprocessing steps are necessary before the GLM analysis of any real or simulated fMRI human brain data.Our simulated data  lso based on a real fMRI experiment of a subject at res-a ting state.In the preprocessing steps, (1) the whole images data were realigned according to the first volume for the correction of head motion, (2) the images were corrected for difference in slice timing, (3) images were normalized to the Montreal Neurological Institute (MNI) template using parameters defined from the normalization of the anatomical scan to the MNI template and finally, (4) images were smoothed with a Gaussian kernel of 8 mm full-width at half-maximum to reduce noise.For the rest and task conditions of the design matrix X, GLM analysis to estimate * β was performed at every voxel of this preprocessed sim ted fMRI data by using ReMlL approach of SPM and NH method.
A volume table of SPM {t} maps ula of the effect β or artificial activity from the whole data were constructed corresponding to the hypothesis "task > rest" at p < 0.001 (uncorrected) level of significance by using the SPM method as shown in Table 2 and by using the NH method as shown in Table 3.A standard volume table of SPM {t} map showed only top three significant voxel of each cluster or region with three digits of p-values and remaining all significant voxels with the full length of p-values (p < 0.001) we showed in the Figure 6.In standard SPM package (SPM8b), any hypothesis must be tested at set-level, cluster-level and voxel-level in the form of adjusted p-values.The columns in volume tables show, from left to right (1) set-level: the chance (p) of finding this (c) or a greater number of clusters in the search

 
, , x y z sents the m : coordinates in the ROIs space instead of epre-MNI space.In the tabular data the bold numbers r ost "significant" voxel within the cluster at the level p < 0.001 (uncorrected) and there are four significant cluster   4 c  which are our ROI3, ROI1, ROI4 and ROI2 respectively.The number of voxels in the activated ROIs clearly shows that the NH method increases substantially the statistical significance of the four activated regions.Finally, in the analysis of simulated data, the NH method detects, over four ROIs, 20%, 19%, 29% and 27% more activated voxels as shown in Table 3 than the SPM method as shown in Table 2.
A comparison of the two methods is also made with all the uncorrected p-values (p < 0.001) o f activated voxel which are tested at voxel-level.The SPM method detected 310 voxels whereas the NH method detected 405 voxels from the total 600 voxels of our four ROIs; the p-values of all the voxels are shown in Figure 6.The figure shows that the NH method has smaller p-values of activated voxels than the SPM.The p-values by the NH method are dispersed over a much wider range than the SPM method which is reflected in the ROC curves described in the Section (3.3) with larger true positive rate and smaller false positive rate for the NH method.

ROC Curve
We will make use characteristics (ROC of the well-known receiver operating ) curve analysis Kim et al. [16] and ctivated voxels and 1 false positiv Fayyaz et al. [15] to compare the SPM and NH methods.A ROC curve is a graphical representation of the true positive rate (TPR) or sensitivity versus the false positive rate (FPR) or (1-specificity) for a binary classifier system over a range of its discrimination threshold.The TPR and FPR values obtained by using the number of true positive (TP) and false positive (FP) activated voxels at specific threshold level in a volume of SPM {t} maps.A volume of SPM {t} maps of both methods corresponding to any p-value or threshold level (T-statistics) can be obtain with the use of SPM package (SPM8b).The TP and FP activated voxels are obtained from the volumes of SPM {t} maps of both methods by applying the lower four sequentially threshold levels as shown in Table 4.Moreover, these obtained TP voxels of simulated data demonstrated those voxels which activated within our predefined areas of ROI whereas FP voxels activated somewhere else in the brain.
In the Table 4, the measures of NH method would yield 496 true positives a e activated voxel, with a true-positive rate (sensitivity) of 496/599 = 0.8280 and a false-positive rate (1-specificity) of 1/33 = 0.0303.Similarly, 49 true positives and 9 false positives had level between 5.0 and 5.99.Thus, if any value of level greater than 4.99 were taken as an indication of NH method, this measure would yield 496 + 49 = 545 true positives activated voxels and 1 + 4 = 5 false positive activated voxels, with a true-positive rate of 547/599 = 0.9098 and a false-positive rate of 5/33 = 0.1515.And similarly so on for SPM and HM method for the other levels, 4.00 to 4.99, and 3.00 to 3.99.
We plotted true positive verses false positive rates and fitted a simple logarithmic curve ( ˆb r r y ax  ) to the four bi e of acti variate pairs of sensitivity and (1-specificity) for each method, with x r = false-positive rat vated voxels, y r = true-positive rate of activated voxels.The fitted logarithmic curves represented as the ROC curves of each method as shown below in Figure 7 where solid line curve is for the NH and dots line curve is for the SPM method.The area under the ROC curve from each of the two methods is 95.11% and 89.25% respectively for the NH and the SPM method.The ROC curve and Table 4 demonstrating the superiority of the NH method because  the NH method w activated voxel; (4) a larger area under the ROC curve by the NH than by the SPM and finally; (5) the closer the curve follows the left-hand border and then the top border of the ROC space, the more accurate the test.

Real fMRI Data
We also compared the SPM real fMRI data obtained periment, which is available in the public domain at http://cnl.web.arizona.edu/spm.htm.The three condition ("study") block images data were acquired on a GE 1 Sigma 5× Wholebody Echospeed Horizon System.The whole brain BOLD/EPI acquisition consisted of 17 slices, each 5 mm thick, with a 1 mm skip (image volume size in voxels: x = 64, y = 64, z = 17 voxel size: 3.44 mm × 3.44 mm × 6 mm; field of view (FOV) is 220).The ac-quisition took 160 seconds, with the scan-to-scan repetition time (TR) set to 2 secconds.A total of 80 acquisitions were made: 8 blocks of 20 seconds each (i.e., 8 blocks of 10 reps each).During each 20 second block, it was presented 4 stimuli of bird pictures, in each for 5 seconds.Successive blocks alternated between hard and easy learned birds and a control condition of familiar birds (crows and swans), starting with learned birds.The pattern was as in Table 5.
A "2D" T1 anatomical (same plane and section as the BOLD image; Series 2) and a "3D" image (124 Sagitta of Real fMRI Data ps of the condition ethod of the above od, we assume that each voxel is highly neighboring voxels partly due k design of a real fMRI data. B Block 8 lslices, Series 4) were also acquired to be used as structural images.

Analysis
We obtained a volume of SPM {t} ma ("easy and hard > control") for each m real fMRI data by using the same preprocessing steps (1) to (5) as discussed in the Section (3.2) corresponding to FWE-corrected p < 0.05 and contrast C = [1 -2 1] value.Results are overlaid on a volume rendered brain.The area of activated voxels of each method clearly observed in views of dorsal, lateral and ventral surface of the brain.Yellow and red areas reflect the "easy and hard" condition causing higher brain activities than the control condition as shown in Figure 8. Finally, our proposed method has more activated regions with substantially increased statistical significance, which makes it possible to decide with more confidence if a certain brain region is activated or not.

Dicussion
In the NH meth spatio-correlated with the to the smoothing of the data images and we make use of the information from neighboring voxels of a voxel, for estimation at the voxel.Whereas in the SPM, the error covariance matrix, V, does not consider any spatio-correlation among the voxels and hyperparameters of, V, also do not vary over voxels [3,4].Due this phenomenon some weak signals can not be detected with the use of ReMLE approach of SPM.On the other hand, we admit more heterogeneity of the covariance matrix or variance-covariance structure in the NH method.The distinction of the ROC curves between the NH and the SPM methods becomes more apparent for simulate data as we discussed in the Section (3.3).Worsley et al. [3] used pre-whitening strategy with a patially varying error covariance s p ated at every voxel to estimate regression coefficient β in (2) without using information from neighboring voxels.Zhang et al. [17] and [18,19] used the matrix

V
to estimate regression coefficient β of whitened model and then applied Durbin-Watson (AR(1) correlation test) on the residuals of whitened model to improve on the accuracy of the autocorrelation model.The variation of the autocorrelation coefficient (ρ) calls for the need for autocorrelation modeling with initial estimates borrowed from neighboring voxels in order to attain more accurate inferences at every voxel.
Finally, In the analysis of simulated data, the NH method detects (405 -310)/405 = 23% more activated voxels over all the ROIs as shown in Tables 2 and 3 than the SPM method, and real data detects 40% more activated voxels.The proposed method detects voxels more accurately of the simulated data even in case of weak but true artificial signals than the SPM with a better ROC performance.The true positive rate and higher activation of voxels of simulated data showed validity of the proposed method and it is more apparent as far as the real data are concerned.The ROC curves and the p-values of activated voxels show that the NH method is superior than the SPM method based on the true simulated activations.

Conclusions
In this paper, we proposed a method in an effort to cope with the r the error variance-cova fMRI data.Under timates do not change abruptly between neighboring voxels, we employed an estimate-transfusion approach between neighboring voxels by using the estimates from neighboring voxels as initial values of the estimates for their new neighboring voxels.Since the intial values may affect the final result of the estimate Wu [20] and Kim [21], it is desirable that we apply the estimate-transfusion approach to obtain the estimates that are affected by the estimates from neighboring voxels.
In both of the experiments, one with artificial data and the other with real data, we showed superiority of the proposed NH method over the traditional SPM method in the context of the ROC curve.In the M, we assume a variance-covariance structure contains hyperparameters which do not vary over voxels.This may possibly deteriorate the detection accuracy of the activated voxels when the signals are relatively small.This kind of undesirable phenomenon can be avoided by applying the NH method.
) + wn based on pooling the sample autocorrelations of the fMRI raw data that create bias result.In the fMRI data set, there were variations in the strength of correlation and SPM's global ϱ cannot adapt to these local changes.This supports the need for local autocorrelation modeling based on neighboring voxels information in order to ensure unbiased result in the form of accurate detection of activation, p-values and valid inferences at every voxel.The estimates of parameters ρ and the variance of the white noise of model (3) of all the voxels of a middle slice from real fMRI data of Section 3.3 are shown in Figure 1.Moreover the variations in the estimated variance values of all the voxel of a same slice in Figure

Figure 1
Figure 1.Estimated variance and 2 ˆξ σ ρ (histogram) val- ues of model (3) of all the voxels of a middle slice of real data.The x-axis represents voxel labels (top panel) and thê values (bottom panel).The estimates are from the NH method.

Figure 2 .
Figure 2. A bird-view of the estimated variance values of model (3) of all the voxels of a middle slice of real data.2 ˆξ σ

Figure 3 .
Figure 3. Three slice views of a volume of fMRI data and the regions inside black circles demonstrated artificial activation areas in the dimension of (x l × y l ).150 × 4 corresponding to the co-ordinates of   , , l l l x y z four ROIs.The fluctuations of weak and strong artificial signals il of Gaussian distribution with means μ = (15, 20) and variance of four different ROI are also shown inFigure4.a

Figure 5 .
Figure 5. Design paradigm of simulation data.

Figure 7 .
Figure 7. ROC curves for the NH method (solid line) and SPM method (dotted line).The markers represented the

Figure 8 .
Figure 8. Regions are rendered in yellow and red colors o the MNI template of SPM o tivated regions of the brain as a result of "higher active and hard conditions control 4 easy control 4 hard contro heterogeneity of the parameters, in particula riance structure of the noise in the the assumption that the parameter es-ReMLE approach of SP im through the National Research Found-26895). .P. Deiber, P. Cope, a, C. Kennard and S. Frackowiak, "The Colour Centre in the Cerebral Cortex of Man," N 32, 1989, pp.386-

Table 1 . Artificial activation in four different regions of interest (ROI).
il Y  is fMRI response or time points matrix of size 150 × 4 and is matrix of artificial signals of size il a Copyright © 2012 SciRes.OJBIPHY