Journal of Neuroscience Methods 240 (2015) 101–115
Contents lists available at ScienceDirect
Journal of Neuroscience Methods journal homepage: www.elsevier.com/locate/jneumeth
Computational Neuroscience
Nonparametric variogram modeling with hole effect structure in analyzing the spatial characteristics of fMRI data Jun Ye a,∗ , Nicole A. Lazar b , Yehua Li c a b c
Department of Statistics, University of Akron, Akron, OH, United States Department of Statistics, University of Georgia, Athens, GA, United States Department of Statistics & Statistical Laboratory, Iowa State University, Ames, IA, United States
h i g h l i g h t s • The new periodic variogram model can well describe the fluctuating spatial structure of fMRI data. • The hole effect model considers both the nonlinear physical relationship between the proximate voxels and the functional relationship between distant voxels.
a r t i c l e
i n f o
Article history: Received 17 August 2014 Received in revised form 10 November 2014 Accepted 11 November 2014 Available online 18 November 2014 Keywords: Bessel model Cross-validation Kriging Mean squared deviation ratio
a b s t r a c t When analyzing functional neuroimaging data, it is particularly important to consider the spatial structure of the brain. Some researchers have applied geostatistical methods in the analysis of functional magnetic resonance imaging (fMRI) data to enhance the detection of activated brain regions. In this paper, we propose a nonparametric variogram model for the complicated spatial characteristics of fMRI data. The new periodic variogram model can well describe the fluctuating spatial structure of fMRI data, considering both the nonlinear physical relationship between the proximate voxels and the functional relationship between distant voxels. We demonstrate the effectiveness of the new variogram model using fMRI data from a saccade study. © 2014 Elsevier B.V. All rights reserved.
1. Introduction The advent of modern neuroimaging techniques has revolutionized our understanding of the working of the human brain. The complicated dependence structures in the data obtained from these imaging methods provide challenges for analysis, but at the same time open a new window for the exploration of the spatial and temporal characteristics of activation patterns in the brain. In recent years there has been explosive growth in the number of neuroimaging studies using functional magnetic resonance imaging (fMRI). In fMRI experiments, researchers record the blood oxygenation level dependent (BOLD) contrast at locations across the subject’s brain to see which areas are involved in performing a given task or are responding to a particular stimulus. BOLD imaging takes advantage of inherent differences in the magnetic properties of oxygenated
∗ Corresponding author. Tel.: +1 3309728008. E-mail addresses: [email protected] (J. Ye), [email protected] (N.A. Lazar), [email protected] (Y. Li). http://dx.doi.org/10.1016/j.jneumeth.2014.11.008 0165-0270/© 2014 Elsevier B.V. All rights reserved.
and deoxygenated hemoglobin, hence fMRI has become a very powerful research tool for studying, albeit indirectly, neuronal activity. According to Kosslyn (1999), there are two questions addressed in neuroimaging: one is when the particular structures and processes are involved in the brain activations over time; the other is how the brain implements the information processing. Solid statistical analysis plays a key role in finding answers to these questions. A brain image is a matrix of numbers corresponding to spatial locations. We generally refer to each volume element in the image as a voxel, which is the three dimensional analogue to a pixel. Whereas structural MR images generally comprise a single three dimensional image, fMRI data are represented as a time series of such images. As a result, it is particularly important to consider the spatial and temporal correlation in the data and see the relationship of different areas in the brain. Geostatistics is a branch of applied statistics that focuses on providing quantitative descriptions of natural variables distributed in space or time. In geostatistics, researchers use variogram modeling to obtain predictions from data on spatially or temporally variable properties. Nowadays geostatistics is widely used in many fields of science and
102
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
industry where there is a need for evaluating spatially or temporally correlated data. Because of the spatial and temporal properties of the fMRI data, researchers have started to apply geostatistical methods to enhance the detection of activated regions of the brain. The variogram is a traditional geostatistical tool that provides quantification of the degree of directional spatial properties in the measured random variables (Rossi et al., 1992). The variogram gives an intuitive interpretation for the variations of the data in a secondorder stationary process (Schabenberger and Pierce, 2002); it only compares the average square difference between locations; and one need not to know the constant population mean (Ye et al., 2009). For these reasons, it is common to work with the variogram rather than the covariance function in geostatistical applications. Spence et al. (2007) use a Gaussian variogram model to find proximate voxels of interest in fMRI data, and show that such spatial analysis can identify regions of the brain exhibiting statistically significant group differences. They claim that close neighbors tend to have similar intensities and distant voxels show less similarity. They further show that neighboring voxels have strong nonlinear correlation so that the sample variogram values of the data are well fitted by the Gaussian variogram model. Their geostatistical analysis method offers important advantages over more traditional voxel-by-voxel approaches. However, their method only considers the close physical neighbors of the voxels. In fMRI data analysis, different regions in the brain may be functionally related, which means that the spatial correlation between voxels may not necessarily decrease with their increasing geometric distances. For example, the different areas involved in language processing will be functionally related, even though they may be physically separated. Bowman (2007) notes that distant voxels may exhibit correlation and proposes a novel way to describe the spatial structure of the fMRI data. He defines two different measures of relationship: one is the physical relationship, which is the geometric separations between two voxel locations and is the regular lag distance used in Spence et al.’s (2007) variogram model; the other is the functional relationship, which is the signal change between two voxels (i.e., the difference in the measured signals between two voxels) and is also called anatomical distance. As the variogram is defined as a function of the lag distance in geostatistics, Bowman uses the anatomical distance instead of the regular geometric distance as the lag distance in his new variogram model. He shows that the proposed variogram can find the functional relationship of the voxels. However, it is worth pointing out that this new defined lag distance in Bowman’s proposed variogram contains unknown measurement errors. A basic assumption in defining the variogram is to ignore an error in the lag distance and to assume the error only exists in the measurement of the signal (Stein, 1999). This means that the lag distance is assumed to be known without error in the measurement process. In fMRI data, the geometric distance between two voxel locations can be measured precisely. However, the signal values of the voxels always contain measurement errors and furthermore the errors may vary in different voxels. Hence it is very hard to well define the anatomical measure between two voxels. Although Bowman states that empirical values from unspecified auxiliary data can correct the anatomical measure, he does not explicitly describe how to deal with the measurement errors among the voxels and his approach might generate other additional problems, e.g., sometimes we do not know if the empirical values from the selected auxiliary data are well-suited for the current analysis. Our proposed model is also related to the random field model in fMRI data analysis. Random field theory is a branch of mathematics that has been extensively used in describing the spatial correlation of brain image data. In work by Worsley and colleagues (Worsley and Friston, 1995; Worsley et al., 1996, 1998, 2002; Friston et al., 2004), Gaussian or other monotonically homologous random field models describe the spatial correlation between neighboring
voxels. However, it is important to consider the spatial correlation among distant voxels too. In this paper, we will build an advanced periodic random field model that can well describe the fluctuant spatial structure of the random filed in the brain image data. Compared with the existing variogram modeling methods in fMRI data analysis (Bowman, 2007; Spence et al., 2007), the main contribution of the paper is that we propose a new nonparametric variogram model for the complicated spatial characteristics of fMRI data. Our hole effect model directly considers both the nonlinear physical relationship between proximate voxels and the functional relationship between distant voxels. Therefore it can effectively capture the main features of the functional relationship among different voxels. The rest of the paper is organized as follows: In Section 2, we introduce basic concepts and definitions of the new geostatistical method. Application of our method to a real fMRI data set is in Section 3. Section 4 discusses the proposed method and gives the conclusion. 2. Methods 2.1. Definition of the variogram Consider a slice of the fMR image at a specific time point; the signals are obtained over a regular grid of voxels, i.e. measured Z(s) : s ∈ D , where s = (x, y) denotes the coordinates of the sample site (x, y); Z(s) is a particular voxel value at location s, D is the set of the region of interest {s1 , · · · , sn }. For a general spatial model Z(s) = + (s), where is an unknown constant population mean and (s) is the zero mean random function at s, Z(s) is second order stationary if E{Z(si )} = E{(si + h)} = and Cov{Z(si + h), Z(si )} = C(h) for any si + h, si ∈ D, where lag h is a directional distance between two locations. Under second order stationarity, the variogram is defined as an alternative measure of spatial dependence (Schabenberger and Gotway, 2005), i.e. (h) = 12 Var{Z(s + h) − Z(s)}. Since (h) = C(0) − C(h), the variogram can reach a limiting maximum value called the sill, i.e., (h0 ) = C(0). Here lag h0 is known as the correlation range with C(h0 ) = 0, because it defines the average distance within which the random function remains correlated spatially. The variogram summarizes spatial relations in the data; it can be calculated directly from the data, in which case it is called the empirical or sample variogram. The computational formula for the empirical variogram (h) ˆ is (Isaaks and Srivastava, 1989). (h) ˆ =
1 2N(h)
{Z(si ) − Z(sj )}2 ,
(1)
(si ,sj )|si −sj =h
where s1 , . . ., sn in D are lying on an already known regular lattice, and the summation is over only the N(h) pairs of voxels which are separated by lag h. Since a relatively large N(h) is associated with a relatively small standard error (Morris, 1991; Rossi et al., 1992), we choose the number of data pairs to be at least 30 (Journel and Huijbregts, 1978). The assumption that the mean of the spatial model is constant is not realistic for fMRI data. Hence before the variogram calculations, we must first remove any trend surface (Schabenberger and Gotway, 2005). After trend removal, we consider Z(s) to be second order stationary (Cressie, 1993). In practice, we usually assume that the model only has a low-order trend, and the error term has the high-order structure. For example, we set up the model as Y (s) = f (s) + + (s) = f (s) + Z(s),
(2)
where f(s) is a trend function. After removing the trend f(s), we render the residuals stationary and use Z(s) to compute the variogram
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
103
for further analysis, and add back the trend f(s) at the end of the analysis. For a specific spatial structure, if the variogram function is direction independent, it is called isotropy (Chiles and Delfiner, 1999). We consider our detrended data as isotropy which is validated by rose diagram (Banerjee et al., 2005; Isaaks and Srivastava, 1989; Waller and Gotway, 2004) of the empirical variogram.
periodicity of variograms, highlighting that the variogram with hole effect structure can be used to describe the complicated spatial structure in the fMRI data too. In the next section, we build a nonparametric variogram model with the geometric lag distance to describe the complex relationship of the different voxels (Tobler, 1970; Gringarten and Deutsch, 2001).
2.2. Properties of the variogram
2.3. Variogram modeling
By the assumption of second order stationarity, the variogram defines the relationship between the variance and the lag distance. Hence variance is a function of lag alone. Understanding the characteristics of different variations of the voxel values at different locations of the brain is very important in fMRI data analysis. Based on the special properties of the fMRI data, we consider two types of relationship in the variogram, which have been described already in Spence et al. (2007) and Bowman (2007).
It is useful to model the empirical variogram because the results in a smooth fit provide estimators at arbitrary lags (Webster and Oliver, 2007; Schabenberger and Gotway, 2005). Current geostatistical practice in selecting a variogram model is often rather subjective, relying on empirical guidelines (Gorsich and Genton, 2000). By the special properties of the fMRI data, we examine two kinds of variogram model approaches: one is parametric, the other is nonparametric. The parametric approach, proposed by Spence et al. (2007), considers the physical relationship of the voxels only; our new nonparametric approach considers both the physical relationship and the functional relationship of the voxels.
2.2.1. Physical relationship and functional relationship In the mining and petroleum industries, as the magnitude of the lag separation vector increases, the variogram increases too, showing some linear properties. This is generally observed in nature and is called the physical relationship since regions close in space tend to have similar signal values. In brain image data, a similar relationship exists and is pointed out by Spence et al. (2007). They find that the spatial correlation between neighboring voxels is much stronger than that in many other data types, therefore the variogram approaches the origin with a somewhat reversed curvature, showing quadratic patterns (Spence et al., 2007). Hence the physical relationship may more properly be described as nonlinear. The finding from Spence et al. (2007) leads to an improved approach in analyzing the regions of interest in fMRI data. We will explore this property in more detail via a Gaussian variogram model. In fMRI data, not only the physical neighbors of the voxels are related, but also distant voxels with similar anatomy or function are related in some ways. Bowman (2007) observes this phenomenon in the fMRI data however Spence et al. (2007) ignore it. Bowman (2007) uses a monotonically increasing variogram model with an anatomical lag to describe the functional relationship of the voxels, where the regular lag distance is changed from geometric distance to anatomical distance. However, as we discussed in the Introduction, the anatomical lag always contains unknown measurement errors. 2.2.2. Hole effect relationship We consider the geostatistical method for model selection and model inference in the analysis of fMRI data. In some instances the variogram may seem to fluctuate more or less periodically, rather than increasing monotonically. That is, the variogram can decrease from its maximum to a local minimum and then increase again, which appears as a “down-hole”. This fluctuant structure of the variogram is called the hole effect (Journel and Huijbregts, 1978; Webster and Oliver, 2007). Hole effect structure of the variogram is usually ignored in geostatistics, because in the mining and petroleum industries from which most geostatistical analysis arose, the fluctuations in the empirical variogram are considered as more or less random noise (Webster and Oliver, 2007). But for some data (e.g., fMRI data), if we have already known that the fluctuations in the variogram are due to correlation among different distant regions, the hole effect structure is not artificial and should not be ignored. Furthermore, the hole effect structure can provide valuable information concerning functional variability of the data (Pyrcz and Deutsch, 2007). Several recent studies from ecological and environmental science (e.g., Radeloff et al., 2000; Chen, 2005) have revealed relations between periodic landscape patterns and
2.3.1. Parametric variogram modeling The parametric variogram is monotonically increasing. Wackernagel (2003) defines the Gaussian-type model in a variogram form c
(h) = 2 [1 − exp{−(ah) }],
(3)
where a is a parameter to determine the weights of the neighbors and 1 ≤ c ≤ 2. As c → 1, the model is more linear at small lags and the spatial correlation between nearby points is weaker. As c → 2, the model approaches the origin with a more quadratic shape and appears more sigmoidal, and the spatial correlation between nearby points is higher. The two extreme cases are: when c = 1, Eq. (3) becomes an exponential model; when c = 2, it becomes a Gaussian model. The pure Gaussian model is unstable, because the quadratic behavior at the origin in the variogram will generate extreme values at the borders of the estimated map in prediction (kriging) (Wackernagel, 2003). Hence, a nugget effect is usually added as a nested structure in the Gaussian model to yield a discontinuity at the origin; this avoids the extreme extrapolation properties and makes the kriging results more stable. In geostatistics, the name of nugget effect was coined by Matheron (1962) in reference to small nuggets of ore distributed throughout a larger body of rock. The nugget effect can be thought of as white noise, usually due to micro-scale variation or measurement error (Schabenberger and Gotway, 2005). It can be expressed in a variogram form as (h) = 2 for h > 0. Since the sum of valid covariance functions is itself a valid function, the Gaussian model used in the parametric approach becomes
2
(h) = 2 + ( 2 − 2 ) 1 − exp{−(ah) } .
(4)
2.3.2. Nonparametric variogram modeling Bessel model: We choose the Bessel model to describe the functional relationship because it can well describe the pronounced hole effect structure in the variogram function (Schabenberger and Gotway, 2005; Webster and Oliver, 2007). The covariance function of the Bessel model is defined as C(h) = 2
d 2
(bh)
−
J (bh),
(5)
where b is the number of sign changes, d is the dimension of the data, = (d/2) − 1 and J (·) is the Bessel function of the first kind of order (Yaglom, 1987; Cressie, 1993; Stein, 1999; Webster and
104
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
(a) Two sign changes (one period)
(b) Three sign changes (one period) 1.5 Variogram
Variogram
1.5
1
0.5
0
0
2 3 Lag (c) Four sign changes (two periods)
1
0.5
0
2
0
1
2 Lag
3
0
1
2 3 Lag (d) Five sign changes (two periods)
k=0
k!2
2
0
1
.
2 Lag
3
(6)
(h) = 2 1 − J0 (bh) ,
(7)
where 2 is the sill. Fig. 1 gives examples of Bessel variogram models (h) = 1 − J0 (bh) with lag h ∈ [0, 3.8] and parameter b = 2, 3, 4, 5, where the parameter b determines the number of sign changes in the model function. The number of sign changes (periods) monotonically increases as b increases. Nonparametric model approach: To increase the flexibility in modeling the variogram, we choose a linear combination of valid covariance functions. Treating the valid covariance functions as though they are basis functions yields a nonparametric model fitting approach (Ecker and Gelfand, 1997; Schabenberger and Gotway, 2005; Shapiro and Botha, 1991). The nonparametric variogram model is then defined as
2
2
2
(h) = + ( − )
1−
p
(9)
2
p−1
2
k=1
0.5
The Bessel variogram model used in our analysis is
,
1 1 c (h) = + ( − ) 1 − J0 (kbh) − exp{−(ah) } . (10) p1 p2
Oliver, 2007). When d = 1, C(h) = cos(bh); when d = 2, C(h) = J0 (bh); when d = 3, C(h) = (1/bh) sin(bh); when d→ ∞, C(h) = exp { − (bh)2 } (Stein, 1999; Schabenberger and Gotway, 2005). The periodic function has a weaker hole effect structure as d increases, and it becomes a Gaussian function when d→ ∞. Detailed proofs are in Stein (1999). We consider random field is of dimension two in what follows (i.e., d = 2). As J is the Bessel function of the first kind of order and = (d/2) − 1, the order is fixed as = (d/2) − 1 =0 in the data analysis (Webster and Oliver, 2007; Schabenberger and Gotway, 2005). Hence C(h) = J0 (bh), which can be expressed as ∞ 2k (−1)k bh
p
where 2 is the nugget effect, 2 is the sill, b is the number of sign changes of the basis function, p is the number of basis functions. Because of the nonlinear property of the fMRI data, we further replace one Bessel basis function with a Gaussian-type basis function to get more flexibility in the model fitting, i.e., 2
Fig. 1. Examples of Bessel variogram models (h) = 1 − J0 (bh). Panels (a)–(d) are the different variograms with parameter b = 2, 3, 4, 5 under lag h ∈ [0, 3.8]. The hole effect structure of the model changes with b.
J0 (bh) =
1 1− J0 (kbh) p
2
k=1
1
0
2
(h) = + ( − )
0.5
1.5 Variogram
Variogram
1.5
1
0
1
is also suggested that each basis function is equally weighted, i.e., ωk = 1/p, and k is chosen as k to simplify model fitting. Now the variogram model becomes
ωk J0 ( k bh)
,
(8)
k=1
where p is the number of basis functions, and ωk and k are weights. Theoretically the more basis functions we use, the better the fit (Cherry et al., 1996). However in practice it is too computationally intensive for large data sets to greatly increase the number of basis functions. Also, as the number of basis functions increases and the fit gets better, the model is less generalizable and less interpretable. When the Bessel function J0 (·) is considered as the basis function, Ecker and Gelfand (1997) suggest that p should be at most five. It
The present model is a hybrid between the Bessel basis functions and the Gaussian-type basis function. In the Gaussian-type basis function, the parameter a determines the weights of neighbors for the target voxel, and takes account of the physical relations in the voxels; in the Bessel basis functions, the parameter b determines the different sign changes of the target voxel, and hence considers the functional relations in the voxels; the parameters p1 and p2 satisfy 4/p1 + 1/p2 = 1, which allows for the different weights of the Gaussian-type model and Bessel model in this nonparametric fitting approach. Variogram fitting method: Methods for variogram fitting include maximum likelihood and least squares. Maximum likelihood estimation has an assumption of Gaussian distribution of the data and the fit has to be based on all data; by contrast, least squares does not have the distribution assumption and one can restrict the lag distance for the data. We choose least squares here because of its flexibility (Schabenberger and Gotway, 2005; Schabenberger and Pierce, 2002). Since fitting models in this way is a form of nonlinear regression, we use the Levenberg-Marquardt algorithm to obtain parameter estimates of the variogram (Marquardt, 1963). It is a challenging problem to well fit the empirical variogram because it is hard to match all its peaks and troughs by the variogram model. Manual fitting is simple and is done most commonly in geostatistics (Ma and Jones, 2001). However, it is preferable to use automatic fitting if we have to model many variograms at the same time. In our setting, we take the approach of “automatic fitting at first, and manual fitting afterwards”. We first fit the empirical variogram automatically by using plausible starting values and then check the value of the coefficient of multiple determination (R2 ) in the least squares method. If the R2 value is low (e.g. less than 0.8), we adjust the starting values of the initial parameters in the Levenberg–Marquardt algorithm and see if the fit is improved. If the fit is still not satisfactory, we manually fit the variogram by a procedure that embodies both visual inspection and statistical fitting, and follow these three rules (Webster and Oliver, 2007; Ma and Jones, 2001): (i) ignore the point-to-point fluctuation and concentrate on general trends; (ii) estimate the variogram to be accurate at short lags, with less accuracy at longer lags; (iii) match the cyclic pattern of the variogram at least to the first peak or trough. 2.4. Kriging To compare variograms we must work in the predictive space (Ecker and Gelfand, 1997). Kriging is a standard tool in geostatistics for spatial prediction. We choose the best variogram model by cross-validation in kriging. For each observation in the validation data set, the mean square error (MSE) could be calculated after kriging. However, MSE is insensitive to the variogram used for kriging, and kriging variance is sensitive to the variogram (Lark, 2000,
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
2009). Hence mean squared deviation ratio (MSDR) is more suitable to assess the performance of the different variogram models in estimating the functional or physical relations among the voxels (Ripley, 1981; Davis, 1987; Webster and Oliver, 2007). MSDR is the MSE divided by the kriging variance, hence it can be considered a weighted MSE. Dividing the error by MSE allows one to compare the magnitudes of both the actual and the predicted error in the cross-validation. If the model for the variogram is accurate then the MSE should be equivalent to the kriging variance; and so the MSDR should be 1. 2.4.1. General characteristics of kriging For predicting the value of the random variable Z(s0 ) at an unsampled site s0 , from the data Z(s1 ), . . ., Z(sn ) at sampled sites s1 , . . ., sn , kriging is given as the weighted sum of the values of the neighbors. The weights are calculated by minimizing the error variance of an assumed variogram model for the data with regard to ˆ 0 ) at the spatial distribution of the observed data points. When Z(s an unsampled site s0 is estimated from Z(s1 ), . . ., Z(sn ) at sampled sites s1 , . . ., sn , Z(si ) is decomposed as a geostatistical model Z(si ) = (si ) + e(si ), i = 0, . . ., n,
(11)
where (si ) is an unknown function, e(si ) is the zero mean secondorder stationary random function at spatial position si . n ˆ 0) = Ordinary kriging takes Z(s w Z(si ), subject i=1 i to
n
i=1
wi = 1.
Letting
2
n
m
be
the
Lagrange
parameter,
ˆ 0 )} − 2m{
− 1} is minimized by setting its E{Z(s0 ) − Z(s i=1 i partial derivatives with respect to wi and m equal to zero. Hence, (si , s0 ) =
n
wj (si , sj ) − m,
subject to
j=1
n
wi = 1.
(12)
i=1
The ordinary kriging variance is 2 OK =
n
wi (xi − x0 ) + m.
(13)
i=1
Here we first estimate the linear trend, then analyze the residuals of the trend by ordinary kriging, and finally add back the trend after finishing the analysis (Cressie, 1993). 2.4.2. Model selection by cross-validation We use K−fold cross-validation to choose between different variogram models based on their kriging results. In our setting, for voxel values Z(s1 ), . . ., Z(sN ) at locations s1 , . . ., sN , K−fold crossvalidation splits the N voxels into K folds, where the kth fold has k ) for k = 1, . . ., K. The voxel values in voxel values Z(s1k ), . . ., Z(sN/K each fold are removed in turn and new values are predicted at those locations by the voxel values in other folds using kriging. ˆ k ) be the kriging values predicted at locations ˆ k ), . . ., Z(s Let Z(s 1
N/K
k 2 (sk ), . . ., 2 (sk , and OK ) be the kriging variance. s1k , . . ., sN/K OK N/K 1
Then
2 1 k ˆ k ) / 2 (sk ). Z(si ) − Z(s OK i i N K N/K
MSDR =
(14)
k=1 i=1
When MSDR >1, the kriging variance underestimates the true estimation variance; when MSDR 8 s). We find that the mean values of MSDR in the four time
1.4
1.2
1
0.8
0.6
0.4
0.2
Gau−100
Gau−150
Gau−200
Three−basis
Four−basis
Fig. 8. MSDR values in the 28 time points under different approaches.
Five−basis
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
111
Fig. 9. Panel (a) is the 32 × 26 rectangular map and its 16 equally sized pieces. Panel (b) is a resampled map in bootstrap under Anti task. Panel (c) is a resampled map in bootstrap under Pro task.
Table 2 In the different time groups, MSDR mean values under the Gaussian-type model approach and the nonparametric model approach. Time points
Gau-100
Gau-150
Gau-200
Three-basis
Four-basis
Five-basis
61, 65, 69, 73, 77, 81, 85 62, 66, 70, 74, 78, 82, 86 63, 67, 71, 75, 79, 83, 87 64, 68, 72, 76, 80, 84, 88
0.4447 0.2909 0.3085 0.2959
0.4736 0.2914 0.3625 0.3774
0.7480 0.5017 0.6795 0.5881
0.5757 0.6133 0.6203 0.7149
0.9714 0.6772 0.7918 0.6506
0.9879 0.8783 0.8990 0.8162
Mean
0.3350
0.3762
0.6293
0.6311
0.7727
0.8954
groups are quite similar (Table 2), and the results are consistent with those from the 28 time points. That is, in every time group, the MSDRs are different under the following three comparisons: Gaussian-type versus nonparametric approach; Gau-200 model fitting versus Five-basis functions model fitting; and Gau-100 model fitting versus Gau-200 model fitting.
3.6. Further analysis in a bootstrap resampling method To further assess the differences of the models, we performed a bootstrap analysis. In this method, we cut every 32 × 26 rectangular map into 16 equally sized pieces (panel (a) in Fig. 9). Since there are 156 images, one for each time point, where 73 of them are under Anti task and 83 of them are under Pro task, we randomly combine the different pieces of the images according to the two different tasks, where there are 1400 resamples for each task (panels (b) and (c) in Fig. 9). The variogarm modeling and kriging are executed for the resampled images. We compare our Five-basis model with the Gau-200 model. The averages of MSDR and its standard error are 0.9374 and 0.0342 for the Five-basis model, and are 0.8168 and 0.0594 for the Gau-200 model (Table 3). The bootstrap results confirm the efficacy of our proposed method.
Table 3 In the bootstrap resampling method under Anti and Pro tasks, MSDR and its standard error (stder) with the Gau-200 model approach and the Five-basis model approach, where there are 1400 resamples for each task. Task
MSDR of Gau-200
Stder of Gau-200
MSDR of Five-basis
Stder of Five-basis
Pro Anti
0.7212 0.9124
0.0449 0.0739
0.9386 0.9361
0.0352 0.0331
Mean
0.8168
0.0594
0.9374
0.0342
4. Conclusions and future work In many geostatistical studies, the variogram is monotonically increasing; any hole effect structure is considered to be artificial and is ignored. However in fMRI data, high correlation exists between a given voxel and neighboring voxels as well as some relatively distant voxels. Hence the cyclical pattern of the variogram is very important for fMRI data analysis and it cannot be ignored. This is a major difference between geostatistical and fMRI data. In geostatistics, the physical relationship is most relevant since close regions are more related than distant regions, and distant regions will tend to be unrelated. But in fMRI, voxels may be functionally related even though they are not neighbors. This claim is similar to that in Bowman (2007), but our approach is totally different. Instead of using the anatomical measure in the variogram, we use the regular geometric measure to describe the functional relationship between voxels. The advantage of our method is that it does not need auxiliary data, and considers the original fMRI data directly. When we use the hole effect variogram model in model fitting, the presence of such a hole effect provides valuable information concerning complicated spatial variability of the fMRI data. In summary, the use of the hole effect structure in the variogram model provides a promising approach in the analysis of brain imaging data. It shows a good example of applying geostatistical ideas in a new area. However, most of the variogram models with hole effect structure are with dampening, which means that the functional relations are always weaker than the physical relations. But for a target voxel in fMRI data, its distant voxels may exhibit higher correlation than some closer neighbors (Bowman, 2007). How to well define a suitable hole effect variogram model and keep the balance between the physical relationship and the functional relationship in the model fitting will be considered in future work. For applying the method to the whole brain, conceptually, the procedure of the analysis would be the same as the analysis for a
112
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
single slice. Computationally, the extension to the whole brain is more of a challenge, since we quickly run into memory and storage issues. We will address efficient computation in future work.
Appendix A.
Acknowledgement
As suggested by one referee, we also try the fifth slice in the data analysis. After the data preprocessing steps (Fig. A.1), we consider the different variogram models for the empirical variograms at time points 61–88 (Fig. A.2). The best model will be chosen by
We wish to express our appreciation to Dr. Rebecca L. McNamee for kindly providing the saccade data set.
A.1. Data analysis in the fifth slice
Fig. A.1. Image maps at the fifth slice. Panel (a) is the mean map of all 156 time points. Panel (b) is the masked mean map. Panel (c) is the map at a randomly selected time point. Panel (d) is the mean subtracted map. Panel (e) is the linear trend map. Panel (f) is the trend removed map. The same processing procedure is carried out at every time point.
Fig. A.2. Panel (a) is the empirical variogram map in 0◦ , 90◦ and 270◦ directions at a randomly selected time point. Panel (b) is the empirical variogram maps in 0◦ , 90◦ and 180◦ directions at the same time point.
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
113
Fig. A.3. Maps for 5-fold cross-validation. Panels (a)–(e) are the five folds where the voxels denoted by white squares (144 voxels for (1), (2) and (3); 145 voxels for (4) and (5)) are removed in each fold for estimation. Panel (f) is the original map with 722 voxels.
Table A.1 MSDR for different time points (61–88) under the Gaussian-type model approach and the nonparametric model approach. Time point
Gau-100
Gau-150
Gau-200
Three-basis
Four-basis
Five-basis
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
0.5387 0.0631 0.4386 0.3762 0.7270 0.7289 0.3787 0.6260 0.3242 0.1525 0.6497 0.6355 0.3040 0.4838 1.1951 0.5643 0.6891 0.4665 0.5140 0.0758 0.2059 0.4237 0.7218 0.2859 0.8263 0.3955 0.0430 0.2303
0.7064 0.0456 0.5244 0.5389 0.8670 0.9909 0.4783 0.7154 0.3338 0.1110 0.7592 0.9643 0.3365 0.6369 1.1367 0.8192 0.6080 0.6127 0.6246 0.0539 0.2312 0.5833 0.0383 0.2891 1.2053 0.5584 0.0280 0.2840
0.5388 0.5097 1.0137 1.1245 1.0280 1.1017 0.6970 0.5122 0.3751 0.3099 1.0733 0.6384 0.4626 1.0047 0.6214 0.7779 1.1955 0.6061 0.9579 0.4168 0.2875 0.9995 0.2709 0.2477 0.1777 0.5315 0.2184 0.3682
0.6105 0.8059 0.2914 0.8072 0.5815 0.5850 0.8182 0.2191 0.8270 0.3281 1.0299 0.4657 0.7225 0.6977 0.8812 0.9769 0.5769 0.1933 0.2752 0.4512 0.6611 0.9451 0.4399 0.6937 0.4994 0.5505 0.9559 0.4822
0.4997 0.6191 0.3082 0.4753 0.4691 0.6505 0.7905 0.9759 1.2151 0.4163 0.7412 0.4701 0.7582 0.8377 0.6007 1.1788 0.9498 0.9141 1.2208 0.8534 0.6175 0.8721 0.4314 0.8536 0.7149 0.5987 1.0430 0.6214
0.6829 0.9862 1.0175 0.6117 0.6281 1.1862 1.1842 0.9334 0.6272 0.8674 0.8248 0.6855 0.6347 1.1562 1.1175 1.0726 1.1221 0.8923 0.6726 0.9266 0.8332 0.9964 0.6474 0.7938 0.8884 0.8881 0.8543 0.6215
Mean
0.4662
0.5386
0.6452
0.6204
0.7497
0.8697
Table A.2 In the different time groups, MSDR mean values under the Gaussian-type model approach and the nonparametric model approach. Time points
Gau-100
Gau-150
Gau-200
Three-basis
Four-basis
Five-basis
61, 65, 69, 73, 77, 81, 85 62, 66, 70, 74, 78, 82, 86 63, 67, 71, 75, 79, 83, 87 64, 68, 72, 76, 80, 84, 88
0.5150 0.3877 0.5630 0.3991
0.6216 0.5056 0.5128 0.5236
0.5807 0.7233 0.6932 0.5837
0.6398 0.5865 0.6702 0.5851
0.7463 0.7476 0.7292 0.7755
0.7738 0.9961 0.9026 0.8064
Mean
0.4662
0.5386
0.6452
0.6204
0.7497
0.8697
114
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
1.2
1
0.8
0.6
0.4
0.2
0
Gau−100
Gau−150
Gau−200
Three−basis
Four−basis
Five−basis
Fig. A.4. MSDR values in the 28 time points under different approaches.
5-fold cross-validation in kriging (Fig. A.3). The MSDR results in Tables A.1 and A.2 and Fig. A.4 show that the variogram model with hole effect structure has better fit. The conclusion is consistent with results obtained using the fourth slice.
References Banerjee S, Carlin BP, Gelfand AE. Hierarchical modeling and analysis for spatial data. Boca Raton, FL: Chapman & Hall CRC Press; 2005. Bowman FD. Spatiotemporal models for region of interest analyses of functional neuroimaging data. J Am Stat Assoc 2007;102:442–53. Carroll MK, Cecchi GA, Rish I, Garg R, Rao AR. Prediction and interpretation of distributed neural activity with sparse models. NeuroImage 2009;44:112–22. Chen X. Statistical and geostatistical featuures of streambed hydraulic conductivities in the Platte River. Nebraska, Environ Geol 2005;48:693–701. Cherry S, Banfield J, Quimby WF. An evaluation of a nonparametric method of estimating semi-variograms of isotropic spatial processes. J Appl Stat 1996;23:435–49. Chiles JP, Delfiner P. Geostatistics: modeling spatial uncertainty. New York, NY: John Wiley & Sons, Inc; 1999. Cressie NAC. Statistics for spatial data. revised ed. New York, NY: John Wiley & Sons, Inc; 1993. Davis BM. Use and abuse of cross-validation in geostatistics. Math Geol 1987;19:241–8. Ecker MD, Gelfand AE. Bayesian variogram modeling for an isotropic spatial process. J Agric Biol Environ Stat 1997;4:347–69. Friston KJ, Frith CD, Dolan RJ, Price CJ, Zeki S, Ashburner JT, et al. Human brain function. 2nd ed. London UK: Elsevier; 2004. Gorsich DJ, Genton MG. Variogram model selection via nonparametric derivative estimation. Math Geol 2000;3:249–70. Goutte C, Hansen LK, Liptrot MG, Rostrup E. Feature space clustering for fMRI meta analysis. Hum Brain Mapp 2001;13:165–83. Goutte C, Toft P, Rostrup E, Nielsen FA, Hansen LK. On clustering fMRI time series. NeuroImage 1999;9:298–310. Gringarten E, Deutsch CV. Teacher’s aide variogram interpolation and modeling. Math Geol 2001;33:507–34. Guitton D, Buchtel HA, Douglas RM. Frontal lobe lesions in man cause difficulties in suppressing reflexive glances and generating goal-directed saccades. Exp Brain Res 1985;58:455–72. Haas TC. Lognormal and moving window methods of estimating acid deposition. J Am Stat Assoc 1990;85:950–63.
Hastie T, Tibshirani R, Friedman J. The elements of statistical learning: data mining, inference, and prediction. New York, NY: Springer Science; 2001. Huettel SC, Song AW, McCarthy G. Functional Magnetic Resonance Imaging. Sunderland, MA: Sinauer Associates Inc; 2004. Isaaks EH, Srivastava RM. An introduction to applied geostatistics. New York, NY: Oxford University Press, Inc; 1989. Journel A, Huijbregts ChJ. Mining Geostat. New York, NY: Academic Press, Inc; 1978. Kosslyn SM. If neuroimaging is the answer what is the question. Philos Trans R Soc Lond B: Biol Sci 1999;354:1283–94. Lark RM. A comparison of some robust estimators of the variogram for use in soil survey. Eur J Soil Sci 2000;51:137–57. Lark RM. Kriging a soil variable with a simple nonstationary variance model. J Agric Biol Environ Stat 2009;14:301–21. Lindquist MA. The statistical analysis of fMRI data. Stat Sci 2008;23:439–64. Ma YZ, Jones TA. Teacher’s aide modeling hole-effect variogram of lithologyindicator variables. Math Geol 2001;33:631–48. Marquardt DW. An algorithm for least-squares estimation of nonlinear parameters. J Soc Ind Appl Math 1963;11:431–41. Matheron G. Traite de Geostatistique Appliquee, Tome I, Memoires du Bureau de Recherches Geologiques et Minieres No. 14. Paris: Editions Technip; 1962. Morris MD. On counting the number of data pairs for semivariogram estimation. Math Geol 1991;23:929–43. Pham TD, Eisenblatter U, Baune BT, Berger K. Preprocessing film-copied MRI for studying morphological brain changes. J Neurosci Methods 2009;180: 352–62. Pyrcz MJ, Deutsch CV. The whole story on the hole effect; 2007 www. sciencedirect.com/science/article/pii/S016890020802010X Radeloff VC, Miller TF, He HS, Mladenoff DJ. Periodicity in spatial data and geostatistical models: autocorrelation between patches. Ecography 2000;23:81–91. Ripley BD. Spatial statistics. New York, NY: John Wiley & Sons, Inc; 1981. Rossi RE, Mulla DJ, Journel AG, Franz EH. Geostatistical tools for modeling and interpreting ecological spatial dependence. Ecol Monogr 1992;62:277–314. Schabenberger O, Gotway CA. Statistical methods for spatial data analysis. Boca Raton, FL: Chapman & Hall/ CRC Press; 2005. Schabenberger O, Pierce FJ. Contemporary statistical models for the plant and soil sciences. Boca Raton, FL: Chapman & Hall/CRC Press; 2002. Shapiro A, Botha JD. Variogram fitting with a general class of conditionally nonnegative definite functions. Computational statistics and data analysis 1991;11:87–96. Sommer FT, Wichert A. Exploratory analysis and data modeling in functional neuroimaging. Cambridge, MA: The MIT Press; 2003. Spence JS, Carmack PS, Gunst RF, Schucany WR, Woodward WA, Haley RW. The accounting for spatial dependence in the analysis of SPECT brain imaging data. J Am Stat Assoc 2007;102:464–73.
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115 Stanberry L, Nandy R, Cordes D. Cluster analysis of fMRI data using dendrogram sharpening. Hum Brain Mapp 2003;20:201–19. Stein ML. Interpolation of spatial data, some theory for Kriging. New York, NY: Springer-Verlag; 1999. Tobler W. A computer movie simulating urban growth in the Detroit region. Econ Geogr 1970;46:234–40. Wackernagel H. Multivariate geostatistics, an introduction with applications. 3rd ed. New York, NY: Springer-Verlag; 2003. Waller LA, Gotway CA. Applied spatial statistics for public health data. Hoboken, NI: John Wiley & Sons, Inc; 2004. Wang K, Jiang T, Yu C, Tian L, Li J, Liu Y, et al. Spontaneous activity associated with primary visual cortex: a resting-state fMRI study. Cerebr Cortex 2008;18:697–704. Webster R, Oliver MA. Geostatistics for environmental scientists. 2nd ed. Chichester, England: John Wiley & Sons, Ltd; 2007.
115
Worsley KJ, Friston KJ. Analysis of fMRI time series revisited – again. Neuroimage 1995;2:173–81. Worsley KJ, Marrett S, Neelin P, Vandal AC, Friston KJ, Evans AC. A unified statistical approach for determining significant signals in images of cerebral activation. Hum Brain Mapp 1996;4:58–73. Worsley KJ, Cao J, Paus T, Petrides M, Evans AC. Applications of random field theory to functional connectivity. Hum Brain Mapp 1998;6:364–7. Worsley K, Liao C, Aston JAD, Petre V, Duncan G, Morales F, et al. A general statistical analysis for fMRI data. NeuroImage 2002;15:1–15. Yaglom A. Correlation Theory of stationary and related random functions I, basic results. New York, NY: Springer-Verlag; 1987. Ye J, Lazar NA, Li Y. Geostatistical analysis in clustering fMRI time series. Stat Med 2009;28:2490–508.
Contents lists available at ScienceDirect
Journal of Neuroscience Methods journal homepage: www.elsevier.com/locate/jneumeth
Computational Neuroscience
Nonparametric variogram modeling with hole effect structure in analyzing the spatial characteristics of fMRI data Jun Ye a,∗ , Nicole A. Lazar b , Yehua Li c a b c
Department of Statistics, University of Akron, Akron, OH, United States Department of Statistics, University of Georgia, Athens, GA, United States Department of Statistics & Statistical Laboratory, Iowa State University, Ames, IA, United States
h i g h l i g h t s • The new periodic variogram model can well describe the fluctuating spatial structure of fMRI data. • The hole effect model considers both the nonlinear physical relationship between the proximate voxels and the functional relationship between distant voxels.
a r t i c l e
i n f o
Article history: Received 17 August 2014 Received in revised form 10 November 2014 Accepted 11 November 2014 Available online 18 November 2014 Keywords: Bessel model Cross-validation Kriging Mean squared deviation ratio
a b s t r a c t When analyzing functional neuroimaging data, it is particularly important to consider the spatial structure of the brain. Some researchers have applied geostatistical methods in the analysis of functional magnetic resonance imaging (fMRI) data to enhance the detection of activated brain regions. In this paper, we propose a nonparametric variogram model for the complicated spatial characteristics of fMRI data. The new periodic variogram model can well describe the fluctuating spatial structure of fMRI data, considering both the nonlinear physical relationship between the proximate voxels and the functional relationship between distant voxels. We demonstrate the effectiveness of the new variogram model using fMRI data from a saccade study. © 2014 Elsevier B.V. All rights reserved.
1. Introduction The advent of modern neuroimaging techniques has revolutionized our understanding of the working of the human brain. The complicated dependence structures in the data obtained from these imaging methods provide challenges for analysis, but at the same time open a new window for the exploration of the spatial and temporal characteristics of activation patterns in the brain. In recent years there has been explosive growth in the number of neuroimaging studies using functional magnetic resonance imaging (fMRI). In fMRI experiments, researchers record the blood oxygenation level dependent (BOLD) contrast at locations across the subject’s brain to see which areas are involved in performing a given task or are responding to a particular stimulus. BOLD imaging takes advantage of inherent differences in the magnetic properties of oxygenated
∗ Corresponding author. Tel.: +1 3309728008. E-mail addresses: [email protected] (J. Ye), [email protected] (N.A. Lazar), [email protected] (Y. Li). http://dx.doi.org/10.1016/j.jneumeth.2014.11.008 0165-0270/© 2014 Elsevier B.V. All rights reserved.
and deoxygenated hemoglobin, hence fMRI has become a very powerful research tool for studying, albeit indirectly, neuronal activity. According to Kosslyn (1999), there are two questions addressed in neuroimaging: one is when the particular structures and processes are involved in the brain activations over time; the other is how the brain implements the information processing. Solid statistical analysis plays a key role in finding answers to these questions. A brain image is a matrix of numbers corresponding to spatial locations. We generally refer to each volume element in the image as a voxel, which is the three dimensional analogue to a pixel. Whereas structural MR images generally comprise a single three dimensional image, fMRI data are represented as a time series of such images. As a result, it is particularly important to consider the spatial and temporal correlation in the data and see the relationship of different areas in the brain. Geostatistics is a branch of applied statistics that focuses on providing quantitative descriptions of natural variables distributed in space or time. In geostatistics, researchers use variogram modeling to obtain predictions from data on spatially or temporally variable properties. Nowadays geostatistics is widely used in many fields of science and
102
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
industry where there is a need for evaluating spatially or temporally correlated data. Because of the spatial and temporal properties of the fMRI data, researchers have started to apply geostatistical methods to enhance the detection of activated regions of the brain. The variogram is a traditional geostatistical tool that provides quantification of the degree of directional spatial properties in the measured random variables (Rossi et al., 1992). The variogram gives an intuitive interpretation for the variations of the data in a secondorder stationary process (Schabenberger and Pierce, 2002); it only compares the average square difference between locations; and one need not to know the constant population mean (Ye et al., 2009). For these reasons, it is common to work with the variogram rather than the covariance function in geostatistical applications. Spence et al. (2007) use a Gaussian variogram model to find proximate voxels of interest in fMRI data, and show that such spatial analysis can identify regions of the brain exhibiting statistically significant group differences. They claim that close neighbors tend to have similar intensities and distant voxels show less similarity. They further show that neighboring voxels have strong nonlinear correlation so that the sample variogram values of the data are well fitted by the Gaussian variogram model. Their geostatistical analysis method offers important advantages over more traditional voxel-by-voxel approaches. However, their method only considers the close physical neighbors of the voxels. In fMRI data analysis, different regions in the brain may be functionally related, which means that the spatial correlation between voxels may not necessarily decrease with their increasing geometric distances. For example, the different areas involved in language processing will be functionally related, even though they may be physically separated. Bowman (2007) notes that distant voxels may exhibit correlation and proposes a novel way to describe the spatial structure of the fMRI data. He defines two different measures of relationship: one is the physical relationship, which is the geometric separations between two voxel locations and is the regular lag distance used in Spence et al.’s (2007) variogram model; the other is the functional relationship, which is the signal change between two voxels (i.e., the difference in the measured signals between two voxels) and is also called anatomical distance. As the variogram is defined as a function of the lag distance in geostatistics, Bowman uses the anatomical distance instead of the regular geometric distance as the lag distance in his new variogram model. He shows that the proposed variogram can find the functional relationship of the voxels. However, it is worth pointing out that this new defined lag distance in Bowman’s proposed variogram contains unknown measurement errors. A basic assumption in defining the variogram is to ignore an error in the lag distance and to assume the error only exists in the measurement of the signal (Stein, 1999). This means that the lag distance is assumed to be known without error in the measurement process. In fMRI data, the geometric distance between two voxel locations can be measured precisely. However, the signal values of the voxels always contain measurement errors and furthermore the errors may vary in different voxels. Hence it is very hard to well define the anatomical measure between two voxels. Although Bowman states that empirical values from unspecified auxiliary data can correct the anatomical measure, he does not explicitly describe how to deal with the measurement errors among the voxels and his approach might generate other additional problems, e.g., sometimes we do not know if the empirical values from the selected auxiliary data are well-suited for the current analysis. Our proposed model is also related to the random field model in fMRI data analysis. Random field theory is a branch of mathematics that has been extensively used in describing the spatial correlation of brain image data. In work by Worsley and colleagues (Worsley and Friston, 1995; Worsley et al., 1996, 1998, 2002; Friston et al., 2004), Gaussian or other monotonically homologous random field models describe the spatial correlation between neighboring
voxels. However, it is important to consider the spatial correlation among distant voxels too. In this paper, we will build an advanced periodic random field model that can well describe the fluctuant spatial structure of the random filed in the brain image data. Compared with the existing variogram modeling methods in fMRI data analysis (Bowman, 2007; Spence et al., 2007), the main contribution of the paper is that we propose a new nonparametric variogram model for the complicated spatial characteristics of fMRI data. Our hole effect model directly considers both the nonlinear physical relationship between proximate voxels and the functional relationship between distant voxels. Therefore it can effectively capture the main features of the functional relationship among different voxels. The rest of the paper is organized as follows: In Section 2, we introduce basic concepts and definitions of the new geostatistical method. Application of our method to a real fMRI data set is in Section 3. Section 4 discusses the proposed method and gives the conclusion. 2. Methods 2.1. Definition of the variogram Consider a slice of the fMR image at a specific time point; the signals are obtained over a regular grid of voxels, i.e. measured Z(s) : s ∈ D , where s = (x, y) denotes the coordinates of the sample site (x, y); Z(s) is a particular voxel value at location s, D is the set of the region of interest {s1 , · · · , sn }. For a general spatial model Z(s) = + (s), where is an unknown constant population mean and (s) is the zero mean random function at s, Z(s) is second order stationary if E{Z(si )} = E{(si + h)} = and Cov{Z(si + h), Z(si )} = C(h) for any si + h, si ∈ D, where lag h is a directional distance between two locations. Under second order stationarity, the variogram is defined as an alternative measure of spatial dependence (Schabenberger and Gotway, 2005), i.e. (h) = 12 Var{Z(s + h) − Z(s)}. Since (h) = C(0) − C(h), the variogram can reach a limiting maximum value called the sill, i.e., (h0 ) = C(0). Here lag h0 is known as the correlation range with C(h0 ) = 0, because it defines the average distance within which the random function remains correlated spatially. The variogram summarizes spatial relations in the data; it can be calculated directly from the data, in which case it is called the empirical or sample variogram. The computational formula for the empirical variogram (h) ˆ is (Isaaks and Srivastava, 1989). (h) ˆ =
1 2N(h)
{Z(si ) − Z(sj )}2 ,
(1)
(si ,sj )|si −sj =h
where s1 , . . ., sn in D are lying on an already known regular lattice, and the summation is over only the N(h) pairs of voxels which are separated by lag h. Since a relatively large N(h) is associated with a relatively small standard error (Morris, 1991; Rossi et al., 1992), we choose the number of data pairs to be at least 30 (Journel and Huijbregts, 1978). The assumption that the mean of the spatial model is constant is not realistic for fMRI data. Hence before the variogram calculations, we must first remove any trend surface (Schabenberger and Gotway, 2005). After trend removal, we consider Z(s) to be second order stationary (Cressie, 1993). In practice, we usually assume that the model only has a low-order trend, and the error term has the high-order structure. For example, we set up the model as Y (s) = f (s) + + (s) = f (s) + Z(s),
(2)
where f(s) is a trend function. After removing the trend f(s), we render the residuals stationary and use Z(s) to compute the variogram
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
103
for further analysis, and add back the trend f(s) at the end of the analysis. For a specific spatial structure, if the variogram function is direction independent, it is called isotropy (Chiles and Delfiner, 1999). We consider our detrended data as isotropy which is validated by rose diagram (Banerjee et al., 2005; Isaaks and Srivastava, 1989; Waller and Gotway, 2004) of the empirical variogram.
periodicity of variograms, highlighting that the variogram with hole effect structure can be used to describe the complicated spatial structure in the fMRI data too. In the next section, we build a nonparametric variogram model with the geometric lag distance to describe the complex relationship of the different voxels (Tobler, 1970; Gringarten and Deutsch, 2001).
2.2. Properties of the variogram
2.3. Variogram modeling
By the assumption of second order stationarity, the variogram defines the relationship between the variance and the lag distance. Hence variance is a function of lag alone. Understanding the characteristics of different variations of the voxel values at different locations of the brain is very important in fMRI data analysis. Based on the special properties of the fMRI data, we consider two types of relationship in the variogram, which have been described already in Spence et al. (2007) and Bowman (2007).
It is useful to model the empirical variogram because the results in a smooth fit provide estimators at arbitrary lags (Webster and Oliver, 2007; Schabenberger and Gotway, 2005). Current geostatistical practice in selecting a variogram model is often rather subjective, relying on empirical guidelines (Gorsich and Genton, 2000). By the special properties of the fMRI data, we examine two kinds of variogram model approaches: one is parametric, the other is nonparametric. The parametric approach, proposed by Spence et al. (2007), considers the physical relationship of the voxels only; our new nonparametric approach considers both the physical relationship and the functional relationship of the voxels.
2.2.1. Physical relationship and functional relationship In the mining and petroleum industries, as the magnitude of the lag separation vector increases, the variogram increases too, showing some linear properties. This is generally observed in nature and is called the physical relationship since regions close in space tend to have similar signal values. In brain image data, a similar relationship exists and is pointed out by Spence et al. (2007). They find that the spatial correlation between neighboring voxels is much stronger than that in many other data types, therefore the variogram approaches the origin with a somewhat reversed curvature, showing quadratic patterns (Spence et al., 2007). Hence the physical relationship may more properly be described as nonlinear. The finding from Spence et al. (2007) leads to an improved approach in analyzing the regions of interest in fMRI data. We will explore this property in more detail via a Gaussian variogram model. In fMRI data, not only the physical neighbors of the voxels are related, but also distant voxels with similar anatomy or function are related in some ways. Bowman (2007) observes this phenomenon in the fMRI data however Spence et al. (2007) ignore it. Bowman (2007) uses a monotonically increasing variogram model with an anatomical lag to describe the functional relationship of the voxels, where the regular lag distance is changed from geometric distance to anatomical distance. However, as we discussed in the Introduction, the anatomical lag always contains unknown measurement errors. 2.2.2. Hole effect relationship We consider the geostatistical method for model selection and model inference in the analysis of fMRI data. In some instances the variogram may seem to fluctuate more or less periodically, rather than increasing monotonically. That is, the variogram can decrease from its maximum to a local minimum and then increase again, which appears as a “down-hole”. This fluctuant structure of the variogram is called the hole effect (Journel and Huijbregts, 1978; Webster and Oliver, 2007). Hole effect structure of the variogram is usually ignored in geostatistics, because in the mining and petroleum industries from which most geostatistical analysis arose, the fluctuations in the empirical variogram are considered as more or less random noise (Webster and Oliver, 2007). But for some data (e.g., fMRI data), if we have already known that the fluctuations in the variogram are due to correlation among different distant regions, the hole effect structure is not artificial and should not be ignored. Furthermore, the hole effect structure can provide valuable information concerning functional variability of the data (Pyrcz and Deutsch, 2007). Several recent studies from ecological and environmental science (e.g., Radeloff et al., 2000; Chen, 2005) have revealed relations between periodic landscape patterns and
2.3.1. Parametric variogram modeling The parametric variogram is monotonically increasing. Wackernagel (2003) defines the Gaussian-type model in a variogram form c
(h) = 2 [1 − exp{−(ah) }],
(3)
where a is a parameter to determine the weights of the neighbors and 1 ≤ c ≤ 2. As c → 1, the model is more linear at small lags and the spatial correlation between nearby points is weaker. As c → 2, the model approaches the origin with a more quadratic shape and appears more sigmoidal, and the spatial correlation between nearby points is higher. The two extreme cases are: when c = 1, Eq. (3) becomes an exponential model; when c = 2, it becomes a Gaussian model. The pure Gaussian model is unstable, because the quadratic behavior at the origin in the variogram will generate extreme values at the borders of the estimated map in prediction (kriging) (Wackernagel, 2003). Hence, a nugget effect is usually added as a nested structure in the Gaussian model to yield a discontinuity at the origin; this avoids the extreme extrapolation properties and makes the kriging results more stable. In geostatistics, the name of nugget effect was coined by Matheron (1962) in reference to small nuggets of ore distributed throughout a larger body of rock. The nugget effect can be thought of as white noise, usually due to micro-scale variation or measurement error (Schabenberger and Gotway, 2005). It can be expressed in a variogram form as (h) = 2 for h > 0. Since the sum of valid covariance functions is itself a valid function, the Gaussian model used in the parametric approach becomes
2
(h) = 2 + ( 2 − 2 ) 1 − exp{−(ah) } .
(4)
2.3.2. Nonparametric variogram modeling Bessel model: We choose the Bessel model to describe the functional relationship because it can well describe the pronounced hole effect structure in the variogram function (Schabenberger and Gotway, 2005; Webster and Oliver, 2007). The covariance function of the Bessel model is defined as C(h) = 2
d 2
(bh)
−
J (bh),
(5)
where b is the number of sign changes, d is the dimension of the data, = (d/2) − 1 and J (·) is the Bessel function of the first kind of order (Yaglom, 1987; Cressie, 1993; Stein, 1999; Webster and
104
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
(a) Two sign changes (one period)
(b) Three sign changes (one period) 1.5 Variogram
Variogram
1.5
1
0.5
0
0
2 3 Lag (c) Four sign changes (two periods)
1
0.5
0
2
0
1
2 Lag
3
0
1
2 3 Lag (d) Five sign changes (two periods)
k=0
k!2
2
0
1
.
2 Lag
3
(6)
(h) = 2 1 − J0 (bh) ,
(7)
where 2 is the sill. Fig. 1 gives examples of Bessel variogram models (h) = 1 − J0 (bh) with lag h ∈ [0, 3.8] and parameter b = 2, 3, 4, 5, where the parameter b determines the number of sign changes in the model function. The number of sign changes (periods) monotonically increases as b increases. Nonparametric model approach: To increase the flexibility in modeling the variogram, we choose a linear combination of valid covariance functions. Treating the valid covariance functions as though they are basis functions yields a nonparametric model fitting approach (Ecker and Gelfand, 1997; Schabenberger and Gotway, 2005; Shapiro and Botha, 1991). The nonparametric variogram model is then defined as
2
2
2
(h) = + ( − )
1−
p
(9)
2
p−1
2
k=1
0.5
The Bessel variogram model used in our analysis is
,
1 1 c (h) = + ( − ) 1 − J0 (kbh) − exp{−(ah) } . (10) p1 p2
Oliver, 2007). When d = 1, C(h) = cos(bh); when d = 2, C(h) = J0 (bh); when d = 3, C(h) = (1/bh) sin(bh); when d→ ∞, C(h) = exp { − (bh)2 } (Stein, 1999; Schabenberger and Gotway, 2005). The periodic function has a weaker hole effect structure as d increases, and it becomes a Gaussian function when d→ ∞. Detailed proofs are in Stein (1999). We consider random field is of dimension two in what follows (i.e., d = 2). As J is the Bessel function of the first kind of order and = (d/2) − 1, the order is fixed as = (d/2) − 1 =0 in the data analysis (Webster and Oliver, 2007; Schabenberger and Gotway, 2005). Hence C(h) = J0 (bh), which can be expressed as ∞ 2k (−1)k bh
p
where 2 is the nugget effect, 2 is the sill, b is the number of sign changes of the basis function, p is the number of basis functions. Because of the nonlinear property of the fMRI data, we further replace one Bessel basis function with a Gaussian-type basis function to get more flexibility in the model fitting, i.e., 2
Fig. 1. Examples of Bessel variogram models (h) = 1 − J0 (bh). Panels (a)–(d) are the different variograms with parameter b = 2, 3, 4, 5 under lag h ∈ [0, 3.8]. The hole effect structure of the model changes with b.
J0 (bh) =
1 1− J0 (kbh) p
2
k=1
1
0
2
(h) = + ( − )
0.5
1.5 Variogram
Variogram
1.5
1
0
1
is also suggested that each basis function is equally weighted, i.e., ωk = 1/p, and k is chosen as k to simplify model fitting. Now the variogram model becomes
ωk J0 ( k bh)
,
(8)
k=1
where p is the number of basis functions, and ωk and k are weights. Theoretically the more basis functions we use, the better the fit (Cherry et al., 1996). However in practice it is too computationally intensive for large data sets to greatly increase the number of basis functions. Also, as the number of basis functions increases and the fit gets better, the model is less generalizable and less interpretable. When the Bessel function J0 (·) is considered as the basis function, Ecker and Gelfand (1997) suggest that p should be at most five. It
The present model is a hybrid between the Bessel basis functions and the Gaussian-type basis function. In the Gaussian-type basis function, the parameter a determines the weights of neighbors for the target voxel, and takes account of the physical relations in the voxels; in the Bessel basis functions, the parameter b determines the different sign changes of the target voxel, and hence considers the functional relations in the voxels; the parameters p1 and p2 satisfy 4/p1 + 1/p2 = 1, which allows for the different weights of the Gaussian-type model and Bessel model in this nonparametric fitting approach. Variogram fitting method: Methods for variogram fitting include maximum likelihood and least squares. Maximum likelihood estimation has an assumption of Gaussian distribution of the data and the fit has to be based on all data; by contrast, least squares does not have the distribution assumption and one can restrict the lag distance for the data. We choose least squares here because of its flexibility (Schabenberger and Gotway, 2005; Schabenberger and Pierce, 2002). Since fitting models in this way is a form of nonlinear regression, we use the Levenberg-Marquardt algorithm to obtain parameter estimates of the variogram (Marquardt, 1963). It is a challenging problem to well fit the empirical variogram because it is hard to match all its peaks and troughs by the variogram model. Manual fitting is simple and is done most commonly in geostatistics (Ma and Jones, 2001). However, it is preferable to use automatic fitting if we have to model many variograms at the same time. In our setting, we take the approach of “automatic fitting at first, and manual fitting afterwards”. We first fit the empirical variogram automatically by using plausible starting values and then check the value of the coefficient of multiple determination (R2 ) in the least squares method. If the R2 value is low (e.g. less than 0.8), we adjust the starting values of the initial parameters in the Levenberg–Marquardt algorithm and see if the fit is improved. If the fit is still not satisfactory, we manually fit the variogram by a procedure that embodies both visual inspection and statistical fitting, and follow these three rules (Webster and Oliver, 2007; Ma and Jones, 2001): (i) ignore the point-to-point fluctuation and concentrate on general trends; (ii) estimate the variogram to be accurate at short lags, with less accuracy at longer lags; (iii) match the cyclic pattern of the variogram at least to the first peak or trough. 2.4. Kriging To compare variograms we must work in the predictive space (Ecker and Gelfand, 1997). Kriging is a standard tool in geostatistics for spatial prediction. We choose the best variogram model by cross-validation in kriging. For each observation in the validation data set, the mean square error (MSE) could be calculated after kriging. However, MSE is insensitive to the variogram used for kriging, and kriging variance is sensitive to the variogram (Lark, 2000,
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
2009). Hence mean squared deviation ratio (MSDR) is more suitable to assess the performance of the different variogram models in estimating the functional or physical relations among the voxels (Ripley, 1981; Davis, 1987; Webster and Oliver, 2007). MSDR is the MSE divided by the kriging variance, hence it can be considered a weighted MSE. Dividing the error by MSE allows one to compare the magnitudes of both the actual and the predicted error in the cross-validation. If the model for the variogram is accurate then the MSE should be equivalent to the kriging variance; and so the MSDR should be 1. 2.4.1. General characteristics of kriging For predicting the value of the random variable Z(s0 ) at an unsampled site s0 , from the data Z(s1 ), . . ., Z(sn ) at sampled sites s1 , . . ., sn , kriging is given as the weighted sum of the values of the neighbors. The weights are calculated by minimizing the error variance of an assumed variogram model for the data with regard to ˆ 0 ) at the spatial distribution of the observed data points. When Z(s an unsampled site s0 is estimated from Z(s1 ), . . ., Z(sn ) at sampled sites s1 , . . ., sn , Z(si ) is decomposed as a geostatistical model Z(si ) = (si ) + e(si ), i = 0, . . ., n,
(11)
where (si ) is an unknown function, e(si ) is the zero mean secondorder stationary random function at spatial position si . n ˆ 0) = Ordinary kriging takes Z(s w Z(si ), subject i=1 i to
n
i=1
wi = 1.
Letting
2
n
m
be
the
Lagrange
parameter,
ˆ 0 )} − 2m{
− 1} is minimized by setting its E{Z(s0 ) − Z(s i=1 i partial derivatives with respect to wi and m equal to zero. Hence, (si , s0 ) =
n
wj (si , sj ) − m,
subject to
j=1
n
wi = 1.
(12)
i=1
The ordinary kriging variance is 2 OK =
n
wi (xi − x0 ) + m.
(13)
i=1
Here we first estimate the linear trend, then analyze the residuals of the trend by ordinary kriging, and finally add back the trend after finishing the analysis (Cressie, 1993). 2.4.2. Model selection by cross-validation We use K−fold cross-validation to choose between different variogram models based on their kriging results. In our setting, for voxel values Z(s1 ), . . ., Z(sN ) at locations s1 , . . ., sN , K−fold crossvalidation splits the N voxels into K folds, where the kth fold has k ) for k = 1, . . ., K. The voxel values in voxel values Z(s1k ), . . ., Z(sN/K each fold are removed in turn and new values are predicted at those locations by the voxel values in other folds using kriging. ˆ k ) be the kriging values predicted at locations ˆ k ), . . ., Z(s Let Z(s 1
N/K
k 2 (sk ), . . ., 2 (sk , and OK ) be the kriging variance. s1k , . . ., sN/K OK N/K 1
Then
2 1 k ˆ k ) / 2 (sk ). Z(si ) − Z(s OK i i N K N/K
MSDR =
(14)
k=1 i=1
When MSDR >1, the kriging variance underestimates the true estimation variance; when MSDR 8 s). We find that the mean values of MSDR in the four time
1.4
1.2
1
0.8
0.6
0.4
0.2
Gau−100
Gau−150
Gau−200
Three−basis
Four−basis
Fig. 8. MSDR values in the 28 time points under different approaches.
Five−basis
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
111
Fig. 9. Panel (a) is the 32 × 26 rectangular map and its 16 equally sized pieces. Panel (b) is a resampled map in bootstrap under Anti task. Panel (c) is a resampled map in bootstrap under Pro task.
Table 2 In the different time groups, MSDR mean values under the Gaussian-type model approach and the nonparametric model approach. Time points
Gau-100
Gau-150
Gau-200
Three-basis
Four-basis
Five-basis
61, 65, 69, 73, 77, 81, 85 62, 66, 70, 74, 78, 82, 86 63, 67, 71, 75, 79, 83, 87 64, 68, 72, 76, 80, 84, 88
0.4447 0.2909 0.3085 0.2959
0.4736 0.2914 0.3625 0.3774
0.7480 0.5017 0.6795 0.5881
0.5757 0.6133 0.6203 0.7149
0.9714 0.6772 0.7918 0.6506
0.9879 0.8783 0.8990 0.8162
Mean
0.3350
0.3762
0.6293
0.6311
0.7727
0.8954
groups are quite similar (Table 2), and the results are consistent with those from the 28 time points. That is, in every time group, the MSDRs are different under the following three comparisons: Gaussian-type versus nonparametric approach; Gau-200 model fitting versus Five-basis functions model fitting; and Gau-100 model fitting versus Gau-200 model fitting.
3.6. Further analysis in a bootstrap resampling method To further assess the differences of the models, we performed a bootstrap analysis. In this method, we cut every 32 × 26 rectangular map into 16 equally sized pieces (panel (a) in Fig. 9). Since there are 156 images, one for each time point, where 73 of them are under Anti task and 83 of them are under Pro task, we randomly combine the different pieces of the images according to the two different tasks, where there are 1400 resamples for each task (panels (b) and (c) in Fig. 9). The variogarm modeling and kriging are executed for the resampled images. We compare our Five-basis model with the Gau-200 model. The averages of MSDR and its standard error are 0.9374 and 0.0342 for the Five-basis model, and are 0.8168 and 0.0594 for the Gau-200 model (Table 3). The bootstrap results confirm the efficacy of our proposed method.
Table 3 In the bootstrap resampling method under Anti and Pro tasks, MSDR and its standard error (stder) with the Gau-200 model approach and the Five-basis model approach, where there are 1400 resamples for each task. Task
MSDR of Gau-200
Stder of Gau-200
MSDR of Five-basis
Stder of Five-basis
Pro Anti
0.7212 0.9124
0.0449 0.0739
0.9386 0.9361
0.0352 0.0331
Mean
0.8168
0.0594
0.9374
0.0342
4. Conclusions and future work In many geostatistical studies, the variogram is monotonically increasing; any hole effect structure is considered to be artificial and is ignored. However in fMRI data, high correlation exists between a given voxel and neighboring voxels as well as some relatively distant voxels. Hence the cyclical pattern of the variogram is very important for fMRI data analysis and it cannot be ignored. This is a major difference between geostatistical and fMRI data. In geostatistics, the physical relationship is most relevant since close regions are more related than distant regions, and distant regions will tend to be unrelated. But in fMRI, voxels may be functionally related even though they are not neighbors. This claim is similar to that in Bowman (2007), but our approach is totally different. Instead of using the anatomical measure in the variogram, we use the regular geometric measure to describe the functional relationship between voxels. The advantage of our method is that it does not need auxiliary data, and considers the original fMRI data directly. When we use the hole effect variogram model in model fitting, the presence of such a hole effect provides valuable information concerning complicated spatial variability of the fMRI data. In summary, the use of the hole effect structure in the variogram model provides a promising approach in the analysis of brain imaging data. It shows a good example of applying geostatistical ideas in a new area. However, most of the variogram models with hole effect structure are with dampening, which means that the functional relations are always weaker than the physical relations. But for a target voxel in fMRI data, its distant voxels may exhibit higher correlation than some closer neighbors (Bowman, 2007). How to well define a suitable hole effect variogram model and keep the balance between the physical relationship and the functional relationship in the model fitting will be considered in future work. For applying the method to the whole brain, conceptually, the procedure of the analysis would be the same as the analysis for a
112
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
single slice. Computationally, the extension to the whole brain is more of a challenge, since we quickly run into memory and storage issues. We will address efficient computation in future work.
Appendix A.
Acknowledgement
As suggested by one referee, we also try the fifth slice in the data analysis. After the data preprocessing steps (Fig. A.1), we consider the different variogram models for the empirical variograms at time points 61–88 (Fig. A.2). The best model will be chosen by
We wish to express our appreciation to Dr. Rebecca L. McNamee for kindly providing the saccade data set.
A.1. Data analysis in the fifth slice
Fig. A.1. Image maps at the fifth slice. Panel (a) is the mean map of all 156 time points. Panel (b) is the masked mean map. Panel (c) is the map at a randomly selected time point. Panel (d) is the mean subtracted map. Panel (e) is the linear trend map. Panel (f) is the trend removed map. The same processing procedure is carried out at every time point.
Fig. A.2. Panel (a) is the empirical variogram map in 0◦ , 90◦ and 270◦ directions at a randomly selected time point. Panel (b) is the empirical variogram maps in 0◦ , 90◦ and 180◦ directions at the same time point.
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
113
Fig. A.3. Maps for 5-fold cross-validation. Panels (a)–(e) are the five folds where the voxels denoted by white squares (144 voxels for (1), (2) and (3); 145 voxels for (4) and (5)) are removed in each fold for estimation. Panel (f) is the original map with 722 voxels.
Table A.1 MSDR for different time points (61–88) under the Gaussian-type model approach and the nonparametric model approach. Time point
Gau-100
Gau-150
Gau-200
Three-basis
Four-basis
Five-basis
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
0.5387 0.0631 0.4386 0.3762 0.7270 0.7289 0.3787 0.6260 0.3242 0.1525 0.6497 0.6355 0.3040 0.4838 1.1951 0.5643 0.6891 0.4665 0.5140 0.0758 0.2059 0.4237 0.7218 0.2859 0.8263 0.3955 0.0430 0.2303
0.7064 0.0456 0.5244 0.5389 0.8670 0.9909 0.4783 0.7154 0.3338 0.1110 0.7592 0.9643 0.3365 0.6369 1.1367 0.8192 0.6080 0.6127 0.6246 0.0539 0.2312 0.5833 0.0383 0.2891 1.2053 0.5584 0.0280 0.2840
0.5388 0.5097 1.0137 1.1245 1.0280 1.1017 0.6970 0.5122 0.3751 0.3099 1.0733 0.6384 0.4626 1.0047 0.6214 0.7779 1.1955 0.6061 0.9579 0.4168 0.2875 0.9995 0.2709 0.2477 0.1777 0.5315 0.2184 0.3682
0.6105 0.8059 0.2914 0.8072 0.5815 0.5850 0.8182 0.2191 0.8270 0.3281 1.0299 0.4657 0.7225 0.6977 0.8812 0.9769 0.5769 0.1933 0.2752 0.4512 0.6611 0.9451 0.4399 0.6937 0.4994 0.5505 0.9559 0.4822
0.4997 0.6191 0.3082 0.4753 0.4691 0.6505 0.7905 0.9759 1.2151 0.4163 0.7412 0.4701 0.7582 0.8377 0.6007 1.1788 0.9498 0.9141 1.2208 0.8534 0.6175 0.8721 0.4314 0.8536 0.7149 0.5987 1.0430 0.6214
0.6829 0.9862 1.0175 0.6117 0.6281 1.1862 1.1842 0.9334 0.6272 0.8674 0.8248 0.6855 0.6347 1.1562 1.1175 1.0726 1.1221 0.8923 0.6726 0.9266 0.8332 0.9964 0.6474 0.7938 0.8884 0.8881 0.8543 0.6215
Mean
0.4662
0.5386
0.6452
0.6204
0.7497
0.8697
Table A.2 In the different time groups, MSDR mean values under the Gaussian-type model approach and the nonparametric model approach. Time points
Gau-100
Gau-150
Gau-200
Three-basis
Four-basis
Five-basis
61, 65, 69, 73, 77, 81, 85 62, 66, 70, 74, 78, 82, 86 63, 67, 71, 75, 79, 83, 87 64, 68, 72, 76, 80, 84, 88
0.5150 0.3877 0.5630 0.3991
0.6216 0.5056 0.5128 0.5236
0.5807 0.7233 0.6932 0.5837
0.6398 0.5865 0.6702 0.5851
0.7463 0.7476 0.7292 0.7755
0.7738 0.9961 0.9026 0.8064
Mean
0.4662
0.5386
0.6452
0.6204
0.7497
0.8697
114
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115
1.2
1
0.8
0.6
0.4
0.2
0
Gau−100
Gau−150
Gau−200
Three−basis
Four−basis
Five−basis
Fig. A.4. MSDR values in the 28 time points under different approaches.
5-fold cross-validation in kriging (Fig. A.3). The MSDR results in Tables A.1 and A.2 and Fig. A.4 show that the variogram model with hole effect structure has better fit. The conclusion is consistent with results obtained using the fourth slice.
References Banerjee S, Carlin BP, Gelfand AE. Hierarchical modeling and analysis for spatial data. Boca Raton, FL: Chapman & Hall CRC Press; 2005. Bowman FD. Spatiotemporal models for region of interest analyses of functional neuroimaging data. J Am Stat Assoc 2007;102:442–53. Carroll MK, Cecchi GA, Rish I, Garg R, Rao AR. Prediction and interpretation of distributed neural activity with sparse models. NeuroImage 2009;44:112–22. Chen X. Statistical and geostatistical featuures of streambed hydraulic conductivities in the Platte River. Nebraska, Environ Geol 2005;48:693–701. Cherry S, Banfield J, Quimby WF. An evaluation of a nonparametric method of estimating semi-variograms of isotropic spatial processes. J Appl Stat 1996;23:435–49. Chiles JP, Delfiner P. Geostatistics: modeling spatial uncertainty. New York, NY: John Wiley & Sons, Inc; 1999. Cressie NAC. Statistics for spatial data. revised ed. New York, NY: John Wiley & Sons, Inc; 1993. Davis BM. Use and abuse of cross-validation in geostatistics. Math Geol 1987;19:241–8. Ecker MD, Gelfand AE. Bayesian variogram modeling for an isotropic spatial process. J Agric Biol Environ Stat 1997;4:347–69. Friston KJ, Frith CD, Dolan RJ, Price CJ, Zeki S, Ashburner JT, et al. Human brain function. 2nd ed. London UK: Elsevier; 2004. Gorsich DJ, Genton MG. Variogram model selection via nonparametric derivative estimation. Math Geol 2000;3:249–70. Goutte C, Hansen LK, Liptrot MG, Rostrup E. Feature space clustering for fMRI meta analysis. Hum Brain Mapp 2001;13:165–83. Goutte C, Toft P, Rostrup E, Nielsen FA, Hansen LK. On clustering fMRI time series. NeuroImage 1999;9:298–310. Gringarten E, Deutsch CV. Teacher’s aide variogram interpolation and modeling. Math Geol 2001;33:507–34. Guitton D, Buchtel HA, Douglas RM. Frontal lobe lesions in man cause difficulties in suppressing reflexive glances and generating goal-directed saccades. Exp Brain Res 1985;58:455–72. Haas TC. Lognormal and moving window methods of estimating acid deposition. J Am Stat Assoc 1990;85:950–63.
Hastie T, Tibshirani R, Friedman J. The elements of statistical learning: data mining, inference, and prediction. New York, NY: Springer Science; 2001. Huettel SC, Song AW, McCarthy G. Functional Magnetic Resonance Imaging. Sunderland, MA: Sinauer Associates Inc; 2004. Isaaks EH, Srivastava RM. An introduction to applied geostatistics. New York, NY: Oxford University Press, Inc; 1989. Journel A, Huijbregts ChJ. Mining Geostat. New York, NY: Academic Press, Inc; 1978. Kosslyn SM. If neuroimaging is the answer what is the question. Philos Trans R Soc Lond B: Biol Sci 1999;354:1283–94. Lark RM. A comparison of some robust estimators of the variogram for use in soil survey. Eur J Soil Sci 2000;51:137–57. Lark RM. Kriging a soil variable with a simple nonstationary variance model. J Agric Biol Environ Stat 2009;14:301–21. Lindquist MA. The statistical analysis of fMRI data. Stat Sci 2008;23:439–64. Ma YZ, Jones TA. Teacher’s aide modeling hole-effect variogram of lithologyindicator variables. Math Geol 2001;33:631–48. Marquardt DW. An algorithm for least-squares estimation of nonlinear parameters. J Soc Ind Appl Math 1963;11:431–41. Matheron G. Traite de Geostatistique Appliquee, Tome I, Memoires du Bureau de Recherches Geologiques et Minieres No. 14. Paris: Editions Technip; 1962. Morris MD. On counting the number of data pairs for semivariogram estimation. Math Geol 1991;23:929–43. Pham TD, Eisenblatter U, Baune BT, Berger K. Preprocessing film-copied MRI for studying morphological brain changes. J Neurosci Methods 2009;180: 352–62. Pyrcz MJ, Deutsch CV. The whole story on the hole effect; 2007 www. sciencedirect.com/science/article/pii/S016890020802010X Radeloff VC, Miller TF, He HS, Mladenoff DJ. Periodicity in spatial data and geostatistical models: autocorrelation between patches. Ecography 2000;23:81–91. Ripley BD. Spatial statistics. New York, NY: John Wiley & Sons, Inc; 1981. Rossi RE, Mulla DJ, Journel AG, Franz EH. Geostatistical tools for modeling and interpreting ecological spatial dependence. Ecol Monogr 1992;62:277–314. Schabenberger O, Gotway CA. Statistical methods for spatial data analysis. Boca Raton, FL: Chapman & Hall/ CRC Press; 2005. Schabenberger O, Pierce FJ. Contemporary statistical models for the plant and soil sciences. Boca Raton, FL: Chapman & Hall/CRC Press; 2002. Shapiro A, Botha JD. Variogram fitting with a general class of conditionally nonnegative definite functions. Computational statistics and data analysis 1991;11:87–96. Sommer FT, Wichert A. Exploratory analysis and data modeling in functional neuroimaging. Cambridge, MA: The MIT Press; 2003. Spence JS, Carmack PS, Gunst RF, Schucany WR, Woodward WA, Haley RW. The accounting for spatial dependence in the analysis of SPECT brain imaging data. J Am Stat Assoc 2007;102:464–73.
J. Ye et al. / Journal of Neuroscience Methods 240 (2015) 101–115 Stanberry L, Nandy R, Cordes D. Cluster analysis of fMRI data using dendrogram sharpening. Hum Brain Mapp 2003;20:201–19. Stein ML. Interpolation of spatial data, some theory for Kriging. New York, NY: Springer-Verlag; 1999. Tobler W. A computer movie simulating urban growth in the Detroit region. Econ Geogr 1970;46:234–40. Wackernagel H. Multivariate geostatistics, an introduction with applications. 3rd ed. New York, NY: Springer-Verlag; 2003. Waller LA, Gotway CA. Applied spatial statistics for public health data. Hoboken, NI: John Wiley & Sons, Inc; 2004. Wang K, Jiang T, Yu C, Tian L, Li J, Liu Y, et al. Spontaneous activity associated with primary visual cortex: a resting-state fMRI study. Cerebr Cortex 2008;18:697–704. Webster R, Oliver MA. Geostatistics for environmental scientists. 2nd ed. Chichester, England: John Wiley & Sons, Ltd; 2007.
115
Worsley KJ, Friston KJ. Analysis of fMRI time series revisited – again. Neuroimage 1995;2:173–81. Worsley KJ, Marrett S, Neelin P, Vandal AC, Friston KJ, Evans AC. A unified statistical approach for determining significant signals in images of cerebral activation. Hum Brain Mapp 1996;4:58–73. Worsley KJ, Cao J, Paus T, Petrides M, Evans AC. Applications of random field theory to functional connectivity. Hum Brain Mapp 1998;6:364–7. Worsley K, Liao C, Aston JAD, Petre V, Duncan G, Morales F, et al. A general statistical analysis for fMRI data. NeuroImage 2002;15:1–15. Yaglom A. Correlation Theory of stationary and related random functions I, basic results. New York, NY: Springer-Verlag; 1987. Ye J, Lazar NA, Li Y. Geostatistical analysis in clustering fMRI time series. Stat Med 2009;28:2490–508.
