mandymejia / fMRIscrub / 199

pending completion

Build # 199

Build Type

Pull #31

travis-ci-com

Committed by

web-flow

Commit Message

<a href="https://github.com/mandymejia/fMRIscrub/commit/<a class=hub.com/mandymejia/fMRIscrub/commit/<a class="double-link" href="https://git"><a class=hub.com/mandymejia/fMRIscrub/commit/<a class="double-link" href="https://git"><a class=hub.com/mandymejia/fMRIscrub/commit/5116b0778a69b1279d4ab39561c03d4905f6848b">5116b0778">&lt;a href=&quot;https://github.com/mandymejia/fMRIscrub/commit/</a><a class="double-link" href="https://github.com/mandymejia/fMRIscrub/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/mandymejia/fMRIscrub/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/mandymejia/fMRIscrub/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/mandymejia/fMRIscrub/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/mandymejia/fMRIscrub/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/mandymejia/fMRIscrub/commit/5116b0778a69b1279d4ab39561c03d4905f6848b">5116b0778</a><a href="https://github.com/mandymejia/fMRIscrub/commit/5116b0778a69b1279d4ab39561c03d4905f6848b">&quot;&gt;&amp;lt;a href=&amp;quot;https://github.com/mandymejia/fMRIscrub/commit/&lt;/a&gt;&lt;a class=&quot;double-link&quot; href=&quot;https://github.com/mandymejia/fMRIscrub/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/mandymejia/fMRIscrub/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/mandymejia/fMRIscrub/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/mandymejia/fMRIscrub/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/mandymejia/fMRIscrub/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/mandymejia/fMRIscrub/commit/5116b0778a69b1279d4ab39561c03d4905f6848b&quot;&gt;5116b0778&lt;/a&gt;&lt;a href=&quot;https://github.com/mandymeji... (continued)

Pull Request Pull Request #31: 7.0

Run Details

1117 of 1117 new or added lines in 19 files covered. (100.0%)

375 of 1257 relevant lines covered (29.83%)

69966.36 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/R/CompCor.R

#' 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)
}

#' CompCor
#'
#' The aCompCor algorithm (Behzadi et. al., 2007) 10.1016/j.neuroimage.2007.04.042
#' 
#' First, the principal components (PCs) of each noise 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_clever_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 \times 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"} and \code{"noise"}
#'
#'  The entry \code{"data"} will be a \code{V x T} matrix where each row is a 
#'  data voxel (if it was originally an array, the voxels will be in spatial
#'  order) time series with each noise PC regressed from it. Otherwise, this 
#'  entry will be \code{NULL}.
#'
#'  The entry \code{"noise"} is a list of \code{T} \eqn{x} \code{noise_nPC} PC 
#'  scores, one for each \code{ROI_noise}.
#'
#' @importFrom robustbase rowMedians
#'
#' @export
CompCor <- function(
  X, ROI_data="infer", ROI_noise=NULL, 
  noise_nPC=5, noise_erosion=NULL, 
  center=TRUE, scale=TRUE, nuisance=NULL
  ){

  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(ROI_data)) {
    # Perform nuisance regression.
    design <- do.call(cbind, out2$noise_comps)
    if (!(is.null(nuisance) || isFALSE(nuisance) || identical(nuisance, 0))) {
      nuisance <- check_design_matrix(nuisance, T_)
      design <- cbind(design, nuisance)
    }
    design <- check_design_matrix(cbind(1, design))
    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)
  }

  list(
    data=out1$X, 
    noise=structure(list(PCs=out2$noise_comps, var=out2$noise_var, ROI_noise=out1$ROI_noise), class="clever_CompCor")
  )
}

1	#' CompCor: get noise components
2	#'
3	#' Get noise components for aCompCor.
4	#'
5	#' @param X_noise The noise ROIs data
6	#' @param center,scale Center & scale robustly
7	#' @param noise_nPC Number of PCs to obtain for each noise ROI
8	#'
9	#' @return A list with components X, X_noise, ROI_data, ROI_noise, noise_nPC,
10	#' noise_erosion, noise_comps, and noise_var.
11	#'
12	#' @keywords internal
13	CompCor.noise_comps <- function(X_noise, center, scale, noise_nPC){
14	TOL <- 1e-8	×
15
16	N <- length(X_noise)	×
17	noise_comps <- vector("list", N); names(noise_comps) <- names(X_noise)	×
18	noise_var <- vector("list", N); names(noise_var) <- names(X_noise)	×
19	noise_vartotal <- vector("list", N); names(noise_vartotal) <- names(X_noise)	×
20
21	if (length(noise_nPC) == 1) {	×
22	noise_nPC <- as.list(rep(noise_nPC, N))	×
23	} else {
24	noise_nPC <- as.list(noise_nPC)	×
25	}
26	names(noise_nPC) <- names(X_noise)	×
27
28	for (ii in 1:N) {	×
29	T_ <- nrow(X_noise[[ii]])	×
30	if (ncol(X_noise[[ii]])==0) { next }	×
31
32	# Transpose.
33	X_noise[[ii]] <- t(X_noise[[ii]])	×
34	# Center.
35	if (center) { X_noise[[ii]] <- X_noise[[ii]] - c(rowMedians(X_noise[[ii]], na.rm=TRUE)) }	×
36	# Compute MADs.
37	mad <- 1.4826 * rowMedians(abs(X_noise[[ii]]), na.rm=TRUE)	×
38	X_constant <- mad < TOL	×
39	if (any(X_constant)) {	×
40	if (all(X_constant)) {	×
41	stop(paste0("All data locations in noise ROI ", ii, " are zero-variance.\n"))	×
42	} else {
43	warning(paste0("Warning: removing ", sum(X_constant),	×
44	" constant data locations (out of ", length(X_constant),	×
45	") in noise ROI ", ii,	×
46	".\n"))	×
47	}
48	}
49	mad <- mad[!X_constant]; X_noise[[ii]] <- X_noise[[ii]][!X_constant,]	×
50	# Scale.
51	if (scale) { X_noise[[ii]] <- X_noise[[ii]]/c(mad) }	×
52	# Revert transpose.
53	X_noise[[ii]] <- t(X_noise[[ii]])	×
54	if (ncol(X_noise[[ii]])==0) { next }	×
55
56	# Compute the PC scores.
57	x <- svd(tcrossprod(X_noise[[ii]]))	×
58	noise_var[[ii]] <- x$d	×
59	noise_vartotal[[ii]] <- sum(noise_var[[ii]])	×
60	if (noise_nPC[[ii]] < 1) {	×
61	# Use enough PCs to explain the desired proportion of variance.
62	noise_nPC[[ii]] <- min(	×
63	length(x$d),	×
64	sum(cumsum(noise_var[[ii]]) < noise_vartotal[[ii]]*noise_nPC[[ii]]) + 1	×
65	)
66	}
67	noise_comps[[ii]] <- x$u[,seq(noise_nPC[[ii]]),drop=FALSE]	×
68	noise_var[[ii]] <- noise_var[[ii]][seq(noise_nPC[[ii]])]	×
69	}
70
71	list(noise_comps=noise_comps, noise_var=noise_var, noise_vartotal=noise_vartotal)	×
72	}
73
74	#' CompCor
75	#'
76	#' The aCompCor algorithm (Behzadi et. al., 2007) 10.1016/j.neuroimage.2007.04.042
77	#'
78	#' First, the principal components (PCs) of each noise ROI are calculated. For each ROI,
79	#' voxels are centered and scaled (can be disabled with the arguments
80	#' \code{center} and \code{scale}), and then the PCs are calculated via the
81	#' singular value decomposition.
82	#'
83	#' Next, aCompCor is performed to remove the shared variation between the noise ROI
84	#' PCs and each location in the data. This is accomplished by a nuisance regression
85	#' using a design matrix with the noise ROI PCs, any additional regressors specified
86	#' by \code{nuisance}, and an intercept term. (To detrend the data and perform aCompCor
87	#' in the same regression, \code{nuisance} can be set to DCT bases obtained with
88	#' the function \code{\link{dct_bases}}.)
89	#'
90	#' @inheritParams data_clever_CompCor_Params
91	#' @inheritParams noise_Params
92	#' @param center,scale Center the columns of the noise ROI data by their medians,
93	#' and scale by their MADs? Default: \code{TRUE} for both. Note that this argument
94	#' affects the noise ROI data and not the data that is being cleaned with aCompCor.
95	#' Centering and scaling of the data being cleaned can be done after this function call.
96	#' @param nuisance Nuisance signals to regress from each data column in addition
97	#' to the noise ROI PCs. Should be a \eqn{T \times N} numeric matrix where
98	#' \eqn{N} represents the number of nuisance signals. To not perform any nuisance
99	#' regression set this argument to \code{NULL}, \code{0}, or \code{FALSE}.
100	#' Default: \code{NULL}.
101	#' @return A list with entries \code{"data"} and \code{"noise"}
102	#'
103	#' The entry \code{"data"} will be a \code{V x T} matrix where each row is a
104	#' data voxel (if it was originally an array, the voxels will be in spatial
105	#' order) time series with each noise PC regressed from it. Otherwise, this
106	#' entry will be \code{NULL}.
107	#'
108	#' The entry \code{"noise"} is a list of \code{T} \eqn{x} \code{noise_nPC} PC
109	#' scores, one for each \code{ROI_noise}.
110	#'
111	#' @importFrom robustbase rowMedians
112	#'
113	#' @export
114	CompCor <- function(
115	X, ROI_data="infer", ROI_noise=NULL,
116	noise_nPC=5, noise_erosion=NULL,
117	center=TRUE, scale=TRUE, nuisance=NULL
118	){
119
120	out1 <- format_data(	×
121	X=X, ROI_data=ROI_data, ROI_noise=ROI_noise,	×
122	noise_nPC=noise_nPC, noise_erosion=noise_erosion	×
123	)
124
125	out2 <- CompCor.noise_comps(	×
126	X_noise=out1$X_noise, center=center, scale=scale, noise_nPC=out1$noise_nPC	×
127	)
128
129	T_ <- nrow(out1$X)	×
130
131	if (any(ROI_data)) {	×
132	# Perform nuisance regression.
133	design <- do.call(cbind, out2$noise_comps)	×
134	if (!(is.null(nuisance) \|\| isFALSE(nuisance) \|\| identical(nuisance, 0))) {	×
135	nuisance <- check_design_matrix(nuisance, T_)	×
136	design <- cbind(design, nuisance)	×
137	}
138	design <- check_design_matrix(cbind(1, design))	×
139	out1$X <- nuisance_regression(out1$X, design)	×
140
141	# # Normalize data.
142	# out1$X <- t(out1$X)
143	# if (center) { out1$X <- out1$X - c(rowMedians(out1$X, na.rm=TRUE)) }
144	# if (scale) {
145	# mad <- 1.4826 * rowMedians(abs(out1$X), na.rm=TRUE)
146	# mad_inv <- ifelse(mad < 1e-8, 0, 1/mad)
147	# out1$X <- out1$X * mad_inv
148	# }
149	# out1$X <- t(out1$X)
150	} else {
151	out1["X"] <- list(NULL)	×
152	}
153
154	list(	×
155	data=out1$X,	×
156	noise=structure(list(PCs=out2$noise_comps, var=out2$noise_var, ROI_noise=out1$ROI_noise), class="clever_CompCor")	×
157	)
158	}

mandymejia / fMRIscrub / 199

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous