Spatial permutations for significance testing

This example shows how to perform spatial permutations tests (a.k.a spin-tests; Alexander-Bloch et al., 2018, NeuroImage) to assess whether two brain patterns are correlated above and beyond what would be expected from a spatially-autocorrelated null model.

While the original spin-tests were designed for comparing surface maps, we generally work with parcellated data in our lab. Using parcellations presents some novel difficulties in effectively implementing the spin-test, so this example demonstrates three spin-test methods.

For the original MATLAB toolbox published alongside the paper by Alexander-Bloch and colleagues refer to https://github.com/spin-test/spin-test.

An example dataset

The spin-test assumes that we have two sets of correlated data and that we are interested in assessing the degree to which this correlation exceeds a spatially-autocorrelated null model. We generate this null by permuting the original data by assuming it can be represented on the surface of the sphere and “spinning” the sphere.

First, let’s get some parcellated spatial maps that we can compare:

from netneurotools import datasets as nndata

data = nndata.fetch_vazquez_rodriguez2019(verbose=0)
rsq, grad = data['rsquared'], data['gradient']

The above function returns the \(R^{2}\) values of a structure-function linear regression model for each parcel (rsq) as well as the scores of each parcel along the first gradient computed from diffusion map embedding of a functional connectivity matrix (grad). (Refer to Vazquez-Rodriguez et al., 2019 for more information on these variables.)

These two vectors contain values for 1000 brain regions (a high-resolution sub-division of the Desikian-Killiany atlas; Cammoun et al., 2012). We’re interested in assessing the degree to which these two vectors are correlated; that is, how does the strength of the structure-function relationship in a brain region relate to its position along the first diffusion gradient?

from scipy.stats import pearsonr

r, p = pearsonr(rsq, grad)
print('r = {:.2f}, p = {:.4g}'.format(r, p))

The p-value suggests that our data are, indeed, highly correlated. Unfortunately, when doing this sort of correlation with brain data the null model used to generate the p-value does not take into account that the data are constrained by a spatial toplogy (i.e., the brain) and are therefore spatially auto-correlated. The p-values will be “inflated” because our true degrees of freedom are less than the number of samples we have.

To address this we can use a spatial permutation test (called a “spin test”), formally introduced by Alexander-Bloch et al., 2018, NeuroImage. This test works by considering the brain as a sphere and using random rotations of this sphere to construct a null distribution. If we rotate our data and resample datapoints based on their rotated values, we can generate a null that is more appropriate to our spatially auto-correlated data.

The original spin test was designed for working with vertex-level data; however, since we have parcellated data there are a few different options we have to choose between when performing a spin test.

Option 1: The “original” spin test

The original spin test assumes that you are working with two relatively high- resolution surface maps. It uses the coordinates of the vertices of these surfaces and applies random angular rotations, re-assigning values to the closest vertex (i.e., having the minimum Euclidean distance) after the rotation.

However, there are instances when two vertices may be assigned the same value because their closest rotated vertex is identical. When working with surfaces that are sampled at a sufficiently high resolution this will occur less frequently, but does still happen with some frequency. To demonstrate we can grab the coordinates of the fsaverage6 surface and generate a few rotations.

First, we’ll grab the spherical projections of the fsaverage6 surface and extract the vertex coordinates:

import nibabel as nib

# if you have FreeSurfer installed on your computer this will simply grab the
# relevant files from the $SUBJECTS_DIR directory; otherwise, it will download
# them to the $HOME/nnt-data/tpl-fsaverage directory
lhsphere, rhsphere = nndata.fetch_fsaverage('fsaverage6', verbose=0)['sphere']

lhvert, lhface = nib.freesurfer.read_geometry(lhsphere)
rhvert, rhface = nib.freesurfer.read_geometry(rhsphere)

Then, we’ll provide these to the function for generating the spin-based resampling array. We also need an indicator array designating which coordinates belong to which hemisphere so we’ll create that first:

from netneurotools import stats as nnstats
import numpy as np

