Source Code for Module mvpa.measures.glm

  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  """The general linear model (GLM).""" 
 10   
 11  __docformat__ = 'restructuredtext' 
 12   
 13  import numpy as N 
 14   
 15  from mvpa.measures.base import FeaturewiseDatasetMeasure 
 16  from mvpa.misc.state import StateVariable 
 17   
 18 -class GLM(FeaturewiseDatasetMeasure): 
 19      """General linear model (GLM). 
 20   
 21      Regressors can be defined in a design matrix and a linear fit of the data 
 22      is computed univariately (i.e. indepently for each feature). This measure 
 23      can report 'raw' parameter estimates (i.e. beta weights) of the linear 
 24      model, as well as standardized parameters (z-stat) using an ordinary 
 25      least squares (aka fixed-effects) approach to estimate the parameter 
 26      estimate. 
 27   
 28      The measure is reported in a (nfeatures x nregressors)-shaped array. 
 29      """ 
 30   
 31      pe = StateVariable(enabled=False, 
 32          doc="Parameter estimates (nfeatures x nparameters).") 
 33   
 34      zstat = StateVariable(enabled=False, 
 35          doc="Standardized parameter estimates (nfeatures x nparameters).") 
 36   
 37 -    def __init__(self, design, voi='pe', **kwargs): 
 38          """ 
 39          :Parameters: 
 40            design: array(nsamples x nregressors) 
 41              GLM design matrix. 
 42            voi: 'pe' | 'zstat' 
 43              Variable of interest that should be reported as feature-wise 
 44              measure. 'beta' are the parameter estimates and 'zstat' returns 
 45              standardized parameter estimates. 
 46          """ 
 47          FeaturewiseDatasetMeasure.__init__(self, **kwargs) 
 48          # store the design matrix as a such (no copying if already array) 
 49          self._design = N.asmatrix(design) 
 50   
 51          # what should be computed ('variable of interest') 
 52          if not voi in ['pe', 'zstat']: 
 53              raise ValueError, \ 
 54                    "Unknown variable of interest '%s'" % str(voi) 
 55          self._voi = voi 
 56   
 57          # will store the precomputed Moore-Penrose pseudo-inverse of the 
 58          # design matrix (lazy calculation) 
 59          self._inv_design = None 
 60          # also store the inverse of the inner product for beta variance 
 61          # estimation 
 62          self._inv_ip = None 
 63   
 64   
 65 -    def _call(self, dataset): 
 66          # just for the beauty of it 
 67          X = self._design 
 68   
 69          # precompute transformation is not yet done 
 70          if self._inv_design is None: 
 71              self._inv_ip = (X.T * X).I 
 72              self._inv_design = self._inv_ip * X.T 
 73   
 74          # get parameter estimations for all features at once 
 75          # (betas x features) 
 76          betas = self._inv_design * dataset.samples 
 77   
 78          # charge state 
 79          self.states.pe = pe = betas.T.A 
 80   
 81          # if betas and no z-stats are desired return them right away 
 82          if self._voi == 'pe' and not self.states.isEnabled('zstat'): 
 83              # return as (feature x beta) 
 84              return pe 
 85   
 86          # compute residuals 
 87          residuals = X * betas 
 88          residuals -= dataset.samples 
 89   
 90          # estimates of the parameter variance and compute zstats 
 91          # assumption of mean(E) == 0 and equal variance 
 92          # XXX next lines ignore off-diagonal elements and hence covariance 
 93          # between regressors. The humble being writing these lines asks the 
 94          # god of statistics for forgives, because it knows not what it does 
 95          diag_ip = N.diag(self._inv_ip) 
 96          # (features x betas) 
 97          beta_vars = N.array([ r.var() * diag_ip for r in residuals.T ]) 
 98          # (parameter x feature) 
 99          zstat = pe / N.sqrt(beta_vars) 
100   
101          # charge state 
102          self.states.zstat = zstat 
103   
104          if self._voi == 'pe': 
105              # return as (feature x beta) 
106              return pe 
107          elif self._voi == 'zstat': 
108              # return as (feature x zstat) 
109              return zstat 
110   
111          # we shall never get to this point 
112          raise ValueError, \ 
113                "Unknown variable of interest '%s'" % str(self._voi) 
114