Package 'brainGraph'

Threshold additional set of matrices

Description

apply_thresholds thresholds an additional set of matrices (e.g., FA-weighted matrices for DTI tractography) based on the matrices that have been returned from create_mats. This ensures that the same connections are present in both sets of matrices.

Usage

apply_thresholds(sub.mats, group.mats, W.files, inds)

Arguments

sub.mats	List (length equal to number of thresholds) of numeric arrays (3-dim) for all subjects
group.mats	List (length equal to number of thresholds) of numeric arrays (3-dim) for group-level data
W.files	Character vector of the filenames of the files with connectivity matrices
inds	List (length equal to number of groups) of integers; each list element should be a vector of length equal to the group sizes

Details

The argument W.files accepts the same formats as A.files; see create_mats for details.

Value

List containing:

W	A 3-d array of the raw connection matrices
W.norm.sub	List of 3-d arrays of the normalized connection matrices for all given thresholds
W.norm.mean	List of 3-d arrays of the normalized connection matrices averaged for each group

Author(s)

Christopher G. Watson, [email protected]

Examples

## Not run: 
  W.mats <- apply_thresholds(A.norm.sub, A.norm.mean, f.W, inds)

## End(Not run)

---

Atlas helper functions

Description

guess_atlas tries to determine which atlas is being used based on the data; i.e., the number of vertices/regions.

as_atlas and create_atlas converts/coerces an object to a a data.table, or creates one, that is compatible with brainGraph.

Usage

guess_atlas(x)

as_atlas(object)

create_atlas(regions, coords, lobes, hemis, other = NULL)

Arguments

x, object	An object to test or convert to an atlas data.table
regions	Character vector of region names
coords	Numeric matrix of spatial coordinates; must have 3 columns
lobes	Character or factor vector of lobe membership
hemis	Character or factor vector of hemisphere membership. There should probably not be more than 3 unique elements (for left, right, and bi-hemispheric regions)
other	A named list of vectors with other data. The names of the list will become column names in the return object.

Value

guess_atlas - Character string; either the matched atlas or NA

as_atlas and create_atlas return a data.table that conforms to other atlases in the package, or exits with an error.

Guessing the atlas from an object

There are several valid inputs to guess_atlas:

data.table

The atlas will be guessed based on the number of columns (subtracting by 1 if a “Study ID” column is present). This is the same behavior as for data.frame objects, as well.

igraph

The vertex count

brainGraph

If there is a atlas graph-level attribute, it will return that. Otherwise, the vertex count.

matrix,array

The number of rows, which should equal the number of columns if the input is a connectivity matrix.

Note that this will only work properly for atlases that are currently in the package. If you are using a custom atlas and you receive errors, please open an issue on GitHub.

Coercing to an atlas

There are several things as_atlas tries to do to make it work without error:

• Coerce the object to data.table

• Add a column of integers named index

• Change columns named 'x', 'y', or 'z' to have .mni at the end

• Convert the lobe and hemi columns to be factors

Examples

my_atlas <- data.frame(name=paste('Region', 1:10), x.mni=rnorm(10),
  y.mni=rnorm(10), z.mni=rnorm(10),
  lobe=rep(c('Frontal', 'Parietal', 'Temporal', 'Occipital', 'Limbic'), 2),
  hemi=c(rep('L', 5), rep('R', 5)))
my_atlas2 <- as_atlas(my_atlas)
str(my_atlas)
str(my_atlas2)
regions <- paste('Region', 1:10)
xyz <- matrix(rnorm(30), nrow=10, ncol=3)
lobe <- rep(c('Frontal', 'Parietal', 'Temporal', 'Occipital', 'Limbic'), 2)
hemi <- c(rep('L', 5), rep('R', 5))
other <- list(network=rep(c('Default mode', 'Task positive'), 5))
my_atlas <- create_atlas(regions, xyz, lobe, hemi, other)
str(my_atlas)

---

Set graph, vertex, and edge attributes common in MRI analyses

Description

set_brainGraph_attr is a convenience function that sets a number of graph, vertex, and edge attributes for a given graph object. Specifically, it calculates measures that are common in MRI analyses of brain networks.

Usage

set_brainGraph_attr(g, type = c("observed", "random"),
  use.parallel = TRUE, A = NULL, xfm.type = c("1/w", "-log(w)",
  "1-w", "-log10(w/max(w))", "-log10(w/max(w)+1)"),
  clust.method = "louvain")

xfm.weights(g, xfm.type = c("1/w", "-log(w)", "1-w", "-log10(w/max(w))",
  "-log10(w/max(w)+1)"), invert = FALSE)

Arguments

g	A graph object
type	Character string indicating the type of graphs. Default: observed
use.parallel	Logical indicating whether to use foreach. Default: TRUE
A	Numeric matrix; the (weighted) adjacency matrix, which can be used for faster calculation of local efficiency. Default: NULL
xfm.type	Character string specifying how to transform the weights. Default: 1/w
clust.method	Character string indicating which method to use for community detection. Default: 'louvain'
invert	Logical indicating whether or not to invert the transformation. Default: FALSE

Details

Including type='random' in the function call will reduce the number of attributes calculated. It will only add graph-level attributes for: clustering coefficient, characteristic path length, rich club coefficient, global efficiency, and modularity.

Value

A graph object with the following attributes:

Graph-level	Density, connected component sizes, diameter, # of triangles, transitivity, average path length, assortativity, global & local efficiency, modularity, vulnerability, hub score, rich-club coefficient, # of hubs, edge asymmetry
Vertex-level	Degree, strength; betweenness, eigenvector, and leverage centralities; hubs; transitivity (local); k-core, s-core; local & nodal efficiency; color (community, lobe, component); membership (community, lobe, component); gateway and participation coefficients, within-module degree z-score; vulnerability; and coordinates (x, y, and z)
Edge-level	Color (community, lobe, component), edge betweenness, Euclidean distance (in mm), weight (if weighted)

xfm.weights returns the same graph object, with transformed edge weights plus a graph attribute (xfm.type) recording the method of transformation

Negative edge weights

If there are any negative edge weights in the graph, several of the distance-based metrics will not be calculated, because they can throw errors which is undesirable when processing a large dataset. The metrics are: local and nodal efficiency, diameter, characteristic path length, and hubness.

Transforming edge weights

For distance-based measures, it is important to transform the edge weights so that the strongest connections are re-mapped to having the lowest weights. Then you may calculate e.g., the shortest path length which will include the strongest connections.

xfm.type allows you to choose from 5 options for transforming edge weights when calculating distance-based metrics (e.g., shortest paths). There is no “best-practice” for choosing one over the other, but the reciprocal is probably most common.

1/w

reciprocal (default)

-log(w)

the negative (natural) logarithm

1-w

subtract weights from 1

-log10(w/max(w))

negative (base-10) log of normalized weights

-log10(w/max(w)+1)

same as above, but add 1 before taking the log

To transform the weights back to original values, specify invert=TRUE.

Community detection

clust.method allows you to choose from any of the clustering (community detection) functions available in igraph. These functions begin with cluster_; the function argument should not include this leading character string. There are a few possibilities, depending on the value and the type of input graph:

1. By default, louvain is used, calling cluster_louvain

2. Uses spinglass if there are any negative edges and/or the selected method is spinglass

3. Uses walktrap if there are any negative edge weights and any other method (besides spinglass) is selected

4. Automatically transforms the edge weights if edge_betweenness is selected and the graph is weighted, because the algorithm considers edges as distances

Author(s)

Christopher G. Watson, [email protected]

See Also

components, diameter, centr_betw, betweenness, centr_eigen, transitivity, distances, assortativity, coreness, communities, knn

---

Bootstrapping for global graph measures

Description

Perform bootstrapping to obtain groupwise standard error estimates of a global graph measure.

The plot method returns two ggplot objects: one with shaded regions based on the standard error, and the other based on confidence intervals (calculated using the normal approximation).

Usage

brainGraph_boot(densities, resids, R = 1000, measure = c("mod",
  "E.global", "Cp", "Lp", "assortativity", "strength", "mod.wt",
  "E.global.wt"), conf = 0.95, .progress = getOption("bg.progress"),
  xfm.type = c("1/w", "-log(w)", "1-w", "-log10(w/max(w))",
  "-log10(w/max(w)+1)"))

## S3 method for class 'brainGraph_boot'
summary(object, ...)

## S3 method for class 'brainGraph_boot'
plot(x, ..., alpha = 0.4)

Arguments

densities	Numeric vector of graph densities to loop through
resids	An object of class brainGraph_resids (the output from get.resid)
R	Integer; the number of bootstrap replicates. Default: 1e3
measure	Character string of the measure to test. Default: mod
conf	Numeric; the level for calculating confidence intervals. Default: 0.95
.progress	Logical indicating whether or not to show a progress bar. Default: getOption('bg.progress')
xfm.type	Character string specifying how to transform the weights. Default: 1/w
object, x	A brainGraph_boot object
...	Unused
alpha	A numeric indicating the opacity for the confidence bands

Details

The confidence intervals are calculated using the normal approximation at the $100 \times conf$% level (by default, 95%).

For getting estimates of weighted global efficiency, a method for transforming edge weights must be provided. The default is to invert them. See xfm.weights.

Value

brainGraph_boot – an object of class brainGraph_boot containing some input variables, in addition to a list of boot objects (one for each group).

plot – list with the following elements:

se	A ggplot object with ribbon representing standard error
ci	A ggplot object with ribbon representing confidence intervals

Author(s)

Christopher G. Watson, [email protected]

See Also

boot, boot.ci

Other Group analysis functions: GLM, Mediation, NBS, brainGraph_permute, mtpc

Other Structural covariance network functions: IndividualContributions, Residuals, brainGraph_permute, corr.matrix, import_scn, plot_volumetric

Examples

## Not run: 
boot.E.global <- brainGraph_boot(densities, resids.all, 1e3, 'E.global')

## End(Not run)

---

Coordinates for data from brain atlases

Description

Datasets containing spatial coordinates for: the original AAL atlases, the newer AAL2 atlases, Freesurfer atlases, Brainsuite, Craddock200, Dosenbach160, Harvard-Oxford, and LONI probabilistic brain atlas. In addition to coordinates, there are indices for the major lobes and hemispheres of the brain, the class variable (for Destrieux atlases), functional networks (for Dosenbach, Power, and Gordon atlases; plus the Yeo network labels for the Brainnetome atlas).

Usage

aal116

aal90

aal2.120

aal2.94

destrieux

destrieux.scgm

dk

dk.scgm

dkt

dkt.scgm

brainsuite

craddock200

dosenbach160

hoa112

lpba40

hcp_mmp1.0

power264

brainnetome

gordon333

Format

A data frame with 90 or 116 (for the original AAL atlases), 94 or 120 (for the newer AAL2 atlases), 148 or 162 (for Destrieux), 68 or 82 (for DK), 62 or 76 (for DKT), 74 (Brainsuite), 200 (Craddock), 160 (Dosenbach), 112 (Harvard-Oxford), 40 (LONI), 246 (Brainnetome), 360 (HCP), 264 (Power), or 333 (Gordon) observations on (some of) the following 19 variables:

name

a character vector of region names

x.mni

a numeric vector of x-coordinates (in MNI space)

y.mni

a numeric vector of y-coordinates (in MNI space)

z.mni

a numeric vector of z-coordinates (in MNI space)

lobe

a factor with some of levels Frontal Parietal Temporal Occipital Insula Limbic Cingulate SCGM Cerebellum (for aal116 and aal2.120) and Brainstem (for craddock200)

hemi

a factor with levels L R and B (for dosenbach160)

index

a numeric vector

name.full

a character vector of full region names, for the DK and DKT atlases

class

a factor with levels G G_and_S S, for the Destrieux atlases

network

(dosenbach160) a factor with levels default fronto-parietal cingulo-opercular sensorimotor cerebellum occipital

gyrus

(brainnetome) Abbreviated names of gyri/regions (including subcortical), with 24 unique values

gyrus.full

(brainnetome) Full names of gyrus

subregion

(brainnetome) Abbreviated names of subregions (including subdivisions of subcortical gray matter)

