Insight and inference for DVARS

1Functional connectivity obtained with resting state functional magnetic resonance imaging (rs-fMRI) is typically computed by correlation coefficients between different brain regions, or with a multivariate decomposition like Independent Components Analysis. Both approaches can be corrupted by artifacts due to head motion or physiological effects, and much effort is dedicated to preprocessing rs-fMRI data and using diagnostic measure to identify bad scans.

2proposed and popularized a measure to characterize the quality of fMRI data, an image-wide summary that produces a time series that can detect problem scans. They called their measure DVARS, defined as the spatial standard deviation of successive difference images. In fact, DVARS can be linked to old statistical methods developed to estimate noise variance in the presence of drift (see Appendix A for DVARS history).

3While DVARS appears to perform well at the task of detecting bad scans — bad pairs of scans — it does not have any absolute units nor a reference null distribution from which to obtain p-values. In particular, the typical “good” values of DVARS varies over sites and protocols which makes it difficult to create comparable summaries of data quality across data sets. The emergence of the large scale data sets such as the Human Connectome Project (HCP) (1k subjects) and the UK Biobank (10k subjects) further motivates the need for automated, yet reliable, quantitative techniques to control and improve the data quality.

4The purpose of this work is to provide a formal description of DVARS as part of a sum of squares (SS) decomposition of the data, propose more interpretable standardized versions of DVARS, and compute DVARS p-values for a null hypothesis of homogeneity.

5The remainder of this work is organized as follows. We first describe the sum of squares decomposition for the 4D data and how this relates to traditional DVARS, and other new diagnostic measures it suggests. Then we describe a sampling distribution for DVARS under the null hypothesis, and mechanisms for estimating the parameters of this null distribution. We establish the validity and sensitivity of the DVARS test with simulations, and use two different fRMI cohorts to demonstrate how both the DVARS test and our‘DSE’ decomposition are useful to identify problem subjects and diagnose the source of artifacts within individual subjects.

6Here we state our results concisely relegating full derivations to Appendices.

7For T time-points and I voxels, let the original raw rs-fMRI data at voxel i and t be. Denote the mean at voxel i as, and by some type of overall mean value (i.e. a summary of the mean image, like median or mean). We take as our starting point for all calculations the centered and scaled data:

8The scaling ensures that typical brain values before centering are around 100 and are comparable across datasets.

9Denote the total (“all”) variability at scan t as

10We interchangeably use the terms variability and sum of squares, not to be confused with variance. Define two mean squared terms, one for fast (“differenced”) variabilitythe half difference between time t and at each voxel, squared and averaged over space, and one for slow variabilitythe average between t and at each voxel, squared and averaged over space.

11We then have the following decomposition of the average SS at time points t and,.for. This has a particularly intuitive graphical interpretation: If we plot and at, they sum to the midpoint between and found at (see Fig. 1).

12Sincewe now have a concrete interpretation for DVARS, with being the “fast” DSE component in the average SS at t and.

13This also leads to a decomposition of the total sum of squares over all scans: With averages A, D, S and E defined in Table S1 (row 1) we have the following “DSE” decomposition

14That is, the total variability (“A-var”) in the 4D dataset is the sum of terms attributable to fast (“D-var”), slow (“S-var”) and edge variability (“E-var”). D is also the average mean squared difference (MSSD; see Appendix A). Each term in the “DSE” decomposition can be split into global and non-global components, as shown in Table S1, rows 2–3 (as also noted by for ).

15Elements of the DSE decomposition can be visualized as time series (see Table 1) or as images. For example, just as a variability image with voxels is useful, we find that a D-var image, and a S-var image, offer more informative views of the noise structure.

16We see the arrangement of DSE values in Table S1 as a variant of an Analysis of Variance (ANOVA) table that summarizes contributions from fast, slow, end, global and non-global components to the total mean-squares in a 4D dataset. Traditionally ANOVA tables use sum-of-squares to partition variability, but we instead focus on root mean squared (RMS) or mean squared (MS) values to leverage intuition on typical noise standard deviation (or variance) of resting state fMRI data. To avoid any confusion, we call this variant ‘DSE table’.

