test_mne_anova_against_r_3way_1way.py

import numpy as np
import pandas as pd
from rpy2 import robjects as rob

# test data as found in MNE-Python
data = np.random.RandomState(42).randn(20, 8)

# prepare categorial factors for R long table format
var = data.flatten()
aa = ['a1', 'a1', 'a1', 'a1', 'a2', 'a2', 'a2', 'a2'] * 20
bb = ['b1', 'b1', 'b2', 'b2', 'b1', 'b1', 'b2', 'b2'] * 20
cc = ['c1', 'c2', 'c1', 'c2', 'c1', 'c2', 'c1', 'c2'] * 20
dd = ['d1', 'd2', 'd3', 'd4', 'd5', 'd6', 'd7', 'd8'] * 20
subject = np.repeat(np.arange(20), 8).astype('S2')


# compose data frame and save in file
df = pd.DataFrame(dict(var=var, A=aa, B=bb, C=cc, D=dd, subject=subject))
df.to_csv('test_r_mne_anova.csv', index=False)

# quick and dirty setup R
rob.r("""
df <- read.csv('test_r_mne_anova.csv')
df$subject <- as.factor(df$subject)
options(digits=20)
""")

# our repeated measures anova formula for 3 ways
formula = 'var ~ A * B * C + Error(subject / (A * B * C))'

# now print results
print rob.r('(summary(ano <- aov(%s, df)))' % formula)

""" Result (8 /5/2013)
Error: Subject
          Df             Sum Sq             Mean Sq F value Pr(>F)
Residuals 19 10.842172197658172 0.57064064198200903

Error: Subject:A
          Df              Sum Sq             Mean Sq            F value  Pr(>F)
A          1  2.1505188758333156 2.15051887583331558 2.5676199999999998 0.12557
Residuals 19 15.9135046176327215 0.83755287461224848

Error: Subject:B
          Df               Sum Sq             Mean Sq F value  Pr(>F)
B          2  0.40948309086464102 0.20474154543232040 0.24006 0.78776
Residuals 38 32.40868450052717975 0.85286011843492571

Error: Subject:A:B
          Df              Sum Sq             Mean Sq            F value Pr(>F)
A:B        2  3.3986864423342231 1.69934322116711156 1.7563800000000001 0.1864
Residuals 38 36.7660584898879179 0.96752785499705052
"""


""" Result (27 /5/2015)
Error: subject
          Df             Sum Sq             Mean Sq F value Pr(>F)
Residuals 19 16.718626365490039 0.87992770344684412

Error: subject:A
          Df               Sum Sq             Mean Sq             F value
A          1  0.59677877912566002 0.59677877912566002 0.74783999999999995
Residuals 19 15.16206330175759653 0.79800333167145243
           Pr(>F)
A         0.39795
Residuals

Error: subject:B
          Df              Sum Sq             Mean Sq F value  Pr(>F)
B          1  0.1967705112514305 0.19677051125143052 0.20895 0.65278
Residuals 19 17.8928578387976600 0.94172935993671891

Error: subject:C
          Df               Sum Sq             Mean Sq             F value
C          1  0.62588369501616603 0.62588369501616603 0.99404000000000003
Residuals 19 11.96303588242783000 0.62963346749620153
           Pr(>F)
C         0.33128
Residuals

Error: subject:A:B
          Df               Sum Sq             Mean Sq F value  Pr(>F)
A:B        1  0.23625737084348761 0.23625737084348764 0.21378 0.64907
Residuals 19 20.99771363935110102 1.10514282312374212

Error: subject:A:C
          Df               Sum Sq             Mean Sq              F value
A:C        1  0.13534296835637619 0.13534296835637624 0.094039999999999999
Residuals 19 27.34417366322780296 1.43916703490672648
           Pr(>F)
A:C       0.76243
Residuals

Error: subject:B:C
          Df               Sum Sq              Mean Sq F value  Pr(>F)
B:C        1 0.055428413815978697 0.055428413815978662 0.11685 0.73622
Residuals 19 9.012490297394636585 0.474341594599717697

Error: subject:A:B:C
          Df              Sum Sq             Mean Sq F value Pr(>F)
A:B:C      1  2.3684998313995598 2.36849983139955977 2.78749 0.1114
Residuals 19 16.1440923060018058 0.84968906873693717
"""

# our repeated measures anova formula for 1 way
formula = 'var ~ D + Error(subject / (D))'

# now print results
print rob.r('(summary(ano <- aov(%s, df)))' % formula)

"""
Error: subject
          Df             Sum Sq             Mean Sq F value Pr(>F)
Residuals 19 16.718626365490039 0.87992770344684412

Error: subject:D
           Df               Sum Sq             Mean Sq             F value
D           7   4.2149615698086684 0.60213736711552401 0.67571999999999999
Residuals 133 118.5164269289584382 0.89110095435307102
           Pr(>F)
D         0.69232
Residuals
"""