From 9c14ce67a1a7bb691e62ca531ee8b002cb6df14f Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Tue, 24 Jan 2023 17:17:36 +0700 Subject: [PATCH 01/14] Use `fMRItools` where possible. Move some things to that package. --- DESCRIPTION | 7 +- NAMESPACE | 26 +- R/CompCor.R | 175 ------------- R/CompCor_HCP.R | 191 -------------- R/DVARS.R | 57 +---- R/FD.R | 1 + R/artifact_images.R | 49 +--- R/carpetplot.R | 168 ------------- R/format_data.R | 328 ------------------------- R/fsl_bptf.R | 53 ---- R/leverage.R | 1 + R/nuisance_regression.R | 48 ---- R/pscrub_multi.R | 21 +- R/rob_stabilize.R | 8 +- R/rox_args_docs.R | 2 +- R/utils.R | 210 ---------------- README.Rmd | 2 +- README.md | 5 +- man/CompCor.Rd | 148 ----------- man/CompCor.noise_comps.Rd | 23 -- man/CompCor_HCP.Rd | 123 ---------- man/DVARS.Rd | 2 +- man/Mode.Rd | 19 -- man/as.matrix2.Rd | 24 -- man/carpetplot.Rd | 66 ----- man/carpetplot_stack.Rd | 45 ---- man/check_design_matrix.Rd | 18 -- man/dct_bases.Rd | 19 -- man/dct_convert.Rd | 39 --- man/detrend.Rd | 24 -- man/erode_vol.Rd | 32 --- man/fft_detrend.Rd | 21 -- man/format_data.Rd | 98 -------- man/fsl_bptf.Rd | 32 --- man/get_NIFTI_ROI_masks.Rd | 41 ---- man/hat_matrix.Rd | 18 -- man/infer_BOLD_format.Rd | 21 -- man/is_constant.Rd | 21 -- man/is_integer.Rd | 20 -- man/neighbor_counts.Rd | 22 -- man/nuisance_regression.Rd | 20 -- man/pct_sig.Rd | 26 -- man/pscrub.Rd | 2 +- man/pscrub_Params.Rd | 2 +- man/pscrub_multi.Rd | 2 +- man/read_nifti.Rd | 19 -- man/scale_med.Rd | 28 --- man/unmask_vol.Rd | 26 -- tests/testthat.R | 1 + tests/testthat/test-test_all_methods.R | 6 +- tests/testthat/test-test_scrubbing.R | 6 +- vignettes/projection_scrubbing.rmd | 6 +- 52 files changed, 53 insertions(+), 2319 deletions(-) delete mode 100644 R/CompCor.R delete mode 100644 R/CompCor_HCP.R delete mode 100644 R/carpetplot.R delete mode 100644 R/format_data.R delete mode 100644 R/fsl_bptf.R delete mode 100644 R/nuisance_regression.R delete mode 100644 R/utils.R delete mode 100644 man/CompCor.Rd delete mode 100644 man/CompCor.noise_comps.Rd delete mode 100644 man/CompCor_HCP.Rd delete mode 100644 man/Mode.Rd delete mode 100644 man/as.matrix2.Rd delete mode 100644 man/carpetplot.Rd delete mode 100644 man/carpetplot_stack.Rd delete mode 100644 man/check_design_matrix.Rd delete mode 100644 man/dct_bases.Rd delete mode 100644 man/dct_convert.Rd delete mode 100644 man/detrend.Rd delete mode 100644 man/erode_vol.Rd delete mode 100644 man/fft_detrend.Rd delete mode 100644 man/format_data.Rd delete mode 100644 man/fsl_bptf.Rd delete mode 100644 man/get_NIFTI_ROI_masks.Rd delete mode 100644 man/hat_matrix.Rd delete mode 100644 man/infer_BOLD_format.Rd delete mode 100644 man/is_constant.Rd delete mode 100644 man/is_integer.Rd delete mode 100644 man/neighbor_counts.Rd delete mode 100644 man/nuisance_regression.Rd delete mode 100644 man/pct_sig.Rd delete mode 100644 man/read_nifti.Rd delete mode 100644 man/scale_med.Rd delete mode 100644 man/unmask_vol.Rd diff --git a/DESCRIPTION b/DESCRIPTION index a01de9f..69c31b5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: fMRIscrub Type: Package Title: Scrubbing and Other Data Cleaning Routines for fMRI -Version: 0.11.0 +Version: 0.11.1 Authors@R: c( person(given = "Amanda", family = "Mejia", @@ -30,7 +30,7 @@ Description: Data-driven fMRI denoising with projection scrubbing (Pham et al (2012) ), aCompCor (anatomical Components Correction) (Muschelli et al (2014) ), detrending, and nuisance - regression. Projection scrubbing and DVARS are also applicable to other + regression. Projection scrubbing is also applicable to other outlier detection tasks involving high-dimensional data. Depends: R (>= 3.5.0) License: GPL-3 @@ -39,6 +39,7 @@ Imports: expm, MASS, e1071, + fMRItools (>= 0.2.2), pesel, robustbase, stats, @@ -61,7 +62,7 @@ Suggests: testthat (>= 3.0.0), covr Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.0 +RoxygenNote: 7.2.3 URL: https://github.com/mandymejia/fMRIscrub BugReports: https://github.com/mandymejia/fMRIscrub/issues LazyData: true diff --git a/NAMESPACE b/NAMESPACE index 12c5663..62bd41e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,40 +12,30 @@ S3method(print,summary.scrub_projection) S3method(summary,scrub_DVARS) S3method(summary,scrub_FD) S3method(summary,scrub_projection) -export(CompCor) -export(CompCor_HCP) export(DVARS) export(FD) export(artifact_images) -export(bandstop_filter) -export(carpetplot) -export(carpetplot_stack) -export(dct_bases) -export(dct_convert) -export(detrend) -export(erode_vol) -export(fsl_bptf) export(fusedPCA) export(high_kurtosis) export(leverage) -export(nuisance_regression) -export(pct_sig) export(pscrub) export(rob_stabilize) export(scrub) export(scrub_xifti) -export(unmask_vol) importFrom(MASS,mvrnorm) importFrom(e1071,kurtosis) -importFrom(grDevices,dev.off) -importFrom(grDevices,gray.colors) -importFrom(grDevices,pdf) -importFrom(grDevices,png) +importFrom(fMRItools,as.matrix_ifti) +importFrom(fMRItools,dct_bases) +importFrom(fMRItools,hat_matrix) +importFrom(fMRItools,is_constant) +importFrom(fMRItools,is_integer) +importFrom(fMRItools,nuisance_regression) +importFrom(fMRItools,scale_med) +importFrom(fMRItools,validate_design_mat) importFrom(pesel,pesel) importFrom(robustbase,rowMedians) importFrom(stats,mad) importFrom(stats,median) -importFrom(stats,mvfft) importFrom(stats,pchisq) importFrom(stats,qnorm) importFrom(stats,quantile) diff --git a/R/CompCor.R b/R/CompCor.R deleted file mode 100644 index 08c93d5..0000000 --- a/R/CompCor.R +++ /dev/null @@ -1,175 +0,0 @@ -#' CompCor: get noise components -#' -#' Get noise components for aCompCor. -#' -#' @param X_noise The noise ROIs data -#' @param center,scale Center & scale robustly -#' @param noise_nPC Number of PCs to obtain for each noise ROI -#' -#' @return A list with components X, X_noise, ROI_data, ROI_noise, noise_nPC, -#' noise_erosion, noise_comps, and noise_var. -#' -#' @keywords internal -CompCor.noise_comps <- function(X_noise, center, scale, noise_nPC){ - TOL <- 1e-8 - - N <- length(X_noise) - noise_comps <- vector("list", N); names(noise_comps) <- names(X_noise) - noise_var <- vector("list", N); names(noise_var) <- names(X_noise) - noise_vartotal <- vector("list", N); names(noise_vartotal) <- names(X_noise) - - if (length(noise_nPC) == 1) { - noise_nPC <- as.list(rep(noise_nPC, N)) - } else { - noise_nPC <- as.list(noise_nPC) - } - names(noise_nPC) <- names(X_noise) - - for (ii in 1:N) { - T_ <- nrow(X_noise[[ii]]) - if (ncol(X_noise[[ii]])==0) { next } - - # Transpose. - X_noise[[ii]] <- t(X_noise[[ii]]) - # Center. - if (center) { X_noise[[ii]] <- X_noise[[ii]] - c(rowMedians(X_noise[[ii]], na.rm=TRUE)) } - # Compute MADs. - mad <- 1.4826 * rowMedians(abs(X_noise[[ii]]), na.rm=TRUE) - X_constant <- mad < TOL - if (any(X_constant)) { - if (all(X_constant)) { - stop(paste0("All data locations in noise ROI ", ii, " are zero-variance.\n")) - } else { - warning(paste0("Warning: removing ", sum(X_constant), - " constant data locations (out of ", length(X_constant), - ") in noise ROI ", ii, - ".\n")) - } - } - mad <- mad[!X_constant]; X_noise[[ii]] <- X_noise[[ii]][!X_constant,] - # Scale. - if (scale) { X_noise[[ii]] <- X_noise[[ii]]/c(mad) } - # Revert transpose. - X_noise[[ii]] <- t(X_noise[[ii]]) - if (ncol(X_noise[[ii]])==0) { next } - - # Compute the PC scores. - x <- svd(tcrossprod(X_noise[[ii]])) - noise_var[[ii]] <- x$d - noise_vartotal[[ii]] <- sum(noise_var[[ii]]) - if (noise_nPC[[ii]] < 1) { - # Use enough PCs to explain the desired proportion of variance. - noise_nPC[[ii]] <- min( - length(x$d), - sum(cumsum(noise_var[[ii]]) < noise_vartotal[[ii]]*noise_nPC[[ii]]) + 1 - ) - } - noise_comps[[ii]] <- x$u[,seq(noise_nPC[[ii]]),drop=FALSE] - noise_var[[ii]] <- noise_var[[ii]][seq(noise_nPC[[ii]])] - } - - list(noise_comps=noise_comps, noise_var=noise_var, noise_vartotal=noise_vartotal) -} - -#' Anatomical CompCor -#' -#' The aCompCor algorithm for denoising fMRI data using noise ROIs data -#' -#' First, the principal components (PCs) of each noise region of interest (ROI) -#' are calculated. For each ROI, voxels are centered and scaled -#' (can be disabled with the arguments \code{center} and \code{scale}), -#' and then the PCs are calculated via the singular value decomposition. -#' -#' Next, aCompCor is performed to remove the shared variation between the noise ROI -#' PCs and each location in the data. This is accomplished by a nuisance regression -#' using a design matrix with the noise ROI PCs, any additional regressors specified -#' by \code{nuisance}, and an intercept term. (To detrend the data and perform aCompCor -#' in the same regression, \code{nuisance} can be set to DCT bases obtained with -#' the function \code{\link{dct_bases}}.) -#' -#' @inheritParams data_CompCor_Params -#' @inheritParams noise_Params -#' @param center,scale Center the columns of the noise ROI data by their medians, -#' and scale by their MADs? Default: \code{TRUE} for both. Note that this argument -#' affects the noise ROI data and not the data that is being cleaned with aCompCor. -#' Centering and scaling of the data being cleaned can be done after this function call. -#' @param nuisance Nuisance signals to regress from each data column in addition -#' to the noise ROI PCs. Should be a \eqn{T} by \eqn{N} numeric matrix where -#' \eqn{N} represents the number of nuisance signals. To not perform any nuisance -#' regression set this argument to \code{NULL}, \code{0}, or \code{FALSE}. -#' Default: \code{NULL}. -#' @return A list with entries \code{"data"}, \code{"noise"}, and potentially -#' \code{"ROI_data"}. -#' -#' The entry \code{"data"} will be a \code{V x T} matrix where each row corresponds to a -#' data location (if it was originally an array, the locations will be voxels in spatial -#' order). Each row will be a time series with each noise PC regressed from it. This entry -#' will be \code{NULL} if there was no data. -#' -#' The entry \code{"noise"} is a list of noise PC scores, their corresponding variance, -#' and their ROI mask, for each noise ROI. -#' -#' If the data ROI is not all \code{TRUE}, the entry \code{"ROI_data"} will have -#' the ROI mask for the data. -#' -#' @importFrom robustbase rowMedians -#' -#' @section References: -#' \itemize{ -#' \item{Behzadi, Y., Restom, K., Liau, J. & Liu, T. T. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage 37, 90-101 (2007).} -#' \item{Muschelli, J. et al. Reduction of motion-related artifacts in resting state fMRI using aCompCor. NeuroImage 96, 22-35 (2014).} -#' } -#' -#' @export -#' @seealso CompCor_HCP -CompCor <- function( - X, ROI_data="infer", ROI_noise=NULL, - noise_nPC=5, noise_erosion=NULL, - center=TRUE, scale=TRUE, nuisance=NULL - ){ - - if (is.null(ROI_noise)) {stop("`ROI_noise` must be provided. (It is not inferred.)")} - - out1 <- format_data( - X=X, ROI_data=ROI_data, ROI_noise=ROI_noise, - noise_nPC=noise_nPC, noise_erosion=noise_erosion - ) - - out2 <- CompCor.noise_comps( - X_noise=out1$X_noise, center=center, scale=scale, noise_nPC=out1$noise_nPC - ) - - T_ <- nrow(out1$X) - - if (any(out1$ROI_data)) { - # Perform nuisance regression. - design <- do.call(cbind, out2$noise_comps) - if (!(is.null(nuisance) || isFALSE(nuisance) || identical(nuisance, 0))) { - if (nrow(nuisance) != T_) { stop("`nuisance` does not have the same number of timepoints as `X`.") } - nuisance <- check_design_matrix(nuisance) - design <- cbind(design, nuisance) - } - design <- check_design_matrix(cbind(1, design), T_) - out1$X <- nuisance_regression(out1$X, design) - - # # Normalize data. - # out1$X <- t(out1$X) - # if (center) { out1$X <- out1$X - c(rowMedians(out1$X, na.rm=TRUE)) } - # if (scale) { - # mad <- 1.4826 * rowMedians(abs(out1$X), na.rm=TRUE) - # mad_inv <- ifelse(mad < 1e-8, 0, 1/mad) - # out1$X <- out1$X * mad_inv - # } - # out1$X <- t(out1$X) - } else { - out1["X"] <- list(NULL) - } - - z <- list( - data = out1$X, - noise = list(PCs=out2$noise_comps, var=out2$noise_var, ROI_noise=out1$ROI_noise) - ) - if (!all(out1$ROI_data)) { z$ROI_data <- out1$ROI_data } - class(z) <- "CompCor" - z -} \ No newline at end of file diff --git a/R/CompCor_HCP.R b/R/CompCor_HCP.R deleted file mode 100644 index 0511777..0000000 --- a/R/CompCor_HCP.R +++ /dev/null @@ -1,191 +0,0 @@ -#' Get NIFTI ROI masks -#' -#' Get NIFTI ROI masks -#' -#' @param nii_labels \eqn{I} by \eqn{J} by \eqn{K} -#' NIFTI object or array (or file path to the NIFTI) which -#' contains the corresponding labels to each voxel in \code{nii}. Values should -#' be according to this table: -#' https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/AnatomicalROI/FreeSurferColorLUT . -#' In the HCP, the corresponding file is "ROIs/Atlas_wmparc.2.nii.gz". -#' @param ROI_noise A list of numeric vectors. Each entry should represent labels -#' in \code{nii_labels} belonging to a single noise ROI, named by that entry's -#' name. Or, this can be a character vector of at least one of the following: -#' \code{"wm_cort"} (cortical white matter), \code{"wm_cblm"} (cerebellar white -#' matter), \code{"csf"} (cerebrospinal fluid). In the latter case, these labels -#' will be used: -#' -#' \describe{ -#' \item{\code{"wm_cort"}}{\code{c(3000:4035, 5001, 5002)}} -#' \item{\code{"wm_cblm"}}{\code{c(7, 46)}} -#' \item{\code{"csf"}}{\code{c(4, 5, 14, 15, 24, 31, 43, 44, 63, 250, 251, 252, 253, 254, 255))}} -#' } -#' -#' These default ROIs are based on this forum post: -#' https://www.mail-archive.com/hcp-users@humanconnectome.org/msg00931.html -#' -#' Default: \code{c("wm_cort", "csf")} -#' @return The ROIs -#' @keywords internal -#' -get_NIFTI_ROI_masks <- function(nii_labels, ROI_noise=c("wm_cort", "csf")){ - # Read NIFTI labels file. - if (is.character(nii_labels)) { - nii_labels <- read_nifti(nii_labels) - } - stopifnot(length(dim(nii_labels))==3) - - # Get ROI masks. - ROI_noise.default <- list( - wm_cort = c(3000:4035, 5001, 5002), - wm_cblm = c(7, 46), - csf = c(4, 5, 14, 15, 24, 31, 43, 44, 63, 250, 251, 252, 253, 254, 255) - ) - if (is.null(ROI_noise)) { - ROI_noise <- ROI_noise.default[c("wm_cort", "csf")] - } else if (is.character(ROI_noise)) { - stopifnot(all(ROI_noise %in% names(ROI_noise.default))) - ROI_noise <- ROI_noise.default[unique(ROI_noise)] - } else if (is.list(ROI_noise)) { - ROI_noise <- ROI_noise - stopifnot(is.numeric(do.call(c, ROI_noise))) - } else { - stop("`ROI_noise` argument is not in a recognized form.") - } - lapply( - ROI_noise, - function(x){array(nii_labels %in% x, dim=dim(nii_labels))} - ) -} - -#' Anatomical CompCor for HCP NIFTI and CIFTI data -#' -#' Wrapper to \code{\link{CompCor}} for HCP-format data. Can be used to clean -#' the surface-based CIFTI data with aCompCor using the noise PCs and ROIs -#' calculated from the NIFTI fMRI data and NIFTI mask. Can also be used to just -#' obtain the noise PCs and ROIs without performing aCompCor, if the CIFTI -#' data is not provided. -#' -#' @inheritParams get_NIFTI_ROI_masks -#' @param nii \eqn{I} by \eqn{J} by \eqn{K} by \eqn{T} -#' NIFTI object or array (or file path to the NIFTI) which contains -#' whole-brain data, including the noise ROIs. In the HCP, the corresponding -#' file is e.g. "../Results/rfMRI_REST1_LR/rfMRI_REST1_LR.nii.gz" -#' @param noise_nPC The number of principal components to compute for each noise -#' ROI. Alternatively, values between 0 and 1, in which case they will -#' represent the minimum proportion of variance explained by the PCs used for -#' each noise ROI. The smallest number of PCs will be used to achieve this -#' proportion of variance explained. -#' -#' Should be a list or numeric vector with the same length as \code{ROI_noise}. -#' It will be matched to each ROI based on the name of each entry, or if the -#' names are missing, the order of entries. If it is an unnamed vector, its -#' elements will be recycled. Default: \code{5} (compute the top 5 PCs for -#' each noise ROI). -#' @param noise_erosion The number of voxel layers to erode the noise ROIs by. -#' Should be a list or numeric vector with the same length as \code{ROI_noise}. -#' It will be matched to each ROI based on the name of each entry, or if the -#' names are missing, the order of entries. If it is an unnamed vector, its -#' elements will be recycled. Default: \code{NULL}, which will use a value of -#' 0 (do not erode the noise ROIs). -#' @param brainstructures Choose among "left", "right", and "subcortical". -#' Default: \code{c("left", "right")} (cortical data only) -#' @param idx A numeric vector indicating the timepoints to use, or -#' \code{NULL} (default) to use all idx. (Indexing begins with 1, so the -#' first timepoint has index 1 and the last has the same index as the length of -#' the scan.) -#' @param cii \code{"xifti"} (or file path to the CIFTI) from which the noise -#' ROI components will be regressed. In the HCP, the corresponding file is e.g. -#' "../Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll.dtseries.nii". If not -#' provided, only the noise components will be returned (no data will be cleaned). -#' @param center,scale Center the columns of the data by median, and scale the -#' columns of the data by MAD? Default: \code{TRUE} for both. Affects both -#' \code{X} and the noise data. \code{center} also applies to \code{nuisance_too} -#' so if it is \code{FALSE}, \code{nuisance_too} must already be centered. -#' @param DCT Add DCT bases to the nuisance regression? Use an integer to -#' indicate the number of cosine bases. Use \code{0} (default) to forgo detrending. -#' -#' The data must be centered, either before input or with \code{center}. -#' @param nuisance_too A matrix of nuisance signals to add to the nuisance -#' regression. Should have \eqn{T} rows. \code{NULL} to not add additional -#' nuisance regressors (default). -#' @param verbose Should occasional updates be printed? Default: \code{FALSE}. -#' -#' @return The noise components, and if \code{cii} is provided, the cleaned -#' surface-based data as a \code{"xifti"} object. -#' -#' @section References: -#' \itemize{ -#' \item{Behzadi, Y., Restom, K., Liau, J. & Liu, T. T. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage 37, 90-101 (2007).} -#' \item{Muschelli, J. et al. Reduction of motion-related artifacts in resting state fMRI using aCompCor. NeuroImage 96, 22-35 (2014).} -#' } -#' -#' @export -#' -#' @seealso CompCor -CompCor_HCP <- function( - nii, nii_labels, - ROI_noise=c("wm_cort", "csf"), noise_nPC=5, noise_erosion=NULL, - idx=NULL, cii=NULL, brainstructures=c("left", "right"), - center = TRUE, scale = TRUE, DCT = 0, nuisance_too = NULL, - verbose=FALSE){ - - # Get NIFTI data. - if (is.character(nii)) { - if (verbose) { cat("Reading data NIFTI.\n") } - nii <- read_nifti(nii) - } - stopifnot(length(dim(nii))==4) - - # Get NIFTI ROI masks. - if (verbose) { cat("Reading labels NIFTI.\n") } - ROI_noise <- get_NIFTI_ROI_masks(nii_labels, ROI_noise) - stopifnot(all( dim(ROI_noise[[1]]) == dim(nii)[1:3] )) - - # Drop idx in NIFTI. - T_ <- dim(nii)[4] - if (!is.null(idx)) { - stopifnot(all(idx %in% seq(T_))) - nii <- nii[,,,idx] - } - T_ <- dim(nii)[4] - if (T_ < 10) { - warning("There are very few `idx`.\n") - } - - if (verbose) { cat("Computing noise components.\n") } - out <- CompCor( - nii, ROI_data=NULL, ROI_noise=ROI_noise, - noise_erosion=noise_erosion, noise_nPC=noise_nPC, - center=center, scale=scale, - # We do nuisance regression with aCompCor, not prior. - nuisance=NULL - ) - - if (is.null(cii)) { return(out$noise) } - - if (is.character(cii)) { - if (requireNamespace("ciftiTools", quietly = TRUE)) { - cii <- ciftiTools::read_cifti(cii, brainstructures=brainstructures, verbose=verbose) - } else { - stop("Package `ciftiTools` required to read the CIFTI file. Please install it.") - } - } - stopifnot(all(names(cii) == c("data", "surf", "meta"))) - - # Drop idx. - if (!is.null(idx)) { cii <- ciftiTools::select_xifti(cii, idx) } - - # Make design matrix. - design <- do.call(cbind, out$noise$PCs) - if (DCT > 0) { design <- cbind(design, dct_bases(T_, DCT)) } - if (!is.null(nuisance_too)) { - nuisance_too <- check_design_matrix(nuisance_too) - design <- cbind(design, nuisance_too) - } - design <- check_design_matrix(cbind(1, design)) - - # Return data in \code{"xifti"} format. - out$data <- ciftiTools::newdata_xifti(cii, nuisance_regression(cii, design)) - out -} \ No newline at end of file diff --git a/R/DVARS.R b/R/DVARS.R index aba8a2d..c1e5456 100644 --- a/R/DVARS.R +++ b/R/DVARS.R @@ -22,58 +22,6 @@ sd_hIQR <- function(x, d=1){ return(as.numeric(out)) } -#' Mode of data vector -#' -#' Get mode of a data vector. But use the median instead of the mode if all -#' data values are unique. -#' -#' @param x The data vector -#' -#' @return The mode -#' -#' @keywords internal -#' -Mode <- function(x) { - q <- unique(x) - # Use median instead of the mode if all data values are unique. - if (length(q) == length(x)) { return(median(x)) } - q[which.max(tabulate(match(x, q)))] -} - -#' Convert data values to percent signal. -#' -#' Convert data values to percent signal. -#' -#' @param X a \eqn{T} by \eqn{N} numeric matrix. The columns will be normalized to -#' percent signal. -#' @param center A function that computes the center of a numeric vector. -#' Default: \code{median}. Other common options include \code{mean} and -#' \code{mode}. -#' @param by Should the center be measured individually for each \code{"column"} -#' (default), or should the center be the same across \code{"all"} columns? -#' -#' @return \code{X} with its columns normalized to percent signal. (A value of -#' 85 will represent a -15% signal change.) -#' -#' @export -pct_sig <- function(X, center=median, by=c("column", "all")){ - stopifnot(is.numeric(X)) - stopifnot(length(dim(X))==2) - stopifnot(is.function(center)) - by <- match.arg(by, c("column", "all")) - - T_ <- nrow(X); N_ <- ncol(X) - X <- t(X) - - if (by=="column") { - m <- apply(X, 1, center) - } else { - m <- center(as.numeric(X)) - } - - t(X / m * 100) -} - #' DVARS #' #' Computes the DSE decomposition and DVARS-related statistics. Based on code @@ -108,6 +56,7 @@ pct_sig <- function(X, center=median, by=c("column", "all")){ #' } #' @export #' @importFrom stats median pchisq qnorm +#' @importFrom fMRItools as.matrix_ifti #' #' @section References: #' \itemize{ @@ -117,12 +66,12 @@ pct_sig <- function(X, center=median, by=c("column", "all")){ DVARS <- function( X, normalize=TRUE, cutoff_DPD=5, - cutoff_ZD=qnorm(1 - .05 / nrow(as.matrix2(X))), + cutoff_ZD=qnorm(1 - .05 / nrow(as.matrix_ifti(X))), verbose=FALSE){ cutoff_DVARS <- NULL - X <- as.matrix2(X, verbose=verbose) + X <- as.matrix_ifti(X, verbose=verbose) T_ <- nrow(X); N_ <- ncol(X) cutoff <- list(DVARS=cutoff_DVARS, DPD=cutoff_DPD, ZD=cutoff_ZD) diff --git a/R/FD.R b/R/FD.R index 57ff28b..75a8300 100644 --- a/R/FD.R +++ b/R/FD.R @@ -47,6 +47,7 @@ #' } #' #' @importFrom utils read.table +#' @importFrom fMRItools nuisance_regression dct_bases #' @export #' #' @section References: diff --git a/R/artifact_images.R b/R/artifact_images.R index b9670ed..9576d78 100644 --- a/R/artifact_images.R +++ b/R/artifact_images.R @@ -22,6 +22,7 @@ #' matrix give the index of the top component ("V" in PCA and "S" in ICA) at #' each timepoint. #' +#' @importFrom fMRItools scale_med #' @export artifact_images <- function(psx, idx=NULL, use_dt=TRUE){ @@ -55,7 +56,7 @@ artifact_images <- function(psx, idx=NULL, use_dt=TRUE){ } V <- psx$fusedPCA$V } else if ("ICA" %in% names(psx)) { - U <- scale_med(psx$ICA$M) + U <- fMRItools::scale_med(psx$ICA$M) if (!("S" %in% names (psx$ICA))) { stop("No directions. Run `pscrub` again with `get_dirs=TRUE`.") } @@ -90,50 +91,4 @@ artifact_images <- function(psx, idx=NULL, use_dt=TRUE){ colnames(lev_imgs$top) <- NULL lev_imgs -} - -#' Undo a volumetric mask -#' -#' Un-applies a mask to vectorized data to yield its volumetric representation. -#' The mask and data should have compatible dimensions: the number of rows in -#' \code{dat} should equal the number of locations within the \code{mask}. -#' -#' @param dat Data matrix with locations along the rows and measurements along -#' the columns. If only one set of measurements were made, this may be a -#' vector. -#' @param mask Volumetric binary mask. \code{TRUE} indicates voxels inside the -#' mask. -#' @param fill The value for locations outside the mask. Default: \code{NA}. -#' -#' @return The 3D or 4D unflattened volume array -#' -#' @export -#' -unmask_vol <- function(dat, mask, fill=NA) { - - # Check that dat is a vector or matrix. - if (is.vector(dat) || is.factor(dat)) { dat <- matrix(dat, ncol=1) } - stopifnot(length(dim(dat)) == 2) - - # Check that mask is numeric {0, 1} or logical, and is 3D. - if (is.numeric(mask)) { - mask_vals <- unique(as.vector(mask)) - stopifnot(length(mask_vals) <= 2) - stopifnot(all(mask_vals %in% c(0,1))) - mask <- array(as.logical(mask), dim=dim(mask)) - } - stopifnot(length(dim(mask)) == 3) - - # Other checks. - stopifnot(is.vector(fill) && length(fill)==1) - stopifnot(sum(mask) == nrow(dat)) - - # Make volume and fill. - vol <- array(fill, dim=c(dim(mask), ncol(dat))) - for (ii in seq_len(ncol(dat))) { - vol[,,,ii][mask] <- dat[,ii] - } - if (ncol(dat)==1) { vol <- vol[,,,1] } - - vol } \ No newline at end of file diff --git a/R/carpetplot.R b/R/carpetplot.R deleted file mode 100644 index 69bf367..0000000 --- a/R/carpetplot.R +++ /dev/null @@ -1,168 +0,0 @@ -#' Carpetplot -#' -#' Plot a matrix with \code{graphics::image}. For fMRI data, this is the "carpetplot" -#' or grayplot coined by (Power, 2017). The \code{graphics} package is required. -#' -#' @param x The \eqn{T x V} numeric data matrix, or a \code{"xifti"} object. -#' In the plot, the \eqn{T} index will increase from left to right, and the -#' \eqn{V} will increase from top to bottom. -#' @param qcut Sets blackpoint at the \code{qcut} quantile, and the -#' whitepoint at the \code{1-qcut} quantile. Default: \code{.1}. This is -#' equivalent to setting the color range between the 10% and 90% quantiles. -#' The quantiles are computed across the entire data matrix after any -#' centering or scaling. -#' -#' Must be between 0 and .49. If \code{0} or \code{NULL} (default), do not -#' clamp the data values. -#' @param fname A \code{.pdf} (highly recommended) or \code{.png} file path -#' to write the carpetplot to. If \code{NULL} (default), return the plot directly -#' instead of writing a file. -#' @param center,scale Center and scale the data? If \code{x} is fMRI data -#' which has not otherwise been centered or scaled, it is recommended to center -#' but not scale it (default). -#' @param colors \code{"gray255"} (default) will use a grayscale color ramp -#' from black to white. Otherwise, this should be a character vector of -#' color names to use. -#' -#' Colors will be assigned from the lowest to the highest data value, after -#' any clamping of the data values by \code{qcut}. -#' @param sortSub If \code{x} is a \code{"xifti"} object with subcortical data, -#' should the voxels be sorted by structure alphabetically? Default: \code{TRUE}. -#' @param ... Additional arguments to \code{pdf} or \code{png}, such as width -#' and height. -#' -#' @importFrom grDevices dev.off pdf png gray.colors -#' -#' @return The image or \code{NULL, invisibly} if a file was written. -#' -#' @section References: -#' \itemize{ -#' \item{Power, J. D. A simple but useful way to assess fMRI scan qualities. NeuroImage 154, 150-158 (2017).} -#' } -#' -#' @export -#' -carpetplot <- function( - x, qcut=.1, fname=NULL, center=TRUE, scale=FALSE, colors="gray255", sortSub=TRUE, ...){ - - if (!requireNamespace("graphics", quietly = TRUE)) { - stop("Package \"graphics\" needed since `svd` failed. Please install it.", call. = FALSE) - } - - # Get T x V matrix. - x <- as.matrix2(x, sortSub=sortSub) - - # Center or scale. - if (center | scale) { x <- scale(x, center=center, scale=scale) } - - # Orient correctly for `image` - x <- x[,seq(ncol(x), 1)] - - # Clamp to quantile cutoffs. - if (!is.null(qcut)) { - qcut <- as.numeric(qcut) - if (qcut > 1e-8) { - stopifnot(qcut < .49) - x[x < quantile(x, qcut)] <- quantile(x, qcut) - x[x > quantile(x, 1-qcut)] <- quantile(x, 1-qcut) - } - } - - # Check colors. - if (identical(colors, "gray255")) { - colors <- gray.colors(255, start=0, end=1) - } - - # Check file name - if (!is.null(fname)) { - fname <- as.character(fname) - ftype <- ifelse(endsWith(fname, ".png"), "png", "pdf") - if (ftype == "pdf" && !endsWith(fname, ".pdf")) { - fname <- paste0(fname, ".pdf") - } - - switch(ftype, png=png, pdf=pdf)(fname, ...) - } - - graphics::image( - x, - axes=FALSE, frame.plot = FALSE, xlab="", ylab="", ann=FALSE, - useRaster=TRUE, col=colors - ) - - if (!is.null(fname)) { - dev.off() - } - - invisible(NULL) -} - -#' Stacked carpetplot -#' -#' Stacks carpetplots on top of one anotherr by rbinding the matrices. -#' -#' @param x_list List of data matrices -#' @param center,scale Center and scale the data? If \code{x} is fMRI data -#' which has not otherwise been centered or scaled, it is recommended to center -#' but not scale it (default). -#' @param qcut Sets blackpoint at the \code{qcut} quantile, and the -#' whitepoint at the \code{1-qcut} quantile. Default: \code{.1}. This is -#' equivalent to setting the color range between the 10% and 90% quantiles. -#' The quantiles are computed across the entire data matrix after any -#' centering or scaling. -#' -#' Must be between 0 and .49. If \code{0} or \code{NULL} (default), do not -#' clamp the data values. -#' @param match_scale Match the scales of the carpetplots? Default: \code{TRUE}. -#' @param nsep Equivalent number of data locations for size of gap between -#' carpetplots. Default: zero (no gap). -#' @param ... Additional arguments to \code{carpetplot} -#' @return \code{NULL}, invisibly -#' @export -carpetplot_stack <- function(x_list, center=TRUE, scale=FALSE, qcut=.1, match_scale=TRUE, nsep=0, ...){ - - # [TO DO] check args - args <- list(...) - sortSub <- ifelse("sortSub" %in% names(args), args$sortSub, TRUE) - - x_list <- lapply(x_list, as.matrix2, sortSub=sortSub) - - # Center or scale. - if (center | scale) { - x_list <- lapply(x_list, scale, center=center, scale=scale) - } - - # Clamp to quantile cutoffs. - if (!is.null(qcut)) { - qcut <- as.numeric(qcut) - if (qcut > 1e-8) { - stopifnot(qcut < .49) - for (ii in seq(length(x_list))) { - x_list[[ii]][x_list[[ii]] < quantile(x_list[[ii]], qcut)] <- quantile(x_list[[ii]], qcut) - x_list[[ii]][x_list[[ii]] > quantile(x_list[[ii]], 1-qcut)] <- quantile(x_list[[ii]], 1-qcut) - } - } - } - - if (match_scale) { - scale_cp <- function(x){ q <- min(x[]); x[] <- (x[]-q)/(max(x[]) - q); x } - x_list <- lapply(x_list, scale_cp) - } - xmax <- max(vapply(x_list, max, 0, na.rm=TRUE), na.rm=TRUE) - - if (nsep > 0) { - x2 <- vector("list", length(x_list)*2 - 1) - for (ii in seq(length(x2))) { - if (ii %% 2 == 1) { - x2[[ii]] <- x_list[[ceiling(ii/2)]] - } else { - x2[[ii]] <- matrix(xmax, nrow=nrow(x_list[[1]]), ncol=nsep) - } - } - x2 <- do.call(cbind, x2) - } else { - x2 <- do.call(cbind, x_list) - } - - carpetplot(x2, center=FALSE, scale=FALSE, qcut=0, ...) -} \ No newline at end of file diff --git a/R/format_data.R b/R/format_data.R deleted file mode 100644 index 8a75851..0000000 --- a/R/format_data.R +++ /dev/null @@ -1,328 +0,0 @@ -#' Count each voxel's neighbors with value \code{TRUE} -#' -#' For each voxel in a 3D logical array, count the number of immediate neighbors -#' with value \code{TRUE}. -#' -#' @param arr The 3D logical array. -#' @param pad Pad value for edge. -#' -#' @return An array with the same dimensions as \code{arr}. Each voxel value -#' will be the number of its immediate neighbors (0 to 6) which are \code{TRUE}. -#' -#' @keywords internal -neighbor_counts <- function(arr, pad=FALSE){ - stopifnot(length(dim(arr)) == 3) - stopifnot(all(unique(arr) %in% c(TRUE, FALSE))) - arrPad <- ciftiTools:::pad_vol(arr, matrix(1, 3, 2), fill=pad) - dPad <- dim(arrPad) - # Look up, down, left, right, forward, behind (not diagonally) - arrPad[1:(dPad[1]-2),2:(dPad[2]-1),2:(dPad[3]-1)] + - arrPad[3:(dPad[1]),2:(dPad[2]-1),2:(dPad[3]-1)] + - arrPad[2:(dPad[1]-1),1:(dPad[2]-2),2:(dPad[3]-1)] + - arrPad[2:(dPad[1]-1),3:(dPad[2]),2:(dPad[3]-1)] + - arrPad[2:(dPad[1]-1),2:(dPad[2]-1),1:(dPad[3]-2)] + - arrPad[2:(dPad[1]-1),2:(dPad[2]-1),3:(dPad[3])] -} - -#' Erode volumetric mask -#' -#' Erode a volumetric mask by a certain number of voxel layers. For each layer, -#' any in-mask voxel adjacent to at least one out-of-mask voxel is removed -#' from the mask. -#' -#' Diagonal voxels are not considered adjacent, i.e. the voxel at (0,0,0) is not -#' adjacent to the voxel at (1,1,0) or (1,1,1), although it is adjacent to -#' (1,0,0). -#' -#' @param vol The volume to erode. Out-of-mask voxels should be indicated by a -#' value in \code{out_of_mask_val}. -#' @param n_erosion The number of layers to erode the mask by. -#' @param out_of_mask_val A voxel is not included in the mask if and only if its -#' value is in this vector. The first value in this vector will be used to replace -#' the eroded voxels. Default: \code{NA}. -#' -#' @return The eroded \code{vol}. It is the same as \code{vol} but with eroded -#' voxels replaced with the value \code{out_of_mask_val[1]}. -#' -#' @export -erode_vol <- function(vol, n_erosion=1, out_of_mask_val=NA){ - stopifnot(is_integer(n_erosion, nneg=TRUE)) - stopifnot(length(dim(vol)) == 3) - if (n_erosion==0) { return(vol) } - mask <- !array(vol %in% out_of_mask_val, dim=dim(vol)) - for (ii in seq(n_erosion)) { - to_erode <- mask & neighbor_counts(mask, pad=TRUE) < 6 - mask[to_erode] <- FALSE - vol[to_erode] <- out_of_mask_val[1] - } - vol -} - -#' Format data for pscrub and CompCor -#' -#' @inheritParams data_CompCor_Params -#' @param noise_nPC The number of principal components to compute for each noise -#' ROI. Alternatively, values between 0 and 1, in which case they will -#' represent the minimum proportion of variance explained by the PCs used for -#' each noise ROI. The smallest number of PCs will be used to achieve this -#' proportion of variance explained. -#' -#' Should be a list or numeric vector with the same length as \code{ROI_noise}. -#' It will be matched to each ROI based on the name of each entry, or if the -#' names are missing, the order of entries. If it is an unnamed vector, its -#' elements will be recycled. Default: \code{5} (compute the top 5 PCs for -#' each noise ROI). -#' @param noise_erosion The number of voxel layers to erode the noise ROIs by. -#' Should be a list or numeric vector with the same length as \code{ROI_noise}. -#' It will be matched to each ROI based on the name of each entry, or if the -#' names are missing, the order of entries. If it is an unnamed vector, its -#' elements will be recycled. Default: \code{NULL}, which will use a value of -#' 0 (do not erode the noise ROIs). -#' -#' @return A list with components "X", "X_noise", "ROI_data", and "ROI_noise" -#' -#' @keywords internal -format_data <- function(X, ROI_data="infer", ROI_noise=NULL, noise_nPC=5, noise_erosion=NULL){ - - # TO DO: explain how to use ROI_data - - # if X is a file, read it. - if (is.character(X)) { - if (endsWith(X, ".dtseries.nii") | endsWith(X, ".dscalar.nii")) { - if (!requireNamespace("ciftiTools", quietly = TRUE)) { - stop("Package \"ciftiTools\" needed to read `X`. Please install it", call. = FALSE) - } - if (identical(ROI_data, "infer")) { ROI_data <- "all" } - X <- ciftiTools::read_cifti(X, brainstructures=ROI_data) - } else if (endsWith(X, ".nii") | endsWith(X, ".nii.gz")) { - X <- read_nifti(X) - } else { - stop("The argument `X`, '", X, "', does not look like a CIFTI or NIFTI file name.") - } - } - - # X must be a matrix, array, or "xifti" - if (is.matrix(X)) { - T_ <- nrow(X); V_ <- ncol(X) - if (T_ > V_) { - warning( - "Data matrix has more rows than columns. Check that observations\ - are in rows and variables are in columns." - ) - } - X_type <- "vector" - } else if (is.array(X)) { - if (length(dim(X))==3) { X <- array(X, dim=c(dim(X)[1:2], 1, dim(X)[3])) } - stopifnot(length(dim(X))==4) - T_ <- dim(X)[4] - X_type <- "volume" - } else if (inherits(X, "xifti")) { - xifti_meta <- X$meta - X <- t(as.matrix(X)) - T_ <- nrow(X); V_ <- ncol(X) - X_type <- "xifti" - } else { - stop("`X` must be a matrix, array, NIFTI, path to a NIFTI, CIFTI, or path to a CIFTI.") - } - - if (T_ < 2) { stop("There are less than two timepoints.") } - - # ROI_data - if (is.null(ROI_data)) { - if (X_type == "volume") { - ROI_data <- X[,,,1] > Inf - } else { - ROI_data <- rep(FALSE, V_) - } - ROI_data_infer <- FALSE - - } else { - ROI_data_infer <- identical(ROI_data, "infer") - if (X_type == "vector") { - if (ROI_data_infer) { ROI_data <- rep(TRUE, V_) } - ROI_data <- as.vector(ROI_data) - if (length(ROI_data) != V_) { - stop("The `ROI_data` must be a logical vector with the same length as columns in the data.") - } - ROI_data <- as.logical(ROI_data) - } else if (X_type == "volume") { - if (ROI_data_infer) { ROI_data <- c(0, NA, NaN) } - if (is.character(ROI_data)) { - ROI_data <- read_nifti(ROI_data) - } - if (is.vector(ROI_data)) { - if (length(ROI_data) != length(unique(ROI_data))) { - stop( - "`X` is a volume, and `ROI_data` is a vector, so `ROI_data` should be\ - values which out-of-mask voxels take on. But, the values of `ROI_data`\ - are not unique." - ) - } - ROI_data <- apply(array(X %in% ROI_data, dim=dim(X)), 1:3, function(x){!all(x)}) - } else if (is.array(ROI_data)) { - if (length(dim(ROI_data))==2) { ROI_data <- array(ROI_data, dim=c(dim(ROI_data), 1)) } - if (all(dim(ROI_data) != dim(X)[1:3])) { - stop("The `ROI_data` must have the same dimensions as the first three dimensions of `X`.") - } - mode(ROI_data) <- "logical" - } - } else if (X_type == "xifti") { - ROI_data <- rep(TRUE, V_) - } else { stop("Internal error: unrecognized `X_type`") } - if (sum(ROI_data) == 0) { warning("The data ROI was empty.\n") } - } - - # ROI_noise - if (!is.null(ROI_noise)) { #&& length(ROI_noise)>0) { - if (!is.list(ROI_noise)) { ROI_noise <- list(Noise1=ROI_noise) } - if (is.null(names(ROI_noise))) { names(ROI_noise) <- paste0("Noise", 1:length(ROI_noise)) } - if (length(names(ROI_noise)) != length(unique(names(ROI_noise)))) { - stop("The `ROI_noise` names must be unique.") - } - stopifnot(!any(names(ROI_noise) == "data")) - - # noise_nPC - if (length(noise_nPC) == 1) { - noise_nPC <- as.list(rep(noise_nPC, length(ROI_noise))) - } else { - noise_nPC <- as.list(noise_nPC) - } - names(noise_nPC) <- names(ROI_noise) - - if (is.null(names(noise_nPC))) { - noise_nPC <- noise_nPC[rep(1:length(noise_nPC), length(ROI_noise))[1:length(ROI_noise)]] - names(noise_nPC) <- names(ROI_noise) - } - else { - if (all(sort(names(noise_nPC)) != sort(names(ROI_noise)))) { - stop("The names of `noise_nPC` do not match those of `ROI_noise`.") - } - } - noise_nPC <- noise_nPC[names(ROI_noise)] - - # noise_erosion - if (X_type == "volume") { - if (is.null(noise_erosion)) { - noise_erosion = 0 - } else { - if (!all(noise_erosion==0) & !any(sapply(ROI_noise, is.array))) { - warning("`noise_erosion` was provided, but there are no array/NIFTI noise ROIs to erode.\n") - } - } - noise_erosion <- as.list(noise_erosion) - if (is.null(names(noise_erosion))) { - noise_erosion <- noise_erosion[rep(1:length(noise_erosion), length(ROI_noise))[1:length(ROI_noise)]] - names(noise_erosion) <- names(ROI_noise) - } else { - if (all(sort(names(noise_erosion)) != sort(names(ROI_noise)))) { - stop("The names of `noise_erosion` do not match those of `ROI_noise`.") - } - } - noise_erosion <- noise_erosion[names(ROI_noise)] - - } else { - if (!is.null(noise_erosion)) { - warning( - "Erosion requires volumetric data, but the data is not volumetric.\ - No erosion will happen.\n" - ) - } - } - - X_noise <- vector("list", length(ROI_noise)); names(X_noise) <- names(ROI_noise) - for (ii in 1:length(ROI_noise)) { - if (is.null(ROI_noise[[ii]])) { ROI_noise[ii] <- list(NULL); next } - if (X_type == "vector") { - if (is.vector(ROI_noise[[ii]])) { - stopifnot(length(ROI_noise[[ii]]) == V_) - ROI_noise[[ii]] <- as.logical(ROI_noise[[ii]]) - X_noise[[ii]] <- X[,ROI_noise[[ii]]] - } else if (is.matrix(ROI_noise[[ii]])) { - stopifnot(nrow(ROI_noise[[ii]]) == T_) - X_noise[[ii]] <- ROI_noise[[ii]]; ROI_noise[ii] <- list(NULL) - } else { - stop( - "Each entry in `ROI_noise` must be a logical vector, or matrix\ - with the same number of rows as timepoints in `X`." - ) - } - } else if (X_type == "volume") { - if (is.character(ROI_noise[[ii]])) { - if (!file.exists(ROI_noise[[ii]])) { - stop(paste( - "The `ROI_noise` entry", ROI_noise[[ii]], "is not an existing file." - )) - } - ROI_noise[[ii]] <- read_nifti(ROI_noise[[ii]]) - } - if (is.matrix(ROI_noise[[ii]])) { - ROI_noise[[ii]] <- array(ROI_noise[[ii]], dim=c(dim(ROI_noise[[ii]]), 1)) - } - if (is.array(ROI_noise[[ii]])) { - stopifnot(all(dim(ROI_noise[[ii]]) == dim(X)[1:3])) - ROI_noise[[ii]][,,] <- as.logical(ROI_noise[[ii]]) * 1 - ROI_noise[[ii]] <- erode_vol(ROI_noise[[ii]], noise_erosion[[ii]], c(-1, 0, NA, NaN)) - X_noise[[ii]] <- t(matrix(X[ROI_noise[[ii]] > 0], ncol=T_)) - - } else if (is.matrix(ROI_noise[[ii]])) { - stopifnot(nrow(ROI_noise[[ii]]) == T_) - X_noise[[ii]] <- ROI_noise[[ii]]; ROI_noise[ii] <- list(NULL) - } else { - stop( - "Each entry in `ROI_noise` must be a logical array, or matrix\ - with the same number of rows as timepoints in `X`." - ) - } - } else if (X_type == "xifti") { - stopifnot(is.matrix(ROI_noise[[ii]])) - stopifnot(nrow(ROI_noise[[ii]]) == T_) - X_noise[[ii]] <- ROI_noise[[ii]]; ROI_noise[ii] <- list(NULL) - } else { stop("Internal error: unrecognized `X_type`") } - if (ncol(X_noise[[ii]]) == 0) { - warning(paste("The noise ROI", names(ROI_noise)[ii], "is empty.")) - } - } - ROI_noise <- ROI_noise[!vapply(ROI_noise, is.null, FALSE)] - - if (!all(sapply(ROI_noise, is.null))) { - # check that ROI are mutually exclusive - all_ROI_noises <- apply(do.call(rbind, ROI_noise) > 0, 2, sum) - if (!all(all_ROI_noises < 2)) { - stop("The noise ROIs must all be mutually exclusive.") - } - all_ROI_noises <- all_ROI_noises > 0 - if (ROI_data_infer) { - ROI_data[all_ROI_noises] <- FALSE - } else { - if (any(all_ROI_noises & as.vector(ROI_data))) { - warning("The noise ROIs overlapped with the data ROI. Labeling overlapped voxels as noise.") - if (X_type == "volume") { - ROI_data[,,] <- ROI_data & !all_ROI_noises - } else { - ROI_data <- ROI_data & !all_ROI_noises - } - } - } - } - if (length(ROI_noise) == 0) { ROI_noise <- NULL } - - } else { - X_noise <- NULL; noise_nPC <- NULL - } - - # make sure nPCs greater than the rank of the noise ROIs! - # similarly for data - # ... - - if (X_type == "volume") { - X <- t(matrix(X[ROI_data], ncol=dim(X)[4])) - } else { - X <- X[,ROI_data] - } - - list( - X=X, X_noise=X_noise, - ROI_data=ROI_data, ROI_noise=ROI_noise, - noise_nPC=noise_nPC, noise_erosion=noise_erosion - ) -} \ No newline at end of file diff --git a/R/fsl_bptf.R b/R/fsl_bptf.R deleted file mode 100644 index 30fd86f..0000000 --- a/R/fsl_bptf.R +++ /dev/null @@ -1,53 +0,0 @@ -#' \code{bptf} function from FSL -#' -#' Copy of highpass filter as implemented by \code{bptf} in FSL. The results are very -#' similar but not identical. -#' -#' Sources: -#' https://cpb-us-w2.wpmucdn.com/sites.udel.edu/dist/7/4542/files/2016/09/fsl_temporal_filt-15sywxn.m -#' https://github.com/rordenlab/niimath/blob/master/src/core32.c -#' -#' @param orig_data \eqn{T} by \eqn{V} data matrix whose columns will be detrended -#' @param HP_sigma The frequency parameter for the highpass filter -#' -#' @return The data with detrended columns -#' -#' @export -#' -#' @section References: -#' \itemize{ -#' \item{Jenkinson, M., Beckmann, C. F., Behrens, T. E. J., Woolrich, M. W. & Smith, S. M. FSL. NeuroImage 62, 782-790 (2012).} -#' } -#' -fsl_bptf <- function(orig_data, HP_sigma=2000) { - T_ <- nrow(orig_data) - - orig_data <- nuisance_regression(orig_data, cbind(1, seq(T_))) - - HP_filt_size <- ceiling(HP_sigma*3)#round(HP_sigma*8) - HP_lin <- seq(-HP_filt_size/2, HP_filt_size/2, length.out=HP_filt_size) - HP_gfilt <- exp( -(HP_lin^2) / (2*(HP_sigma^2)) ) - HP_gfilt <- HP_gfilt/sum(HP_gfilt) - - filt_data <- matrix(NA, nrow=T_, ncol=ncol(orig_data)) - back <- floor((HP_filt_size-1)/2) - front <- ceiling((HP_filt_size-1)/2) - for (t in seq(T_)) { - if ((t-back < 1) && (t+front > T_)) { - trunc_HP_gfilt <- HP_gfilt[seq(back-t+2, HP_filt_size-(t+front-T_))] - trunc_HP_gfilt <- trunc_HP_gfilt/sum(trunc_HP_gfilt) - filt_data[t,] <- trunc_HP_gfilt %*% orig_data - } else if (t-back < 1) { - trunc_HP_gfilt <- HP_gfilt[seq(back-t+2, HP_filt_size)] - trunc_HP_gfilt <- trunc_HP_gfilt/sum(trunc_HP_gfilt) - filt_data[t,] <- trunc_HP_gfilt %*% orig_data[seq(t+front),] - } else if (t+front > T_) { - trunc_HP_gfilt <- HP_gfilt[seq(HP_filt_size-(t+front-T_))] - trunc_HP_gfilt <- trunc_HP_gfilt/sum(trunc_HP_gfilt) - filt_data[t,] <- trunc_HP_gfilt %*% orig_data[seq(t-back, T_),] - } else { - filt_data[t,] <- HP_gfilt %*% orig_data[seq(t-back, t+front),] - } - } - orig_data - filt_data -} \ No newline at end of file diff --git a/R/leverage.R b/R/leverage.R index 51160ec..cfb30b4 100644 --- a/R/leverage.R +++ b/R/leverage.R @@ -16,6 +16,7 @@ #' \code{!is.null(median_cutoff)}, \code{"cut"} and \code{"flag"} are omitted. #' #' @importFrom stats median +#' @importFrom fMRItools hat_matrix #' #' @export leverage <- function(Comps, are_norm=FALSE, median_cutoff=NULL){ diff --git a/R/nuisance_regression.R b/R/nuisance_regression.R deleted file mode 100644 index d856288..0000000 --- a/R/nuisance_regression.R +++ /dev/null @@ -1,48 +0,0 @@ -#' Hat matrix -#' -#' Get the hat matrix from a design matrix -#' -#' @param design The \eqn{T} by \eqn{Q} design matrix -#' -#' @return The \eqn{T} by \eqn{T} hat matrix -#' -#' @keywords internal -hat_matrix <- function(design){ - design <- as.matrix(design) - # https://stackoverflow.com/questions/19100600/extract-maximal-set-of-independent-columns-from-a-matrix - # https://stackoverflow.com/questions/39167204/in-r-how-does-one-extract-the-hat-projection-influence-matrix-or-values-from-an - qrd <- qr(design) - design <- design[, qrd$pivot[seq_len(qrd$rank)], drop=FALSE] - qrd <- qr(design) - Qd <- qr.Q(qrd) - tcrossprod(Qd) -} - -#' Nuisance regression -#' -#' Performs nuisance regression. The data and design matrix must both be -#' centered, or an intercept must be included in the design matrix! -#' -#' @param Y The \eqn{T} by \eqn{V} or \eqn{V} by \eqn{T} data. -#' @param design The \eqn{T} by \eqn{Q} matrix of nuisance regressors. -#' -#' @return The data after nuisance regression. -#' -#' @export -nuisance_regression <- function(Y, design){ - # Z <- design - # if(nrow(Y) != nrow(Z)) stop('Y and Z must have same number of rows') - # invZtZ <- solve(t(Z) %*% Z) #(Z'Z)^{-1} - # betahat <- invZtZ %*% t(Z) %*% Y #(Z'Z)^{-1} Z'Y - # return(Y - Z %*% betahat) - - Y <- as.matrix(Y); design <- as.matrix(design) - I_m_H <- diag(nrow(design)) - hat_matrix(design) - if (nrow(Y)==nrow(design)) { - return(I_m_H %*% Y) - } else if (ncol(Y)==nrow(design)) { - return(Y %*% I_m_H) - } else { - stop("Y and design are not of compatible dimensions.") - } -} \ No newline at end of file diff --git a/R/pscrub_multi.R b/R/pscrub_multi.R index f42fbc7..ced947e 100644 --- a/R/pscrub_multi.R +++ b/R/pscrub_multi.R @@ -82,7 +82,8 @@ #' @importFrom pesel pesel #' @importFrom robustbase rowMedians #' @importFrom stats mad qnorm var setNames -#' +#' @importFrom fMRItools is_integer is_constant nuisance_regression as.matrix_ifti dct_bases validate_design_mat +#' #' @keywords internal #' pscrub_multi = function( @@ -100,7 +101,7 @@ pscrub_multi = function( # `X` ------------------------------------------------------------------------ if (verbose) { cat("Checking for missing, infinite, and constant data.\n") } - X <- as.matrix2(X, verbose=verbose); class(X) <- "numeric" + X <- as.matrix_ifti(X, verbose=verbose); class(X) <- "numeric" V0_ <- ncol(X) X_NA_mask <- apply(X, 2, function(x){any(x %in% c(NA, NaN, -Inf, Inf))}) if (any(X_NA_mask)) { @@ -111,7 +112,7 @@ pscrub_multi = function( ) X <- X[,!X_NA_mask,drop=FALSE] } - X_const_mask <- apply(X, 2, is_constant) + X_const_mask <- apply(X, 2, fMRItools::is_constant) if (any(X_const_mask)) { if (all(X_const_mask)) { stop("All data columns are constant.\n") } warning( @@ -172,8 +173,8 @@ pscrub_multi = function( } do_nuisance <- !(is.null(nuisance) || isFALSE(nuisance) || identical(nuisance, 0)) if (do_nuisance) { - nuisance <- check_design_matrix(nuisance, T_) - design_const_mask <- apply(nuisance, 2, is_constant) + nuisance <- validate_design_mat(nuisance, T_) + design_const_mask <- apply(nuisance, 2, fMRItools::is_constant) if (!any(design_const_mask)) { if (!any(abs(apply(X, 2, mean)) < 1e-8)) { warning("No intercept detected in `design`, yet the data are not centered.") @@ -190,7 +191,7 @@ pscrub_multi = function( comps_mean_dt <- 0 } else { comps_mean_dt <- as.numeric(comps_mean_dt) - stopifnot(is_integer(comps_mean_dt, nneg=TRUE)) + stopifnot(fMRItools::is_integer(comps_mean_dt, nneg=TRUE)) } if (isTRUE(comps_var_dt)) { comps_var_dt <- 4 @@ -198,7 +199,7 @@ pscrub_multi = function( comps_var_dt <- 0 } else { comps_var_dt <- as.numeric(comps_var_dt) - stopifnot(is_integer(comps_var_dt, nneg=TRUE)) + stopifnot(fMRItools::is_integer(comps_var_dt, nneg=TRUE)) } comps_dt <- (comps_mean_dt > 0) || (comps_var_dt > 0) kurt_quantile <- as.numeric(kurt_quantile) @@ -214,7 +215,7 @@ pscrub_multi = function( get_dirs <- as.logical(get_dirs); stopifnot(isTRUE(get_dirs) || isFALSE(get_dirs)) full_PCA <- as.logical(full_PCA); stopifnot(isTRUE(full_PCA) || isFALSE(full_PCA)) get_outliers <- as.logical(get_outliers); stopifnot(isTRUE(get_outliers) || isFALSE(get_outliers)) - cutoff <- as.numeric(cutoff); stopifnot(is_integer(cutoff, nneg=TRUE)) + cutoff <- as.numeric(cutoff); stopifnot(fMRItools::is_integer(cutoff, nneg=TRUE)) verbose <- as.logical(verbose); stopifnot(isTRUE(verbose) || isFALSE(verbose)) # ---------------------------------------------------------------------------- @@ -232,7 +233,7 @@ pscrub_multi = function( cat(action, "the data columns.\n") } - # Center & scale here, rather than calling `scale_med`, to save memory. + # Center & scale here, rather than calling `fMRItools::scale_med`, to save memory. X <- t(X) if (center) { X <- X - c(rowMedians(X)) } if (scale) { X <- X / (1.4826 * rowMedians(abs(X))) } @@ -350,7 +351,7 @@ pscrub_multi = function( out$fusedPCA$PC_exec_times <- NULL; out$fusedPCA$nItes <- NULL # V matrix from PCA no longer needed. if(!get_dirs){ out$PCA$V <- NULL } - + tf_const_mask <- apply(out$fusedPCA$u, 2, is_constant) if(any(tf_const_mask)){ warning( diff --git a/R/rob_stabilize.R b/R/rob_stabilize.R index d6c3121..4592597 100644 --- a/R/rob_stabilize.R +++ b/R/rob_stabilize.R @@ -11,13 +11,14 @@ #' #' @return The output of \code{robustbase::lmrob} #' +#' @importFrom fMRItools is_integer dct_bases #' @keywords internal rob_trend <- function(x, nDCT=4, lmrob_method="MM", seed=0) { x <- as.vector(x) T_ <- length(x) nDCT <- as.numeric(nDCT) - stopifnot(is_integer(nDCT)); stopifnot(nDCT >= 0) + stopifnot(fMRItools::is_integer(nDCT, nneg=TRUE)) if (nDCT == 0) { mat <- data.frame(rep(1, T_)) colnames(mat) <- "x_int" @@ -53,6 +54,7 @@ rob_trend <- function(x, nDCT=4, lmrob_method="MM", seed=0) { #' #' @return the timeseries with its center and scale stabilized #' +#' @importFrom fMRItools is_integer #' @export #' rob_stabilize <- function(x, center=TRUE, scale=TRUE, lmrob_method="MM", rescale=TRUE) { @@ -69,7 +71,7 @@ rob_stabilize <- function(x, center=TRUE, scale=TRUE, lmrob_method="MM", rescale center <- 0 } else { center <- as.numeric(center) - stopifnot(is_integer(center, nneg=TRUE)) + stopifnot(fMRItools::is_integer(center, nneg=TRUE)) } if (isTRUE(scale)) { @@ -78,7 +80,7 @@ rob_stabilize <- function(x, center=TRUE, scale=TRUE, lmrob_method="MM", rescale scale <- 0 } else { scale <- as.numeric(scale) - stopifnot(is_integer(scale, nneg=TRUE)) + stopifnot(fMRItools::is_integer(scale, nneg=TRUE)) } if (center > 0) { diff --git a/R/rox_args_docs.R b/R/rox_args_docs.R index 4244b4d..c2815f4 100644 --- a/R/rox_args_docs.R +++ b/R/rox_args_docs.R @@ -21,7 +21,7 @@ #' non-time-series data because the observations are not temporally related. #' #' Additional nuisance regressors can be specified like so: -#' \code{cbind(1, dct_bases(nrow(x), 4), more_nuisance)}. +#' \code{cbind(1, fMRItools::dct_bases(nrow(x), 4), more_nuisance)}. #' @param 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: \code{TRUE}. #' Centering is necessary for computing the projections, so if \code{center} is diff --git a/R/utils.R b/R/utils.R deleted file mode 100644 index 191fd88..0000000 --- a/R/utils.R +++ /dev/null @@ -1,210 +0,0 @@ -#' Is this an integer? -#' -#' @param x The putative integer -#' @param nneg Require \code{x>=0} (non-negative) too? -#' @return Logical indicating whether \code{x} is an integer -#' -#' @keywords internal -is_integer <- function(x, nneg=FALSE){ - y <- FALSE - if (length(x)==1 && is.numeric(x)) { - if (x%%1==0) { - if (x>=0 || !nneg) { y <- TRUE } - } - } - y -} - -#' Is this numeric vector constant? -#' -#' @param x The numeric vector -#' @param TOL minimum range of \code{x} to be considered non-constant. -#' Default: \code{1e-8} -#' -#' @return Is \code{x} constant? -#' -#' @keywords internal -is_constant <- function(x, TOL=1e-8) { - abs(max(x) - min(x)) < TOL -} - -#' Check design matrix -#' -#' @param design The design matrix -#' -#' @return The (modified) design matrix -#' -#' @keywords internal -check_design_matrix <- function(design, T_=nrow(design)) { - class(design) <- "numeric" - if (identical(design, 1)) { design <- matrix(1, nrow=T_) } - design <- as.matrix(design) - stopifnot(nrow(design) == T_) - # Set constant columns (intercept regressor) to 1, and scale the other columns. - design_const_mask <- apply(design, 2, is_constant) - design[,design_const_mask] <- 1 - design[,!design_const_mask] <- scale(design[,!design_const_mask]) - design -} - -#' Robust scaling -#' -#' Centers and scales the columns of a matrix robustly -#' -#' Centers each column on its median, and scales each column by its median -#' absolute deviation (MAD). Constant-valued columns are set to \code{NA} -#' (or removed if \code{drop_const}) and a warning is raised. If all -#' MADs are zero, an error is raised. -#' -#' @param mat A numerical matrix. -#' @param TOL minimum MAD to consider a column non-constant. -#' Default: \code{1e-8} -#' @param drop_const Drop -#' -#' @return The input matrix with its columns centered and scaled. -#' -#' @importFrom robustbase rowMedians -scale_med <- function(mat, TOL=1e-8, drop_const=TRUE){ - # Transpose. - mat <- t(mat) - - # Center. - mat <- mat - c(rowMedians(mat, na.rm=TRUE)) - - # Scale. - mad <- 1.4826 * rowMedians(abs(mat), na.rm=TRUE) - mad <- as.numeric(mad) - const_mask <- mad < TOL - if (any(const_mask)) { - if (all(const_mask)) { - stop("All columns are zero-variance.\n") - } else { - warning(paste0( - "Warning: ", sum(const_mask), - " constant columns (out of ", length(const_mask), - " ). These will be removed.\n" - )) - } - } - mad <- mad[!const_mask] - mat[const_mask,] <- NA - mat[!const_mask,] <- mat[!const_mask,] / mad - - # Revert transpose. - mat <- t(mat) - - if (drop_const) { mat <- mat[!const_mask,] } - - mat -} - -#' Infer fMRI data format -#' -#' @param BOLD The fMRI data -#' @param verbose Print the format? Default: \code{FALSE}. -#' @return The format: \code{"CIFTI"} file path, \code{"xifti"} object, -#' \code{"NIFTI"} file path, \code{"nifti"} object, or \code{"data"}. -#' @keywords internal -infer_BOLD_format <- function(BOLD, verbose=FALSE){ - - # Character vector: CIFTI or NIFTI - if (is.character(BOLD)) { - stopifnot(length(BOLD)==1) - if (endsWith(BOLD, ".dtseries.nii") | endsWith(BOLD, ".dscalar.nii")) { - format <- "CIFTI" - } else if (endsWith(BOLD, "gii")) { - format <- "GIFTI" - } else { - format <- "NIFTI" - } - - } else if (inherits(BOLD, "xifti")) { - format <- "xifti" - } else if (inherits(BOLD, "niftiExtension") | inherits(BOLD, "niftiImage") | inherits(BOLD, "nifti")) { - format <- "nifti" - } else if (inherits(BOLD, "gifti")) { - format <- "gifti" - - } else { - if (is.list(BOLD)) { stop("Unknown `BOLD` format.") } - if (is.matrix(BOLD)) { - format <- "data" - } else if (length(dim(BOLD))==4) { - format <- "nifti" - } else { - stop("Unknown `BOLD` format.") - } - } - if (verbose) { cat("Inferred input format:", format, "\n") } - format -} - -#' Wrapper to common functions for reading NIFTIs -#' -#' Tries \code{RNifti::readNifti}, then \code{oro.nifti::readNIfTI}. If -#' neither package is available an error is raised. -#' -#' @param nifti_fname The file name of the NIFTI. -#' @return The NIFTI -#' @keywords internal -read_nifti <- function(nifti_fname){ - if (requireNamespace("RNifti", quietly = TRUE)) { - return(RNifti::readNifti(nifti_fname)) - } else if (requireNamespace("oro.nifti", quietly = TRUE)) { - return(oro.nifti::readNIfTI(nifti_fname, reorient=FALSE)) - } else { - stop( - "Package \"RNifti\" or \"oro.nifti\" needed to read", nifti_fname, - ". Please install at least one", call. = FALSE - ) - } -} - -#' Convert CIFTI, NIFTI, or GIFTI file input to \eqn{T} by \eqn{V} matrix -#' -#' Convert input data to \eqn{T} by \eqn{V} matrix. Reads in a CIFTI, NIFTI, -#' or GIFTI file. Transposes CIFTI and xifti object input. Masks NIFTI and -#' nifti object input. -#' -#' @param x The object to coerce to a matrix -#' @param sortSub Sort subcortex by labels? Default: \code{FALSE} -#' @param verbose Print updates? Default: \code{FALSE} -#' @return x as a matrix. -#' @keywords internal -as.matrix2 <- function(x, sortSub=FALSE, verbose=FALSE) { - - format <- infer_BOLD_format(x) - if (format %in% c("CIFTI", "xifti")) { - if (format == "CIFTI") { - if (!requireNamespace("ciftiTools", quietly = TRUE)) { - stop("Package \"ciftiTools\" needed to read input data. Please install it", call. = FALSE) - } - x <- ciftiTools::read_cifti(x, brainstructures=ciftiTools::info_cifti(x)$cifti$brainstructures) - } - if (sortSub && !is.null(x$data$subcort)) { - x$data$subcort <- x$data$subcort[order(x$meta$subcort$labels),] - } - x <- t(as.matrix(x)) - } else if (format %in% c("GIFTI", "gifti")) { - if (format == "GIFTI") { - if (!requireNamespace("gifti", quietly = TRUE)) { - stop("Package \"gifti\" needed to read `X`. Please install it", call. = FALSE) - } - x <- gifti::read_gifti(x) - } - x <- t(do.call(cbind, x$data)) - if (verbose) {cat("GIFTI dimensions:\n"); print(dim(x))} - } else if (format %in% c("NIFTI", "nifti")) { - if (format == "NIFTI") { - x <- read_nifti(x) - } - if (verbose) {cat("Masking NIFTI by removing locations with constant zero, NA, or NaN.\n")} - z <- array(x %in% c(0, NA, NaN), dim=dim(x)) - mask <- !apply(z, seq(3), all) - x <- matrix(x[rep(mask, dim(x)[4])], ncol=dim(x)[4]) - x <- t(x) - if (verbose) {cat("NIFTI dimensions:\n"); print(dim(x))} - } - - x -} \ No newline at end of file diff --git a/README.Rmd b/README.Rmd index aff8dd6..d6c3238 100644 --- a/README.Rmd +++ b/README.Rmd @@ -12,7 +12,7 @@ output: github_document [![Codecov test coverage](https://codecov.io/gh/mandymejia/fMRIscrub/branch/master/graph/badge.svg)](https://app.codecov.io/gh/mandymejia/fMRIscrub?branch=master) -`fMRIscrub` is a collection of routines for data-driven scrubbing (projection scrubbing and DVARS), motion scrubbing, and other fMRI denoising strategies such as anatomical CompCor, detrending, and nuisance regression. The data-driven scrubbing methods are also applicable to other outlier detection tasks involving high-dimensional data. +`fMRIscrub` is a collection of routines for data-driven scrubbing (projection scrubbing and DVARS), motion scrubbing, and other fMRI denoising strategies such as anatomical CompCor, detrending, and nuisance regression. Projection scrubbing is also applicable to other outlier detection tasks involving high-dimensional data. ## Installation diff --git a/README.md b/README.md index 6771893..694b0bd 100644 --- a/README.md +++ b/README.md @@ -15,9 +15,8 @@ coverage](https://codecov.io/gh/mandymejia/fMRIscrub/branch/master/graph/badge.s `fMRIscrub` is a collection of routines for data-driven scrubbing (projection scrubbing and DVARS), motion scrubbing, and other fMRI denoising strategies such as anatomical CompCor, detrending, and -nuisance regression. The data-driven scrubbing methods are also -applicable to other outlier detection tasks involving high-dimensional -data. +nuisance regression. Projection scrubbing is also applicable to other +outlier detection tasks involving high-dimensional data. ## Installation diff --git a/man/CompCor.Rd b/man/CompCor.Rd deleted file mode 100644 index 970cc2a..0000000 --- a/man/CompCor.Rd +++ /dev/null @@ -1,148 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/CompCor.R -\name{CompCor} -\alias{CompCor} -\title{Anatomical CompCor} -\usage{ -CompCor( - X, - ROI_data = "infer", - ROI_noise = NULL, - noise_nPC = 5, - noise_erosion = NULL, - center = TRUE, - scale = TRUE, - nuisance = NULL -) -} -\arguments{ -\item{X}{Wide numeric data matrix (\eqn{T observations} by \eqn{V variables}, \eqn{T << V}). -For example, if \code{X} represents an fMRI run, \eqn{T} should be the number -of timepoints and \eqn{V} should be the number of brainordinate vertices/voxels. - -Or, a 4D array or NIFTI or file path to a NIFTI (\eqn{I} by \eqn{J} by \eqn{K} by \eqn{T} -observations), in which case \code{ROI_data} must be provided. -(The vectorized data will be \eqn{T timepoints} by \eqn{V_{in-mask} voxels}) - -Or, a \code{ciftiTools} \code{"xifti"} object or a file path to a CIFTI -(The vectorized data will be \eqn{T timepoints} by \eqn{V_{left+right+sub} greyordinates}).} - -\item{ROI_data}{Indicates the data ROI. Allowed arguments depend on \code{X}: - -If \code{X} is a matrix, this must be a length \eqn{V} logical vector, where -the data ROI is indicated by \code{TRUE} values. If \code{"infer"} (default), all -columns of \code{X} will be included in the data ROI (\code{rep(TRUE, V)}). - -If \code{X} is an array or NIFTI, this must be either a vector of values -to expect for out-of-mask voxels in \code{X}, or a (file path to a) 3D NIFTI. -In the latter case, each of the volume dimensions should match the first -three dimensions of \code{X}. Voxels in the data ROI should be indicated by -\code{TRUE} and all other voxels by \code{FALSE}. If \code{"infer"} (default), -will be set to \code{c(0, NA, NaN)} (include all voxels which are not constant -\code{0}, \code{NA}, or \code{NaN}). - -If \code{X} is a \code{"xifti"} this must be the \code{brainstructures} -argument to \code{\link[ciftiTools]{read_cifti}}. If \code{"infer"} (default), -\code{brainstructures} will be set to \code{"all"} (use both left and right -cortex vertices, and subcortical voxels). - -If \code{NULL}, the data ROI will be empty. This is useful for obtaining just -the noise ROI, if the data and noise are located in separate files.} - -\item{ROI_noise}{Indicates the noise ROIs for aCompCor. Should be a list where -each entry corresponds to a distinct noise ROI. The names of the list should -be the ROI names, e.g. \code{"white_matter"} and \code{"csf"}. The expected -formats of the list entries depends on \code{X}: - -For all types of \code{X}, \code{ROI_noise} entries can be a matrix of noise -ROI data. The matrix should have \eqn{T} rows, with each column being a -data location's timeseries. - -If \code{X} is a matrix, entries can also indicate a noise ROI within \code{X}. -These entries must be a length \eqn{V} logical vector with \code{TRUE} values -indicating locations in \code{X} within that noise ROI. Since the ROIs must -not overlap, the masks must be mutually exclusive with each other, and with -\code{ROI_data}. - -If \code{X} is an array or NIFTI, entries can also indicate a noise ROI within \code{X}. -These entries must be a logical array or (file path to) a 3D NIFTI with the -same spatial dimensions as \code{X}, and with \code{TRUE} values indicating -voxels inside the noise ROI. Since the ROIs must not overlap, the masks must -be mutually exclusive with each other, and with \code{ROI_data}. - -(If \code{X} is a \code{"xifti"}, entries must be data matrices, since no -greyordinate locations in \code{X} are appropriate noise ROIs).} - -\item{noise_nPC}{The number of principal components to compute for each noise -ROI. Alternatively, values between 0 and 1, in which case they will -represent the minimum proportion of variance explained by the PCs used for -each noise ROI. The smallest number of PCs will be used to achieve this -proportion of variance explained. - -Should be a list or numeric vector with the same length as \code{ROI_noise}. -It will be matched to each ROI based on the name of each entry, or if the -names are missing, the order of entries. If it is an unnamed vector, its -elements will be recycled. Default: \code{5} (compute the top 5 PCs for -each noise ROI).} - -\item{noise_erosion}{The number of voxel layers to erode the noise ROIs by. -Should be a list or numeric vector with the same length as \code{ROI_noise}. -It will be matched to each ROI based on the name of each entry, or if the -names are missing, the order of entries. If it is an unnamed vector, its -elements will be recycled. Default: \code{NULL}, which will use a value of -0 (do not erode the noise ROIs). Note that noise erosion can only be -performed if the noise ROIs are volumetric.} - -\item{center, scale}{Center the columns of the noise ROI data by their medians, -and scale by their MADs? Default: \code{TRUE} for both. Note that this argument -affects the noise ROI data and not the data that is being cleaned with aCompCor. -Centering and scaling of the data being cleaned can be done after this function call.} - -\item{nuisance}{Nuisance signals to regress from each data column in addition -to the noise ROI PCs. Should be a \eqn{T} by \eqn{N} numeric matrix where -\eqn{N} represents the number of nuisance signals. To not perform any nuisance -regression set this argument to \code{NULL}, \code{0}, or \code{FALSE}. -Default: \code{NULL}.} -} -\value{ -A list with entries \code{"data"}, \code{"noise"}, and potentially -\code{"ROI_data"}. - -The entry \code{"data"} will be a \code{V x T} matrix where each row corresponds to a -data location (if it was originally an array, the locations will be voxels in spatial -order). Each row will be a time series with each noise PC regressed from it. This entry -will be \code{NULL} if there was no data. - -The entry \code{"noise"} is a list of noise PC scores, their corresponding variance, -and their ROI mask, for each noise ROI. - -If the data ROI is not all \code{TRUE}, the entry \code{"ROI_data"} will have -the ROI mask for the data. -} -\description{ -The aCompCor algorithm for denoising fMRI data using noise ROIs data -} -\details{ -First, the principal components (PCs) of each noise region of interest (ROI) -are calculated. For each ROI, voxels are centered and scaled -(can be disabled with the arguments \code{center} and \code{scale}), -and then the PCs are calculated via the singular value decomposition. - -Next, aCompCor is performed to remove the shared variation between the noise ROI -PCs and each location in the data. This is accomplished by a nuisance regression -using a design matrix with the noise ROI PCs, any additional regressors specified -by \code{nuisance}, and an intercept term. (To detrend the data and perform aCompCor -in the same regression, \code{nuisance} can be set to DCT bases obtained with -the function \code{\link{dct_bases}}.) -} -\section{References}{ - -\itemize{ -\item{Behzadi, Y., Restom, K., Liau, J. & Liu, T. T. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage 37, 90-101 (2007).} -\item{Muschelli, J. et al. Reduction of motion-related artifacts in resting state fMRI using aCompCor. NeuroImage 96, 22-35 (2014).} -} -} - -\seealso{ -CompCor_HCP -} diff --git a/man/CompCor.noise_comps.Rd b/man/CompCor.noise_comps.Rd deleted file mode 100644 index b6018da..0000000 --- a/man/CompCor.noise_comps.Rd +++ /dev/null @@ -1,23 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/CompCor.R -\name{CompCor.noise_comps} -\alias{CompCor.noise_comps} -\title{CompCor: get noise components} -\usage{ -CompCor.noise_comps(X_noise, center, scale, noise_nPC) -} -\arguments{ -\item{X_noise}{The noise ROIs data} - -\item{center, scale}{Center & scale robustly} - -\item{noise_nPC}{Number of PCs to obtain for each noise ROI} -} -\value{ -A list with components X, X_noise, ROI_data, ROI_noise, noise_nPC, -noise_erosion, noise_comps, and noise_var. -} -\description{ -Get noise components for aCompCor. -} -\keyword{internal} diff --git a/man/CompCor_HCP.Rd b/man/CompCor_HCP.Rd deleted file mode 100644 index c4638dd..0000000 --- a/man/CompCor_HCP.Rd +++ /dev/null @@ -1,123 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/CompCor_HCP.R -\name{CompCor_HCP} -\alias{CompCor_HCP} -\title{Anatomical CompCor for HCP NIFTI and CIFTI data} -\usage{ -CompCor_HCP( - nii, - nii_labels, - ROI_noise = c("wm_cort", "csf"), - noise_nPC = 5, - noise_erosion = NULL, - idx = NULL, - cii = NULL, - brainstructures = c("left", "right"), - center = TRUE, - scale = TRUE, - DCT = 0, - nuisance_too = NULL, - verbose = FALSE -) -} -\arguments{ -\item{nii}{\eqn{I} by \eqn{J} by \eqn{K} by \eqn{T} -NIFTI object or array (or file path to the NIFTI) which contains -whole-brain data, including the noise ROIs. In the HCP, the corresponding -file is e.g. "../Results/rfMRI_REST1_LR/rfMRI_REST1_LR.nii.gz"} - -\item{nii_labels}{\eqn{I} by \eqn{J} by \eqn{K} -NIFTI object or array (or file path to the NIFTI) which -contains the corresponding labels to each voxel in \code{nii}. Values should -be according to this table: -https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/AnatomicalROI/FreeSurferColorLUT . -In the HCP, the corresponding file is "ROIs/Atlas_wmparc.2.nii.gz".} - -\item{ROI_noise}{A list of numeric vectors. Each entry should represent labels -in \code{nii_labels} belonging to a single noise ROI, named by that entry's -name. Or, this can be a character vector of at least one of the following: -\code{"wm_cort"} (cortical white matter), \code{"wm_cblm"} (cerebellar white -matter), \code{"csf"} (cerebrospinal fluid). In the latter case, these labels -will be used: - -\describe{ -\item{\code{"wm_cort"}}{\code{c(3000:4035, 5001, 5002)}} -\item{\code{"wm_cblm"}}{\code{c(7, 46)}} -\item{\code{"csf"}}{\code{c(4, 5, 14, 15, 24, 31, 43, 44, 63, 250, 251, 252, 253, 254, 255))}} -} - -These default ROIs are based on this forum post: -https://www.mail-archive.com/hcp-users@humanconnectome.org/msg00931.html - -Default: \code{c("wm_cort", "csf")}} - -\item{noise_nPC}{The number of principal components to compute for each noise -ROI. Alternatively, values between 0 and 1, in which case they will -represent the minimum proportion of variance explained by the PCs used for -each noise ROI. The smallest number of PCs will be used to achieve this -proportion of variance explained. - -Should be a list or numeric vector with the same length as \code{ROI_noise}. -It will be matched to each ROI based on the name of each entry, or if the -names are missing, the order of entries. If it is an unnamed vector, its -elements will be recycled. Default: \code{5} (compute the top 5 PCs for -each noise ROI).} - -\item{noise_erosion}{The number of voxel layers to erode the noise ROIs by. -Should be a list or numeric vector with the same length as \code{ROI_noise}. -It will be matched to each ROI based on the name of each entry, or if the -names are missing, the order of entries. If it is an unnamed vector, its -elements will be recycled. Default: \code{NULL}, which will use a value of -0 (do not erode the noise ROIs).} - -\item{idx}{A numeric vector indicating the timepoints to use, or -\code{NULL} (default) to use all idx. (Indexing begins with 1, so the -first timepoint has index 1 and the last has the same index as the length of -the scan.)} - -\item{cii}{\code{"xifti"} (or file path to the CIFTI) from which the noise -ROI components will be regressed. In the HCP, the corresponding file is e.g. -"../Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_MSMAll.dtseries.nii". If not -provided, only the noise components will be returned (no data will be cleaned).} - -\item{brainstructures}{Choose among "left", "right", and "subcortical". -Default: \code{c("left", "right")} (cortical data only)} - -\item{center, scale}{Center the columns of the data by median, and scale the -columns of the data by MAD? Default: \code{TRUE} for both. Affects both -\code{X} and the noise data. \code{center} also applies to \code{nuisance_too} -so if it is \code{FALSE}, \code{nuisance_too} must already be centered.} - -\item{DCT}{Add DCT bases to the nuisance regression? Use an integer to -indicate the number of cosine bases. Use \code{0} (default) to forgo detrending. - -The data must be centered, either before input or with \code{center}.} - -\item{nuisance_too}{A matrix of nuisance signals to add to the nuisance -regression. Should have \eqn{T} rows. \code{NULL} to not add additional -nuisance regressors (default).} - -\item{verbose}{Should occasional updates be printed? Default: \code{FALSE}.} -} -\value{ -The noise components, and if \code{cii} is provided, the cleaned -surface-based data as a \code{"xifti"} object. -} -\description{ -Wrapper to \code{\link{CompCor}} for HCP-format data. Can be used to clean -the surface-based CIFTI data with aCompCor using the noise PCs and ROIs -calculated from the NIFTI fMRI data and NIFTI mask. Can also be used to just -obtain the noise PCs and ROIs without performing aCompCor, if the CIFTI -data is not provided. -} -\section{References}{ - -\itemize{ -\item{Behzadi, Y., Restom, K., Liau, J. & Liu, T. T. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage 37, 90-101 (2007).} -\item{Muschelli, J. et al. Reduction of motion-related artifacts in resting state fMRI using aCompCor. NeuroImage 96, 22-35 (2014).} -} -} - -\seealso{ -CompCor -} diff --git a/man/DVARS.Rd b/man/DVARS.Rd index 718c392..4a3dd7e 100644 --- a/man/DVARS.Rd +++ b/man/DVARS.Rd @@ -8,7 +8,7 @@ DVARS( X, normalize = TRUE, cutoff_DPD = 5, - cutoff_ZD = qnorm(1 - 0.05/nrow(as.matrix2(X))), + cutoff_ZD = qnorm(1 - 0.05/nrow(as.matrix_ifti(X))), verbose = FALSE ) } diff --git a/man/Mode.Rd b/man/Mode.Rd deleted file mode 100644 index 86f55c7..0000000 --- a/man/Mode.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DVARS.R -\name{Mode} -\alias{Mode} -\title{Mode of data vector} -\usage{ -Mode(x) -} -\arguments{ -\item{x}{The data vector} -} -\value{ -The mode -} -\description{ -Get mode of a data vector. But use the median instead of the mode if all -data values are unique. -} -\keyword{internal} diff --git a/man/as.matrix2.Rd b/man/as.matrix2.Rd deleted file mode 100644 index 259c2aa..0000000 --- a/man/as.matrix2.Rd +++ /dev/null @@ -1,24 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{as.matrix2} -\alias{as.matrix2} -\title{Convert CIFTI, NIFTI, or GIFTI file input to \eqn{T} by \eqn{V} matrix} -\usage{ -as.matrix2(x, sortSub = FALSE, verbose = FALSE) -} -\arguments{ -\item{x}{The object to coerce to a matrix} - -\item{sortSub}{Sort subcortex by labels? Default: \code{FALSE}} - -\item{verbose}{Print updates? Default: \code{FALSE}} -} -\value{ -x as a matrix. -} -\description{ -Convert input data to \eqn{T} by \eqn{V} matrix. Reads in a CIFTI, NIFTI, -or GIFTI file. Transposes CIFTI and xifti object input. Masks NIFTI and -nifti object input. -} -\keyword{internal} diff --git a/man/carpetplot.Rd b/man/carpetplot.Rd deleted file mode 100644 index ee0ad24..0000000 --- a/man/carpetplot.Rd +++ /dev/null @@ -1,66 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/carpetplot.R -\name{carpetplot} -\alias{carpetplot} -\title{Carpetplot} -\usage{ -carpetplot( - x, - qcut = 0.1, - fname = NULL, - center = TRUE, - scale = FALSE, - colors = "gray255", - sortSub = TRUE, - ... -) -} -\arguments{ -\item{x}{The \eqn{T x V} numeric data matrix, or a \code{"xifti"} object. -In the plot, the \eqn{T} index will increase from left to right, and the -\eqn{V} will increase from top to bottom.} - -\item{qcut}{Sets blackpoint at the \code{qcut} quantile, and the -whitepoint at the \code{1-qcut} quantile. Default: \code{.1}. This is -equivalent to setting the color range between the 10\% and 90\% quantiles. -The quantiles are computed across the entire data matrix after any -centering or scaling. - -Must be between 0 and .49. If \code{0} or \code{NULL} (default), do not -clamp the data values.} - -\item{fname}{A \code{.pdf} (highly recommended) or \code{.png} file path -to write the carpetplot to. If \code{NULL} (default), return the plot directly -instead of writing a file.} - -\item{center, scale}{Center and scale the data? If \code{x} is fMRI data -which has not otherwise been centered or scaled, it is recommended to center -but not scale it (default).} - -\item{colors}{\code{"gray255"} (default) will use a grayscale color ramp -from black to white. Otherwise, this should be a character vector of -color names to use. - -Colors will be assigned from the lowest to the highest data value, after -any clamping of the data values by \code{qcut}.} - -\item{sortSub}{If \code{x} is a \code{"xifti"} object with subcortical data, -should the voxels be sorted by structure alphabetically? Default: \code{TRUE}.} - -\item{...}{Additional arguments to \code{pdf} or \code{png}, such as width -and height.} -} -\value{ -The image or \code{NULL, invisibly} if a file was written. -} -\description{ -Plot a matrix with \code{graphics::image}. For fMRI data, this is the "carpetplot" -or grayplot coined by (Power, 2017). The \code{graphics} package is required. -} -\section{References}{ - -\itemize{ -\item{Power, J. D. A simple but useful way to assess fMRI scan qualities. NeuroImage 154, 150-158 (2017).} -} -} - diff --git a/man/carpetplot_stack.Rd b/man/carpetplot_stack.Rd deleted file mode 100644 index 3153144..0000000 --- a/man/carpetplot_stack.Rd +++ /dev/null @@ -1,45 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/carpetplot.R -\name{carpetplot_stack} -\alias{carpetplot_stack} -\title{Stacked carpetplot} -\usage{ -carpetplot_stack( - x_list, - center = TRUE, - scale = FALSE, - qcut = 0.1, - match_scale = TRUE, - nsep = 0, - ... -) -} -\arguments{ -\item{x_list}{List of data matrices} - -\item{center, scale}{Center and scale the data? If \code{x} is fMRI data -which has not otherwise been centered or scaled, it is recommended to center -but not scale it (default).} - -\item{qcut}{Sets blackpoint at the \code{qcut} quantile, and the -whitepoint at the \code{1-qcut} quantile. Default: \code{.1}. This is -equivalent to setting the color range between the 10\% and 90\% quantiles. -The quantiles are computed across the entire data matrix after any -centering or scaling. - -Must be between 0 and .49. If \code{0} or \code{NULL} (default), do not -clamp the data values.} - -\item{match_scale}{Match the scales of the carpetplots? Default: \code{TRUE}.} - -\item{nsep}{Equivalent number of data locations for size of gap between -carpetplots. Default: zero (no gap).} - -\item{...}{Additional arguments to \code{carpetplot}} -} -\value{ -\code{NULL}, invisibly -} -\description{ -Stacks carpetplots on top of one anotherr by rbinding the matrices. -} diff --git a/man/check_design_matrix.Rd b/man/check_design_matrix.Rd deleted file mode 100644 index 370319b..0000000 --- a/man/check_design_matrix.Rd +++ /dev/null @@ -1,18 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{check_design_matrix} -\alias{check_design_matrix} -\title{Check design matrix} -\usage{ -check_design_matrix(design, T_ = nrow(design)) -} -\arguments{ -\item{design}{The design matrix} -} -\value{ -The (modified) design matrix -} -\description{ -Check design matrix -} -\keyword{internal} diff --git a/man/dct_bases.Rd b/man/dct_bases.Rd deleted file mode 100644 index 4f21b00..0000000 --- a/man/dct_bases.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DCT_FFT_detrending.R -\name{dct_bases} -\alias{dct_bases} -\title{Generate cosine bases for the DCT} -\usage{ -dct_bases(T_, n) -} -\arguments{ -\item{T_}{Length of timeseries} - -\item{n}{Number of cosine bases} -} -\value{ -Matrix with cosine bases along columns -} -\description{ -Generate cosine bases for the DCT -} diff --git a/man/dct_convert.Rd b/man/dct_convert.Rd deleted file mode 100644 index 400312b..0000000 --- a/man/dct_convert.Rd +++ /dev/null @@ -1,39 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DCT_FFT_detrending.R -\name{dct_convert} -\alias{dct_convert} -\alias{dct2Hz} -\alias{Hz2dct} -\title{DCT and frequency conversion} -\usage{ -dct_convert(T_, TR, n = NULL, f = NULL) - -dct2Hz(T_, TR, n) - -Hz2dct(T_, TR, f) -} -\arguments{ -\item{T_}{Length of timeseries (number of timepoints)} - -\item{TR}{TR of the fMRI scan, in seconds (the time between timepoints)} - -\item{n}{Number of cosine bases} - -\item{f}{Hz of highpass filter} -} -\value{ -If \code{n} was provided, the highpass filter cutoff (Hz) is returned. -Otherwise, if \code{f} was provided, the number of cosine bases is returned. -The result should be rounded before passing to \code{\link{dct_bases}} -} -\description{ -Convert between number of DCT bases and Hz of highpass filter -} -\details{ -Provide either \code{n} or \code{f} to calculate the other. - -If only the total length of the scan is known, you can set that to \code{TR} -and use \code{T_=1}. - -\eqn{f = n / (2 * T_ * TR)} -} diff --git a/man/detrend.Rd b/man/detrend.Rd deleted file mode 100644 index 5adf1dc..0000000 --- a/man/detrend.Rd +++ /dev/null @@ -1,24 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DCT_FFT_detrending.R -\name{detrend} -\alias{detrend} -\title{Detrending with DCT or FFT} -\usage{ -detrend(X, TR, f = 0.008, method = c("DCT", "FFT")) -} -\arguments{ -\item{X}{\eqn{T \times V} numeric matrix. Each column is a voxel or vertex -time series.} - -\item{TR}{TR of the fMRI scan, in seconds (the time between timepoints)} - -\item{f}{Hz of highpass filter. Default: \code{.008}} - -\item{method}{\code{"DCT"} (default) or \code{"FFT"}.} -} -\value{ -Detrended \code{X} -} -\description{ -Detrending with DCT or FFT -} diff --git a/man/erode_vol.Rd b/man/erode_vol.Rd deleted file mode 100644 index 786d2b7..0000000 --- a/man/erode_vol.Rd +++ /dev/null @@ -1,32 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/format_data.R -\name{erode_vol} -\alias{erode_vol} -\title{Erode volumetric mask} -\usage{ -erode_vol(vol, n_erosion = 1, out_of_mask_val = NA) -} -\arguments{ -\item{vol}{The volume to erode. Out-of-mask voxels should be indicated by a -value in \code{out_of_mask_val}.} - -\item{n_erosion}{The number of layers to erode the mask by.} - -\item{out_of_mask_val}{A voxel is not included in the mask if and only if its -value is in this vector. The first value in this vector will be used to replace -the eroded voxels. Default: \code{NA}.} -} -\value{ -The eroded \code{vol}. It is the same as \code{vol} but with eroded -voxels replaced with the value \code{out_of_mask_val[1]}. -} -\description{ -Erode a volumetric mask by a certain number of voxel layers. For each layer, -any in-mask voxel adjacent to at least one out-of-mask voxel is removed -from the mask. -} -\details{ -Diagonal voxels are not considered adjacent, i.e. the voxel at (0,0,0) is not -adjacent to the voxel at (1,1,0) or (1,1,1), although it is adjacent to -(1,0,0). -} diff --git a/man/fft_detrend.Rd b/man/fft_detrend.Rd deleted file mode 100644 index b7839dc..0000000 --- a/man/fft_detrend.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DCT_FFT_detrending.R -\name{fft_detrend} -\alias{fft_detrend} -\title{FFT detrending} -\usage{ -fft_detrend(X, N) -} -\arguments{ -\item{X}{\eqn{T \times V} numeric matrix. Each column is a voxel or vertex -time series.} - -\item{N}{Number of FFT entries to remove from each end of the vector} -} -\value{ -Detrended \code{X} -} -\description{ -FFT detrending -} -\keyword{internal} diff --git a/man/format_data.Rd b/man/format_data.Rd deleted file mode 100644 index bd09f0d..0000000 --- a/man/format_data.Rd +++ /dev/null @@ -1,98 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/format_data.R -\name{format_data} -\alias{format_data} -\title{Format data for pscrub and CompCor} -\usage{ -format_data( - X, - ROI_data = "infer", - ROI_noise = NULL, - noise_nPC = 5, - noise_erosion = NULL -) -} -\arguments{ -\item{X}{Wide numeric data matrix (\eqn{T observations} by \eqn{V variables}, \eqn{T << V}). -For example, if \code{X} represents an fMRI run, \eqn{T} should be the number -of timepoints and \eqn{V} should be the number of brainordinate vertices/voxels. - -Or, a 4D array or NIFTI or file path to a NIFTI (\eqn{I} by \eqn{J} by \eqn{K} by \eqn{T} -observations), in which case \code{ROI_data} must be provided. -(The vectorized data will be \eqn{T timepoints} by \eqn{V_{in-mask} voxels}) - -Or, a \code{ciftiTools} \code{"xifti"} object or a file path to a CIFTI -(The vectorized data will be \eqn{T timepoints} by \eqn{V_{left+right+sub} greyordinates}).} - -\item{ROI_data}{Indicates the data ROI. Allowed arguments depend on \code{X}: - -If \code{X} is a matrix, this must be a length \eqn{V} logical vector, where -the data ROI is indicated by \code{TRUE} values. If \code{"infer"} (default), all -columns of \code{X} will be included in the data ROI (\code{rep(TRUE, V)}). - -If \code{X} is an array or NIFTI, this must be either a vector of values -to expect for out-of-mask voxels in \code{X}, or a (file path to a) 3D NIFTI. -In the latter case, each of the volume dimensions should match the first -three dimensions of \code{X}. Voxels in the data ROI should be indicated by -\code{TRUE} and all other voxels by \code{FALSE}. If \code{"infer"} (default), -will be set to \code{c(0, NA, NaN)} (include all voxels which are not constant -\code{0}, \code{NA}, or \code{NaN}). - -If \code{X} is a \code{"xifti"} this must be the \code{brainstructures} -argument to \code{\link[ciftiTools]{read_cifti}}. If \code{"infer"} (default), -\code{brainstructures} will be set to \code{"all"} (use both left and right -cortex vertices, and subcortical voxels). - -If \code{NULL}, the data ROI will be empty. This is useful for obtaining just -the noise ROI, if the data and noise are located in separate files.} - -\item{ROI_noise}{Indicates the noise ROIs for aCompCor. Should be a list where -each entry corresponds to a distinct noise ROI. The names of the list should -be the ROI names, e.g. \code{"white_matter"} and \code{"csf"}. The expected -formats of the list entries depends on \code{X}: - -For all types of \code{X}, \code{ROI_noise} entries can be a matrix of noise -ROI data. The matrix should have \eqn{T} rows, with each column being a -data location's timeseries. - -If \code{X} is a matrix, entries can also indicate a noise ROI within \code{X}. -These entries must be a length \eqn{V} logical vector with \code{TRUE} values -indicating locations in \code{X} within that noise ROI. Since the ROIs must -not overlap, the masks must be mutually exclusive with each other, and with -\code{ROI_data}. - -If \code{X} is an array or NIFTI, entries can also indicate a noise ROI within \code{X}. -These entries must be a logical array or (file path to) a 3D NIFTI with the -same spatial dimensions as \code{X}, and with \code{TRUE} values indicating -voxels inside the noise ROI. Since the ROIs must not overlap, the masks must -be mutually exclusive with each other, and with \code{ROI_data}. - -(If \code{X} is a \code{"xifti"}, entries must be data matrices, since no -greyordinate locations in \code{X} are appropriate noise ROIs).} - -\item{noise_nPC}{The number of principal components to compute for each noise -ROI. Alternatively, values between 0 and 1, in which case they will -represent the minimum proportion of variance explained by the PCs used for -each noise ROI. The smallest number of PCs will be used to achieve this -proportion of variance explained. - -Should be a list or numeric vector with the same length as \code{ROI_noise}. -It will be matched to each ROI based on the name of each entry, or if the -names are missing, the order of entries. If it is an unnamed vector, its -elements will be recycled. Default: \code{5} (compute the top 5 PCs for -each noise ROI).} - -\item{noise_erosion}{The number of voxel layers to erode the noise ROIs by. -Should be a list or numeric vector with the same length as \code{ROI_noise}. -It will be matched to each ROI based on the name of each entry, or if the -names are missing, the order of entries. If it is an unnamed vector, its -elements will be recycled. Default: \code{NULL}, which will use a value of -0 (do not erode the noise ROIs).} -} -\value{ -A list with components "X", "X_noise", "ROI_data", and "ROI_noise" -} -\description{ -Format data for pscrub and CompCor -} -\keyword{internal} diff --git a/man/fsl_bptf.Rd b/man/fsl_bptf.Rd deleted file mode 100644 index 2544138..0000000 --- a/man/fsl_bptf.Rd +++ /dev/null @@ -1,32 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/fsl_bptf.R -\name{fsl_bptf} -\alias{fsl_bptf} -\title{\code{bptf} function from FSL} -\usage{ -fsl_bptf(orig_data, HP_sigma = 2000) -} -\arguments{ -\item{orig_data}{\eqn{T} by \eqn{V} data matrix whose columns will be detrended} - -\item{HP_sigma}{The frequency parameter for the highpass filter} -} -\value{ -The data with detrended columns -} -\description{ -Copy of highpass filter as implemented by \code{bptf} in FSL. The results are very -similar but not identical. -} -\details{ -Sources: -https://cpb-us-w2.wpmucdn.com/sites.udel.edu/dist/7/4542/files/2016/09/fsl_temporal_filt-15sywxn.m -https://github.com/rordenlab/niimath/blob/master/src/core32.c -} -\section{References}{ - -\itemize{ -\item{Jenkinson, M., Beckmann, C. F., Behrens, T. E. J., Woolrich, M. W. & Smith, S. M. FSL. NeuroImage 62, 782-790 (2012).} -} -} - diff --git a/man/get_NIFTI_ROI_masks.Rd b/man/get_NIFTI_ROI_masks.Rd deleted file mode 100644 index 066c0aa..0000000 --- a/man/get_NIFTI_ROI_masks.Rd +++ /dev/null @@ -1,41 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/CompCor_HCP.R -\name{get_NIFTI_ROI_masks} -\alias{get_NIFTI_ROI_masks} -\title{Get NIFTI ROI masks} -\usage{ -get_NIFTI_ROI_masks(nii_labels, ROI_noise = c("wm_cort", "csf")) -} -\arguments{ -\item{nii_labels}{\eqn{I} by \eqn{J} by \eqn{K} -NIFTI object or array (or file path to the NIFTI) which -contains the corresponding labels to each voxel in \code{nii}. Values should -be according to this table: -https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/AnatomicalROI/FreeSurferColorLUT . -In the HCP, the corresponding file is "ROIs/Atlas_wmparc.2.nii.gz".} - -\item{ROI_noise}{A list of numeric vectors. Each entry should represent labels -in \code{nii_labels} belonging to a single noise ROI, named by that entry's -name. Or, this can be a character vector of at least one of the following: -\code{"wm_cort"} (cortical white matter), \code{"wm_cblm"} (cerebellar white -matter), \code{"csf"} (cerebrospinal fluid). In the latter case, these labels -will be used: - -\describe{ -\item{\code{"wm_cort"}}{\code{c(3000:4035, 5001, 5002)}} -\item{\code{"wm_cblm"}}{\code{c(7, 46)}} -\item{\code{"csf"}}{\code{c(4, 5, 14, 15, 24, 31, 43, 44, 63, 250, 251, 252, 253, 254, 255))}} -} - -These default ROIs are based on this forum post: -https://www.mail-archive.com/hcp-users@humanconnectome.org/msg00931.html - -Default: \code{c("wm_cort", "csf")}} -} -\value{ -The ROIs -} -\description{ -Get NIFTI ROI masks -} -\keyword{internal} diff --git a/man/hat_matrix.Rd b/man/hat_matrix.Rd deleted file mode 100644 index 3593cb5..0000000 --- a/man/hat_matrix.Rd +++ /dev/null @@ -1,18 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/nuisance_regression.R -\name{hat_matrix} -\alias{hat_matrix} -\title{Hat matrix} -\usage{ -hat_matrix(design) -} -\arguments{ -\item{design}{The \eqn{T} by \eqn{Q} design matrix} -} -\value{ -The \eqn{T} by \eqn{T} hat matrix -} -\description{ -Get the hat matrix from a design matrix -} -\keyword{internal} diff --git a/man/infer_BOLD_format.Rd b/man/infer_BOLD_format.Rd deleted file mode 100644 index 5771fc3..0000000 --- a/man/infer_BOLD_format.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{infer_BOLD_format} -\alias{infer_BOLD_format} -\title{Infer fMRI data format} -\usage{ -infer_BOLD_format(BOLD, verbose = FALSE) -} -\arguments{ -\item{BOLD}{The fMRI data} - -\item{verbose}{Print the format? Default: \code{FALSE}.} -} -\value{ -The format: \code{"CIFTI"} file path, \code{"xifti"} object, -\code{"NIFTI"} file path, \code{"nifti"} object, or \code{"data"}. -} -\description{ -Infer fMRI data format -} -\keyword{internal} diff --git a/man/is_constant.Rd b/man/is_constant.Rd deleted file mode 100644 index 8f4e9a9..0000000 --- a/man/is_constant.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{is_constant} -\alias{is_constant} -\title{Is this numeric vector constant?} -\usage{ -is_constant(x, TOL = 1e-08) -} -\arguments{ -\item{x}{The numeric vector} - -\item{TOL}{minimum range of \code{x} to be considered non-constant. -Default: \code{1e-8}} -} -\value{ -Is \code{x} constant? -} -\description{ -Is this numeric vector constant? -} -\keyword{internal} diff --git a/man/is_integer.Rd b/man/is_integer.Rd deleted file mode 100644 index f0549a0..0000000 --- a/man/is_integer.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{is_integer} -\alias{is_integer} -\title{Is this an integer?} -\usage{ -is_integer(x, nneg = FALSE) -} -\arguments{ -\item{x}{The putative integer} - -\item{nneg}{Require \code{x>=0} (non-negative) too?} -} -\value{ -Logical indicating whether \code{x} is an integer -} -\description{ -Is this an integer? -} -\keyword{internal} diff --git a/man/neighbor_counts.Rd b/man/neighbor_counts.Rd deleted file mode 100644 index 50f8cdf..0000000 --- a/man/neighbor_counts.Rd +++ /dev/null @@ -1,22 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/format_data.R -\name{neighbor_counts} -\alias{neighbor_counts} -\title{Count each voxel's neighbors with value \code{TRUE}} -\usage{ -neighbor_counts(arr, pad = FALSE) -} -\arguments{ -\item{arr}{The 3D logical array.} - -\item{pad}{Pad value for edge.} -} -\value{ -An array with the same dimensions as \code{arr}. Each voxel value -will be the number of its immediate neighbors (0 to 6) which are \code{TRUE}. -} -\description{ -For each voxel in a 3D logical array, count the number of immediate neighbors -with value \code{TRUE}. -} -\keyword{internal} diff --git a/man/nuisance_regression.Rd b/man/nuisance_regression.Rd deleted file mode 100644 index 4dd8c1d..0000000 --- a/man/nuisance_regression.Rd +++ /dev/null @@ -1,20 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/nuisance_regression.R -\name{nuisance_regression} -\alias{nuisance_regression} -\title{Nuisance regression} -\usage{ -nuisance_regression(Y, design) -} -\arguments{ -\item{Y}{The \eqn{T} by \eqn{V} or \eqn{V} by \eqn{T} data.} - -\item{design}{The \eqn{T} by \eqn{Q} matrix of nuisance regressors.} -} -\value{ -The data after nuisance regression. -} -\description{ -Performs nuisance regression. The data and design matrix must both be -centered, or an intercept must be included in the design matrix! -} diff --git a/man/pct_sig.Rd b/man/pct_sig.Rd deleted file mode 100644 index b4aabc9..0000000 --- a/man/pct_sig.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DVARS.R -\name{pct_sig} -\alias{pct_sig} -\title{Convert data values to percent signal.} -\usage{ -pct_sig(X, center = median, by = c("column", "all")) -} -\arguments{ -\item{X}{a \eqn{T} by \eqn{N} numeric matrix. The columns will be normalized to -percent signal.} - -\item{center}{A function that computes the center of a numeric vector. -Default: \code{median}. Other common options include \code{mean} and -\code{mode}.} - -\item{by}{Should the center be measured individually for each \code{"column"} -(default), or should the center be the same across \code{"all"} columns?} -} -\value{ -\code{X} with its columns normalized to percent signal. (A value of -85 will represent a -15\% signal change.) -} -\description{ -Convert data values to percent signal. -} diff --git a/man/pscrub.Rd b/man/pscrub.Rd index 522d9b5..6c2fa81 100644 --- a/man/pscrub.Rd +++ b/man/pscrub.Rd @@ -49,7 +49,7 @@ 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: -\code{cbind(1, dct_bases(nrow(x), 4), more_nuisance)}.} +\code{cbind(1, fMRItools::dct_bases(nrow(x), 4), more_nuisance)}.} \item{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: \code{TRUE}. diff --git a/man/pscrub_Params.Rd b/man/pscrub_Params.Rd index 02cdcf8..5d2f083 100644 --- a/man/pscrub_Params.Rd +++ b/man/pscrub_Params.Rd @@ -26,7 +26,7 @@ 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: -\code{cbind(1, dct_bases(nrow(x), 4), more_nuisance)}.} +\code{cbind(1, fMRItools::dct_bases(nrow(x), 4), more_nuisance)}.} \item{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: \code{TRUE}. diff --git a/man/pscrub_multi.Rd b/man/pscrub_multi.Rd index ddd5492..ed6dac5 100644 --- a/man/pscrub_multi.Rd +++ b/man/pscrub_multi.Rd @@ -68,7 +68,7 @@ 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: -\code{cbind(1, dct_bases(nrow(x), 4), more_nuisance)}.} +\code{cbind(1, fMRItools::dct_bases(nrow(x), 4), more_nuisance)}.} \item{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: \code{TRUE}. diff --git a/man/read_nifti.Rd b/man/read_nifti.Rd deleted file mode 100644 index 7c2c912..0000000 --- a/man/read_nifti.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{read_nifti} -\alias{read_nifti} -\title{Wrapper to common functions for reading NIFTIs} -\usage{ -read_nifti(nifti_fname) -} -\arguments{ -\item{nifti_fname}{The file name of the NIFTI.} -} -\value{ -The NIFTI -} -\description{ -Tries \code{RNifti::readNifti}, then \code{oro.nifti::readNIfTI}. If -neither package is available an error is raised. -} -\keyword{internal} diff --git a/man/scale_med.Rd b/man/scale_med.Rd deleted file mode 100644 index b15dc2e..0000000 --- a/man/scale_med.Rd +++ /dev/null @@ -1,28 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{scale_med} -\alias{scale_med} -\title{Robust scaling} -\usage{ -scale_med(mat, TOL = 1e-08, drop_const = TRUE) -} -\arguments{ -\item{mat}{A numerical matrix.} - -\item{TOL}{minimum MAD to consider a column non-constant. -Default: \code{1e-8}} - -\item{drop_const}{Drop} -} -\value{ -The input matrix with its columns centered and scaled. -} -\description{ -Centers and scales the columns of a matrix robustly -} -\details{ -Centers each column on its median, and scales each column by its median -absolute deviation (MAD). Constant-valued columns are set to \code{NA} -(or removed if \code{drop_const}) and a warning is raised. If all -MADs are zero, an error is raised. -} diff --git a/man/unmask_vol.Rd b/man/unmask_vol.Rd deleted file mode 100644 index 00d5900..0000000 --- a/man/unmask_vol.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/artifact_images.R -\name{unmask_vol} -\alias{unmask_vol} -\title{Undo a volumetric mask} -\usage{ -unmask_vol(dat, mask, fill = NA) -} -\arguments{ -\item{dat}{Data matrix with locations along the rows and measurements along -the columns. If only one set of measurements were made, this may be a -vector.} - -\item{mask}{Volumetric binary mask. \code{TRUE} indicates voxels inside the -mask.} - -\item{fill}{The value for locations outside the mask. Default: \code{NA}.} -} -\value{ -The 3D or 4D unflattened volume array -} -\description{ -Un-applies a mask to vectorized data to yield its volumetric representation. -The mask and data should have compatible dimensions: the number of rows in -\code{dat} should equal the number of locations within the \code{mask}. -} diff --git a/tests/testthat.R b/tests/testthat.R index ce55c84..2625306 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -5,6 +5,7 @@ my_wb <- "../../workbench" library(testthat) library(ciftiTools) if (interactive()) { ciftiTools.setOption("wb_path", my_wb) } +library(fMRItools) library(fMRIscrub) diff --git a/tests/testthat/test-test_all_methods.R b/tests/testthat/test-test_all_methods.R index 5d4bc79..ead9064 100644 --- a/tests/testthat/test-test_all_methods.R +++ b/tests/testthat/test-test_all_methods.R @@ -18,7 +18,7 @@ test_that("pscrub works", { plot(psx) psx <- pscrub( - matrix(rnorm(4000), nrow=30), "fusedPCA", 1, center=FALSE, PESEL=FALSE, kurt_quantile=.8, + matrix(rnorm(9000), nrow=30), "fusedPCA", 1, center=FALSE, PESEL=FALSE, kurt_quantile=.8, full_PCA = TRUE, cutoff=5, verbose=TRUE ) plot(psx) @@ -28,11 +28,11 @@ test_that("pscrub works", { )) psx <- testthat::expect_warning(pscrub( - matrix(rnorm(10000), nrow=100) + 100, nuisance=dct_bases(100, 2) + matrix(rnorm(10000), nrow=100) + 100, nuisance=fMRItools::dct_bases(100, 2) )) psx <- testthat::expect_warning(pscrub( - Dat2, projection="PCA", nuisance=cbind(1, dct_bases(nrow(Dat2), 12)), + Dat2, projection="PCA", nuisance=cbind(1, fMRItools::dct_bases(nrow(Dat2), 12)), comps_mean_dt=2, comps_var_dt=2, get_dirs=TRUE, get_outliers=FALSE )) plot(psx) diff --git a/tests/testthat/test-test_scrubbing.R b/tests/testthat/test-test_scrubbing.R index d07ebad..6b4298b 100644 --- a/tests/testthat/test-test_scrubbing.R +++ b/tests/testthat/test-test_scrubbing.R @@ -28,11 +28,11 @@ test_that("pscrub works", { )) psx <- testthat::expect_warning(pscrub( - matrix(rnorm(10000), nrow=100) + 100, nuisance=dct_bases(100, 2) + matrix(rnorm(10000), nrow=100) + 100, nuisance=fMRItools::dct_bases(100, 2) )) psx <- testthat::expect_warning(pscrub( - Dat2, projection="fusedPCA", nuisance=cbind(1, dct_bases(nrow(Dat2), 12)), + Dat2, projection="fusedPCA", nuisance=cbind(1, fMRItools::dct_bases(nrow(Dat2), 12)), comps_mean_dt=2, comps_var_dt=2, get_dirs=TRUE, get_outliers=FALSE )) plot(psx) @@ -70,5 +70,5 @@ test_that("ciftiTools-related functions work", { }) test_that("Miscellaneous functions work", { - testthat::expect_warning(summary(pscrub(fsl_bptf(Dat2)))) + testthat::expect_warning(summary(pscrub(fMRItools::fsl_bptf(Dat2)))) }) diff --git a/vignettes/projection_scrubbing.rmd b/vignettes/projection_scrubbing.rmd index b585653..601996e 100644 --- a/vignettes/projection_scrubbing.rmd +++ b/vignettes/projection_scrubbing.rmd @@ -91,7 +91,7 @@ Let's see what the problem was in the first scan. We can reconstruct the origina fname = system.file("extdata", "Dat1_mask.nii.gz", package = "fMRIscrub") Mask1 = readNIfTI(fname) > 0 #Pitt_0050048 (full of artifacts) Mask1 = array(Mask1, dim=c(dim(Mask1), 1)) # 2D --> 3D slice -Img1 = fMRIscrub::unmask_vol(t(Dat1), Mask1) +Img1 = fMRItools::unvec_vol(t(Dat1), Mask1) ``` We will compare the timepoint of median leverage (left) to the timepoint of maximum leverage (right). @@ -126,8 +126,8 @@ par(mfrow=c(1,2)) # Constant voxels are deleted during the `pscrub` algorithm, so the artifact images will have # missing values where the constant voxels were. -artImg1.mean = unmask_vol(t(artImg1$mean), Mask1) -artImg1.top = unmask_vol(t(artImg1$top), Mask1) +artImg1.mean = fMRItools::unvec_vol(t(artImg1$mean), Mask1) +artImg1.top = fMRItools::unvec_vol(t(artImg1$top), Mask1) idx = which(which(psx$outlier_flag) == t_max) image(artImg1.mean[,,1,idx], main=paste0('Lev image, mean (T=',t_max,')')) From 4909a0a70829336b4cf786679c3563182344b90a Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Tue, 24 Jan 2023 17:17:43 +0700 Subject: [PATCH 02/14] belong in prev commit --- R/DCT_FFT_detrending.R | 110 ----------------------------------------- R/bandstop_filter.R | 38 -------------- man/bandstop_filter.Rd | 33 ------------- 3 files changed, 181 deletions(-) delete mode 100644 R/DCT_FFT_detrending.R delete mode 100644 R/bandstop_filter.R delete mode 100644 man/bandstop_filter.Rd diff --git a/R/DCT_FFT_detrending.R b/R/DCT_FFT_detrending.R deleted file mode 100644 index afa19d0..0000000 --- a/R/DCT_FFT_detrending.R +++ /dev/null @@ -1,110 +0,0 @@ -#' Generate cosine bases for the DCT -#' -#' @param T_ Length of timeseries -#' @param n Number of cosine bases -#' -#' @return Matrix with cosine bases along columns -#' -#' @export -dct_bases <- function(T_, n){ - b <- matrix(NA, T_, n) - idx <- (seq(T_)-1)/(T_-1) - for (ii in seq(n)) { b[,ii] <- cos(idx*pi*ii) } - b -} - -#' DCT and frequency conversion -#' -#' Convert between number of DCT bases and Hz of highpass filter -#' -#' Provide either \code{n} or \code{f} to calculate the other. -#' -#' If only the total length of the scan is known, you can set that to \code{TR} -#' and use \code{T_=1}. -#' -#' \eqn{f = n / (2 * T_ * TR)} -#' -#' @param T_ Length of timeseries (number of timepoints) -#' @param TR TR of the fMRI scan, in seconds (the time between timepoints) -#' @param n Number of cosine bases -#' @param f Hz of highpass filter -#' -#' @return If \code{n} was provided, the highpass filter cutoff (Hz) is returned. -#' Otherwise, if \code{f} was provided, the number of cosine bases is returned. -#' The result should be rounded before passing to \code{\link{dct_bases}} -#' -#' @export -dct_convert <- function(T_, TR, n=NULL, f=NULL){ - stopifnot(xor(is.null(n), is.null(f))) - if (is.null(n)) { - return(f * 2 * T_ * TR) - } else if (is.null(f)) { - return(n / (2 * T_ * TR)) - } else { stop() } -} - -#' @rdname dct_convert -dct2Hz <- function(T_, TR, n){ - dct_convert(T_, TR, n=n) -} - -#' @rdname dct_convert -Hz2dct <- function(T_, TR, f){ - dct_convert(T_, TR, f=f) -} - -#' FFT detrending -#' -#' @param X \eqn{T \times V} numeric matrix. Each column is a voxel or vertex -#' time series. -#' @param N Number of FFT entries to remove from each end of the vector -#' -#' @return Detrended \code{X} -#' -#' @importFrom stats mvfft -#' -#' @keywords internal -fft_detrend <- function(X, N) { - T_ <- nrow(X) - Y <- mvfft(X) - Y[seq(N),] <- 0 - Y[seq(T_-N+1, T_),] <- 0 - Re(stats::mvfft(Y, inverse=TRUE)) -} - -#' Detrending with DCT or FFT -#' -#' @param X \eqn{T \times V} numeric matrix. Each column is a voxel or vertex -#' time series. -#' @param TR TR of the fMRI scan, in seconds (the time between timepoints) -#' @param f Hz of highpass filter. Default: \code{.008} -#' @param method \code{"DCT"} (default) or \code{"FFT"}. -#' -#' @return Detrended \code{X} -#' -#' @export -detrend <- function(X, TR, f=.008, method=c("DCT", "FFT")) { - X <- as.matrix(X) - is_pos_num <- function(q){ is.numeric(q) && length(q)==1 && q > 0 } - stopifnot(is_pos_num(TR)) - stopifnot(is_pos_num(f)) - method <- match.arg(method, c("DCT", "FFT")) - - T_ <- nrow(X) - N <- round(f * TR * T_) - - N <- switch(method, - DCT = round(f * TR * T_ * 2), - FFT = round(f * TR * T_) - ) - - if (N < 1) { return(t( scale(t(X), scale=FALSE)) ) } - if (N > ifelse(method=="DCT", T_, T_/2)) { - stop("Maximum `f` for this data is: ", round(1/(TR*2), digits=5), " Hz.") - } - - switch(method, - DCT = nuisance_regression(X, cbind(1, dct_bases(T_, N))), - FFT = fft_detrend(X, N) - ) -} \ No newline at end of file diff --git a/R/bandstop_filter.R b/R/bandstop_filter.R deleted file mode 100644 index 2f4c073..0000000 --- a/R/bandstop_filter.R +++ /dev/null @@ -1,38 +0,0 @@ -#' Bandstop filter -#' -#' Filter out frequencies within a given range using a Chebyshev Type II -#' stopband. Basically a convenience wrapper for the \code{gsignal::cheby2} -#' function. -#' -#' @param x A numeric matrix, with each column being a timeseries to apply -#' the stopband filter -#' @param TR the time difference between rows in x, in seconds. -#' @param f1,f2 The frequency limits for the filter, in Hz. \code{f1 < f2}. -#' @param Rs The amount of attenuation of the stopband ripple, in dB -#' -#' @return The filtered data -#' @export -#' -#' @examples -#' library(gsignal) -#' n_voxels = 1e4 -#' n_timepoints = 100 -#' X = cbind(arima.sim(n=100, list(ar=.6)), arima.sim(n=100, list(ar=.6))) -#' Y = bandstop_filter(X, .72, .31, .43) -bandstop_filter <- function(x, TR, f1, f2, Rs=20){ - - if (!requireNamespace("gsignal", quietly = TRUE)) { - stop("Package \"gsignal\" needed for bandstop filter.", call. = FALSE) - } - - if (is.vector(x)) { x <- as.matrix(x) } - stopifnot(nrow(x)>5) - is_pos_num <- function(q){ length(q)==1 && is.numeric(q) && q>0 } - stopifnot(is_pos_num(TR)) - stopifnot(is_pos_num(f1)); stopifnot(is_pos_num(f2)) - stopifnot(f1 < f2) - - nf <- c(f1, f2) * TR * 2 - nfilt <- gsignal::cheby2(2, Rs=Rs, w=nf, type="stop") - gsignal::filtfilt(nfilt, x) -} \ No newline at end of file diff --git a/man/bandstop_filter.Rd b/man/bandstop_filter.Rd deleted file mode 100644 index 1160c4d..0000000 --- a/man/bandstop_filter.Rd +++ /dev/null @@ -1,33 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/bandstop_filter.R -\name{bandstop_filter} -\alias{bandstop_filter} -\title{Bandstop filter} -\usage{ -bandstop_filter(x, TR, f1, f2, Rs = 20) -} -\arguments{ -\item{x}{A numeric matrix, with each column being a timeseries to apply -the stopband filter} - -\item{TR}{the time difference between rows in x, in seconds.} - -\item{f1, f2}{The frequency limits for the filter, in Hz. \code{f1 < f2}.} - -\item{Rs}{The amount of attenuation of the stopband ripple, in dB} -} -\value{ -The filtered data -} -\description{ -Filter out frequencies within a given range using a Chebyshev Type II -stopband. Basically a convenience wrapper for the \code{gsignal::cheby2} -function. -} -\examples{ -library(gsignal) -n_voxels = 1e4 -n_timepoints = 100 -X = cbind(arima.sim(n=100, list(ar=.6)), arima.sim(n=100, list(ar=.6))) -Y = bandstop_filter(X, .72, .31, .43) -} From 10dbc6575f3186b76092cb85e43f0b09a6ce1237 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Tue, 24 Jan 2023 18:57:12 +0700 Subject: [PATCH 03/14] Grey --> Gray --- R/rox_args_docs.R | 4 ++-- man/data_CompCor_Params.Rd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/rox_args_docs.R b/R/rox_args_docs.R index c2815f4..3b80143 100644 --- a/R/rox_args_docs.R +++ b/R/rox_args_docs.R @@ -99,7 +99,7 @@ NULL #' (The vectorized data will be \eqn{T timepoints} by \eqn{V_{in-mask} voxels}) #' #' Or, a \code{ciftiTools} \code{"xifti"} object or a file path to a CIFTI -#' (The vectorized data will be \eqn{T timepoints} by \eqn{V_{left+right+sub} greyordinates}). +#' (The vectorized data will be \eqn{T timepoints} by \eqn{V_{left+right+sub} grayordinates}). #' @param ROI_data Indicates the data ROI. Allowed arguments depend on \code{X}: #' #' If \code{X} is a matrix, this must be a length \eqn{V} logical vector, where @@ -143,7 +143,7 @@ NULL #' be mutually exclusive with each other, and with \code{ROI_data}. #' #' (If \code{X} is a \code{"xifti"}, entries must be data matrices, since no -#' greyordinate locations in \code{X} are appropriate noise ROIs). +#' grayordinate locations in \code{X} are appropriate noise ROIs). #' @name data_CompCor_Params #' @keywords internal NULL diff --git a/man/data_CompCor_Params.Rd b/man/data_CompCor_Params.Rd index 65a661f..0655f69 100644 --- a/man/data_CompCor_Params.Rd +++ b/man/data_CompCor_Params.Rd @@ -13,7 +13,7 @@ observations), in which case \code{ROI_data} must be provided. (The vectorized data will be \eqn{T timepoints} by \eqn{V_{in-mask} voxels}) Or, a \code{ciftiTools} \code{"xifti"} object or a file path to a CIFTI -(The vectorized data will be \eqn{T timepoints} by \eqn{V_{left+right+sub} greyordinates}).} +(The vectorized data will be \eqn{T timepoints} by \eqn{V_{left+right+sub} grayordinates}).} \item{ROI_data}{Indicates the data ROI. Allowed arguments depend on \code{X}: @@ -59,7 +59,7 @@ voxels inside the noise ROI. Since the ROIs must not overlap, the masks must be mutually exclusive with each other, and with \code{ROI_data}. (If \code{X} is a \code{"xifti"}, entries must be data matrices, since no -greyordinate locations in \code{X} are appropriate noise ROIs).} +grayordinate locations in \code{X} are appropriate noise ROIs).} } \description{ fMRI data for \code{scrub} and \code{CompCor} From cf42a132d0d3bc3b3c7276c134bb9f1b89ca3e4e Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 15:26:54 +0700 Subject: [PATCH 04/14] app veyor rem and cleanup --- .github/workflows/R-CMD-check.yaml | 53 +++++----------------------- README.Rmd | 1 - appveyor.yml | 56 ------------------------------ cran-comments.md | 36 +++---------------- 4 files changed, 13 insertions(+), 133 deletions(-) delete mode 100644 appveyor.yml diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index a5418a5..1d8c68c 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -1,4 +1,4 @@ -# Workflow derived from https://github.com/r-lib/actions/tree/master/examples +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: @@ -8,57 +8,22 @@ on: name: R-CMD-check -jobs: +jobs: R-CMD-check: - runs-on: ${{ matrix.config.os }} - - name: ${{ matrix.config.os }} (${{ matrix.config.r }}) - - strategy: - fail-fast: false - matrix: - config: - - {os: macOS-latest, r: 'release'} - - {os: windows-latest, r: 'release'} - - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} - - {os: ubuntu-latest, r: 'release'} - - {os: ubuntu-latest, r: 'oldrel-1'} - + runs-on: ubuntu-latest env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} R_KEEP_PKG_SOURCE: yes - steps: - - uses: actions/checkout@v2 - - - uses: r-lib/actions/setup-pandoc@v1 + - uses: actions/checkout@v3 - - uses: r-lib/actions/setup-r@v1 + - uses: r-lib/actions/setup-r@v2 with: - r-version: ${{ matrix.config.r }} - http-user-agent: ${{ matrix.config.http-user-agent }} use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@v1 + - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: rcmdcheck + extra-packages: any::rcmdcheck + needs: check - - name: Install dependencies - run: | - install.packages(c("remotes","testthat"),dependencies=TRUE) - remotes::install_github("statsmaths/glmgen/R_pkg/glmgen") - shell: Rscript {0} - - - uses: r-lib/actions/check-r-package@v1 - - - name: Show testthat output - if: always() - run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true - shell: bash - - - name: Upload check results - if: failure() - uses: actions/upload-artifact@main - with: - name: ${{ runner.os }}-r${{ matrix.config.r }}-results - path: check + - uses: r-lib/actions/check-r-package@v2 diff --git a/README.Rmd b/README.Rmd index d6c3238..d11ef28 100644 --- a/README.Rmd +++ b/README.Rmd @@ -8,7 +8,6 @@ output: github_document [![R-CMD-check](https://github.com/mandymejia/fMRIscrub/workflows/R-CMD-check/badge.svg)](https://github.com/mandymejia/fMRIscrub/actions) -[![AppVeyor build status](https://ci.appveyor.com/api/projects/status/github/mandymejia/fMRIscrub?branch=master&svg=true)](https://ci.appveyor.com/project/mandymejia/fMRIscrub) [![Codecov test coverage](https://codecov.io/gh/mandymejia/fMRIscrub/branch/master/graph/badge.svg)](https://app.codecov.io/gh/mandymejia/fMRIscrub?branch=master) diff --git a/appveyor.yml b/appveyor.yml deleted file mode 100644 index 9d152bb..0000000 --- a/appveyor.yml +++ /dev/null @@ -1,56 +0,0 @@ -# DO NOT CHANGE the "init" and "install" sections below - -# Download script file from GitHub -init: - ps: | - $ErrorActionPreference = "Stop" - Invoke-WebRequest http://raw.github.com/krlmlr/r-appveyor/master/scripts/appveyor-tool.ps1 -OutFile "..\appveyor-tool.ps1" - Import-Module '..\appveyor-tool.ps1' - -install: - ps: Bootstrap - -cache: - - C:\RLibrary - -environment: - global: - NOT_CRAN: true - USE_RTOOLS: true - R_REMOTES_NO_ERRORS_FROM_WARNINGS: true - WARNINGS_ARE_ERRORS: 1 - -# Adapt as necessary starting from here - -build_script: - - travis-tool.sh install_deps - -test_script: - - travis-tool.sh run_tests - -on_failure: - - 7z a failure.zip *.Rcheck\* - - appveyor PushArtifact failure.zip - -artifacts: - - path: '*.Rcheck\**\*.log' - name: Logs - - - path: '*.Rcheck\**\*.out' - name: Logs - - - path: '*.Rcheck\**\*.fail' - name: Logs - - - path: '*.Rcheck\**\*.Rout' - name: Logs - - - path: '\*_*.tar.gz' - name: Bits - - - path: '\*_*.zip' - name: Bits - -branches: - only: - - master \ No newline at end of file diff --git a/cran-comments.md b/cran-comments.md index bd8c960..2bd9367 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,48 +1,20 @@ ## Test environments -* Linux, R 4.0.4 -* Windows 10 x64, R 4.1.1 +* Windows x86_64-w64-mingw32/x64, R 4.2.2 +* Mac x86_64-apple-darwin17.0, R 4.1.1 ## R CMD check results -> checking compiled code ... NOTE - Note: information on .o files for x64 is not available - File 'C:/Users/damon/Desktop/fMRI/fMRIscrub.Rcheck/fMRIscrub/libs/x64/fMRIscrub.dll': - Found 'abort', possibly from 'abort' (C), 'runtime' (Fortran) - Found 'exit', possibly from 'exit' (C), 'stop' (Fortran) - Found 'printf', possibly from 'printf' (C) - - Compiled code should not call entry points which might terminate R nor - write to stdout/stderr instead of to the console, nor use Fortran I/O - nor system RNGs. The detected symbols are linked into the code but - might come from libraries and not actually be called. - - See 'Writing portable packages' in the 'Writing R Extensions' manual. - WARNING - 'qpdf' is needed for checks on size reduction of PDFs - -> checking for old-style vignette sources ... NOTE - Vignette sources only in 'inst/doc': - 'projection_scrubbing.rmd' - A 'vignettes' directory is required as from R 3.1.0 - and these will not be indexed nor checked - -0 errors v | 1 warning x | 2 notes x +0 errors ✔ | 0 warnings ✔ | 0 notes ✔ s -------------------------------------------------------------------------------- We indeed use a "vignettes" directory. -`abort`, `exit` and `printf` is never called by our code. - ## Downstream dependencies -`templateICAr` works fine with this new version. +None. ## Tests Passes all the tests in `tests/testthat.R`. - -## First submisison - -This is the first submission since `fMRIscrub` was removed from CRAN on May 2, 2022 due to a package in Suggests being no longer available. The current submission resolves this issue. \ No newline at end of file From 44ffe2c167b6c11f85fe487c59892d85588ee1f2 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 17:16:28 +0700 Subject: [PATCH 05/14] Tidy up for CRAN --- DESCRIPTION | 2 +- NEWS.md | 5 +++++ README.md | 2 -- cran-comments.md | 6 +----- 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 69c31b5..211c3cf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: fMRIscrub Type: Package Title: Scrubbing and Other Data Cleaning Routines for fMRI -Version: 0.11.1 +Version: 0.12.0 Authors@R: c( person(given = "Amanda", family = "Mejia", diff --git a/NEWS.md b/NEWS.md index e61e428..a9524e2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# 0.12.0 + +* Use `fMRItools` where possible and remove duplicated functions from here. +* Remove AppVeyor + # 0.8.0 * Renamed `clever` to `projection scrubbing` diff --git a/README.md b/README.md index 694b0bd..46fe576 100644 --- a/README.md +++ b/README.md @@ -6,8 +6,6 @@ [![R-CMD-check](https://github.com/mandymejia/fMRIscrub/workflows/R-CMD-check/badge.svg)](https://github.com/mandymejia/fMRIscrub/actions) -[![AppVeyor build -status](https://ci.appveyor.com/api/projects/status/github/mandymejia/fMRIscrub?branch=master&svg=true)](https://ci.appveyor.com/project/mandymejia/fMRIscrub) [![Codecov test coverage](https://codecov.io/gh/mandymejia/fMRIscrub/branch/master/graph/badge.svg)](https://app.codecov.io/gh/mandymejia/fMRIscrub?branch=master) diff --git a/cran-comments.md b/cran-comments.md index 2bd9367..5d6688e 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -5,11 +5,7 @@ ## R CMD check results -0 errors ✔ | 0 warnings ✔ | 0 notes ✔ s - --------------------------------------------------------------------------------- - -We indeed use a "vignettes" directory. +0 errors ✔ | 0 warnings ✔ | 0 notes ✔ ## Downstream dependencies From 8bf481bcad605c5025ff2de814c779cb08bf3722 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 17:17:18 +0700 Subject: [PATCH 06/14] Create CRAN-SUBMISSION --- CRAN-SUBMISSION | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 CRAN-SUBMISSION diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION new file mode 100644 index 0000000..31867f6 --- /dev/null +++ b/CRAN-SUBMISSION @@ -0,0 +1,3 @@ +Version: 0.12.0 +Date: 2023-01-25 10:17:14 UTC +SHA: 44ffe2c167b6c11f85fe487c59892d85588ee1f2 From 0a70a8fb4a47dd064c0b2292bc35783a69fd3859 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 21:44:10 +0700 Subject: [PATCH 07/14] rejected from CRAN (fusedPCA) --- CRAN-SUBMISSION | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 CRAN-SUBMISSION diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION deleted file mode 100644 index 31867f6..0000000 --- a/CRAN-SUBMISSION +++ /dev/null @@ -1,3 +0,0 @@ -Version: 0.12.0 -Date: 2023-01-25 10:17:14 UTC -SHA: 44ffe2c167b6c11f85fe487c59892d85588ee1f2 From dae688522d890bda5ed03457af6e79ff5a052acf Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 21:51:57 +0700 Subject: [PATCH 08/14] get rid of unnecessary Imports and Suggests --- DESCRIPTION | 4 ---- 1 file changed, 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 211c3cf..0ffc77b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -36,7 +36,6 @@ Depends: R (>= 3.5.0) License: GPL-3 Encoding: UTF-8 Imports: - expm, MASS, e1071, fMRItools (>= 0.2.2), @@ -54,11 +53,8 @@ Suggests: rmarkdown, RNifti, ggplot2, - ggpubr, - gsignal, fastICA, oro.nifti, - gridExtra, testthat (>= 3.0.0), covr Roxygen: list(markdown = TRUE) From 1de35d00342e548fb82f51005fa4ba9304afd5d0 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 22:07:12 +0700 Subject: [PATCH 09/14] Remove fusedPCA for CRAN --- .Rbuildignore | 1 + NAMESPACE | 1 - NEWS.md | 3 +- R/artifact_images.R | 12 +- R/pscrub.R | 44 ++++--- R/pscrub_from_multi.R | 26 ++-- R/pscrub_multi.R | 160 +++++++++++++------------ R/rox_args_docs.R | 18 +-- man/fusedPCA.Rd | 73 ----------- man/fusedPCA_check_kwargs.Rd | 28 ----- man/pscrub.Rd | 27 +---- man/pscrub_Params.Rd | 10 -- man/pscrub_multi.Rd | 26 ---- tests/testthat/test-test_all_methods.R | 10 +- tests/testthat/test-test_scrubbing.R | 17 +-- vignettes/projection_scrubbing.rmd | 6 +- 16 files changed, 164 insertions(+), 298 deletions(-) delete mode 100644 man/fusedPCA.Rd delete mode 100644 man/fusedPCA_check_kwargs.Rd diff --git a/.Rbuildignore b/.Rbuildignore index a19808e..e12b870 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -16,3 +16,4 @@ robdist_new.*\.R ^codecov\.yml$ ^\.github$ ^CRAN-SUBMISSION$ +^R/fusedPCA\.R$ diff --git a/NAMESPACE b/NAMESPACE index 62bd41e..c377a06 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,7 +15,6 @@ S3method(summary,scrub_projection) export(DVARS) export(FD) export(artifact_images) -export(fusedPCA) export(high_kurtosis) export(leverage) export(pscrub) diff --git a/NEWS.md b/NEWS.md index a9524e2..e5b59b0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,7 +1,8 @@ -# 0.12.0 +# 0.12.1 * Use `fMRItools` where possible and remove duplicated functions from here. * Remove AppVeyor +* Remove fusedPCA (glmgen is not on CRAN, and C++ has some problems) # 0.8.0 diff --git a/R/artifact_images.R b/R/artifact_images.R index 9576d78..57b6af6 100644 --- a/R/artifact_images.R +++ b/R/artifact_images.R @@ -49,12 +49,12 @@ artifact_images <- function(psx, idx=NULL, use_dt=TRUE){ stop("No directions. Run `pscrub` again with `get_dirs=TRUE`.") } V <- psx$PCA$V - } else if ("fusedPCA" %in% names(psx)) { - U <- psx$fusedPCA$U - if (!("V" %in% names(psx$PCA))) { - stop("No directions. Run `pscrub` again with `get_dirs=TRUE`.") - } - V <- psx$fusedPCA$V + # } else if ("fusedPCA" %in% names(psx)) { + # U <- psx$fusedPCA$U + # if (!("V" %in% names(psx$PCA))) { + # stop("No directions. Run `pscrub` again with `get_dirs=TRUE`.") + # } + # V <- psx$fusedPCA$V } else if ("ICA" %in% names(psx)) { U <- fMRItools::scale_med(psx$ICA$M) if (!("S" %in% names (psx$ICA))) { diff --git a/R/pscrub.R b/R/pscrub.R index 26f5006..523d85a 100644 --- a/R/pscrub.R +++ b/R/pscrub.R @@ -18,8 +18,8 @@ #' } #' #' @inheritParams pscrub_Params -#' @param projection One of the following: \code{"ICA"} (default), \code{"PCA"}, -#' or \code{"fusedPCA"}. +#' @param projection One of the following: \code{"ICA"} (default) or \code{"PCA"}. +# or \code{"fusedPCA"}. # @param R_true The \eqn{T} by \eqn{T} correlation matrix, if known. Used for the bootstrap # robust distance measure. #' @param PESEL Use \code{\link[pesel]{pesel}} to select the number of @@ -54,17 +54,17 @@ #' If PCA was not used, all entries except \code{nPCs_PESEL} and/or \code{nPCs_avgvar} will not be included, depending on which #' method(s) was used to select the number of PCs. #' } -#' \item{fusedPCA}{ -#' If fusedPCA was used, this will be a list with components: -#' \describe{ -#' \item{U}{The \eqn{T} by \eqn{Q} PC score matrix.} -#' \item{D}{The standard deviation of each PC.} -#' \item{V}{The \eqn{P} by \eqn{Q} PC directions matrix. Included only if \code{get_dirs}} -#' \item{highkurt}{The length \code{Q} logical vector indicating scores of high kurtosis.} -#' \item{U_dt}{Detrended components of \code{U}. Included only if components were mean- or variance-detrended.} -#' \item{highkurt}{The length \code{Q} logical vector indicating detrended scores of high kurtosis. Included only if components were mean- or variance-detrended.} -#' } -#' } +# \item{fusedPCA}{ +# If fusedPCA was used, this will be a list with components: +# \describe{ +# \item{U}{The \eqn{T} by \eqn{Q} PC score matrix.} +# \item{D}{The standard deviation of each PC.} +# \item{V}{The \eqn{P} by \eqn{Q} PC directions matrix. Included only if \code{get_dirs}} +# \item{highkurt}{The length \code{Q} logical vector indicating scores of high kurtosis.} +# \item{U_dt}{Detrended components of \code{U}. Included only if components were mean- or variance-detrended.} +# \item{highkurt}{The length \code{Q} logical vector indicating detrended scores of high kurtosis. Included only if components were mean- or variance-detrended.} +# } +# } #' \item{ICA}{ #' If ICA was used, this will be a list with components: #' \describe{ @@ -87,15 +87,24 @@ #' #' psx = pscrub(X) pscrub = function( - X, projection=c("ICA", "fusedPCA", "PCA"), + X, projection=c( + "ICA", + #"fusedPCA", + "PCA" + ), nuisance="DCT4", center=TRUE, scale=TRUE, comps_mean_dt=FALSE, comps_var_dt=FALSE, - PESEL=TRUE, kurt_quantile=.99, fusedPCA_kwargs=NULL, + PESEL=TRUE, kurt_quantile=.99, + #fusedPCA_kwargs=NULL, get_dirs=FALSE, full_PCA=FALSE, get_outliers=TRUE, cutoff=4, seed=0, verbose=FALSE){ - projection <- match.arg(projection, c("ICA", "fusedPCA", "PCA")) + projection <- match.arg(projection, c( + "ICA", + #"fusedPCA", + "PCA" + )) if (!PESEL) { projection <- paste0(projection, "2") } if (kurt_quantile > 0) { projection <- paste0(projection, "_kurt") @@ -108,7 +117,8 @@ pscrub = function( X=X, projection=projection, nuisance=nuisance, center=center, scale=scale, comps_mean_dt=comps_mean_dt, comps_var_dt=comps_var_dt, - kurt_quantile=kurt_quantile, fusedPCA_kwargs=fusedPCA_kwargs, + kurt_quantile=kurt_quantile, + #fusedPCA_kwargs=fusedPCA_kwargs, get_dirs=get_dirs, full_PCA=full_PCA, get_outliers=get_outliers, cutoff=cutoff, seed=seed, verbose=verbose diff --git a/R/pscrub_from_multi.R b/R/pscrub_from_multi.R index 97682f4..3c57cfe 100644 --- a/R/pscrub_from_multi.R +++ b/R/pscrub_from_multi.R @@ -36,19 +36,19 @@ pscrub_from_multi <- function(psx) { } } psx$PCA$nPCs_avgvar <- psx$PCA$nPCs_PESEL <- NULL - # For fusedPCA - } else if (!is.null(psx$fusedPCA)) { - if (nrow(psx$fusedPCA$U) != ncol(psx$fusedPCA$U)) { - psx$fusedPCA$U <- psx$fusedPCA$U[, seq(nComps), drop=FALSE] - psx$fusedPCA$D <- psx$fusedPCA$D[seq(nComps), drop=FALSE] - if ("V" %in% names(psx$fusedPCA)) { - psx$fusedPCA$V <- psx$fusedPCA$V[, seq(nComps), drop=FALSE] - } - if ("highkurt" %in% names(psx$fusedPCA)) { - psx$fusedPCA$highkurt <- psx$fusedPCA$highkurt[seq(nComps)] - } - } - psx$PCA <- NULL + # # For fusedPCA + # } else if (!is.null(psx$fusedPCA)) { + # if (nrow(psx$fusedPCA$U) != ncol(psx$fusedPCA$U)) { + # psx$fusedPCA$U <- psx$fusedPCA$U[, seq(nComps), drop=FALSE] + # psx$fusedPCA$D <- psx$fusedPCA$D[seq(nComps), drop=FALSE] + # if ("V" %in% names(psx$fusedPCA)) { + # psx$fusedPCA$V <- psx$fusedPCA$V[, seq(nComps), drop=FALSE] + # } + # if ("highkurt" %in% names(psx$fusedPCA)) { + # psx$fusedPCA$highkurt <- psx$fusedPCA$highkurt[seq(nComps)] + # } + # } + # psx$PCA <- NULL # For ICA } else if (!is.null(psx$ICA)) { if (nrow(psx$ICA$M) != ncol(psx$ICA$M)) { diff --git a/R/pscrub_multi.R b/R/pscrub_multi.R index ced947e..4900e06 100644 --- a/R/pscrub_multi.R +++ b/R/pscrub_multi.R @@ -12,10 +12,10 @@ #' \item{\code{"PCA_kurt"}}{PCA using the high-kurtosis PCs among the top \eqn{k}.} #' \item{\code{"PCA2"}}{PCA using the top \eqn{k2} PCs.} #' \item{\code{"PCA2_kurt"}}{PCA using the high-kurtosis PCs among the top \eqn{k2}.} -#' \item{\code{"fusedPCA"}}{fusedPCA using the top \eqn{k} fused PCs.} -#' \item{\code{"fusedPCA_kurt"}}{fusedPCA using the high-kurtosis fused PCs among the top \eqn{k}.} -#' \item{\code{"fusedPCA2"}}{fusedPCA using the top \eqn{k2} fused PCs.} -#' \item{\code{"fusedPCA2_kurt"}}{fusedPCA using the high-kurtosis fused PCs among the top \eqn{k2}.} +# \item{\code{"fusedPCA"}}{fusedPCA using the top \eqn{k} fused PCs.} +# \item{\code{"fusedPCA_kurt"}}{fusedPCA using the high-kurtosis fused PCs among the top \eqn{k}.} +# \item{\code{"fusedPCA2"}}{fusedPCA using the top \eqn{k2} fused PCs.} +# \item{\code{"fusedPCA2_kurt"}}{fusedPCA using the high-kurtosis fused PCs among the top \eqn{k2}.} #' \item{\code{"ICA"}}{ICA using the top \eqn{k} ICs.} #' \item{\code{"ICA_kurt"}}{ICA using the high-kurtosis ICs among the top \eqn{k}.} #' \item{\code{"ICA2"}}{ICA using the top \eqn{k2} ICs.} @@ -56,17 +56,17 @@ #' If PCA was not used, all entries except \code{nPCs_PESEL} and/or \code{nPCs_avgvar} will not be included, depending on which #' method(s) was used to select the number of PCs. #' } -#' \item{fusedPCA}{ -#' If fusedPCA was used, this will be a list with components: -#' \describe{ -#' \item{U}{The \eqn{T} by \eqn{Q} PC score matrix.} -#' \item{D}{The standard deviation of each PC.} -#' \item{V}{The \eqn{P} by \eqn{Q} PC directions matrix. Included only if \code{get_dirs}} -#' \item{highkurt}{The length \code{Q} logical vector indicating scores of high kurtosis.} -#' \item{U_dt}{Detrended components of \code{U}. Included only if components were mean- or variance-detrended.} -#' \item{highkurt}{The length \code{Q} logical vector indicating detrended scores of high kurtosis. Included only if components were mean- or variance-detrended.} -#' } -#' } +# \item{fusedPCA}{ +# If fusedPCA was used, this will be a list with components: +# \describe{ +# \item{U}{The \eqn{T} by \eqn{Q} PC score matrix.} +# \item{D}{The standard deviation of each PC.} +# \item{V}{The \eqn{P} by \eqn{Q} PC directions matrix. Included only if \code{get_dirs}} +# \item{highkurt}{The length \code{Q} logical vector indicating scores of high kurtosis.} +# \item{U_dt}{Detrended components of \code{U}. Included only if components were mean- or variance-detrended.} +# \item{highkurt}{The length \code{Q} logical vector indicating detrended scores of high kurtosis. Included only if components were mean- or variance-detrended.} +# } +# } #' \item{ICA}{ #' If ICA was used, this will be a list with components: #' \describe{ @@ -90,7 +90,7 @@ pscrub_multi = function( X, projection = "ICA_kurt", nuisance="DCT4", center=TRUE, scale=TRUE, comps_mean_dt=FALSE, comps_var_dt=FALSE, - kurt_quantile=.99, fusedPCA_kwargs=NULL, + kurt_quantile=.99, #fusedPCA_kwargs=NULL, get_dirs=FALSE, full_PCA=FALSE, get_outliers=TRUE, cutoff=4, seed=0, verbose=FALSE){ @@ -137,7 +137,7 @@ pscrub_multi = function( outlier_flag = list(), mask = rep(1, V0_), PCA = NULL, - fusedPCA = NULL, + #fusedPCA = NULL, ICA = NULL ) @@ -145,8 +145,16 @@ pscrub_multi = function( out$mask[!X_NA_mask][X_const_mask] <- -2 # `projection`---------------------------------------------------------------- - valid_projection_PESEL <- c("PCA", "PCA_kurt", "fusedPCA", "fusedPCA_kurt", "ICA", "ICA_kurt") - valid_projection_avgvar <- c("PCA2", "PCA2_kurt", "fusedPCA2", "fusedPCA2_kurt", "ICA2", "ICA2_kurt") + valid_projection_PESEL <- c( + "PCA", "PCA_kurt", + #"fusedPCA", "fusedPCA_kurt", + "ICA", "ICA_kurt" + ) + valid_projection_avgvar <- c( + "PCA2", "PCA2_kurt", + #"fusedPCA2", "fusedPCA2_kurt", + "ICA2", "ICA2_kurt" + ) valid_projection <- c(valid_projection_PESEL, valid_projection_avgvar) if ("all" %in% projection) { projection <- valid_projection @@ -204,14 +212,14 @@ pscrub_multi = function( comps_dt <- (comps_mean_dt > 0) || (comps_var_dt > 0) kurt_quantile <- as.numeric(kurt_quantile) stopifnot(kurt_quantile >= 0 && kurt_quantile <= 1) - if(!identical(fusedPCA_kwargs, NULL)){ - names(fusedPCA_kwargs) <- match.arg( - names(fusedPCA_kwargs), c("lambda", "niter_max", "TOL", "verbose"), - several.ok=TRUE) - if(length(names(fusedPCA_kwargs)) != length(unique(names(fusedPCA_kwargs)))){ - stop("Duplicate fusedPCA_kwargs were given.") - } - } + # if(!identical(fusedPCA_kwargs, NULL)){ + # names(fusedPCA_kwargs) <- match.arg( + # names(fusedPCA_kwargs), c("lambda", "niter_max", "TOL", "verbose"), + # several.ok=TRUE) + # if(length(names(fusedPCA_kwargs)) != length(unique(names(fusedPCA_kwargs)))){ + # stop("Duplicate fusedPCA_kwargs were given.") + # } + # } get_dirs <- as.logical(get_dirs); stopifnot(isTRUE(get_dirs) || isFALSE(get_dirs)) full_PCA <- as.logical(full_PCA); stopifnot(isTRUE(full_PCA) || isFALSE(full_PCA)) get_outliers <- as.logical(get_outliers); stopifnot(isTRUE(get_outliers) || isFALSE(get_outliers)) @@ -265,10 +273,10 @@ pscrub_multi = function( } # Compute PCA. - if ("PCA" %in% base_projection || "fusedPCA" %in% base_projection) { + if ("PCA" %in% base_projection){# || "fusedPCA" %in% base_projection) { if (verbose) { cat("Computing PCA.\n") } # Compute the SVD of X or XXt. - if (get_dirs || "fusedPCA" %in% base_projection) { + if (get_dirs) {# || "fusedPCA" %in% base_projection) { out$PCA <- tryCatch( { svd(X)[c("u", "d", "v")] }, error = function(cond) { @@ -329,44 +337,44 @@ pscrub_multi = function( } } - # Compute fusedPCA. - if ("fusedPCA" %in% base_projection) { - maxK_fusedPCA <- max(as.numeric(list( - fusedPCA = out$PCA$nPCs_PESEL, - fusedPCA_kurt = out$PCA$nPCs_PESEL, - fusedPCA2 = out$PCA$nPCs_avgvar, - fusedPCA2_kurt = out$PCA$nPCs_avgvar - )[projection[grepl("fusedPCA", projection)]])) - if (verbose) { cat("Computing fusedPCA.\n") } - out$fusedPCA <- do.call( - fusedPCA, - c( - list( - X=X, X.svd=out$PCA[c("U", "D", "V")], - K=maxK_fusedPCA, solve_directions=get_dirs - ), - fusedPCA_kwargs - ) - ) - out$fusedPCA$PC_exec_times <- NULL; out$fusedPCA$nItes <- NULL - # V matrix from PCA no longer needed. - if(!get_dirs){ out$PCA$V <- NULL } + # # Compute fusedPCA. + # if ("fusedPCA" %in% base_projection) { + # maxK_fusedPCA <- max(as.numeric(list( + # fusedPCA = out$PCA$nPCs_PESEL, + # fusedPCA_kurt = out$PCA$nPCs_PESEL, + # fusedPCA2 = out$PCA$nPCs_avgvar, + # fusedPCA2_kurt = out$PCA$nPCs_avgvar + # )[projection[grepl("fusedPCA", projection)]])) + # if (verbose) { cat("Computing fusedPCA.\n") } + # out$fusedPCA <- do.call( + # fusedPCA, + # c( + # list( + # X=X, X.svd=out$PCA[c("U", "D", "V")], + # K=maxK_fusedPCA, solve_directions=get_dirs + # ), + # fusedPCA_kwargs + # ) + # ) + # out$fusedPCA$PC_exec_times <- NULL; out$fusedPCA$nItes <- NULL + # # V matrix from PCA no longer needed. + # if(!get_dirs){ out$PCA$V <- NULL } - tf_const_mask <- apply(out$fusedPCA$u, 2, is_constant) - if(any(tf_const_mask)){ - warning( - "Warning: ", sum(tf_const_mask), " out of ", length(tf_const_mask), - " fused PC scores are zero-variance.\n" - ) - } - names(out$fusedPCA)[names(out$fusedPCA) %in% c("u", "d", "v")] <- toupper( - names(out$fusedPCA)[names(out$fusedPCA) %in% c("u", "d", "v")] - ) - } - # Identify which fused PCs have high kurtosis. - if (any(c("fusedPCA_kurt", "fusedPCA2_kurt") %in% projection)) { - out$fusedPCA$highkurt <- high_kurtosis(out$fusedPCA$U, kurt_quantile=kurt_quantile) - } + # tf_const_mask <- apply(out$fusedPCA$u, 2, is_constant) + # if(any(tf_const_mask)){ + # warning( + # "Warning: ", sum(tf_const_mask), " out of ", length(tf_const_mask), + # " fused PC scores are zero-variance.\n" + # ) + # } + # names(out$fusedPCA)[names(out$fusedPCA) %in% c("u", "d", "v")] <- toupper( + # names(out$fusedPCA)[names(out$fusedPCA) %in% c("u", "d", "v")] + # ) + # } + # # Identify which fused PCs have high kurtosis. + # if (any(c("fusedPCA_kurt", "fusedPCA2_kurt") %in% projection)) { + # out$fusedPCA$highkurt <- high_kurtosis(out$fusedPCA$U, kurt_quantile=kurt_quantile) + # } # Remove extra PCA information. if (!is.null(out$PCA$U) && !full_PCA) { @@ -420,12 +428,12 @@ pscrub_multi = function( ) out$PCA$highkurt_dt <- high_kurtosis(out$PCA$U_dt, kurt_quantile=kurt_quantile) } - if (!is.null(out$fusedPCA$U)) { - out$fusedPCA$U_dt <- apply( - out$fusedPCA$U, 2, rob_stabilize, center=comps_mean_dt, scale=comps_var_dt - ) - out$fusedPCA$highkurt_dt <- high_kurtosis(out$fusedPCA$U_dt, kurt_quantile=kurt_quantile) - } + # if (!is.null(out$fusedPCA$U)) { + # out$fusedPCA$U_dt <- apply( + # out$fusedPCA$U, 2, rob_stabilize, center=comps_mean_dt, scale=comps_var_dt + # ) + # out$fusedPCA$highkurt_dt <- high_kurtosis(out$fusedPCA$U_dt, kurt_quantile=kurt_quantile) + # } if (!is.null(out$ICA$M)) { out$ICA$M_dt <- apply( out$ICA$M, 2, rob_stabilize, center=comps_mean_dt, scale=comps_var_dt @@ -457,10 +465,10 @@ pscrub_multi = function( PCA_kurt = which(out$PCA[[highkurt_ii]][seq(out$PCA$nPCs_PESEL)]), PCA2 = seq(out$PCA$nPCs_avgvar), PCA2_kurt = which(out$PCA[[highkurt_ii]][seq(out$PCA$nPCs_avgvar)]), - fusedPCA = seq(out$PCA$nPCs_PESEL), - fusedPCA_kurt = which(out$fusedPCA[[highkurt_ii]][seq(out$PCA$nPCs_PESEL)]), - fusedPCA2 = seq(out$PCA$nPCs_avgvar), - fusedPCA2_kurt = which(out$fusedPCA[[highkurt_ii]][seq(out$PCA$nPCs_avgvar)]), + #fusedPCA = seq(out$PCA$nPCs_PESEL), + #fusedPCA_kurt = which(out$fusedPCA[[highkurt_ii]][seq(out$PCA$nPCs_PESEL)]), + #fusedPCA2 = seq(out$PCA$nPCs_avgvar), + #fusedPCA2_kurt = which(out$fusedPCA[[highkurt_ii]][seq(out$PCA$nPCs_avgvar)]), ICA = seq(out$PCA$nPCs_PESEL), ICA_kurt = which(out$ICA[[highkurt_ii]][seq(out$PCA$nPCs_PESEL)]), ICA2 = seq(out$PCA$nPCs_avgvar), diff --git a/R/rox_args_docs.R b/R/rox_args_docs.R index 3b80143..2afd6f1 100644 --- a/R/rox_args_docs.R +++ b/R/rox_args_docs.R @@ -63,15 +63,15 @@ #' We model each component as a length \eqn{T} vector of Normal iid random variables, #' for which the distribution of kurtosis values can be approximated. The #' quantile is estimated based on this distribution. -#' @param fusedPCA_kwargs Arguments to \code{\link{fusedPCA}} in list form. Valid -#' entries are: -#' -#' \describe{ -#' \item{lambda}{Trend-filtering parameter. Default: \code{5}.} -#' \item{niter_max}{Maximum number of iterations. Default: \code{1000}.} -#' \item{TOL}{Convergence tolerance parameter. Default: \code{1e-8}.} -#' \item{verbose}{Print updates? Default: \code{FALSE}.} -#' } +# @param fusedPCA_kwargs Arguments to \code{\link{fusedPCA}} in list form. Valid +# entries are: +# +# \describe{ +# \item{lambda}{Trend-filtering parameter. Default: \code{5}.} +# \item{niter_max}{Maximum number of iterations. Default: \code{1000}.} +# \item{TOL}{Convergence tolerance parameter. Default: \code{1e-8}.} +# \item{verbose}{Print updates? Default: \code{FALSE}.} +# } #' @param get_dirs Should the projection directions be returned? This is the #' \eqn{V} matrix in PCA and \eqn{S} matrix in ICA. The default is \code{FALSE} #' to save memory. However, \code{get_dirs==TRUE} is required for \code{\link{artifact_images}}. diff --git a/man/fusedPCA.Rd b/man/fusedPCA.Rd deleted file mode 100644 index 3496153..0000000 --- a/man/fusedPCA.Rd +++ /dev/null @@ -1,73 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/fusedPCA.R -\name{fusedPCA} -\alias{fusedPCA} -\title{Fused PCA} -\usage{ -fusedPCA( - X, - X.svd = NULL, - solve_directions = TRUE, - K = NULL, - lambda = 0.5, - niter_max = 1000, - TOL = 1e-08, - verbose = FALSE -) -} -\arguments{ -\item{X}{A numerical data matrix (observations x variables).} - -\item{X.svd}{(Optional) The svd decomposition of X. Save time by providing -this argument if the svd has already been computed. Default NULL.} - -\item{solve_directions}{Should the principal directions be solved for? These -will be needed to display the artifact images for outlying observations.} - -\item{K}{(Optional) The number of fused PCs to solve for. If not -provided, it will be set to the number of regular PCs with variance above -the mean, up to 100 PCs.} - -\item{lambda}{The trend filtering parameter; roughly, the filtering intensity. -Default is 0.5 . Can be NULL (lets algorithm decide).} - -\item{niter_max}{The number of iterations to use for approximating the PC.} - -\item{TOL}{The maximum 2-norm between iterations to accept as convergence.} - -\item{verbose}{Print statements about convergence?} -} -\value{ -SVD The fused SVD decomposition of X (list with u, d, v). -} -\description{ -From: https://github.com/Lei-D/fusedPCA . Requires the \code{glmgen} package. -} -\section{References}{ - -\itemize{ -\item{Kim, S.-J., Koh, K., Boyd, S. & Gorinevsky, D. \$\\ell_1\$ Trend Filtering. SIAM Rev. 51, 339-360 (2009).} -\item{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).} -\item{Tibshirani, R. J. Adaptive piecewise polynomial estimation via trend filtering. The Annals of Statistics 42, 285-323 (2014).} -} -} - -\examples{ -U = matrix(rnorm(100*3),ncol=3) -U[20:23,1] = U[20:23,1] + 3 -U[40:43,2] = U[40:43,2] - 2 -U = svd(U)$u -D = diag(c(10,5,1)) -V = svd(matrix(rnorm(3*20),nrow=20))$u -X = U \%*\% D \%*\% t(V) -out3 = fusedPCA(X, K=3, lambda=.75) -matplot(out3$u, ty='l') -out3$d -plot(rowSums(out3$u^2), ty='l') - -# Orthonormalized -out3_svd = svd(out3$u) -matplot(out3_svd$u, ty='l') -out3_svd$d -plot(rowSums(out3_svd$u^2), ty='l') -} diff --git a/man/fusedPCA_check_kwargs.Rd b/man/fusedPCA_check_kwargs.Rd deleted file mode 100644 index daa8e96..0000000 --- a/man/fusedPCA_check_kwargs.Rd +++ /dev/null @@ -1,28 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/fusedPCA.R -\name{fusedPCA_check_kwargs} -\alias{fusedPCA_check_kwargs} -\title{Check fusedPCA args} -\usage{ -fusedPCA_check_kwargs( - X, - X.svd, - solve_directions, - K, - lambda, - niter_max, - TOL, - verbose -) -} -\arguments{ -\item{X, X.svd, solve_directions, K, lambda, niter_max, TOL, verbose}{See respective -functions.} -} -\value{ -\code{NULL}, invisibly. -} -\description{ -Check fusedPCA arguments for both R-core and cpp-core versions -} -\keyword{internal} diff --git a/man/pscrub.Rd b/man/pscrub.Rd index 6c2fa81..aaf2cd4 100644 --- a/man/pscrub.Rd +++ b/man/pscrub.Rd @@ -6,7 +6,7 @@ \usage{ pscrub( X, - projection = c("ICA", "fusedPCA", "PCA"), + projection = c("ICA", "PCA"), nuisance = "DCT4", center = TRUE, scale = TRUE, @@ -14,7 +14,6 @@ pscrub( comps_var_dt = FALSE, PESEL = TRUE, kurt_quantile = 0.99, - fusedPCA_kwargs = NULL, get_dirs = FALSE, full_PCA = FALSE, get_outliers = TRUE, @@ -29,8 +28,7 @@ pscrub( number of timepoints and \eqn{V} should be the number of vertices/voxels. Projection scrubbing will measure the outlyingness of each row in \code{X}.} -\item{projection}{One of the following: \code{"ICA"} (default), \code{"PCA"}, -or \code{"fusedPCA"}.} +\item{projection}{One of the following: \code{"ICA"} (default) or \code{"PCA"}.} \item{nuisance}{Nuisance signals to regress from each column of \code{X}. Should be specified as a design matrix: a \eqn{T} by \eqn{N} numeric matrix @@ -99,16 +97,6 @@ We model each component as a length \eqn{T} vector of Normal iid random variable for which the distribution of kurtosis values can be approximated. The quantile is estimated based on this distribution.} -\item{fusedPCA_kwargs}{Arguments to \code{\link{fusedPCA}} in list form. Valid -entries are: - -\describe{ -\item{lambda}{Trend-filtering parameter. Default: \code{5}.} -\item{niter_max}{Maximum number of iterations. Default: \code{1000}.} -\item{TOL}{Convergence tolerance parameter. Default: \code{1e-8}.} -\item{verbose}{Print updates? Default: \code{FALSE}.} -}} - \item{get_dirs}{Should the projection directions be returned? This is the \eqn{V} matrix in PCA and \eqn{S} matrix in ICA. The default is \code{FALSE} to save memory. However, \code{get_dirs==TRUE} is required for \code{\link{artifact_images}}.} @@ -156,17 +144,6 @@ where \code{Q} is the number of PCs selected by PESEL or of above-average varian If PCA was not used, all entries except \code{nPCs_PESEL} and/or \code{nPCs_avgvar} will not be included, depending on which method(s) was used to select the number of PCs. } -\item{fusedPCA}{ -If fusedPCA was used, this will be a list with components: -\describe{ -\item{U}{The \eqn{T} by \eqn{Q} PC score matrix.} -\item{D}{The standard deviation of each PC.} -\item{V}{The \eqn{P} by \eqn{Q} PC directions matrix. Included only if \code{get_dirs}} -\item{highkurt}{The length \code{Q} logical vector indicating scores of high kurtosis.} -\item{U_dt}{Detrended components of \code{U}. Included only if components were mean- or variance-detrended.} -\item{highkurt}{The length \code{Q} logical vector indicating detrended scores of high kurtosis. Included only if components were mean- or variance-detrended.} -} -} \item{ICA}{ If ICA was used, this will be a list with components: \describe{ diff --git a/man/pscrub_Params.Rd b/man/pscrub_Params.Rd index 5d2f083..1df8778 100644 --- a/man/pscrub_Params.Rd +++ b/man/pscrub_Params.Rd @@ -72,16 +72,6 @@ We model each component as a length \eqn{T} vector of Normal iid random variable for which the distribution of kurtosis values can be approximated. The quantile is estimated based on this distribution.} -\item{fusedPCA_kwargs}{Arguments to \code{\link{fusedPCA}} in list form. Valid -entries are: - -\describe{ -\item{lambda}{Trend-filtering parameter. Default: \code{5}.} -\item{niter_max}{Maximum number of iterations. Default: \code{1000}.} -\item{TOL}{Convergence tolerance parameter. Default: \code{1e-8}.} -\item{verbose}{Print updates? Default: \code{FALSE}.} -}} - \item{get_dirs}{Should the projection directions be returned? This is the \eqn{V} matrix in PCA and \eqn{S} matrix in ICA. The default is \code{FALSE} to save memory. However, \code{get_dirs==TRUE} is required for \code{\link{artifact_images}}.} diff --git a/man/pscrub_multi.Rd b/man/pscrub_multi.Rd index ed6dac5..c5ae44d 100644 --- a/man/pscrub_multi.Rd +++ b/man/pscrub_multi.Rd @@ -13,7 +13,6 @@ pscrub_multi( comps_mean_dt = FALSE, comps_var_dt = FALSE, kurt_quantile = 0.99, - fusedPCA_kwargs = NULL, get_dirs = FALSE, full_PCA = FALSE, get_outliers = TRUE, @@ -36,10 +35,6 @@ likely to contain outlier information. Choose at least one of the following: \item{\code{"PCA_kurt"}}{PCA using the high-kurtosis PCs among the top \eqn{k}.} \item{\code{"PCA2"}}{PCA using the top \eqn{k2} PCs.} \item{\code{"PCA2_kurt"}}{PCA using the high-kurtosis PCs among the top \eqn{k2}.} -\item{\code{"fusedPCA"}}{fusedPCA using the top \eqn{k} fused PCs.} -\item{\code{"fusedPCA_kurt"}}{fusedPCA using the high-kurtosis fused PCs among the top \eqn{k}.} -\item{\code{"fusedPCA2"}}{fusedPCA using the top \eqn{k2} fused PCs.} -\item{\code{"fusedPCA2_kurt"}}{fusedPCA using the high-kurtosis fused PCs among the top \eqn{k2}.} \item{\code{"ICA"}}{ICA using the top \eqn{k} ICs.} \item{\code{"ICA_kurt"}}{ICA using the high-kurtosis ICs among the top \eqn{k}.} \item{\code{"ICA2"}}{ICA using the top \eqn{k2} ICs.} @@ -114,16 +109,6 @@ We model each component as a length \eqn{T} vector of Normal iid random variable for which the distribution of kurtosis values can be approximated. The quantile is estimated based on this distribution.} -\item{fusedPCA_kwargs}{Arguments to \code{\link{fusedPCA}} in list form. Valid -entries are: - -\describe{ -\item{lambda}{Trend-filtering parameter. Default: \code{5}.} -\item{niter_max}{Maximum number of iterations. Default: \code{1000}.} -\item{TOL}{Convergence tolerance parameter. Default: \code{1e-8}.} -\item{verbose}{Print updates? Default: \code{FALSE}.} -}} - \item{get_dirs}{Should the projection directions be returned? This is the \eqn{V} matrix in PCA and \eqn{S} matrix in ICA. The default is \code{FALSE} to save memory. However, \code{get_dirs==TRUE} is required for \code{\link{artifact_images}}.} @@ -172,17 +157,6 @@ where \code{Q} is the number of PCs selected by PESEL or of above-average varian If PCA was not used, all entries except \code{nPCs_PESEL} and/or \code{nPCs_avgvar} will not be included, depending on which method(s) was used to select the number of PCs. } -\item{fusedPCA}{ -If fusedPCA was used, this will be a list with components: -\describe{ -\item{U}{The \eqn{T} by \eqn{Q} PC score matrix.} -\item{D}{The standard deviation of each PC.} -\item{V}{The \eqn{P} by \eqn{Q} PC directions matrix. Included only if \code{get_dirs}} -\item{highkurt}{The length \code{Q} logical vector indicating scores of high kurtosis.} -\item{U_dt}{Detrended components of \code{U}. Included only if components were mean- or variance-detrended.} -\item{highkurt}{The length \code{Q} logical vector indicating detrended scores of high kurtosis. Included only if components were mean- or variance-detrended.} -} -} \item{ICA}{ If ICA was used, this will be a list with components: \describe{ diff --git a/tests/testthat/test-test_all_methods.R b/tests/testthat/test-test_all_methods.R index ead9064..398d35c 100644 --- a/tests/testthat/test-test_all_methods.R +++ b/tests/testthat/test-test_all_methods.R @@ -17,11 +17,11 @@ test_that("pscrub works", { psx <- testthat::expect_warning(pscrub(Dat1)) plot(psx) - psx <- pscrub( - matrix(rnorm(9000), nrow=30), "fusedPCA", 1, center=FALSE, PESEL=FALSE, kurt_quantile=.8, - full_PCA = TRUE, cutoff=5, verbose=TRUE - ) - plot(psx) + # psx <- pscrub( + # matrix(rnorm(9000), nrow=30), "fusedPCA", 1, center=FALSE, PESEL=FALSE, kurt_quantile=.8, + # full_PCA = TRUE, cutoff=5, verbose=TRUE + # ) + # plot(psx) psx <- testthat::expect_warning(pscrub( matrix(rnorm(10000), ncol=50) diff --git a/tests/testthat/test-test_scrubbing.R b/tests/testthat/test-test_scrubbing.R index 6b4298b..2c30cc4 100644 --- a/tests/testthat/test-test_scrubbing.R +++ b/tests/testthat/test-test_scrubbing.R @@ -3,11 +3,14 @@ test_that("pscrub works", { Dat1, projection = "all" )) - testthat::expect_warning(fMRIscrub:::plot.scrub_projection_multi(psx)) + fMRIscrub:::plot.scrub_projection_multi(psx) # expect_warning w/ fusedPCA? psx <- testthat::expect_warning(fMRIscrub:::pscrub_multi( Dat2, - projection = c("fusedPCA", "ICA_kurt"), + projection = c( + #"fusedPCA", + "ICA_kurt" + ), kurt_quantile = .90, cutoff = 5, verbose = TRUE @@ -31,11 +34,11 @@ test_that("pscrub works", { matrix(rnorm(10000), nrow=100) + 100, nuisance=fMRItools::dct_bases(100, 2) )) - psx <- testthat::expect_warning(pscrub( - Dat2, projection="fusedPCA", nuisance=cbind(1, fMRItools::dct_bases(nrow(Dat2), 12)), - comps_mean_dt=2, comps_var_dt=2, get_dirs=TRUE, get_outliers=FALSE - )) - plot(psx) + # psx <- testthat::expect_warning(pscrub( + # Dat2, projection="fusedPCA", nuisance=cbind(1, fMRItools::dct_bases(nrow(Dat2), 12)), + # comps_mean_dt=2, comps_var_dt=2, get_dirs=TRUE, get_outliers=FALSE + # )) + # plot(psx) }) test_that("DVARS works", { diff --git a/vignettes/projection_scrubbing.rmd b/vignettes/projection_scrubbing.rmd index 601996e..3ee3b84 100644 --- a/vignettes/projection_scrubbing.rmd +++ b/vignettes/projection_scrubbing.rmd @@ -166,7 +166,11 @@ For testing purposes, we can try different projections at the same time using th # the default (ICA + kurtosis), fusedPCA + kurtosis, and ICA ps.Dat1.3 <- fMRIscrub:::pscrub_multi( - Dat1, projection=c("ICA_kurt", "fusedPCA_kurt", "ICA"), verbose=TRUE, comps_mean_dt=1, comps_var_dt=1 + Dat1, projection=c( + "ICA_kurt", + #"fusedPCA_kurt", + "ICA" + ), verbose=TRUE, comps_mean_dt=1, comps_var_dt=1 ) fMRIscrub:::plot.scrub_projection_multi(ps.Dat1.3, legend.position="bottom") ``` From cbbda93be454bdc46cfba98dd6c8ac4051b02b60 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 22:31:51 +0700 Subject: [PATCH 10/14] Update CRAN comments; iterate version --- DESCRIPTION | 2 +- cran-comments.md | 5 +++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0ffc77b..72b9e3d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: fMRIscrub Type: Package Title: Scrubbing and Other Data Cleaning Routines for fMRI -Version: 0.12.0 +Version: 0.12.1 Authors@R: c( person(given = "Amanda", family = "Mejia", diff --git a/cran-comments.md b/cran-comments.md index 5d6688e..0d39dd3 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -14,3 +14,8 @@ None. ## Tests Passes all the tests in `tests/testthat.R`. + +# Resubmission #1 + +The previous submission failed automatic checks because a suggested package is +not on CRAN (`glmgen`). Any references to this package have now been removed. \ No newline at end of file From 23625400304755e56b72f0233781a77b484d6da6 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 22:36:19 +0700 Subject: [PATCH 11/14] rm glmgen from DESC... lol --- DESCRIPTION | 1 - 1 file changed, 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 72b9e3d..63e790d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -48,7 +48,6 @@ Suggests: cowplot, ciftiTools, gifti, - glmgen, knitr, rmarkdown, RNifti, From 3dfcb041a0b941214ab064ad646aad9744e17314 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 22:37:26 +0700 Subject: [PATCH 12/14] make example data smaller --- R/pscrub.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/pscrub.R b/R/pscrub.R index 523d85a..29a844a 100644 --- a/R/pscrub.R +++ b/R/pscrub.R @@ -81,8 +81,8 @@ #' #' @examples #' library(fastICA) -#' n_voxels = 5e3 -#' n_timepoints = 70 +#' n_voxels = 2e3 +#' n_timepoints = 35 #' X = matrix(rnorm(n_timepoints*n_voxels), ncol = n_voxels) #' #' psx = pscrub(X) From f47ea59f215f7df7e60d29676018efdb5d282e67 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 22:37:57 +0700 Subject: [PATCH 13/14] rox --- man/pscrub.Rd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/man/pscrub.Rd b/man/pscrub.Rd index aaf2cd4..4fe360b 100644 --- a/man/pscrub.Rd +++ b/man/pscrub.Rd @@ -179,8 +179,8 @@ outline of the algorithm: \code{vignette("projection_scrubbing", package="fMRIsc \examples{ library(fastICA) -n_voxels = 5e3 -n_timepoints = 70 +n_voxels = 2e3 +n_timepoints = 35 X = matrix(rnorm(n_timepoints*n_voxels), ncol = n_voxels) psx = pscrub(X) From 1a45c47839569c231f33bff75b759210b91adce2 Mon Sep 17 00:00:00 2001 From: Damon Pham Date: Wed, 25 Jan 2023 22:43:05 +0700 Subject: [PATCH 14/14] Create CRAN-SUBMISSION --- CRAN-SUBMISSION | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 CRAN-SUBMISSION diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION new file mode 100644 index 0000000..2cbe70a --- /dev/null +++ b/CRAN-SUBMISSION @@ -0,0 +1,3 @@ +Version: 0.12.1 +Date: 2023-01-25 15:43:00 UTC +SHA: f47ea59f215f7df7e60d29676018efdb5d282e67