Working with Adjacency¶

Adjacency holds square matrices over a set of nodes: similarity and distance matrices, functional connectivity, directed graphs. A symmetric matrix is stored as its strict upper triangle and rebuilt on demand, and one object can hold a stack of matrices — one per subject, region or timepoint. Most of its methods mirror BrainData's. Representational similarity analysis is the archetypal use of these matrices and has its own tutorial.

Create one¶

The constructor takes a square matrix, a polars or pandas DataFrame, or a .csv or .h5 path, and a list of those stacks them — except .h5 paths, which load one at a time. matrix_type declares what the values mean: 'similarity' (higher is more alike), 'distance' (higher is further apart), or 'directed' (asymmetric, stored in full). Append '_flat' to pass values that are already vectorized; a bare vector with no matrix_type is read as a distance. labels names the nodes and travels with the object. Distance and similarity matrices store the strict upper triangle, so an input diagonal is dropped, and a matrix declared symmetric has to be symmetric — hence the symmetrized noise below. The example has three groups of four nodes, connected within group at strengths 1, 2 and 3 and unconnected across groups:

import numpy as np
from scipy.linalg import block_diag

from nltools.data import Adjacency

def symmetric_noise(rng, scale):
    noise = rng.standard_normal((12, 12)) * scale
    return (noise + noise.T) / 2

m1 = block_diag(np.ones((4, 4)), np.zeros((4, 4)), np.zeros((4, 4)))
m2 = block_diag(np.zeros((4, 4)), np.ones((4, 4)), np.zeros((4, 4)))
m3 = block_diag(np.zeros((4, 4)), np.zeros((4, 4)), np.ones((4, 4)))

blocks = Adjacency(
    m1 + 2 * m2 + 3 * m3 + symmetric_noise(np.random.default_rng(0), 0.1),
    matrix_type="similarity",
    labels=[f"{group}.{node}" for group in ("C1", "C2", "C3") for node in range(4)],
)

print(blocks)

nltools.data.adjacency.Adjacency(shape=(12, 12), Y=(0, 0), is_symmetric=True, matrix_type=similarity)

Printing reports the logical shape, the size of the metadata frame Y, symmetry, and the matrix type. .shape is that logical shape — (n_matrices, n_nodes, n_nodes) for a stack — while .vector_shape is what is stored, here the n(n-1)/2 entries above the diagonal. squareform() rebuilds the square matrix as a plain numpy array with a zero diagonal, or a list of them for a stack.

BrainData.distance() returns an Adjacency too: every image against every other under a scipy metric ('euclidean' by default). Twenty-one images from the pain dataset — seven subjects at low, medium and high intensity — give a 21-node distance matrix, and joblib keeps the build from redownloading them:

print(f"shape:        {blocks.shape}  ({blocks.n_nodes} nodes)")
print(f"vector_shape: {blocks.vector_shape}")
print(f"squareform(): {blocks.squareform().shape} array")

shape:        (12, 12)  (12 nodes)
vector_shape: (66,)
squareform(): (12, 12) array

from joblib import Memory

from nltools.datasets import fetch_pain

memory = Memory(".tutorial-cache", verbose=0)
subset = memory.cache(fetch_pain)()[:21]
pain_distance = subset.distance(metric="correlation")
pain_distance.labels = [f"s{subject}" for subject in subset.X["SubjectID"]]

print(pain_distance)

nltools.data.adjacency.Adjacency(shape=(21, 21), Y=(0, 0), is_symmetric=True, matrix_type=distance)

Plot¶

plot() draws a seaborn heatmap of the square form, using labels for the ticks when they are set. Keywords the signature does not name go to seaborn.heatmap, so cmap, vmin and vmax work. The three planted groups are the brighter diagonal squares below. The pain matrix holds distances, so its blocks are dark instead: each image is close to the same subject's other two images and far from everyone else's.

blocks.plot()

pain_distance.plot()

Threshold and transform¶

