In tdeenes/eegR: Analyze EEG Signals

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 8, fig.height = 8)

The eegR package has been developed to analyze EEG (or MEG) signals. It provides all basic functionalities such as data import, filtering, ICA, channel interpolation, re-referencing, segmenting, time-frequency decomposition, complex visualization and statistical analysis of single-trial and averaged ERPs. However, instead of being a comprehensive tool, it has been developed to serve mainly as a general frame which is highly extendable by other, more specialized analytic methods.

ERP analysis

The eegR package can perform various analytic methods on averaged event-related potential (ERP) data:

point-to-point topographical ANOVA (see Koenig & Melie-Garcia, 2009)
point-to-point t-test and ANOVA with TFCE correction (see Smith & Nichols, 2009; Mensen & Khatami, 2013)
peak detection

The package also provides several convenience functions to handle complex ERP datasets. In the following, we demonstrate the use of the package on a real dataset.

library(eegR)
data(erps)
?erps

We will need the electrode positions and the group memberships, so we extract them.

id_dat <- attributes(erps)$id
chan_pos <- attributes(erps)$chan

Plot all channel curves for each condition/participant separately

First we will check the ERPs for each participant and condition separately. This way outliers can be spotted and excluded. As an example, plot the curves from the participant "03" and for the stimulus classes "A" and "B".

# subset erps 
dat03 <- subsetArray(erps, list(id = "03", stimclass = c("A", "B")))

# plot curves
plotERParray(dat03, lty = 1, col = "grey70", grid_dim = c(3, 2))

Plot GFP curves for each group and each condition

You can use lattice or ggplot2 to plot whatever aspect of your data. Before calling the graphing functions, you can use transformArray to arrange your data in appropriate format, and perform some preliminary computations.

Now we will compute the Global Field Power on the grand averages of dyslexic and control participants, and compare the curves.

# prepare data
tempdat <- transformArray(
    compGfp(amplitude) ~ . | id, 
    erps, list(readgroup = id_dat$group))

# create plot
library(ggplot2)
ggplot(tempdat, aes(x = time, y = amplitude, col = readgroup, lty = pairtype)) + 
    geom_line() + 
    facet_grid(stimclass ~ .)

Plot topographies

If you want to produce topoplots (2D maps of channel amplitudes), use plot2dview for single maps and complexplot2dview to compare topographies of different groups or conditions side-by-side.

Now we will take the topographies in the identical condition, and compare the reading groups and stimulus classes at 300 ms.

# subset the data (pairtype = "ident") and average across id based on 
# reading group membership
tempdat <- transformArray(y ~ . | id, erps, 
                          group = list(readgroup = id_dat$group), 
                          subset = list(pairtype = "ident"), 
                          datfr = FALSE)

# split the array along the readgroup and stimclass dimensions
tempdat <- splitArray(tempdat, c("readgroup", "stimclass"), drop = TRUE)

# compute GFP
gfpdat <- lapply(tempdat, compGfp)

# subset tempdat to include only the time samples at 300 ms
# use [] to keep the attributes of tempdat
plotdat <- tempdat
plotdat[] <- lapply(tempdat, subsetArray, list(time = "300"))

# plot the topographies and the centroids at 300 ms, and add the GFP curves 
# of the whole segments; do not display the channel positions and labels
complexplot2dview(
    plotdat, chan_pos, 300, 
    gfp = gfpdat, plot_ch = FALSE)

Note that using this approach, you can easily create an animation for the whole segment: just install the animation package and loop through all time points when calling complexplot2dview.

Perform a point-to-point ANOVA with TFCE correction

Threshold-free cluster enhancement (TFCE) is an excellent method to compare ERPs because it takes into account the dependencies of individual measurement points in the spatial and time domain.

To be able to use TFCE, you have to create a channel neigbourhood matrix. The package provides an interactive function (chanNb) for this purpose.

# first call chanNb without the alpha parameter
#chanNb(chan_pos)

# now that you know which alpha parameter gives good results, call the 
# function again with the given value
ChN <- chanNb(chan_pos, alpha = 0.7)

Now you can start the analysis. For serious purposes, set the number of permutations to 999 or higher (preferably 4999). Because the computation is resource intensive, use parallelization if possible. In the present example, we will compute bin averages (collapse 12 ms long time bins) to decrease the computational burden, and use only two CPU-cores.

# do bin averaging
tempdat <- avgBin(erps, "time", bin_length = 6)

# perform the analysis
results <- arrayAnova(
    tempdat, 
    factordef = list(between = "group",
                     w_id = "id", 
                     within = c("stimclass", "pairtype"),
                     observed = "group"),
    bwdat = id_dat, 
    perm = .(n = 99), 
    tfce = .(ChN = ChN),
    parallel = .(ncores = 2))

# plot the results (for now, only the p-values of the TFCE method)
modelplot(results, "p-value")

# print the summary
summary(results)