Source Code for Module mvpa.measures.anova

  1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 
  2  # vi: set ft=python sts=4 ts=4 sw=4 et: 
  3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  4  # 
  5  #   See COPYING file distributed along with the PyMVPA package for the 
  6  #   copyright and license terms. 
  7  # 
  8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  9  """FeaturewiseDatasetMeasure performing a univariate ANOVA.""" 
 10   
 11  __docformat__ = 'restructuredtext' 
 12   
 13  import numpy as N 
 14   
 15  from mvpa.measures.base import FeaturewiseDatasetMeasure 
 16   
 17  # TODO: Extend with access to functionality from scipy.stats? 
 18  # For binary: 
 19  #  2-sample kolmogorov-smirnof might be interesting 
 20  #   (scipy.stats.ks_2samp) to judge if two conditions are derived 
 21  #   from different distributions (take it as 'activity' vs 'rest'), 
 22  # 
 23  # For binary+multiclass: 
 24  #  kruskal-wallis H-test (scipy.stats.kruskal) 
 25  # 
 26  # and may be some others 
 27   
 28 -class OneWayAnova(FeaturewiseDatasetMeasure): 
 29      """`FeaturewiseDatasetMeasure` that performs a univariate ANOVA. 
 30   
 31      F-scores are computed for each feature as the standard fraction of between 
 32      and within group variances. Groups are defined by samples with unique 
 33      labels. 
 34   
 35      No statistical testing is performed, but raw F-scores are returned as a 
 36      sensitivity map. As usual F-scores have a range of [0,inf] with greater 
 37      values indicating higher sensitivity. 
 38      """ 
 39   
 40 -    def _call(self, dataset, labels=None): 
 41          """Computes featurewise f-scores using simple comparisons.""" 
 42          # group means 
 43          means = [] 
 44          # with group variance 
 45          vars_ = [] 
 46   
 47          # if no labels were provided -- use the ones within the dataset 
 48          if labels is None: 
 49              labels = dataset.labels 
 50   
 51          # split by groups -> [groups x [samples x features]] 
 52          for ul in N.unique(labels): 
 53              ul_samples = dataset.samples[labels == ul] 
 54              means.append(ul_samples.mean(axis=0)) 
 55              vars_.append(ul_samples.var(axis=0)) 
 56   
 57          # mean of within group variances 
 58          mvw = N.array(vars_).mean(axis=0) 
 59          # variance of group means 
 60          vgm = N.array(means).var(axis=0) 
 61   
 62          # compute f-scores (in-place to save some cycles) 
 63          # XXX may cause problems when there are features with no variance in 
 64          # some groups. One could deal with them here and possibly assign a 
 65          # zero f-score to throw them out, but at least theoretically zero 
 66          # variance is possible. Another possiblilty could be to apply 
 67          # N.nan_to_num(), but this might hide the problem. 
 68          # Michael therefore thinks that it is best to let the user deal with 
 69          # it prior to any analysis. 
 70   
 71          # for features where there is no variance between the groups, 
 72          # we should simply leave 0 as is, and avoid that way NaNs for 
 73          # invariance features 
 74          vgm0 = vgm.nonzero() 
 75          vgm[vgm0] /= mvw[vgm0] 
 76   
 77          return vgm 
 78   
 79   
 80 -class CompoundOneWayAnova(OneWayAnova): 
 81      """Compound comparisons via univariate ANOVA. 
 82   
 83      Provides F-scores per each label if compared to the other labels. 
 84      """ 
 85   
 86 -    def _call(self, dataset): 
 87          """Computes featurewise f-scores using compound comparisons.""" 
 88   
 89          orig_labels = dataset.labels 
 90          labels = orig_labels.copy() 
 91   
 92          results = [] 
 93          for ul in dataset.uniquelabels: 
 94              labels[orig_labels == ul] = 1 
 95              labels[orig_labels != ul] = 2 
 96              results.append(OneWayAnova._call(self, dataset, labels)) 
 97   
 98          # features x labels 
 99          return N.array(results).T 
100