PCA leverage: outlier detection for high-dimensional functional magnetic resonance imaging data

# PCA leverage: outlier detection for high-dimensional functional magnetic resonance imaging data

AMANDA F. MEJIA, MARY BETH NEBEL, ANI ELOYAN, BRIAN CAFFO and MARTIN A. LINDQUIST
Department of Biostatistics, Johns Hopkins University, USA
Center for Neurodevelopmental and Imaging Research, Kennedy Krieger Institute, USA
Department of Biostatistics, Brown University, USA
mlindqui@jhsph.edu
###### Abstract

Outlier detection for high-dimensional (HD) data is a popular topic in modern statistical research. However, one source of HD data that has received relatively little attention is functional magnetic resonance images (fMRI), which consists of hundreds of thousands of measurements sampled at hundreds of time points. At a time when the availability of fMRI data is rapidly growing—primarily through large, publicly available grassroots datasets—automated quality control and outlier detection methods are greatly needed. We propose PCA leverage and demonstrate how it can be used to identify outlying time points in an fMRI run. Furthermore, PCA leverage is a measure of the influence of each observation on the estimation of principal components, which are often of interest in fMRI data. We also propose an alternative measure, PCA robust distance, which is less sensitive to outliers and has controllable statistical properties. The proposed methods are validated through simulation studies and are shown to be highly accurate. We also conduct a reliability study using resting-state fMRI data from the Autism Brain Imaging Data Exchange and find that removal of outliers using the proposed methods results in more reliable estimation of subject-level resting-state networks using ICA.High-dimensional statistics; fMRI; Image analysis; Leverage; Outlier detection; Principal component analysis; Robust statistics

footnotetext: To whom correspondence should be addressed.\@sect

section1[Introduction]Introduction

Outliers in high-dimensional (HD) settings, such as genetics, medical imaging and chemometrics, are a common problem in modern statistics and have been the focus of much recent research (Hubert and others, 2005; Filzmoser and others, 2008; Hadi and others, 2009; Shieh and Hung, 2009; Fritsch and others, 2012; Ro and others, 2015). One such source of HD data is functional magnetic resonance imaging (fMRI). An fMRI run usually contains 100,000-200,000 volumetric elements or “voxels” within the brain, which are sampled at hundreds of time points. Here, we consider voxels to be variables and time points to be observations, in which case the outlier problem is to identify time points that contain high levels of noise or artifacts.

Multiple noise sources related to the hardware and the participant (Lindquist and others, 2008) can corrupt fMRI data, including magnetic field instabilities, head movement, and physiological effects, such as heartbeat and respiration. Noise sources appear as high-frequency “spikes”, image artifacts distortions, and signal drift. fMRI data also undergo a series of complex preprocessing steps before being analyzed; errors during any one of these steps could introduce additional artifacts. Thus, performing adequate quality control prior to statistical analysis is critical.

In recent years, the availability of fMRI data has increased rapidly. The emergence of a number of publicly available fMRI databases, often focusing on a specific disease or disorder, presents a great opportunity to study brain function and organization. However, these datasets are usually collected from multiple sites with varying methods for acquisition, processing and quality control, resulting in widely varying levels of quality and high rates of artifacts. In the absence of automated outlier detection methods appropriate for fMRI data, quality inspection often takes place in a manual or semi-automated manner by individual research groups. This presents a timely opportunity for statisticians to develop more automated methods.

Here we propose an HD outlier detection method based on dimension reduction through principal components analysis (PCA) and established measures of outlyingness, namely leverage and robust distances. While leverage has not typically been employed for outlier identification outside of a regression framework, we argue for leverage as a meaningful measure when the principal components (PCs) are themselves of interest, which is often true for fMRI data.

Several outlier detection methods for standard and HD data use PCA, including PCA influence functions and other PC sensitivity measures (Brooks, 1994; Gao and others, 2005). However, these methods are often computationally demanding as they rely on re-estimating the PCs with each observation left out. Similarly, methods that depend on robust covariance estimation (see Hadi and others (2009) for a review) are usually not suited for HD settings. One such method, the minimum covariance determinant (MCD) estimator, identifies the observation subset with the smallest sample covariance matrix determinant (Rousseeuw, 1985). Hubert and others (2005) proposed ROBPCA, a robust PCA method for HD data that can also identify outliers, which lie far from the robust PCs space. Filzmoser and others (2008) proposed PCOut and Sign, two computationally efficient methods that perform standard PCA after robustly scaling the data and looking for outliers within the principal directions explaining 99% of the variance. Ro and others (2015) proposed the minimum diagonal product estimator, which is related to the MCD but ignores off-diagonal elements and is identifiable when there are more variables than observations. Fritsch and others (2012) proposed an HD adaptation of the MCD through regularization and applied the method to neuroimaging summary statistics.

However, such methods are often validated using only moderately sized data containing more observations than variables. One exception comes from Shieh and Hung (2009) who proposed identifying outlying genes in microarray data by performing PCA dimension reduction prior to robust distance computation on the reduced data. The method was validated on a dataset of approximately 100 observations and 2000 variables and shown to result in fewer false positives and false negatives than ROBPCA.

Existing methods for fMRI artifact identification have focused on head motion and ad-hoc measures of quality. While the removal of affected time points (“scrubbing” or “spike regression”) using these methods appears beneficial (Satterthwaite and others, 2013; Power and others, 2014), a more unified outlier detection framework is needed, as motion is only one potential artifact source in fMRI data. In addition, existing methods result in a collection of measures that must somehow be combined. We propose a single measure of outlyingness related to the influence of each time point on PC estimation, which is the basis of several common brain connectivity measures (see Section PCA leverage: outlier detection for high-dimensional functional magnetic resonance imaging data).

The remainder of this paper is organized as follows. We begin with a description of our statistical methodology. We then present a simulation study, which is used to assess the sensitivity and specificity of the proposed methods. Next, we present a reliability analysis employing the Autism Brain Imaging Data Exchange (ABIDE) dataset. We conclude with a brief discussion.

\@sect

section1[Methods]Methods

As described in detail below, we propose two PCA-based measures of outlyingness, PCA leverage and PCA robust distance, and develop thresholding rules to label outliers using either measure. For both measures, we begin with PCA dimension reduction. All computations are performed in the R statistical environment version 3.1.1 (R Core Team, 2014).

\@sect

subsection2[Dimension Reduction]Dimension Reduction

Let be the number of 3-dimensional “volumes” collected over time in an fMRI run, and let be the number of voxels in the brain. Let () represent an fMRI run, where each row of is a vectorized volume. We first center and scale each column of relative to its median and median absolute deviation (Hampel and others, 1986), respectively, to avoid the influence of outliers. The singular value decomposition (SVD) (Golub and Reinsch, 1970) of is given by , where is diagonal with elements and . Here denotes the transpose of matrix . The rows of contain the PCs or eigenimages of , and the columns of contain the corresponding PC scores. Note that to avoid memory limitations, rather than compute the SVD of directly, one generally computes the SVD of to obtain and and then solves for .

We retain principal components, so that the “reduced data” are given by the submatrices of and corresponding to the first principal components. For ease of notation, we redefine and to represent these submatrices and . To choose the model order , we retain only components with a greater-than-average eigenvalue. While more sophisticated cutoff methods exist, we find that this simple cutoff rule works well in practice (Jackson, 1993). To avoid extreme solutions, we require .

\@sect

subsection2[Principal components analysis leverage]Principal components analysis leverage

In regression, leverage is defined as the diagonals of the “hat matrix” , where is a matrix of explanatory variables (Neter and others, 1996). The hat matrix projects the outcome variable(s) onto the column space of , yielding the projected data . Leverage, bounded between and , is often used to assess the potential influence of an observation on the regression fit, as it is the change in due to a 1-unit change in and is proportional to the uncertainty in the estimate, since . Particularly relevant to our context, leverage is also a measure of outlyingness among the explanatory variables, as it is related to the Mahalanobis distance.

