One of the reviewers of our paper An antidote to the imager's fallacy, or how to identify brain areas that are in limbo suggested a simpler approach for the in-limbo test: simply subtract the contrast size of the comparison voxel and do the level 2 analysis again.
This test is implemented in this notebook.
First we import import libraries
import nipype.algorithms.modelgen as model
import pandas
import nipype.pipeline.engine as pe
from nipype.interfaces.base import Bunch
import nipype.interfaces.utility as util
import nipype.interfaces.io as nio
import nipype.interfaces.fsl as fsl
We find the copes, varcopes and dof_files of the first level FSL analysis
import glob, os
copes = sorted(glob.glob('/home/gdholla1/data/in_limbo/birte/copes_mni/*/cope1.nii.gz'))
varcopes = sorted(glob.glob('/home/gdholla1/data/in_limbo/birte/varcopes_mni/*/varcope1.nii.gz'))
dof_files = [os.path.abspath('dof')] * len(copes)
Set up different thresholds
thresholds = [2.3, 3.1, 3.7]
This is a completely standard level 2 analysis pipeline with FSL FLAME. Because the SPMs of the first level analysis were unsmoothed, we have to do that first.
from nipype.workflows.fmri.fsl import create_fixed_effects_flow
# Create mixed effects work flow
mixedfx = create_fixed_effects_flow('in_limbo_level2_birte_fsl')
mixedfx.base_dir = '/home/gdholla1/workflow_folders/'
# Set up smoothing nodes
smoother_copes = pe.MapNode(fsl.IsotropicSmooth(),
iterfield=['in_file'],
name='smoother_copes')
smoother_copes.inputs.in_file = copes
smoother_copes.inputs.fwhm = 8.0
smoother_varcopes = pe.MapNode(fsl.IsotropicSmooth(),
iterfield=['in_file'],
name='smoother_varcopes')
smoother_varcopes.inputs.in_file = varcopes
smoother_varcopes.inputs.fwhm = 8.0
mixedfx.connect(smoother_copes, 'out_file', mixedfx.get_node('inputspec'), 'copes')
mixedfx.connect(smoother_varcopes, 'out_file', mixedfx.get_node('inputspec'), 'varcopes')
mixedfx.inputs.inputspec.dof_files = dof_files
mixedfx.inputs.flameo.mask_file = '/usr/share/fsl/data/standard/MNI152_T1_2mm_brain_mask.nii.gz'
mixedfx.inputs.l2model.num_copes = len(copes)
### Mixed effects (flame1) instead of fixed effects (fe)
mixedfx.inputs.flameo.run_mode = 'flame1'
# Set up a datasink to write data to
ds = pe.Node(nio.DataSink(), name='datasink')
ds.inputs.base_directory = '/home/gdholla1/data/in_limbo/birte/level2_fsl'
mixedfx.connect(mixedfx.get_node('outputspec'), 'zstats', ds, 'zstats')
# Estimate smoothness for Gaussian Random Field Theory
smooth_est = pe.Node(fsl.SmoothEstimate(), name='smooth_est')
smooth_est.inputs.mask_file = '/usr/share/fsl/data/standard/MNI152_T1_2mm_brain_mask.nii.gz'
get_one = lambda x: x[0]
mixedfx.connect(mixedfx.get_node('outputspec'), ('zstats', get_one), smooth_est, 'zstat_file')
cluster = pe.Node(fsl.Cluster(), name='cluster')
cluster.inputs.out_threshold_file = True
# Iterate over different thresholds
cluster.iterables = [('threshold',
thresholds)]
# Set up GRF and write threhsolded z-map to datasink
mixedfx.connect(smooth_est, 'dlh', cluster, 'dlh')
mixedfx.connect(mixedfx.get_node('outputspec'), ('zstats', get_one), cluster, 'in_file')
mixedfx.connect(cluster, 'threshold_file', ds, 'thresholded_z')
We set up a function argmin that gets the index of the least-significantly activated voxel in the level 2 result.
def argmin(image, threshold=0):
import nibabel as nb
import numpy as np
image = nb.load(image)
data = np.ma.masked_less_equal(image.get_data(), threshold)
return np.unravel_index(np.ma.argmin(data), data.shape)
get_comparison_voxel = pe.Node(util.Function(function=argmin,
input_names=['image'],
output_names=['index']),
name='get_comparison_voxel')
We set up a function get_min_value that gets the contrast sizes in this least-significantly activated voxel for all subjects
get_min_value = pe.MapNode(util.Function(function=get_value_at,
input_names=['image', 'index'],
output_names=['value']),
iterfield=['image'],
name='get_min_value')
def get_value_at(image, index):
import nibabel as nb
return nb.load(image).get_data()[index]
Now we connect the standard level 2 analysis to these nodes
mixedfx.connect(cluster, 'threshold_file', get_comparison_voxel, 'image')
mixedfx.connect(get_comparison_voxel, 'index', get_min_value, 'index')
mixedfx.connect(smoother_copes, 'out_file', get_min_value, 'image')
We subtract the contrasts sizes in the comparison voxel from the individual contrast maps. Then we multiply everything with -1.
def subtract_and_inverse(image, value):
import nibabel as nb
import os
image = nb.load(image)
data = image.get_data()
data -= value
data *= -1
fn = os.path.abspath('cope_in_limbo.nii.gz')
nb.save(nb.Nifti1Image(data, image.get_affine()), fn)
return fn
subtract_cope = pe.MapNode(util.Function(function=subtract_and_inverse,
input_names=['image', 'value'],
output_names=['in_limbo_cope'],
),
iterfield=['image',
'value'],
name='subtract_from_cope')
mixedfx.connect(get_min_value, 'value', subtract_cope, 'value')
mixedfx.connect(smoother_copes, 'out_file', subtract_cope, 'image')
Now we create a second standard level 2 analysis pipeline. We input the 'subtracted-and-inversed' copes from the first level 2 analysis.
in_limbo_mixedfx = create_fixed_effects_flow('in_limbo_mixedfx')
mixedfx.connect(smoother_varcopes, 'out_file', in_limbo_mixedfx, 'inputspec.varcopes')
in_limbo_mixedfx.inputs.inputspec.dof_files = dof_files
in_limbo_mixedfx.inputs.flameo.mask_file = '/usr/share/fsl/data/standard/MNI152_T1_2mm_brain_mask.nii.gz'
in_limbo_mixedfx.inputs.l2model.num_copes = len(copes)
mixedfx.connect(subtract_cope, 'in_limbo_cope', in_limbo_mixedfx, 'inputspec.copes')
Now we have to sort all voxels according to the outcome of the two pipelines. We have three categories and apply the following logic:
def get_in_limbo_area(thresholded_z, in_limbo_t, threshold):
import nibabel as nb
import numpy as np
import os
# Load the two level 2 results
z_image = nb.load(thresholded_z)
z = z_image.get_data()
in_limbo_t = nb.load(in_limbo_t).get_data()
# Set up outputs
in_limbo = np.zeros_like(z)
unactivated = np.zeros_like(z)
# Apply logic
in_limbo[(z == 0) & (in_limbo_t < threshold) & (in_limbo_t != 0)] = 1
unactivated[(z == 0) & (in_limbo_t > threshold)] = 1
nb.save(nb.Nifti1Image(in_limbo, z_image.get_affine(),),
os.path.abspath('in_limbo.nii.gz'))
nb.save(nb.Nifti1Image(unactivated, z_image.get_affine(),),
os.path.abspath('in_limbo_unactivated.nii.gz'))
return os.path.abspath('in_limbo.nii.gz'), os.path.abspath('in_limbo_unactivated.nii.gz')
in_limbo_mapper = pe.Node(util.Function(function=get_in_limbo_area,
input_names=['thresholded_z',
'in_limbo_t',
'threshold'],
output_names=['in_limbo', 'unactivated']),
name='in_limbo_mapper')
We have to set an $\alpha$-value for the in-limbo test and calculate the corresponding t-value. Here we do a one-sided alpha level of .05 and calculate the corresponding t-threshold. We put that into the get_in_limbo_area function.
import scipy as sp
from scipy import stats
dof = len(copes) - 1
in_limbo_mapper.inputs.threshold = sp.stats.t(dof).ppf(0.95)
mixedfx.connect(in_limbo_mixedfx.get_node('outputspec'), ('tstats', get_one), in_limbo_mapper, 'in_limbo_t')
mixedfx.connect(cluster, 'threshold_file', in_limbo_mapper, 'thresholded_z')
Now we can make a graphical representation of the workflow
from IPython.display import Image
mixedfx.write_graph()
Image('/home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/graph.dot.png')
INFO:workflow:Converting dotfile: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/graph.dot to png format
And run it!
mixedfx.run()
INFO:workflow:['check', 'execution', 'logging'] INFO:workflow:Running serially. INFO:workflow:Executing node l2model in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/l2model INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node smoother_copes in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/smoother_copes INFO:workflow:Executing node l2model in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/l2model INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node copemerge in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/copemerge INFO:workflow:Executing node gendofvolume in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/gendofvolume INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node smoother_varcopes in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/smoother_varcopes INFO:workflow:Executing node varcopemerge in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/varcopemerge INFO:workflow:Executing node flameo in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/flameo INFO:workflow:Executing node smooth_est in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/smooth_est INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node cluster.aI.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_2.3/cluster INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node get_comparison_voxel.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_2.3/get_comparison_voxel INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node get_min_value.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_2.3/get_min_value INFO:workflow:Executing node cluster.aI.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.7/cluster INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node get_comparison_voxel.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.7/get_comparison_voxel INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node get_min_value.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.7/get_min_value INFO:workflow:Executing node datasink.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_2.3/datasink INFO:workflow:Executing node datasink.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.7/datasink INFO:workflow:Executing node cluster.aI.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.1/cluster INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node datasink.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.1/datasink INFO:workflow:Executing node get_comparison_voxel.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.1/get_comparison_voxel INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node subtract_from_cope.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_2.3/subtract_from_cope INFO:workflow:Executing node copemerge.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_2.3/copemerge INFO:workflow:Executing node subtract_from_cope.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.7/subtract_from_cope INFO:workflow:Executing node copemerge.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_3.7/copemerge INFO:workflow:Executing node gendofvolume.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_3.7/gendofvolume INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node gendofvolume.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_2.3/gendofvolume INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node get_min_value.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.1/get_min_value INFO:workflow:Executing node subtract_from_cope.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.1/subtract_from_cope INFO:workflow:Executing node copemerge.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_3.1/copemerge INFO:workflow:Executing node gendofvolume.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_3.1/gendofvolume INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node varcopemerge in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/varcopemerge INFO:workflow:Executing node flameo.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_2.3/flameo INFO:workflow:Executing node in_limbo_mapper.a0 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_2.3/in_limbo_mapper INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node flameo.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_3.7/flameo INFO:workflow:Executing node in_limbo_mapper.a2 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.7/in_limbo_mapper INFO:workflow:Collecting precomputed outputs INFO:workflow:Executing node flameo.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/in_limbo_mixedfx/_threshold_3.1/flameo INFO:workflow:Executing node in_limbo_mapper.a1 in dir: /home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_3.1/in_limbo_mapper INFO:workflow:Collecting precomputed outputs
<networkx.classes.digraph.DiGraph at 0x7e7b090>
Load Seaborn for nice plots
import seaborn as sns
palette = sns.blend_palette(["seagreen", "lightblue"]);
sns.set_style('white')
Plot for the thresholds
# Don't put warnings in my plots, please, Python
import warnings
warnings.filterwarnings('ignore')
x_slice = 41
y_slice = 49
z_slice = 31
for threshold in [2.3, 3.1, 3.7]:
## STANDARD FSL AR(1)-model
in_limbo = nb.load('/home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_%s/in_limbo_mapper/in_limbo.nii.gz' % threshold).get_data()
z_map = nb.load('/home/gdholla1/workflow_folders/in_limbo_level2_birte_fsl/_threshold_%s/cluster/zstat1_threshold.nii.gz' % threshold).get_data()
brain = nb.load('/usr/share/fsl/data/standard/MNI152_T1_2mm_brain.nii.gz').get_data()
z_map = np.ma.masked_less_equal(z_map, 0)
in_limbo = np.ma.masked_less_equal(in_limbo, 0)
plt.figure(figsize=(20, 5))
plt.suptitle('AR(1) estimation - z > %s' % threshold, fontsize=25)
plt.subplot(131)
plt.imshow(brain[x_slice, :, :].T, origin='lower', cmap=plt.cm.gray)
plt.imshow(z_map[x_slice, :, :].T, origin='lower', cmap=plt.cm.hot)
plt.imshow(in_limbo[x_slice, :, :].T, origin='lower', cmap=plt.cm.Greens_r)
plt.axis('off')
plt.subplot(132)
plt.imshow(brain[:, y_slice, :].T, origin='lower', cmap=plt.cm.gray)
plt.imshow(z_map[:, y_slice, :].T, origin='lower', cmap=plt.cm.hot)
plt.imshow(in_limbo[:, y_slice, :].T, origin='lower', cmap=plt.cm.Greens_r)
plt.axis('off')
plt.subplot(133)
plt.imshow(brain[:, :, z_slice].T, origin='lower', cmap=plt.cm.gray)
plt.imshow(z_map[:, :, z_slice].T, origin='lower', cmap=plt.cm.hot)
plt.imshow(in_limbo[:, :, z_slice].T, origin='lower', cmap=plt.cm.Greens_r)
plt.axis('off')
## Sandwich approach
plt.figure(figsize=(20, 5))
in_limbo = nb.load('/home/gdholla1/projects/in_limbo/final_results/%s_in_limbo_mask_brain_masked.nii.gz' % threshold).get_data()
z_map = nb.load('/home/gdholla1/projects/in_limbo/final_results/%s_z_threshold.nii.gz' % threshold).get_data()
z_map = np.ma.masked_less_equal(z_map, 0)
in_limbo = np.ma.masked_less_equal(in_limbo, 0)
plt.suptitle('Sandwich- estimation z > %s' % threshold, fontsize=25)
plt.subplot(131)
plt.imshow(brain[x_slice, :, :].T, origin='lower', cmap=plt.cm.gray)
plt.imshow(z_map[x_slice, :, :].T, origin='lower', cmap=plt.cm.hot)
plt.imshow(in_limbo[x_slice, :, :].T, origin='lower', cmap=plt.cm.Greens_r)
plt.axis('off')
plt.subplot(132)
plt.imshow(brain[:, y_slice, :].T, origin='lower', cmap=plt.cm.gray)
plt.imshow(z_map[:, y_slice, :].T, origin='lower', cmap=plt.cm.hot)
plt.imshow(in_limbo[:, y_slice, :].T, origin='lower', cmap=plt.cm.Greens_r)
plt.axis('off')
plt.subplot(133)
plt.imshow(brain[:, :, z_slice].T, origin='lower', cmap=plt.cm.gray)
plt.imshow(z_map[:, :, z_slice].T, origin='lower', cmap=plt.cm.hot)
plt.imshow(in_limbo[:, :, z_slice].T, origin='lower', cmap=plt.cm.Greens_r)
plt.axis('off')
(Due to standard coding this is radiological convention) with right hemisphere on right side of the screen.