17We calculate expected values for each of the DSE values for artifact-free data using different null models. In Appendix D we detail the most arbitrary version of this model, based only on time-constant spatial covariance,. Another model is based on time-space-separable correlation; this noise model assumes data with arbitrary spatial covariance but a common temporal autocorrelation for all voxels with a constant lag-1 autocorrelation ρ. While this is a less restrictive time series model than AR(1), in practice temporal autocorrelation varies widely over space, and we stress we only consider this as a working model to gain intuition on the DSE table. (Our null model for DVARS p-values, below, does not assume time-space separability). We also consider the idealized model of “perfect” data with completely independent and identically distributed (IID) 4D data.

18Table 2 shows three sets of reference values for the DSE table.1 The first pair of rows shows the expected value of the MS for each component for the separable model. This shows that all DSE components scale with the average voxel-wise variance, and as temporal autocorrelation ρ increases D-var shrinks and S-var grows. The global components are seen to depend on, the average of the elements of. This indicates, intuitively, that the greater the spatial structure in the data the more variability that is explained by the global.

19The next pair of rows in Table 2 show the expected MS values normalized to the expected A-var term. The A-var-normalized D-var and S-var diverge from 1/2 exactly depending on ρ, specifically. The global terms here depend on the ratio of average spatial covariance and average variance,.

