Skip to content

Instantly share code, notes, and snippets.

Created Apr 16, 2010
Embed
What would you like to do?
#Now I download file into pymvpa
from mvpa.suite import *
#filename of source fMRI timeseries image for each run
fmri_src1 = os.path.join(pymvpa_dataroot, '/home/cabiubuntu/temiadelore/s101CS_12runs/4D_s101
CS.nii.gz')
# load the samples attributes as usual and preserve the
# literal labels
attr = SampleAttributes(
os.path.join(pymvpa_dataroot,
'/home/cabiubuntu/temiadelore/s101CS_12runs/attributes.txt'),
literallabels=True)
#setting up timeseries events
evs = [ev for ev in attr.toEvents()
if ev['label'] in ['AI', 'VI', 'AC', 'VC', 'NULL']]
#Since we might want to look at the sensitivity profile ranging from just a little before to a little after each event,
#I will extract starting from ove volume prior to the event and one
#volume after the events i.e. two consecutive volumes are extracted a as sample for each event
verbose(1, "Segmenting timeseries into events")
fmri_src_ev = ERNiftiDataset(samples=fmri_src,
events=evs,
mask=None,
evconv = True,
enforce_dim=None,
labels_map={'AI': 1,
'VI': 2,
'AC': 3,
'VC': 4,
'NULL' : 0})
#Begin Preprocessing
# detrend on full timeseries
detrend(fmri_src_ev, perchunk=True, model='linear')
#We perform normalization by Z-scoring the data
# using rest/NULL as baseline
zscore(fmri_src_ev, perchunk=True)
#Using meta-classifier
#using k nearest neighbor as clasifier
from mvpa.clfs.knn import kNN
from mvpa.datasets.splitters import NFoldSplitter, OddEvenSplitter
from mvpa.clfs.meta import SplitClassifier, FeatureSelectionClassifier
#Performing a cross-validated sensitivity analysis.
#NFoldSplitter used to split the dataset prior to training
#cvtype = 1 is the leave one out cross validation split
clf = FeatureSelectionClassifier(
kNN(k=5),
SensitivityBasedFeatureSelection(OneWayAnova(transformer=N.abs),
FractionTailSelector(0.05, mode='select', tail='upper')),
enable_states = ['feature_ids'],descr="kNN on 5%(ANOVA)")
sclf = SplitClassifier(clf, NFoldSplitter(cvtype = 1),
enable_states=['confusion', 'training_confusion','feature_ids','trained_dataset'])
#Now train data
sclf.train(fmri_src_ev)
#The features chosen to predict test data
final_dataset = fmri_src_ev.selectFeatures(sclf.feature_ids)
#Checking the performance of the classifier
print sclf.confusion.asstring(description=True) \
#or plot the performance results
import pylab as P
sclf.confusion.plot() \
P.show() \
#Calculate the cross validated transfer error of the trained classifier
from mvpa.algorithms.cvtranserror import CrossValidatedTransferError
from mvpa.datasets.splitters import NFoldSplitter
clf = KNN(k=5) # you can change classifier type
terr = TransferError(clf)
cvterr = CrossValidatedTransferError(terr,
NFoldSplitter(cvtype=1),
enable_states=['confusion'])
error = cvterr(fmri_src_ev)
#Now I download file into pymvpa
from mvpa.suite import *
#filename of source fMRI timeseries image for each run
fmri_src = os.path.join(pymvpa_dataroot, '/home/cabiubuntu/temiadelore/s101CS_12runs/4D_s101CS.nii.gz')
mask = os.path.join(pymvpa_dataroot, '/home/cabiubuntu/temiadelore/ROIS/AudVis_ROI.nii')
T1 = os.path.join(pymvpa_dataroot, '/home/cabiubuntu/temiadelore/s101CS_12runs/TI.nii.gz')
# load the samples (each scan) attributes as usual and preserve the
# literal labels
attr = SampleAttributes(
os.path.join(pymvpa_dataroot,
'/home/cabiubuntu/temiadelore/s101CS_12runs/attributes.txt'),
literallabels=True)
# actual labels are not important here, could be 'labels=1'
pre_fmri_src = NiftiDataset(samples=fmri_src, labels=attr.labels, chunks=attr.chunks)
# convert to floats
pre_fmri_src.setSamplesDType('float32')
#Begin Preprocessing
# detrend on full timeseries
detrend(pre_fmri_src, perchunk=True, model='linear')
# zscore dataset relative to baseline ('rest') mean
zscore(pre_fmri_src, perchunk=True, baselinelables='NULL', targetdtype='float32')
#load the reference function
from mvpa.misc.io.base import *
ev = SampleAttributes(os.path.join(pymvpa_dataroot,'/home/cabiubuntu/temiadelore/s101CS_12runs/events.txt'),header=['labels','onset','chunks'])
evs = ev.toEvents()
#Segmenting timeseries into events
fmri_src_ev = ERNiftiDataset(samples=pre_fmri_src.map2Nifti(),
events=evs,
mask=mask,
evconv = False,
enforce_dim=None,
labels_map={'AI': 1,
'VI': 2,
'AC': 3,
'VC': 4,
'NULL' : 0})
# remove baseline 'rest' samples from dataset for analysis
fmri_src_ev = fmri_src_ev.selectSamples(N.array([l != 0 for l in fmri_src_ev.labels],
dtype='bool'))
#get dataset summary
fmri_src_ev.summary()
#Using meta-classifier to select features and choose classifier type
#using k nearest neighbor as classifier (k = 5 , number of neighboring voxels used)
from mvpa.clfs.knn import kNN
from mvpa.datasets.splitters import NFoldSplitter, OddEvenSplitter
from mvpa.clfs.meta import SplitClassifier, FeatureSelectionClassifier
#Performing a cross-validated sensitivity analysis.
#NFoldSplitter used to split the dataset prior to training
#cvtype = 1 is the leave one out cross validation split
clf = FeatureSelectionClassifier(
kNN(k=5),
SensitivityBasedFeatureSelection(OneWayAnova(transformer=N.abs),
FractionTailSelector(0.05, mode='select', tail='upper')),
enable_states = ['feature_ids'],descr="kNN on 5%(ANOVA)")
sclf = SplitClassifier(clf, NFoldSplitter(cvtype = 1),
enable_states=['confusion', 'training_confusion','feature_ids','trained_dataset'])
#need to register pickler/unpickler functions for ufuncs
import copy_reg
import numpy
def ufunc_pickler(ufunc):
return ufunc.__name__
def ufunc_unpickler(name):
import numpy
return getattr(numpy, name)
copy_reg.pickle(numpy.ufunc, ufunc_pickler, ufunc_unpickler)
import cPickle
cPickle.dumps(numpy.add)
'cnumpy.core.umath\nadd\np1\n.'
cPickle.loads(cPickle.dumps(numpy.add))
#Now train data
sclf.train(fmri_src_ev)
#The features chosen to predict test data
final_dataset = fmri_src_ev.selectFeatures(sclf.feature_ids)
#convert to nfti image features selected
nftifeaures = fmri_src_ev.map2Nifti(final_dataset)
#visualizing/p;lotting features space on subject's brain
fig = P.figure(figsize=(12, 4), facecolor='white')
plotMRI(background=T1, background_mask=None, cmap_bg='gray', overlay=nftifeatures, overlay_mask=mask, cmap_overlay='autumn', vlim=(None,None), vlim_type=None, do_stretch_colors=True, add_info=True, add_hist=True, add_colorbar=True, fig=fig, interactive=None, nrows=None, ncolumns=None)
P.title('Feature Selection map')
#Checking the performance of the classifier
print sclf.confusion.asstring(description=True) \
#or plot the performance results
import pylab as P
sclf.confusion.plot() \
P.show() \
#Calculate the cross validated transfer error of the trained classifier
from mvpa.algorithms.cvtranserror import CrossValidatedTransferError
from mvpa.datasets.splitters import NFoldSplitter
clf = KNN(k=5)
terr = TransferError(clf)
cvterr = CrossValidatedTransferError(terr,
NFoldSplitter(cvtype=1),
enable_states=['confusion'])
error = cvterr(fmri_src_ev)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment