Impact of sampling rate on statistical significance for single subject fMRI connectivity analysis.

A typical time series in functional magnetic resonance imaging (fMRI) exhibits autocorrelation, that is, the samples of the time series are dependent. In addition, temporal filtering, one of the crucial steps in preprocessing of functional magnetic resonance images, induces its own autocorrelation. While performing connectivity analysis in fMRI, the impact of the autocorrelation is largely ignored. Recently, autocorrelation has been addressed by variance correction approaches, which are sensitive to the sampling rate. In this article, we aim to investigate the impact of the sampling rate on the variance correction approaches. Toward this end, we first derived a generalized expression for the variance of the sample Pearson correlation coefficient (SPCC) in terms of the sampling rate and the filter cutoff frequency, in addition to the autocorrelation and cross-covariance functions of the time series. Through simulations, we illustrated the importance of the variance correction for a fixed sampling rate. Using the real resting state fMRI data sets, we demonstrated that the data sets with higher sampling rates were more prone to false positives, in agreement with the existing empirical reports. We further demonstrated with single subject results that for the data sets with higher sampling rates, the variance correction strategy restored the integrity of true connectivity.

INTRODUCTION

The voxel time series from functional magnetic resonance imaging (fMRI) is widely used to investigate the functional connectivity between various brain regions. The connectivity analysis in resting state fMRI studies widely employs sample Pearson correlation coefficient (SPCC) to measure the statistical association between two time series (Birn et al., 2013; Biswal, Yetkin, Haughton, & Hyde, 1995; Murphy, Birn, & Bandettini, 2013). When the time series are white (no autocorrelation), the variance of the SPCC, called nominal variance, is used for the inference on connectivity. However, by its very nature, the time series in fMRI exhibits sample‐to‐sample dependence or an inherent autocorrelation (Bollmann, Puckett, Cunnington, & Bartha, 2018; Fiecas, Cribben, Bahktiari, & Cummine, 2017; Friston, Jezzard, & Turner, 1994; Weisskoff et al., 1993). The autocorrelation affects the variance of the SPCC and hence the subsequent inference (Afyouni, Smith, & Nichols, 2018; Arbabshirani et al., 2014; Fiecas et al., 2017; Hart, Cribben, & Fiecas, 2018). The origin of the inherent autocorrelation is attributed to multiple physiological noise sources (Aguirre, Zarahn, & Esposito, 1997; Friston et al., 2000; Liu, 2016; Lund, Madsen, Sidaros, Luo, & Nicholsd, 2006; Purdon & Weisskoff, 1998; Zarahn, Aguirre, & Esposito1, 1997), repetition time (TR; Arbabshirani et al., 2014; Bollmann et al., 2018), and the low‐frequency fluctuations from the hemodynamic response function (Friston et al., 1995; Rajapakse, Kruggel, Maisog, von Cramon, & Yves, 1998).

The primary sources of the autocorrelation due to the physiological noise in the fMRI include pseudo‐periodic heartbeat (0.6–4.2 Hz; Cordes et al., 2001; Glover, Li, & Ress, 2000), subject respiration noise (0.1–0.3 Hz; Lund et al., 2006; Wise, Ide, Poulin, & Tracey, 2004), and equipment‐related noises such as the thermal drift in shims (less than 0.01 Hz; Purdon & Weisskoff, 1998; Weisskoff, 1996). The noise from subject initiated motion is usually impulsive in nature and may occur at any frequency (Friston, Williams, Howard, Frackowiak, & Turner, 1996; Patel et al., 2014; Power, Schlaggar, & Petersen, 2015). The signal of interest that is believed to contain the neural component in fMRI, usually spans the frequencies between 0.01 Hz and 0.1 Hz (Bharat, Edgar, & James, 1996; Cordes et al., 2001). Thus, the time series are usually denoised with a temporal band pass filter in order to reduce the level of autocorrelation and to provide a better signal‐to‐noise ratio. Owing to its importance, temporal band pass filtering is employed in various fMRI connectivity initiatives such as the Human Connectome Project (HCP; Essen et al., 2012).

A temporal filter, by its operational principle, is a linear system that convolves an input time series with the filter impulse response to produce a filtered time series (Hayes, 1996). While temporal filtering reduces the level of the physiological noise, it is well‐known that the shape of the filter introduces an autocorrelation artifact in the filtered time series (Bright, Tench, & Murphy, 2017; Davey, Graydena, Egan, & Johnston, 2013). This type of autocorrelation is called the filter‐induced autocorrelation as opposed to the inherent signal autocorrelation. Davey et al. (2013) showed that filtering increases the variance of the SPCC, thereby introduces false positives in the seed‐voxel‐based connectivity analysis. Based on this observation, they suggested a corrective measure constructed using the filter frequency response to reduce the fraction of the false positives. However, the corrective measure is derived without accounting for the inherent autocorrelation.

Recently, investigation of the impact of the inherent autocorrelation on connectivity has become increasingly important in the resting state fMRI community (Afyouni et al., 2018; Arbabshirani et al., 2014; Fiecas et al., 2017; Hart et al., 2018). In these studies, the approach is to compute the variance of the SPCC accounting for the inherent autocorrelation and use the variance for the reliable inference on connectivity. In this article, we refer to these approaches as variance correction. Arbabshirani et al. (2014) provided a heuristic variance correction for the time series that follow an autoregressive (AR) model of the first order and investigated its impact on brain networks. Using the classical result of Roy (1989) on the co‐variance structure of the autocorrelation, Fiecas et al. (2017) proposed a variance component model that accounts for the inherent autocorrelation and the heterogeneity across subjects on the variance of the SPCC. This approach was later extended by Hart et al. (2018) to investigate the impact of longitudinal autocorrelation in the functional connectivity networks of Alzheimer's patients. The unpublished work by Afyouni et al. (2018) used the approach in Roy (1989) to derive the effective degrees of freedom (DOF) of the SPCC and investigated the impact of the autocorrelation on the graph‐theoretic measures. We noted that the sampling rate was fixed and low (except in Afyouni et al. (2018)) in these variance correction approaches.