Extending leverage to the PCA context, we treat as an estimated design matrix in the estimation of . With and fixed, is equivalent to the least squares estimate in the multivariate regression model . We therefore define PCA leverage as , where . Note that is simply a scaling factor applied to each variable and therefore has no effect on leverage. Continuing the regression analogy, in PCA, projects onto the column space of , the principal directions, as . Furthermore, PCA leverage is a measure of outlyingness among the PCA scores and within , since . While in reality and are not fixed and PCA leverage therefore only approximately represents the influence of each observation on the PCs and fitted values in , we find this approximation to be quite close in practice, as illustrated in Figure 1. Note that dimension reduction is essential for PCA leverage to be informative, since when all components are retained.

In regression, leverage only represents the potential influence of an observation on regression coefficient estimation; influence points must be outliers in the explanatory variables (“leverage points”) as well as in the response variable(s). In contrast, PCA leverage is a more direct measure of influence, as PCA leverage points are outliers in both and the original data . Thus, PCA leverage points are also influence points. Furthermore, while in regression we discern “good” from “bad” leverage points, fMRI observations with high PCA leverage are unlikely to represent true signal, since the signal change associated with neuronal sources is very small compared with noise and artifacts. Therefore, we assume that all observations with high PCA leverage are “bad” influence points in the fMRI context.

Moreover, the interpretation of PCA leverage as the influence of each observation on PC estimation is particularly relevant for resting-state fMRI (rs-fMRI). PC estimation is a preprocessing step for one of the most common types of analysis for rs-fMRI data: estimation of spatially independent brain networks and the functional connectivity between those networks. In such analyses, PCA leverage is both a measure of influence on the quantity of interest and of outlyingness.

In setting a leverage threshold to identify outliers, it is important to recognize that the sum of leverages for a set of observations equals the number of variables in the design matrix. The mean leverage of all observations is therefore fixed at . Outliers wield a large amount of leverage, and their presence reduces the leverage of all remaining observations, such that the mean may be significantly greater than the leverage of normal observations. Thus, the median is a more appropriate reference quantity than the mean for normal observations. If the leverage of observation exceeds times the median , it is labeled a “leverage outlier”. In the simulations and experimental data analysis described below, we consider .

While such practical rules may work well in the absence of convenient statistical properties, a formal statistical test for outliers with known and controllable properties is desirable. In the following section, we propose an alternative robust distance measure based on minimum covariance determinant (MCD) estimators (Rousseeuw, 1985).

\@sect

subsection2[Principal components robust distance]Principal components robust distance

For a design matrix with an intercept or centered variables, leverage is related to the squared empirical Mahalanobis distance (Mahalanobis, 1936), defined for an matrix and observation as , where and are the sample mean and covariance matrix of , respectively. The Mahalanobis distance is known to be sensitive to outliers due to their influence on the sample mean and covariance, which may lead to “masking”, a phenomenon in which truly outlying observations appear normal due to the presence of more extreme outliers (Rousseeuw and Van Zomeren, 1990; Rousseeuw and Hubert, 2011).

As an alternate measure, we adopt the minimum covariance determinant (MCD) distance proposed by Rousseeuw (1985). For a general dataset, let be the number of observations and be the number of variables. The MCD estimators of location, , and scale, , are obtained by computing the sample mean and covariance within a subset of the data of size for which the confidence ellipsoid determined by and centered at has minimal volume. The maximum breakdown point of MCD estimators is obtained by setting and approaches as . The MCD distance is then computed as a Mahalanobis distance using and in place of and . For ease of notation, let .

Let , and let , , be the indices of the observations selected to compute and . Let be the indices of the remaining observations, among which we look for outliers. For Gaussian data, approximately follow a distribution (Hubert and others, 2005; Shieh and Hung, 2009), while for ,

 ~d2i:=c(m−p+1)pmd2i∼Fp,m−p+1,

where and can be estimated asymptotically or through simulation. (While some previous work has simply assumed a distribution for , we find this to result in many false positives.) To estimate we use the asymptotic form, , while to estimate we use the small sample-corrected asymptotic form given in Hardin and Rocke (2005). To improve the F-distribution fit, Maronna and Zamar (2002) and Filzmoser and others (2008) scale the distances to match the median of the theoretical distribution. However, as contains at most half of the original observations, the median within represents the th or greater quantile within and may be contaminated with outliers. Therefore, we let , where is the th sample quantile of . We label a “distance outlier” any observation in with greater than the th quantile of the theoretical F distribution. In our simulations and experimental data analysis, we consider ranging from to .

