-
Notifications
You must be signed in to change notification settings - Fork 1
/
sensor_anova_Baleen_gk_peripheral_acsc.R
172 lines (110 loc) · 5.24 KB
/
sensor_anova_Baleen_gk_peripheral_acsc.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
sensor_anova_Baleen_gk_peripheral_acsc <-function(acfilePrefix,scfilePrefix,t1,t2){
library('ez')
options(digits=2)
filePath <- "/cluster/kuperberg/SemPrMM/MEG/results/sensor_level/R/"
load(paste(filePath,acfilePrefix,".BaleenLP_All.gk_peripheral.",t1,"-",t2,".df",sep=""))
erpData.LP.ac.all <- erpData.all
load(paste(filePath,acfilePrefix,".BaleenHP_All.gk_peripheral.",t1,"-",t2,".df",sep=""))
erpData.HP.ac.all <- erpData.all
load(paste(filePath,scfilePrefix,".BaleenLP_All.gk_peripheral.",t1,"-",t2,".df",sep=""))
erpData.LP.sc.all <- erpData.all
load(paste(filePath,scfilePrefix,".BaleenHP_All.gk_peripheral.",t1,"-",t2,".df",sep=""))
erpData.HP.sc.all <- erpData.all
#combine LP and HP datasets for SC and AC
erpData.all <-rbind(erpData.LP.ac.all,erpData.HP.ac.all,erpData.LP.sc.all,erpData.HP.sc.all)
#more transparent factor labels
colnames(erpData.all)[3] <- "prop"
colnames(erpData.all)[4] <- "prime"
#get rid of bad channels T9, T10, TP9 and TP10
erpData.all <- subset(erpData.all,elec !=29 & elec != 39 & elec !=40 & elec !=50)
outfile <-paste(filePath,"acsc.Baleen_All.",t1,"-",t2,"_anova_gk_peripheral.txt",sep="")
sink(outfile);
################################################################
#############QUADRANT ANALYSIS##################################
################################################################
#Just get those electrodes that are in the defined groups
erpData.quad <-subset(erpData.all, hemCode != 'XX')
erpData.quad <- droplevels(erpData.quad)
#####Omnibus ANOVA######
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad,dv = .(amp),wid=.(subj),within=.(prime,prop,hemCode,antCode),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("All")
print(eztest)
#####Within Low Proportion ANOVA####
erpData.quad.lo <- subset(erpData.quad, prop == 'BaleenLP_All')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.lo,dv = .(amp),wid=.(subj),within=.(prime,hemCode,antCode),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("LowProp - Related vs. Unrelated")
print(eztest)
#####Within High Proportion ANOVA####
erpData.quad.hi <- subset(erpData.quad, prop == 'BaleenHP_All')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.hi,dv = .(amp),wid=.(subj),within=.(prime,hemCode,antCode),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("HighProp - Related vs. Unrelated")
print(eztest)
#####Within Anterior ANOVA####
erpData.quad.frontal <- subset(erpData.quad, antCode == 'frontal')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.frontal,dv = .(amp),wid=.(subj),within=.(prime,prop,hemCode),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("frontal")
print(eztest)
#####Within Posterior ANOVA####
erpData.quad.parietal <- subset(erpData.quad, antCode == 'parietal')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.parietal,dv = .(amp),wid=.(subj),within=.(prime,prop,hemCode),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("parietal")
print(eztest)
#####Within Left ANOVA####
erpData.quad.left <- subset(erpData.quad, hemCode == 'left')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.left,dv = .(amp),wid=.(subj),within=.(prime,prop,antCode),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("Left")
print(eztest)
#####Within Right ANOVA####
erpData.quad.right <- subset(erpData.quad, hemCode == 'right')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.right,dv = .(amp),wid=.(subj),within=.(prime,prop,antCode),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("Right")
print(eztest)
#####Within Left, Within Anterior ANOVA####
erpData.quad.left.frontal <- subset(erpData.quad, antCode == 'frontal' & hemCode == 'left')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.left.frontal,dv = .(amp),wid=.(subj),within=.(prime,prop),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("Left frontal")
print(eztest)
#####Within Left, Within Posterior ANOVA####
erpData.quad.left.parietal <- subset(erpData.quad, antCode == 'parietal' & hemCode == 'left')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.left.parietal,dv = .(amp),wid=.(subj),within=.(prime,prop),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("Left parietal")
print(eztest)
#####Within Right, Within Anterior ANOVA####
erpData.quad.right.frontal <- subset(erpData.quad, antCode == 'frontal' & hemCode == 'right')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.right.frontal,dv = .(amp),wid=.(subj),within=.(prime,prop),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("Right frontal")
print(eztest)
#####Within Right, Within Posterior ANOVA####
erpData.quad.right.parietal <- subset(erpData.quad, antCode == 'parietal' & hemCode == 'right')
#compute the ANOVA
eztest <- ezANOVA(data=erpData.quad.right.parietal,dv = .(amp),wid=.(subj),within=.(prime,prop),between=.(subjGroup),type=3,detailed=TRUE)
#print results to file
print("Right parietal")
print(eztest)
########
#Print the marginal means
erpData.quad.aov <- aov(amp ~ prime * prop * hemCode * antCode + Error(subj/(prime * prop * hemCode * antCode)),data=erpData.quad)
print("Marginal Means")
print(model.tables(erpData.quad.aov,"means"),digits=5)
sink()
}