threshold zeroes part of the range and keeps the rest. upper= keeps values at or above the cutoff, lower= keeps values at or below it, and giving both keeps the two tails outside the band. Read each as the edge of the region you keep, not the region you drop: the interesting edges of a similarity matrix are the high ones, so reach for upper=, and of a distance matrix the low ones, so reach for lower=. A string ending in % is a percentile of the stored values; binarize=True turns whatever survives into ones.

strong = blocks.threshold(upper=0.5)
top_decile = blocks.threshold(upper="90%")
binary = blocks.threshold(upper=0.5, binarize=True)

print(f"upper=0.5    keeps {(strong.data != 0).sum()} of {blocks.vector_shape[0]}")
print(f"upper='90%'  keeps {(top_decile.data != 0).sum()}")
print(f"binarized values: {np.unique(binary.data)}")

upper=0.5    keeps 18 of 66
upper='90%'  keeps 7
binarized values: [0. 1.]

distance_to_similarity converts a distance matrix into a similarity one: metric='correlation' returns 1 - d, undoing the correlation distance above, and metric='euclidean' returns exp(-beta * d / sd), where sd is the standard deviation of the square matrix. There is no method for the other direction — a correlation similarity goes back as 1 - s, with matrix_type set on the constructor. Correlations are also bounded at ±1 and their sampling variance shrinks near those ends, so averages and parametric tests misbehave there. Fisher's r-to-z is the usual fix: r_to_z() is arctanh and z_to_r() is tanh, applied elementwise to the stored triangle:

pain_similarity = pain_distance.distance_to_similarity(metric="correlation")
pain_z = pain_similarity.r_to_z()

_dist, _sim, _z = pain_distance.data, pain_similarity.data, pain_z.data

print(f"distance   [{_dist.min():.2f}, {_dist.max():.2f}]")
print(f"similarity [{_sim.min():.2f}, {_sim.max():.2f}]")
print(f"z          [{_z.min():.2f}, {_z.max():.2f}]")
print(f"round trip: {np.allclose(pain_z.z_to_r().data, _sim)}")

distance   [0.08, 1.42]
similarity [-0.42, 0.92]
z          [-0.45, 1.58]
round trip: True

Arithmetic and statistics¶

+, -, * and / work elementwise against another Adjacency or a scalar. Two matrices have to agree on node count, matrix type, and node labels and their order — nothing else would catch nodes stacked in different orders. mean, median, std and sum collapse the object: on a single matrix they return one number over the stored edges; on a stack the default axis=0 averages across matrices into an Adjacency, and axis=1 returns one number per matrix. Subtracting the first group's mask leaves the other two groups and the noise:

residual = blocks - Adjacency(m1, matrix_type="similarity", labels=blocks.labels)

print(f"blocks:              mean {blocks.mean():.2f}, sd {blocks.std():.2f}")
print(f"first group removed: mean {residual.mean():.2f}, sd {residual.std():.2f}")

blocks:              mean 0.55, sd 1.00
first group removed: mean 0.46, sd 1.00

A list of matrices makes a stack, and append adds to an existing one. Fifteen matrices stand in for fifteen subjects: the first group is connected in five of them, absent in the next five and connected again in the last five, under noise five times larger than above. Averaging across the stack recovers the planted group at two thirds of its strength, since a third of the matrices do not carry it. ttest then tests every edge against zero, returning 'mean', 't', 'z' and 'p' as four Adjacency maps; permutation=True swaps the parametric p-value for a sign-flip one and leaves 't' parametric. The maps are unthresholded: 66 edges is 66 tests, so correct for them.

stack_rng = np.random.default_rng(1)
stack = Adjacency(
    [m1 + symmetric_noise(stack_rng, 0.5) for _ in range(5)]
    + [symmetric_noise(stack_rng, 0.5) for _ in range(5)]
    + [m1 + symmetric_noise(stack_rng, 0.5) for _ in range(5)],
    matrix_type="similarity",
    labels=blocks.labels,
)

print(stack)

nltools.data.adjacency.Adjacency(shape=(15, 12, 12), Y=(0, 0), is_symmetric=True, matrix_type=similarity)

stack.mean().plot()