For time series data, where the assumption of independence among observations is violated, the distributional results given in Hardin and Rocke (2005) may be invalid. As the autocorrelation in fMRI time series is often modeled as an AR(1) process with a coefficient of , we divide each fMRI time series into three subsets, each consisting of every third observation. The autocorrelation within each subset is negligible at approximately . To obtain the MCD distance for each observation, we use the MCD estimates of center and scale within each subset, averaged across subsets, and we find that this significantly improves the distributional fit.

\@sect

section1[Simulation Study]Simulation Study

\@sect

subsection2[Construction of baseline scans]Construction of baseline scans

Our simulated dataset is based on fMRI scans from three subjects collected as part of the ABIDE dataset (described in Section PCA leverage: outlier detection for high-dimensional functional magnetic resonance imaging data). For generalizability, each subject was chosen from a different data collection site. For each scan, we identify a contiguous subset of volumes containing no detectable artifacts, resulting in 141, 171 and 89 volumes, respectively. We reduce dimensionality by using only the 45th axial (horizontal) slice, corresponding roughly to the center of the brain. For scan , let be the resulting length of the scan and be the resulting number of voxels in the brain mask, so that scan is represented by the matrix . We can separate into an anatomical baseline, the mean image (), and the residual , representing primarily functional information. Then , where is a vector of s of length .

We then use independent components analysis (ICA), a blind-source separation algorithm, to decompose the intrinsic activity in into a number of spatially independent sources (McKeown and others (1997)) (described in Section PCA leverage: outlier detection for high-dimensional functional magnetic resonance imaging data). Let be the number of sources of neuronal signal identified for scan . Then , where () contains the spatial maps of each source, and () contains the time courses of each source. The residual contains structured (spatially and temporally correlated) noise. Let .

\@sect

subsection2[Artifact-free images]Artifact-free images

For each scan , we construct three simulation setups: baseline image () plus white noise (setup 1); baseline image plus functional signal () plus white noise (setup 2); and baseline image plus functional signal plus structured noise () (setup 3).

To test the specificity of each outlier detection method in the artifact-free setting, we generate images with varying signal-to-noise ratio (SNR) in the following way. For scan , we have true signal variance and true noise variance . Defining SNR as the ratio of signal variance to noise variance, let 0.025, 0.050, 0.075, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0 be the desired SNR of the simulated scans. For setups 1 and 2, we generate the white noise matrix for scan as independent, mean-zero Gaussian noise with variance . For setup 3, we generate the structured noise matrix , where is the baseline SNR of scan , equal to , and , respectively. Therefore, the simulated artifact-free data at SNR is for setup 1; for setup 2; and for setup 3. For setups 1 and 2, we randomly generate times; for setup 3, the noise is fixed.

The specificity, or percentage of observations not labeled as outliers that are truly non-outliers, in this case is simply the percentage of volumes in each scan not labeled as outliers. Figure 3 shows the mean specificity across iterations, where each line represents a scan and SNR level. The lines in red correspond to SNR of , which is close to the observed SNR of each scan. Specificity is nearly for both leverage and robust distance methods across all thresholds and SNR levels considered. In the presence of structured noise, the specificity of the robust distance method is approximately - in some cases, unless the th quantile threshold is used.

\@sect

subsection2[Images with artifacts]Images with artifacts

We generate four common types of fMRI artifacts: spikes, motion, banding, and ghosting. Spikes are created by increasing the intensity of an entire volume by a given percentage. Motion artifacts are created by rotating a volume by a given angle. Banding artifacts are generated by changing an intensity in the Fourier transform of the image, resulting in a striped appearance. Ghosting artifacts are created by superimposing a figure or “ghost” moving through space over time.

At each of iterations, one simulated fMRI scan is generated for each subject, SNR level, artifact type and simulation setup. For spike, motion and banding artifacts, 10 volumes are randomly selected from each scan, and the artifact intensity for each volume is generated from a uniform distribution (range percent intensity increase for spike artifacts; rotation for motion artifacts; times change at location of the Fourier transformed image). For ghosting artifacts, 9 sequential volumes are randomly selected, and the mean intensity of the ghost, relative to the mean intensity of the image, is randomly generated from a uniform distribution (range ). An example of each artifact type is displayed in Figure 2.

