[pymvpa] basic GLM

[Michael Hanke]

Hi,

On Tue, Jun 09, 2009 at 06:09:29PM -0700, Daqiang Sun wrote:
> 
> Hi,
> 
> Could anybody show me the right way to run GLM?
> I have a dataset and a design array (mydesign) as:
> array([[1, 1, 1, 1, 1, ..., 0, 0, 0, 0, 0],
>        [0, 1, 2, 3, 4, 5, ..., 6, 7, 8, 9]])
> 
> I tried:
> 
> myglm = measures.glm.GLM(N.transpose(mydesign), voi='pe')
> myglm(dataset)
> 
> but the the number of elements in output array is the same as feature number, 
> not nfeatures * nregressors. I tried more designs with more regressors too, 
> but got the same issue. I'm not sure what's going wrong there.

The reason for this is a default combiner that collapses the betas
across all regressors into a single score -- which doesn't make much
sense here (Note that this is a candidate for change in a future major
release, but we want to accumulate all the changes that change the
current behavior (in maybe backward-incompatible way) into a single
release with a big fat warning sign). If you change the above line to

  myglm = measures.glm.GLM(N.transpose(mydesign), voi='pe', combiner=None)

you should get what you desire -- a matrix (nfeatures x nregressors).

> Also should the design array have dummy coding for the intercept?

Yes! This implementation is very simple, no automatic demeaning is done. If you
expect an intercept it needs a constant column in the design matrix.

If you take a look at mvpa/tests/test_stats.py you can find a more
elaborate example at the end...

> I tried scipy.stats.glm by the way, but it seems to me that it only deals with 
> dichotomous regressors. 
> 
> I just found that rather than resorting to other packages, running non-classification 
> things within is such a convenience. I really appreciate your extra efforts to include 
> so many useful features in the package.

That is true. However, adding that kind of special stuff always implies
the risk of having bugs in it that are never detected since only very
few people use it. With the rise of NiPy we plan to remove our custom
implementations and integrate with NiPy to allow for a seamless
combination of the available functionality.

I'm very grateful that you are testing this code. However, keep being
suspicious, since you are probably one of the first people to touch it!

Thanks,

Michael

-- 
GPG key:  1024D/3144BE0F Michael Hanke
http://apsy.gse.uni-magdeburg.de/hanke
ICQ: 48230050