group = stack.ttest()
_p = group["p"].data

print(f"{(_p < 0.05).sum()} of {len(_p)} edges at p < .05")
print(f"the planted group is {int(m1[np.triu_indices(12, 1)].sum())} edges")

8 of 66 edges at p < .05
the planted group is 6 edges

Regression¶

regress answers two different questions and the design decides which. Both return a dict of beta, sigma, t, p, df and residual, and neither adds an intercept for you. tail is keyword-only and takes 2 (two-tailed, the default) or 1.

Pass an Adjacency and this matrix is decomposed into a weighted sum of the predictor matrices: edges are the observations, matrices the predictors, and the coefficients come back as plain arrays. Only a single response matrix works, and predictor and response must agree on node ordering, so the design carries the same labels. The three group masks recover the strengths the data was built with. Pass a DesignMatrix and each edge is regressed across the stack instead — the adjacency analogue of a mass-univariate imaging analysis, with the same multiple-comparisons problem. Each coefficient is then an Adjacency holding one matrix per predictor, and thresholding the t map leaves the on-off-on group:

design = Adjacency([m1, m2, m3], matrix_type="similarity", labels=blocks.labels)
block_fit = blocks.regress(design)

print(f"beta: {block_fit['beta'].round(2)}")
print(f"t:    {block_fit['t'].round(1)}")
print(f"df:   {block_fit['df']}")

beta: [0.97 2.05 3.02]
t:    [ 33.   69.6 102.5]
df:   63

from nltools.data import DesignMatrix

on_off_on = DesignMatrix(
    np.array([1] * 5 + [0] * 5 + [1] * 5).reshape(-1, 1),
    columns=["on"],
    sampling_freq=1.0,
)
edge_fit = stack.regress(on_off_on)

edge_fit["t"].threshold(upper=2).plot()

Multidimensional scaling¶

plot_mds lays a single distance matrix out in two or three dimensions (n_components); metric scaling is the default and metric_mds=False asks for non-metric. Node names come from labels, and labels_color takes one color per node, here each image's intensity.

The images land in seven tight clusters, one per subject; the manipulated variable, intensity, does not organize the picture at all. The distances agree: two images of one subject are more than three times closer than two of different subjects, and sharing an intensity across subjects buys nothing:

intensity_color = {"low": "#4c72b0", "medium": "#dd8452", "high": "#c44e52"}

pain_distance.plot_mds(
    labels_color=[intensity_color[level] for level in subset.X["PainIntensity"]]
)

_rows, _cols = np.triu_indices(pain_distance.n_nodes, 1)
_edges = pain_distance.squareform()[_rows, _cols]
_subject = subset.X["SubjectID"].to_numpy()
_level = subset.X["PainIntensity"].to_numpy()
_across = _subject[_rows] != _subject[_cols]
_same_level = _level[_rows] == _level[_cols]

print(f"same subject:         {_edges[~_across].mean():.2f}")
print(f"different subject:    {_edges[_across].mean():.2f}")
print(f"  and same intensity: {_edges[_across & _same_level].mean():.2f}")

same subject:         0.25
different subject:    0.87
  and same intensity: 0.85

Graphs¶

to_graph hands a single matrix to networkx — a Graph, or a DiGraph for a directed matrix — with labels as the node names, so give the nodes distinct labels or two of them will merge into one. Zero entries become non-edges, so the binarized matrix from earlier arrives as three disconnected cliques of four nodes and every node has degree 3. Every networkx metric and layout applies from there. A node-level measure is one number per node, and when the nodes are brain regions roi_to_brain writes those numbers back into the parcellation as an image; the BrainData tutorial runs that round trip end to end.

import networkx as nx

clique_graph = binary.to_graph()
_degrees = sorted(set(dict(clique_graph.degree()).values()))

print(f"{len(clique_graph)} nodes, {clique_graph.number_of_edges()} edges")
print(f"degrees: {_degrees}")

12 nodes, 18 edges
degrees: [3]

nx.draw_circular(clique_graph, node_color="lightsteelblue", with_labels=True)