Analyzing EEG time-frequency data in a factorial design I am trying to find the best strategy to analyze a set of EEG time-frequency data from 24 subjects, using 64 electrodes. This is an exploratory study as I have no a-priori hypothesis on where or when changes in two frequency bands might be observed. In this case I have found that using non-parametrical statistical testing might be the best soulution: in particular, cluster-based permutation testing is widely used in cognitive electrophysiology.

However this approach is not feasible as I have a 2X2X2 within subject design, and it is not possible to look at interaction with permutation testing. On the other hand, running an ANOVA on the whole data-set might well result in an high number of false positives.

A solution to my proble might be to define ROI(s), so to perform ANOVA(s) on a limited number of electrodes within a time window. However I have no clue on how to define such ROI(s).

What is the best strategy according to you?
It is possible to look at within-subject interactions using cluster permutation tests. 
I.e. Eric Maris here: https://mailman.science.ru.nl/pipermail/fieldtrip/2011-January/003447.html

2x2x2 with factors A,B,C has 4 interactions AxB, AxC, BxC and AxBxC and 8 cells

+--------+-------------+-------------+
| C = 0  |    B = 0    |    B = 1    |
+--------+-------------+-------------+
| A = 0  | A0 B0 C0    | A0 B1 C0    |
| A = 1  | A1 B0 C0    | A1 B1 C0    |
+--------+-------------+-------------+

+--------+-------------+-------------+
| C = 1  |    B = 0    |    B = 1    |
+--------+-------------+-------------+
| A = 0  | A0 B0 C1    |  A0 B1 C1   |
| A = 1  | A1 B0 C1    |  A1 B1 C1   |
+--------+-------------+-------------+



2-way: You calculate the difference of differences at the average of missing factor:

$$
interaction_{AB} = \frac{interaction_{AB_{c=0}}  + interaction_{AB_{c=1}}}{2}  = \frac{(A_0B_0C_0-A_1B_0C_0)- (A_0B_1C_0-A_1B_1C_0) + (A_0B_0C_1-A_1B_0C_1)- (A_0B_1C_1-A_1B_1C_1)}{2}
$$

and for the three way interaction is the difference of interactions
$$
interaction_{AB_{c=1}} - interaction_{AB_{c=0}}
$$
(or any of the other pairs, same result)

I wrote a small script that confirms this:




library(tidyverse)
library("ANOVApower")

# simulate simple 2x2x2 design
design_result &lt;- ANOVA_design(design = "2w*2w*2w",
                              n = 10, 
                              mu = c(1:8), 
                              sd = 1.0, 
                              labelnames = c("A",0,1,"B",0,1,"C",0,1))

d = design_result$dataframe
contrasts(d$A) = c(-.5,.5)
contrasts(d$B) = c(-.5,.5)
contrasts(d$C) = c(-.5,.5)

# Cool, the ANOVApower script already gives us a unique column for each cell
# we can forget about subject because everything is balanced :)
d_cell = d%&gt;%group_by(cond) %&gt;%
  summarise(y = mean(y))


d_interactions = d_cell %&gt;% summarise(
                     AB_C0 = (y[cond=='A_0_B_0_C_0']-y[cond=='A_0_B_1_C_0']) - (y[cond=='A_1_B_0_C_0']-y[cond=='A_1_B_1_C_0']),
                     AB_C1 = (y[cond=='A_0_B_0_C_1']-y[cond=='A_0_B_1_C_1']) - (y[cond=='A_1_B_0_C_1']-y[cond=='A_1_B_1_C_1']),
                     AC_B0 = (y[cond=='A_0_B_0_C_0']-y[cond=='A_0_B_0_C_1']) - (y[cond=='A_1_B_0_C_0']-y[cond=='A_1_B_0_C_1']),
                     AC_B1 = (y[cond=='A_0_B_1_C_0']-y[cond=='A_0_B_1_C_1']) - (y[cond=='A_1_B_1_C_0']-y[cond=='A_1_B_1_C_1']),
                     CB_A0 = (y[cond=='A_0_B_0_C_0']-y[cond=='A_0_B_1_C_0']) - (y[cond=='A_0_B_0_C_1']-y[cond=='A_0_B_1_C_1']),
                     CB_A1 = (y[cond=='A_1_B_0_C_0']-y[cond=='A_1_B_1_C_0']) - (y[cond=='A_1_B_0_C_1']-y[cond=='A_1_B_1_C_1'])
)

d_2way = d_interactions %&gt;% summarise(AB = (AB_C0 + AB_C1)/2,
                             AC = (AC_B0 + AC_B1)/2,
                             BC = (CB_A0 + CB_A1)/2)


d_interactions$AB_C0 - d_interactions$AB_C1
d_interactions$AB_C1 - d_interactions$AB_C0

# confirm with linear model results
lm(y~A*B*C,d)



Which returns:

&gt; d_2way
# A tibble: 1 x 3
      AB     AC    BC
   &lt;dbl&gt;  &lt;dbl&gt; &lt;dbl&gt;
1 -0.683 -0.827 0.129


&gt; d_interactions$AB_C1 - d_interactions$AB_C0
[1] 1.307831
&gt; d_interactions$AB_C1 - d_interactions$AB_C0
[1] 1.307831


&gt; lm(y~A*B*C,d)

Call:
lm(formula = y ~ A * B * C, data = d)

Coefficients:
(Intercept)           A1           B1           C1        A1:B1        A1:C1        B1:C1     A1:B1:C1  
     4.3853       3.8473       1.8056       1.1670      -0.6832      -0.8272       0.1288       1.3078  


Identical results. 

This is what I would do