Skip to content

Instantly share code, notes, and snippets.

@effigies
Created July 24, 2018 14:25
  • Star 1 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save effigies/0267812074c8eb49492ac0584d6730c1 to your computer and use it in GitHub Desktop.
Markiewicz and Bohland, 2016 permutation testing

Partially-parametric permutation testing

Cluster thresholding in multi-voxel pattern analysis (MVPA) is an open problem, as many figures of merit are possible, few (if any) of which have been sufficiently analyzed to permit a parametric solution or a guarantee of compatibility with pre-computed simulations. Stelzer, et al. 2012 represents probably the most conservative approach, constructing voxel-wise and then cluster-wise null distributions at the group level, based on permuting the training labels at the individual level.

In Mapping the cortical representation of speech sounds in a syllable repetition task, we adapted this approach to skip the voxel-wise null distribution, instead relying on a parametric approach to cluster definition derived from Lee, et al., 2012. In brief, the approach takes the spatial mean within subject, and performs a one-sided t-test across subjects to find regions where the activation is consistently above 0 (i.e., the spatial mean).

I have recently recovered the code that performed this analysis, but it is currently in a repository context that is not fully abstracted from the details of the system it had been running on. This Gist is a curation of the code used to perform this statistical analysis. It is not meant to be complete, but to help fill in the pieces that are not specific to my studies.

This code is presented as-is. I have not attempted to clean it up, and inconsistencies in style and poor names are the detritus of 5 years of iterative work on the same few problems.

The code is in Python. We will be using these libraries (and one utility function):

import os
import numpy as np
import nibabel as nib
import scipy as sp
import scipy.stats
import mne
import networkx as nx
from collections import Counter

def prep_dirs(targets):
    """Create all necessary directories in a path"""
    if isinstance(targets, basestring):
        targets = [targets]

    for dirname in map(os.path.dirname, targets):
        if not os.path.isdir(dirname):
            try:
                os.makedirs(dirname)
            except OSError as e:
                # Already exists errors are common when running parallel
                # jobs
                if e[0] != 17:
                    raise e

The following functions are useful for performing some simple operations on nibabel images:

def image_like(img, data):
    return img.__class__(data, img.affine, img.header,
                         extra=img.extra.copy())

def mapImage(func, img):
    return image_like(img, func(img.get_data()))

def tTestImages(images, popmean=0.0):
    """Perform per-entry t-test on nibabel spatial images"""
    stack = nib.concat_images(images, check_affines=False)

    tstats, pvalues = sp.stats.ttest_1samp(stack.get_data(), popmean, axis=3)
    reject, corrected = mne.stats.fdr_correction(pvalues)

    return (image_like(stack, tstats), image_like(stack, pvalues),
            image_like(stack, corrected))

def leeTTestMGHs(overlays):
    """Perform group-level t-tests with each session's spatial mean removed.

    Alternative to non-parametric group statistics from Lee, et al. 2012"""
    demean = lambda x: mapImage(lambda y: y - np.mean(y), x)
    return tTestImages(map(demean, map(nib.load, overlays)), 0.0)

These are functions that I used to work with FreeSurfer directories/surface meshes:

def subjFromSess(sess, sessions_dir=os.environ['FUNCTIONALS_DIR']):
    with open(os.path.join(sessions_dir, sess, 'subjectname'), 'r') as f:
        return f.readline().rstrip('\n')

def getSurface(subject, hemi, surface='white',
               subject_dir=os.environ['SUBJECTS_DIR']):
    surf_path = os.path.join(subject_dir, subject, 'surf',
                             hemi + '.' + surface)
    return nib.freesurfer.read_geometry(surf_path)

def buildSurfaceGraph(subject, hemi, surface='white',
                      subject_dir=os.environ['SUBJECTS_DIR']):
    coords, faces = getSurface(subject, hemi, surface, subject_dir)
    graph = nx.Graph()

    for i, j, k in faces:
        graph.add_edge(i, j)
        graph.add_edge(i, k)
        graph.add_edge(j, k)

    return graph
  
def getLabelSurf(sess, hemi, parcellation, results_dir=params.results_dir):
    return nib.freesurfer.read_annot(os.path.join(results_dir, sess, 'label',
                                     '{}.{}.annot'.format(hemi, parcellation)))
  
def filterLabelsFromGraph(graph, label, session, hemi, parcellation,
                          results_dir=params.results_dir):
    labels, ctab, names = getLabelSurf(session, hemi, parcellation,
                                       results_dir)
    indicator = (np.asarray(names)[labels] != label) & (labels > 0)
    return graph.subgraph(np.nonzero(indicator)[0])

These last two functions assume that a non-standard FreeSurfer label file is being used, and it can be found in a results directory that mirrors the FreeSurfer directory structure. I believe that for the purposes used here, it is only used to remove the medial wall, and so could use pretty much any label file in FreeSurfer.

In the following, I assume the existence of the function:

def getPermutationInputs(session, hemi, attribute, mapid, ds_path,
                         balanced, permutations):
    # session = FreeSurfer session
    # hemi = hemisphere
    # attribute = distingusihing name for statistic
    # mapid = permutation number
    # ds_path = base path for dataset results
    # balanced = number of times training set was rebalanced
    # permutations = total number of permutations

This should be replaced with some way of getting a random set of individual level chance maps. For computational purposes, I would recommend finding some way of returning the same set of inputs for the same set of parameters. I solved this by simply saving the set in a text file when first generated, and only regenerating if the text file was deleted.

The following function calculates cluster sizes from a single group permutation. It follows a similar protocol to the above, with some defaults pre-filled (graphs is a dictionary that may contain pre-built left/right hemisphere meshes).

