A semiparametric PCA approach to fMRI data analysis

A SEMIPARAMETRIC PCA APPROACH TO FMRI DATA ANALYSIS M.O. Ulfarsson† and V. Solo‡ †

‡

University of Iceland, Dept. Electrical Eng., Reykjavik, ICELAND University of New South Wales, School of Electrical Eng., Sydney, AUSTRALIA ABSTRACT

Functional Magnetic Resonance (fMRI) data is most often analyzed using linear regression type methods that consider each voxel separately or by using exploratory methods such as Principal Component Analysis (PCA) or Independent Component Analysis (ICA). In this paper we introduce a model, which we call XnPCA, that combines regression with PCA. Unlike the linear regression methods XnPCA allows for non-stationary noise. Additionally, since XnPCA is based on the maximum likelihood framework the Bayesian information criterion (BIC) can be used for model selection and comparison. We compare XnPCA to a regression model commonly used in fMRI research using real data from a combined visual-motor experiment. Index Terms— Principal Component Analysis (PCA), functional Magnetic Resonance Imaging (fMRI), Regression. 1. INTRODUCTION Functional Magnetic Resonance Imaging (fMRI) has established itself to be an invaluable tool for Neuroscientists and Psychologists to understand human brain organization and function. Blood Oxygen Dependent Level (BOLD) fMRI is the dominant fMRI technique. BOLD fMRI is based on the fact that the magnetic susceptibilities of oxyhemoglobin and deoxyhemoglobin differ slightly, i.e., the signal decay rate of deoxyhemoglobin is more rapid than its oxygenated counterpart. Following neural activation, oxygen rich blood flows to the area of activation which leads to a local changes in the magnetic resonance decay parameter T2∗ . This subsequently leads to a local rise in the intensity of the observed signal. The intensity rise, which is called the BOLD response, can be detected [1] and used to generate a BOLD-weighted image. There are two main approaches to analyze fMRI data: univoxel and multivoxel approaches. In univoxel approaches the time series recorded at each brain volume element (voxel) are analyzed separately and are typically modeled as a sum of the BOLD response, drift terms, and stationary noise such as AR(p), or ARMA(p,q). The drift terms along with the stationary noise are supposed to account for low-frequency noise, head movement and various other effects of physiological or hardware related origin.

Multivoxel methods such as Principal Component Analysis (PCA) and Independent Component Analysis (ICA) are usually used in fMRI as exploratory tools [2, 3]. However, there is an underlying model for PCA which we call noisy PCA [4]. The papers [5–7] demonstrate application of nPCA for fMRI analysis. In the following, we propose a semiparametric model based on nPCA which we call XnPCA that combines the traditional unvoxel regression model with nPCA. The most important feature of this model is that is allows for nonstationary noise. We propose to use BIC to select the number of nPCs in the model. Finally, we compare, using the Bayesian information criterion (BIC), XnPCA to a typical univoxel regression model. A model somewhat similar to XnPCA was discussed in [8], but used in completely different way, i.e., in a multi-subject setting. The paper is constructed as follows. In Section 2 we describe the univoxel regression model. Section 3 introduces the XnPCA model and section 4 describes the maximum likelihood estimation of the XnPCA parameters. Section 5 treats model selection and section 6 model comparison. Then in section 7 we apply the new model on a real fMRI data set. Finally, in section 8 conclusions are drawn. Throughout the paper, boldface denotes vectors or matrices. I r denotes the r × r identity matrix; N (0, Ω) denotes multivariate Gaussian density with zero mean and covariance Ω; i.i.d. denotes independent and identically distributed; |Ω| is the determinant of Ω. 2. THE UNIVOXEL REGRESSION MODEL A typical univoxel regression model used for fMRI signal analysis is given by y v = Xβv + η v ,

v = 1, ..., M.

Here y v is a T × 1 vector of the voxel time series observed at voxel number v, X is a T × p regression matrix that contains a baseline, third order polynomial bases that accounts for drift and the BOLD response. The noise η v is assumed to be zeromean Gaussian AR(2) process. Details regarding estimation, inference, and model selection can for example be found in [9].

3. THE XNPCA MODEL

It can be shown that this cyclic descent algorithm is globally convergent. We can use the result of the algorithm to estimate the the sample nPCs and the corresponding covariance:

The XnPCA model is given by y v = Xβ v + Guv + ǫv ,

v = 1, ..., M.

Here y v and X have the same meaning as in the univoxel iid regression model above, G is a T × r loading matrix, uv ∼ iid N (0, I r ) are the nPCs, ǫv ∼ N (0, σ 2 I T ) represents noise, uv and ǫv are mutually uncorrelated. The log-likelihood function for the XnPCA model is given by lθ (y) = −