We are interested in both the specificity and the sensitivity, or the percentage of true outliers identified as outliers. Figure 4 shows the mean sensitivity and specificity for each outlier detection method, simulation setup, and artifact type, where each line represents a scan and SNR level. The realistic SNR of is shown in red. As the simulation setup becomes more realistic, the sensitivity to outliers tends to decrease, while the specificity is relatively stable. The robust distance method has nearly specificity in all scenarios and tends to display higher sensitivity than the leverage method, particularly for banding and spike artifacts. While differences across artifact types are apparent, these are likely driven by the range of intensities chosen.

\@sect

section1[Experimental Data Results]Experimental Data Results

Using a large, multi-site fMRI dataset, we assess the result of outlier removal on the scan-rescan reliability of a common type of analysis. This section is organized as follows. We begin with a description of the dataset employed and show an example. We then describe the reliability analysis. Finally, we quantify the improvement to reliability with the proposed outlier detection methods using a linear mixed model to account for subject and site effects.

\@sect

subsection2[fMRI Dataset]fMRI Dataset

ABIDE is a publicly available resource of neuroimaging and phenotypic information from 1112 subjects consisting of 20 datasets collected at 16 sites (Di Martino and others, 2014). Fully anonymized data from 91 children collected at Kennedy Krieger Institute after the ABIDE release were also included. Image acquisition parameters and demographic information are available at http://fcon_1000.projects.nitrc.org/indi/abide/. For each subject, a -weighted MPRAGE volume and one or more rs-fMRI sessions were collected on the same day. Details of data pre-processing and quality control are provided in Supplementary Materials Appendix A, where Table 1 lists the number of subjects in each dataset.

For a single example scan, Figure 5 shows the leverage and robust distance functions, along with 6 motion parameters (roll, pitch, yaw, and translation in each direction) and their derivatives, which are commonly used for artifact detection. Below the plot, the volumes corresponding to the spikes at time points 60, 90, 134 and 150 are shown. Three of the spikes are leverage and distance outliers using any of the thresholds considered ( for leverage; for robust distance), while the spike at time point 90 is only a leverage outlier at . Obvious banding artifacts are seen at time points 60 and 150, a moderate banding artifact is seen at time point 134, and no visible artifact is apparent at time point 90. While the artifact at time point 150 would be detected using motion measures, the other spikes would likely go undetected.

\@sect

subsection2[Estimation of subject-level brain networks and connectivity]Estimation of subject-level brain networks and connectivity

Resting-state brain networks represent regions of the brain that act in a coordinated manner during rest. While such networks have traditionally been identified at the group level, there is growing interest in estimating these networks at the subject level, where the higher levels of noise make accurate estimation difficult. There is also interest in estimating the subject-level “functional connectivity”, or temporal dependence of neuronal activation, between these different networks (van den Heuvel and Pol, 2010). We assess the benefits of outlier removal on the reliability of these networks and their functional connectivity. Details of the estimation of subject-level resting-state networks are provided in Supplementary Materials Appendix B. Here we briefly describe the procedure. We begin by performing group ICA (GICA) separately for each of the ABIDE datasets . The result of GICA is the matrix , where is the number of voxels in the group-level brain mask and is the number of independent components (ICs). Each row of may represent a source of noise (e.g. motion, respiration) or a resting-state network. After identification of those ICs corresponding to resting-state networks, let denote the matrix containing the resting-state networks identified for dataset . Using dual regression (Beckmann and others, 2009), we obtain , where is the matrix whose rows contain the estimated resting-state brain networks for subject , and is the “mixing matrix” representing the activation of each network over time. We are interested in reliable estimation of two quantities: and the matrix , which represents the functional connectivity between each pair of networks.

\@sect

subsection2[Measuring reliability of subject-level brain networks and functional connectivity]Measuring reliability of subject-level brain networks and functional connectivity

