[pymvpa] Pattern localization

[Matthias Ekman]

Hi,

I am trying to create a "discriminative map" which shows which features 
are driving the classifiers performance. The procedure (below) will give 
me one map per Nfold. My idea was, to calculate a mean map of those maps 
and (may be) threshold it at an arbitrary level.

The alternative I was thinking about, is a "stability analysis", that 
picks voxels which have consistent clf weights across cross-validation.

What do you think would be an appropriate way (for linear classifiers) 
of visualizing voxels which most strongly affect classifier predictions?


Best regards,
Matthias




from mvpa.suite import *

# ------------------------------------------------------------- #

# load PyMVPA example dataset
attr = SampleAttributes(os.path.join(pymvpa_dataroot, 'attributes.txt'))
dataset = NiftiDataset(samples=os.path.join(pymvpa_dataroot, 'bold.nii.gz'),
                        labels=attr.labels,
                        chunks=attr.chunks,
                        mask=os.path.join(pymvpa_dataroot, 'mask.nii.gz'))
					   				
detrend(dataset, perchunk=True, model='linear')

zscore(dataset, perchunk=True, baselinelabels=[0],
        targetdtype='float32')

dataset = dataset.selectSamples(
                 N.array([l in [1, 2] for l in dataset.labels],
                         dtype='bool'))
						
# ------------------------------------------------------------- #

clf = FeatureSelectionClassifier(
		   LinearCSVMC(),
            		SensitivityBasedFeatureSelection(
               			 OneWayAnova(transformer=N.abs),
				   FixedNElementTailSelector(50, mode='select', tail='upper')),
            	enable_states = ['feature_ids'])	
        			
# ------------------------------------------------------------- #
	
cv = CrossValidatedTransferError(TransferError(clf),
                                  NFoldSplitter(),								 
enable_states=['confusion'], 
harvest_attribs=['transerror.clf.getSensitivityAnalyzer(transformer=N.abs, 
force_training=False)()'])
								 								
# ------------------------------------------------------------- #	

error = cv(dataset)
print 'Error:', error

sensitivities = cv.harvested.values()[0]
dataset.map2Nifti(sensitivities).save('sensitivities.nii.gz')