M M tr(S e Ω−1 ) − log |Ω|, 2 2

ˆv u

ˆ ˆ T (y v − X β) ˆ −1 G = W

S uˆ

ˆ −1 , = σ ˆ2 W

ˆ +σ ˆ =G ˆTG ˆ2 I r . where W 5. MODEL SELECTION FOR XNPCA To finish the model fitting we need to choose the number of nPCs r. Here we use the Bayesian information criterion (BIC) [11]. In its general form BIC is given by

where θ = (β, G, σ 2 ), Ω = GGT + σ 2 I T , and

1 BIC = −lθˆ (y) + dim(θ) log M, 2

M 1 X Se = (y − Xβ v )(y v − Xβv )T , M v=1 v

ˆ is the maximum likelihood estimate of the parameter where θ vector θ, and dim(θ) is the number of free parameters. The BIC for the XnPCA model (scaled and irrelevant parameters excluded) is given by

is the data covariance matrix. 4. ML ESTIMATION FOR XNPCA

1X T −r 1 log lj + log σ ˆ2 + dim(θ) log M, 2 j=1 2 2M r

The task at hand is to estimate θ. Differentiating the likelihood and solving the Euler equations leads to ˆ β v ˆ G σ ˆ

2

= = =

(X T Ω−1 X)−1 X T Ω−1 y v , 2

v = 1, ..., M (1)

1/2

P r (Lr − σ ˆ Ir) R tr(S e ) − tr(Lr ) , T −r

(2) (3)

where P r is the T × r matrix of unit eigenvectors of S e , Lr = diag(l1 , ..., lr ) contains the corresponding eigenvalues, R is an arbitrary r × r rotation matrix. The Euler equations (1)-(3) are linked nonlinearly together so the solution is not trivial. We suggest the following cyclic ascent algorithm [10] to maximize the likelihood: Algorithm 1 (XnPCA). The XnPCA algorithm consists of following steps: Initialization: Set β v0 = [¯ y v , 01×p−1 ]T , v = 1, ..., M and 2 compute G0 and σ0 according to (2) and (3).

where dim(θ) = T r − r(r − 1)/2 + r + 1 + M p. The number of nPCs r is selected according to the minimum of BIC (or close to it). 6. MODEL COMPARISON The BIC statistic can be used to compare different models given that they are scaled the same way and same irrelevant terms are excluded (e.g., the 2π factor in Gaussian models). We are interested in comparing XnPCA with the univoxel regression model of section 2. The BIC for the univoxel regression model has the following form BIC =

Iteration: For k = 1, ... compute 2 β-step: Given Gk−1 and σk−1 compute β k using (1).

Ω-step: Given β k compute Gk and (3).

BICr =

σk2

using (2) and

Termination: Stop when absolute relative difference between lθk−1 (y) and lθk (y) is below 10−6 and call the ˆ ˆ σ final estimates G, ˆ 2 , and β.

M 1 X 1 log |Σv | + (T p + 1) log(M ). 2M v=1 2

where Σv is the AR(2) covariance matrix. This expression is used to compute the BIC in the results section below. 7. RESULTS The fMRI data comes from a motor-visual experiment [12]. The subject was instructed to do right hand thumb opposition

268

315

267.5 310 267 305

266.5

BICr

BICr

266 265.5 265 264.5

300 295 290

264 285 263.5 263 0

5

10

15

20

25

280 0

30

5

10

15

r

Fig. 1. BIC for brain slice number 5.

X = [1T , d1 , d2 , d3 , α1 , α2 , α3 ] where the terms d1,t = t, d2,t = t2 and d3,t = t3 , t = 1, ..., T are included to capture drift terms that increase with time. They were normalized and mean corrected to avoid numerical instabilities. The terms α1 , α2 and α3 are models for the BOLD response corresponding to the stimulus pattern for the finger-thumb opposition, annular, and anti-annular checkerboards, respectively. The model used for the BOLD response is a simple parametric model from [13]. The XnPCA method was applied to each of the 21 brain slices and for each brain slice BIC was used to select the number of nPCs r. Fig. 1 shows BIC for brain slice #5 and Fig. 2 depicts BIC for brain slice #19. Their minimum occurs at r = 8 and r = 17, respectively. Fig. 3 displays the number of nPCs for each of the 21 brain slices. It is interesting that there is a significant difference in the number of nPCs between brain slices. The BIC statistic for the univoxel regression model presented in section 2 are compared to the XnPCA BIC values in Table 1. For all the brain slices the BIC for XnPCA is lower and therefore XnPCA should be considered the better method for this example. An activation map for XnPCA can be constructed by computing the signal to noise ratio

v

25

30

Fig. 2. BIC for brain slice number 19. 22 20 18