subregion.full

(brainnetome) Full names of subregion

Yeo_7network

(brainnetome) Factor with 8 levels consisting of SCGM plus the 7 networks from Yeo et al.

Yeo_17network

(brainnetome) Factor with 18 levels consisting of SCGM plus the 17 networks from Yeo et al.

area

(HCP) a factor with 23 cortical areas

Anatomy

(power264) Full region/gyrus names for the Power atlas; contains 53 unique regions

Brodmann

(power264) Integer values for Brodmann areas

Note

Use of the HCP parcellation is subject to the terms at https://balsa.wustl.edu/WN56. In particular: "I will acknowledge the use of WU-Minn HCP data and data derived from WU-Minn HCP data when publicly presenting any results or algorithms that benefitted from their use."

Region names in the gordon333 atlas were chosen to match those of the hcp_mmp1.0 atlas. Many were determined from the coordinates (using FSL's atlasquery), while the rest were entered manually by me. The lobe values were matched to the HCP atlas, as well.

Source

https://neuroimaging-core-docs.readthedocs.io/en/latest/pages/atlases.html

References

Tzourio-Mazoyer, N. and Landeau, B. and Papathanassiou, D. and Crivello, F. and Etard, O. and Delcroix, N. and Mazoyer, B. and Joliot, M. (2002) Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage, 15(1), 273–289. doi:10.1006/nimg.2001.0978

Rolls, E.T. and Joliot, M. and Tzourio-Mazoyer, N. (2015) Implementation of a new parcellation of the orbitofrontal cortex in the automated anatomical labelling atlas. NeuroImage, 122, 1–5. doi:10.1016/j.neuroimage.2015.07.075

Destrieux, C. and Fischl, B. and Dale, A. and Halgren E. (2010) Automatic parcellation of human cortical gyri and sulci using standard anatomic nomenclature. NeuroImage, 53(1), 1–15. doi:10.1016/j.neuroimage.2010.06.010

Desikan, R.S. and Segonne, F. and Fischl, B. et al. (2006) An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. NeuroImage, 31, 968–980. doi:10.1016/j.neuroimage.2006.01.021

Klein, A. and Tourville, J. (2012) 101 labeled brain images and a consistent human cortical labeling protocol. Front Neurosci, 6. doi:10.3389/fnins.2012.00171

Shattuck, D.W. and Leahy, R.M. (2002) BrainSuite: an automated cortical surface identification tool. Medical Image Analysis, 8(2), 129–142.

Pantazis, D. and Joshi, A.A. and Jintao, J. and Shattuck, D.W. and Bernstein, L.E. and Damasio, H. and Leahy, R.M. (2009) Comparison of landmark-based and automatic methods for cortical surface registration. NeuroImage, 49(3), 2479–2493.

Craddock, R.C. and James, G.A. and Holtzheimer, P.E. and Hu, X.P. and Mayberg, H.S. (2012) A whole brain fMRI atlas generated via spatially constrained spectral clustering. Human Brain Mapping, 33, 1914–1928. doi:10.1002/hbm.21333

Dosenbach, N.U. and Nardos, B. and Cohen, A.L. and Fair, D.A. and Power, J.D. and Church, J.A. and Nelson, S.M. and Wig, G.S. and Vogel, A.C. and Lessov-Schlaggar, C.N. and Barnes, K.A. (2010) Prediction of individual brain maturity using fMRI. Science, 329(5997), 1358–1361.

Makris, N. and Goldstein, J.M. and Kennedy, D. et al. (2006) Decreased volume of left and total anterior insular lobule in schizophrenia. Schizophr Res, 83(2-3), 155–171.

Shattuck, D.W. and Mirza, M. and Adisetiyo, V. and Hojatkashani, C. and Salamon, G. and Narr, K.L. and Poldrack, R.A. and Bilder, R.M. and Toga, A.W. (2008) Construction of a 3D probabilistic atlas of human cortical structures. NeuroImage, 39(3), 1064–1080. doi:10.1016/j.neuroimage.2007.09.031

Glasser, M.F. and Coalson, T.S. and Robinson, E.C. and Hacker, C.D. and Harwell, J. and Yacoub, E. and Ugurbil, K. and Andersson, J. and Beckmann, C.F. and Jenkinson, M. and Smith, S.M. and van Essen, D.C. (2016) A multi-modal parcellation of human cerebral cortex. Nature, 536, 171–178. doi:10.1038/nature18933. PMID: 27437579.

Power, J.D. and Cohen, A.L. and Nelson, S.M. and Wig, G.S. and Barnes, K.A. and Church, J.A. and Vogel, A.C. and Laumann, T.O. and Miezin, F.M. and Schlaggar, B.L. and Petersen, S.E. (2011) Functional network organization of the human brain. Neuron, 72(4), 665–678. doi:10.1016/j.neuron.2011.09.006

Fan, L. and Li, H. and Zhuo, J. and Zhang, Y. and Wang, J. and Chen, L. and Yang, Z. and Chu, C. and Xie, S. and Laird, A.R. and Fox, P.T. and Eickhoff, S.B. and Yu, C. and Jiang, T (2016) The Human Brainnetome Atlas: A New Brain Atlas Based on Connectional Architecture. Cerebral Cortex, 26(8), 3508–3526. doi:10.1093/cercor/bhw157

Gordon, E.M. and Laumann, T.O. and Adeyemo, B. and Huckins, J.F. and Kelley, W.M. and Petersen, S.E. (2014) Generation and Evaluation of a Cortical Area Parcellation from Resting-State Correlations. Cerebral Cortex, 26(1), 288–303. doi:10.1093/cercor/bhu239

---

Default options for brainGraph

Description

brainGraph is a package for performing graph theory analysis of brain MRI data.

Package options

brainGraph uses the following options to configure behavior:

• bg.subject_id: character string specifying the name your project/study uses as a subject identifier. All imported data (e.g., covariates tables) MUST have a column matching this. One possible alternative is 'participant_id', recommended by BIDS. Default: 'Study.ID'

• bg.group: character string specifying the name your project/study uses as a group identifier. All imported data (e.g., covariates tables) MUST have a column matching this. One possible alternative is 'group', recommended by BIDS. Default: 'Group'

• bg.session: character string specifying the name your project/study uses as a “time” or session identifier, in the case of longitudinal studies. All imported data (e.g., covariates tables) MUST have a column matching this. One possible alternative is 'session_id', recommended by BIDS. Default: 'Time'

• bg.progress: logical indicating whether to show progress bars for functions that provide the option. Default: TRUE

• bg.ncpus: integer indicating the number of cores to use for parallel operations. Only used if you have not already registered a parallel backend (see Chapter 5 of the User Guide or https://github.com/cwatson/brainGraph/blob/master/README.md for examples). Default: 2L

---

Permutation test for group difference of graph measures

Description

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.

Usage

brainGraph_permute(densities, resids, N = 5000, 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)

## S3 method for class 'brainGraph_permute'
summary(object, measure = object$measure,
  alternative = c("two.sided", "less", "greater"), alpha = 0.05,
  p.sig = c("p", "p.fdr"), ...)

## S3 method for class 'brainGraph_permute'
plot(x, measure = x$measure,
  alternative = c("two.sided", "less", "greater"), alpha = 0.05,
  p.sig = c("p", "p.fdr"), ptitle = NULL, ...)

Arguments

densities	Numeric vector of graph densities
resids	An object of class brainGraph_resids (the output from get.resid)
N	Integer; the number of permutations (default: 5e3)
perms	Numeric matrix of permutations, if you would like to provide your own (default: NULL)
auc	Logical indicating whether or not to calculate differences in the area-under-the-curve of metrics (default: FALSE)
level	A character string for the attribute “level” to calculate differences (default: graph)
measure	A character string specifying the vertex-level metric to calculate, only used if level='vertex' (default: btwn.cent). For the summary method, this is to focus on a single graph-level measure (since multiple are calculated at once).
.function	A custom function you can pass if level='other'
object, x	A brainGraph_permute object (output by brainGraph_permute).
alternative	Character string, whether to do a two- or one-sided test. Default: 'two.sided'
alpha	Numeric; the significance level. Default: 0.05
p.sig	Character string specifying which p-value to use for displaying significant results (default: p)
...	Unused
ptitle	Character string specifying a title for the plot (default: NULL)

Details

If you would like to calculate differences in the area-under-the-curve (AUC) across densities, then specify auc=TRUE.

There are three possible “levels”:

1. graph Calculate modularity (Louvain algorithm), clustering coefficient, characteristic path length, degree assortativity, and global efficiency.

2. 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.

3. 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: g (a list of lists of igraph graph objects); and densities (numeric vector).

Value

An object of class brainGraph_permute with input arguments in addition to:

DT	A data table with permutation statistics
obs.diff	A data table of the observed group differences
Group	Group names

The plot method returns a list of ggplot objects

Author(s)

Christopher G. Watson, [email protected]

See Also

Other Group analysis functions: Bootstrapping, GLM, Mediation, NBS, mtpc

Other Structural covariance network functions: Bootstrapping, IndividualContributions, Residuals, corr.matrix, import_scn, plot_volumetric

Examples

## Not run: 
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)

## End(Not run)

---

brainGraph generic methods

Description

These functions are S3 generics for various brainGraph-defined objects.

groups returns the “Group” graph attribute for each graph or observation in the object.

region.names is a generic method for extracting region names from various brainGraph objects. These are generally convenience functions.

nregions is a generic method for extracting the number of regions from various brainGraph objects.

Usage

## S3 method for class 'brainGraphList'
groups(x)

## S3 method for class 'corr_mats'
groups(x)

region.names(object)

## S3 method for class 'data.table'
region.names(object)

nregions(object)

Arguments

x, object

An object

Details

For a data.table, region.names assumes that it contains a factor column named region.

---

Create a list of brainGraph graphs

Description

make_brainGraphList creates a brainGraphList object, a list containing a set of graphs for all subjects (or group-average graphs) in a study at a specific threshold (or density), in addition to some graph-level attributes common to those graphs.

The [ method will let you subset/slice the graphs for individual subjects and/or groups.

as_brainGraphList coerces a list of graphs to a brainGraphList object. It is assumed that certain metadata attributes – threshold, package version, atlas, imaging modality, edge weighting, and whether they are random graphs – are identical for all graphs in the list.

Usage

make_brainGraphList(x, atlas, type = c("observed", "random"),
  level = c("subject", "group", "contrast"), set.attrs = TRUE,
  modality = NULL, weighting = NULL, threshold = NULL,
  gnames = NULL, ...)

## S3 method for class 'array'
make_brainGraphList(x, atlas, type = c("observed",
  "random"), level = c("subject", "group", "contrast"),
  set.attrs = TRUE, modality = NULL, weighting = NULL,
  threshold = NULL, gnames = NULL, grpNames = NULL, subnet = NULL,
  mode = "undirected", weighted = NULL, diag = FALSE,
  .progress = getOption("bg.progress"), ...)

## S3 method for class 'corr_mats'
make_brainGraphList(x, atlas = x$atlas,
  type = "observed", level = "group", set.attrs = TRUE,
  modality = NULL, weighting = NULL, threshold = x$densities,
  gnames = names(x$r.thresh), grpNames = gnames, mode = "undirected",
  weighted = NULL, diag = FALSE,
  .progress = getOption("bg.progress"), ...)

## S3 method for class 'brainGraphList'
x[i, g = NULL, drop = TRUE]

## S3 method for class 'brainGraphList'
print(x, ...)

is.brainGraphList(x)

## S3 method for class 'brainGraphList'
nobs(object, ...)

as_brainGraphList(g.list, type = c("observed", "random"),
  level = c("subject", "group", "contrast"))

Arguments

x	3-D numeric array of all subjects' connectivity matrices (for a single threshold) or a corr_mats object
atlas	Character string specifying the brain atlas
type	Character string indicating the type of graphs. Default: observed
level	Character string indicating whether the graphs are subject-, group-, or contrast-specific. Default: 'subject'
set.attrs	Logical indicating whether to assign all graph-, vertex-, and edge-level attributes (via set_brainGraph_attr). Default: TRUE
modality	Character string indicating imaging modality (e.g. 'dti'). Default: NULL
weighting	Character string indicating how the edges are weighted (e.g., 'fa', 'pearson', etc.). Default: NULL
threshold	Integer or number indicating the threshold used when “sparsifying” the connectivity matrix (if any). Default: NULL
gnames	Character vector of graph names (e.g., study IDs if level='subject'). Default: NULL
...	Other arguments passed to set_brainGraph_attr
grpNames	Character (or factor) vector of group names. If level == 'group', then you do not need to include this argument (the group names will be the same as gnames). Default: NULL)
subnet	Integer or character vector indicating the vertices to keep, if you are interested in working with a subset of an atlas. By default, all vertices are used.
mode	Character string defining how the matrix should be interpreted. Default: 'undirected'
weighted	Logical specifying whether to create a weighted network
diag	Logical indicating whether to include the diagonal of the connectivity matrix. Default: FALSE
.progress	Logical indicating whether to print a progress bar. Default: getOption('bg.progress')
i	Integer, character, or logical vector for subsetting by subject, or by group (if x$level='group')
g	Integer, character, or logical vector for subsetting by group (if x$level='subject')
drop	If TRUE (the default), then return only the list of graphs; otherwise, subset the graphs and return the entire object
object	A brainGraphList object
g.list	List of graph objects

Details

In addition to creating the initial igraph graphs from the connectivity matrices, then attributes will be calculated and assigned for each graph via set_brainGraph_attr if set.attrs=TRUE. Other arguments can be passed to that function. You may display a progress bar by setting .progress=TRUE.

This object can be considered comparable to a 4-D NIfTI file, particularly that returned by FSL's TBSS “prestats” step since that file contains the FA volumes for all study subjects.

To convert an object with 3 “levels” (i.e., subject-level lists from an older brainGraph version), see the code in the Examples below.

Value

make_brainGraphList returns an object of class brainGraphList with elements:

threshold	The specified threshold/density
version	The versions of R, igraph, and brainGraph used when creating the graphs
atlas	The atlas common to all the graphs
modality	The imaging modality (if supplied)
weighting	A string indicating what edge weights represent (if applicable)
graphs	A named list of brainGraph graphs; the names correspond to the individual graphs' Study IDs

[ – A brainGraphList object (if drop=FALSE) or a list of graphs

Subsetting/extracting

The first index is for subsetting the individual graphs. The second index is for subsetting by group membership and requires that the graphs have a Group graph attribute. When both are included, the first index cannot have length or numeric value greater than the number of remaining subjects after subsetting by group.

If the indexing vector(s) is (are) character, the vector(s) must contain one (or more) of the subject or group names. If logical, its length must equal the number of subjects or groups.

Note

If the input is a corr_mats object, and the extent of the 3-D array is greater than 1, then only the first will be converted to a graph.

Author(s)

Christopher G. Watson, [email protected]

See Also

Other Graph creation functions: Creating_Graphs_GLM, Creating_Graphs, make_ego_brainGraph

Examples

## Not run: 
# Create a list, one for each threshold
g <- vector('list', length(thresholds))
for (i in seq_along(thresholds)) {
  g[[i]] <- make_brainGraphList(A.norm.sub[[i]], thresholds[i], atlas,
      covars.dti$Study.ID, covars.dti$Group, modality='dti', weighting='fa')
}

## End(Not run)
## Not run: 
# Subset the first 10 subjects, irrespective of group
my.bgl[1:10]

# Return object for only 'Control' subjects
my.bgl[, 'Control']

# Return object with graphs from groups 1 and 3
my.bgl[g=c(1, 3), drop=FALSE]

# Subset the first 10 subjects of group 2
my.bgl[1:10, 2]

## End(Not run)
## Not run: 
## Convert old version single-subject graph lists
## g[[1]] is group 1, g[[1]][[1]] is threshold 1, g[[1]][[1]][[1]] is subj. 1
kNumThresholds <- length(g[[1]])
g.l <- vector('list', kNumThresholds)
for (i in seq_len(kNumThresholds)) {
  g.l[[i]] <- as_brainGraphList(do.call(Map, c(c, g))[[i]])
}

## End(Not run)

---

Calculate communicability betweenness centrality

Description

centr_betw_comm calculates the communicability betweenness of the vertices of a graph. The centrality for vertex r is

$\omega_r = \frac{1}{C} \sum_p \sum_q \frac{(e^{\mathbf{A}})_{pq} - (e^{\mathbf{A} + \mathbf{E}(r)})_{pq}}{(e^{\mathbf{A}})_{pq}}$

where $C = (n - 1)^2 - (n - 1)$ is a normalization factor.

Usage

centr_betw_comm(g, A = NULL)

Arguments

g	An igraph graph object
A	Numeric matrix; the adjacency matrix of the input graph. Default: NULL

Value

A numeric vector of the centrality for each vertex

Author(s)

Christopher G. Watson, [email protected]

References

Estrada, E. and Higham, D.J. and Hatano N. (2009) Communicability betweenness in complex networks. Physica A, 388, 764–774. doi:10.1016/j.physa.2008.11.011

See Also

Other Centrality functions: centr_lev

---

Calculate a vertex's leverage centrality

Description

Calculates the leverage centrality of each vertex in a graph.

Usage

centr_lev(g, A = NULL)

Arguments

g	An igraph graph object
A	Numeric matrix; the adjacency matrix of the input graph. Default: NULL

Details

The leverage centrality relates a vertex's degree with the degree of its neighbors. The equation is:

$l_i = \frac{1}{k_i} \sum_{j \in N_i} \frac{k_i - k_j}{k_i + k_j}$

where $k_i$ is the degree of the $i^{th}$ vertex and $N_i$ is the set of neighbors of i. This function replaces NaN with NA (for functions that have the argument na.rm).

Value

A vector of the leverage centrality for all vertices.

Author(s)

Christopher G. Watson, [email protected]

References

Joyce, K.E. and Laurienti P.J. and Burdette J.H. and Hayasaka S. (2010) A new measure of centrality for brain networks. PLoS One, 5(8), e12200. doi:10.1371/journal.pone.0012200

See Also

Other Centrality functions: centr_betw_comm

---

Test if an object is a character vector of numbers

Description

check_sID is a convenience function to test if a vector (typically the subject ID column in a data.table) is a character vector of numbers, a factor vector of numbers, or a numeric vector. If so, it will zero-pad the variable to have equal width.

pad_zeros pads a vector with zeros to avoid issues with ordering a column of integers or integers converted to character.

Usage

check_sID(x)

pad_zeros(x)

Arguments

x	pad_zeros accepts either a vector (numeric or character) or a single integer. check_sID accepts a character, numeric, or factor vector

Details

This function is meant to avoid issues that arise when sorting a vector of numbers that have been converted to character. For example, import_scn automatically reads in the first column (with FreeSurfer outputs this is the column of subject IDs) as a character variable. If the subject IDs had been all numbers/integers, then sorting (i.e., setting the key in a data.table) would be incorrect: e.g., it might be '1', '10', '2', ....

If “x” is a numeric vector, then the resultant string width will be determined by max(x) or x itself if the input is a single integer. For example, if x=10, it will return '01', '02', ..., '10'. If “x” is a character vector, then the output's string width will be max(nchar(x)). For example, if x includes both '1' and '1000', it will return '0001', etc.

Value

check_sID returns either the input vector or a character vector padded with 0

A character vector with zero-padded values

Examples

pad_zeros(10)  # '01' '02' ... '10'
x <- c(1, 10, 100)
pad_zeros(x)   # '001' '010' '100'
x <- as.character(x)
pad_zeros(x)   # '001' '010' '100'

---

Calculate coefficient of variation

Description

coeff_var is a S3 generic that calculates the coefficient of variation, defined as

$CV(x) = \frac{sd(x)}{mean(x)}$

Usage

coeff_var(x, na.rm = FALSE, ...)

## Default S3 method:
coeff_var(x, na.rm = FALSE, ...)

Arguments

x	Numeric vector, matrix, or array
na.rm	Logical indicating whether NA values should be stripped when calculating sums. Default: FALSE
...	Unused

Details

If x is a matrix, it will calculate the CV for each column. If x is a 3D array, it will calculate the coefficient of variation for each row-column combination. If the input dimensions are $n \times n \times r$, a matrix with size $n \times n$ will be returned.

Value

A numeric vector or matrix

---

Calculate communicability

Description

communicability calculates the communicability of a network, a measure which takes into account all possible paths (including non-shortest paths) between vertex pairs.

Usage

communicability(g, weights = NULL)

Arguments

g	An igraph graph object
weights	Numeric vector of edge weights; if NULL (the default), and if the graph has edge attribute weight, then that will be used. To avoid using weights, this should be NA.

Details

The communicability $G_{pq}$ is a weighted sum of the number of walks from vertex p to q and is calculated by taking the exponential of the adjacency matrix A:

$G_{pq} = \sum_{k=0}^{\infty} \frac{(\mathbf{A}^k)_{pq}}{k!} = (e^{\mathbf{A}})_{pq}$

where $k$ is walk length.

For weighted graphs with $D = diag(d_i)$ a diagonal matrix of vertex strength,

$G_{pq} = (e^{\mathbf{D}^{-1/2} \mathbf{A} \mathbf{D}^{-1/2}})_{pq}$

Value

A numeric matrix of the communicability

Author(s)

Christopher G. Watson, [email protected]

References

Estrada, E. and Hatano, N. (2008) Communicability in complex networks. Physical Review E. 77, 036111. doi:10.1103/PhysRevE.77.036111

Crofts, J.J. and Higham, D.J. (2009) A weighted communicability measure applied to complex brain networks. J. R. Soc. Interface. 6, 411–414. doi:10.1098/rsif.2008.0484

---

Contract graph vertices based on brain lobe and hemisphere

Description

Create a new graph after merging vertices within specified groups. By default, groups are brain lobe and hemisphere membership.

Usage

contract_brainGraph(g, vgroup = "lobe.hemi")

Arguments

g	A brainGraph graph object
vgroup	Character string; the name of the vertex attribute to use when contracting the graph. Default: 'lobe.hemi'

Details

The size vertex-level attribute of the resultant graph is equal to the number of vertices in each group. The x-, y-, and z-coordinates of the new graph are equal to the mean coordinates of the vertices per group. The new edge weights are equal to the number of inter-group connections of the original graph.

Value

A new brainGraph graph object with vertex-level attributes representing the mean spatial coordinates, and vertex- and edge-level attributes of color names

Author(s)

Christopher G. Watson, [email protected]

See Also

contract

---

Calculate the p-value for differences in correlation coefficients

Description

Given two sets of correlation coefficients and sample sizes, this function calculates and returns the z-scores and p-values associated with the difference between correlation coefficients. This function was adapted from https://stackoverflow.com/a/14519007/3357706.

Usage

cor.diff.test(r1, r2, n, alternative = c("two.sided", "less", "greater"))

Arguments

r1, r2	Numeric (vector or matrix) of correlation coefficients for both groups
n	Integer vector; number of observations for both groups
alternative	Character string, whether to do a two- or one-sided test. Default: 'two.sided'

Value

A list with elements p and z, the p-values and z-scores for the difference in correlations.

Author(s)

Christopher G. Watson, [email protected]

Examples

## Not run: 
kNumSubjs <- summary(covars$Group)
corr.diffs <- cor.diff.test(corrs$R[, , 1], corrs$R[, , 2], kNumSubjs)
edge.diffs <- t(sapply(which(corr.diffs$p < .05), function(x)
                       mapply('[[',
                              dimnames(corr.diffs$p),
                              arrayInd(x, dim(corr.diffs$p)))
                              ))

## End(Not run)

---

Calculate correlation matrix and threshold

Description

corr.matrix calculates the correlation between all column pairs of a given data frame, and thresholds the resultant correlation matrix based on a given density (e.g., 0.1 if you want to keep only the 10% strongest correlations). If you want to threshold by a specific correlation coefficient (via the thresholds argument), then the densities argument is ignored.

The plot method will plot “heat maps” of the correlation matrices.

Usage

corr.matrix(resids, densities, thresholds = NULL, what = c("resids",
  "raw"), exclude.reg = NULL, type = c("pearson", "spearman"),
  rand = FALSE)

## S3 method for class 'corr_mats'
x[i, g = NULL]

## S3 method for class 'corr_mats'
plot(x, mat.type = c("thresholded", "raw"),
  thresh.num = 1L, ordered = TRUE, order.by = "lobe",
  graphs = NULL, grp.names = NULL, legend.title = NULL, ...)

## S3 method for class 'corr_mats'
region.names(object)

## S3 method for class 'corr_mats'
nregions(object)

Arguments

resids	An object of class brainGraph_resids (the output from get.resid)
densities	Numeric vector indicating the resultant network densities; keeps the top X% of correlations
thresholds	Numeric; absolute correlation value to threshold by (default: NULL)
what	Character string indicating whether to correlate the residuals or the raw structural MRI values (default: 'resids')
exclude.reg	Character vector of regions to exclude (default: NULL)
type	Character string indicating which type of correlation coefficient to calculate (default: 'pearson')
rand	Logical indicating whether the function is being called for permutation testing; not intended for general use (default: FALSE)
x, object	A corr_mats object
i	Integer for subsetting by density/threshold
g	Integer, character, or logical for subsetting by group
mat.type	Character string indicating whether to plot raw or thresholded (binarized) matrices. Default: 'raw'
thresh.num	Integer specifying which threshold to plot (if mat.type='thresholded'). Default: 1L
ordered	Logical indicating whether to order the vertices by some grouping. Default: TRUE
order.by	Character string indicating how to group vertices. Default: 'lobe'
graphs	A brainGraphList object containing graphs with the vertex-level attribute of interest. Default: NULL
grp.names	Character vector specifying the names of each group of vertices. Default: NULL
legend.title	Character string for the legend title. Default is to leave blank
...	Unused

Details

If you wish to exclude regions from your analysis, you can give the indices of their columns with the exclude.reg argument.

By default, the Pearson correlation coefficients are calculated, but you can return Spearman by changing the type argument.

Value

A corr_mats object containing the following components:

R, P	Numeric arrays of correlation coefficients and P-values. The length of the 3rd dimension equals the number of groups
r.thresh	A list of 3-d binary arrays indicating correlations that are above a certain threshold. The length of the list equals the number of groups, and the length of the 3rd dimension equals the number of thresholds/densities.
thresholds	Numeric matrix of the thresholds supplied. The number of columns equals the number of groups.
what	Residuals or raw values
exclude.reg	Excluded regions (if any)
type	Pearson or Spearman
atlas	The brain atlas used
densities	Numeric vector; the densities of the resulting graphs, if you chose to threshold each group to have equal densities.

Plotting correlation matrices

There are several ways to control the plot appearance. First, you may plot the “raw” correlations, or only those of the thresholded (binarized) matrices. Second, you may order the vertices by a given vertex attribute; by default, they will be ordered by lobe, but you may also choose to order by, e.g., network (for the dosenbach160 atlas) or by community membership. In the latter case, you need to pass a brainGraphList object to the graphs argument; each graph in the object must have a vertex attribute specified in order.by. Finally, you can control the legend text with grp.names.

Author(s)

Christopher G. Watson, [email protected]

See Also

rcorr

Other Structural covariance network functions: Bootstrapping, IndividualContributions, Residuals, brainGraph_permute, import_scn, plot_volumetric

Examples

## Not run: 
myResids <- get.resid(lhrh, covars)
corrs <- corr.matrix(myResids, densities=densities)))