def getLeeChanceClusters(dataset, hemi, attribute, mapid, vthresh=0.05,
                         balanced=2, permutations=100, graphs=GRAPHS,
                         session='SRavgstc',
                         results_dir=params.results_dir):
    timestamp('getLeeChanceClusters({}, {}, {}, {:d}, {:d}, {:d}, {})'.format(
              dataset, hemi, attribute, mapid, balanced, permutations,
              session))
    ds_path = os.path.join(results_dir, session, dataset)

    inputs = getPermutationInputs(session, hemi, attribute, mapid, ds_path,
                                  balanced, permutations)

    if not all(map(os.path.exists, inputs)):
        print 'Cannot find all of: ' + str(inputs)
        return []

    tstats, pvals, fdr = leeTTestMGHs(inputs)
    if hemi not in graphs.setdefault(session, {}):
        base = buildSurfaceGraph(subjFromSess(session), hemi)
        graphs[session][hemi] = filterLabelsFromGraph(base, 'None', session,
                                                      hemi, 'SLaparc17')

    # One-sided threshold from two-sided test
    thresholded = (pvals.get_data() < 2 * vthresh) & (tstats.get_data() > 0)
    subgraphs = nx.algorithms.connected_component_subgraphs(
        graphs[session][hemi].subgraph(np.nonzero(thresholded)[0]))

    return map(len, subgraphs)

The following gets you a cumulative distribution function from a list (hist) of sizes:

def cdfFromHist(hist):
    # Count instances of each cluster size
    c = Counter(hist)

    # Cumulative sum gives clusters of a given size or lower
    accumulated = np.cumsum([count for size, count in sorted(c.items())])
    # print (accumulated[-1], len(hist))
    assert accumulated[-1] == len(hist)

    # Map cluster sizes to proportion of smaller clusters in histogram
    propSmaller = dict(zip((size + 1 for size in sorted(c.keys())),
                           accumulated / float(accumulated[-1])))

    # Cumulative distribution function CDF(x) = P(X < x)
    # Using strict inequality, there are no clusters smaller than 1, so
    # cdf[0] = cdf[1] = 0
    #
    # By using proportion smaller, strict probabilities are maintained
    # Additionally, ensure that the region between jumps in the CDF is
    # equal to the start of the step
    last = None
    cdf = np.zeros(max(propSmaller.keys()))
    for i in sorted(propSmaller.keys()):
        if last is not None:
            cdf[last:i] = propSmaller[last]
        last = i

    return cdf

With these, we can now construct a CDF of chance clusters for any number of group-level chance maps we want. Again, the directory structure is leaking into the code.

def getLeeCDF(dataset, attribute, nmaps, vthresh=0.05, session='SRavgstc',
              results_dir=params.results_dir):
    cdf_path = os.path.join(results_dir, session, dataset,
                            'averageAccuracyMaps', attribute, 'cluster_sizes',
                            'lee{:06d}.{:g}.npy'.format(nmaps, 100*vthresh))

    if os.path.exists(cdf_path):
        return np.load(cdf_path)

    lengths = []
    for i in range(nmaps):
        for hemi in ('lh', 'rh'):
            lengths.extend(getLeeChanceClusters(dataset, hemi, attribute, i,
                                                vthresh=vthresh,
                                                session=session,
                                                results_dir=results_dir))

    cdf = cdfFromHist(lengths)

    prep_dirs(cdf_path)
    np.save(cdf_path, cdf)

    return cdf

We can now find the percentile ranks of clusters from the empirical group map, using the CDF of the null distribution defined by the above functions. The following is somewhat study-specific, in that the directory structure and file names are fairly directly coded in. This will need some abstraction from these details to be usable.

def leeClusters(vpercentile, dataset, attribute, nmaps=10000, graphs=GRAPHS,
                session='SRavgstc', results_dir=params.results_dir):
    vthresh = 1 - vpercentile
    cdf = getLeeCDF(dataset, attribute, nmaps, vthresh=vthresh,
                    session=session, results_dir=results_dir)
    ds_path = os.path.join(results_dir, session, dataset)

    # ``analysis`` needs to be defined

    sessions = subsessions(session, results_dir=params.results_dir)

    in_fmt = os.path.join(ds_path, '{}',
                          '{{}}.{}_{}.mgh'.format(attribute, analysis)).format
    out_fmt = os.path.join(results_dir, session, 'overlays',
                           '{{}}.{}_{}_{}.lee{:g}.mgh'.format(
                               dataset, attribute, analysis,
                               vpercentile)).format

    for hemi in ('lh', 'rh'):
        tstats, pvals, fdr = leeTTestMGHs(in_fmt(sess, hemi)
                                          for sess in sessions)
        ranks = np.zeros(tstats.shape, '>f4')
        if hemi not in graphs.setdefault(session, {}):
            base = buildSurfaceGraph(subjFromSess(session), hemi)
            graphs[session][hemi] = filterLabelsFromGraph(base, 'None',
                                                          session, hemi,
                                                          'SLaparc17')

        # pvals are two-tailed, so take positive only and double threshold
        # for one-tailed
        superthreshold = np.nonzero((pvals.get_data() < 2 * vthresh) &
                                    (tstats.get_data() > 0))[0]

        subgraphs = nx.algorithms.connected_component_subgraphs(
            graphs[session][hemi].subgraph(superthreshold))

        for sub in subgraphs:
            ranks[sub.nodes()] = cdf[len(sub)] if len(sub) < len(cdf) else 1

        image_like(tstats, ranks).to_filename(out_fmt(hemi))

Note that these percentile ranks are 1 - p for uncorrected p-values. They will need to be Bonferroni corrected for the number of analyses that you're running.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment