Nothing
#' Permutation test for group difference of graph measures
#'
#' \code{brainGraph_permute} draws permutations from linear model residuals to
#' determine the significance of between-group differences of a global or
#' vertex-wise graph measure. It is intended for structural covariance networks
#' (in which there is only one graph per group), but can be extended to other
#' types of data.
#'
#' If you would like to calculate differences in the area-under-the-curve (AUC)
#' across densities, then specify \code{auc=TRUE}.
#'
#' There are three possible \dQuote{levels}:
#' \enumerate{
#' \item \emph{graph} Calculate modularity (Louvain algorithm), clustering
#' coefficient, characteristic path length, degree assortativity, and global
#' efficiency.
#' \item \emph{vertex} Choose one of: centrality metrics (betweenness,
#' closeness, communicability, eigenvector, leverage, pagerank, subgraph);
#' k-core; degree; eccentricity; nodal or local efficiency; k-nearest neighbor
#' degree; shortest path length; transitivity; or vulnerability.
#' \item \emph{other} Supply your own function. This is useful if you want to
#' calculate something that I haven't hard-coded. It must take as its own
#' arguments: \code{g} (a list of lists of \code{igraph} graph objects); and
#' \code{densities} (numeric vector).
#' }
#'
#' @param densities Numeric vector of graph densities
#' @param resids An object of class \code{brainGraph_resids} (the output from
#' \code{\link{get.resid}})
#' @param N Integer; the number of permutations (default: 5e3)
#' @param perms Numeric matrix of permutations, if you would like to provide
#' your own (default: \code{NULL})
#' @param auc Logical indicating whether or not to calculate differences in the
#' area-under-the-curve of metrics (default: \code{FALSE})
#' @param level A character string for the attribute \dQuote{level} to calculate
#' differences (default: \code{graph})
#' @param measure A character string specifying the vertex-level metric to
#' calculate, only used if \code{level='vertex'} (default: \code{btwn.cent}).
#' For the \code{summary} method, this is to focus on a single
#' \emph{graph-level} measure (since multiple are calculated at once).
#' @param .function A custom function you can pass if \code{level='other'}
#' @export
#' @importFrom permute shuffleSet
#' @importFrom foreach getDoParRegistered
#' @importFrom doParallel registerDoParallel
#'
#' @return An object of class \code{brainGraph_permute} with input arguments in
#' addition to:
#' \item{DT}{A data table with permutation statistics}
#' \item{obs.diff}{A data table of the observed group differences}
#' \item{Group}{Group names}
#'
#' @family Group analysis functions
#' @family Structural covariance network functions
#'
#' @author Christopher G. Watson, \email{cgwatson@@bu.edu}
#' @examples
#' \dontrun{
#' myResids <- get.resid(lhrh, covars)
#' myPerms <- shuffleSet(n=nrow(myResids$resids.all), nset=1e3)
#' out <- brainGraph_permute(densities, m, perms=myPerms)
#' out <- brainGraph_permute(densities, m, perms=myPerms, level='vertex')
#' out <- brainGraph_permute(densities, m, perms=myPerms,
#' level='other', .function=myFun)
#' }
brainGraph_permute <- function(densities, resids, N=5e3, perms=NULL, auc=FALSE,
level=c('graph', 'vertex', 'other'),
measure=c('btwn.cent', 'coreness', 'degree', 'eccentricity',
'clo.cent', 'communicability', 'ev.cent', 'lev.cent',
'pagerank', 'subg.cent', 'E.local', 'E.nodal',
'knn', 'Lp', 'transitivity', 'vulnerability'),
.function=NULL) {
stopifnot(inherits(resids, 'brainGraph_resids'))
gID <- getOption('bg.group')
measure <- match.arg(measure)
level <- match.arg(level)
if (level == 'other') {
if (!is.function(.function)) stop('".function" must be a function!')
measure <- 'other'
} else if (level == 'graph') {
measure <- NULL
}
if (is.null(perms)) perms <- shuffleSet(n=nobs(resids), nset=N)
dims <- dim(perms)
N <- dims[1L]
perms <- rbind(perms, seq_len(dims[2L])) # last row is observed metrics
grps <- as.numeric(resids$resids.all[, get(gID)])
if (!getDoParRegistered()) {
cl <- makeCluster(getOption('bg.ncpus'))
registerDoParallel(cl)
}
# Should return a "data.table" (graph-level) or numeric matrix (vertex-level)
# When "auc=FALSE", should have a "densities" column
res.perm <- switch(level,
vertex=permute_vertex_foreach(perms, densities, resids, grps, auc, measure),
other=permute_other_foreach(perms, densities, resids, grps, .function),
graph=permute_graph_foreach(perms, densities, resids, grps, auc))
if (length(densities) == 1L) res.perm <- cbind(densities=densities, res.perm)
if (level == 'vertex') {
res.perm <- as.data.table(res.perm)
start <- if (isTRUE(auc)) 1L else 2L
setnames(res.perm, seq.int(start, dim(res.perm)[2L]), region.names(resids))
}
obs.ind <- if (isTRUE(auc)) N + 1L else (N + 1L) * seq_along(densities)
obs.diff <- res.perm[obs.ind]
res.perm <- res.perm[-obs.ind]
out <- list(atlas=resids$atlas, auc=auc, N=N, level=level, measure=measure, densities=densities,
resids=resids, DT=res.perm, obs.diff=obs.diff, Group=resids$Group)
class(out) <- c('brainGraph_permute', class(out))
return(out)
}
#==============================================================================
# Helper functions
#==============================================================================
make_graphs_perm <- function(densities, resids, inds, grps) {
corrs <- lapply(unique(grps), function(x)
corr.matrix(resids[which(grps[inds] == x)],
densities=densities, rand=TRUE))
sapply(corrs, lapply, function(x)
apply(x$r.thresh, 3L, graph_from_adjacency_matrix, mode='undirected', diag=FALSE))
}
# Graph level
graph_attr_perm <- function(g, densities) {
mod <- sapply(g, sapply, function(x) modularity(cluster_louvain(x)))
Cp <- sapply(g, sapply, function(x) transitivity(x, type='localaverage'))
Lp <- sapply(g, sapply, mean_distance)
assort <- sapply(g, sapply, assortativity_degree)
E.global <- sapply(g, sapply, efficiency, 'global')
list(mod=mod, Cp=Cp, Lp=Lp, assort=assort, E.global=E.global)
}
permute_graph_foreach <- function(perms, densities, resids, grps, auc) {
i <- NULL
N <- dim(perms)[1L]
if (isTRUE(auc)) {
res.perm <- foreach(i=seq_len(N), .combine='rbind') %dopar% {
g <- make_graphs_perm(densities, resids, perms[i, ], grps)
meas.list <- graph_attr_perm(g, densities)
t(vapply(meas.list, function(y) auc_diff_perm(densities, y), numeric(1L)))
}
res.perm <- as.data.table(res.perm)
} else {
# Returns a list of (Ndensities * Nperms) X 2 matrices for all metrics
res.perm <- foreach(i=seq_len(N), .combine=function(a, b) Map(rbind, a, b)) %dopar% {
g <- make_graphs_perm(densities, resids, perms[i, ], grps)
meas.list <- graph_attr_perm(g, densities)
}
res.perm <- data.table(densities=rep.int(densities, N),
sapply(res.perm, function(x) as.numeric(-diff(t(x)))))
setkey(res.perm, densities)
}
return(res.perm)
}
# Vertex-level
vertex_attr_funs <- function(measure) { # Move the "switch" outside the "foreach" loop
switch(measure,
coreness=coreness,
degree=degree,
eccentricity=eccentricity,
clo.cent=function(x) centr_clo(x)$res,
communicability=centr_betw_comm,
ev.cent=function(x) centr_eigen(x)$vector,
lev.cent=function(x) centr_lev(x)$res,
pagerank=function(x) page_rank(x, weights=NA)$vector,
subg.cent=subgraph_centrality,
E.local=function(x) efficiency(x, type='local', weights=NA, use.parallel=FALSE),
E.nodal=function(x) efficiency(x, type='nodal', weights=NA),
knn=function(x) knn(x)$knn,
Lp=function(x) mean_distance_wt(x, level='vertex', weights=NA),
transitivity=function(x) transitivity(x, type='local', isolates='zero'),
vulnerability=function(x) vulnerability(x, use.parallel=FALSE),
function(x) centr_betw(x)$res)
}
permute_vertex_foreach <- function(perms, densities, resids, grps, auc, measure) {
i <- NULL
if (isTRUE(auc)) {
diffFun <- function(densities, meas.list) {
sapply(seq_len(dim(meas.list[[1L]])[2L]), function(x)
auc_diff(densities, cbind(meas.list[[1L]][, x], meas.list[[2L]][, x])))
}
} else {
diffFun <- function(densities, meas.list) cbind(densities, meas.list[[1L]] - meas.list[[2L]])
}
fun <- vertex_attr_funs(measure)
res.perm <- foreach(i=seq_len(dim(perms)[1L]), .combine='rbind') %dopar% {
g <- make_graphs_perm(densities, resids, perms[i, ], grps)
meas.list <- lapply(g, function(x) t(sapply(x, fun)))
diffFun(densities, meas.list)
}
}
# Other-level
permute_other_foreach <- function(perms, densities, resids, grps, .function) {
i <- NULL
res.perm <- foreach(i=seq_len(dim(perms)[1L]), .combine='rbind') %dopar% {
g <- make_graphs_perm(densities, resids, perms[i, ], grps)
.function(g, densities)
}
}
#==============================================================================
# Methods
#==============================================================================
#' Print a summary from a permutation analysis
#'
#' @param object,x A \code{brainGraph_permute} object (output by
#' \code{\link{brainGraph_permute}}).
#' @param p.sig Character string specifying which p-value to use for displaying
#' significant results (default: \code{p})
#' @param ... Unused
#' @inheritParams GLM
#' @export
#' @rdname brainGraph_permute
summary.brainGraph_permute <- function(object, measure=object$measure,
alternative=c('two.sided', 'less', 'greater'),
alpha=0.05, p.sig=c('p', 'p.fdr'), ...) {
perm.diff <- p <- p.fdr <- region <- obs.diff <- ..measure <- NULL
gID <- getOption('bg.group')
if (object$level == 'other') { # Hack to figure out which level it is when level="other"
object$level <- if (dim(object$DT)[2L] > 6L) 'vertex' else 'graph'
}
if (object$level == 'graph' && is.null(measure)) measure <- 'mod'
group_str <- paste0(measure, '.', object$Group)
permDT <- copy(object$DT)
g <- with(object, make_graphs_perm(densities, resids, seq_len(nobs(resids)),
resids$resids.all[, as.numeric(get(gID))]))
densities <- object$densities
# VERTEX-LEVEL
#-------------------------------------
if (object$level == 'vertex') {
obsDT <- copy(object$obs.diff)
fun <- vertex_attr_funs(measure)
obs <- lapply(g, function(x) t(sapply(x, fun)))
if (isTRUE(object$auc)) {
obs <- lapply(obs, apply, 2L, function(y) -auc_diff(densities, y))
permDT[, densities := 1]
obsDT[, densities := 1]
densities <- 1
}
sum.dt <- data.table(densities=densities,
region=rep(V(g[[1L]][[1L]])$name, each=length(densities)),
g1=c(obs[[1L]]), g2=c(obs[[2L]]),
key=c('densities', 'region'))
setnames(sum.dt, c('g1', 'g2'), group_str)
obs.m <- melt(obsDT, id.vars='densities', variable.name='region', value.name='obs.diff')
sum.dt <- merge(sum.dt, obs.m)
permDT <- melt(permDT, id.vars='densities', variable.name='region', value.name=measure)
setkeyv(permDT, key(sum.dt))
permdiff <- permDT[, list(perm.diff=mean(get(measure))), by=key(sum.dt)]
sum.dt <- merge(sum.dt, permdiff, by=key(sum.dt))
# GRAPH-LEVEL
#-------------------------------------
} else if (object$level == 'graph') {
stopifnot(hasName(permDT, measure))
permDT[, region := 'graph']
obs <- graph_attr_perm(g, densities)[[measure]]
if (isTRUE(object$auc)) {
obs <- t(apply(obs, 2L, function(y) -auc_diff(densities, y)))
permDT[, densities := 1]
densities <- 1
}
sum.dt <- data.table(densities=densities, region='graph', obs)
setnames(sum.dt, c('V1', 'V2'), group_str)
sum.dt[, obs.diff := object$obs.diff[[measure]]]
sum.dt[, perm.diff := permDT[, mean(get(measure)), by=densities]$V1]
}
result.dt <- merge(permDT[, c('densities', 'region', ..measure)],
sum.dt[, c('densities', 'region', 'obs.diff')],
by=c('densities', 'region'))
alt <- match.arg(alternative)
compfun <- switch(alt,
two.sided=function(perm, obs) sum(abs(perm) >= abs(unique(obs))),
less=function(perm, obs) sum(perm <= unique(obs)),
greater=function(perm, obs) sum(perm >= unique(obs)))
result.dt[, p := (compfun(get(measure), obs.diff) + 1L) / (.N + 1L), by=key(result.dt)]
CI <- switch(alt, two.sided=c(alpha / 2, 1 - (alpha / 2)), less=c(alpha, 1), greater=c(1 / object$N, 1 - alpha))
result.dt[, c('ci.low', 'ci.high') := as.list(sort(get(measure))[ceiling(.N * CI)]), by=key(result.dt)]
result.dt <- result.dt[, .SD[1L], by=key(result.dt)]
sum.dt <- merge(sum.dt, result.dt[, !c(measure, 'obs.diff'), with=FALSE], by=key(result.dt))
setcolorder(sum.dt, c('densities', 'region', group_str, 'obs.diff', 'perm.diff', 'ci.low', 'ci.high', 'p'))
if (isFALSE(object$auc)) {
groupby <- if (object$level == 'graph') NULL else 'densities'
sum.dt[, p.fdr := p.adjust(p, 'fdr'), by=groupby]
}
meas.full <- print_full_metric(measure)
p.sig <- match.arg(p.sig)
perm.sum <- c(object, list(DT.sum=sum.dt, meas.full=meas.full, alt=alt, alpha=alpha, p.sig=p.sig))
perm.sum$measure <- measure
class(perm.sum) <- c('summary.brainGraph_permute', class(perm.sum))
perm.sum
}
#' @aliases summary.brainGraph_permute
#' @method print summary.brainGraph_permute
#' @export
print.summary.brainGraph_permute <- function(x, ...) {
print_title_summary('Permutation analysis')
cat('# of permutations:', prettyNum(x$N, ','), '\n')
cat('Level: ', x$level, '\n')
cat('Graph metric: ', x$meas.full, '\n')
if (isTRUE(x$auc)) cat('Area-under-the-curve (AUC) calculated across', length(x$densities), 'densities:\n', x$densities, '\n')
symb <- switch(x$alt, two.sided='!=', greater='>', less='<')
alt <- sprintf('%s - %s %s 0', x$Group[1L], x$Group[2L], symb)
cat('Alternative hypothesis: ', '\t', alt, '\n')
cat('Alpha: ', x$alpha, '\n\n')
if (with(x, dim(DT.sum[get(p.sig) < alpha])[1L]) == 0L) {
cat('No significant results!\n')
} else {
clp <- 100 * (1 - x$alpha)
setnames(x$DT.sum, c('ci.low', 'ci.high'), paste0(clp, '% CI ', c('low', 'high')))
with(x, print(DT.sum[get(p.sig) < alpha]))
}
invisible(x)
}
#' Plot results from permutation testing
#'
#' @param ptitle Character string specifying a title for the plot (default:
#' \code{NULL})
#' @export
#' @return The \code{plot} method returns a \emph{list} of \code{ggplot} objects
#' @rdname brainGraph_permute
plot.brainGraph_permute <- function(x, measure=x$measure,
alternative=c('two.sided', 'less', 'greater'),
alpha=0.05, p.sig=c('p', 'p.fdr'), ptitle=NULL, ...) {
densities <- Group <- sig <- trend <- yloc <- obs <- mylty <- ci.low <- ci.high <-
variable <- value <- reg.num <- region <- perm.diff <- obs.diff <- plot2 <- reg.num2 <- NULL
p.sig <- match.arg(p.sig)
perm.sum <- summary(x, measure=measure, alternative=alternative, alpha=alpha)
measure <- perm.sum$measure
sum.dt <- perm.sum$DT.sum
if (is.null(ptitle)) ptitle <- perm.sum$meas.full
ylabel2 <- sprintf('Observed and permutation difference (%s - %s)', x$Group[1L], x$Group[2L])
# GRAPH LEVEL
#-------------------------------------
if (x$level == 'graph') {
plot.dt <- melt(sum.dt, id.vars=setdiff(names(sum.dt), paste0(measure, '.', x$Group)),
variable.name='Group', value.name='obs')
plot.dt[, Group := factor(Group, labels=x$Group)]
idvars <- c('densities', 'region', 'p', 'Group', 'obs')
if (isFALSE(x$auc)) idvars <- c(idvars, 'p.fdr')
plot.dt <- melt(plot.dt, id.vars=idvars)
plot.dt[, c('sig', 'trend') := '']
plot.dt[get(p.sig) < alpha, sig := '*']
plot.dt[get(p.sig) >= alpha & get(p.sig) < 2 * alpha, trend := '*']
plot.dt[, yloc := extendrange(obs, f=0.07)[1L]]
plot.dt[, mylty := 0]
plot.dt[variable == 'obs.diff', mylty := 1]
plot.dt[variable %in% c('ci.low', 'ci.high'), mylty := 2]
plot.dt[variable == 'perm.diff', mylty := 3]
# 'base' plotting
if (!requireNamespace('ggplot2', quietly=TRUE)) {
# 1. Line plot of observed group values
gID <- getOption('bg.group')
grps <- paste(perm.sum$measure, x$Group, sep='.')
ymin <- plot.dt[1L, yloc]
ymax <- plot.dt[, max(obs)]
par(mar=c(8.6, 4.1, 4.1, 2.1), xpd=TRUE)
plot.dt[variable == 'obs.diff' & get(gID) == grps[1L],
plot(densities, obs, type='l', col=plot.cols[1L],
main=perm.sum$meas.full, xlab='Density', ylab='Observed values')]
plot.dt[variable == 'obs.diff' & get(gID) == grps[2L],
lines(densities, obs, col=plot.cols[2L])]
plot.dt[variable == 'obs.diff' & sig == '*',
text(densities, yloc, sig, cex=1.75, col='red')]
plot.dt[variable == 'obs.diff' & trend == '*',
text(densities, yloc, trend, cex=1.75, col='blue')]
legend('bottom', title=gID, sub(paste0(perm.sum$measure, '.'), '', grps),
fill=plot.cols[1L:2L], inset=c(0, -0.35), horiz=TRUE)
# 2. Line plot of observed difference and CI's
ymin <- plot.dt[variable == 'ci.low', min(value)]
ymax <- plot.dt[variable == 'ci.high', max(value)]
par(mar=c(9.1, 4.1, 4.1, 2.1), xpd=TRUE)
plot.dt[variable == 'obs.diff' & get(gID) == grps[1L],
plot(densities, value, type='l', col='red', ylim=c(ymin, ymax),
main=perm.sum$meas.full, xlab='Density', ylab=ylabel2)]
plot.dt[variable == 'obs.diff' & get(gID) == grps[1L],
points(densities, value, col='red', pch=19)]
plot.dt[variable == 'ci.high' & get(gID) == grps[1L],
lines(densities, value, lty=2)]
plot.dt[variable == 'ci.low' & get(gID) == grps[1L],
lines(densities, value, lty=2)]
plot.dt[variable == 'perm.diff' & get(gID) == grps[1L],
lines(densities, value, lty=3)]
legend('bottom', c('Observed diff.', '95% CI', 'Mean permutation diff.'),
lty=1L:3L, inset=c(0, -0.4))
return(invisible(x))
# 'ggplot2' plotting
} else {
# 1. Line plot of observed group values
p <- ggplot2::ggplot(plot.dt, ggplot2::aes(x=densities))
plot1 <- p +
ggplot2::geom_line(ggplot2::aes(y=obs, col=Group)) +
ggplot2::geom_text(ggplot2::aes(y=yloc, label=sig), col='red', size=8) +
ggplot2::geom_text(ggplot2::aes(y=yloc, label=trend), col='blue', size=8) +
ggplot2::theme(legend.position='bottom')
# 2. Line plot of observed difference and CI's
plot2 <- p +
ggplot2::geom_line(data=plot.dt[variable =='obs.diff'], ggplot2::aes(y=value, lty=factor(mylty)), col='red') +
ggplot2::geom_point(data=plot.dt[variable == 'obs.diff'], ggplot2::aes(y=value), col='red', size=3) +
ggplot2::geom_line(data=plot.dt[variable == 'perm.diff'], ggplot2::aes(y=value, lty=factor(mylty))) +
ggplot2::geom_line(data=plot.dt[variable == 'ci.low'], ggplot2::aes(y=value, lty=factor(mylty))) +
ggplot2::geom_line(data=plot.dt[variable == 'ci.high'], ggplot2::aes(y=value, lty=factor(mylty))) +
ggplot2::scale_linetype_manual(name='',
labels=c('Observed diff.', '95% CI', 'Mean permutation diff.'),
values=1:3) +
ggplot2::theme(legend.position='bottom', legend.spacing.x=ggplot2::unit(9, 'pt'),
legend.key=ggplot2::element_rect(fill='white'),
legend.background=ggplot2::element_rect(fill='gray79'))
ylabel1 <- 'Observed values'
plot1 <- plot1 +
ggplot2::labs(title=ptitle, y=ylabel1, x='Density') +
ggplot2::theme(plot.title=ggplot2::element_text(hjust=0.5, face='bold'))
plot2 <- plot2 +
ggplot2::labs(title=ptitle, y=ylabel2, x='Density') +
ggplot2::theme(plot.title=ggplot2::element_text(hjust=0.5, face='bold'))
return(list(plot1, plot2))
}
# VERTEX LEVEL
#-------------------------------------
} else {
plot.dt <- droplevels(sum.dt[get(p.sig) < alpha])
if (dim(plot.dt)[1L] == 0L) stop('No significant results!')
# 'base' plotting
if (!requireNamespace('ggplot2', quietly=TRUE)) {
ymin <- plot.dt[, min(c(obs.diff, ci.low))]
ymax <- plot.dt[, max(c(obs.diff, ci.high))]
plot.dt[, reg.num2 := as.numeric(as.factor(region))]
dotplot(perm.diff ~ region | as.factor(densities),
xlab='Region', ylab=ylabel2, main=perm.sum$meas.full,
data=plot.dt, ylim=extendrange(c(ymin, ymax)),
x0=plot.dt$reg.num2, y0=plot.dt$ci.low, y1=plot.dt$ci.high, yobs=plot.dt$obs.diff,
scales=list(x=list(relation='free')),
panel=function(x, y, x0, y0, y1, yobs, subscripts) {
lpoints(x, y, pch=19)
panel.abline(h=0, lty=3)
larrows(x0=x0[subscripts], y0=y0[subscripts], x1=x0[subscripts], y1=y1[subscripts],
code=3, angle=90, length=0.1)
lsegments(x0[subscripts] - 0.1, yobs[subscripts],
x0[subscripts] + 0.1, yobs[subscripts],
col='red', lwd=3)
})
# 'ggplot2' plotting
} else {
plot.dt[, reg.num := seq_len(.N), by=densities]
plot1 <- ggplot2::ggplot(plot.dt, ggplot2::aes(x=region)) +
ggplot2::geom_point(ggplot2::aes(y=perm.diff)) +
ggplot2::geom_errorbar(ggplot2::aes(ymin=ci.low, ymax=ci.high)) +
ggplot2::geom_segment(ggplot2::aes(x=reg.num - .25, xend=reg.num + .25, y=obs.diff, yend=obs.diff),
col='red', size=1.25) +
ggplot2::facet_wrap(~ densities, scales='free') +
ggplot2::labs(title=ptitle, y=ylabel2, x='Region') +
ggplot2::theme(plot.title=ggplot2::element_text(hjust=0.5, face='bold'))
return(list(plot1, plot2))
}
}
}
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.