tutorial_searchlight.ipynb

{
  "metadata": {
    "name": "Looking here and there -- Searchlights"
  }, 
  "nbformat": 3, 
  "nbformat_minor": 0, 
  "worksheets": [
    {
      "cells": [
        {
          "cell_type": "heading", 
          "level": 1, 
          "metadata": {}, 
          "source": [
            "Looking here and there -- Searchlights"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "In ", 
            "*chap_tutorial_classifiers* we have seen how we can implement a\nclassification analysis, but we still have no clue about where in the brain (or\nour chosen ROIs) our signal of interest is located.  And that is despite the\nfact that we have analyzed the data repeatedly, with different classifiers and\ninvestigated error rates and confusion matrices. So what can we do?\n\n", 
            "Ideally, we would like to have some way of estimating a score for each feature\nthat indicates how important that particular feature (most of the time a\nvoxel) is in the context of a certain classification task. There are various\npossibilities to get a vector of such per-feature scores in PyMVPA. We could\nsimply compute an ", 
            "[ANOVA](http://en.wikipedia.org/wiki/Analysis_of_variance) F-score per each feature, yielding scores that would\ntell us which features vary significantly between any of the categories in our\ndataset.\n\n", 
            "Before we can take a look at the implementation details, let's first recreate\nour preprocessed demo dataset. The code is very similar to that from\n", 
            "*chap_tutorial_classifiers* and should raise no questions. We get a\ndataset with one sample per category per run."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "from mvpa2.tutorial_suite import *\n", 
            "ds = get_haxby2001_data(roi='vt')\n", 
            "ds.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "Measures"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now that we have the dataset, computing the desired ANOVA F-scores is\nrelatively painless:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "aov = OneWayAnova()\n", 
            "f = aov(ds)\n", 
            "print f"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "If the code snippet above is of no surprise then you probably got the basic\nidea. We created an object instance `aov` being a\n", 
            "[OneWayAnova](http://pymvpa.org/generated/mvpa2.measures.anova.OneWayAnova.html#mvpa2-measures-anova-onewayanova). This instance is subsequently\n", 
            "*called* with a dataset and yields the F-scores wrapped into a\n", 
            "[Dataset](http://pymvpa.org/generated/mvpa2.datasets.base.Dataset.html#mvpa2-datasets-base-dataset). Where have we seen this before?  Right!\nThis one differs little from a call to\n", 
            "[CrossValidation](http://pymvpa.org/generated/mvpa2.measures.base.CrossValidation.html#mvpa2-measures-base-crossvalidation).  Both are objects that get\ninstantiated (potentially with some custom arguments) and yield the results in\na dataset when called with an input dataset. This is called a ", 
            "[processing](http://pymvpa.org/glossary.html#term-processing) and is a common concept in PyMVPA.\n\n", 
            "However, there is a difference between the two processing objects.\n", 
            "[CrossValidation](http://pymvpa.org/generated/mvpa2.measures.base.CrossValidation.html#mvpa2-measures-base-crossvalidation) returns a dataset with a single\nfeature -- the accuracy or error rate, while\n", 
            "[OneWayAnova](http://pymvpa.org/generated/mvpa2.measures.anova.OneWayAnova.html#mvpa2-measures-anova-onewayanova) returns a vector with one value per\nfeature. The latter is called a\n", 
            "[FeaturewiseMeasure](http://pymvpa.org/generated/mvpa2.measures.base.FeaturewiseMeasure.html#mvpa2-measures-base-featurewisemeasure). But other than the number of\nfeatures in the returned dataset there is not much of a difference. All\nmeasures in PyMVPA, for example, support an optional post-processing step.\nDuring instantiation of a measure an arbitrary mapper can be specified to be\ncalled internally to forward-map the results before they are returned. If, for\nsome reason, the F-scores need to be scaled into the interval [0,1], an\n", 
            "[FxMapper](http://pymvpa.org/generated/mvpa2.mappers.fx.FxMapper.html#mvpa2-mappers-fx-fxmapper) can be used to achieve that:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "aov = OneWayAnova(\n        postproc=FxMapper('features',\n                          lambda x: x / x.max(),\n                          attrfx=None))\n", 
            "f = aov(ds)\n", 
            "print f.samples.max()"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n**Exercise**", 
            "\n\n", 
            "Map the F-scores back into a brain volume and look at their distribution\nin the ventral temporal ROI."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# you can use this cell for this exercise"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n", 
            "\n\n", 
            "Now that we know how to compute feature-wise F-scores we can start worrying\nabout them. Our original goal was to decipher information that is encoded\nin the multivariate pattern of brain activation. But now we are using an\nANOVA, a ", 
            "**univariate** measure, to localize important voxels? There must\nbe something else -- and there is!"
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "Searching, searching, searching, ..."
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "*Kriegeskorte et al. (2006)* suggested an algorithm that takes\na small, sphere-shaped neighborhood of brain voxels and computes a\nmultivariate measure to quantify the amount of information encoded in its\npattern (e.g.  ", 
            "[mutual information](http://en.wikipedia.org/wiki/Mutual_information)). Later on this ", 
            "[searchlight](http://pymvpa.org/glossary.html#term-searchlight)\napproach has been extended to run a full classifier cross-validation in\nevery possible sphere in the brain. Since that, numerous studies have\nemployed this approach to localize relevant information in a locally\nconstraint fashion.\n\n", 
            "We know almost all pieces to implement a searchlight analysis in\nPyMVPA. We can load and preprocess datasets, we can set up a\ncross-validation procedure."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "clf = kNN(k=1, dfx=one_minus_correlation, voting='majority')\n", 
            "cv = CrossValidation(clf, HalfPartitioner())"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "The only thing left to do is that we have to split the dataset into all\npossible sphere neighborhoods that intersect with the brain. To achieve this,\nwe can use ", 
            "[sphere_searchlight()](http://pymvpa.org/generated/mvpa2.measures.searchlight.sphere_searchlight.html#mvpa2-measures-searchlight-sphere-searchlight):"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "sl = sphere_searchlight(cv, radius=3, postproc=mean_sample())"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "This single line configures a searchlight analysis that runs a full\ncross-validation in every possible sphere in the dataset. Each sphere has a\nradius of three voxels. The algorithm uses the coordinates (by default\n`voxel_indices`) stored in a feature attribute of the input dataset to\ndetermine local neighborhoods. From the `postproc` argument you might\nhave guessed that this object is also a measure -- and your are right. This\nmeasure returns whatever value is computed by the basic measure (here this\nis a cross-validation) and assigns it to the feature representing the\ncenter of the sphere in the output dataset. For this initial example we are\nnot interested in the full cross-validation output (error per each fold),\nbut only in the mean error, hence we are using an appropriate mapper for\npost-processing. As with any other ", 
            "[processing object](http://pymvpa.org/glossary.html#term-processing-object) we have to\ncall it with a dataset to run the actual analysis:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "res = sl(ds)\n", 
            "print res"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "That was it. However, this was just a toy example with only our ventral\ntemporal ROI. Let's now run it on a much larger volume, so we can actually\nlocalize something (even loading and preprocessing will take a few seconds).\nWe will reuse the same searchlight setup and run it on this data as well.\nDue to the size of the data it might take a few minutes to compute the\nresults, depending on the number of CPUs in the system."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "ds = get_haxby2001_data_alternative(roi=0)\n", 
            "print ds.nfeatures"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "res = sl(ds)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now let's see what we got. Since a vector with 35k elements is a little\nhard to comprehend we have to resort to some statistics."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "sphere_errors = res.samples[0]\n", 
            "res_mean = np.mean(res)\n", 
            "res_std = np.std(res)\n", 
            "chance_level = 1.0 - (1.0 / len(ds.uniquetargets))"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "As you'll see, the mean empirical error is just barely below the chance level.\nHowever, we would not expect a signal for perfect classification\nperformance in all spheres anyway. Let's see for how many spheres the error\nis more the two standard deviations lower than chance."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "frac_lower = np.round(np.mean(sphere_errors < chance_level - 2 * res_std), 3)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "So in almost 10% of all spheres the error is substantially lower than what\nwe would expect for random guessing of the classifier -- that is more than\n3000 spheres!"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n**Exercise**", 
            "\n\n", 
            "Look at the distribution of the errors\n(hint: `hist(sphere_errors, bins=np.linspace(0, 1, 18))`.\nIn how many spheres do you think the classifier actually picked up\nreal signal? What would be a good value to threshold the errors to\ndistinguish false from true positives? Think of it in the context of\nstatistical testing of fMRI data results. What problems are we facing\nhere?\n\n", 
            "Once you are done thinking about that -- and only ", 
            "*after* you're done,\nproject the sphere error map back into the fMRI volume and look at it as\na brain overlay in your favorite viewer (hint: you might want to store\naccuracies instead of errors, if your viewer cannot visualize the lower\ntail of the distribution:\n`map2nifti(ds, 1.0 - sphere_errors).to_filename('sl.nii.gz')`).\nDid looking at the image change your mind?"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# you can use this cell for this exercise"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n"
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "For real!"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now that we have an idea of what can happen in a searchlight analysis,\nlet's do another one, but this time on a more familiar ROI -- the full brain."
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n**Exercise**", 
            "\n\n", 
            "Load the dataset with `get_haxby2001_data_alternative(roi='brain')`\nthis will apply any required preprocessing for you. Now run a searchlight\nanalysis for radii 0, 1 and 3. For each resulting error map look at the\ndistribution of values, project them back into the fMRI volume and\ncompare them. How does the distribution change with radius and how does\nit compare to results of the previous exercise? What would be a good\nchoice for the threshold in this case?"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# you can use this cell for this exercise"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n", 
            "\n\n", 
            "You have now performed a number of searchlight analyses, investigated the\nresults and probably tried to interpret them. What conclusions did you draw\nfrom these analyses in terms of the neuroscientific aspects? What have you\nlearned about object representation in the brain? In this case we have run\n8-way classification analyses and have looked at the average error rate across\nall conditions in thousands of sphere-shaped ROIs in the brain. In some spheres the\nclassifier could perform well, i.e. it could predict all\nsamples equally well. However, this only applies to a handful of over 30k\nspheres we have tested, and does not reveal whether the classifier was capable of\nclassifying ", 
            "*all* of the conditions or just some.  For the vast majority\nwe observe errors somewhere\nbetween the theoretical chance level and zero and we don't know what caused\nthe error to decrease. We don't even know which samples get misclassified.\n\n", 
            "From ", 
            "*chap_tutorial_classifiers* we know that there is a way out of this\ndilemma. We can look at the confusion matrix of a classifier to get a lot more\ninformation that is otherwise hidden. However, we cannot reasonably do this for\nthousands of searchlight spheres (Note that this is not completely true. See\ne.g. ", 
            "*Connolly et al., 2012* for some creative use-cases for\nsearchlights).  It becomes obvious that a searchlight analysis is probably not\nthe end of a data exploration but rather a crude take off, as it raises more\nquestions than it answers.\n\n", 
            "Moreover, a searchlight cannot detect signals that extend beyond a small\nlocal neighborhood. This property effectively limits the scope of analyses\nthat can employ this strategy. A study looking a global brain circuitry\nwill hardly restrict the analysis to patches of a few cubic millimeters of\nbrain tissue. As we have seen before, searchlights also have another nasty\naspect. Although they provide us with a multivariate localization measure,\nthey also inherit the curse of univariate fMRI data analysis --\n", 
            "[multiple comparisons](http://en.wikipedia.org/wiki/Multiple_comparisons). PyMVPA comes with an algorithm that can help to\ncope with the problem in the context of group analyses:\n", 
            "[GroupClusterThreshold](http://pymvpa.org/generated/mvpa2.algorithms.group_clusterthr.GroupClusterThreshold.html#mvpa2-algorithms-group-clusterthr-groupclusterthreshold).\n\n", 
            "Despite these limitations a searchlight analysis can be a valuable exploratory\ntool if used appropriately. The capabilities of PyMVPA's searchlight\nimplementation go beyond what we looked at in this tutorial. It is not only\npossible to run ", 
            "*spatial* searchlights, but multiple spaces can be considered\nsimultaneously. This is further illustrated in\n", 
            "*chap_tutorial_eventrelated_searchlight*."
          ]
        }
      ], 
      "metadata": {}
    }
  ]
}