## End(Not run)
## Not run: 
corrs <- corr.matrix(myResids, densities)
plot(corrs, order.by='comm', graphs=g.list, grp.names='Community')

## End(Not run)

---

Count number of edges of a brain graph

Description

count_homologous counts the number of edges between homologous regions in a brain graph (e.g. between L and R superior frontal).

count_inter counts the number of edges between and within all vertices in one group (e.g. lobe, hemi, network, etc.).

Usage

count_homologous(g)

count_inter(g, group = c("lobe", "hemi", "network", "class", "gyrus",
  "Yeo_7network", "Yeo_17network", "area", "Brodmann"))

Arguments

g	A brainGraph graph object
group	Character string specifying which grouping to calculate edge counts for. Default: 'lobe'

Value

count_homologous - a named vector of the edge ID's connecting homologous regions

count_inter - a data.table of total, intra-, and inter-group edge counts

Author(s)

Christopher G. Watson, [email protected]

Examples

## Not run: 
g1.lobecounts <- count_inter(g[[1]][[N]], 'lobe')

## End(Not run)

---

Create connection matrices for tractography or fMRI data

Description

create_mats creates arrays from connection matrices (e.g., fdt_network_matrix from FSL or ROICorrelation.txt from DPABI). You may choose to normalize these matrices by the waytotal or region size (tractography), or not at all.

