Title: | Scrubbing and Other Data Cleaning Routines for fMRI |
---|---|
Description: | Data-driven fMRI denoising with projection scrubbing (Pham et al (2022) <doi:10.1016/j.neuroimage.2023.119972>). Also includes routines for DVARS (Derivatives VARianceS) (Afyouni and Nichols (2018) <doi:10.1016/j.neuroimage.2017.12.098>), motion scrubbing (Power et al (2012) <doi:10.1016/j.neuroimage.2011.10.018>), aCompCor (anatomical Components Correction) (Muschelli et al (2014) <doi:10.1016/j.neuroimage.2014.03.028>), detrending, and nuisance regression. Projection scrubbing is also applicable to other outlier detection tasks involving high-dimensional data. |
Authors: | Amanda Mejia [aut, cre], John Muschelli [aut] , Damon Pham [aut] , Daniel McDonald [ctb], Fatma Parlak [ctb] |
Maintainer: | Amanda Mejia <[email protected]> |
License: | GPL-3 |
Version: | 0.13.0 |
Built: | 2024-10-29 04:18:09 UTC |
Source: | https://github.com/mandymejia/fmriscrub |
Visualize artifact patterns from the results of pscrub
.
Requires pscrub(..., get_dirs=TRUE)
.
artifact_images(psx, idx = NULL, use_dt = TRUE)
artifact_images(psx, idx = NULL, use_dt = TRUE)
psx |
A |
idx |
The timepoints or column indexes for which to compute artifact
images. If |
use_dt |
If detrended components are available (the "U" matrix of PCA
or "M" matrix of ICA), should they be used to compute the artifact images?
Otherwise, use the non-detrended components. Default: |
Computes two types: "mean" artifact images based on a weighted sum of the projection directions, with weights determined by the scores for each component at the flagged timepoint, and "top" artifact images based on the projection direction with the greatest score at the flagged timepoint.
A list of three: idx
, the timepoints for which the artifact images
were computed; mean
, the "mean" artifact images; and top
, the
"top" artifact images. The row names of the top
artifact images
matrix give the index of the top component ("V" in PCA and "S" in ICA) at
each timepoint.
A sagittal slice from the fMRI time series for subject 0050048. The scan was obtained at the University of Pittsburgh School of Medicine. The scan has been pre-processed with slice time correction, rigid body realignment estimation, spatial normalization to MNI space, and linear detrending. Subject 0050048 was a typically-developing 11-year-old male. The scan has many artifacts. A mask was applied to vectorize the spatial dimensions.
Dat1
Dat1
A numeric matrix of 193 time points by 4675 voxels
Source: http://fcon_1000.projects.nitrc.org/indi/abide/abide_I.html
1.Di Martino, A. et al. The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism. Mol Psychiatry 19, 659–667 (2014).
A sagittal slice from the fMRI time series for subject 0051479. The scan was obtained at the California Institute of Technology. The scan has been pre-processed with slice time correction, rigid body realignment estimation, spatial normalization to MNI space, and linear detrending. Subject 0051479 was a typically-developing 20-year-old female. The scan has few visible artifacts. A mask was applied to vectorize the spatial dimensions.
Dat2
Dat2
A numeric matrix of 145 time points by 4679 voxels
Source: http://fcon_1000.projects.nitrc.org/indi/abide/abide_I.html
1.Di Martino, A. et al. The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism. Mol Psychiatry 19, 659–667 (2014).
Computes the DSE decomposition and DVARS-related statistics. Based on code from github.com/asoroosh/DVARS .
DVARS( X, normalize = TRUE, cutoff_DPD = 5, cutoff_ZD = qnorm(1 - 0.05/nrow(as.matrix_ifti(X))), verbose = FALSE )
DVARS( X, normalize = TRUE, cutoff_DPD = 5, cutoff_ZD = qnorm(1 - 0.05/nrow(as.matrix_ifti(X))), verbose = FALSE )
X |
a |
normalize |
Normalize the data? Default: To replicate Afyouni and Nichols' procedure for the HCP MPP data, since the HCP scans are already normalized to 10,000, just divide the data by 100 and center the voxels on their means:
Note that while voxel centering doesn't affect DVARS, it does affect DPD and ZD. |
cutoff_DPD , cutoff_ZD
|
Numeric outlier cutoffs. Timepoints exceeding these cutoffs will be flagged as outliers. |
verbose |
Should occasional updates be printed? Default is |
A list with components
A data.frame with rows, each column being a different variant of DVARS.
"DVARS"
The outlier cutoff value(s).
A logical data.frame with rows, where
TRUE
indicates suspected outlier presence.
Afyouni, S. & Nichols, T. E. Insight and inference for DVARS. NeuroImage 172, 291-312 (2018).
Calculate Framewise Displacement (FD)
FD( X, trans_units = c("mm", "cm", "in"), rot_units = c("deg", "rad", "mm", "cm", "in"), brain_radius = NULL, detrend = FALSE, lag = 1, cutoff = 0.3 )
FD( X, trans_units = c("mm", "cm", "in"), rot_units = c("deg", "rad", "mm", "cm", "in"), brain_radius = NULL, detrend = FALSE, lag = 1, cutoff = 0.3 )
X |
An Alternatively, this can be the file path to an |
trans_units |
|
rot_units |
|
brain_radius |
If If |
detrend |
Detrend each RP with the DCT before computing FD?
Default: |
lag |
The difference of indices between which to calculate change in
position. Default: |
cutoff |
FD values higher than this will be flagged. Default: |
The FD formula is taken from Power et. al. (2012):
where is the timepoint;
,
and
are the
translational realignment parameters (RPs);
,
and
are the rotational RPs;
and
(and similarly for the other RPs).
A list with components
A length vector of FD values in
trans_units
.
"FD"
cutoff
A length-N logical vector, where TRUE
indicates suspected outlier presence.
Power, J. D., Barnes, K. A., Snyder, A. Z., Schlaggar, B. L. & Petersen, S. E. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. Neuroimage 59, 2142-2154 (2012).
Convert flagged volumes to corresponding one-hot encoded vectors which can be used for nuisance regression.
flags_to_nuis_spikes(flags, n_time)
flags_to_nuis_spikes(flags, n_time)
flags |
Numeric vector of integers indicating the indices of
the flagged volumes. Or, a logical vector of length |
n_time |
The length of the vectors to obtain. For nuisance regression,
this is the length of the BOLD data. The highest index in |
A numeric matrix of ones and zeroes. The number of rows will be
n_time
and the number of columns will be the number of flags. Each
column will have a 1
at the flag index, and 0
elsewhere.
See help(package="fMRIscrub")
for a list of functions.
The kurtosis cutoff is a high quantile (default 0.99) of the sampling distribution of kurtosis for Normal iid data of the same length as the components; it is estimated by simulation or calculated from the theoretical asymptotic distribution if the components are long enough.
high_kurtosis(Comps, kurt_quantile = 0.99, n_sim = 5000, min_1 = FALSE)
high_kurtosis(Comps, kurt_quantile = 0.99, n_sim = 5000, min_1 = FALSE)
Comps |
A matrix; each column is a component. For PCA, this is the U matrix. For ICA, this is the M matrix. |
kurt_quantile |
components with kurtosis of at least this quantile are kept. |
n_sim |
The number of simulation data to use for estimating the sampling
distribution of kurtosis. Only used if a new simulation is performed. (If
|
min_1 |
Require at least one component to be selected? In other words, if
no components meet the quantile cutoff, should the component with the highest
kurtosis be returned? Default: |
The components should not have any strong low-frequency trends, because trends can affect kurtosis in unpredictable ways unrelated to outlier presence.
A logical vector indicating whether each component has high kurtosis.
Computes the leverage of each observation in the PC score (U) or IC mixing (M) matrix for projection scrubbing. Can threshold to flag potential outliers.
leverage(Comps, are_norm = FALSE, median_cutoff = NULL)
leverage(Comps, are_norm = FALSE, median_cutoff = NULL)
Comps |
The |
are_norm |
Assume the columns of |
median_cutoff |
The outlier cutoff, in multiples of the median leverage.
Default: |
A list with entries "meas"
(the leverage values),
"cut"
(the leverage cutoff value) and
"flag"
(logical vector indicating the outliers). If
!is.null(median_cutoff)
, "cut"
and "flag"
are omitted.
"scrub_DVARS"
objectPlot a "scrub_DVARS"
object
## S3 method for class 'scrub_DVARS' plot(x, title = NULL, ...)
## S3 method for class 'scrub_DVARS' plot(x, title = NULL, ...)
x |
The |
title |
(Optional) If provided, will add a title to the plot. |
... |
Additional arguments to ggplot, e.g. |
A ggplot
"scrub_FD"
objectPlot a "scrub_FD"
object
## S3 method for class 'scrub_FD' plot(x, title = NULL, ...)
## S3 method for class 'scrub_FD' plot(x, title = NULL, ...)
x |
The |
title |
(Optional) If provided, will add a title to the plot. |
... |
Additional arguments to ggplot, e.g. |
A ggplot
"scrub_projection"
objectPlot a "scrub_projection"
object
## S3 method for class 'scrub_projection' plot(x, title = NULL, ...)
## S3 method for class 'scrub_projection' plot(x, title = NULL, ...)
x |
The |
title |
(Optional) If provided, will add a title to the plot. |
... |
Additional arguments to ggplot, e.g. |
A ggplot
Projection scrubbing is a data-driven method for identifying artifact-contaminated volumes in fMRI. It works by identifying component directions in the data likely to represent patterns of burst noise, and then computing a composite measure of outlyingness based on leverage within these directions, at each volume. The projection can be PCA, ICA, or "fused PCA." Projection scrubbing can also be used for other outlier detection tasks involving high-dimensional data.
pscrub( X, projection = c("ICA", "PCA"), nuisance = "DCT4", center = TRUE, scale = TRUE, comps_mean_dt = FALSE, comps_var_dt = FALSE, PESEL = TRUE, kurt_quantile = 0.99, get_dirs = FALSE, full_PCA = FALSE, get_outliers = TRUE, cutoff = 4, seed = 0, verbose = FALSE )
pscrub( X, projection = c("ICA", "PCA"), nuisance = "DCT4", center = TRUE, scale = TRUE, comps_mean_dt = FALSE, comps_var_dt = FALSE, PESEL = TRUE, kurt_quantile = 0.99, get_dirs = FALSE, full_PCA = FALSE, get_outliers = TRUE, cutoff = 4, seed = 0, verbose = FALSE )
X |
Wide numeric data matrix ( |
projection |
One of the following: |
nuisance |
Nuisance signals to regress from each column of Detrending is highly recommended for time-series data, especially if there are many time points or evolving circumstances affecting the data. Additionally, if kurtosis is being used to select the projection directions, trends can induce positive or negative kurtosis, contaminating the connection between high kurtosis and outlier presence. Detrending should not be used with non-time-series data because the observations are not temporally related. Additional nuisance regressors can be specified like so:
|
center , scale
|
Center the columns of the data by their medians, and scale the
columns of the data by their median absolute deviations (MADs)? Default: Note that centering and scaling occur after nuisance regression, so even if
|
comps_mean_dt , comps_var_dt
|
Stabilize the mean and variance of each
projection component's timecourse prior to computing kurtosis and leverage?
These arguments should be Slow-moving mean and variance patterns in the components will interfere with
the roles of kurtosis and leverage in identifying outliers. While
Overall, for fMRI we recommend enabling |
PESEL |
Use |
kurt_quantile |
What quantile cutoff should be used to select the
components? Default: We model each component as a length |
get_dirs |
Should the projection directions be returned? This is the
|
full_PCA |
Only applies to the PCA projection. Return the full SVD?
Default: |
get_outliers |
Should outliers be flagged based on |
cutoff |
Median leverage cutoff value. Default: |
seed |
Set a seed right before the call to |
verbose |
Should occasional updates be printed? Default: |
Refer to the projection scrubbing vignette for a demonstration and an
outline of the algorithm: vignette("projection_scrubbing", package="fMRIscrub")
A "pscrub"
object, i.e. a list with components
A numeric vector of leverage values.
The numeric outlier cutoff value (cutoff
times the median leverage).
A logical vector where TRUE
indicates where leverage exceeds the cutoff, signaling suspected outlier presence.
A length numeric vector corresponding to the data locations in
X
. Each value indicates whether the location was masked:
The data location was not masked out.
The data location was masked out, because it had at least one NA
or NaN
value.
The data location was masked out, because it was constant.
This will be a list with components:
The by
PC score matrix.
The standard deviation of each PC.
The by
PC directions matrix. Included only if
get_dirs
.
The length Q
logical vector indicating scores of high kurtosis.
Detrended components of U
. Included only if components were mean- or variance-detrended.
The length Q
logical vector indicating detrended scores of high kurtosis.
The number of PCs selected by PESEL.
The number of above-average variance PCs.
where Q
is the number of PCs selected by PESEL or of above-average variance (or the greater of the two if both were used).
If PCA was not used, all entries except nPCs_PESEL
and/or nPCs_avgvar
will not be included, depending on which
method(s) was used to select the number of PCs.
If ICA was used, this will be a list with components:
The by
source signals matrix. Included only if
get_dirs
The by
mixing matrix.
The length Q
logical vector indicating mixing scores of high kurtosis.
Detrended components of M
. Included only if components were mean- or variance-detrended.
The length Q
logical vector indicating detrended mixing scores of high kurtosis. Included only if components were mean- or variance-detrended.
Mejia, A. F., Nebel, M. B., Eloyan, A., Caffo, B. & Lindquist, M. A. PCA leverage: outlier detection for high-dimensional functional magnetic resonance imaging data. Biostatistics 18, 521-536 (2017).
Pham, D., McDonald, D., Ding, L., Nebel, M. B. & Mejia, A. Less is more: balancing noise reduction and data retention in fMRI with projection scrubbing. (2022).
library(fastICA) n_voxels = 2e3 n_timepoints = 35 X = matrix(rnorm(n_timepoints*n_voxels), ncol = n_voxels) psx = pscrub(X)
library(fastICA) n_voxels = 2e3 n_timepoints = 35 X = matrix(rnorm(n_timepoints*n_voxels), ncol = n_voxels) psx = pscrub(X)
Stabilize the center and scale of a timeseries using robust regression of DCT bases on the first and second moments.
rob_stabilize( x, center = TRUE, scale = TRUE, lmrob_method = "MM", rescale = TRUE )
rob_stabilize( x, center = TRUE, scale = TRUE, lmrob_method = "MM", rescale = TRUE )
x |
The timeseries to stabilize. |
center , scale
|
Center and scale? Default: |
lmrob_method |
The |
rescale |
After stabilizing |
the timeseries with its center and scale stabilized
Scrubbing with robust distance.
robdist( X, RD_cutoff = 4, RD_quantile = 0.99, trans = c("none", "robust-YJ", "SHASH"), bootstrap_n = 1000, bootstrap_alpha = 0.01, projection = c("ICA", "PCA"), nuisance = "DCT4", center = TRUE, scale = TRUE, comps_mean_dt = FALSE, comps_var_dt = FALSE, PESEL = TRUE, kurt_quantile = 0.99, get_dirs = FALSE, full_PCA = FALSE, get_outliers = TRUE, cutoff = 4, seed = 0, skip_dimred = FALSE, verbose = FALSE )
robdist( X, RD_cutoff = 4, RD_quantile = 0.99, trans = c("none", "robust-YJ", "SHASH"), bootstrap_n = 1000, bootstrap_alpha = 0.01, projection = c("ICA", "PCA"), nuisance = "DCT4", center = TRUE, scale = TRUE, comps_mean_dt = FALSE, comps_var_dt = FALSE, PESEL = TRUE, kurt_quantile = 0.99, get_dirs = FALSE, full_PCA = FALSE, get_outliers = TRUE, cutoff = 4, seed = 0, skip_dimred = FALSE, verbose = FALSE )
X |
Wide numeric data matrix ( |
RD_cutoff |
Default: |
RD_quantile |
Quantile cutoff...? |
trans |
Apply a transformation prior to univariate outlier detection?
Three options: |
bootstrap_n |
Use bootstrapping to estimate the robust distance null
distribution? If so, set this to the number of bootstraps. Default:
|
bootstrap_alpha |
If using bootstrap ( |
projection |
One of the following: |
nuisance |
Nuisance signals to regress from each column of Detrending is highly recommended for time-series data, especially if there are many time points or evolving circumstances affecting the data. Additionally, if kurtosis is being used to select the projection directions, trends can induce positive or negative kurtosis, contaminating the connection between high kurtosis and outlier presence. Detrending should not be used with non-time-series data because the observations are not temporally related. Additional nuisance regressors can be specified like so:
|
center , scale
|
Center the columns of the data by their medians, and scale the
columns of the data by their median absolute deviations (MADs)? Default: Note that centering and scaling occur after nuisance regression, so even if
|
comps_mean_dt , comps_var_dt
|
Stabilize the mean and variance of each
projection component's timecourse prior to computing kurtosis and leverage?
These arguments should be Slow-moving mean and variance patterns in the components will interfere with
the roles of kurtosis and leverage in identifying outliers. While
Overall, for fMRI we recommend enabling |
PESEL |
Use |
kurt_quantile |
What quantile cutoff should be used to select the
components? Default: We model each component as a length |
get_dirs |
Should the projection directions be returned? This is the
|
full_PCA |
Only applies to the PCA projection. Return the full SVD?
Default: |
get_outliers |
Should outliers be flagged based on |
cutoff |
Median leverage cutoff value. Default: |
seed |
Set a seed right before the call to |
skip_dimred |
Skip dimension reduction? Default: |
verbose |
Should occasional updates be printed? Default: |
A "robdist"
object, i.e. a list with components
...
...
...
library(fastICA) n_voxels = 2e3 n_timepoints = 35 X = matrix(rnorm(n_timepoints*n_voxels), ncol = n_voxels) rdx = robdist(X)
library(fastICA) n_voxels = 2e3 n_timepoints = 35 X = matrix(rnorm(n_timepoints*n_voxels), ncol = n_voxels) rdx = robdist(X)
Performs projection scrubbing or DVARS scrubbing, and optionally thresholds to identify artifactual time points.
scrub(X, method = c("projection", "DVARS"), ...)
scrub(X, method = c("projection", "DVARS"), ...)
X |
A |
method |
|
... |
Additional arguments to the specific scrubbing function: see
|
A list with components
A length vector or data.frame with
rows, giving the outlyingness measure(s)
Describes the outlyingness measure(s)
The outlier cutoff value(s).
A length vector or data.frame with
rows, where
TRUE
indicates suspected outlier presence.
Performs projection scrubbing or DVARS scrubbing, and optionally thresholds
to identify artifactual time points. Requires ciftiTools
and the
Connectome Workbench.
scrub_xifti( X, method = c("projection", "DVARS"), brainstructures = c("left", "right"), ... )
scrub_xifti( X, method = c("projection", "DVARS"), brainstructures = c("left", "right"), ... )
X |
Path to a CIFTI file, or a |
method |
|
brainstructures |
Character vector indicating which brain structure(s)
to use: |
... |
Additional arguments to each specific scrubbing function:
|
A list with components
A length vector or data.frame with
rows, giving the outlyingness measure(s)
Describes the outlyingness measure(s)
The outlier cutoff value(s).
A length vector or data.frame with
rows, where
TRUE
indicates suspected outlier presence.
A robust outlier detection based on modeling the data as coming from a SHASH distribution.
SHASH_out(x, maxit = 100, out_lim = 4, weight_init = NULL)
SHASH_out(x, maxit = 100, out_lim = 4, weight_init = NULL)
x |
The numeric vector in which to detect outliers. |
maxit |
The maximum number of iterations. Default: |
out_lim |
SD threshold for outlier flagging. Default: |
weight_init |
Initial weights. Default: |
A "SHASH_out"
object, i.e. a list with components
Indices of the detected outliers.
The normalized data.
Coefficients for the SHASH-to-normal transformation.
TRUE for the detected outliers for each itertation.
Last iteration number.
Logical indicating whether the convergence criteria was satisfied or not.
x <- rnorm(100) + (seq(100)/200) x[77] <- 13 SHASH_out(x)
x <- rnorm(100) + (seq(100)/200) x[77] <- 13 SHASH_out(x)
Transform SHASH-distributed data to normal-distributed data.
SHASH_to_normal(x, mu, sigma, nu, tau, inverse = FALSE)
SHASH_to_normal(x, mu, sigma, nu, tau, inverse = FALSE)
x |
Numeric vector of data to transform. |
mu |
Parameter that modulates the mean of |
sigma |
Parameter that modulates the variance of |
nu |
Parameter that modulates the skewness of |
tau |
Parameter that modulates the tailweight of |
inverse |
Transform normal data to SHASH instead? Default: |
The transformed data.
"scrub_DVARS"
objectSummary method for class "scrub_DVARS"
## S3 method for class 'scrub_DVARS' summary(object, ...) ## S3 method for class 'summary.scrub_DVARS' print(x, ...) ## S3 method for class 'scrub_DVARS' print(x, ...)
## S3 method for class 'scrub_DVARS' summary(object, ...) ## S3 method for class 'summary.scrub_DVARS' print(x, ...) ## S3 method for class 'scrub_DVARS' print(x, ...)
object |
Object of class |
... |
further arguments passed to or from other methods. |
x |
Object of class |
A plot of the scrubbing results
"scrub_FD"
objectSummary method for class "scrub_FD"
## S3 method for class 'scrub_FD' summary(object, ...) ## S3 method for class 'summary.scrub_FD' print(x, ...) ## S3 method for class 'scrub_FD' print(x, ...)
## S3 method for class 'scrub_FD' summary(object, ...) ## S3 method for class 'summary.scrub_FD' print(x, ...) ## S3 method for class 'scrub_FD' print(x, ...)
object |
Object of class |
... |
further arguments passed to or from other methods. |
x |
Object of class |
A plot of the scrubbing results
"scrub_projection"
objectSummary method for class "scrub_projection"
## S3 method for class 'scrub_projection' summary(object, ...) ## S3 method for class 'summary.scrub_projection' print(x, ...) ## S3 method for class 'scrub_projection' print(x, ...)
## S3 method for class 'scrub_projection' summary(object, ...) ## S3 method for class 'summary.scrub_projection' print(x, ...) ## S3 method for class 'scrub_projection' print(x, ...)
object |
Object of class |
... |
further arguments passed to or from other methods. |
x |
Object of class |
A plot of the scrubbing results