Let and be two sets of estimated resting-state networks for subject obtained by performing dual regression separately for two different scanning sessions of subject in dataset . There is no need to match components between and , since the ICs in each correspond to the same group-level ICs in . To assess reliability, for each subject and component we compute the number of overlapping voxels between and after both have been thresholded (as described in Supplementary Materials Appendix B). We then average over all networks to obtain the average scan-rescan overlap for each subject per network, denoted for subject in dataset using outlier removal method . Outlier removal methods include no outlier removal, leverage-based outlier removal with , and robust distance-based outlier removal with . Similarly, to assess reliability of functional connectivity between networks, let and be the mixing matrices corresponding to and . We compute the mean squared error (MSE) between the upper triangles of and and denote the result for subject in dataset using outlier removal method .

Although most subjects in the ABIDE dataset have only a single scanning session, we can simulate scan-rescan data by splitting each subject’s data into two contiguous subsets consisting of the first and second half of time points, respectively. We use the resulting pseudo scan-rescan data to obtain , , and , then compute our reliability measures and . While this approach may produce an optimistic estimate of the true scan-rescan reliability, this is not a concern as we are primarily interested in the change in reliability due to outlier removal.

To test for changes in reliability of resting-state networks or functional connectivity due to outlier removal, we fit a linear mixed effects model with a fixed effect for each outlier removal method, a fixed effect for each dataset, and a random effect for each subject. We employ this model for its ability to test several groups and methods simultaneously and to account for within-subject correlation across methods. Using the subjects for whom at least one outlier was identified using any method, we estimate the following model for :

 Mikm =bi0+γk+αmIm>0+ϵikm, ϵikm∼N(0,σ2), bi0∼N(0,τ2),

where indicates no outlier removal. Here, represents the baseline reliability for subjects in dataset when no outlier removal is performed, and represents the change in reliability when outlier removal method is used. To obtain coefficient estimates, we fit this model using the lme function from the nlme package (Pinheiro and others, 2014). Since we have a large sample size, we compute Normal 95% confidence intervals.

\@sect

subsection2[The effect of outlier removal on reliability]The effect of outlier removal on reliability

Figure 6 displays estimates and 95% confidence intervals for the coefficients of the models for reliability of resting state networks (a) and functional connectivity (b). For (a), larger overlap values represent greater reliability; for (b), smaller MSE values represent greater reliability. The left-hand panels of (a) and (b) display the fixed effects for each dataset () and illustrate the heterogeneity in baseline reliability across ABIDE datasets before outlier removal. This reflects the substantial differences in acquisition, processing and quality control methods across the data collection sites contributing to ABIDE. The middle and right-hand panels of (a) and (b) display the coefficients for each outlier removal method (). The average percentage of volumes labeled as outliers using each method is also displayed in gray. Both leverage-based and distance-based methods significantly improve reliability of estimates of subject-level brain networks and functional connectivity. Improvement in reliability is maximized using a threhold of or for leverage-based outlier removal and for distance-based outlier removal. The maximum improvement in overlap of resting-state networks is approximately voxels using either method, while the maximum reduction in MSE of functional connectivity is and is achieved using leverage-based outlier removal with .

We also stratify the model by those subjects who passed quality inspection and those who failed. Figure 1 of Supplementary Materials Appendix C shows estimates and 95% confidence intervals for the model coefficients after stratification. While subjects who failed quality inspection tend to improve more than those who passed quality inspection, differences are not statistically significant, and both groups of subjects benefit substantially from outlier removal.

\@sect

section1[Discussion]Discussion

We have proposed a method to detect outlying time points in an fMRI scan by drawing on the traditional statistical ideas of PCA, leverage, and outlier detection. The proposed methods have been validated through simulated data and a large, diverse fMRI dataset. We have demonstrated that the proposed methods are accurate and result in improved reliability of two common types of analysis for resting-state fMRI data, namely identification of resting-state networks through ICA and estimation of functional connectivity between these networks.

