Last active
December 17, 2015 02:58
-
-
Save dengemann/5539403 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
""" |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment