Nothing
#' Multi-threshold permutation correction
#'
#' Applies the \emph{multi-threshold permutation correction (MTPC)} method to
#' perform inference in graph theory analyses of brain MRI data.
#'
#' This is a multi-step procedure: (steps 3-4 are the time-consuming steps)
#' \enumerate{
#' \item Apply thresholds \eqn{\tau} to the networks, and compute network
#' metrics for all networks and thresholds. (already done beforehand)
#' \item Compute test statistics \eqn{S_{obs}} for each threshold. (done by
#' \code{\link{brainGraph_GLM}})
#' \item Permute group assignments and compute test statistics for each
#' permutation and threshold. (done by \code{\link{brainGraph_GLM}})
#' \item Build a null distribution of the maximum statistic across thresholds
#' (and across brain regions) for each permutation. (done by
#' \code{\link{brainGraph_GLM}})
#' \item Determine the critical value, \eqn{S_{crit}} from the null
#' distribution of maximum statistics.
#' \item Identify clusters where \eqn{S_{obs} > S_{crit}} and compute the AUC
#' for these clusters (denoted \eqn{A_{MTPC}}).
#' \item Compute a critical AUC (\eqn{A_{crit}}) from the mean of the
#' supra-critical AUC's for the permuted tests.
#' \item Reject \eqn{H_0} if \eqn{A_{MTPC} > A_{crit}}.
#' }
#'
#' @param g.list A list of \code{brainGraphList} objects for all thresholds
#' @param thresholds Numeric vector of the thresholds applied to the raw
#' connectivity matrices.
#' @param clust.size Integer indicating the size of \dQuote{clusters} (i.e.,
#' consecutive thresholds for which the observed statistic exceeds the null)
#' (default: \code{3L})
#' @param res.glm A list of \code{bg_GLM} objects, as output by a previous run
#' of \code{mtpc}. Useful if you want to change the cluster size without
#' re-running all of the GLM's and permutations (default: \code{NULL})
#' @param ... Other arguments passed to \code{\link{brainGraph_GLM}} and/or
#' \code{\link{brainGraph_GLM_design}}
#' @inheritParams GLM
#' @export
#' @importFrom permute shuffleSet
#'
#' @return An object of class \code{mtpc} with some input arguments plus the
#' following elements:
#' \item{X,qr,cov.unscaled}{Design matrix, QR decomposition, and unscaled
#' covariance matrix, if the design is the same across thresholds}
#' \item{contrasts}{The contrast matrix or list of matrices}
#' \item{con.name}{Contrast names}
#' \item{removed.subs}{Named integer vector of subjects with incomplete data}
#' \item{atlas}{The atlas of the input graphs}
#' \item{rank,df.residual}{The model rank and residual degrees of freedom}
#' \item{res.glm}{List with length equal to the number of thresholds; each
#' list element is the output from \code{\link{brainGraph_GLM}}}
#' \item{DT}{A \code{data.table} for all thresholds, combined from the outputs
#' of \code{\link{brainGraph_GLM}}}
#' \item{stats}{A data.table containing \code{S.mtpc} (the max. observed
#' statistic), \code{tau.mtpc} (the threshold of the max. observed
#' statistic), \code{S.crit} (the critical statistic value), and
#' \code{A.crit} (the critical AUC)}
#' \item{null.dist}{Numeric array with \code{N} columns and number of rows
#' equal to the number of thresholds. The 3rd dimension is for each
#' contrast. Each element of the array is the maximum statistic
#' for that permutation, threshold, and contrast combination.}
#' \item{perm.order}{Numeric matrix; the permutation set applied for all
#' thresholds (each row is a separate permutation)}
#'
#' @family Group analysis functions
#' @family GLM functions
#' @author Christopher G. Watson, \email{cgwatson@@bu.edu}
#' @references Drakesmith, M. and Caeyenberghs, K. and Dutt, A. and Lewis, G. and
#' David, A.S. and Jones, D.K. (2015) Overcoming the effects of false
#' positives and threshold bias in graph theoretical analyses of neuroimaging
#' data. \emph{NeuroImage}, \bold{118}, 313--333.
#' \doi{10.1016/j.neuroimage.2015.05.011}
#' @examples
#' \dontrun{
#' diffs.mtpc <- mtpc(g.list=g.norm, thresholds=thresholds, N=N,
#' covars=covars.dti, measure='E.nodal.wt', coding='effects',
#' contrasts=c(0, 0, 0, 0, -2), alt='greater',
#' binarize=c('Sex', 'Scanner'), con.name='Group 1 > Group 2')
#' sig.regions <- diffs.mtpc$DT[A.mtpc > A.crit]
#' }
mtpc <- function(g.list, thresholds, covars, measure, contrasts, con.type=c('t', 'f'),
outcome=NULL, con.name=NULL, level=c('vertex', 'graph'),
clust.size=3L, perm.method=c('freedmanLane', 'terBraak', 'smith',
'draperStoneman', 'manly', 'stillWhite'),
part.method=c('beckmann', 'guttman', 'ridgway'), N=500L, perms=NULL,
alpha=0.05, res.glm=NULL, long=TRUE, ...) {
S.mtpc <- tau.mtpc <- A.crit <- A.mtpc <- Contrast <- S.crit <- region <- stat <-
threshold <- values <- null.out <- NULL
# Check if components are 'brainGraphList' objects
matches <- vapply(g.list, is.brainGraphList, logical(1L))
if (any(!matches)) stop("Input must be a list of 'brainGraphList' objects.")
stopifnot(length(g.list) == length(thresholds))
if (is.null(res.glm)) {
# If 'outcome != measure' make sure the same perm order is used throughout
if (is.null(perms)) perms <- shuffleSet(n=length(g.list[[1L]]$graphs), nset=N)
if (!is.null(outcome) && (outcome != measure)) {
n <- sum(complete.cases(covars))
if (dim(perms)[1L] != n) perms <- shuffleSet(n=n, nset=N)
}
res.glm <- lapply(g.list, function(z)
brainGraph_GLM(z, covars, measure, contrasts, con.type, outcome, con.name=con.name, N=N,
level=level, permute=TRUE, perm.method=perm.method, part.method=part.method,
perms=perms, alpha=alpha, long=TRUE, ...))
}
alt <- res.glm[[1L]]$alt
for (i in seq_along(thresholds)) res.glm[[i]]$DT[, threshold := thresholds[i]]
mtpc.all <- rbindlist(lapply(res.glm, function(x) x$DT))
setkey(mtpc.all, Contrast, region)
mtpc.all[, c('S.crit', 'A.mtpc', 'A.crit') := 0]
conNames <- res.glm[[1L]]$con.name
myMax <- switch(alt, two.sided=colMaxAbs, less=colMin, greater=colMax)
mySort <- sortfun(alt)
index <- floor((1 - alpha) * N) + 1L
# Dimensions will be "thresholds X Nperms X contrasts"
null.max.all <- aperm(sapply(res.glm, function(x) x$perm$null.max, simplify='array'),
c(3L, 1L, 2L))
dimnames(null.max.all) <- list(NULL, NULL, conNames)
null.max.perms <- apply(null.max.all, 3L, myMax)
Scrit <- structure(apply(null.max.perms, 2L, mySort, index), names=conNames)
for (x in conNames) mtpc.all[Contrast == x, S.crit := Scrit[x]]
rlefun <- rle_comp(alt)
rle_obs <- mtpc.all[, rlefun(stat, S.crit), by=key(mtpc.all)]
crit.regions <- rle_obs[(values)][lengths >= clust.size, list(region=unique(region)), by=Contrast]
setkeyv(crit.regions, key(mtpc.all))
mtpc.all[crit.regions,
A.mtpc := auc_rle(stat, S.crit, thresholds, alt, clust.size),
by=key(mtpc.all)]
null.crit <- setNames(vector('list', length(conNames)), conNames)
Acrit <- setNames(rep.int(0, length(conNames)), conNames)
for (x in conNames) {
null.crit[[x]] <- which(apply(null.max.all[, , x], 2L, function(y)
any(with(rlefun(y, Scrit[x]), values == TRUE & lengths >= clust.size))))
if (length(null.crit) > 0L) {
null.crits <- abind::adrop(null.max.all[, null.crit[[x]], x, drop=FALSE], drop=3L)
Acrit[x] <- mean(apply(null.crits, 2L, auc_rle, Scrit[x], thresholds, alt, clust.size))
mtpc.all[Contrast == x, A.crit := Acrit[x]]
}
}
maxobsfun <- switch(alt,
two.sided=function(x) max(is.finite(abs(x)) * abs(x), na.rm=TRUE),
less=function(x) min(is.finite(x) * x, na.rm=TRUE),
greater=function(x) max(is.finite(x) * x, na.rm=TRUE))
mtpc.stats <- mtpc.all[, list(S.mtpc=unique(maxobsfun(stat))), by=Contrast]
mtpc.stats[, tau.mtpc := mtpc.all[mtpc.stats, .SD[S.mtpc == stat, threshold]]]
mtpc.stats[, S.crit := Scrit]
mtpc.stats[, A.crit := Acrit]
setcolorder(mtpc.stats, c('Contrast', 'tau.mtpc', 'S.mtpc', 'S.crit', 'A.crit'))
glm.attr <- res.glm[[1L]]
#TODO: may have to remove also 'var.covar'
#TODO: is it possible for 'removed.subs' and 'df.residual' to be different?
glm.attr[c('y', 'DT.Xy', 'DT', 'permute', 'runX', 'runY', 'coefficients', 'residuals',
'sigma', 'fitted.values', 'se', 'perm', 'perm.order')] <- NULL
if (glm.attr$outcome != glm.attr$measure) glm.attr[c('X', 'qr', 'cov.unscaled')] <- NULL
for (i in seq_along(thresholds)) {
res.glm[[i]]$perm[c('null.dist', 'null.max')] <- NULL
res.glm[[i]]$perm.order <- NULL
}
if (isTRUE(long)) null.out <- null.max.all
out <- c(glm.attr, list(res.glm=res.glm, clust.size=clust.size, DT=mtpc.all, stats=mtpc.stats,
null.dist=null.out, perm.order=perms))
class(out) <- c('mtpc', class(out))
return(out)
}
#' Run length encoding of the comparison of 2 vectors
#'
#' Returns a function that will compare 2 vectors and return the run length
#' encoding (RLE) of the logical vector from that comparison.
#' @noRd
rle_comp <- function(alt) {
switch(alt,
two.sided=function(x, y) rle(abs(x) > abs(y)),
less=function(x, y) rle(x < y),
greater=function(x, y) rle(x > y))
}
#' Calculate the AUC for the observed T- or F-statistic
#'
#' @noRd
auc_rle <- function(t.stat, S.crit, thresholds, alt, clust.size) {
inds <- get_rle_inds(clust.size, alt, t.stat, S.crit, thresholds)
x <- Map(function(a, b) thresholds[a:b], inds$starts, inds$ends)
y <- Map(function(a, b) t.stat[a:b], inds$starts, inds$ends)
auc <- mapply(auc_diff, x, y)
return(sum(abs(auc) * is.finite(auc), na.rm=TRUE))
}
#' Get the starting and ending points for significant runs
#'
#' Given a vector of observed statistics and a critical value, calculate the
#' \dQuote{runs} (using \code{\link{rle}}) based on comparing these values. Then
#' return the starting and ending points of the runs.
#' @noRd
get_rle_inds <- function(clust.size, alt, t.stat, S.crit, thresholds) {
compfun <- rle_comp(alt)
runs <- compfun(t.stat, S.crit)
sigruns <- which(runs$values == TRUE & runs$lengths >= clust.size)
ends <- cumsum(runs$lengths)
starts <- c(1L, head(ends, -1L) + 1L)
# ends <- runs.len.cumsum[myruns]
# newindex <- ifelse(sigruns > 1L, sigruns - 1L, 0L)
# if (0L %in% newindex) starts <- c(1L, starts)
return(list(starts=starts[sigruns], ends=ends[sigruns]))
}
################################################################################
# S3 METHODS
################################################################################
#' Print a summary of MTPC results
#'
#' @param object,x A \code{mtpc} object
#' @inheritParams summary.bg_GLM
#' @export
#' @rdname mtpc
summary.mtpc <- function(object, contrast=NULL, digits=max(3L, getOption('digits') - 2L), print.head=TRUE, ...) {
Contrast <- A.mtpc <- A.crit <- stat <- region <- S.mtpc <- NULL
object$printCon <- contrast
object$digits <- digits
object$print.head <- print.head
# Summary table
whichmaxfun <- switch(object$alt, two.sided=function(x) which.max(abs(x)), less=which.min, greater=which.max)
DT.sum <- object$DT[A.mtpc > A.crit, .SD[whichmaxfun(is.finite(stat) * stat)], by=list(Contrast, region)]
DT.sum <- DT.sum[, c('region', 'stat', 'Contrast', 'threshold', 'S.crit', 'A.mtpc', 'A.crit'), with=FALSE]
setcolorder(DT.sum, c('Contrast', 'region', 'threshold', 'stat', 'S.crit', 'A.mtpc', 'A.crit'))
setnames(DT.sum, c('region', 'stat', 'threshold'), c('Region', 'S.mtpc', 'tau.mtpc'))
setorder(DT.sum, Contrast, -S.mtpc)
object$DT.sum <- DT.sum
class(object) <- c('summary.mtpc', class(object))
return(object)
}
#' @aliases summary.mtpc
#' @method print summary.mtpc
#' @export
print.summary.mtpc <- function(x, ...) {
print_title_summary('MTPC results (', x$level, '-level)')
print_model_summary(x)
print_contrast_type_summary(x)
print_subs_summary(x)
cat('# of thresholds: ', length(x$res.glm), '\n')
cat('Cluster size (across thresholds): ', x$clust.size, '\n')
print_permutation_summary(x)
message('\n', 'Statistics', '\n', rep.int('-', getOption('width') / 4))
clp <- 100 * (1 - x$alpha)
cat('tau.mtpc: threshold of the maximum observed statistic\n')
cat('S.mtpc: maximum observed statistic\n')
cat('S.crit: the critical', paste0('(', clp, 'th percentile)'), 'value of the null max. statistic\n')
cat('A.crit: critical AUC from the supra-critical null AUCs\n\n')
print(x$stats)
cat('\n')
# Print results for each contrast
print_contrast_stats_summary(x)
invisible(x)
}
#' Plot statistics from an MTPC analysis
#'
#' Plot the statistics from an MTPC analysis, along with the maximum permuted
#' statistics. The output is similar to Figure 11 in Drakesmith et al. (2015).
#'
#' @param only.sig.regions Logical indicating whether to plot only significant
#' regions (default: \code{TRUE})
#' @param show.null Logical indicating whether to plot points of the maximum
#' null statistics (per permutation)
#' @param caption.stats Logical indicating whether to print the MTPC statistics
#' in the caption of the plot. Default: \code{FALSE}
#' @inheritParams plot.bg_GLM
#' @export
#' @rdname mtpc
#'
#' @return The \code{plot} method returns a \code{trellis} object or a list of
#' \code{ggplot} objects
#' @examples
#' \dontrun{
#' mtpcPlots <- plot(mtpc.diffs)
#'
#' ## Arrange plots into 3x3 grids
#' ml <- marrangeGrob(mtpcPlots, nrow=3, ncol=3)
#' ggsave('mtpc.pdf', ml)
#' }
plot.mtpc <- function(x, contrast=1L, region=NULL, only.sig.regions=TRUE,
show.null=TRUE, caption.stats=FALSE, ...) {
Contrast <- stat_ribbon <- stat <- threshold <- S.crit <- A.mtpc <- A.crit <- panel.num <-
xstart <- starts <- xend <- NULL
stopifnot(inherits(x, 'mtpc'))
conNames <- x$con.name
mycontrast <- conNames[contrast]
DT <- x$DT[Contrast == mycontrast, !'p.fdr']
DT[, stat_ribbon := stat] # For filling in supra-threshold areas
thresholds <- DT[, unique(threshold)]
DT$nullthresh <- x$stats[Contrast == mycontrast, unique(S.crit)]
whichmaxfun <- switch(x$alt, two.sided=function(y) which.max(abs(y)), less=which.min, greater=which.max)
myMax <- switch(x$alt, two.sided=colMaxAbs, less=colMin, greater=colMax)
thr <- apply(x$null.dist[, , mycontrast], 2L, whichmaxfun)
thr.y <- myMax(x$null.dist[, , mycontrast])
nullcoords <- data.table(threshold=thresholds[thr], y=thr.y)
# 'ggplot2'; Local function to plot for a single region
plot_single <- function(x, DT, nullcoords, show.null) {
stat <- S.crit <- threshold <- stat_ribbon <- nullthresh <- y <- A.mtpc <- A.crit <- NULL
myMax <- switch(x$alt, two.sided=function(y) max(abs(y)), less=min, greater=max)
inds <- DT[, get_rle_inds(x$clust.size, x$alt, stat, S.crit, threshold)]
n <- dim(inds)[1L]
if (n == 0L) {
DT[, stat_ribbon := NA]
} else {
if (n == 1L) {
all.inds <- c(apply(inds, 1L, function(z) z[1L]:z[2L]))
} else if (n > 1L) {
all.inds <- c(unlist(apply(inds, 1L, function(z) z[1L]:z[2L])))
} else {
all.inds <- 1L:dim(DT)[1L]
}
DT[-all.inds, stat_ribbon := NA]
}
lineplot <- ggplot2::ggplot(data=DT, mapping=ggplot2::aes(x=threshold)) +
ggplot2::geom_line(ggplot2::aes(y=stat), col='red4', size=1.25, na.rm=TRUE) +
ggplot2::geom_hline(ggplot2::aes(yintercept=nullthresh), lty=2)
if (n > 0L) {
if (x$alt == 'less') {
lineplot <- lineplot +
ggplot2::geom_ribbon(ggplot2::aes(ymax=stat_ribbon, ymin=nullthresh), fill='red4', alpha=0.3)
} else {
lineplot <- lineplot +
ggplot2::geom_ribbon(ggplot2::aes(ymin=stat_ribbon, ymax=nullthresh), fill='red4', alpha=0.3)
}
lineplot <- lineplot + ggplot2::geom_point(ggplot2::aes(y=stat_ribbon), col='red4', size=2, na.rm=TRUE)
}
if (isTRUE(all(diff(thresholds) < 0))) lineplot <- lineplot + ggplot2::scale_x_reverse()
p.region <- if (x$level == 'graph') 'graph-level' else DT[, as.character(unique(region))]
p.title <- paste0('Region: ', p.region)
p.subtitle <- paste0('Outcome: ', x$outcome)
if (isTRUE(show.null)) {
lineplot <- lineplot + ggplot2::geom_point(data=nullcoords, ggplot2::aes(y=y), col='darkgreen', alpha=0.4, na.rm=TRUE)
}
if (isTRUE(caption.stats)) {
Smtpc <- bquote('S'['mtpc']*' = '~ .(DT[, format(myMax(stat))]))
Scrit <- bquote('S'['crit']*' = '~ .(DT[, format(unique(S.crit))]))
Amtpc <- bquote('A'['mtpc']*' = '~ .(DT[, format(max(A.mtpc))]))
Acrit <- bquote('A'['crit']*' = '~ .(DT[, format(unique(A.crit))]))
statslabel <- bquote(atop(.(Scrit)~ ';'~ .(paste('\t'))~ .(Acrit),
.(Smtpc)~ ';'~ .(paste('\t'))~ .(Amtpc)))
if (DT[, unique(A.mtpc) > unique(A.crit)]) statslabel <- bquote(.(statslabel)~ .(paste0('\t(p < ', x$alpha, ')')))
lineplot <- lineplot +
ggplot2::labs(caption=statslabel) +
ggplot2::theme(plot.caption=ggplot2::element_text(hjust=0))
}
lineplot <- lineplot +
ggplot2::labs(title=p.title, subtitle=p.subtitle, x='Threshold', y=paste0(toupper(x$con.type), '-statistic')) +
ggplot2::theme(plot.title=ggplot2::element_text(hjust=0.5, face='bold'),
plot.subtitle=ggplot2::element_text(hjust=0.5))
return(lineplot)
}
# 'base' plotting
#-------------------------------------------------------
if (!requireNamespace('ggplot2', quietly=TRUE)) {
if (x$level == 'vertex') {
if (is.null(region)) {
if (isFALSE(only.sig.regions)) {
myregions <- DT[, levels(region)]
} else {
myregions <- droplevels(DT[A.mtpc > A.crit])[, levels(region)]
}
} else {
myregions <- region
}
}
DT <- droplevels(DT[region %in% myregions])
DT[, panel.num := as.numeric(region)]
inds <- droplevels(DT[, get_rle_inds(x$clust.size, x$alt, stat, S.crit, threshold), by=region])
inds[, panel.num := as.numeric(region)]
inds[, xstart := thresholds[starts]]
inds[, xend := thresholds[ends]]
if (isTRUE(show.null)) {
panelfun <- function(x, y, ...) {
panel.num <- starts <- nullthresh <- NULL
panel.points(-nullcoords$threshold, nullcoords$y, pch=19, fill='darkgreen')
for (i in seq_len(inds[panel.num == panel.number(), .N])) {
xcoords <- -thresholds[inds[panel.num == panel.number()][i, seq(starts, ends)]]
ycoords <- DT[panel.num == panel.number() & threshold %in% -xcoords, c(stat, nullthresh)]
panel.polygon(x=c(xcoords, rev(xcoords)), y=ycoords, col='lightgrey', border=NULL)
}
panel.abline(h=DT[, nullthresh], lty=2)
panel.xyplot(x, y, ..., col='red', lwd=2)
}
} else {
panelfun <- function(x, y, ...) {
panel.num <- starts <- nullthresh <- NULL
for (i in seq_len(inds[panel.num == panel.number(), .N])) {
xcoords <- thresholds[inds[panel.num == panel.number()][i, seq(starts, ends)]]
ycoords <- DT[panel.num == panel.number() & threshold %in% xcoords, c(stat, nullthresh)]
panel.polygon(x=c(xcoords, rev(xcoords)), y=ycoords, col='lightgrey', border=NULL)
}
panel.abline(h=DT[, nullthresh], lty=2)
panel.xyplot(x, y, ..., col='red', lwd=2)
}
}
lineplot <- xyplot(stat ~ -threshold | region, data=DT, type=c('l', 'p'), pch=19,
xlab='Threshold', ylab=paste0(toupper(x$con.type), '-statistic'),
main=paste0('Outcome: ', x$outcome),
scales=list(y=list(relation='free')),
panel=panelfun)
return(lineplot)
# 'ggplot2' plotting
#-------------------------------------------------------
} else {
if (x$level == 'graph') {
lineplot <- plot_single(x, DT, nullcoords, show.null)
return(lineplot)
} else if (x$level == 'vertex') {
if (is.null(region)) {
if (isFALSE(only.sig.regions)) {
region <- DT[, levels(region)]
} else {
region <- droplevels(DT[A.mtpc > A.crit])[, levels(region)]
}
}
lineplots <- setNames(vector('list', length(region)), region)
for (z in region) {
lineplots[[z]] <- plot_single(x, DT[region == z], nullcoords, show.null)
}
return(lineplots)
}
}
}
#' @export
#' @rdname mtpc
nobs.mtpc <- function(object, ...) nobs(object$res.glm[[1L]])
#' @export
#' @rdname mtpc
terms.mtpc <- function(x, ...) terms(x$res.glm[[1]])
#' @export
#' @rdname mtpc
formula.mtpc <- formula.bg_GLM
#' @export
#' @rdname mtpc
labels.mtpc <- labels.bg_GLM
#' @method case.names mtpc
#' @export
#' @rdname mtpc
case.names.mtpc <- function(object, ...) case.names(object$res.glm[[1]])
#' @method variable.names mtpc
#' @export
#' @rdname mtpc
variable.names.mtpc <- function(object, ...) variable.names(object$res.glm[[1]])
#' @method df.residual mtpc
#' @export
#' @rdname mtpc
df.residual.mtpc <- df.residual.bg_GLM
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.