The proposed techniques are, to the best of our knowledge, the first to provide a single measure of outlyingness for time points in an fMRI scan, which can be easily thresholded to identify outliers. Unlike motion-based outlier detection methods for fMRI, they are agnostic to the source of artifact and therefore may be used as a general method to detect artifacts, including those unrelated to motion. Furthermore, PCA leverage is directly related to the estimation of principal components, which are an important quantity in the analysis of resting-state fMRI, as they are used as input to ICA for the identification of resting-state networks.

One limitation of our approach is that we perform validation on a single dataset, the ABIDE. However, this dataset is in fact a diverse collection of 20 datasets from 16 international sites, which strengthens the generalizability of our results. Another limitation of the proposed methods is that they may be sensitive to the number of principal components retained. However, we have found that the method performs well with different model orders (e.g. or ), and we propose an automated method of selecting model order in order to provide a fully automated approach.

A limitation of PCA leverage-based outlier removal is that the proposed thresholding rule is, as in regression, somewhat ad-hoc. However, the use of the median leverage across observations as a benchmark is a reasonable approach, and we have tested a range of threholds. Based on our reliability analysis we expect a threshold of of times the median leverage to work well in practice for fMRI, but for different types of data the researcher may wish to re-evaluate this choice. In particular, fMRI volumes containing artifacts tend to be very different from those volumes free of outliers, so a relatively high threshold for leverage and robust distance tends to work well; in other contexts, outliers may be more subtle.

While the proposed methods have been designed and validated for resting-state fMRI data, they may be easily extended to other types of medical imaging data, such as task fMRI and EEG data, as well as other types of HD data. Furthermore, they may also be extended to group analyses. Future work should focus on exploring these directions.

As the availability of large fMRI datasets continues to grow, automated outlier detection methods are becoming essential for the effective use of such data. In particular, the reliability of analyses employing these diverse datasets may be negatively impacted by the presence of poor quality data. The outlier detection methods we propose have the potential to improve the quality of such datasets, thus enhancing the possibilities to use these data to understand neurological diseases and brain function in general.

\@ssect

Supplementary Materials The reader is referred to the on-line Supplementary Materials for additional technical details. Appendix A contains details regarding the fMRI data acquisition, pre-processing, and quality control procedures. It also contains information about each dataset forming the ABIDE. Appendix B details the estimation of subject-level resting-state networks through ICA and dual regression. Appendix C displays the results of reliability analysis after stratifying by subjects whose scans passed quality inspection and those whose did not.

\@ssect

Funding This work was supported by the National Institutes of Health [R01 EB016061], the National Institute of Biomedical Imaging and Bioengineering [P41 EB015909], and the National Institute of Mental Health [R01 MH095836]. Conflict of Interest: None declared.

## References