r (# nPCs)

and observe 8Hz annular and anti-annular flashing checkerboard according to predetermined stimulus pattern. A 3T scanner was used to obtain T = 100 observations on 21 brain slices with sample interval TR=2s. The voxel size is 3.5 × 3.5 × 5mm. The regression matrix used for the XnPCA model and the univoxel regression model is the following

ˆ ˆ T X T X αβ β α,v α,v α , SNRv = ˆ k2 ky − X β

20

r

16 14 12 10 8 6 0

5

10

15

20

25

# brain slice Fig. 3. The number of nPCs for each brain slice. 2.5 10 2 20

30

1.5

40

1

50 0.5 60 10

20

30

40

50

60

Fig. 4. An activation map for brain slice # 5. v = 1, ..., M,

v

ˆ where X α = [α1 , α2 , α3 ], and β α,v are the β parameters

corresponding to X α . Fig.4 shows an activation plot for brain slice #5 clearly showing activation in the motor cortices.

Table 1. Model comparison between the univoxel regression method and XnPCA using BIC. # brain slice Regression XnPCA 1

300.97

273.68

2

283.30

264.19

3

278.88

263.24

4

275.90

262.40

5

276.47

263.04

6

275.01

262.26

7

274.33

262.66

8

274.41

262.46

9

275.42

263.75

10

275.08

263.43

11

277.26

264.53

12

282.19

268.93

13

282.79

270.31

14

287.43

275.56

15

290.51

277.85

16

291.36

276.61

17

297.40

281.68

18

304.54

288.72

19

299.04

284.65

20

292.15

279.95

21

292.46

283.98

8. CONCLUSIONS In this paper we presented a new model, which we call XnPCA, that combines the univoxel regression method and nPCA. The main difference between a standard univoxel regression model and XnPCA is that XnPCA allows for nonstationary noise. We compared a typical univoxel regression model with XnPCA by comparing the BIC statistic for real fMRI data showing that XnPCA outperformed the univoxel regression method. 9. REFERENCES [1] S. Ogawa, D.W. Tank, R. Menon, J.M. Ellerman, S.G. Kim, and H. Merkle, “Intrinsic signal changes accompanying sensory stimulation: Functional brain mapping with magnetic resonance imaging,” Proc. Natl. Acad. Sci. USA, vol. 89, pp. 5951–5955, 1992. [2] K. J. Friston, C. D. Frith, P. F. Liddle, and R. S. Frackowiak, “Functional Connectivity: The principal compo-

nent analysis of large data sets,” J. Cereb. Blood Flow Metab., vol. 13, pp. 5–14, 1993. [3] W. Xiong, Y. Li, T. Li, T. Adali, and V.D. Calhoun, “On ica of complex valued fMRI: Advantages and order selection,” in Proc. IEEE International Conference on Aucoustics, Speech and Signal Processing (ICASSP’01), Las Vegas, Nevada, USA, 2008. [4] D.N. Lawley, “A modified method of estimation in factor analysis and some large sample results,” in Uppsala symposium on psychological factor analysis, Uppsala, Sweden, 1953, pp. 35–42. [5] M.O. Ulfarsson and V. Solo, “Smooth principal component analysis with application to functional magnetic resonance imaging,” in Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP’06), Toulouse, France, May 2006, vol. 2, pp. II–993 – II–996. [6] M.O. Ulfarsson and V. Solo, “Spatially local and temporally smooth PCA for fMRI,” in Proc. IEEE International Conference on Image Processing (ICIP’06), Atlanta, Georgia, USA, 2006, pp. 2853–2856. [7] M.O. Ulfarsson and V. Solo, “Sparse variable principal component analysis with application to fMRI,” in Proc. IEEE International Symposium on Biomedical Imaging (ISBI’07), Washington D.C., 2007, pp. 460–463. [8] C.J. Long, E.N. Brown, C. Triantafyllou, I. Aharon, L.L. Wald, and V. Solo, “Nonstationary noise estimation in functional MRI,” Neuroimage, vol. 28, pp. 890–903, 2005. [9] R.H. Shumway and D.S. Stoffer, Time Series Analysis and its applications, Springer, New York, NY, second edition, 2006. [10] D.P. Bertsekas, Nonlinear Programming, Athena Scientific, Belmont, MA, second edition, 1999. [11] G. Schwartz, “Estimating the dimension of a model,” Ann. Stat., vol. 6, no. 2, pp. 461–464, 1978. [12] Robert Cox, “AFNI,” Bethesda, Md, National Institute of Mental Health NIMH, http://afni.nimh.nih.gov/afni/. [13] M.S. Cohen, “Parametric analysis of fMRI data using linear system methods.,” NeuroImage, vol. 6, pp. 93– 103, 1997.