20The final pair of rows shows expected values under the most restrictive case of IID noise. Here D-var and S-var are exactly equal, about, and we see that the global variability explained should be tiny,. This suggests that normalized global variability relative to the nominal IID value, i.e., an estimate of, can be used as a unitless index of the strength of spatial structure in the data. (This particular result doesn't depend on the separable model; see Appendix D).

21The handy result on the approximating ρ generalizes beyond the time-space-separable model: For an arbitrary model, both and normalized to A estimate a weighted average of the lag-1 temporal autocorrelations (see Appendix D.8). Hence, the convergence of D-var and S-var we observe as data is cleaned up has the specific interpretation of reduction in the average lag-1 autocorrelation.

22These reference models provide a means to provide DSE values in three useful forms. For each A-var, D-var, S-var and E-var term we present:RMS, the square root of the mean squared quantity,%A-var, a variability as a percentage of total mean-square A, andRelative IID, A-var-normalized values in ratio to nominal IID values.

23RMS, the square root of the mean squared quantity,

24%A-var, a variability as a percentage of total mean-square A, and

25Relative IID, A-var-normalized values in ratio to nominal IID values.

26For example, for A-var we have (1) RMS is, (2) %A-var is 100% and (3) relative IID is 1.0. For D-var, (1) RMS is, (2) %A-var is and (3) relative IID is

27For -var, (2) RMS is, (2) %A-var is and (3) relative IID isnoting that we normalize to A and not.

28We note that the fast and slow components can be defined as responses of linear time-invariant filters. The slow component corresponds to an integrator filter with power transfer function and the fast component corresponds to a differentiator filter with power transfer function, where ω is angular frequency and is the repetition time (TR). In other words, in time domain, can be interpreted as average of convolved BOLD signals with a rectangular window of [1 1] and with a [1 -1] window.

29We seek a significance test for the null hypothesiswhere is the mean under artifact-free conditions. Note this is equivalent to a null of homogeneity for or. If we further assume that the null data are normally distributed, we can create a test statisticapproximately following a distribution with degrees of freedom, where is the null variance (see Appendix E).

30What remains is finding estimates of and. The null mean of is the average differenced data variance,where is the variance of the differenced time series at voxel i. To avoid sensitivity to outliers, we robustly estimate each via the interquartile range (IQR) of the differenced data,where is the IQR of a standard normal, and is the inverse cumulative distribution function of the standard normal. Below we evaluate alternate estimates of, including the median of and directly as the median of.

31The variance of unfortunately depends on the full spatial covariance, and thus we're left to robustly estimating sample variance of directly. We consider several estimates based on IQR and evaluate each with simulations below. Since the IQR-to-standard deviation ratio depends on a normality assumption, and we consider various power transformations before IQR-based variance estimation (see Appendix F). We also consider a “half IQR” estimate of variancewhere hIQR is the difference between the median and first quartile, and. This provides additional robustness against contamination of the variance estimate from upward spikes.

32Finally, the values can be converted to p-values with reference to a distribution, and subsequently converted into equivalent Z scores,

33Note that for extremely large values of numerical underflow will result in p-values of zero; in such cases an approximate Z score can be obtained directly as.

34Under complete spatial independence the degrees of freedom will equal the number of voxels I, and so ν can be thought of an effective number of spatial elements; large scale structure will decrease ν while larger ν should be found with cleaner data. Though we caution that estimates of ν will be very sensitive to the particular estimators used for and.

35For intra-cohort investigation of corruptions, we propose that our D-var time series,, is a more interpretable variant of DVARS, as it represents a particular “fast” portion of noise variance, and when added to “slow” mean-square,, gives the total mean-square of the 4D data. However, these components are not suitable for inter-cohort comparisons, as the variance characteristics may vary with acquisition or scanner differences. In this section we propose a set of transformations which makes the inter-cohort comparison of the DSE components (including DVARS) possible.

36First consider the percent D-var explained at a single time point. Eqn. (5) could be used to find, in sums-of-squares units, the percent variability attributable to D-var at t,:However, problem scans can inflate and could mask spikes. Hence we instead propose to replace with its average A and compute percent D-var at time t asThis has the merit of being interpretable across datasets, regardless of total sum of squares. This is just percent normalization to A as discussed above.

37While %D-var can be more interpretable than unnormalized D-var, its overall mean is still influenced by the temporal autocorrelation. For example, if %D-var is overall around 30% and at one point there is a spike up to 50%, what is interesting is the 20 percentage point change, not 30% or 50% individually. Hence another useful alternative is change in percent D-var from baselineinterpretable as the excess fast variability as a percentage of average sum of squares. Later in Section 4.2.1, we show how is used as measure of “practical significance” to complement DVARS p-values.

38We previously have proposed scaling DVARS relative to its null mean,

39(While we had called this “Standardized DVARS”, a better label is “Relative DVARS.”) This gives a positive quantity that is near 1 for good scans and substantially larger than one for bad ones. However, there is no special interpretation “how large” as the units (multiples of ) are arbitrary; as noted above, DVARS falls with increased temporal correlation, making the comparison of these values between datasets difficult.

40Finally the Z-score or may be useful summaries of evidence for anomalies.

41To validate our null distribution and p-values for we simulate 4D data as completely independent 4D normally distributed noisefor drawn uniformly between and for each i,.

42We manipulate two aspects in our simulations, time series length and heterogeneity of variance over voxels. We consider T of 100, 200, 600 and 1200 data-points, reflecting typical lengths as well as those in the Human Connectome Project. We use three variance scenarios, homogeneous with, low heterogeneity and, and high heterogeneity and.

43We consider four estimates of. First is the very non-robust sample mean of, denoted, considered only for comparative purposes. Next we compute the mean and median of (Eqn. (13)), the robust IQR-based estimates of differenced data variance at each voxel. Finally we also consider the empirical median of,. For, all estimates were based directly on; for comparative purposes we considered the (non-robust) sample variance of,, and IQR-based and hIQR-based estimates of variance with power transformations d of 1, 1/2, 1/3 and 1/4, denoted generically; note is theoretically optimal for (see Appendix F).

44For p-value evaluations, we only evaluate the most promising null moment estimators, and for, and with, and,. We measure the bias our estimators in percentage terms, as and, where the true value are and (as per Appendix E when ). For each method we obtain P-values and create log P-P plots (probability-probability plots) and histograms of equivalent Z-scores.

45Similar simulation settings are used to evaluate the power of the DVARS hypothesis test, except we consider 4 different autocorrelation levels. This range is chosen to reflect observed estimates of lag-1 autocorrelation coefficients in the HCP cohort. Inferences are assessed in terms of sensitivity and specificity.

46All simulations use 1000 realisations.

47We evaluate the impact of the DVARS test as a tool for “scrubbing” (scan deletion) on functional connectivity (FC) measued with Pearson's correlation coefficient. We consider FC between all possible pairs of Region of Interests (ROI) in each subject for a given ROI atlas. The mean time series of each ROI is obtained by averaging all the time series within a ROI. To parcellate the brain, we use two data-driven atlases; Power Atlas which is constructed of 264 non-neighboring cortical and sub-cortical ROIs and each ROIs has 81 voxels (is case of 2 mm isotropic volumes) and Gordon Atlas which is constructed of 333 cortical regions of interests with different sizes.

48We use two popular methods to evaluate the effect of the DVARS inference on functional connectivity. First, we use the QC-FC analysis which begins by creating per-edge, intersubject scores, the correlation of the number of removed volumes and FC; these scores are plotted against the inter-ROI distance (in mm). We then use LOESS smoothing method (with span window of %1) to summarize the association for each method. For further details about QC-FC method, see,,. We use QC-FC to compare our DVARS hypothesis test to four other scan scrubbing methods. From we use two FD thresholds, lenient (0.2 mm) and conservative (0.5 mm), and a DVARS threshold of 5. From FSL's fsl_motion_outliers tool, we use a DVARS threshold corresponding to box-plot right-outliers, 1.5 IQRs above the 75%ile. Note that the first three approaches used a fixed threshold, while the FSL approach gives a run-specific threshold.

49The objective of this FC analysis is to investigate whether DVARS inference test performs as well as the available thresholding methods (such as arbitrary thresholding of FD and DVARS) and if so, whether it delivers the optimal results while sacrificing the fewest temporal degree of freedom as possible. Therefore, we only present the results for the Minimally pre-processed data sets.

50We use two publicly available data-sets to demonstrate the results of methods proposed in this paper on real-data. First, we use 100 subjects from ”100 Unrelated” package in the Human Connectome Project (HCP, S1200 release). We chose this dataset due to the high quality and long sessions of the data. Second, we used first 25 healthy subjects from the New York University (NYU) cohort of the Autism Brain Imaging Data Exchange (ABIDE) consortium via Preprocessed Connectome Project (PCP). We selected this cohort for its high signal-to-noise ratio and the more typical (shorter) time series length.

51For full details see; in brief, 15 min eyes-open resting acquisitions were taken on a Siemens customized Connectome 3T scanner with a gradient-echo EPI sequence, TR = 720ms, TE = 33.1 ms, flip angle = 52° and 2 mm3 isotropic voxels. For each subject, we used the first session, left to right phase encoding direction (See Table S2 for full details of subjects). We considered each subject's data in three states of pre-processing: unprocessed, minimally pre-processed and ICA-FIXed processed. Unprocessed refers to the raw data as acquired from the machine without any pre-processing step performed, useful as a reference to see how the DSE components change with preprocessing steps. Minimally pre-processed (MPP) data have undergone a range of conventional pre-processing steps such as correction of gradient-nonlinearity-induced distortion, realignment aiming to correct the head movements, registration of the scans to the structural (T1w) images, modest (2000s) high pass filtering and finally transformation of the images to the MNI standard space.

52Finally, after regressing out the 24-motion parameters, an ICA-based clean up algorithm called ICA-FIX is applied, where artifactual ICA components, such as movement, physiological noises of the heart beat and respiration, are regressed out non-aggressively. Due to extent of the FIX denoising and an ongoing debate regarding the nature of the global signal, we did not consider global signal regression with the HCP data. From now on, we call this stage ’fully pre-processed (FPP)’ to be consistent with the ABIDE-NYU cohort we describe in the following.

53We use use 20 healthy subjects of New York University (NYU) data-set. For full details visit Pre-processed Connectome Project website http://preprocessed-connectomes-project.org/; in brief, 6 min eyes-closed resting acquisitions were taken on an Allegra 3T scanner with a gradient echo EPI sequence, TR = 2000ms, TE = 15ms, flip angle = 90°, and 3 mm isotropic voxels (See Table S2 for full details of subjects). In this study, each subject was analyzed using Configurable Pipeline for the Analysis of Connectomes (C-PAC) pipeline, in three stages; unprocessed, minimally pre-processed and fully pre-processed. The unprocessed data are raw except for brain extraction with FSL's BET. Minimally pre-processed data were only corrected for slice timing, motion by realignment and then the data were transformed into a template with 3 mm3 isotropic voxels. Fully pre-processed data additionally had residualisation with respect to 24-motion-parameters, signals from white matter (WM) and cerebrospinal fluid (CSF), and linear and quadratic low-frequency drifts. Conventionally this pipeline deletes the first three volumes to account for T1 equilibration effects, but we examine the impact of omitting this step for the raw data.

54Further, we also use all healthy subject of ABIDE (530 subjects) to show how DSE decomposition can be used to compare the data-sets, cohorts and pipelines.

55Fig. 2 shows the percentage bias for the null expected value (left panel) and variance (right panel) for different levels of variance heterogeneity and time series length.

56The direct estimates of the based on the time series perform best on this clean, artifact-free data, while estimated on variance of the differenced data ( and ) degrades with increasing heterogeneity. The estimates of variance have relatively less bias but it is difficult to identify one particular best method, save for IQR often (but not always) having less bias than hIQR, and lower d generally associated with less bias.

57On balance, given the generally equivocal results and concerns about robustness, for further consideration we focus on (median of ) and (median of ) as promising candidates for, and hIQR with and hIQR with for.

58Fig. 3 shows log P-P plots for p-values and histograms of approximate Z scores,; values above the identity in the P-P plot correspond to valid behavior. While all methods have good performance under homogeneous data, (panels A & C) is not robust to variance heterogeneity and results in inflated significance. In contrast, (panels B & D) has good performance over all, for variance estimated with either or (top and bottom panels, respectively), and also yields good approximate Z-scores. On the basis of these results, we elected to use as the only reliable option for the mean, and hIQR, as a variance estimate, and use these settings going forward.

59Fig. 4 shows the results of the power simulation. For all sample sizes and autocorrelation parameters, and for the 1% and 10% artifact rates, power was always above 80% and often 100%. Increased autocorrelation resulted in improvements in power, while higher artifact rates reduced power. For the 20% artifact rate power was adequate ( 80%), but falls to zero for the 30% artifact rate. These results suggest that, at the highest spike rate, the artifacts start to be become indistinguishable from the overall noise (see Fig. S2 for one realization). However, the distribution of DVARS values (Fig. S1) suggest that the constituent null and artifact components are distinguishable even at the highest spike rate, but would require yet more robust methods for estimating the null component than we have employed.

60We first focus on selected results of two HCP subjects, then later summarize results for all HCP and ABIDE subjects.

61Fig. 5 shows different standardized DVARS measures, as introduced in section 2.5, as well as the other DSE components for subject 118730 of the HCP cohort (See Figs. S4, S5 and S6 for more results.). The first six plots corresponds to the variants listed in Table 3; the bottom two plots show “DSE plots,” plots of,, and components, upper plot with minimal pre-processing, lower with full pre-processing. The gray and magenta stripes indicate 19 data points identified as having significant DVARS after Bonferroni correction, with magenta indicating time-points that are additionally practically significant by the criterion. In Fig. 5, the largest occurs at index 7 (i.e. 7th and 8th data points) and has = 4.07, large in terms of being -var = 70.16% of average variability, Z = 36.33 indicating extreme evidence for a spike, and having %D-var = 41.20% more sum-of-squares variability than expected. The least significant occurs at index 726, with = 2.83; while its Z = 4.36 is not a small Z-score, with just %D-var = 4.95% excess variation, it is a relatively modest disturbance. In contrast, we find that the values of original DVARS or relative DVARS do not offer a meaningful interpretation. Table S4 shows values for all significant scans.

62The bottom panel of Fig. 5 shows the DSE plot for fully pre-processed data. This data now exhibits the idealized behavior of IID data, with D-var and S-var components converging at 50% of average variability (see right-hand y-axis). However, interestingly, the change is not similar for all DSE components. Note how is around 2.6 before clean up, and 2.5 after clean up, while falls dramatically with cleaning, indicating that nuisance variance removed was largely of a “slow” variability component. Also observe that cleaned up results drops in total variability component, where artifacts were found, indicating variance is removed by FIX.

63Finally, Table 4 explores the use of the estimated degrees of freedom ν as an index of spatial effective degrees of freedom. Raw data, exhibiting substantial spatial structure, has, which increases to for fully preprocessed data, still only about 5% of the actual number of voxels.

64FC evaluations based on 55,278 unique edges from the Gordon atlas are shown in Fig. 6 (see Fig. S9 for Power Atlas results). Panel A shows the QC-FC analysis of five thresholding methods, compared to unscrubbed QC-FC. The results from the DVARS test appear comparable to the other methods, but Panel B of Fig. 6 show that the DVARS test removes many fewer scans on average, preserving temporal degrees of freedom. A related evaluation, comparing DVARS hypothesis test scrubbing to random scrubbing, finds that FC is significant impacted by the DVARS scrubbing (Figs. S10 and S11).

65We note that the sole purpose of preceeding QC-FC analysis is to ensure that the DVARS inference test outperforms other arbitrary thresholds available in literature, and therefore we do not show the similar results for FPP data.

66Fig. 7, Fig. 8 shows the minimally and fully pre-processed DSE decompositions, respectively, of HCP subject 115320.

67Fig. 7, upper panel, shows that if the strict FD threshold, 0.2 mm, were used 47% of scans would be flagged, while the lenient threshold, 0.5 mm, appears to miss several important events. For example, around scans 775 and 875 there are two surges in, rising to about 60% and 40% average sum-of-squares (excesses of 30% and 10%, respectively, from a baseline of about 30%) while FD remains low. The lower panel's pie chart shows that S-var explains just under 75% of total, and almost all of global sum-of-squares. The Edge component is also 1.5 above its expectation.

68In Fig. 8, the fully preprocessed data-set shows roughly equal fast and slow components, as reflected in the overlapping and sum-of-squares time series (blue and yellow, respectively) and the pie and bar charts for total sum-of-squares. Edge component E-var has also dropped to fall in line with IID expectations. However, this convergence is not homogeneous over scans and excursions of S-var are still found after scan 650. However, these are much reduced relative to MPP data (no more than 75% of average sum-of-squares, compared to over 150% in Fig. 7).

69Note that while significant DVARS are found in the FPP data, they are small in magnitude: Table 5 lists the 10 significant tests, none with %D-var greater than 6%. If we used a %D-var of 5% we would still mark 4 of these 10 significant; while we might hope for better performance from the FIX method, note the severe problems detected towards the end of the scan (Fig. 7).

70The smallest significant %D-var is 2.66%, which is smaller than the least significant scan detected in the minimally preprocessed data, 3.78%. This indicates the increased sensitivity in our procedure as the background noise in the data is reduced. Note that the majority of the spikes detected in Fig. 7 has been removed by ICA-FIX (Fig. 8), however the algorithm has left down-spikes which could be detected via a two-sided version of the test explained in section 2.4.

71Temporal diagnostics of before and after clean-up for three other subjects (HCP subject 118730, NYU-ABIDE subjects 51050 and 51050) also reported in Supplementary Materials. See Figs. S12 and S13 for HCP subject 118730, Figs. S14 and S15 for NYU-ABIDE subject 51050 and Figs. S16 and S17 for NYU-ABIDE 51055.

72The DSE tables for minimally and fully preprocessed (Table 6) gives concise summaries of the data quality. The RMS values provide concrete values that can be used to build intuition for data from a given scanner or protocol. The total noise standard deviation falls from 5.015 to 3.437 with clean-up, but it is notable that the fast component, D-var, falls only slightly from 2.598 to 2.406 (in RMS units), while slow variability falls dramatically from about 4.287 to 2.454. This indicates that much of the variance reduction in “cleaning” comes from removal of low frequency drifts and other slowly-varying effects. The magnitude of temporally structured noise is reflected by S-var explaining 73% of total sum-of-squares, and after clean-up S-var and D-var fall into line around 50%. A measure of the spatially structured noise is the global -var that, while small as a percentage, is seen to be about 1500 that expected with IID before preprocessing, and falling to about 275 relative to IID after preprocessing. That the majority of -var is due to -var indicates that the global signal is generally low frequency in nature.

73We also show the DSE tables for three other subjects; HCP subject 118730 in Table S5, NYU-ABIDE subject 51050 and 51055 in Tables S6 and S7, respectively.

74We observe that the cleaned data has, which implies that the average lag-1 autocorrelation is close to zero (Sec. Appendix D.8). However, temporal autocorrelation is a ubiquitous feature of fMRI data, suggesting a contradiction. To address this, Fig. 9 shows maps of the lag-1 temporal autocorrelation across the pre-processing steps. For raw data, the autocorrelation coefficient is between 0.4 and 0.6, but with successive pre-processing steps, the autocorrelation coefficient decreases until the FFP level where the median of voxel-wise autocorrelation coefficients is approximately zero. (See Fig. S18 for similar results on 20 HCP subjects).

75Thus, while temporal autocorrelation is present in the data, we find that the lag-1 autocorrelation coefficients do get close to zero with cleaned data, indicating that the heuristic is correctly indicating negligible average autocorrelation.

76Fig. 10 illustrates the use of the DSE decomposition to summarize the DSE components of 100 unrelated subjects in the HCP cohort, normalized as a percentage of total variability, A-var, to be maximally comparable across subjects. (See Fig. S22 for same results for ABIDE-NYU cohort). A non-normalized version of this plot (Fig. S23) is useful for viewing absolute changes, showing that S-var dramatically drops with preprocessing while -var is relatively stable.

77For the raw data, %D-var ranges from just over 5%–40%, and S-var varies between 60% and 96%; the var only ever explains a negligible portion of the sum-of-squares, 0.027 to across all three pre-processing levels. For all but two subjects the %D-var and %S-var components successively converge to 50% 5% for FPP data.

78Considering only the global variability, the slow %-var is small, usually falling well below 1%, and fast %-var is negligible, never exceeding 0.1%. This reflects the low frequency nature of the global signal. Similarly, the global edge component, %-var only explains a small proportion of the global variabilities.

79To demonstrate the utility of the DSE decomposition in data quality control, we isolate four subjects and observe how their DSE values change with successive preprocessing.

80Subject 151627, marked with a square, is one of the most extreme subjects for S-var and D-var in raw and MPP data, but has one of the smallest %S-var %D-var differences for FPP data. This dramatic reduction in autocorrelation is confirmed in Fig. 11-A, showing the cumulative distribution of lag-1 autocorrelation, and is likely linked to physiological noise around brain stem and other inferior regions (Fig. 12-A1) successfully removed by ICA-FIX (Fig. 12-A2).

81Subject 122620, marked with an triangle, has small %S-var %D-var differences for all versions of the data, also reflected in its distribution of autocorrelation (Fig. 11-B). However, there is still some notable spatial structure in the S-var and D-var images even after clean up (Fig. 12-B2). This illustrates that if a small portion of the image possess problems, it may not be detected in any simple summary.

82Subject 135932, marked with circle, has absolutely typical S-var and D-var among the 100 subjects in the raw and MPP data, but in the FPP data it has one of worst %S-var %D-var differences. The distribution of autocorrelation coefficients reflects this (Fig. 11-C), with FPP (red line) having more large values of ρ than the other subjects. Inspection of the raw data S-var map (Fig. reffig:VarImg-C1) shows evidence of substantial structured noise that is, by in large, mostly removed by ICA-FIX correction (Fig. 12-C2). However the FPP S-var map shows vascular structure, likely a branch of the posterior cerebral artery near the lingual gryus; this is likely an element of physiological noise that ICA-FIX would have ideally removed but missed. Note also that this subject has low movement as measured by median FD (Fig. S21), eliminating motion as the likely source of the problem.

83Finally, subject 101107, marked as a diamond, has the worst quality as measured by divergent %S-var and %D-var across preprocessing levels, with FPP level having S-var and D-var, and reflected in the largest autocorrelation values among the four subjects (Fig. 11-D). Images of S-var show substantial structured variability that remains even in the FPP data (Fig. 12-D), while the D-var image is improves notably with ICA-FIX. (This was a high-motion subject; note loss of ventromedial prefrontal cortex).

84DSE time series plots of these four subjects confirm these findings, with 122620 and 151627 having flat and converged S-var and D-var time series, while 135932 and especially 101107 have structured and diverged S-var time series (Figs. S19 and S20).

85To demonstrate the value of the S-var time series, Fig. 13 explores time points where is particularly large and small for subject 135932. Four “ images” are shown, for voxel i, the constituents of (Eqn. (4)). Panel A of Fig. 13 shows a ’clean’ time point, with a minimum of structured noise apparent, while panels B–D all show a similar vascular pattern. Examination of the ICA components fed into FIX finds 3 components that reflect this vascular structure that were classified as ’good’ (Fig. S24). This demonstrates the value of the DSE decomposition to identify subtle structured noise in the data.

86Finally, in addition to using DSE plots to investigate the quality of scans across pre-processing levels, they can also be used as a universal measure to compare the quality of scans across cohorts, data-sets and pipelines. We computed the DSE decomposition of 530 healthy subjects across 20 acquisition sites in the ABIDE dataset (Fig. S25), identifying particular sites (e.g. NYU & OHSU) and the CPAC preprocessing pipeline generally to have minimal temporal autocorrelation as reflected in S-var/D-var divergence.

87We have provided a formal context for the diagnostic measure DVARS, showing to be part of a decomposition of sum-of-squares at each successive scan pair and over the whole 4D data. We have proposed a significance test for the DVARS measure which, when detected scans are removed based on p-values, we found to address corruptions of FC while preserving temporal degrees of freedom better than other arbitrary approaches. We have also proposed the DSE decomposition which is particularly useful for summarizing data quality via DSE plots and DSE tables. These tools concisely summarize the interplay of the fast, slow, total and global sum-of-squares, and our derived nominal expected values for each table entry facilitates the identification spatial and temporal artifacts.

88Our analysis shows that D-var (and DVARS) scales with overall noise variance, and is deflated by temporal autocorrelation. We observe that as data becomes cleaner, and the background noise falls, we have greater power to identify spikes. Therefore, to avoid ’over-cleaning’ the data we complement the statistical significance of DVARS p-values with the practical significance of %D-var, a standardized measure of the excess variance explained by a spike as a percentage of average variance. Consequently, the final candidate time-points to be scrubbed is a conjunction of statistical and practical significance; we choose a 5% familywise error rate significance level via Bonferroni and a 5% %D-var cut-off; this practical significance threshold worked adequately in the HCP data we examined but may need to be recalibrated for other data sources.

89Yet one more advantage of using tests, proposed in this work, is that we can estimate the effective spatial degrees of freedom which may prove to be a useful index of spatial structure in the data, but we stress this particular degrees-of-freedom ν is specific to this setting and is unlikely to be useful in other contexts (e.g. as a Bonferroni correction over space).

90Besides %D-var, we have introduced two standardised measures which facilitate the inter-cohort comparison of the fast (or DVARS) component regardless of intensity normalisation used in the pre-processing pipelines. For example, standardised measure -var shows the proportion of variability which can be explained via fast component while -var shows the similar proportion for the slow variability in data.

91The DSE plots allow D-var to be judged relative to S-var, checking for convergence to approximately 50% of A-var as data approaches temporal independence, and consequently the level of autocorrelation as measure of corruption can be tightly monitored across pre-processing steps.

92Using DSE plots we found two HCP subjects (101107 & 136932) where the motion-parameter regression and further ICA-FIX algorithm failed to clean the data and clearly stand out from others in the 100 unrelated subject cohort. We have used the DSE variability images to temporally and spatially locate the corruptions. It is important to note that the DSE decomposition technique should only be used before any form of resting-state bandpass filtering (such as 0.01Hz-0.1Hz) and autocorrelation modelling (such as FILM pre-whitening techniques).

93Finally, we stress that we don not believe there is any one strategy can address all fMRI artifacts. Each method used in this work has it is merits and pitfalls. For example, while scrubbing was shown to be useful to remove the head motion induced spikes, it fails to remove the nuisance due to physiological signals on it is own and requires alternatives like ICA-based methods. Regardless of method, we still see value of using DSE plots and images throughout the analysis to choose a right combination of methods; see for a recent comparison of various combinations of artifact methods.

94Our DVARS p-values depend critically on accurate estimates of and. Despite finding exact expressions for the null mean and variance, we found the most practical and reliable estimates to be based on the sample time series itself, using median for and hIQR to find.Of course this indicates that our inference procedure can only infer relative to the background noise level of the data, picking out extreme values that are inconsistent with our approximating approximation.

95There are two essentials avenues as continuation of this work. First is to study the effect of global signal regression via DSE decompositions. As regressing out the global mean deflates the global segment of each variability component, the DSE decomposition can be used to investigate whether global signal regression is helpful to suppress the spatial artifacts. Second, both cleaning algorithms used in this work, scrubbing and ICA-FIX, leave down-spikes (or dips) after regressing out the nuisance. These down-spikes may also affect the FC and could be detected with a two-sided variant of our hypothesis test.

96In this work majority of the analysis have been done on MATLAB versions 2015b and 2016b, supported by FSL 5.0.9 for neuroimaging analysis.

97Inference on DVARS as well as DSE decomposition techniques proposed in this paper are available via MATLAB scripts, found at http://www.github.com/asoroosh/DVARS. Also, a dedicated web page, http://www.nisox.org/Software/DSE/, present the DSE decompositions of HCP and ABIDE cohort and is regularly updated with new publicly available resting-state data sets.

98Results and figure scripts presented in this work are publicly available for reproducibility purposes on http://www.github.com/asoroosh/DVARS_Paper17.