• Beckmann and others (2009) Beckmann, Christian F, Mackay, Clare E, Filippini, Nicola and Smith, Stephen M. (2009). Group comparison of resting-state fMRI data using multi-subject ICA and dual regression. NeuroImage 47(Suppl 1), S148.
• Brooks (1994) Brooks, Stephen P. (1994). Diagnostics for principal components: influence functions as diagnostic tools. The Statistician, 483–494.
• Di Martino and others (2014) Di Martino, Adriana, Yan, Chao-Gan, Li, Qingyang, Denio, Erin, Castellanos, Francisco X, Alaerts, Kaat, Anderson, Jeffrey S, Assaf, Michal and others. (2014). The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism. Molecular Psychiatry 19(6), 659–667.
• Filzmoser and others (2008) Filzmoser, Peter, Maronna, Ricardo and Werner, Mark. (2008). Outlier identification in high dimensions. Computational Statistics & Data Analysis 52(3), 1694–1711.
• Fritsch and others (2012) Fritsch, Virgile, Varoquaux, Gaël, Thyreau, Benjamin, Poline, Jean-Baptiste and Thirion, Bertrand. (2012). Detecting outliers in high-dimensional neuroimaging datasets with robust covariance estimators. Medical image analysis 16(7), 1359–1370.
• Gao and others (2005) Gao, Shaogen, Li, Guoying and Wang, Dongqian. (2005). A new approach for detecting multivariate outliers. Communications in Statistics: Theory and Methods 34(8), 1857–1865.
• Golub and Reinsch (1970) Golub, Gene H and Reinsch, Christian. (1970). Singular value decomposition and least squares solutions. Numerische mathematik 14(5), 403–420.
• Hadi and others (2009) Hadi, Ali S, Imon, AHM and Werner, Mark. (2009). Detection of outliers. Wiley Interdisciplinary Reviews: Computational Statistics 1(1), 57–70.
• Hampel and others (1986) Hampel, Frank R, Ronchetti, Elvezio M, Rousseeuw, Peter and Stahel, Werner A. (1986). Robust Statistics: the Approach based on Influence Functions. New York: John Wiley.
• Hardin and Rocke (2005) Hardin, Johanna and Rocke, David M. (2005). The distribution of robust distances. Journal of Computational and Graphical Statistics 14(4).
• Hubert and others (2005) Hubert, Mia, Rousseeuw, Peter J and Vanden Branden, Karlien. (2005). ROBPCA: a new approach to robust principal component analysis. Technometrics 47(1), 64–79.
• Jackson (1993) Jackson, Donald A. (1993). Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches. Ecology, 2204–2214.
• Lindquist and others (2008) Lindquist, Martin A and others. (2008). The statistical analysis of fMRI data. Statistical Science 23(4), 439–464.
• Mahalanobis (1936) Mahalanobis, Prasanta Chandra. (1936). On the generalized distance in statistics. Proceedings of the National Institute of Sciences (Calcutta) 2, 49–55.
• Maronna and Zamar (2002) Maronna, Ricardo A and Zamar, Ruben H. (2002). Robust estimates of location and dispersion for high-dimensional datasets. Technometrics 44(4).
• McKeown and others (1997) McKeown, Martin J, Makeig, Scott, Brown, Greg G, Jung, Tzyy-Ping, Kindermann, Sandra S, Bell, Anthony J and Sejnowski, Terrence J. (1997). Analysis of fmri data by blind separation into independent spatial components. Technical Report, DTIC Document.
• Neter and others (1996) Neter, John, Kutner, Michael H, Nachtsheim, Christopher J and Wasserman, William. (1996). Applied linear statistical models, Volume 4. Irwin Chicago.
• Pinheiro and others (2014) Pinheiro, Jose, Bates, Douglas, DebRoy, Saikat, Sarkar, Deepayan and R Core Team. (2014). nlme: Linear and nonlinear mixed effects models. R package version 3.1-118.
• Power and others (2014) Power, Jonathan D, Mitra, Anish, Laumann, Timothy O, Snyder, Abraham Z, Schlaggar, Bradley L and Petersen, Steven E. (2014). Methods to detect, characterize, and remove motion artifact in resting state fMRI. NeuroImage 84, 320–341.
• R Core Team (2014) R Core Team. (2014). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
• Ro and others (2015) Ro, Kwangil, Zou, Changliang, Wang, Zhaojun and Yin, Guosheng. (2015). Outlier detection for high-dimensional data. Biometrika.
• Rousseeuw (1985) Rousseeuw, Peter J. (1985). Multivariate estimation with high breakdown point. Mathematical statistics and applications 8, 283–297.
• Rousseeuw and Hubert (2011) Rousseeuw, Peter J and Hubert, Mia. (2011). Robust statistics for outlier detection. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery 1(1), 73–79.
• Rousseeuw and Van Zomeren (1990) Rousseeuw, Peter J and Van Zomeren, Bert C. (1990). Unmasking multivariate outliers and leverage points. Journal of the American Statistical Association 85(411), 633–639.
• Satterthwaite and others (2013) Satterthwaite, Theodore D, Elliott, Mark A, Gerraty, Raphael T, Ruparel, Kosha, Loughead, James, Calkins, Monica E, Eickhoff, Simon B and others. (2013). An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data. NeuroImage 64, 240–256.
• Shieh and Hung (2009) Shieh, Albert D and Hung, Yeung Sam. (2009). Detecting outlier samples in microarray data. Statistical applications in genetics and molecular biology 8(1), 1–24.
• van den Heuvel and Pol (2010) van den Heuvel, Martin P and Pol, Hilleke E Hulshoff. (2010). Exploring the brain network: a review on resting-state fMRI functional connectivity. European Neuropsychopharmacology 20(8), 519–534.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters