Statistics and Its Interface Volume 3 (2010) 33–43 Adaptive statistical parametric mapping for fMRI Ping Bai † , Haipeng Shen ‡ , Jianhua Z. Huang § and Young K. Truong †,∗ Brain activity is accompanied by changes in cerebral blood flow (CBF) and the differential blood oxygenation that are detectable using functional magnetic resonance imaging (fMRI). The process of identifying brain activation regions can be facilitated by estimating the hemodynamic response function (HRF). There have been some remarkable new developments in statistics to handle this problem. In this paper, we introduce a novel procedure which is capable of adapting itself to any of the existing methods by improving its performance through the application of a penalized smoothing technique. Using a computer experiment and a real fMRI data set, the proposed procedure is assessed by comparing its performance very favorably to the popular SPM based method. AMS 2000 subject classifications: Primary 62H35, 62G08; secondary 62H12. 1. INTRODUCTION Functional magnetic resonance imaging (fMRI) is a non- invasive medical method for detecting the temporal hemody- namic response to a given stimulus. This is possible because of the differential magnetic susceptibility of oxygenated and deoxygenated hemoglobins. Thus the response is also re- ferred to as the blood oxygenation level dependent (BOLD) effect. The data is described in terms of 3 dimensional grids of M volume elements (voxels), which is scanned repeat- edly over a period of time so that each voxel will possess (say) N time point BOLD effects. In practice, we denote the BOLD effect of the j th voxel recorded at time i by y ij , where 1 ≤ i ≤ N and 1 ≤ j ≤ M . Here M is a scanner res- olution specific constant. Typically, M ∼ 10 5 and N ∼ 10 2 . Based on this vast amount of data, one of the objectives of fMRI is to identify which voxel has responded to the exper- imental stimulus. Model-driven strategies have been playing a dominant role in relating the stimulus to the hemodynamic response. Among them, the most widely used approach is based on the general linear model (GLM) given by (1.1) Y = Xβ + ǫ, ∗ Corresponding author. † Research partially supported by NSF DMS-0707090 ‡ Research partially supported by NSF DMS-0606577 § Research partially supported by NSF DMS-0606580 where Y = (y ij ) N×M is the BOLD effect data, X = (x ik ) N×(K+1) is related to the stimuli, β =(β kj ) (K+1)×M is a matrix of parameters, and ǫ =(ǫ ij ) N×M whose columns consist of noise series with mean zero and covariance σ 2 Σ. Suppose there are K different conditions involved in the experiment. For example, these conditions correspond to four different hand movements included in one of our ex- periments of the fMRI study [15]. Then X is an N × (K +1) design matrix, where the first column is usually constant representing the condition when the subject is at rest with no experimental activities performed, and each of the re- maining columns contains the time component for an ex- periment condition. The (K + 1) × M matrix β contains the related parameters, and due to the special formulation of the design matrix, the second to the last rows naturally give the contrast between each experiment condition and the benchmark when the subject is at rest. The inferences on each row of β will generate a spatial map related to one experiment condition. The construction of the spatial map is carried out in two steps. In step one, a time series regression method is carried out at each voxel in order to obtain the estimate of the effect of the stimulus on the BOLD signals. That is, we regress each voxel time series Y j (a column vector of Y) on X via (1.2) Y j = Xβ j + ǫ j , ǫ j ∼ (0,σ 2 Σ), j =1, 2,...,M. Here β j is the j th column of β. Once this model has been fitted at each voxel, inferences of the model parameters are then made according to the experiment hypothesis. The re- sulting statistics from all the voxels are assembled spatially into an image, which is the so-called statistical parametric map. The second step then focuses on the analysis of the statistical map in order to identify those areas of the brain that are activated by the stimuli [4]. For more detailed dis- cussions, see statistical parametric mapping (SPM [3]), FM- RIB Software Library (FSL [16]) and AFNI [2]. The nature of the BOLD response implies that in areas of activation there is a delayed and blurred version of the stimulus sequence (Figure 1). Hence each column of the de- sign matrix X, which represents the temporal characteristic of the expected response, is commonly modeled through the