Package mvpa :: Package measures :: Module glm

[hide private]
[frames] | no frames]

# 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