Usage

create_mats(A.files, modality = c("dti", "fmri"), divisor = c("none",
  "waytotal", "size", "rowSums"), div.files = NULL,
  threshold.by = c("consensus", "density", "mean", "consistency", "raw"),
  mat.thresh = 0, sub.thresh = 0.5, inds = list(seq_along(A.files)),
  algo = c("probabilistic", "deterministic"), P = 5000, ...)

Arguments

A.files	Character vector of the filenames with connection matrices
modality	Character string indicating data modality (default: dti)
divisor	Character string indicating how to normalize the connection matrices; either 'none' (default), 'waytotal', 'size', or 'rowSums' (ignored if modality equals fmri)
div.files	Character vector of the filenames with the data to normalize by (e.g. a list of waytotal files) (default: NULL)
threshold.by	Character string indicating how to threshold the data; choose density, mean, or consistency if you want all resulting matrices to have the same densities (default: consensus)
mat.thresh	Numeric (vector) for thresholding connection matrices (default: 0)
sub.thresh	Numeric (between 0 and 1) for thresholding by subject numbers (default: 0.5)
inds	List (length equal to number of groups) of integers; each list element should be a vector of length equal to the group sizes
algo	Character string of the tractography algorithm used (default: 'probabilistic'). Ignored if modality is fmri.
P	Integer; number of samples per seed voxel (default: 5000)
...	Arguments passed to symmetrize

Value

A list containing:

A	A 3-d array of the raw connection matrices
A.norm	A 3-d array of the normalized connection matrices
A.bin	A list of 3-d arrays of binarized connection matrices, one array for each threshold
A.bin.sums	A list of 3-d arrays of connection matrices, with each entry signifying the number of subjects with a connection present; the number of list elements equals the length of mat.thresh, and the extent of the arrays equals the number of groups
A.inds	A list of arrays of binarized connection matrices, containing 1 if that entry is to be included
A.norm.sub	List of 3-d arrays of the normalized connection matrices for all given thresholds
A.norm.mean	List of 3-d arrays of connection matrices averaged for each group

Connection matrix files

The A.files argument is mandatory and may be specified in a few ways:

1. A character vector of the filenames (preferably with full path).

2. A single character string specifying the directory in which all connectivity matrices are located. This will load all files in the directory.

3. A named list in which the names match the arguments to list.files. This will load all files in path that match the pattern argument, if present, and will load all files in child directories if recursive=TRUE. See examples below.

The same options apply to div.files as well.

Thresholding methods

The argument threshold.by has 5 options:

1. consensus Threshold based on the raw (normalized, if selected) values in the matrices. If this is selected, it uses the sub.thresh value to perform “consensus” thresholding.

2. density Threshold the matrices to yield a specific graph density (given by the mat.thresh argument).

3. mean Keep only connections for which the cross-subject mean is at least 2 standard deviations higher than the threshold (specified by mat.thresh)

4. consistency Threshold based on the coefficient of variation to yield a graph with a specific density (given by mat.thresh). The edge weights will still represent those of the input matrices. See Roberts et al. (2017) for more on “consistency-based” thresholding.

5. raw Threshold each subject's matrix individually, irrespective of group membership. Ignores sub.thresh.

The argument mat.thresh allows you to choose a numeric threshold, below which the connections will be replaced with 0; this argument will also accept a numeric vector. The argument sub.thresh will keep only those connections for which at least X% of subjects have a positive entry (the default is 0.5, or 50%).

Author(s)

Christopher G. Watson, [email protected]

References

Roberts, JA and Perry, A and Roberts, G and Mitchell, PB and Breakspear, M (2017) Consistency-based thresholding of the human connectome. NeuroImage. 145, 118–129. doi:10.1016/j.neuroimage.2016.09.053

Examples

## Not run: 
thresholds <- seq(from=0.001, to=0.01, by=0.001)
fmri.mats <- create_mats(f.A, modality='fmri', threshold.by='consensus',
  mat.thresh=thresholds, sub.thresh=0.5, inds=inds)
dti.mats <- create_mats(f.A, divisor='waytotal', div.files=f.way,
  mat.thresh=thresholds, sub.thresh=0.5, inds=inds)

# Specify a directory and filename pattern
conn_files <- list(path='~/data', pattern='.*fdt_network_matrix')
dti.mats <- create_mats(conn_files, ...)

## End(Not run)

---

Create a brainGraph object

Description

make_brainGraph is the main creation function for creating a brainGraph graph object. This is simply an igraph graph object with additional attributes (at all levels). Several of the graph-level attributes serve the purpose of providing metadata on how the connectivity matrices/networks were created.

make_brainGraph.bg_mediate creates a graph only for vertex-level analyses.

make_empty_brainGraph creates an empty undirected brainGraph object with vertex count equal to the atlas specified; i.e., it creates a graph with 0 edges. Typically used to present results from an analysis in which edges don't make sense (e.g., GLM comparing differences in a vertex-level attribute).

Usage

make_brainGraph(x, atlas, type = c("observed", "random"),
  level = c("subject", "group", "contrast"), set.attrs = TRUE,
  modality = NULL, weighting = NULL, threshold = NULL, ...)

## S3 method for class 'igraph'
make_brainGraph(x, atlas, type = c("observed",
  "random"), level = c("subject", "group", "contrast"),
  set.attrs = TRUE, modality = NULL, weighting = NULL,
  threshold = NULL, name = NULL, Group = NULL, subnet = NULL, ...)

## S3 method for class 'matrix'
make_brainGraph(x, atlas, type = c("observed",
  "random"), level = c("subject", "group", "contrast"),
  set.attrs = TRUE, modality = NULL, weighting = NULL,
  threshold = NULL, name = NULL, Group = NULL, subnet = NULL,
  mode = "undirected", weighted = NULL, diag = FALSE, ...)

## S3 method for class 'bg_mediate'
make_brainGraph(x, atlas = x$atlas,
  type = "observed", level = "contrast", set.attrs = FALSE,
  modality = NULL, weighting = NULL, threshold = NULL, ...)

is.brainGraph(x)

## S3 method for class 'brainGraph'
summary(object, print.attrs = c("all", "graph",
  "vertex", "edge", "none"), ...)

make_empty_brainGraph(atlas, type = c("observed", "random"),
  level = c("subject", "group", "contrast"), modality = NULL,
  weighting = NULL, threshold = NULL, name = NULL, Group = NULL,
  ...)

Arguments

x	An igraph graph object, numeric matrix, or bg_mediate object
atlas	Character string specifying the brain atlas
type	Character string indicating the type of graphs. Default: observed
level	Character string indicating whether the graphs are subject-, group-, or contrast-specific. Default: 'subject'
set.attrs	Logical indicating whether to assign all graph-, vertex-, and edge-level attributes (via set_brainGraph_attr). Default: TRUE
modality	Character string indicating imaging modality (e.g. 'dti'). Default: NULL
weighting	Character string indicating how the edges are weighted (e.g., 'fa', 'pearson', etc.). Default: NULL
threshold	Integer or number indicating the threshold used when “sparsifying” the connectivity matrix (if any). Default: NULL
...	Arguments passed to set_brainGraph_attr
name	Character string indicating subject ID or group/contrast name, depending on the level. Default: NULL
Group	Character string indicating group membership. Default: NULL
subnet	Integer or character vector indicating the vertices to keep, if you are interested in working with a subset of an atlas. By default, all vertices are used.
mode	Character string defining how the matrix should be interpreted. Default: 'undirected'
weighted	Logical specifying whether to create a weighted network
diag	Logical indicating whether to include the diagonal of the connectivity matrix. Default: FALSE
object	A brainGraph object
print.attrs	Character string indicating whether or not to list the object's attributes (default: all)

