Effective degrees of freedom of the Pearson's correlation coefficient under autocorrelation

Abstract

The dependence between pairs of time series is commonly quantified by Pearson's correlation. However, if the time series are themselves dependent (i.e. exhibit temporal autocorrelation), the effective degrees of freedom (EDF) are reduced, the standard error of the sample correlation coefficient is biased, and Fisher's transformation fails to stabilise the variance. Since fMRI time series are notoriously autocorrelated, the issue of biased standard errors - before or after Fisher's transformation - becomes vital in individual-level analysis of resting-state functional connectivity (rsFC) and must be addressed anytime a standardised Z-score is computed. We find that the severity of autocorrelation is highly dependent on spatial characteristics of brain regions, such as the size of regions of interest and the spatial location of those regions. We further show that the available EDF estimators make restrictive assumptions that are not supported by the data, resulting in biased rsFC inferences that lead to distorted topological descriptions of the connectome on the individual level. We propose a practical "xDF" method that accounts not only for distinct autocorrelation in each time series, but instantaneous and lagged cross-correlation. We find the xDF correction varies substantially over node pairs, indicating the limitations of global EDF corrections used previously. In addition to extensive synthetic and real data validations, we investigate the impact of this correction on rsFC measures in data from the Young Adult Human Connectome Project, showing that accounting for autocorrelation dramatically changes fundamental graph theoretical measures relative to no correction.

Keywords: Autocorrelation; Cross correlation; Functional connectivity; Graph theory; Pearson correlation coefficient; Quadratic covariance; Resting state; Serial correlation; Time-series; Toeplitz matrix; Variance; fMRI.

Fig. 1

Analysis of null resting state functional connectivity to illustrate the problem of inflated correlation coefficient significance. Panel A shows standardised BOLD data for the Left Dorsal Prefrontal Cortex (PFCd; 421 voxles) of HCP one subject (HCP-1). Panel B illustrates the standardised BOLD time series of R-Insula (red; 35 voxels) and L-SomMotCent (blue; 773 voxels), illustrating the dramatically different degree of autocorrelation. Panel C maps the Z-scores of correlation between this PCFd region and time series from a different HCP subject (HCP-2), parcellated with the Yeo's atlas and overlaid on an MNI standard volume. Panel D compares Z-scores accounting for autocorrelation vs. naive Z-scores, showing apparent significance (in this null data) with naive Z-scores and expected chance significance with xDF-adjusted Z-scores. On the horizontal axis are naive Z-scores that ignore autocorrelation, while on the vertical axis are Z-scores adjusted according to xDF. Uncorrected critical values ($\pm 1.96$) are plotted in dashed lines. Panel E shows autocorrelation of each time series (left). Horizontal solid lines indicate the confidence intervals calculated as described in Section 2.3. The difference in magnitude and form of autocorrelation among the three time series is evident, with PFCd exhibiting strong, long-range autocorrelation and R-Insula showing virtually no autocorrelation. Also shown is the cross-correlation (right panels) between HCP-1's PCFd and HCP-2's Left Central SomatoMotor Cortex (L-SomMotCent) (top), and HCP-1's PCFd and HCP-2's Right Insula (R-Insula) (bottom).

Fig. 2

Variation in strength of autocorrelation over space within an atlas, and between atlases. Panel A maps the autocorrelation index (ACI) voxelwise and for 3 different atlases, averaged over subjects (15 for voxel-wise, 100 for ROIs); variation is particularly evident for Yeo and Gordon; Power atlas is more homogeneous (but see Panel D). Panel B shows the impact of averaging within ROIs on autocorrelation. Left, shows ACI of individual voxels (blue dots) of a single subject across three regions of interests (ROIs) from the Yeo atlas. Right panel illustrates the ACI of ROI-averaged time series (blue dots) for 100 subjects, showing dramatic increase in ACI; red lines indicates the median. ROIs are Left Posterior Cingulate (LH-PCC), Left Somatosensory Motor (LH-SomMot) and Left Dorsal Prefrontal Cortex (LH-PFC). Panel C plots the ACI, averaged over subjects of the HCP 100 unrelated-subjects package, vs. region size for ACI time series and three atlases, where ICA and one atlas (MMP) are surface-based. There is a strong relationship between ACI and ROI size. The “ROI size” for ICA is defined as number of voxels in each component above an arbitrary threshold of 50. For MMP, the ROI size is defined as number of vertices comprising an ROI. Panel D considers the Power atlas, which has identically sized spherical ROIs, plotting ACI vs. distance to a voxel in the thalamus. Cortical ROIs have systematically larger ACI than subcortical ROIs. Panel E shows variance explained by inter-subject and inter-node ACI profiles for the Gordon, ICA200, Power and Yeo atlases; the large variance explained by inter-subject mean indicates substantial consistency in ACI over subjects.

Fig. 3

Evaluation of false positive rate control for testing $\rho =0$ with different autocorrelation correction methods. Panel A shows results using real data and inter-subject scrambling of HCP 100 unrelated subjects with the Yeo atlas ROIs, comprising 235,500 distinct Z-scores (see Fig. S6 for same results with other atlases). Left shows the QQ plot of Z-scores of each method, top right shows the $-{\text{log}}_{10}$ KS statistics (larger is better, more similar to Gaussian), and bottom right the FPR, all of which show that Naive and B35 have very poor performance. Panel B depicts a similar evaluation with simulated data, where a single ACF is used to simulate all time series with identical autocorrelation (see Section S3.5), again under the null; we additionally consider two “global” correction methods that assume common ACF between the nodes, G-Q47 and AR1MCPS. Here the Naive and the two global methods have poor false positive control. Panel C shows the FPR at the nominal 5% α level across five methods (columns) for identical (top row) and different (bottom row) ACFs, over a range of time series lengths. Naive (note different y-axis limits) and B35 have poor FPR control, while BH, Q47 and xDF all have good performance for long time series, with xDF having some inflation for the most severe autocorrelation structures with short time series. The setting of each simulation is coded by plotting symbol and colour, as shown at the bottom of the figure.

Fig. 4

Percentage bias of estimated standard deviation of $\hat{\rho }$ for different autocorrelation correction methods. Panel A plots the bias of the B35 method for $T=100$ (top) and $T=1200$ (bottom) for equal (left) and unequal (right) ACF's. Panel B plots the same for BH, and Panel C for Q47. Panel D plots the same information for a wider range of time series lengths T. These results show the dramatic standard error bias in BH35, BH and Q47 with increasing ρ. All results here are for our adaptive truncation method; see Figs. S8 and S9 for percent bias of different tapering methods. The setting of each simulation is coded by plotting symbol and colour, as shown at the bottom of the figure. Details of simulations and bias computation are found in Supplementary Materials; see Algorithm S2 and Eq. S(10). We exclude the results for biases of Naive standard error as they often exceed up to $\text{\%}60$ for autocorrelated time series; see Fig. S2.

Fig. 5

Performance of testing $\rho =0$ at level $\alpha =0.05$ on 5000 simulated correlation matrices ($114\times 114$, matching Yeo atlas) with 15% non-null edges (see Section S3.5). From top to bottom, specificity, sensitivity and accuracy (sum of detections at non-null edges and non-detections at null edges) are shown. Specificity (i.e. FPR control) is good for xDF, BH, Q47 and G-Q47, and sensitivity increases with time series length; accuracy is best for xDF, closely followed by BH and Q47.

Fig. 6

Impact of Naive, xDF and BH corrections on rsFC in one HCP subject parcellated with the Yeo atlas. Panel A plots rsFC Z-scores of xDF-corrected connectivity vs. Naive, showing that the significance of edges with Naive computation of $\mathbb{V}\left(\hat{\rho }\right)$ is almost always inflated, but to varying degrees. Solid lines are the critical values corresponding to the cost-efficient (CE) density. Dashed lines illustrates the critical values of FDR-corrected q-values. Taking xDF as reference, edges that are incorrectly detected with Naive are coloured green (FDR but not CE) and blue (CE). The black point marks edge (37,94) and red point (13,25), discussed in body text. Panel B plots rsFC Z-scores of xDF-corrected connectivity vs. BH, same conventions as Panel A, showing deflated significance of Z-scores computed with the BH method. The green point marks edge (103,104). Panel C pp-plot of p-values of Z-scores from xDF (green), BH (blue) and Naive (red) corrections. Dashed line is %5 Bonferroni threshold for 6441 edges. Panel D shows the differences in mean functional connectivity (mFC) of each correction method across statistical (FDR) and proportional (CE) thresholding. See Fig. S10 for a similar plot for a different HCP Subject.

Fig. 7

Overall changes in global and local graph theoretical measures with the 100 unrelated HCP package parcellated by Yeo atlas. Panel A, Bland Altman plots of xDF vs. Naive for weighted degree (top), betweenness (middle) and local efficiency (bottom) computed with a cost-efficient threshold. There is one point for each of 114 nodes, the particular measure averaged over subjects, and the nodes are colour coded according to their resting-state network assignment. Panel B Shows the same graph measures, but with statistical thresholding (corrected via FDR correction). Panel C shows the differences in weighted CE density (left) and Global efficiency (right), and Panel D illustrates the same results statistical FDR thresholding. There is a dramatic impact of correction method on all graph metrics considered. Similar results for Gordon (Fig. S12), Power (Fig. S13) and ICA200 (Fig. S14) is available in Supplementary Materials.

Fig. 8

Overall changes in global and local graph theoretical measures with the 100 unrelated HCP package parcellated by Yeo atlas. Panel A, Bland Altman plots of xDF vs. Naive for binary degree (top), betweenness (middle) and local efficiency (bottom) computed with a cost-efficient threshold. There is one point for each of 114 nodes, the particular measure averaged over subjects, and the nodes are colour coded according to their resting-state network assignment. Panel B Shows the same graph measures, but with statistical thresholding (corrected via FDR correction). Panel C shows the differences in weighted CE density (left) and Global efficiency (right), and Panel D illustrates the same results statistical FDR thresholding. There is a dramatic impact of correction method on all graph metrics considered. Similar results for Gordon (Fig. S15), Power (Fig. S16) and ICA200 (Fig. S17) is available in Supplementary Materials.