We will explain and apply in R
the
permutation-based cluster-Mass method proposed by Maris and
Oostenveld, 2007 and developed in the R
package
permuco4brain
by Frossard
and Renaud, 2018, using electroencephalography (EEG) data. The
cluster-mass is computed considering:
Finally, the All-Resolution Inference (ARI) from Rosenblatt et
al. 2018 is applied to compute the lower bound for the true
discovery proportion inside the clusters computed. Here, we will use the
ARIeeg
and hommel
R
packages.
First of all, you need to install and load the following packages:
#devtools::install_github("angeella/ARIeeg")
#devtools::install_github("bnicenboim/eeguana")
#devtools::install_github("jaromilfrossard/permuco4brain")
#devtools::install_github("jaromilfrossard/permuco")
#devtools::install_github("livioivl/eegusta")
library(ARIeeg) #to compute ARI for spatial-temporal clusters
library(eeguana) #to manage eeg data
library(eegusta) #to manage eeg data
library(permuco4brain) #to compute the spatial-temporal clusters
library(plotly)
library(tidyverse)
library(permuco) #to compute the temporal clusters
library(hommel) #to compute ARI for temporal clusters
library(abind)
The Dataset from the package ARIeeg
is an ERP
experiment composed by:
You can load it using:
load(system.file("extdata", "data_eeg_emotion.RData", package = "ARIeeg"))
We transform the data as eeg_lst
class object from the
R
package eeguana
using the function
eegUtils2eeguana
from the R
package
eegusta
:
dati_lst = eegUtils2eeguana(data = dati)
is_eeg_lst(dati_lst) #check
## [1] TRUE
and we drop off the final five channels:
chan_to_rm <- c("RM" , "EOGvo" ,"EOGvu"
, "EOGhl", "EOGhr")
dati_lst <-
dati_lst %>%
dplyr::select(-one_of(chan_to_rm))
Finally, we segment the data and select two conditions, i.e., disgust face(number \(3\)) and object (number \(5\)):
data_seg <- dati_lst %>%
eeg_segment(.description %in% c(3,5),
.lim = c(min(dati$timings$time), max(dati$timings$time))
) %>% eeg_baseline() %>%
mutate(
condition =
description
) %>%
dplyr::select(-c(type,description))
Some plots to understand the global mean difference between the two conditions:
A<-data_seg %>%
select(Fp1,Fp2, F3, F4) %>%
ggplot(aes(x = .time, y = .value)) +
geom_line(aes(group = condition)) +
stat_summary(
fun = "mean", geom = "line", alpha = 1, size = 1,
aes(color = condition),show.legend = TRUE
) +
facet_wrap(~.key) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_vline(xintercept = .17, linetype = "dotted") +
theme(legend.position = "bottom")+
scale_color_manual(labels = c("Disgust", "Object"), values = c("#80bfff", "#ff8080"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
A
if you want an interactive plot, you can use the function
ggplotly
from the package plotly
:
plotly::ggplotly(A)
The aim is to test if the difference in brain signal during the two conditions is different from \(0\) for each time point, i.e., \(500\). If the full set of channels is considered, we also have tests for each channel, i.e., \(27\), returning a total number of tests equals \(500 \cdot 27\). Therefore, we have \(500\) or \(500 \cdot 27\) statistical tests to perform at group-level, so considering the random subject effect. The multiple testing problem is obvious, and correction methods such as Bonferroni or similar do not capture the time(-spatial) correlation structure of the statistical tests; it will be a conservative method.
The cluster mass method is then used, proposed by Maris and Oostenveld, 2007. It is based on permutation theory; it gains some power with respect to other procedures correcting at the (spatial-)temporal cluster level instead of at the level of single tests. It is similar to the cluster mass approach in the fMRI framework, but in this case, the voxels, i.e., the single object of the analysis, are expressed in terms of time-points or combination time-points/channels. The method can then gain some power with respect to some traditional conservative FWER correction methods exploiting the (spatial-)temporal structure of the data.
The cluster mass method is based on the repeated measures ANOVA model, i.e.,
\[ y = 1_{N \times 1} \mu + X^{\eta} \eta + X^{\pi}\pi + X^{\eta \pi}\eta \pi + \epsilon \]
where \(1_{N \times 1}\) is a matrix with ones and
Therefore, \(y \sim (1\mu + X^{\eta} \eta, \Sigma)\), \(\pi \sim (0, \sigma^2_{\pi} I_{n_{subj}})\) and \(\eta \pi \sim (0,\text{cov}(\eta \pi))\). (N.B: \(\eta \pi\) is not the product between \(\eta\) and \(\pi\) but refers to the interaction effects between subjects and conditions).
We want to make inference on \(\eta\), such that \(H_0: \eta = 0\) vs \(H_1: \eta \ne 0\). We do that using the F statistic, i.e.,
\[ F = \dfrac{y^\top H_{X^{\eta}} y / (n_{stimuli} - 1)}{ y^\top H_{X^{\eta \pi}}y/(n_{stimuli} -1)(n_{subj} -1)} \] where \(H_{X}\) is the projection matrix, i.e., \(H_{X} = X(X^\top X)^{-1} X^\top\). In order to compute this test, we use an alternative definition of \(F\) based on the residuals:
\[ F_r = \dfrac{r^\top H_{X^{\eta}} r / (n_{stimuli} - 1)}{ r^\top H_{X^{\eta \pi}}r/(n_{stimuli} -1)(n_{subj} -1)} \]
where \(r = (H_{X^{\eta}} + H_{X^{\eta\pi}})y\). For further details, see Kherad Pajouh and Renaud, 2014.
So, let the group of permutation, including the identity transformation, \(\mathcal{P}\), we use \(r^\star = P r\), where \(P \in \mathcal{P}\) to compute the null distribution of our test, i.e., \(\mathcal{R}\), and then the p-value, i.e.,
\[ \text{p-value} = \dfrac{1}{B} \sum_{F^\star_r \in \mathcal{R}} \mathbb{I}(|F^\star_r| \ge |F_r|) \]
if the two-tailed is considered, where \(F^\star_r = f(r^\star)\).
We have this model for each time point \(t \in \{1, \dots, 500\}\) and each channel, so finally we will have \(n_{\text{time-points}} \times n_{\text{channels}}\) statistical tests/p-values (raw).
This method has been proposed by Maris and Oostenveld, 2007. It assumes that an effect will appear in clusters of adjacent time frames. Having statistics for each time point, we form these clusters using a threshold \(\tau\) as follows:
All contiguous time points with statistics above this threshold define a single cluster \(C_i\) with \(i \in \{1, \dots, n_C\}\), where \(n_C\) is the number of clusters found. For each time point in the same cluster \(C_i\), we assign the same cluster mass statistic \(m_i = f(C_i)\), where \(f\) is the function that summarizes the statistics of the entire cluster. Typically, it is the sum of the \(F\) statistics. The null distribution of the cluster mass \(\mathcal{M}\) is computed by iterating the above process for each permutation. The contribution of a permutation to the cluster-mass null distribution is the maximum overall cluster mass of that permutation. To check the significance of the cluster \(C_i\) of interest, we compare its cluster mass \(m_i = f(C_i)\) with the cluster mass null distribution \(\mathcal{M}\). Therefore, for each cluster \(C_i\), we have the associated p-values computed as
\[ p_i = \dfrac{1}{n_P} \sum_{m^\star \in \mathcal{M}} I\{m^\star \ge m_i\} \]
where \(m^\star \in \mathcal{M}\) is then calculated given permutation statistics. This method makes sense when analysing EEG data because if a difference in brain activity is thought to occur at time \(s\) for a given factor, then it is very likely that this difference will also occur at time \(s + 1\) (or \(s - 1\)).
In this case, we use graph theory, where the vertices represent the channels and the edges represent the adjacency relationships between two channels. The adjacency must be defined using prior information, so the three-dimensional Euclidean distance between channels is used. Two channels are defined as adjacent if their Euclidean distance is less than the threshold \(\delta\), where \(\delta\) is the smallest Euclidean distance that yields a connected graph Cheval, et al.,). This follows from the fact that there is no unconnected subgraph for a connected graph. The existence of subgraphs implies that some tests cannot be enter in the same cluster, which is not a useful assumption for the present analysis (Frossard and Renaud, 2018; Frossard, 2019).
Once we have a definition of spatial contiguity, we need to define temporal contiguity. We reproduce this graph \(n_{\text{time-points}}\) times, and we have edges between pairs of two vertices associated to the same electrode if they are temporally adjacent. The final graph has a total number of vertices, i.e., number of tests, equals (\(n_{\text{channels}} \times n_{\text{time-points}}\)). The following figure shows an example with \(64\) channels and \(3\) time measures:
We then delete all the vertices in which statistics are below a threshold, e.g., the \(95\) percentile of the null distribution of the \(F\) statistics. So, we have a new graph composed of multiple connected components, where each connected component defines the spatial-temporal cluster. We compute for each spatial-temporal cluster the cluster-mass statistic as before.
The cluster-mass null distribution is calculated using permutations that preserve the spatial-temporal correlation structure of the statistical tests, i.e., no changing the position of the electrodes and mixing the time points. We construct a three-dimensional array, where the first dimension represents the design of our experiments (subjects of \(\times\) stimuli), the second one the time points, and the third one the electrodes. So, we apply permutations only in the first dimension using the method proposed by Kherad Pajouh and Renaud, 2014.
In R, all of this is possible thanks to the permuco
and
permuco4brain
packages developed by Frossard
and Renaud, 2018.
So, we select one channel from our dataset, e.g. the Fp1:
Fp1 <- data_seg %>% select(Fp1)
signal_Fp1 <- Fp1%>%
signal_tbl()%>%
group_by(.id)%>%
nest()%>% #creates a list-column of 40 data frames having dim 500 x 2
mutate(data = map(data,~as.matrix(.x[-1])))%>% #drop off the first column of each dataframe, i.e., we take the column oof the signal for channel Fp1
pull(data)%>% #takes data created in mutate
invoke(abind,.,along = 2)%>% #we merge each dataframe in order to have one db having dimension 500 x 40
aperm(c(2,1)) #traspose
## Warning: `invoke()` was deprecated in purrr 1.0.0.
## ℹ Please use `exec()` instead.
dim(signal_Fp1)
## [1] 40 500
So, signal_Fp1
is a data frame that expresses the
channel Fp1 signals under two conditions in \(500\) time points for \(20\) subjects.
design <-
segments_tbl(Fp1)%>%
select(participant_id, condition)
dim(design)
## [1] 40 2
f <- signal_Fp1 ~ condition + Error(participant_id/(condition))
In the formula, we need to specify the Error(.)
term
since we are dealing with a repeated measures design. We specify a
subject-level random effect and a condition fixed effect nested within
subjects.
Thanks to the permuco
package, we can apply the temporal
cluster-Mass for the channel Fp1:
lm_Fp1 <- clusterlm(f,data = design)
summary(lm_Fp1)
## Effect: condition.
## Alternative Hypothesis: two.sided.
## Statistic: fisher(1, 38).
## Resampling Method: Rd_kheradPajouh_renaud.
## Type of Resampling: permutation.
## Number of Dependant Variables: 500.
## Number of Resamples: 5000.
## Multiple Comparisons Procedure: clustermass.
## Threshold: 4.098172.
## Mass Function: the sum.
## Table of clusters.
##
## start end cluster mass P(>mass)
## 1 47 48 10.49306 0.8418
## 2 170 210 424.79648 0.0288
## 3 217 225 55.56859 0.3328
## 4 378 381 21.22319 0.7046
Here we can see:
For example, considering the first significant cluster, we can compute the cluster mass as:
sum(lm_Fp1$multiple_comparison$condition$uncorrected$main[c(170:210), "statistic"])
## [1] 424.7965
We can also plot the temporal clusters:
plot(lm_Fp1)
The red dots represent the significant temporal cluster for the channel Fp1 composed by the time points from \(170\) to \(210\) using a threshold equals \(4.098\).
However, our significant cluster says only that at least one test is different from \(0\), we don’t know how many tests/time-points are significant (spatial specificity paradox). So, we can apply ARI to understand the lower bound of the number of true discovery proportion. The cluster is composed by the time points from \(170\) to \(210\), i.e., the size of the cluster is equal to \(41\).
praw <- lm_Fp1$multiple_comparison$condition$uncorrected$main[,2]
cluster <- c(170:210)
discoveries(hommel(praw), ix = cluster)/length(cluster)*100
## [1] 24.39024
Therefore, we have at least 24.39\(\%\) of true active time points in the cluster computed.
signal <-
data_seg%>%
signal_tbl()%>%
group_by(.id)%>%
nest()%>%
mutate(data = map(data,~as.matrix(.x[-1])))%>%
pull(data)%>%
invoke(abind,.,along = 3)%>%
aperm(c(3,1,2))
dim(signal)
## [1] 40 500 27
design <-
segments_tbl(data_seg)%>%
select(participant_id, condition)
dim(design)
## [1] 40 2
position_to_graph
from the permuco4brain
package:graph <- position_to_graph(channels_tbl(data_seg), name = .channel, delta = 53,
x = .x, y = .y, z = .z)
plot(graph)
f <- signal ~ condition + Error(participant_id/(condition))
Finally, run the main function:
model <- permuco4brain::brainperm(formula = f,
data = design,
graph = graph,
np = 1000,
multcomp = "clustermass",
return_distribution = TRUE)
## Computing Effect:
## 1 (condition) of 1. Start at 2023-02-19 17:05:08.
where np indicates the number of permutation.
Then, we can analyze the output:
print(model)
## Effect: condition.
## Alternative Hypothesis: two.sided.
## Statistic: fisher(1, 38).
## Resampling Method: Rd_kheradPajouh_renaud.
## Type of Resampling: permutation.
## Number of Dependant Variables: 13500.
## Number of Resamples: 1000.
## Multiple Comparisons Procedure: clustermass.
## Threshold: 4.098172.
## Mass Function: the sum.
## Table of clusters.
##
## Cluster id First sample Last sample N. chan. Main chan. Main chan. length
## 1 1 41 44 1 T7 4
## 2 2 47 48 1 Fp1 2
## 3 3 164 487 25 P8 298
## 4 4 260 265 1 P7 6
## 5 5 266 327 8 CPz 62
## 6 6 274 277 1 P7 4
## 7 7 285 500 8 P7 183
## 8 8 373 374 1 CP2 2
## 9 9 378 381 1 Fp1 4
## 10 10 498 500 1 P8 3
## N. test Clustermass P(>mass)
## 1 4 22.922758 1.000
## 2 2 10.493061 1.000
## 3 1322 15901.330605 0.001
## 4 6 25.774934 1.000
## 5 222 1085.206144 0.449
## 6 4 16.838841 1.000
## 7 666 5368.015395 0.037
## 8 2 8.321318 1.000
## 9 4 21.223185 1.000
## 10 3 12.861519 1.000
We have only two significant clusters. The first one is composed by \(25\) channels while the second one by \(8\) channels, with main channels P7. You can see in details the components of this cluster in
head(names(model$multiple_comparison$condition$clustermass$cluster$membership[which(as.vector(model$multiple_comparison$condition$clustermass$cluster$membership)==3)]))
## [1] "P7_164" "Pz_164" "P7_165" "Pz_165" "P7_166" "CPz_166"
You can see the significant cluster (in red) at fixed time points (e.g. 300) using plot:
plot(model, samples = 300)
and the significant cluster over time and over channels using:
image(model)
where the significant clusters are represented in a colour-scale and the non-significant one in grey. The white pixels are tests which statistic are below the threshold.
However, our significant clusters say only that at least one combination channels/time-points is different from \(0\), we do not know how many combinations are significant (spatial specificity paradox). So, we can apply ARI to understand the lower bound of the number of true discovery proportion:
ARIeeg::ARIeeg(model = model)
## ID Total clustermass pvalue False Null True Null Active Proportion
## [1,] 1 4 22.922758 1.000 0 4 0
## [2,] 2 2 10.493061 1.000 0 2 0
## [3,] 3 1322 15901.330605 0.001 0 1322 0
## [4,] 4 6 25.774934 1.000 0 6 0
## [5,] 5 222 1085.206144 0.449 0 222 0
## [6,] 6 4 16.838841 1.000 0 4 0
## [7,] 7 666 5368.015395 0.037 0 666 0
## [8,] 8 2 8.321318 1.000 0 2 0
## [9,] 9 4 21.223185 1.000 0 4 0
## [10,] 10 3 12.861519 1.000 0 3 0
Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of EEG-and MEG-data. Journal of neuroscience methods, 164(1), 177-190.
Kherad-Pajouh, S., & Renaud, O. (2015). A general permutation approach for analyzing repeated measures ANOVA and mixed-model designs. Statistical Papers, 56(4), 947-967.
Frossard, J. (2019). Permutation tests and multiple comparisons in the linear models and mixed linear models, with extension to experiments using electroencephalography. DOI: 10.13097/archive-ouverte/unige:125617.
Frossard, J. & O. Renaud (2018). Permuco: Permutation Tests for Regression, (Repeated Measures) ANOVA/ANCOVA and Comparison of Signals. R Packages.
Cheval, Boris, et al. “Avoiding sedentary behaviors requires more cortical resources than avoiding physical activity: An EEG study.” Neuropsychologia 119 (2018): 68-80.