[Nipy-devel] rm 2-way ANOVA

Christian Brodbeck christianmbrodbeck@gmail....
Mon Nov 16 19:45:27 CST 2009


>> I'd be interested to see how you went about doing this, if you think
>> it is general enough to contribute to the project.  Did you primarily
>> use Josef's code?
>> 
>> I wrote my own code, it depends very much on what your source data looks
>> like and how your factors are coded (as Josef said). It is fairly simple in
>> principle though and could be adapted. If in your data array subjects are
>> alternating inside blocks, the fixed effect of subject can be added to the
>> design matrix as:
>> 
>> (corrected version):
>>     _sub___   _conditions...
>> s1  1  0  0    1 ... \
>> s2  0  1  0    1      \ condition 1
>> s3  0  0  1    1      /
>> s4 -1 -1 -1    1     /
>> s1  1  0  0    0     \
>> s2  0  1  0    0      \ condition 2
>> ...
> 
> I don't think I understand the meaning of the -1 in s4.

I think it's merely a matter of how you set up contrasts between subjects: with -1s predictor variables contrast subjects against the last subject, whereas with 0 the mean for the last subject would be predicted by the intercept whereas the other s have their parameter. I tried substituting 0s for the -1s and got exactly the same results for my factors condition.  


> 
> If I understand this correctly than you have to loop 129 times ?

I am performing the ANOVA for each electrode for each sample point, so 129*150=19350 times


>> 
>> Yes, this is right.  Using a trunk of scipy for block_diag you have
>> 
>> import scikits.statsmodels as sm
>> import numpy as np
>> from scipy import linalg
>> 
>> data1 = sm.datasets.longley.Load()
>> data2 = sm.datasets.longley.Load()
>> data1.exog = sm.add_constant(data1.exog)
>> data2.exog = sm.add_constant(data2.exog)
>> exog = linalg.block_diag(data1.exog, data2.exog)
>> endog = np.vstack((data1.endog[:,None],data1.endog[:,None]))
>> res = sm.OLS(endog,exog).fit()
>> In [17]: res.params
>> Out[17]:
>> array([  1.50618723e+01,  -3.58191793e-02,  -2.02022980e+00,
>>        -1.03322687e+00,  -5.11041057e-02,   1.82915146e+03,
>>        -3.48225863e+06,   1.50618723e+01,  -3.58191793e-02,
>>        -2.02022980e+00,  -1.03322687e+00,  -5.11041057e-02,
>>         1.82915146e+03,  -3.48225863e+06])
>> 
>> res2 = sm.OLS(data1.endog, data1.exog).fit()
>> In [19]: res2.params
>> Out[19]:
>> array([  1.50618723e+01,  -3.58191793e-02,  -2.02022980e+00,
>>        -1.03322687e+00,  -5.11041057e-02,   1.82915146e+03,
>>        -3.48225863e+06])
>> 
>> So if you are using OLS to do the anova, then you can use the block
>> diagonal.  You should probably use sparse matrices though if your
>> design is very big.  There is an example of building a block diagonal
>> matrix in our sandbox (haven't had time to clean this up, write tests,
>> though the results are correct via the 'eyeball' test, or get a
>> polished results class yet).  It's in the __init__ method for the SUR
>> class.  http://bazaar.launchpad.net/~scipystats/statsmodels/trunk/annotate/head%3A/scikits/statsmodels/sandbox/sysreg.py
> 
> 
> For the case when all design matrices are the same, you could save on
> replicating them in the block diagonal matrix.
> OLS parameter estimation is easy for this case, linalg.leastsq allows
> for matrix endog, the work is in getting all the anova and regression
> result statistics.

Yes, all the design matrices are the same, so you are saying with the present sm.OLS I could estimate all anovas at once, using my model as exog and putting the different endogs into a certain matrix form?



Cheers,

--Christian
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/nipy-devel/attachments/20091117/eece3daf/attachment.html 


More information about the Nipy-devel mailing list