Generate p-value histograms from AFNI statistical datasets

pval_hist.py

#!/usr/bin/env python3

"""
Usage:
    pval_hist.py DSET

Display histograms of the p-values contained in an AFNI statistical
dataset as output from 3dDeconvolve or 3dREMLfit.  DSET should contain
no sub-brick selectors.

If your contrast names are so long that AFNI has decided to truncate
them, you will have problems.
"""

from __future__ import print_function
# This is tested with python-3.4 but is intended to work with python-2.7.
# It won't work with python-2.6 or earlier.  Upgrade!

import sys

try:
    infile = sys.argv[1]
except IndexError:
    print(__doc__, file=sys.stderr)
    sys.exit(0)

import collections
import os.path
import re
import subprocess as sp
import tempfile

import matplotlib.pyplot as plt
import nibabel as nib
import numpy as np
import pandas as pd
import seaborn as sns

subbriks = (sp.check_output(['3dinfo', '-sb_delim', '\n', '-label', infile],
                            universal_newlines=True)
                .rstrip('\n').split('\n'))

labels = collections.OrderedDict()
for label in subbriks:
    match = re.match(r'(.*?)(#\d+)?_([^_]*)$', label)
    if match:
        stim, sub, stattype = match.groups('')
        # Betas aren't stats!
        if stattype == 'Coef':
            continue
        # Throw out duplicate tests with different statistics
        labels.setdefault(stim, collections.OrderedDict())[sub] = stattype

keeplabels = []
header = []
for stim, stimdict in labels.items():
    if '#1' not in stimdict:
        try:
            # Throw out t-test when corresponding F-test is equivalent
            del stimdict['#0']
        except KeyError:
            pass
    for sub, stattype in stimdict.items():
        keeplabels.append('{}{}_{}'.format(stim, sub, stattype))
        header.append((stim, stattype, sub))

# I have this set by default because I hate myself
os.environ['AFNI_COMPRESSOR'] = 'NONE'
with tempfile.TemporaryDirectory() as tempdir:
    tempname = os.path.join(tempdir, 'pval.nii')
    sp.check_call(
        ['3dPval', '-prefix', tempname,
         '{}[{}]'.format(infile, ','.join(keeplabels))])
    dset = nib.load(tempname).get_data()

mask = (dset < 1).any(axis=-1)
header = pd.Index(header, names=['model', 'test', 'subtest'])
pval_data = pd.DataFrame(dset[mask,:], columns=header)
del dset

ncol = np.ceil(np.sqrt(len(header.get_loc_level('', level=-1)[1])))
# Might try `sharey=False` instead of `ylim=(0, 5)`
(sns.FacetGrid(pval_data.unstack()
                   .reset_index(level=[0, 1, 2])
                   .rename_axis({0: 'pval'}, axis=1),
               ylim=(0, 5), col='model', col_wrap=ncol)
     .map(plt.hist, 'pval', bins=20, range=(0, 1), normed=True))

# From the Seaborn docs, because plt.hexbin doesn't play nice
def hexbin(x, y, color, **kwargs):
    cmap = sns.light_palette(color, as_cmap=True)
    plt.hexbin(x, y, gridsize=20, cmap=cmap, **kwargs)

for model, test in (header[~header.get_loc_level('', level=-1)[0]]
                        .droplevel(-1).unique()):
    # parameters might need tweaking for your display/backend
    (sns.PairGrid(pval_data.xs([model, test], level=[0, 1], axis=1))
         .map_offdiag(hexbin)
         .map_diag(plt.hist)
         .fig.suptitle('{}_{}'.format(model, test), y=0.965, fontsize=18))
    plt.subplots_adjust(top=0.9)

plt.show()