[pymvpa] searchlight analysis

Fri Sep 8 15:52:15 UTC 2017

- It's not our model's fault:-)! More seriously, we have checked the basic
contrasts in SPM, and I had also performed previously "default"
classification (linear SVM over all voxels in a relevant ROI vs.
classification in a "non-relevant" area as sanity check, permutation tests
of classification labels etc.). I have also picked a very simple design for
starters (it is a localizer, two kinds of sensory stimulation delivered in
random order, and jittered rest in between)

- Regarding chunks: We actually have two runs (separated by the "7" labels
in the middle - so
15 x "1" for class 1, 15 x "2" for class 2, then 6 x "7" indicating
head-movement parameters, then the next run: again 15 x "1" for class 1, 15
x "2" for class 2, 8 x "7" indicating head-movement parameters and
regression intercepts). Though I understood that chunks can also be defined
if there is only one run, through sensible partitioning the data (without
obviously losing the mapping to the labels/targets etc.)

- Below the code (I am new to python, hope it's +/- readable), I have also
attached it (in the last step the plotting, which is also not yet
functional, is added as well, all the files necessary are in the link
<https://www.dropbox.com/sh/6qnrt2l2othc83g/AADBpc-eJSK5893Vrz_A8BS1a?dl=0> as
also already previously shared):

#

from glob import glob
import os
import numpy as np

from mvpa2.suite import *

%matplotlib inline

# enable debug output for searchlight call
if __debug__:
    debug.active += ["SLC"]

# change working directory to 'WB'
os.chdir('mypath/WB')

# use glob to get the filenames of .nii data into a list
nii_fns = glob('beta*.nii')

# read data

labels = [
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
    7, 7, 7, 7, 7, 7, 7,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
    7, 7, 7, 7, 7, 7, 7
]
grps = np.repeat([0, 1], 37, axis=0) # used for `chuncks`

db = mvpa2.datasets.mri.fmri_dataset(
    nii_fns, targets=labels, chunks=grps, mask=None, sprefix='vxl',
tprefix='tpref', add_fa=None
)

# use only the samples of which labels are 1 or 2
db12 = db[np.array([label in [1, 2] for label in labels], dtype='bool')]

# in-place z-score normalization
zscore(db12)

# choose classifier
clf = LinearNuSVMC()

# setup measure to be computed by Searchlight
# cross-validated mean transfer using an N-fold dataset splitter
cv = CrossValidation(clf, NFoldPartitioner())

# define searchlight methods
radius_ = 1
sl = sphere_searchlight(
    cv, radius=radius_, space='vxl_indices',
    postproc=mean_sample()
)

# stripping all attributes from the dataset that are not required for the
searchlight analysis
db12s = db12.copy(
    deep=False,
    sa=['targets', 'chunks'],
    fa=['vxl_indices'],
    a=['mapper']
)

# run searchlight
sl_map = sl(db12s)

# toggle between error rate (searchlight output) and accuracy (for plotting)
sl_map.samples *= -1
sl_map.samples += 1

# The result dataset is fully aware of the original dataspace.
# Using this information we can map the 1D accuracy maps back into
“brain-space” (using NIfTI image header information from the original input
timeseries.
niftiresults = map2nifti(sl_map, imghdr=db12.a.imghdr)

# find the 3-d coordinates of extrema of error rate or accuracy
extremum_i = np.argmax(sl_map.samples[0]) # max
extremum_i = np.argmin(sl_map.samples[0]) # min
coord = db12s.fa[list(db.fa.keys())[0]][extremum_i]

# plotting
plot_args = {
#     'background' : 'Subject2_Struct.nii',
#     'background_mask' : 'brain.nii.gz',
#     'overlay_mask' : 'S1_mask.nii',
    'do_stretch_colors' : False,
    'cmap_bg' : 'gray',
    'cmap_overlay' : 'autumn', # YlOrRd_r # pl.cm.autumn
    'interactive' : cfg.getboolean('examples', 'interactive', True)
}
fig = pl.figure(figsize=(12, 4), facecolor='white')
subfig = plot_lightbox(
    overlay=niftiresults,
    vlim=(0.5, None), slices=range(23,31),
    fig=fig,
    background='Subject2_Struct.nii',
#     background_mask='brain.nii.gz',
#     overlay_mask='S1_mask.nii',
    **plot_args
)
pl.title('Accuracy distribution for radius %i' % radius_)

#

On Fri, Sep 8, 2017 at 2:53 PM, Yaroslav Halchenko <debian at onerussian.com>
wrote:

>
> On Fri, 08 Sep 2017, Pegah Kassraian Fard wrote:
>
> >    That would be great.
> >    - Files (whole brain - not yet masked) are here
> >    - Targets/labels are in the labels.mat file, it is binary
> classification,
> >    so one classes' samples are labeled by "1", the other classes by "2",
> the
> >    "7" are just the head movement parameters and can be ignored for
> >    classification
> >    - I have included a struct image of the subject
> (subject2_struct.nii), as
> >    well as a mask for primary sensory cortex (S1_mask.nii)
> >    I was attempting to reproduce this:
> >    http://www.pymvpa.org/examples/searchlight.html
> >    including the plots, but my highest accuracies (1-error) were 1.
> outside
> >    of the brain:))
>
> ;)  that is the fun of using all data (not just masked) -- you could
> potentially discover interesting aspects of your data/design/code ;-)
> does the rest of the map look "feasible" given your experiment?
>
> but before we get that excited, let's indeed see the code:
>
> > 2. I was not sure what the necessary background files were
> >    (I tried to guess and use similar ones I found elsewhere, but there
> was a
> >    dimension mismatch - is there a resampling method in the toolbox to
> match
> >    the dimensions?)
> >    Many thanks, let me know if I can provide my (buggy?) script
>
> sure -- just paste it in here (make sure no wrapping) or post it
> somewhere online
>
> --
> Yaroslav O. Halchenko
> Center for Open Neuroscience     http://centerforopenneuroscience.org
> Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755
> Phone: +1 (603) 646-9834                       Fax: +1 (603) 646-1419
> WWW:   http://www.linkedin.com/in/yarik
>
> _______________________________________________
> Pkg-ExpPsy-PyMVPA mailing list
> Pkg-ExpPsy-PyMVPA at lists.alioth.debian.org
> http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-exppsy-pymvpa
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.alioth.debian.org/pipermail/pkg-exppsy-pymvpa/attachments/20170908/ccab0353/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: search_light.py
Type: text/x-python-script
Size: 2966 bytes
Desc: not available
URL: <http://lists.alioth.debian.org/pipermail/pkg-exppsy-pymvpa/attachments/20170908/ccab0353/attachment.bin>