The impact of the sampling rate on modeling the inherent autocorrelation is studied in Corbin et al. (2018) for the general linear model. However, as far as we know, no such investigation has been undertaken to examine the impact of the sampling rate on variance correction. Recent simultaneous multi‐slice acquisition methods provide data sets with sampling rates greater than 1.5 Hz for the whole brain analysis. These methods have advantages such as separating the breathing and the cardiac pulsation from the low frequency fluctuations (Boyacioglu, Schulz, Koopmans, Barth, & Norris, 2015), and producing larger spatial extents of activation in the task (Demetriou et al., 2018; Todd et al., 2016) as well as in the resting‐state fMRI (Preibisch, Buhrer, & Riedl, 2015; Smith et al., 2013). Despite the advantages, data sets with high sampling rates exhibit decreased DOF and increased false positive rates (FPRs) in the fMRI connectivity analysis (Eklund, Andersson, Josephson, Johannesson, & Knutsson, 2012, Figures 1 and 2), (Bright et al., 2017, Figure 2). Thus, it is important to investigate the impact of the sampling rate on the variance correction.

Figure 1

The temporal filtering setup used in fMRI analysis. The time series v(n) and w(n) were processed with a filter of impulse response h(n) to obtain the filtered time series x(n) and y(n), respectively. It is desirable to find the statistics of the SPCC r between the filtered time series in terms of the h(n) and the nontrivial autocorrelation sequences r (k) and r (k) [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 2

The seed‐voxel‐based AR(1) network simulation with 10 nodes. The Node 1 acts a seed and the correlation between the seed and 9 other nodes are given in Table 2. The AR(1) model‐based signal autocorrelation is used for unfiltered time series. A practical band pass filter with TR = 1.4 s and the passband 0.009–0.1 Hz is considered [Color figure can be viewed at http://wileyonlinelibrary.com]

Table 2

The ground truth raw Pearson correlations (unfiltered and filtered), variances and the AR(1) parameters used to simulate the network in Figure 2

	Unfiltered		Filtered
Links	ρ vw	VS	ρ xy	VFS	α	β
1 ↔ 2	−0.004	0.000055	−0.0175	0.0027	−0.99	0.965
1 ↔ 3	0.51	0.0028	0.57	0.0091	−0.49	−0.79
1 ↔ 4	0.05	0.00032	0.35	0.0068	0.586	−0.99
1 ↔ 5	0.4	0.13	0.4	0.0263	−0.986	−0.986
1 ↔ 6	−0.1	0.02	−0.1	0.0175	0.89	0.89
1 ↔ 7	0.46	0.26	0.48	0.125	−0.99	−0.999
1 ↔ 8	−0.645	0.00087	−0.66	0.0031	0.12	0.42
1 ↔ 9	0.29	0.0023	0.3	0.0077	0.11	0.51
1 ↔ 10	0.63	0.1	0.63	0.08	−0.991	−0.991

In this article, we investigated the effect of the sampling rate on the variance correction for the resting state fMRI connectivity analysis. Using the classical result of Roy (1989), we derived a generalized analytical expression for the variance of the SPCC between the filtered time series whose unfiltered version exhibits an unknown autocorrelation structure. We expressed the generalized variance in terms of the filter impulse response in addition to the autocorrelation and the cross‐covariance of the time series. This framework enabled us to discern the filter‐induced autocorrelation from that of the inherent autocorrelation. By relating the filter impulse response to parameters such as the TR (inverse of the sampling rate) and the cutoff frequency of the filter, we showed that it is possible to investigate the impact of the sampling rate on the variance. We then present a variance correction strategy for the two commonly used significance tests for the SPCC, namely, the z‐test and the t‐test. With the correction strategy, we demonstrated that it is possible to mitigate the effect of false positives due to autocorrelation. In particular, we studied the resting state fMRI data set with three different sampling rates. We show that data sets with high sampling rates (short TR) produce a large number of FPR in the single subject results, which can be significantly mitigated by the variance correction. Since the autocorrelation may vary voxel‐wise, we discuss the importance of these corrections on the seed‐voxel‐based connectivity analysis.

MATERIALS AND METHODS

Background

Pearson correlation coefficient is widely used in the resting state fMRI connectivity analysis as it measures the normalized linear association between the time series v(n) and w(n). It is given bywhere γ (0) is the cross‐covariance between the time series at lag 0, σ , and σ are the respective standard deviations and | ρ | ≤ 1 is a constant. The time series are assumed to have zero mean. In practice, the SPCC (Fisher, 1915) is used to estimate the ρ . The SPCC between the two zero‐mean time series v(n) and w(n) is given bywhere, N is the length of the time series. Unlike ρ , the SPCC r is a random variable. For a white bivariate Gaussian time series, the asymptotic probability density function (PDF) of the r is normal with the mean ρ and the variancewhere, [x] denotes the variance of a random variable x. The variance of the SPCC plays a major role in the statistical inference on connectivity via hypothesis testing on the r .

Hypothesis testing

The two most important statistical tests for the SPCC widely employed in fMRI connectivity analysis are the z‐test and the t‐test (Sarty, 2006). Prior to the z‐test, the variance of the raw SPCC in Equation (3) is stabilized using Fisher's z‐transformation given by

The asymptotic variance of z(r ) is given by

z‐test: In the z‐test, the null hypothesis of zero correlation is tested against any of the possible alternative hypotheses by comparing the test statisticto a standard normal distribution. The statistic Z is called the Z‐score.

t‐test: In the t‐test, the null hypothesis of zero correlation is tested by comparing the statistic (t‐score)to a t‐distribution with the mean 0, the variance σ 2 = ν∕(ν − 2), and the DOF ν = N − 2. The variance of the SPCC, used in the above hypothesis tests, for the white time series is independent of the sampling rate. The variance for the autocorrelated time series depends on the sampling rate. Thus, it is important to study the effect of the sampling rate on the variance. We began the study by considering the effect of the inherent autocorrelation and the filter‐induced autocorrelation on variance.

Effect of the inherent autocorrelation and the filter parameters on variance of SPCC

Consider a temporal filter depicted as a linear time‐invariant system in Figure 1 with an impulse response h(n). H represents the complex frequency response and C (n) = h(n) * h(−n) the deterministic autocorrelation of the h(n). We assumed that the unfiltered time series v(n) and w(n) in Figure 1 were jointly stationary, that is, v(n) and w(n) were stationary and the cross‐covariance sequence, γ (n + k, n) = γ (k), depended only on the time difference k. Let r (k) and r (k) represent the respective nontrivial autocorrelation sequences of the unfiltered time series.

Let x(n) and y(n) denote the filtered time series which are related to the unfiltered time series v(n) and w(n) via the convolution sum:

Let r (k) and r (k) represent the autocorrelation sequences of the filtered time series. Our aim was to find the variance of the SPCC r between the filtered time series. Using the classical result from the multivariate time series analysis (Roy, 1989) and the linear system theory (Hayes, 1996), it may be verified that the PDF of the SPCC r in Figure 1, is asymptotically normal with the mean ρ and the variancewhere is a 3 × 1 vector given byand D is a 3 × 3 matrix

In the above equations, r (p) = r (p) * C (p), r (p) = r (p) * C (p), γ (p) = γ (p) * C (−p), and γ (p) = γ (−p). The steps leading to the variance in Equation (9) are given in Section S2. We noted that the variance characterizes the SPCC when the autocorrelation in the filtered time series is partly due to the inherent autocorrelation in the unfiltered time series and in part due to the autocorrelation induced by the filter through its impulse response h(n). We considered a practical band pass filter implemented in the Analysis of Functional NeuroImages (AFNI) Cox (1996) with the impulse responsewhere f is the high cutoff frequency (HCF) and f the low cutoff frequency (LCF) of the filter. Since the h(n) is a function of the parameters such as the TR and the f , the variance in Equation (9) is a function of these parameters as well. Thus, we related the variance of the SPCC to the sampling rate and the HCF. We also considered various other choice of band pass filters such as the Gaussian filters in FMRIB Software Library (FSL; Smith et al., 2004) and the finite impulse response (FIR) filters in MATLAB used in many neuroimaging analyses.

Special cases of the Equation (9) were obtained by setting the autocorrelation sequence of the unfiltered time series to a delta function (to study the effect of only the filter‐induced autocorrelation) or by setting the filter impulse response to a delta function (to study the effect of only the inherent autocorrelation). These two special cases are discussed below.

Effect of filter‐induced autocorrelation on variance

Approach in the literature

This approach is investigated in Davey et al. (2013). In this case, the time series v(n) and w(n) were white, that is, and (in Figure 1). Davey et al. (2013, eq. 8) have shown that in this case, the PDF of the SPCC r is asymptotically normal with the variance

We noted that the filter frequency response affects the variance and hence the subsequent inference on the SPCC. There are two shortcomings to this approach: (a) The nontrivial autocorrelation of the typical (unfiltered) fMRI time series is ignored and (b) The variance in Equation (13) is only an upper bound to the actual variance for this case shown (below) in Equation (14). (See Figure S1).

Approach using Equation (9)

Using the white time series assumption in Equation (9), that is, , we noted that , and γ (0) = ρ σ σ . In this case, the variance term converges to the variance (proof in Section S3):where is the normalized deterministic autocorrelation function of the impulse response. We note that the variance in Equation (14) is a function of the filter impulse response, whereas the variance in Equation (13) is a function of filter frequency response. We showed (in Figure S1) that the variance in Equation (13) is an upper bound to the actual variance in Equation (14). The variance in Equation (14) is well‐known in the area of large sample theory (Shumway & Stoffer, 2011, Theorem A.8). If we further assume h(n) = δ(n) that is, no temporal filtering then ρ (n) = δ(n), in which case the variance in Equation (14) converges to , which is the variance term in Equation (3).

Variance due to inherent signal autocorrelation alone

In this case, h(n) = δ(n) and hence , C (0) = 1. Therefore, the asymptotic variance of the SPCC in this case iswhere the vector and the matrix are the same as in the Equations 10, 11 except that the autocorrelation and cross‐covariance sequences were modified accordingly as follows: r (k) = r (k), r (k) = r (k), γ (k) = γ (k) and γ (k) = γ (k). The result in Equation (15) is a special case of the general result presented in (Roy, 1989, eq. 5).

Stochastic signal models such as the moving average (MA), the AR, and the autoregressive moving average (ARMA) that are widely used to explain the autocorrelation in the fMRI time series are categorized to this case. In Arbabshirani et al. (2014, eq. 16), an approximate variance of the SPCC is heuristically derived by assuming an AR(1) model for the inherent autocorrelation. They show that this variance increases for the different values of the model parameters. We showed that the variance derived in Arbabshirani et al. (2014) is an upper bound to the variance in Equation (15) (Figure S6). The large sample variance for the AR(1) was also derived in (Fiecas et al. (2017), Appendix) which is in agreement with the variance in Equation (15).

Corrections for the significance tests

The variance of the SPCC derived in the previous section may be used to correct the test statistics in Section 2.2 during inference. For the null hypothesis of zero correlation, the corrected z‐ and the t‐scores are listed in Table 1 for various autocorrelation scenarios.

Table 1

The variance uncorrected and the corrected Z‐ and the t‐scores

	Z‐score		t‐score
Autocorrelation type	Uncorrected	Corrected	Uncorrected	Corrected
Filter + signal induced (Section 2.3)	Z=zr_xy1N−3	Z=zr_xyc^TDcN	t=r_xyN−21−rxy2	t=r_xy11−rxy2c^TDcN
Only filter‐induced (Section 2.4.2)		Z=zr_xy1N∑nρh2n		t=r_xy11−rxy21N∑nρh2n
Only signal induced (Section 2.5)		Z=zr_xyc˜^TD˜c˜N		t=r_xy11−rxy2c˜^TD˜c˜N

The values of and D in Table 1 are in given in Equations 10, 11, respectively, and the ρ (n) is computed as in Equation (14). We note that in the last row of Table 1, r  = r (no filtering).

Variance estimation

The variance of the SPCC is a function of the autocorrelation and the cross‐covariance sequences and hence must be estimated from the time series data. In this article, we considered four types of variances:

Nominal variance, : The value of this variance is (Equation (6)) and it does not account for the autocorrelation of the time series.

Filter only variance, : This variance is computed using the filter impulse response as (Equation (14)) and it also does not account for the time series autocorrelation.

Signal only variance, : This variance is a function of the autocorrelation and the cross‐covariance between the unfiltered time series. With these functions estimated (explained below), the variance is computed as (Equation (15)).

Filter + Signal variance, : This variance is a function of the autocorrelation and the cross‐covariance between the filtered time series and is estimated as (Equation (9)).

We observed that the and the are scalar values and hence the same level of variance correction is applied to all the voxels. On the other hand, the and the are specific for a voxel and hence voxel‐wise correction is used. For brevity, in this article, we illustrate the variance correction for the Z‐scores, though it is equally applicable for the t‐scores as well. The autocorrelation and cross‐covariance sequences necessary to obtain the variances are estimated via their respective sample estimators (Shumway & Stoffer, 2011, eq. [1.34]). The sample estimator for the autocorrelation is with . Similarly, the sample cross‐covariance estimator is , and , where we assumed zero‐mean time series. Also, the value of the time series is zero for the indices other than 0, 1, ⋯, N – 1. We replaced the true correlations with the estimated correlations to determine the and the . The computational load required for estimating the over the whole brain is given in Section 4.

Empirical data

We used the statistical package R (R Core Team, 2017) for all analysis in this study.

Simulation 1. Inference on connectivity

Inference in the fMRI is largely based on the Z‐ or the t‐scores. Here we verified the accuracy of the corrected Z‐ and the t‐scores in Table 1 using the variances estimated from the time series. Toward this end, we generated 20,000 pairs of time series. The inherent autocorrelation in the time series is modeled using AR(1) with the parameter α. The value of α is constrained to be | α | < 1 in order to obtain a stationary model. A positive value of α implies signal autocorrelation with low pass filter characteristics and a negative value implies high pass filter characteristics (Hayes, 1996, fig. 3.12). In addition, as the |α| increases, the filter characteristics become sharper. The AR(1) model is widely used in the resting state fMRI literature to represent the autocorrelation of the time series as well as the error terms in the GLM (Arbabshirani et al., 2014; Bollmann et al., 2018; Corbin, Todd, Friston, & Callaghan, 2018). In this article, we used two AR(1) models with parameter value 0.2 and 0.8, respectively, to generate correlated bivariate Gaussian time series v(n) and w(n). The generation of the AR(1) time series from the white noise time series with specific variance is discussed in Section S4. The length of the time series was chosen to be N = 895 that corresponded to the length of one of the real fMRI data sets (TR = 0.645 s; see Section 2.9) evaluated in this article. We considered a realistic band pass filter with a passband 0.009 Hz to 0.1 Hz (Section 2.9). We recorded the raw SPCC between the 20,000 unfiltered and the filtered AR(1) time series. We converted the raw correlation values to Z‐ and t‐scores and variance corrected the scores as per Table 1. We computed the empirical PDFs of the variance corrected and the uncorrected scores (Figure 3).

Figure 3

The empirical PDF of the SPCC for the filtered and the unfiltered data. We set ρ  = 0 that corresponds to null hypothesis of zero correlation. The 20,000 pairs of bivariate Gaussian time series were generated each of length N = 895 that corresponds to one of the real fMRI time series data considered in this article with TR = 0.645 s. The AR(1) model with parameters 0.4 and 0.9 were used for the time series autocorrelation. A realistic band pass filter with the passband 0.009–0.1 Hz was considered. (a) Z‐score (b) t‐score [Color figure can be viewed at http://wileyonlinelibrary.com]

Simulation 2. Simulated brain network

To illustrate the importance of the variance estimation and correction strategy for the connectivity analysis, we simulated a seed‐voxel‐based AR(1) network of 10 nodes (Figure 2). The simulation of each link in the network is discussed in Section S4. Each node represents a voxel. The color of the node denotes the strength of the raw Pearson correlation coefficient between the filtered time series of that node and the filtered time series of the seed node (Node 1). The true Pearson correlation coefficient between the nodes is listed in Table 2. The table also lists the AR(1) parameters used to generate autocorrelation in the time series along with the ground truth variance, . The ranges for the and the were chosen to reflect the variances found in the real data set (Figure 6). We considered TR = 1.4 s, with N = 399; which lead to and (Figure 5a for the HCF = 0.1 Hz). After estimating the variances, we computed the corrected Z‐scores. We used the Benjamini‐Hochberg procedure (Benjamini & Hochberg, 1995; Benjamini & Yekutieli, 2001; Genovese, Lazar, & Nichols, 2002) to control the false discovery rate (FDR) for the corrected Z‐scores. We used a q value threshold of 0.001 throughout the article.

Figure 6

The histogram of the variance (a) (b) (band pass filter with passband 0.009–0.1 Hz) (c) for TR = 0.645 s and various HCF [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 5

(a) The plot of (Equation (14)) for various TR and the HCF of a fixed practical band pass filter from AFNI (Equation (12)). (b) The plot of comparing various filters. Without loss of generality, we set ρ  = 0 in Equation (14) [Color figure can be viewed at http://wileyonlinelibrary.com]

Experimental data and preprocessing

Participants

Resting state data sets were obtained from the Enhanced Nathan Kline Institute Rockland sample, Neuroimaging Data Releases Version 8 (Nooner et al., 2012). During the scans, participants were instructed to keep their eyes closed, relax their minds, and to not move.

MRI acquisition

Imaging was performed with a Siemens 3T TimTrio scanner. T1‐weighted images were acquired using a magnetization‐prepared rapid gradient echo (MPRAGE) sequence with TR = 2,500 ms (1 mm isotropic voxels). T1‐weighted images were subsequently used for spatial normalization. Functional images were acquired with a BOLD contrast‐sensitive gradient echo‐planar sequence. Three resting state fMRI sequences were included in the scan protocol with different TRs (NKI‐RS, 2017). Standard EPI sequence was used to acquire the data with TR = 2,500 ms (3 mm isotropic voxels, 5 min duration) and two other data sets were acquired using multi‐band sequences (Feinberg et al., 2010; Smith et al., 2012) with TR = 645 ms (3 mm isotropic voxels, 10 min duration) and TR = 1,400 ms (2 mm isotropic voxels, 10 min duration).

fMRI preprocessing

We used AFNI (Cox, 1996) to process the images. BEaST approach (Eskildsen et al., 2012) was used to skull strip the anatomical images and the skull stripped images were registered to a 3 mm MNI template (ICBM, 2009). White matter, gray matter, and cerebrospinal fluid were segmented from the images and respective masks were obtained. The functional image processing included the following steps: discarding of first five EPI volumes to allow for the signal to reach equilibrium, slice timing correction, rigid body correction of head movement with Friston 24 parameter model (Friston et al., 1996), coregistration of functional volumes to anatomical images and subsequent normalization of coregistered images to 3 mm MNI template. Nuisance‐signal regression was performed using the AFNI GLM with white matter, cerebrospinal fluid time series and 24 motion parameters as regressors. The nuisance regressed time series were band pass filtered. The LCF was set to 0.009 Hz, and we varied the HCF from 0.1 Hz to Hz in order to explore the effect of increasing the filter bandwidth on the FPR. The linear and the quadratic trends were removed before filtering and no spatial smoothening was performed. We chose 5 participants out of 114 whose frame‐wise displacement (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012) was below 0.4 mm for all the three TRs. The lengths of the time series after removing the first five volumes respectively were 895 (for TR = 0.645 s), 399 (for TR = 1.4 s) and 115 (for TR = 2.5 s).

We estimated the for the other choice of filters mentioned in Section 2.3. The Gaussian filter considered was a combination of the high pass filter (Gaussian‐weighted least squares straight line fitting) with the standard deviation σ = 1/(2 f ) seconds and the low pass filter with σ = 1/(2 f ) seconds (Smith et al., 2004). The FIR filter considered was the rectangular windowed to the inverse Fourier transform of the specified “brick wall” band pass filter. (MATLAB, 2018; signal developers, 2014, fir1 function).

Effect of the sampling rate

To study the effect of the sampling rate on single subject connectivity, we consider the seed‐voxel‐based connectivity analysis (Biswal et al., 1995; Mulders, van Eijndhoven, Schene, Beckmann, & Tendolkar, 2015) in the default mode network (DMN). This network consists of posterior cingulate cortex (PCC) with its positive correlated regions such as the medial prefrontal cortex (mPFC) with the inferior parietal lobule (IPL) being reported consistently (Uddin, Kelly, Biswal, Castellanos, & Milham, 2009, fig. 2b, 4‐th axial slice, Mulders et al., 2015). A peculiar negatively correlated region with the PCC is the anterior cingulate cortex (ACC; Uddin et al., 2009, fig. 2b, 3‐rd axial slice). The seed voxel was taken in the PCC whose MNI coordinates were (−2, − 54, 26; Uddin et al., 2009). The time series of the 19 voxels that are within the sphere of 5 mm radius around the seed voxel were averaged and the averaged time series was used as a seed voxel time series. We computed the raw SPCC between the seed voxel time series and each of the voxel time series covering the whole brain. Since the sampling rate is the inverse of the TR, we use the terms sampling rate and the TR interchangeably.

Effect on variance

We computed the variance for the three different TRs and the HCFs considering the practical band pass filter from AFNI. Since the variances and depend on the voxel pairs, we estimated them for each pair of voxels in the DMN. The set of variance values was then used to determine the histogram of the estimated variance. The values of the histogram were smoothed using a Gaussian kernel density estimator and the smoothed histograms were displayed for various TRs and the HCFs. We used the rule‐of‐thumb bandwidth suggested in Silverman (1986, eq. [3.31]) for the kernel smoothening.

Effect on FPR

The raw SPCC, computed in the DMN, were converted to correct Z‐scores via the estimated variances as per Table 1. The corrected scores were tested for the zero correlation (null hypothesis) against the standard Normal distribution and the p values were adjusted for multiple comparisons using the FDR method. The FPR was then computed as the fraction of voxels whose FDR corrected p values were less than the p‐value threshold (Patel & Bullmore, 2016, Section 2.8). We studied the FPR in terms of the TRs and various HCFs.

Effect on connectivity

We determined the seed‐voxel‐based connectivity in the DMN for the TR of 0.645 s and the HCF of 0.1 Hz. The estimated variances, , and in this network for each subject were used to obtain the corrected Z‐score maps after adjusting them for the FDR.

RESULTS

Empirical results

Simulation 1

The empirical PDF for the Z‐score (Figure 3a) and the t‐score (Figure 3b) are presented. Figure 3a indicates that, for the unfiltered case, the Z‐scores that are corrected with conform to the standard Normal distribution indicating that Z‐scores are properly corrected for the signal autocorrelation. On the other hand, the Z‐scores corrected using the no longer conform to the Normal distribution and hence the inference based on ‐corrected Z‐scores may lead to false positives. For example, the false positive significant correlation increases from 2.4% to 9% if uncorrected for the signal autocorrelation at 95% confidence interval.

Similarly, for the filtered case in Figure 3a, the Z‐scores are corrected with , and . For the filtered case, the Z‐scores corrected using conform to the standard Normal distribution. In this case, the false positive significant correlation increases from 2.4% to 21% if uncorrected for the filter and the signal autocorrelation. The increase in the false positive correlation clearly shows that filter‐induced autocorrelation is higher than that of the inherent autocorrelation in the fMRI time series, which is in agreement with the results and the discussions in Davey et al. (2013). The increase in the false positives when corrected using the (filter alone) compared to the (Filter+ Signal) is just 1.2% (from 2.4 to 3.6). This example shows that the gain obtained by exploiting the signal autocorrelation appears insignificant compared to that of using the filter‐induced autocorrelation alone. The signal autocorrelation, in this example, follows a simple AR(1) model. However, in the realistic fMRI data sets, the signal autocorrelation severely affects the false positives due to high sampling rates (Figure 7a,b).

Figure 7

The false positive rates of various variance correction methods as a function of the TR and the HCF. The FPR is calculated as the fraction of voxels whose q values are less than 0.001. The FPR comparison between (a) and (b) and , (c) and (d) at HCF = 0.1 Hz [Color figure can be viewed at http://wileyonlinelibrary.com]

Simulation 2

The Fisher z‐transformed unfiltered network is shown in Figure 4a. The Z‐score corrected networks, using the and the are shown in Figure 4b,c, respectively. We defined the significant nodes as the set of nodes that survive after correcting for the autocorrelation using the (in Figure 4c) or the (in Figure 4g). Those nodes that do not survive the correction were deemed insignificant. We observed from Figure 4c that only 3 out of 9 nodes survived when correcting with the compared to the nominal case where 6 out of 9 nodes are shown as significant. The node 10, which had a high correlation (0.63) with the seed voxel (Figure 2), was removed from the network due to the high variance (Table 2). This example elucidates the fact that the raw Pearson correlation value alone is not sufficient to show the connectivity strength of a link. The variance of the SPCC that reflects the autocorrelation must be incorporated in order to infer proper connectivity strength.

Figure 4

A 10‐node simulated seed‐based network. The Node 1 acts a seed and the raw Pearson correlation values of the unfiltered ρ and the filtered ρ time series are given in Table 2. (a) The unfiltered Fisher Z network shows the Fisher Z values, that is, z(ρ ; Equation (4)). The Z‐score values for the unfiltered case: (b) Nominal , (c) Signal only . (d) Filtered Fisher Z network, z(ρ ). The Z‐score values for the filtered case: (e) Nominal , (f) Filter only , and (g) filter + signal [Color figure can be viewed at http://wileyonlinelibrary.com]

Similarly, in Figure 4e,f,g, we show Z‐score networks corrected using the variances , and , respectively. Figure 4g indicates that the variance correction with removes those nodes just like the variance correction using the does in Figure 4c. The correction using the removes significant nodes, such as node 9, and includes the insignificant nodes {5,7,10}. This is because the filter‐only correction is blind to the inherent autocorrelation. This illustration clarifies the importance of variance correction for proper inference on connectivity. All the graph analyses in Figure 4 were performed using the igraph package in R (Csardi & Nepusz, 2006).

Experimental results

We report here the results reflecting the impact of sampling rate on the real fMRI resting state data sets. The results are presented in the following order:

The plot of the for the realistic filters as a function of the TR and the filter HCF

The smoothed histograms of the and the for various TRs and the HCF

The plot of the FPR as a function of the TR and the HCF

The axial slices of seed‐based DMN corrected Z‐score maps using , and for a fixed TR and an HCF.

Filter variance,

Figure 5a indicates the for the various HCF and shows that as the HCF increases toward Hz, the approaches the (shown as a triangle). For the widely used HCF of 0.1 Hz, we noted that the relative differences, , expressed in percentage for the TRs are 86%, 70%, and 43%, respectively. Thus, the variances of the short TR data sets are more severely affected (and hence generated more false positives) by filtering than the long TR data sets. This is consistent with the empirical findings from Bright et al. (2017), fig. 2), and Eklund et al. (2012), fig. 2) that shorter the TR, more severe the false positives. Figure 5b indicates the for the different band pass filters mentioned in Section 2.3. The filter in Equation (12) has a lower variance compared to the other filters for all the filter parameters.

Smoothed histograms of and

Figure 6a indicates the histogram of the corresponding to a representative participant for various TRs. The spreads out around the showing the range of variances due to autocorrelation. In addition, as the TR decreases (sampling rate increases), the of most of the voxel pairs appear above the . The relative increases in the mean of the compared to the for the TRs are 44%, 21%, and 3.4%, respectively indicating that short TR data sets are more prone to autocorrelation effects. Similar results were obtained for the other participants.

Figure 6b shows the histogram of the for various TRs. The passband of the band pass filter considered was 0.009 Hz to 0.1 Hz. Filtering increases the variance irrespective of the TRs. The relative increases in the mean of the histogram compared to the nominal variance for the TRs are 87%, 70%, and 34%, respectively, reflecting the impact of filter autocorrelation. Figure 6c shows the histogram of the (for TR = 0.645 s) for various HCFs of the filter. As the HCF increases, the variance decreases due to the reduced filter‐induced autocorrelation and hence the histograms move closer to the unfiltered case (Figure 6a). The bandwidths chosen for the kernel smoothening in Figure 6a–c were 0.00029, 0.00037, and 0.00021, respectively.

FPR

In Figure 7a, the FPR computed in the DMN using the (nominal) and the (filter only) is indicated as a function of the TR and HCF. The results are averaged over the participants. With the correction, the FPR is reduced uniformly for all HCF relative to the correction with the . In particular, for the HCF of 0.1 Hz, the relative reduction in the FPR compared to the for the TRs is 84%, 80%, and 67%, respectively (Figure 7d). As the HCF approaches Hz (no filtering), the FPR of the coincides with that of the . Figure 7b shows the FPR for the variance correction using the and the (filter + signal). As in Figure 7a, the FPR of the reduces uniformly for all HCF. For the HCF of 0.1 Hz, the relative reduction in the FPR compared to the for the TRs is 95%, 94%, and 88%, respectively (Figure 7d). As the HCF approaches Hz, unlike in Figure 7a, the FPR converges to the FPR of the (Figure 7c). This observation confirms that the inherent (unfiltered) signal autocorrelation can be exploited to accurately reduce the FPR. It appears that exploiting both the inherent and the filter‐induced autocorrelation reduces the FPR for all the TRs. Similar results were obtained for other resting state networks (data not shown). Spatial resolution may also affect the autocorrelation of an fMRI time series and hence the FPRs reported in Figure 7 may be altered if the spatial resolution changes.

Seed‐voxel based DMN connectivity

In Figure 8, the seed‐voxel‐based correlation maps in the DMN for a representative participant are shown. The Fisher Z correlation maps for the unfiltered (Figure 8a) and the filtered (Figure 8d) data show four clusters of Z‐scores: one cluster around the PCC, one around the mPFC, and two clusters around the IPL. We transformed the unfiltered Fisher Z map (Figure 8a) to a Z‐score map using the and the , as in Figure 8b,c, respectively. The ‐corrected Z map retains the DMN structure while removing the proportion of voxels that are not connected with the DMN. We make two observations in the corrected map. First, the extent of the DMN region is reduced and second, the values of the Z‐scores are reduced relative to the nominal case but still high enough to identify the DMN clusters reported in the literature. Similarly, the filtered Fisher Z map (Figure 8d) was transformed to Z‐score map using the , the , and the . The ‐corrected Z‐score map reduced the possible number of the false connective voxels that may not be the part of the DMN (Figure 8g). Similar results for the rest of the four participants are presented in Figures S2–S5.

Figure 8

The seed‐based DMN connectivity maps for the Participant 1 at TR = 0.645 s and HCF = 0.1 Hz. (a) The Fisher Z map for the unfiltered case. The Z‐score map for the unfiltered case corrected using (b) and (c) . (d) The Fisher Z map for the filtered case. The Z‐score map for the filtered case corrected using (e) , (f) , and (g) [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 9 indicates the Z‐score maps corrected using the and the at various sampling rates. The Z‐score maps in Figure 9a were obtained from different (prospective) measurements of the same participant. However, the maps in Figure 9b were obtained from retrospective down‐sampled measurements of the actual measurement at TR = 0.645 s.

Figure 9

Axial slices of the Z‐score maps (z = 33 mm) for the participant 1 at HCF = 0.1 Hz and for various TRs with and correction (a) prospective analysis (b) retrospective down‐sampling analysis [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 9a indicates that at TR = 2.5 s the correction excludes the ACC where as the correction retains the ACC though the extent of the ACC is reduced. This observation (at TR = 2.5 s) is for a single measurement of the participant and hence it is difficult to make a conclusion regarding the PCC–ACC network. Repeated measurements of the same participant or a group analysis are necessary to investigate the effect of correction on the networks. On the other hand, the retrospective analysis shows the ACC (Figure 9b) after the correction at all the TRs, though the extent of the regions is reduced compared to the correction.

DISCUSSIONS

The fMRI time series is well‐known for its inherent autocorrelation structure. Besides, the preprocessing steps such as temporal filtering induce additional autocorrelation. Exploring the ramifications of the autocorrelation in the resting state fMRI connectivity studies is of huge interest in the neuroimaging community. The impact of autocorrelation is reflected in the variance of the SPCC, which is used for the inference on the connectivity. This is investigated in Afyouni et al. (2018); Fiecas et al. (2017); Hart et al. (2018) for a fixed sampling rate (or TR) and a filter cutoff frequency. The purpose of the current work was to investigate the effect of the sampling rate and the filter cutoff frequency on the variance of the SPCC in addition to the autocorrelation and cross‐covariance functions of the time series (Equation (9)). The popular inference tests using the nominal variance are notably inflated at the high sampling rates. Thus, we suggest a variance correction for these tests (Table 1).

In our framework, we differentiate the impact of the inherent autocorrelation from that of the filter‐induced autocorrelation. We establish that the variance correction for the temporal filtering (without accounting for the inherent autocorrelation) proposed in (Davey et al., 2013, Equation (8)) is an overestimate (at all sampling rates) to the accurate variance presented by us (Section 2.4.2, Figure S1). In addition, we demonstrate that the variance correction accounting for only the filter‐induced autocorrelation is not sufficient.

How big is the impact of the sampling rate?

One way to quantify the impact of the sampling rate on the seed‐based connectivity analysis is by inspecting the percentage of the links removed after variance correction as a function of the TR. From Section 3.2, we note that variance correction in the DMN removes on an average 97%, 88%, and 75% of the false positive links relative to the current practice of not accounting for autocorrelation (at HCF = 0.1 Hz) at TRs of 0.645 s, 1.4 s, and 2.5 s, respectively. This removal manifested as a reduction in the extent of the DMN region while retaining its integrity. In addition, as the HCF increases, the percentage of links removed drops exponentially. At the HCF =  Hz (corresponding to no filtering), it stays at an average of 87%, 83%, and 57% for the respective TRs. This illustrates the huge impact of the sampling rate in the single subject analysis. One may argue that setting a high threshold on the nominally corrected Z‐score map (Figure 8e) may result in the reduced region DMN map just like the variance corrected map in Figure 8g. It is clear from the analysis of Figure 8e that even a high threshold value such as 10 may not lead to the same result as in Figure 8g. We remind the reader here that the significant impact of the autocorrelation discussed above is valid for the single subject results.

As illustrated in Figure 9, the conclusions regarding the impact of sampling rate on the variance corrected resting state networks using the single subject analysis can be misleading. Thus, group analysis is necessary to investigate the true impact on the networks.

The graph‐theoretic connectivity measures are well‐known for their sensitivity to the change in functional connectivity (Afyouni et al., 2018; van den Heuvel et al., 2017). For example, in Afyouni et al. (2018, Section 3.4), the measures such as the nodal strength and the local efficiencies are shown to experience a larger impact of the autocorrelation in the frontoparietal and the DMN nodes than the nodes from the subcortical regions of the HCP data with a fixed sampling rate (TR = 0.725 s, high pass filter).

The impact of the sampling rate in the group analysis has yet to be explored. Fiecas et al. (2017, Section 3.1) showed that the group difference between the skilled and the dyslexic readers exhibits no significant functional connectivity after correcting for the autocorrelation (TR = 1.97 s, band pass filter passband: [0.01, 0.10] Hz). Similarly, the variance corrected global functional connectivity network between Alzheimer's disease patients and the healthy control group is insignificant (Hart et al., 2018; TR = 2.3 s, low pass filter with cutoff frequency 0.1 Hz). The results reported in these works are for the low sampling rate (less than 0.51 Hz). Based on the observations from our single subject results and the graph theoretic results (Afyouni et al., 2018), we hypothesize that the effect of the high sampling rate (greater than 1.5 Hz) on the group analysis may also be significant.

The effect of sampling rate on the filter‐induced autocorrelation may be illustrated by defining correlation susceptibility as the ratio of the SPCC between filtered time series to that of the unfiltered time series. If the ratio is one, then the filtering does not affect the correlation values. The ratio quantifies the change in correlation values due to filtering. In addition, if the ratio is independent of the sampling rate, then correlation susceptibility should be same for all the sampling rates. In Figure 10, we show the ratio for the various sampling rates. The range of the correlation susceptibility values indicates the impact of sampling rate on the ratio. Data sets with lower sampling rates (TR = 2.5 s) have most of the susceptibility values closer to one than the data sets with high sampling rates.

Figure 10

The correlation susceptibility between the filtered and the unfiltered real fMRI data for the various TRs and for the band pass filter with passband 0.009–0.1 Hz [Color figure can be viewed at http://wileyonlinelibrary.com]

A note on whitening

An alternative approach to the variance correction is whitening. The autocorrelation sequence estimated from an fMRI time series can directly be used to whiten the time series before the connectivity analysis as in Bollmann et al. (2018); Bright et al. (2017); Christova, Lewis, Jerde, Lynch, and Georgopoulos (2011); Olszowy, Aston, Rua, and Williams (2018). Rather, we see a value in using the variance correction approach due to the following reasons:

The whitening is inverse filtering. The filter coefficients are obtained by fitting a model such as the ARMA to the time series by using the estimated autocorrelation. The invertibility of the model depends on the model coefficients which in turn depend on the estimated autocorrelation values. For a given estimated autocorrelation, it may be possible that the inverse (whitening) filter may not exist at all.

The whitening leaves behind a residual autocorrelation especially for the data sets with high sampling rates (Bollmann et al., 2018). The reasons for such a residual autocorrelation include mismatches in the model type, the filter‐order and the autocorrelation estimation. If the residual autocorrelation is not corrected, it may increase the false positives in the connectivity analysis according to Equation (9).

As discussed in Afyouni et al. (2018), whitening modifies the definition of the SPCC of the observed time series to that of the unobserved white time series.

Finally, it is crucial to observe that the method presented in this article is not confined only to the fMRI time series. It may be used for the other modalities such as the MEG and the EEG where the signals are acquired at high sampling rates.

Methodological limitations

The computation of the variance (in Equation (9)) needs a computational load in the seed‐voxel‐based connectivity analysis. The variance computation for a pair of voxels in the DMN (Figure 8) needs a mean CPU time of 13.12 ms in R package installed on a 64‐bit laptop computer with Intel Core i5‐7200U processor (2.50 GHz) and 4 GB RAM. The number of brain only voxels in the DMN for the 3 mm MNI registered images is 70,835. The mean CPU times required for a single subject whole brain analysis is therefore 15.49 min. This may be computationally burdensome for the group analysis.

The methodology proposed in this article is applicable for the pairwise SPCC and thus may not be suitable for the multi‐voxel connectivity analysis, such as those proposed in (Anders, Heussen, Sprenger, Haynes, & Ethofer, 2015). The estimators for the autocorrelation and the cross‐covariances in Section 2.7 behave well for a reasonably long time series. For short‐window time series such as those required in the dynamic functional connectivity, the estimators may not yield reliable auto correlation and cross‐covariance estimates. In such a case, short data record estimators such as those proposed in Afyouni et al. (2018) may be used.

Potential extensions

Dynamic network connectivity

It would be greatly beneficial to extend our approach to the temporal network dynamics domain. Sliding window methods that capture the dynamic, nonstationary network organization along time in resting state fMRI time series have been increasingly popular (Chang & Glover, 2010; Preti, Bolton, & Ville, 2017; Sakoglu et al., 2010). However, validating the results in the network dynamics is challenging. This is largely due to the fact that artifacts from the high sampling rates may manifest as significant events which may diminish or destroy the underlying true dynamics.

Autocorrelation due to other preprocessing steps

In this work, we have exploited the effect of the temporal filtering on the variance of the SPCC. The impact of the other preprocessing steps such as the spatial smoothening (Worsley, 2005) and the interpolation schemes (Power, Plitt, Kundu, Bandettini, & Martin, 2017) might be a useful extension.

Higher‐order autocorrelation

Higher‐order autocorrelation exploits higher‐order dependencies in the fMRI signal rather than the second‐order dependencies explored in this article. It would be interesting to investigate the impact of higher‐order dependencies on the variance of the SPCC.

CONCLUSIONS

The autocorrelation, a crucial feature of the time series in the resting state fMRI connectivity analysis, is sensitive to the sampling rate. The autocorrelation affects the inference on connectivity if not properly accounted for. In this study, we investigated the impact of the sampling rate on the functional connectivity analysis that accounts for the autocorrelation. Using the classical results from multivariate time series, we first showed that the effect of the sampling rate is reflected in the variance of the SPCC. Then with the simulation and the real data sets, we showed that the variance correction is an important step for proper inference of the connectivity analysis. Finally, we also showed that in the single subject analysis, the data sets with higher sampling rates are more prone to the false positive connectivity which may be reduced by the variance correction. Thus, our approach brings together the existing methods in the literature and generalizes the variance correction results for the arbitrary sampling rates.

R CODE AVAILABILITY AND REPRODUCIBILITY

A set of R codes, the necessary preprocessed data, and the other additional R data files to reproduce all the figures in this article has been uploaded in the GitHub repository, https://github.com/OliverMount/AutocorrfMRI.

Supplementary Material

Click here for additional data file.