permutation_test.ipynb

{
  "metadata": {
    "name": "Monte-Carlo testing of Classifier-based Analyses"
  }, 
  "nbformat": 3, 
  "nbformat_minor": 0, 
  "worksheets": [
    {
      "cells": [
        {
          "cell_type": "heading", 
          "level": 1, 
          "metadata": {}, 
          "source": [
            "Monte-Carlo testing of Classifier-based Analyses"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "It is often desirable to be able to make statements like ", 
            "*\"Performance is\nsignificantly above chance-level\"* and to help with that PyMVPA supports ", 
            "*Null*\nhypothesis (aka ", 
            "*H0*) testing for any ", 
            "[Measure](http://pymvpa.org/generated/mvpa2.measures.base.Measure.html#mvpa2-measures-base-measure).\nMeasures take an optional constructor argument `null_dist` that can be used\nto provide an instance of some ", 
            "[NullDist](http://pymvpa.org/generated/mvpa2.clfs.stats.NullDist.html#mvpa2-clfs-stats-nulldist) estimator.\nIf the properties of the expected ", 
            "*Null* distribution are known a-priori, it is\npossible to use any distribution specified in SciPy's `stats` module for this\npurpose (see e.g. ", 
            "[FixedNullDist](http://pymvpa.org/generated/mvpa2.clfs.stats.FixedNullDist.html#mvpa2-clfs-stats-fixednulldist)).\n\n", 
            "However, as with other applications of statistics in classifier-based analyses\nthere is the problem that we typically do not know the distribution of a\nvariable like error or performance under the ", 
            "*Null* hypothesis (i.e. the\nprobability of a result given that there is no signal), hence we cannot easily\nassign the adored p-values. Even worse, the chance-level or guess probability\nof a classifier depends on the content of a validation dataset, e.g. balanced\nor unbalanced number of samples per label and total number of labels.\n\n", 
            "One approach to deal with this situation is to ", 
            "*estimate* the ", 
            "*Null*\ndistribution using permutation testing. The ", 
            "*Null* distribution is then\nestimated by computing the measure of interest multiple times using original\ndata samples but with permuted targets.  Since quite often the exploration of\nall permutations is unfeasible, Monte-Carlo testing (see ", 
            "*Nichols*) allows to obtain stable estimate with only a\nlimited number of random permutations.\n\n", 
            "Given the results computed using permuted targets one can now determine the\nprobability of the empirical result (i.e. the one computed from the\noriginal training dataset) under the ", 
            "*no signal* condition. This is\nsimply the fraction of results from the permutation runs that is\nlarger or smaller than the empirical (depending on whether one is\nlooking at performances or errors).\n\n", 
            "Here is how this looks for a simple cross-validated classification in PyMVPA.\nWe start by generated a dataset with 200 samples and 3 features of which 2 carry\nsome relevant signal."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# lazy import\nfrom mvpa2.suite import *\n\n# enable progress output for MC estimation\nif __debug__:\n    debug.active += [\"STATMC\"]\n\n# some example data with signal\nds = normal_feature_dataset(perlabel=100, nlabels=2, nfeatures=3,\n                            nonbogus_features=[0,1], snr=0.3, nchunks=2)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now we can start collecting the pieces that play a role in this analysis. We\nneed a classifier."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "clf = LinearCSVMC()"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We need a ", 
            "[generator](http://pymvpa.org/glossary.html#term-generator) than will produce partitioned datasets, one for each\nfold of the cross-validation. A partitioned dataset is basically the same as the\noriginal dataset, but has an additional samples attribute that indicates whether\nparticular samples will be the ", 
            "*part* of the data that is used for training the\nclassifier, or for testing it. By default, the\n", 
            "[NFoldPartitioner](http://pymvpa.org/generated/mvpa2.generators.partition.NFoldPartitioner.html#mvpa2-generators-partition-nfoldpartitioner) will create a sample\nattribute `partitions` that will label one ", 
            "[chunk](http://pymvpa.org/glossary.html#term-chunk) in each fold\ndifferently from all others (hence mark it as taken-out for testing)."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "partitioner = NFoldPartitioner()"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We need two pieces for the Monte Carlo shuffling. The first of them is\nan instance of an\n", 
            "[AttributePermutator](http://pymvpa.org/generated/mvpa2.generators.permutation.AttributePermutator.html#mvpa2-generators-permutation-attributepermutator) that will\npermute the target attribute of the dataset for each iteration.  We\nwill instruct it to perform 200 permutations. In a real analysis the\nnumber of permutations should be larger to get stable estimates."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "permutator = AttributePermutator('targets', count=200)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "The second mandatory piece for a Monte-Carlo-style estimation of\nthe ", 
            "*Null* distribution is the actual \"estimator\".\n", 
            "[MCNullDist](http://pymvpa.org/generated/mvpa2.clfs.stats.MCNullDist.html#mvpa2-clfs-stats-mcnulldist) will use the\nconstructed `permutator` to shuffle the targets and later on report\np-value from the left tail of the ", 
            "*Null* distribution, because we are\ngoing to compute errors and are interested in them being ", 
            "*lower* than\nchance. Finally we also ask for all results from Monte-Carlo shuffling\nto be stored for subsequent visualization of the distribution."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "distr_est = MCNullDist(permutator, tail='left', enable_ca=['dist_samples'])"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now we have all pieces and can conduct the actual cross-validation. We assign\na post-processing ", 
            "[node](http://pymvpa.org/glossary.html#term-node) `mean_sample` that will take care of averaging\nerror values across all cross-validation fold. Consequently, the ", 
            "*Null*\ndistribution of ", 
            "*average cross-validated classification error* will be estimated\nand used for statistical evaluation."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "cv = CrossValidation(clf, partitioner,\n                     errorfx=mean_mismatch_error,\n                     postproc=mean_sample(),\n                     null_dist=distr_est,\n                     enable_ca=['stats'])\n# run\nerr = cv(ds)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now we have a usual cross-validation error and `cv` stores\n", 
            "[conditional attribute`s such as confusion matrices](http://pymvpa.org/glossary.html#term-conditional-attribute`s-such-as-confusion-matrices):"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print 'CV-error:', 1 - cv.ca.stats.stats['ACC']"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "However, in addition it also provides the results of the statistical\nevaluation. The ", 
            "[conditional attribute](http://pymvpa.org/glossary.html#term-conditional-attribute) `null_prob` has a\ndataset that contains the p-values representing the likelihood of an\nerror equal or lower to the output one under the ", 
            "*Null* hypothesis,\ni.e. no actual relevant signal in the data. For a reason that will\nappear sensible later on, the p-value is contained in a dataset."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "p = cv.ca.null_prob\n# should be exactly one p-value\nassert(p.shape == (1,1))\nprint 'Corresponding p-value:',  np.asscalar(p)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We can now look at the distribution of the errors under ", 
            "*H0* and plot the\nexpected chance level as well as the empirical error."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# make new figure\npl.figure()\n# histogram of all computed errors from permuted data\npl.hist(np.ravel(cv.null_dist.ca.dist_samples), bins=20)\n# empirical error\npl.axvline(np.asscalar(err), color='red')\n# chance-level for a binary classification with balanced samples\npl.axvline(0.5, color='black', ls='--')\n# scale x-axis to full range of possible error values\npl.xlim(0,1)\npl.xlabel('Average cross-validated classification error')"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We can see that the ", 
            "*Null* or chance distribution is centered around the\nexpected chance-level and the empirical error value is in the far left tail,\nthus unlikely to belong to ", 
            "*Null* distribution, and hence the low p-value.\n\n", 
            "This could be the end, but sometimes one needs to have a closer look. Let's say your\ndata is not that homogeneous. Let's say that some ", 
            "[chunk](http://pymvpa.org/glossary.html#term-chunk) may be very\ndifferent from others. You might want to look at the error value probability for\nspecific cross-validation folds. Sounds complicated? Luckily it is very simple.\nIt only needs a tiny change in the cross-validation setup -- the removal of the\n`mean_sample` post-processing ", 
            "[node](http://pymvpa.org/glossary.html#term-node)."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "cv = CrossValidation(clf, partitioner,\n                     errorfx=mean_mismatch_error,\n                     null_dist=distr_est,\n                     enable_ca=['stats'])\n# run\nerr = cv(ds)\n\nassert (err.shape == (2,1))\nprint 'CV-errors:', np.ravel(err)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now we get two errors -- one for each cross-validation fold and\nmost importantly, we also get the two associated p-values."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "p = cv.ca.null_prob\nassert(p.shape == (2,1))\nprint 'Corresponding p-values:',  np.ravel(p)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "What happened is that a dedicated ", 
            "*Null* distribution has been estimated for\neach element in the measure results. Without `mean_sample` an error is\nreported for each CV-fold, hence a separate distributions are estimated for\neach CV-fold too. And because we have also asked for all distribution samples\nto be reported, we can now plot both distribution and both empirical errors.\nBut how do we figure out with value is which?\n\n", 
            "As mentioned earlier all results are returned in Datasets. All datasets have\ncompatible sample and feature axes, hence corresponding elements."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "assert(err.shape == p.shape == cv.null_dist.ca.dist_samples.shape[:2])\n\n# let's make a function this time\ndef plot_cv_results(cv, err, title):\n    # make new figure\n    pl.figure()\n    colors = ['green', 'blue']\n    # null distribution samples\n    dist_samples = np.asarray(cv.null_dist.ca.dist_samples)\n    for i, e in enumerate(err):\n        # histogram of all computed errors from permuted data per CV-fold\n        pl.hist(np.ravel(dist_samples[i]), bins=20, color=colors[i],\n                label='CV-fold %i' %i, alpha=0.5,\n                range=(dist_samples.min(), dist_samples.max()))\n        # empirical error\n        pl.axvline(np.asscalar(e), color=colors[i])\n\n    # chance-level for a binary classification with balanced samples\n    pl.axvline(0.5, color='black', ls='--')\n    # scale x-axis to full range of possible error values\n    pl.xlim(0,1)\n    pl.xlabel(title)\n\nplot_cv_results(cv, err, 'Per CV-fold classification error')"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We have already seen that the statistical evaluation is pretty flexible.\nHowever, we haven't yet seen whether it is flexible enough. To illustrate that\nthink about what was done in the above Monte Carlo analyses.\n\n", 
            "A dataset was shuffled repeatedly, and for each iteration a full\ncross-validation of classification error was performed. However, the shuffling\nwas done on the ", 
            "*full* dataset, hence target were permuted in both training\n", 
            "*and* testing dataset portions in each CV-fold. This basically means that for\neach Monte Carlo iteration the classifier was tested on a new data/signal.\nHowever, we may be more interested in what the classifier has to say on the\n", 
            "*actual* data, but when it was trained on randomly permuted data.\n\n", 
            "As you can guess this is possible too and goes like this. The most important\ndifference is that we are going to use  now a dedicate measure to estimate the\n", 
            "*Null* distribution. That measure is very similar to the cross-validation we\nhave used before, but differs in an important bit. Let's look at the pieces."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# how often do we want to shuffle the data\nrepeater = Repeater(count=200)\n# permute the training part of a dataset exactly ONCE\npermutator = AttributePermutator('targets', limit={'partitions': 1}, count=1)\n# CV with null-distribution estimation that permutes the training data for\n# each fold independently\nnull_cv = CrossValidation(\n            clf,\n            ChainNode([partitioner, permutator], space=partitioner.get_space()),\n            errorfx=mean_mismatch_error)\n# Monte Carlo distribution estimator\ndistr_est = MCNullDist(repeater, tail='left', measure=null_cv,\n                       enable_ca=['dist_samples'])\n# actual CV with H0 distribution estimation\ncv = CrossValidation(clf, partitioner, errorfx=mean_mismatch_error,\n                     null_dist=distr_est, enable_ca=['stats'])"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "The `repeater` is a simple node that returns any given dataset a\nconfigurable number of times. We use the helper to configure the number of\nMonte Carlo iterations. The new `permutator` is again configured to shuffle\nthe `targets` attribute, but only ", 
            "*once* and only for samples that were\nlabeled as being part of the training set in a particular CV-fold (the\n`partitions` sample attribute will be created by the NFoldPartitioner that we\nhave configured earlier).\n\n", 
            "The most important difference is a new dedicated measure that will be used to\nperform a cross-validation analysis under the ", 
            "*H0* hypotheses. To this end\nwe set up a standard CV procedure with a twist: we use a chained generator\n(comprising of the typical partitioner and the new one-time permutator).\nThis will cause the CV to permute the training set for each CV-fold internally\n(and that is what we wanted).\n\n", 
            "Now we assign the ", 
            "*H0* cross-validation procedure to the distribution\nestimator and use the `repeater` to set the number of iterations. Lastly, we\nplug everything into a standard CV analysis with, again, a non-permuting\n`partitioner` and the pimped ", 
            "*Null* distribution estimator.\n\n", 
            "Now we just need to run it, and plot the results the same way we did before."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "err = cv(ds)\nprint 'CV-errors:', np.ravel(err)\np = cv.ca.null_prob\nprint 'Corresponding p-values:',  np.ravel(p)\n# plot\nplot_cv_results(cv, err,\n                'Per CV-fold classification error (only training permutation)')"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "There a many ways to futher tweak the statistical evaluation. For example, if\nthe family of the distribution is known (e.g. Gaussian/Normal) and provided\nwith the `dist_class` parameter of `MCNullDist`, then permutation tests\ndone by `MCNullDist` allow determining the distribution parameters. Under the\n(strong) assumption of Gaussian distribution, 20-30 permutations should be\nsufficient to get sensible estimates of the distribution parameters.\n\n", 
            "But that would be another story..."
          ]
        }
      ], 
      "metadata": {}
    }
  ]
}