<div dir="ltr"><div><br class="gmail-Apple-interchange-newline"><br></div>- 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)<div><br></div><div>- Regarding chunks: We actually have two runs (separated by the "7" labels in the middle - so</div><div>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.)</div><div><br></div><div>- 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 <a href="https://www.dropbox.com/sh/6qnrt2l2othc83g/AADBpc-eJSK5893Vrz_A8BS1a?dl=0">link</a> as also already previously shared):</div><div><br></div><div><br></div><div><br></div><div><div><br></div><div>#</div><div><br></div><div>from glob import glob</div><div>import os</div><div>import numpy as np</div><div><br></div><div>from mvpa2.suite import *</div><div><br></div><div>%matplotlib inline</div><div><br></div><div><br></div><div># enable debug output for searchlight call</div><div>if __debug__:</div><div> debug.active += ["SLC"]</div><div><br></div><div><br></div><div># change working directory to 'WB'</div><div>os.chdir('mypath/WB')</div><div><br></div><div># use glob to get the filenames of .nii data into a list</div><div>nii_fns = glob('beta*.nii')</div><div><br></div><div># read data</div><div><br></div><div>labels = [</div><div> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, </div><div> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, </div><div> 7, 7, 7, 7, 7, 7, 7, </div><div> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, </div><div> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, </div><div> 7, 7, 7, 7, 7, 7, 7</div><div>]</div><div>grps = np.repeat([0, 1], 37, axis=0) # used for `chuncks`</div><div><br></div><div>db = mvpa2.datasets.mri.fmri_dataset(</div><div> nii_fns, targets=labels, chunks=grps, mask=None, sprefix='vxl', tprefix='tpref', add_fa=None</div><div>)</div><div><br></div><div># use only the samples of which labels are 1 or 2</div><div>db12 = db[np.array([label in [1, 2] for label in labels], dtype='bool')]</div><div><br></div><div># in-place z-score normalization</div><div>zscore(db12)</div><div><br></div><div># choose classifier</div><div>clf = LinearNuSVMC()</div><div><br></div><div># setup measure to be computed by Searchlight</div><div># cross-validated mean transfer using an N-fold dataset splitter</div><div>cv = CrossValidation(clf, NFoldPartitioner())</div><div><br></div><div># define searchlight methods</div><div>radius_ = 1</div><div>sl = sphere_searchlight(</div><div> cv, radius=radius_, space='vxl_indices',</div><div> postproc=mean_sample()</div><div>)</div><div><br></div><div># stripping all attributes from the dataset that are not required for the searchlight analysis</div><div>db12s = db12.copy(</div><div> deep=False,</div><div> sa=['targets', 'chunks'],</div><div> fa=['vxl_indices'],</div><div> a=['mapper']</div><div>)</div><div><br></div><div># run searchlight</div><div>sl_map = sl(db12s)</div><div><br></div><div># toggle between error rate (searchlight output) and accuracy (for plotting)</div><div>sl_map.samples *= -1</div><div>sl_map.samples += 1</div><div><br></div><div># The result dataset is fully aware of the original dataspace. </div><div># Using this information we can map the 1D accuracy maps back into “brain-space” (using NIfTI image header information from the original input timeseries.</div><div>niftiresults = map2nifti(sl_map, imghdr=db12.a.imghdr)</div><div><br></div><div># find the 3-d coordinates of extrema of error rate or accuracy</div><div>extremum_i = np.argmax(sl_map.samples[0]) # max</div><div>extremum_i = np.argmin(sl_map.samples[0]) # min</div><div>coord = db12s.fa[list(db.fa.keys())[0]][extremum_i]</div><div><br></div><div># plotting</div><div>plot_args = {</div><div># 'background' : 'Subject2_Struct.nii',</div><div># 'background_mask' : 'brain.nii.gz',</div><div># 'overlay_mask' : 'S1_mask.nii',</div><div> 'do_stretch_colors' : False,</div><div> 'cmap_bg' : 'gray',</div><div> 'cmap_overlay' : 'autumn', # YlOrRd_r # pl.cm.autumn</div><div> 'interactive' : cfg.getboolean('examples', 'interactive', True)</div><div>}</div><div>fig = pl.figure(figsize=(12, 4), facecolor='white')</div><div>subfig = plot_lightbox(</div><div> overlay=niftiresults,</div><div> vlim=(0.5, None), slices=range(23,31),</div><div> fig=fig, </div><div> background='Subject2_Struct.nii',</div><div># background_mask='brain.nii.gz',</div><div># overlay_mask='S1_mask.nii',</div><div> **plot_args</div><div>)</div><div>pl.title('Accuracy distribution for radius %i' % radius_)</div><div><br></div><div><br></div><div><br></div><div><br></div><div>#</div></div><div><br></div><div><br></div><div><br><div><br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Sep 8, 2017 at 2:53 PM, Yaroslav Halchenko <span dir="ltr"><<a href="mailto:debian@onerussian.com" target="_blank">debian@onerussian.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class=""><br>
On Fri, 08 Sep 2017, Pegah Kassraian Fard wrote:<br>
<br>
> That would be great.<br>
> - Files (whole brain - not yet masked) are here<br>
> - Targets/labels are in the labels.mat file, it is binary classification,<br>
> so one classes' samples are labeled by "1", the other classes by "2", the<br>
> "7" are just the head movement parameters and can be ignored for<br>
> classification<br>
> - I have included a struct image of the subject (subject2_struct.nii), as<br>
> well as a mask for primary sensory cortex (S1_mask.nii)<br>
> I was attempting to reproduce this:<br>
> <a href="http://www.pymvpa.org/examples/searchlight.html" rel="noreferrer" target="_blank">http://www.pymvpa.org/<wbr>examples/searchlight.html</a><br>
> including the plots, but my highest accuracies (1-error) were 1. outside<br>
> of the brain:))<br>
<br>
</span>;) that is the fun of using all data (not just masked) -- you could<br>
potentially discover interesting aspects of your data/design/code ;-)<br>
does the rest of the map look "feasible" given your experiment?<br>
<br>
but before we get that excited, let's indeed see the code:<br>
<span class=""><br>
> 2. I was not sure what the necessary background files were<br>
> (I tried to guess and use similar ones I found elsewhere, but there was a<br>
> dimension mismatch - is there a resampling method in the toolbox to match<br>
> the dimensions?)<br>
> Many thanks, let me know if I can provide my (buggy?) script<br>
<br>
</span>sure -- just paste it in here (make sure no wrapping) or post it<br>
somewhere online<br>
<span class="im HOEnZb"><br>
--<br>
Yaroslav O. Halchenko<br>
Center for Open Neuroscience <a href="http://centerforopenneuroscience.org" rel="noreferrer" target="_blank">http://<wbr>centerforopenneuroscience.org</a><br>
Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755<br>
Phone: <a href="tel:%2B1%20%28603%29%20646-9834" value="+16036469834">+1 (603) 646-9834</a> Fax: <a href="tel:%2B1%20%28603%29%20646-1419" value="+16036461419">+1 (603) 646-1419</a><br>
WWW: <a href="http://www.linkedin.com/in/yarik" rel="noreferrer" target="_blank">http://www.linkedin.com/in/<wbr>yarik</a><br>
<br>
</span><div class="HOEnZb"><div class="h5">______________________________<wbr>_________________<br>
Pkg-ExpPsy-PyMVPA mailing list<br>
<a href="mailto:Pkg-ExpPsy-PyMVPA@lists.alioth.debian.org">Pkg-ExpPsy-PyMVPA@lists.<wbr>alioth.debian.org</a><br>
<a href="http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-exppsy-pymvpa" rel="noreferrer" target="_blank">http://lists.alioth.debian.<wbr>org/cgi-bin/mailman/listinfo/<wbr>pkg-exppsy-pymvpa</a><br>
</div></div></blockquote></div><br></div>