Value

A brainGraph graph object with additional graph-, vertex-, and edge-level attributes (see below).

The method for bg_mediate returns a brainGraph_mediate object, which has extra attributes:

Graph	mediator, treat, outcome, nobs
Vertex	b?.acme, p?.acme, b?.ade, p?.ade, b?.prop, p?.prop, b.tot, p.tot

make_empty_brainGraph – An empty brainGraph graph object

Graph-level attributes

Graph-level attributes added are:

version

The R, brainGraph, and igraph package versions used to create the graph

date

The creation date, from as.POSIXct

atlas

Character string denoting the brain atlas used

type

Character string specifying whether this is an observed or random graph

modality

The imaging modality; you can choose anything you like, but the summary.brainGraph knows about dti, fmri, thickness, area, and volume

weighting

What edge weights represent; you can choose anything you like, but summary.brainGraph knows about fa, sld (streamline density, tractography), pearson, spearman, kendall, and partial (partial correlation coefficient)

threshold

Numeric indicating the threshold used to create the final connectivity matrix (if any)

name

Character string specifying the study ID or group/contrast name, depending on the level argument

Group

Character string specifying the experimental group that the given subject belongs to, or if it is a group-level graph

subnet

Integer vector, if subnet was specified in the call

Vertex attributes

Vertex-level attributes added are:

name

The names of the brain regions in the network

lobe

The names of the major brain lobes for each vertex

hemi

The names of the hemisphere for each vertex (either 'L', 'R', or 'B')

lobe.hemi

The lobe-hemisphere combination (represented as an integer vector)

class

The tissue class (if applicable)

network

The network (if the atlas is dosenbach160)

x,y,z

The spatial coordinates of the (centers-of-mass) brain regions in MNI space

x.mni,y.mni,z.mni

Same as above

color.lobe,color.class,color.network

Colors for vertices of their respective membership

circle.layout

Integer vector indicating the order (going counter-clockwise from the top) for circular layouts

Edge attributes

Edge-level attributes added are:

color.lobe,color.class,color.network

Correspond to the vertex attribute of the same name. Inter-group edges will be colored gray

Specifying a subnetwork

You can create a graph for a subset of an atlas's regions with the subnet argument. This can either be a numeric or character vector. If the input object (either a matrix or an igraph graph) has fewer rows/columns or vertices, respectively, than the atlas then the subnet graph attribute will also be added to the return object. This may occur if, for example, you use make_auc_brainGraph on graphs that were initially created from subnetworks.

See Also

Other Graph creation functions: Creating_Graphs_GLM, brainGraphList, make_ego_brainGraph

Examples

## Not run: 
bg <- make_brainGraph(A, 'dkt', modality='dti', weighting='fa',
  mode='undirected', diag=FALSE, weighted=TRUE)

## End(Not run)

---

Create a graph list with GLM-specific attributes

Description

These methods create a brainGraphList with attributes specific to the results of brainGraph_GLM, mtpc, or NBS. The graphs element of the returned object will contain one graph for each contrast.

Usage

## S3 method for class 'bg_GLM'
make_brainGraphList(x, atlas = x$atlas,
  type = "observed", level = "contrast", set.attrs = FALSE,
  modality = NULL, weighting = NULL, threshold = NULL,
  gnames = x$con.name, ...)

## S3 method for class 'mtpc'
make_brainGraphList(x, atlas = x$atlas,
  type = "observed", level = "contrast", set.attrs = FALSE,
  modality = NULL, weighting = NULL, threshold = NULL,
  gnames = x$con.name, ...)

## S3 method for class 'NBS'
make_brainGraphList(x, atlas, type = "observed",
  level = "contrast", set.attrs = TRUE, modality = NULL,
  weighting = NULL, threshold = NULL, gnames = x$con.name,
  mode = "undirected", weighted = TRUE, diag = FALSE, ...)

Arguments

x	A bg_GLM, mtpc, or NBS object
atlas	Character string specifying the brain atlas to use
type	Character string indicating the type of graphs. Default: observed
level	Character string indicating whether the graphs are subject-, group-, or contrast-specific. Default: 'subject'
set.attrs	Logical indicating whether to assign all graph-, vertex-, and edge-level attributes (via set_brainGraph_attr). Default: TRUE
modality	Character string indicating imaging modality (e.g. 'dti'). Default: NULL
weighting	Character string indicating how the edges are weighted (e.g., 'fa', 'pearson', etc.). Default: NULL
threshold	Integer or number indicating the threshold used when “sparsifying” the connectivity matrix (if any). Default: NULL
gnames	Character vector of graph names (e.g., study IDs if level='subject'). Default: NULL
...	Other arguments passed to set_brainGraph_attr
mode	Character string defining how the matrix should be interpreted. Default: 'undirected'
weighted	Logical specifying whether to create a weighted network
diag	Logical indicating whether to include the diagonal of the connectivity matrix. Default: FALSE

Value

A brainGraphList object, with a graph object for each contrast with additional attributes:

Graph	name (contrast name), outcome (the outcome variable), alpha (the significance level); for MTPC: tau.mtpc, S.mtpc, S.crit, A.crit
Vertex	size2 (t-statistic); size (the t-stat transformed for visualization purposes); p (equal to $1-p$); p.fdr (equal to $1-p_{FDR}$, the FDR-adjusted p-value); effect.size (the contrast of parameter estimates for t-contrasts; the extra sum of squares for F-contrasts); se (the standard error of gamma); A.mtpc, sig (binary indicating whether A.mtpc > A.crit) (for MTPC)

make_brainGraphList.NBS returns graphs with additional attributes:

Vertex	comp (integer vector indicating connected component membership), p.nbs (P-value for each component)
Edge	stat (the test statistic for each connection), p (the P-value)

Note

Only valid for vertex-level and NBS analyses.

See Also

brainGraph_GLM, mtpc, NBS

Other Graph creation functions: Creating_Graphs, brainGraphList, make_ego_brainGraph

---

Calculate an asymmetry index based on edge counts

Description

Calculate an asymmetry index, a ratio of intra-hemispheric edges in the left to right hemisphere of a graph for brain MRI data.

Usage

edge_asymmetry(g, level = c("hemi", "vertex"), A = NULL)

Arguments

g	An igraph graph object
level	Character string indicating whether to calculate asymmetry for each region, or the hemisphere as a whole (default: 'hemi')
A	Numeric matrix; the adjacency matrix of the input graph. Default: NULL

Details

The equation is:

$A = \frac{E_{lh} - E_{rh}}{0.5 \times (E_{lh} + E_{rh})}$

where lh and rh are left and right hemispheres, respectively. The range of this measure is $[-2, 2]$ (although the limits will only be reached if all edges are in one hemisphere), with negative numbers indicating more edges in the right hemisphere, and a value of 0 indicating equal number of edges in each hemisphere.

The level argument specifies whether to calculate asymmetry for each vertex, or for the whole hemisphere.

Value

A data table with edge counts for both hemispheres and the asymmetry index; if level is vertex, the data table will have vcount(g) rows.

Author(s)

Christopher G. Watson, [email protected]

---

Calculate graph global, local, or nodal efficiency

Description

This function calculates the global efficiency of a graph or the local or nodal efficiency of each vertex of a graph.

Usage

efficiency(g, type = c("local", "nodal", "global"), weights = NULL,
  xfm = FALSE, xfm.type = NULL, use.parallel = TRUE, A = NULL,
  D = NULL)

Arguments

g	An igraph graph object
type	Character string; either local, nodal, or global. Default: local
weights	Numeric vector of edge weights; if NULL (the default), and if the graph has edge attribute weight, then that will be used. To avoid using weights, this should be NA.
xfm	Logical indicating whether to transform the edge weights. Default: FALSE
xfm.type	Character string specifying how to transform the weights. Default: 1/w
use.parallel	Logical indicating whether or not to use foreach. Default: TRUE
A	Numeric matrix; the adjacency matrix of the input graph. Default: NULL
D	Numeric matrix; the graph's “distance matrix”

Details

Local efficiency for vertex i is:

$E_{local}(i) = \frac{1}{N} \sum_{i \in G} E_{global}(G_i)$

where $G_i$ is the subgraph of neighbors of i, and N is the number of vertices in that subgraph.

Nodal efficiency for vertex i is:

$E_{nodal}(i) = \frac{1}{N-1} \sum_{j \in G} \frac{1}{d_{ij}}$

Global efficiency for graph G with N vertices is:

$E_{global}(G) = \frac{1}{N(N-1)} \sum_{i \ne j \in G} \frac{1}{d_{ij}}$

where $d_{ij}$ is the shortest path length between vertices i and j. Alternatively, global efficiency is equal to the mean of all nodal efficiencies.

Value

A numeric vector of the efficiencies for each vertex of the graph (if type is local|nodal) or a single number (if type is global).

Author(s)

Christopher G. Watson, [email protected]

References

Latora, V. and Marchiori, M. (2001) Efficient behavior of small-world networks. Phys Rev Lett, 87.19, 198701. doi:10.1103/PhysRevLett.87.198701

Latora, V. and Marchiori, M. (2003) Economic small-world behavior in weighted networks. Eur Phys J B, 32, 249–263. doi:10.1140/epjb/e2003-00095-5

---

Fit General Linear Models at each vertex of a graph

Description

brainGraph_GLM specifies and fits a General Linear Model (GLM) at each vertex for a given vertex measure (e.g. degree) or at the graph-level (e.g., global efficiency). Given a contrast matrix or list of contrast(s), and contrast type (for t- or F-contrast(s), respectively) it will calculate the associated statistic(s) for the given contrast(s).

The summary method prints the results, only for which $p < \alpha$, where alpha comes from the bg_GLM object. “Simple” P-values are used by default, but you may change this to the FDR-adjusted or permutation P-values via the function argument p.sig. You may also choose to subset by contrast.

The plot method plots the GLM diagnostics (similar to that of plot.lm). There are a total of 6 possible plots, specified by the which argument; the behavior is the same as in plot.lm. Please see the help for that function.

The [ method allows you to select observations (i.e., rows of X and y) and independent variables (i.e., columns of X) from a bg_GLM object.

Usage

brainGraph_GLM(g.list, covars, measure, contrasts, con.type = c("t",
  "f"), outcome = NULL, X = NULL, con.name = NULL,
  alternative = c("two.sided", "less", "greater"), alpha = 0.05,
  level = c("vertex", "graph"), permute = FALSE,
  perm.method = c("freedmanLane", "terBraak", "smith", "draperStoneman",
  "manly", "stillWhite"), part.method = c("beckmann", "guttman",
  "ridgway"), N = 5000, perms = NULL, long = FALSE, ...)

## S3 method for class 'bg_GLM'
print(x, ...)

## S3 method for class 'bg_GLM'
summary(object, p.sig = c("p", "p.fdr", "p.perm"),
  contrast = NULL, alpha = object$alpha, digits = max(3L,
  getOption("digits") - 2L), print.head = TRUE, ...)

## S3 method for class 'bg_GLM'
plot(x, region = NULL, which = c(1L:3L, 5L),
  ids = TRUE, ...)

## S3 method for class 'bg_GLM'
x[i, j]

Arguments

g.list	A brainGraphList object
covars	A data.table of covariates
measure	Character string of the graph measure of interest
contrasts	Numeric matrix (for T statistics) or list of matrices (for F statistics) specifying the contrast(s) of interest; if only one contrast is desired, you can supply a vector (for T statistics)
con.type	Character string; either 't' or 'f' (for t or F-statistics). Default: 't'
outcome	Character string specifying the name of the outcome variable, if it differs from the graph metric (measure)
X	Numeric matrix, if you wish to supply your own design matrix. Ignored if outcome != measure.
con.name	Character vector of the contrast name(s); if contrasts has row/list names, those will be used for reporting results
alternative	Character string, whether to do a two- or one-sided test. Default: 'two.sided'
alpha	Numeric; the significance level. Default: 0.05
level	Character string; either vertex (default) or graph
permute	Logical indicating whether or not to permute group labels. Default: FALSE
perm.method	Character string indicating the permutation method. Default: 'freedmanLane'
part.method	Character string; the method of partitioning the design matrix into covariates of interest and nuisance. Default: 'beckmann'
N	Integer; number of permutations to create. Default: 5e3
perms	Matrix of permutations, if you would like to provide your own. Default: NULL
long	Logical indicating whether or not to return all permutation results. Default: FALSE
...	Arguments passed to brainGraph_GLM_design
object, x	A bg_GLM object
p.sig	Character string specifying which P-value to use for displaying significant results (default: p)
contrast	Integer specifying the contrast to plot/summarize; defaults to showing results for all contrasts
digits	Integer specifying the number of digits to display for P-values
print.head	Logical indicating whether or not to print only the first and last 5 rows of the statistics tables (default: TRUE)
region	Character string specifying which region's results to plot; only relevant if level='vertex'. Default: NULL
which	Integer vector indicating which of the 6 plots to print to the plot device. Default: c(1:3, 5)
ids	Logical indicating whether to plot subject ID's for outliers. Otherwise plots the integer index
i	Integer/character vector; the observation number(s) or row names to select or remove
j	Integer/character vector; the design matrix column number(s) or names to select or remove

Details

The measure argument will be the graph- or vertex-level measure of interest. Often, this will serve as the model's outcome (or dependent, or response) variable; i.e., the variable typically denoted by y in GLMs. In other cases, you may wish to choose some other variable as the outcome; e.g., IQ, age, etc. Then you could test for a direct association between the network measure and outcome of interest, or test for another association while adjusting for the network metric. For these applications, you must provide the variable name via the outcome argument. This is analogous to -evperdat in FSL's PALM and to --pvr in FreeSurfer.

Value

An object of class bg_GLM containing some input-specific variables (level, outcome, measure, con.type, contrasts, con.name, alt, alpha, permute, perm.method, part.method, N) in addition to:

DT.Xy	A data table from which the design matrices are created and the outcome variable, for all regions.
X	A named numeric matrix or a 3D array of the design matrix. Rownames are subject IDs, column names are predictor variables, and dimnames along the 3rd dimension are region names (if applicable). This is a 3D array only if outcome != measure and level == 'vertex'.
y	A named numeric matrix of the outcome variable. Rownames are Study IDs and column names are regions. There will be multiple columns only if outcome == measure and level == 'vertex'.
DT	A data table with an entry for each vertex (region) containing statistics of interest
removed.subs	A named integer vector in which the names are subject ID's of those removed due to incomplete data (if any). The integers correspond to the row number in the input covars table.
runX	If outcome != measure and level == 'vertex', this will be a character vector of the regions for which the design matrix is invertible. Otherwise, it is NULL.
runY	Character vector of the regions for which the outcome variable has 0 variability. For example, if level='vertex' and measure='degree', some regions may be disconnected or have the same degree for all subjects.
atlas	Character string of the atlas used (guessed based on the vertex count).
perm	A list containing: null.dist (the null distribution of maximum statistics), thresh (the statistic value corresponding to the $100 \times (1 - \alpha)$th% percentile of the null distribution)

The plot method returns a list of ggplot objects (if installed) or writes the plots to a PDF in the current directory named bg_GLM_diagnostics.pdf

A bg_GLM object with the specified row(s) selected or removed from both X and y, and column(s) selected/removed from X

Design matrix

The GLM's design matrix will often be identical to the model matrix associated with lm objects (if “dummy” coding, the default, is used) and is created from the input data.table and arguments passed to brainGraph_GLM_design. The first column should have the name of getOption('bg.subject_id') and its values must match the name graph-level attribute of the input graphs. The covariates table must be supplied even if you provide your own design matrix X. If level='vertex' and outcome == measure, there will be a single design for all regions but a separate model for each region (since the graph measure varies by region). If level='vertex' and outcome != measure, there will be a separate design (and, therefore, a separate model) for each region even though the outcome is the same in all models.

Contrasts and statistics

Either t- or F-contrasts can be calculated (specified by con.type). Multiple t-contrasts can be specified by passing a multi-row matrix to contrasts. Multiple F-contrasts can be specified by passing a list of matrices; all matrices must have the same number of columns. All F-contrasts are necessarily two-sided; t-contrasts can be any direction, but only one can be chosen per function call. If you choose con.type="f", the calculated effect size is represented by the ESS (“extra sum of squares”), the additional variance explained for by the model parameters of interest (as determined by the contrast matrix). The standard error for F-contrasts is the sum of squared errors of the full model.

Non-parametric permutation tests

You can calculate permutations of the data to build a null distribution of the maximum statistic which corrects for multiple testing. To account for complex designs, the design matrix must be partitioned into covariates of interest and nuisance; the default method is the Beckmann method. The default permutation strategy is that of Freedman & Lane (1983), and is the same as that in FSL's randomise. See randomise.

Note

The [ method is used when calculating studentized residuals and other “leave-one-out” diagnostics, and typically should not be called directly by the user.

Author(s)

Christopher G. Watson, [email protected]

See Also

plot.lm

Other GLM functions: GLM design, GLM fits, mtpc

Other Group analysis functions: Bootstrapping, Mediation, NBS, brainGraph_permute, mtpc

Examples

## Not run: 
conmat <- matrix(c(0, 0, 0, 1), nrow=1)
rownames(conmat) <- 'Control > Patient'

res.lm <- brainGraph_GLM(g[[6]], covars=covars.all[tract == 1],
  measure='strength', contrasts=conmat, alt='greater', permute=TRUE, long=TRUE)

## End(Not run)
## Not run: 
## Save objects and then to multipage PDF
lmPlots <- plot(x)
ggsave('lmPlots.pdf', lmPlots)

## Save all the GLM sub-objects from MTPC analysis
res.mtpc <- mtpc(...)
glmPlots <- lapply(res.mtpc$res.glm, plot, which=1:6)
ml <- marrangeGrob(glmPlots, nrow=1, ncol=1)
ggsave('glmPlots.pdf', ml, width=8.5, height=11)

## End(Not run)

---

Extract basic information from a bg_GLM object

Description

These functions return the terms, term labels, model formula, “case names”, “variable names”, region names, and number of observations for a bg_GLM object. The term labels are used for ANOVA tables.

Usage

## S3 method for class 'bg_GLM'
nobs(object, ...)

## S3 method for class 'bg_GLM'
terms(x, ...)

## S3 method for class 'bg_GLM'
formula(x, ...)

## S3 method for class 'bg_GLM'
labels(object, ...)

## S3 method for class 'bg_GLM'
case.names(object, ...)

## S3 method for class 'bg_GLM'
variable.names(object, ...)

## S3 method for class 'bg_GLM'
region.names(object)

## S3 method for class 'bg_GLM'
nregions(object)

Arguments

...	Unused
x, object	A bg_GLM object

Value

terms returns a named integer list in which the names are the term labels and the list elements are the column(s) of the design matrix for each term. nobs returns an integer. The other functions return character vectors.

Note

formula returns a character string, not a formula object.

---

Create a design matrix for linear model analysis

Description

brainGraph_GLM_design takes a data.table of covariates and returns a design matrix to be used in linear model analysis.

Usage

brainGraph_GLM_design(covars, coding = c("dummy", "effects",
  "cell.means"), factorize = TRUE, binarize = NULL, int = NULL,
  mean.center = FALSE, center.how = c("all", "within-groups"),
  center.by = getOption("bg.group"))

Arguments

covars	A data.table of covariates
coding	Character string indicating how factor variables will be coded. Default: 'dummy'
factorize	Logical indicating whether to convert character columns into factor. Default: TRUE
binarize	Character vector specifying the column name(s) of the covariate(s) to be converted from type factor to numeric. Default: NULL
int	Character vector specifying the column name(s) of the covariate(s) to test for an interaction. Default: NULL
mean.center	Logical indicating whether to mean center non-factor variables. Default: FALSE
center.how	Character string indicating whether to use the grand mean or groupwise means. Default: 'all'
center.by	Character string indicating which grouping variable to use for calculating means (if applicable). Default: 'Group'

Details

There are three different ways to code factors: dummy, effects, or cell-means (chosen by the argument coding). Effects coding is sometimes referred to as deviation coding. Dummy coding is the default when calling lm. To understand the difference between these, see Chapter 8 of the User Guide.

Value

A numeric matrix. Rownames are subject ID's and column names are the variable names. There will be additional attributes recording the coding, factorize, and mean.center function arguments. There will also be attributes for binarize and int if they are not NULL, and center.how and center.by if mean.center=TRUE.

Character variables

The default behavior is to convert all character columns (excluding the Study ID column and any that you list in the binarize argument) to factor variables. To change this, set factorize=FALSE. So, if your covariates include multiple character columns, but you want to convert Scanner to binary instead of a factor, you may still specify binarize='Scanner' and get the expected result. binarize will convert the given factor variable(s) into numeric variable(s), which is performed before centering (if applicable).

Centering

The argument mean.center will mean-center (i.e., subtract the mean of from each variable) any non-factor variables (including any dummy/indicator covariates). This is done after “factorizing” and “binarizing”. If center.how='all', then the “grand mean” will be used; otherwise, the groupwise means will be used. The grouping variable is determined by center.by and is by default 'Group'.

Interactions

int specifies which variables should interact with one another. This argument accepts both numeric/continuous (e.g., Age) and factor variables (e.g., Sex). All interaction combinations will be generated: if you supply 3 variables, all two-way and the single three-way interaction will be generated. This variable must have at least two elements; it is otherwise ignored. It is generally recommended that centering be performed when including interaction terms.

Author(s)

Christopher G. Watson, [email protected]

See Also

Other GLM functions: GLM fits, GLM, mtpc

Examples

## Not run: 
# Recreate design matrix when "outcome == measure"
DT <- res.glm$DT.Xy[region == levels(region)[1L],
                    !c('region', res.glm$outcome),
                    with=FALSE]
X <- do.call(brainGraph_GLM_design, c(list(covars=DT),
                                      attributes(res.glm$X)[-c(1L, 2L)]))
all.equal(X, res.glm$X)

## End(Not run)

---

Fit design matrices to one or multiple outcomes

Description

These are the “base” model-fitting functions that solve the least squares problem to estimate model coefficients, residuals, etc. for brain network data.

fastLmBG_t and fastLmBG_f calculate contrast-based statistics for T or F contrasts, respectively. It accepts any number of contrasts (i.e., a multi-row contrast matrix).

Usage

fastLmBG(X, Y, QR = qr.default(X), Q = qr_Q2(QR, n = n, p = p),
  R = qr_R2(QR, p), n = dim(X)[1L], p = QR$rank, ny = dim(Y)[2L],
  dfR = n - p, XtXinv = inv(QR))

fastLmBG_3d(X, Y, runX, QR = qr(X[, , runX, drop = FALSE]),
  Q = lapply(QR, qr_Q2, n = n, p = p), R = lapply(QR, qr_R2, p),
  n = dim(X)[1L], p = QR[[1L]]$rank, ny = length(runX), dfR = n -
  p, XtXinv = inv(QR))

fastLmBG_3dY(X, Y, runX, QR = qr(X[, , runX, drop = FALSE]),
  Q = lapply(QR, qr_Q2, n = n, p = p), R = lapply(QR, qr_R2, p),
  n = dim(X)[1L], p = QR[[1L]]$rank, ny = length(runX), dfR = n -
  p, XtXinv = inv(QR))

fastLmBG_3dY_1p(X, Y, runX, QR = qr(X[, , runX, drop = FALSE]),
  Q = lapply(QR, qr_Q2, diag(1L, n, 1L), n, 1L), R = lapply(QR,
  function(r) r$qr[1L]), n = dim(X)[1L], p = 1L, ny = length(runX),
  dfR = n - 1L, XtXinv = inv(QR))

fastLmBG_t(fits, contrasts, alternative = c("two.sided", "less",
  "greater"), alpha = NULL)

fastLmBG_f(fits, contrasts, rkC = NULL, nC = length(contrasts))

Arguments

X	Design matrix or 3D array of design matrices
Y	Numeric matrix; there should be 1 column for each outcome variable (so that in a graph-level analysis, this is a column matrix)
QR, Q, R	The QR decomposition(s) and Q and R matrix(es) of the design matrix(es). If X is a 3D array, these should be lists
n, p, ny, dfR	Integers; the number of observations, model rank, number of regions/outcome variables, and residual degrees of freedom
XtXinv	Numeric matrix or array; the inverse of the cross-product of the design matrix(es)
runX	Character vector of the regions for which the design matrix is not singular
fits	List object output by one of the model fitting functions (e.g., fastLmBG)
contrasts	Numeric matrix (for T statistics) or list of matrices (for F statistics) specifying the contrast(s) of interest; if only one contrast is desired, you can supply a vector (for T statistics)
alternative	Character string, whether to do a two- or one-sided test. Default: 'two.sided'
alpha	Numeric; the significance level. Default: 0.05
rkC, nC	Integers; the rank of the contrast matrix and number of contrasts, respectively (for F contrasts)

Value

A list with elements

coefficients	Parameter estimates
rank	Model rank
df.residual	Residual degrees of freedom
residuals	Model residuals
sigma	The residual standard deviation, or root mean square error (RMSE)
fitted.values	Model fitted values
qr	The design matrix QR decomposition(s)
cov.unscaled	The “unscaled covariance matrix”

fastLmBG_t – A multidimensional array with the third dimension equaling the number of contrasts; each matrix contains the contrast of parameter estimates, standard error of the contrast, T-statistics, P-values, FDR-adjusted P-values, and confidence intervals (if alpha is given)

fastLmBG_f – A numeric matrix with columns for the effect size, standard error, F statistic, P-values, and FDR-adjusted P-values

Parameter estimation

These functions use the QR decomposition to calculate the least squares solution which is the same as the base lm function. If we substitute $X = QR$ in the standard normal equations, the equation to be solved reduces to

$X^T X \hat{\beta} = X^T y \Rightarrow R \hat{\beta} = Q^T y$

Since R is an upper-triangular matrix, we can use the backsolve function which is a bit faster than solve. In some cases, the fastLmBG* functions are about as fast or faster (particularly when X is not permuted) as one in which the normal equations are solved directly; additionally, using the QR method affords greater numerical stability.

Different scenarios

There are a few different scenarios for fitting models of the data, with a separate function for each:

fastLmBG

The main function for when there is a single design matrix $X$ and any number of outcome variables $Y$.

fastLmBG_3d

Fits models when there is a different design matrix $X$ for each region and a single outcome variable $Y$, which in this case will be a column matrix.

fastLmBG_3dY

Fits models when there is both a different design matrix $X$ and outcome variable $Y$ for each region. Occurs under permutation for the Freedman-Lane, ter Braak, and Still-White methods.

fastLmBG_3dY_1p

Fits models when there is both a different design and outcome variable for each region, and also when $X$ is a rank-1 matrix (i.e., it has 1 column). Only occurs under permutation with the Still-White method if there is a single regressor of interest.

In the last case above, model coefficients are calculated by simple (i.e., non-matrix) algebra.

Improving speed/efficiency

Speed/efficiency gains will be vast for analyses in which there is a single design matrix $X$ for all regions, there are multiple outcome variables (i.e., vertex-level analysis), and the permutation method chosen does not permute $X$. Specifically, these are Freedman-Lane, ter Braak, and Manly methods. Therefore, the QR decomposition, the $Q$ and $R$ matrices, and the “unscaled covariance matrix” (which is $(X^T X)^{-1}$) only need to be calculated once for the entire analysis. Other functions (e.g., lm.fit) would recalculate these for each permutation.

Furthermore, this (and the other model fitting functions in the package) will likely only work in models with full rank. I sacrifice proper error checking in favor of speed, but hopefully any issues with the model will be identified prior to the permutation step. Finally, the number of observations, model rank, number of outcome variables, and degrees of freedom will not change and therefore do not need to be recalculated (although these probably amount to a negligible speed boost).

In case there are multiple design matrices, or the permutation method permutes the design, then the QR decomposition will need to be calculated each time anyway. For these cases, I use more simplified functions qr_Q2 and qr_R2 to calculate the $Q$ and $R$ matrices, and then the fitted values, residuals, and residual standard deviation are calculated at the same time (whereas lm.fit and others would calculate these each time).

Contrast-based statistics

The contrast of parameter estimates, $\gamma$, for T contrasts is

$\gamma = C \hat{\beta}$

where $C$ is the contrast matrix with size $k \times p$ (where $k$ is the number of contrasts) and $\hat{\beta}$ is the matrix of parameter estimates with size $p \times r$ (where $r$ is the number of regions). For F contrasts, the effect size is the extra sum of squares and is calculated as

$\gamma (C (X^T X)^{-1} C^T)^{-1} \gamma^T$

The standard error of a T contrast is

$\sqrt{\hat{\sigma} (X^T X)^{-1}}$

where $\hat{\sigma}$ is the residual standard deviation of the model and the second term is the unscaled covariance matrix. The standard error for F contrasts is simply the residual sum of squares. P-values and FDR-adjusted P-values (across regions) are also calculated. Finally, if $\alpha$ is provided for T contrasts, confidence limits are calculated.

Author(s)

Christopher G. Watson, [email protected]

See Also

randomise

Other GLM functions: GLM design, GLM, mtpc

---

Influence measures for a bg_GLM object

Description

These functions compute common (leave-one-out) diagnostics for the models in a bg_GLM object.

Usage

## S3 method for class 'bg_GLM'
rstandard(model, type = c("sd.1", "predictive"), ...)

## S3 method for class 'bg_GLM'
rstudent(model, ...)

## S3 method for class 'bg_GLM'
hatvalues(model, ...)

## S3 method for class 'bg_GLM'
cooks.distance(model, ...)

dffits.bg_GLM(model)

## S3 method for class 'bg_GLM'
dfbeta(model, ...)

## S3 method for class 'bg_GLM'
dfbetas(model, ...)

covratio.bg_GLM(model)

## S3 method for class 'bg_GLM'
influence(model, do.coef = TRUE, region = NULL, ...)

Arguments

model	A bg_GLM object
type	The type of standardized residuals. Default: 'sd.1'
...	Unused
do.coef	Logical indicating whether to calculate dfbeta
region	Character string of the region(s) to return results for. Default is to calculate for all regions

Details

The influence method calculates all diagnostics present in lm.influence and influence.measures, consisting of the following functions:

rstandard

Standardized residuals. Choosing type='predictive' returns leave-one-out cross validation residuals. The “PRESS” statistic can be calculated as colSums(resids.p^2)

rstudent

Studentized residuals

hatvalues

The leverage, or the diagonal of the hat/projection matrix

cooks.distance

Cook's distance

dffits.bg_GLM

The change in fitted values when deleting observations

dfbeta

The change in parameter estimates (coefficients) when deleting observations

dfbetas

The scaled change in parameter estimates

covratio.bg_GLM

The covariance ratios, or the change in the determinant of the covariance matrix of parameter estimates when deleting observations

Value

Most influence functions return a numeric matrix in which rownames are Study ID's and column names are regions. dfbeta and dfbetas return a numeric array in which each column is a parameter estimate and the 3rd dimension is for each region. influence returns a list with class infl.bg_GLM and elements:

infmat	Numeric array (like dfbeta) with DFBETAs, DFFITs, covratios, Cook's distance, and hat values
is.inf	Logical array of the same data as infmat; values of TRUE indicate the subject-variable-region combination is an outlier value
f	The model formula
sigma	The leave-one-out residual standard deviation
wt.res	Model residuals

Outlier values

Each variable has a different criterion for determining outliers. In the following: x is the influence variable (for DFBETA, the criterion applies to all DFBETAs); k is the number of columns of the design matrix; dfR is the residual degrees of freedom; and n is the number of observations.

DFBETAs

If $|x| > 1$

DFFITs

If $|x| > 3 \sqrt{k / dfR}$

covratio

If $|1 - x| > (3k / dfR)$

cook

If $F_{k, dfR}(x) > 0.5$

hat

If $x > 3k / n$

The return object of influence has a print method which will list the subjects/variables/regions for which an outlier was detected.

Author(s)

Christopher G. Watson, [email protected]

See Also

GLM

---

Model selection for bg_GLM objects

Description

These functions compute the log-likelihood and Akaike's An Information Criterion (AIC) of a bg_GLM object. See logLik.lm and extractAIC for details.

Usage

## S3 method for class 'bg_GLM'
logLik(object, REML = FALSE, ...)

## S3 method for class 'bg_GLM'
extractAIC(fit, scale = 0, k = 2, ...)

Arguments

object, fit	A bg_GLM object
REML	Logical indicating whether to return the restricted log-likelihood. Default: FALSE
...	Unused
scale	Should be left at its default
k	Numeric; the weight of the equivalent degrees of freedom

Details

The functions AIC and BIC will also work for bg_GLM objects because they each call logLik.

Value

logLik returns an object of class logLik with several attributes. extractAIC returns a numeric vector in which the first element is the equivalent degrees of freedom and the remaining are the AIC's for each region

---

Extract model fit statistics from a bg_GLM object

Description

These functions extract or calculate model fit statistics of a bg_GLM object. These can be found in the output from summary.lm.

Usage

## S3 method for class 'bg_GLM'
coef(object, ...)

## S3 method for class 'bg_GLM'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'bg_GLM'
fitted(object, ...)

## S3 method for class 'bg_GLM'
residuals(object, type = c("response", "partial"), ...)

## S3 method for class 'bg_GLM'
deviance(object, ...)

coeff_determ(object, adjusted = FALSE)

## S3 method for class 'bg_GLM'
df.residual(object, ...)

## S3 method for class 'bg_GLM'
sigma(object, ...)

## S3 method for class 'bg_GLM'
vcov(object, ...)

coeff_table(object, CI = FALSE, level = 0.95)

## S3 method for class 'bg_GLM'
anova(object, region = NULL, ...)

Arguments

object	A bg_GLM object
...	Unused
parm	Vector of parameters to calculate confidence intervals for. Default is to use all parameters
level	The confidence level. Default: 0.95
type	Character string specifying the type of residuals to return. Default: 'response'
adjusted	Logical indicating whether to calculate the adjusted R-squared. Default: FALSE
CI	Logical indicating whether to include confidence intervals of parameter estimates in the coefficient summary table. Default: FALSE
region	Character vector indicating the region(s) to calculate ANOVA statistics for. Default: NULL (use all regions)

Details

These mimic the same functions that operate on lm objects, and include:

coef

Regression coefficients (parameter estimates)

confint

Confidence intervals (by default, 95%) for parameter estimates

fitted

Fitted (mean) values; i.e., the design matrix multiplied by the parameter estimates, $X \hat{\beta}$

residuals

Model residuals; i.e., the response/outcome variable minus the fitted values. Partial residuals can also be calculated

deviance

Model deviance, or the residual sum of squares

coeff_determ

Calculate the coefficient of determination (or $R^2$), adjusted or unadjusted

df.residual

Residual degrees of freedom

sigma

Residual standard deviation, sometimes called the root mean squared error (RMSE)

vcov

Variance-covariance matrix of the model parameters

coeff_table returns model coefficients, standard errors, T-statistics, and P-values for all model terms and regions in a bg_GLM object. This is the same as running summary(x)$coefficients for a lm object.

Value

A named numeric vector, matrix, or array, depending on the function:

coef	Matrix in which rownames are parameter names and column names are regions
fitted, residuals	Matrix in which rownames are Study ID's and column names are regions. If type='partial', an array is returned in which columns are terms and the 3rd dimension are regions
deviance, coeff_determ, sigma	Numeric vector with elements for each region
df.residual	Single integer; the degrees of freedom
confint, vcov, coeff_table	Numeric array; the extent of the third dimension equals the number of regions

anova returns a list of tables of class anova

ANOVA tables

The anova method calculates the so-called Type III test statistics for a bg_GLM object. These standard ANOVA statistics include: sum of squares, mean squares, degrees of freedom, F statistics, and P-values. Additional statistics calculated are: $\eta^2$, partial $\eta^2$, $\omega^2$, and partial $\omega^2$ as measures of effect size.

Note

sigma – The denominator is not the number of observations, but rather the model's residual degrees of freedom.

When calculating partial residuals, the parameter estimates are not re-calculated after removing one of the model terms.

Author(s)

Christopher G. Watson, [email protected]

See Also

GLM, Anova

---

Create a data table with graph global and vertex measures

Description

graph_attr_dt is a helper function that takes a brainGraphList or a list of graphs and creates a data.table of global measures for each graph. Each row will be for a different graph.

vertex_attr_dt is a helper function that creates a data.table in which each row is a vertex and each column is a different network measure (degree, centrality, etc.).

Usage

graph_attr_dt(bg.list)

vertex_attr_dt(bg.list)

Arguments

bg.list

A brainGraphList object, or a list of graph objects

Value

A data.table

See Also

graph_attr, graph_attr_names

vertex_attr, vertex_attr_names, graph_from_data_frame

---

Calculate Euclidean distance of edges and vertices

Description

edge_spatial_dist calculates the Euclidean distance of an igraph graph object's edges. The distances are in mm and based on MNI space. These distances are NOT along the cortical surface, so can only be considered approximations, particularly concerning inter-hemispheric connections. The input graph must have atlas as a graph-level attribute.

vertex_spatial_dist calculates, for each vertex of a graph, the average Euclidean distance across all of that vertex's connections.

Usage

edge_spatial_dist(g)

vertex_spatial_dist(g)

Arguments

g	An igraph graph object

Value

edge_spatial_dist - a numeric vector with length equal to the edge count of the input graph, consisting of the Euclidean distance (in mm) of each edge

vertex_spatial_dist - a named numeric vector with length equal to the number of vertices, consisting of the average distance (in mm) for each vertex

Author(s)

Christopher G. Watson, [email protected]

References

Alexander-Bloch, A.F. and Vertes, P.E. and Stidd, R. et al. (2013) The anatomical distance of functional connections predicts brain network topology in health and schizophrenia. Cerebral Cortex, 23, 127–138. doi:10.1093/cercor/bhr388

---

Calculate vertex hubness

Description

hubness calculates the “hubness” (see reference) of the vertices in a graph. These are vertices which meet at least two of the following four criteria:

1. Have high degree/strength

2. Have high betweenness centrality

3. Have low clustering coefficient

4. Have low average path length

For each criterion, “high” or “low” means “in the top 20%” across all vertices. Vertices meeting any of the criteria get a value of 1 for that metric; these are summed to yield the hubness score which ranges from 0-4. As in the reference article, vertices with a score of 2 or higher are to be considered hubs, although that determination isn't made in this function.

Usage

hubness(g, xfm.type = g$xfm.type, weights = NULL, prop.keep = 0.2)

Arguments

g	An igraph graph object
xfm.type	Character string specifying how to transform the weights. Default: 1/w
weights	Numeric vector of edge weights; if NULL (the default), and if the graph has edge attribute weight, then that will be used. To avoid using weights, this should be NA.
prop.keep	Numeric (between 0 and 1) indicating the proportion of vertices to consider as having a high score. Default: 0.2 (20%)

Value

A numeric vector with the vertices' hubness score

Author(s)

Christopher G. Watson, [email protected]

References

van den Heuvel, M.P. and Mandl, R.C.W. and Stam, C.J. and Kahn, R.S. and Pol, H.E.H. (2010) Aberrant frontal and temporal complex network structure in schizophrenia: a graph theoretical analysis. The Journal of Neuroscience, 30(47), 15915–15926. doi:10.1523/JNEUROSCI.2874-10.2010

---

Import data for structural connectivity analysis

Description

Given a directory, atlas name, and imaging modality/structural metric, this function imports data for structural connectivity analysis. It expects files containing a table of region-wise structural MRI measures (e.g., mean cortical thickness), with one file for each hemisphere. The first column of all files should contain the subject ID; the column name will be changed to the value of getOption('bg.subject_id').

Usage

import_scn(datadir, atlas, modality = "thickness", exclude.subs = NULL,
  custom.atlas = NULL)

Arguments

datadir	The path name of the directory containing the data files
atlas	Character string specifying the atlas in use. For a custom atlas, please specify 'custom', and provide the name to the custom.atlas argument
modality	The structural imaging measure (default: 'thickness')
exclude.subs	Vector indicating the subjects to exclude, if any (default: NULL)
custom.atlas	Character string specifying the name of the R object for the atlas in use, if atlas='custom' was also supplied (default: NULL)

Details

The files should have specific names; the second in the following list is only required for atlases/parcellations that include subcortical gray matter (e.g., dk.scgm).

• ${parcellation}_${hemi}_${modality}.csv for cortical volume, thickness, surface area, or local gyrification index (LGI). Here, ${parcellation} can be aparc, aparc.DKTatlas40, or aparc.a2009s. For example, for cortical thickness with the Desikan-Killiany atlas, the filename should be aparc_lh_thickness.csv. If you are using a custom atlas, see the Note below. The ${hemi} variable is either lh or rh. Finally, ${modality} should be either volume, thickness, area, or lgi.

• asegstats.csv for SCGM volume

Value

A list containing:

atlas	Character string
modality	Character string
lhrh	A data.table of structural MRI measures for both hemispheres
aseg	A data.table of structural MRI measures for subcortical gray matter, if applicable
subs.excluded	Vector of subject ID's that were excluded
subs.missing	Vector of subject ID's that are not present in both the cortical and subcortical tables (if applicable)

Note

When using a custom atlas, the name of the atlas's data.table should match the ${parcellation} portion of the filename (specification shown above). Furthermore, it must conform to the output of Freesurfer's aparcstats2table (and asegstats2table, if applicable). Otherwise, please contact me for inclusion of a different data type.

The subject ID column will be zero-padded (to the left) to avoid issues when the variable is numeric; this ensures that all ID's will have the same number of characters and sorting will be done properly.

Author(s)

Christopher G. Watson, [email protected]

See Also

Other Structural covariance network functions: Bootstrapping, IndividualContributions, Residuals, brainGraph_permute, corr.matrix, plot_volumetric

Examples

## Not run: 
  raw_data <- import_scn('/home/cwatson/data', atlas='dkt',
                         exclude.subs=c('con07', 'con23', 'pat15'))

## End(Not run)

---

Approaches to estimate individual network contribution

Description

loo calculates the individual contribution to group network data for each subject in each group using a “leave-one-out” approach. The residuals of a single subject are excluded, and a correlation matrix is created. This is compared to the original correlation matrix using the Mantel test.

aop calculates the individual contribution using an “add-one-patient” approach. The residuals of a single patient are added to those of a control group, and a correlation matrix is created. This is repeated for all individual patients and each patient group.

The summary method prints the group/region-wise means and standard deviations.

The plot method is only valid for regional contribution estimates, and plots the average regional contribution for each vertex/region.

Usage

loo(resids, corrs, level = c("global", "regional"))

aop(resids, corrs, level = c("global", "regional"), control.value = 1L)

## S3 method for class 'IC'
summary(object, region = NULL, digits = max(3L,
  getOption("digits") - 2L), ...)

## S3 method for class 'IC'
plot(x, plot.type = c("mean", "smooth", "boxplot"),
  region = NULL, ids = TRUE, ...)

Arguments

resids	An object of class brainGraph_resids (the output from get.resid)
corrs	List of lists of correlation matrices (as output by corr.matrix).
level	Character string; the level at which you want to calculate contributions (either global or regional)
control.value	Integer or character string specifying the control group (default: 1L)
object, x	A IC object
region	Character vector specifying which regions' IC's to print. Only relevant if method='Leave one out'
digits	Integer specifying the number of digits to display for P-values
...	Unused
plot.type	Character string indicating the type of plot; the default is to plot the mean (along with standard errors)
ids	Logical indicating whether to plot Study ID's for outliers. Otherwise plots the integer index

Value

A data.table with columns for

Study.ID	Subject identifier
Group	Group membership
region	If level='regional'
IC, RC	The value of the individual/regional contributions

Note

For aop, it is assumed by default that the control group is the first group.

Author(s)

Christopher G. Watson, [email protected]

References

Saggar, M. and Hosseini, S.M.H. and Buno, J.L. and Quintin, E. and Raman, M.M. and Kesler, S.R. and Reiss, A.L. (2015) Estimating individual contributions from group-based structural correlations networks. NeuroImage, 120, 274–284. doi:10.1016/j.neuroimage.2015.07.006

See Also

Other Structural covariance network functions: Bootstrapping, Residuals, brainGraph_permute, corr.matrix, import_scn, plot_volumetric

Examples

## Not run: 
IC <- loo(resids.all, corrs)
RC <- loo(resids.all, corrs, level='regional')

## End(Not run)
## Not run: 
IC <- aop(resids.all, corrs)
RC <- aop(resids.all, corrs, level='regional')

## End(Not run)

---

Calculate the inverse of the cross product of a design matrix

Description

inv is a S3 generic that calculates the inverse of the cross product of a design matrix, also referred to as the “unscaled covariance matrix”.

pinv calculates $M^{+} = (M^T M)^{-1} M^T$ for full (column) rank matrices. However, it does not verify the matrix's rank.

Usage

inv(x, ...)

## S3 method for class 'matrix'
inv(x, y = NULL, transpose = FALSE, ...)

## S3 method for class 'array'
inv(x, y = NULL, transpose = FALSE, ...)

## S3 method for class 'qr'
inv(x, p = x$rank, ...)

## S3 method for class 'list'
inv(x, p = x[[1L]]$rank, r = length(x),
  vnames = dimnames(x[[1L]]$qr)[[2L]], nms = names(x), ...)

pinv(x)

Arguments

x	A numeric matrix or array, a qr object, or a list of qr objects
...	Unused
y	A numeric matrix or vector (for the matrix and array methods). If supplied, this will be multiplied by x before the inverse is calculated. Default: NULL
transpose	Logical. If FALSE (the default), take the cross product of the arguments. If TRUE, use tcrossprod
p	The rank of the original matrix
r	The number of design matrices; i.e., the length of the input list
vnames	Character vector of the design matrix's variable names
nms	The region names; i.e., the names of the input list

Details

If x is a matrix, the Cholesky decomposition of the cross product is calculated (or using tcrossprod if transpose=TRUE), and the inverse is calculated from that result. That is,

$inv(X) = (X^T X)^{-1}$

$inv(X, transpose=TRUE) = (X X^T)^{-1}$

$inv(X, y) = (X^T y)^{-1}$

If x is a 3-dimensional array, then the inverse will be calculated for each matrix along the 3rd dimension, with the same input arguments for each.

Finally, there is a method for objects with class qr, and lists of QR decomposition objects.

Value

A numeric matrix or array

pinv returns the input matrix's pseudoinverse

Note

These methods should only be used on full-rank matrices, as there is no error checking being performed.

---

Calculate the AUC across densities of given attributes

Description

Given a list of brainGraphList objects, this function will calculate the area under the curve (AUC) across all thresholds/densities for each subject or group.

Usage

make_auc_brainGraph(g.list, g.attr = NULL, v.attr = NULL,
  norm = FALSE)

Arguments

g.list	A list of brainGraphList objects
g.attr	A character vector of graph attribute name(s). Default: NULL
v.attr	A character vector of vertex attribute name(s). Default: NULL
norm	Logical indicating whether to normalize threshold values to be between 0 and 1 (inclusive). Default: FALSE

Details

If the elements of the input list do not have a threshold element (or if it is NULL) then the AUC will be calculated on the interval $[0, 1]$ (inclusive). This has the same effect as specifying norm=TRUE in the function call.

Value

A brainGraphList object with one graph for each subject

Examples

## Not run: 
g.auc <- make_auc_brainGraph(g.fa, g.attr='E.global.wt')

## End(Not run)

---