coords = np.vstack([lhvert, rhvert])
hemi = [0] * len(lhvert) + [1] * len(rhvert)
spins, cost = nnstats.gen_spinsamples(coords, hemi, n_rotate=10, seed=1234)
print(spins.shape)
print(spins)

spins is an array that contains the indices that we can use to resample the fsaverage surface according to ten random rotations. The cost array indicates the total cost (in terms of Euclidean distance) of the re-assignments for each rotation.

The fsaverage surface has 81,924 vertices; let’s check how many are re-assigned for each rotation and what the average re-assignment distance is for each vertex:

for rotation in range(10):
    uniq = len(np.unique(spins[:, rotation]))
    print('Rotation {:>2}: {} vertices, {:.2f} mm / vertex'
          .format(rotation + 1, uniq, cost[rotation] / len(spins)))

In this case we can see that, for the first rotation, only 75,380 vertices were re-assigned (meaning that we “lost” data from 6,544 vertices), but these were assigned to vertices that were, on average, about 0.67 mm from the original. While this doesn’t seem too bad, when we lower the resolution of our data down even more (as we do with parcellations), this can become especially problematic.

We can demonstrate this for the 1000-node parcellation that we have for our dataset above. We need to define the spatial coordinates of the parcels on a spherical surface projection. To do this, we’ll fetch the left and right hemisphere FreeSurfer annotation files for the parcellation and then find the centroids of each parcel (defined on the spherical projection of the fsaverage surface):

from netneurotools import freesurfer as nnsurf

# this will download the Cammoun et al., 2012 FreeSurfer annotation files to
# the $HOME/nnt-data/atl-cammoun2012 directory
lhannot, rhannot = nndata.fetch_cammoun2012('surface', verbose=0)['scale500']

# this will find the center-of-mass of each parcel in the provided annotations
coords, hemi = nnsurf.find_fsaverage_centroids(lhannot, rhannot, surf='sphere')

The find_fsaverage_centroids() function return the xyz coordinates (coords) for each parcel defined in lhannot and rhannot, as well as an indicator array identifying to which hemisphere each parcel belongs (hemi):

We’ll use these coordinates to generate a resampling array as we did before for the fsaverage6 vertex coordinates:

# we'll generate 1000 rotations here instead of only 10 as we did previously
spins, cost = nnstats.gen_spinsamples(coords, hemi, n_rotate=1000, seed=1234)
for rotation in range(10):
    uniq = len(np.unique(spins[:, rotation]))
    print('Rotation {:>2}: {} parcels, {:.2f} mm / parcel'
          .format(rotation + 1, uniq, cost[rotation] / len(spins)))

We can see two things from this: (1) we’re getting more parcel duplications (i.e., only 727 out of the 1000 parcels were assigned in the first rotation), and (2) the distance from the original parcels to the re-assigned parcels has increased substantially from the fsaverage6 data.

This latter point makes sense: our parcellation provides a much sparser sampling of the cortical surface, so naturally parcels will be farther away from one another. However, this first issue of parcel re-assignment is a bit more problematic. At the vertex-level, when we’re densely sampling the surface, it may not be as much of a problem that some vertices are re-assigned multiple times. But our parcels are a bit more “independent” and losing up to 300 parcels for each rotation may not be desirable.

Nonetheless, we will use it to generate a spatial permutation-derived p-value for the correlation of our original data:

r, p = nnstats.permtest_pearsonr(rsq, grad, resamples=spins, seed=1234)
print('r = {:.2f}, p = {:.4g}'.format(r, p))

(Note that the maximum p-value from a permutation test is equal to 1 / (n_perm + 1).)

The benefit of generating our resampling array independent of a statistical test is that we can re-use it for any number of applications. If we wanted to conduct a Spearman correlation instead of a Pearson correlation we could easily do that:

from scipy.stats import rankdata

rho, prho = nnstats.permtest_pearsonr(rankdata(rsq), rankdata(grad),
                                      resamples=spins, seed=1234)
print('rho = {:.2f}, p = {:.4g}'.format(rho, prho))

Option 2: Exact matching

#