Introduction - Preprocessing and visualization of EEG data
Bruno Nicenboim
2022-10-10
Source:vignettes/intro.Rmd
intro.RmdThe package eeguana provides a framework for doing simple pre-processing with specialized functions and manipulating EEG data with dplyr verbs (e.g., eeg_mutate, eeg_filter, eeg_summarize) extended to a new class eeg_lst, and ggplot wrapper functions. The new class is inspired by tidyverse principles but it’s not really “tidy” (due to space considerations), it’s a list of (i) a wide table that contains the signal amplitudes at every sample point of the EEG, (ii) an events table with information about markers (or triggers), blinks and other exported information, and (iii) a long table with experimental information, such as participant (recording), conditions, etc.
While it’s possible to transform the eeg_lst to a data.frame, data.table or a tibble (with as.data.frame(), as.data.table() and as_tibble()), the motivation for manipulating the data in the eeg_lst format has to do with size considerations. In this case, the original file was 113 MB, converting it to a long format entails a lot of repetition and generates an object of 556 MB. While this will still work here, a long format quickly becomes prohibitive in real settings with longer recordings.
Working with eeguana
Dplyr-like verbs always return an eeg_lst object, which allows us to use magrittr’s pipe, %>% (see ?`dplyr-eeguana`) or the new R pipe |>. In addition, eeg_/ch_ functions will also return an eeg_lst unless they have a suffix _tbl that indicates they return a data frame. In general, we will work with the eeg_lst, unless we want to modify the channels information (with channels_tbl()) or the events table containing markers, artifacts ans custom annotations (with events_tbl).
A practical example: the N170 effect
Here, I exemplify the use of eeguana with raw EEG data exported in the format of BrainVision 1.0. The data belong to a simple experiment where a participant was presented 100 faces and 100 assorted images in random order. The task of the experiment was to mentally count the number of faces.
First we download the data:
# Run the following or just download the files from raw_faces folder in https://osf.io/tbwvz/
library(httr)
GET("https://osf.io/wv5hp//?action=download",
write_disk("./s1_faces.vhdr", overwrite = TRUE),
progress()
)
GET("https://osf.io/c8wea//?action=download",
write_disk("./s1_faces.vmrk", overwrite = TRUE),
progress()
)
GET("https://osf.io/uhsde//?action=download",
write_disk("./s1_faces.eeg", overwrite = TRUE),
progress()
)BrainVision 1.0 exports three files: s1_faces.vhdr, s1_faces.vmrk, and s1_faces.eeg. The file s1_faces.vhdr contains the metadata and links to the other two files, s1_faces.vmrk contains the triggers and other events in the samples, and s1_faces.eeg contains the signals at every sample for every channel recorded.
library(dplyr)
library(ggplot2)
library(stringr)
library(eeguana)
set.seed(123) # ICA will always find the same componentsWe first need to read the data:
faces <- read_vhdr("s1_faces.vhdr")
#> # Reading file s1_faces.vhdr...
#> # Data from ./s1_faces.eeg was read.
#> # Data from 1 segment(s) and 34 channels was loaded.
#> # Object size in memory 113.4 MbThe function read_vhdr() creates a list with data frames for the signal, events, segments information, and incorporates in its attributes generic EEG information.
faces
#> # EEG data:
#>
#> # Signal table:
#> .id .sample Fp1 Fpz Fp2 F7 F3
#> 1: 1 1 -17878.56 -8921.395 -14028.81 -284.2541 -1378.1467
#> 2: 1 2 -17903.75 -8948.022 -14056.53 -303.9134 -1400.6560
#> 3: 1 3 -17916.55 -8961.906 -14070.91 -307.1683 -1407.8466
#> 4: 1 4 -17918.17 -8963.691 -14068.63 -299.8859 -1401.9985
#> 5: 1 5 -17917.71 -8965.658 -14069.16 -294.7915 -1401.2996
#> ---
#> 424484: 1 424484 -16094.71 -8234.730 -13505.21 609.9043 -236.9552
#> 424485: 1 424485 -16140.11 -8280.301 -13547.56 603.1918 -311.0485
#> 424486: 1 424486 -16145.24 -8275.078 -13528.22 541.3096 -291.0218
#> 424487: 1 424487 -16165.22 -8294.479 -13546.50 467.7130 -302.0740
#> 424488: 1 424488 -16167.39 -8305.199 -13561.12 576.3062 -319.2688
#> Fz F4 F8 FC5 FC1 FC2 FC6
#> 1: -8077.128 1167.489 -1663.4490 -2315.080 -9040.689 -795.0549 -3665.146
#> 2: -8101.660 1145.292 -1685.9034 -2338.859 -9064.137 -816.1298 -3688.133
#> 3: -8114.331 1138.028 -1691.4203 -2343.475 -9071.438 -823.4675 -3690.046
#> 4: -8116.114 1139.702 -1683.0713 -2333.746 -9065.608 -817.9690 -3681.642
#> 5: -8119.442 1136.373 -1684.2664 -2326.648 -9061.397 -817.2882 -3679.563
#> ---
#> 424484: -7991.669 1977.146 -162.8618 -2290.511 -9349.164 -318.8641 -2805.395
#> 424485: -8045.019 1938.122 -183.0542 -2424.133 -9410.679 -364.1038 -2833.090
#> 424486: -8032.237 1965.358 -132.3712 -2291.431 -9379.728 -342.3297 -2811.445
#> 424487: -8052.393 1930.877 -180.4978 -2361.589 -9381.456 -358.2924 -2868.325
#> 424488: -8065.579 1927.474 -188.0929 -2368.595 -9409.171 -363.8277 -2853.797
#> M1 T7 C3 Cz C4 T8 M2
#> 1: 0 -5156.998 -2466.320 -1568.4097 -4427.577 -10351.511 -1845.363
#> 2: 0 -5184.417 -2490.411 -1591.3053 -4450.049 -10372.402 -1850.862
#> 3: 0 -5186.146 -2495.726 -1597.7970 -4454.775 -10378.323 -1845.087
#> 4: 0 -5174.358 -2487.358 -1591.5627 -4449.590 -10373.983 -1836.370
#> 5: 0 -5170.146 -2481.878 -1588.5833 -4446.334 -10370.489 -1841.961
#> ---
#> 424484: 0 -3703.158 -2301.288 -395.4219 -2952.459 -9502.499 -2128.771
#> 424485: 0 -3834.959 -2372.439 -444.3945 -2994.426 -9465.224 -2179.546
#> 424486: 0 -3773.260 -2312.064 -419.6414 -2975.410 -9389.198 -2153.340
#> 424487: 0 -3746.392 -2304.304 -426.3353 -2986.242 -9369.264 -2141.644
#> 424488: 0 -3615.768 -2358.720 -437.9209 -2984.145 -9521.864 -2130.978
#> CP5 CP1 CP2 CP6 P7 P3 Pz
#> 1: -8938.128 -4755.838 -1369.6688 -1729.708 -7513.455 -8512.474 -537.63159
#> 2: -8963.231 -4778.917 -1391.3324 -1750.029 -7535.836 -8533.438 -557.97095
#> 3: -8969.171 -4784.986 -1397.7506 -1757.532 -7541.739 -8537.741 -561.55725
#> 4: -8957.989 -4775.938 -1391.9026 -1754.553 -7530.282 -8529.099 -552.78510
#> 5: -8950.726 -4772.205 -1388.7207 -1755.197 -7524.361 -8523.562 -549.19885
#> ---
#> 424484: -8037.571 -4200.332 -677.7635 -1046.501 -5433.473 -8150.485 -34.94090
#> 424485: -8083.840 -4253.405 -721.7708 -1080.431 -5470.565 -8197.932 -79.35287
#> 424486: -8027.603 -4225.636 -702.5348 -1048.874 -5474.501 -8173.215 -60.83399
#> 424487: -8046.398 -4230.160 -709.1735 -1040.856 -5488.533 -8182.244 -64.40166
#> 424488: -8060.742 -4244.541 -712.0606 -1053.931 -5478.491 -8187.780 -69.53255
#> P4 P8 POz O1 Oz O2 HEOG
#> 1: -99.10343 -875.6397 -852.6708 -4117.080 -1950.149 -817.4540 -5006.420
#> 2: -117.14455 -892.3933 -871.1710 -4136.997 -1965.578 -825.3984 -5027.017
#> 3: -122.82701 -892.9816 -873.9114 -4138.725 -1965.137 -824.5893 -5037.775
#> 4: -117.75130 -886.6740 -866.0771 -4126.606 -1957.450 -823.7434 -5033.031
#> 5: -117.82441 -885.3681 -865.5070 -4127.543 -1967.270 -834.0415 -5031.210
#> ---
#> 424484: 496.89807 -182.1897 -304.2810 -2730.290 -1170.744 -67.8408 -3137.224
#> 424485: 458.84891 -229.4153 -345.6953 -2777.571 -1209.235 -104.1242 -3216.466
#> 424486: 476.24600 -226.5095 -342.8815 -2794.121 -1232.792 -135.4974 -3136.874
#> 424487: 471.66693 -247.6580 -343.4331 -2790.444 -1224.516 -152.9679 -3105.096
#> 424488: 475.47369 -247.5108 -350.5685 -2804.640 -1245.077 -186.8792 -3162.988
#> VEOG
#> 1: -9116.806
#> 2: -9133.780
#> 3: -9138.781
#> 4: -9134.222
#> 5: -9139.996
#> ---
#> 424484: -8591.514
#> 424485: -8672.283
#> 424486: -8648.578
#> 424487: -8658.159
#> 424488: -8660.311
#>
#> # Events table:
#> .id .type .description .initial .final .channel
#> 1: 1 New Segment 1 1 <NA>
#> 2: 1 Stimulus s111 10528 10528 <NA>
#> 3: 1 Stimulus s130 10742 10742 <NA>
#> 4: 1 Stimulus s71 11175 11175 <NA>
#> 5: 1 Stimulus s130 12806 12806 <NA>
#> ---
#> 401: 1 Stimulus s130 411197 411197 <NA>
#> 402: 1 Stimulus s71 411469 411469 <NA>
#> 403: 1 Stimulus s121 412943 412943 <NA>
#> 404: 1 Stimulus s122 412966 412966 <NA>
#> 405: 1 Stimulus s102 412989 412989 <NA>
#>
#> # Segments table:
#> .id .recording segment
#> 1: 1 s1_faces.vhdr 1
summary(faces)
#> # EEG data:
#> # Sampling rate: 512 Hz.
#> # Size in memory: 113.4 Mb.
#> # Total duration: 00:13:49.
#> # Summary of segments
#> .recording n_segments n_incomplete
#> 1: s1_faces.vhdr 1 0
#> # Summary of events
#> .type .description n
#> 1: New Segment 1
#> 2: Stimulus s102 1
#> 3: Stimulus s111 1
#> 4: Stimulus s121 1
#> 5: Stimulus s122 1
#> 6: Stimulus s130 200
#> 7: Stimulus s70 100
#> 8: Stimulus s71 100We see that there is no electrode positions in the object, but since we know that the layout was a standard 10/20, we’ll add this layout to the object, using the dataset layout_32_1020. This is required to be able to create topographic plots.
channels_tbl(faces)
#> .channel number .reference resolution unit radius theta phi .x .y .z
#> 1: Fp1 Ch1 <NA> 1 microvolts NA NA NA NA NA NA
#> 2: Fpz Ch2 <NA> 1 microvolts NA NA NA NA NA NA
#> 3: Fp2 Ch3 <NA> 1 microvolts NA NA NA NA NA NA
#> 4: F7 Ch4 <NA> 1 microvolts NA NA NA NA NA NA
#> 5: F3 Ch5 <NA> 1 microvolts NA NA NA NA NA NA
#> 6: Fz Ch6 <NA> 1 microvolts NA NA NA NA NA NA
#> 7: F4 Ch7 <NA> 1 microvolts NA NA NA NA NA NA
#> 8: F8 Ch8 <NA> 1 microvolts NA NA NA NA NA NA
#> 9: FC5 Ch9 <NA> 1 microvolts NA NA NA NA NA NA
#> 10: FC1 Ch10 <NA> 1 microvolts NA NA NA NA NA NA
#> 11: FC2 Ch11 <NA> 1 microvolts NA NA NA NA NA NA
#> 12: FC6 Ch12 <NA> 1 microvolts NA NA NA NA NA NA
#> 13: M1 Ch13 <NA> 1 microvolts NA NA NA NA NA NA
#> 14: T7 Ch14 <NA> 1 microvolts NA NA NA NA NA NA
#> 15: C3 Ch15 <NA> 1 microvolts NA NA NA NA NA NA
#> 16: Cz Ch16 <NA> 1 microvolts NA NA NA NA NA NA
#> 17: C4 Ch17 <NA> 1 microvolts NA NA NA NA NA NA
#> 18: T8 Ch18 <NA> 1 microvolts NA NA NA NA NA NA
#> 19: M2 Ch19 <NA> 1 microvolts NA NA NA NA NA NA
#> 20: CP5 Ch20 <NA> 1 microvolts NA NA NA NA NA NA
#> 21: CP1 Ch21 <NA> 1 microvolts NA NA NA NA NA NA
#> 22: CP2 Ch22 <NA> 1 microvolts NA NA NA NA NA NA
#> 23: CP6 Ch23 <NA> 1 microvolts NA NA NA NA NA NA
#> 24: P7 Ch24 <NA> 1 microvolts NA NA NA NA NA NA
#> 25: P3 Ch25 <NA> 1 microvolts NA NA NA NA NA NA
#> 26: Pz Ch26 <NA> 1 microvolts NA NA NA NA NA NA
#> 27: P4 Ch27 <NA> 1 microvolts NA NA NA NA NA NA
#> 28: P8 Ch28 <NA> 1 microvolts NA NA NA NA NA NA
#> 29: POz Ch29 <NA> 1 microvolts NA NA NA NA NA NA
#> 30: O1 Ch30 <NA> 1 microvolts NA NA NA NA NA NA
#> 31: Oz Ch31 <NA> 1 microvolts NA NA NA NA NA NA
#> 32: O2 Ch32 <NA> 1 microvolts NA NA NA NA NA NA
#> 33: HEOG Ch33 <NA> 1 microvolts NA NA NA NA NA NA
#> 34: VEOG Ch34 <NA> 1 microvolts NA NA NA NA NA NA
#> .channel number .reference resolution unit radius theta phi .x .y .z
## In case the order of the electrodes is different, we do a left_join instead of replacing the table:
channels_tbl(faces) <- select(channels_tbl(faces), .channel) %>%
left_join(layout_32_1020)
#> Joining, by = ".channel"The plots that eeguana produce are ggplot objects that can be modified like regular ggplot.
plot(faces) +
ggtitle("All the experiment")
#> # Downsampling from 512Hz to 256Hz.
#> # Object size in memory 56.7 Mb
Pre-processing
Rereferencing
In this dataset, the signal from all electrodes is monopolar and referenced to the left mastoid (M1). We want the signal (excluding the EOG channels) to be referenced to linked (left and right) mastoids (M1 and M2).
faces <- eeg_rereference(faces, -VEOG, -HEOG, .ref = c("M1", "M2"))Filtering
We apply a band pass filter of 0.1 to 30 Hz to all the channels except the EOG channels. We don’t segment yet, because discontinuities in the signal create artifacts on the edges.
faces_filt <- eeg_filt_band_pass(faces, -HEOG, -VEOG, .freq = c(.1, 30))
#> Setting up band-pass filter from 0.1 - 30 Hz
#> Width of the transition band at the low cut-off frequency is 0.1 Hz
#> Width of the transition band at the high cut-off frequency is 7.5 HzICA
We want to apply ICA to as much data as possible, but to “representative” data: to data of the experiment and not when the participants were moving or reading the instructions. For the same reason, we want to the ignore artifacts that we are sure are not representing brain activity, because of the extreme amplitudes in the signal.
We first cut a large segment that excludes data before and after the experiment was ran using eeg_segment() (that is, only data from the marker “s111” to “s121”), and then we mark differences of 200 microvolts between peaks with eeg_artif_minmax().
faces_ls <- eeg_segment(faces_filt, .description == "s111", .end = .description == "s121") %>%
eeg_artif_minmax(-HEOG, -VEOG, .threshold = 200, .window = 200, .unit = "ms")
#> # Total of 1 segments found.
#> # Object size in memory 107.5 Mb after segmentation.
#> # Number of intervals with artifacts: 118eeg_artif_minimax() only adds the artifacts in the events table and doesn’t modify the signal. We can have a look at the type of artifacts that were detected, by plotting the signal and the events with annotate_events(). It would be a good idea to just check a couple of seconds of the signal, for that we use the dplyr-like function eeg_filter().
faces_ls %>%
eeg_select(-HEOG, -VEOG)%>%
eeg_filter(as_time(.sample, .unit = "s") %>% between(0, 90)) %>%
plot() +
annotate_events() +
theme(legend.position = "bottom")
#> # Downsampling from 512Hz to 256Hz.
#> # Object size in memory 5.9 Mb
Now we can run ICA, removing the EOG and reference electrodes, and ignoring the artifacts.
## By default, it will ignore artifacts
faces_ica <- faces_ls %>%
eeg_ica(-HEOG, -VEOG, -M1, -M2, .method = adapt_fast_ICA, .ignore = .type == "artifact")
#> # ICA is being done using adapt_fast_ICA...
#> # 97% of the samples will be used.
#> # ICA took 2.4 minsWe can now check the different topographic plots of the ICAs with plot_components(), and their correlation with the EOG channels looking a the summary of the object.
faces_ica %>% plot_components()
eeg_ica_summary_tbl(faces_ica)
#> EOG channels detected as: HEOG, VEOG
#> # Downsampling from 512Hz to 256Hz.
#> # Object size in memory 53.8 Mb
#> .recording EOG .ICA cor var
#> 1: s1_faces.vhdr VEOG ICA5 0.3110616894 3.757060e-01
#> 2: s1_faces.vhdr HEOG ICA5 0.0060600201 3.757060e-01
#> 3: s1_faces.vhdr VEOG ICA14 0.0594109202 7.014111e-02
#> 4: s1_faces.vhdr HEOG ICA14 0.0115211860 7.014111e-02
#> 5: s1_faces.vhdr VEOG ICA25 0.0191137136 5.700123e-02
#> 6: s1_faces.vhdr HEOG ICA25 0.0079409656 5.700123e-02
#> 7: s1_faces.vhdr VEOG ICA28 0.0215901453 5.238948e-02
#> 8: s1_faces.vhdr HEOG ICA28 -0.0026613292 5.238948e-02
#> 9: s1_faces.vhdr VEOG ICA13 -0.0733915974 4.883685e-02
#> 10: s1_faces.vhdr HEOG ICA13 0.0028042746 4.883685e-02
#> 11: s1_faces.vhdr VEOG ICA16 0.0189808192 3.920688e-02
#> 12: s1_faces.vhdr HEOG ICA16 -0.0052139432 3.920688e-02
#> 13: s1_faces.vhdr VEOG ICA27 -0.0146958686 3.894467e-02
#> 14: s1_faces.vhdr HEOG ICA27 -0.0005737728 3.894467e-02
#> 15: s1_faces.vhdr HEOG ICA26 -0.0127623329 3.356956e-02
#> 16: s1_faces.vhdr VEOG ICA26 -0.0103929696 3.356956e-02
#> 17: s1_faces.vhdr VEOG ICA21 -0.0110379588 3.209204e-02
#> 18: s1_faces.vhdr HEOG ICA21 0.0041385932 3.209204e-02
#> 19: s1_faces.vhdr VEOG ICA30 -0.0318988532 3.195676e-02
#> 20: s1_faces.vhdr HEOG ICA30 -0.0049492927 3.195676e-02
#> 21: s1_faces.vhdr HEOG ICA11 -0.0080689883 1.531702e-02
#> 22: s1_faces.vhdr VEOG ICA11 -0.0022034332 1.531702e-02
#> 23: s1_faces.vhdr VEOG ICA8 -0.0061078333 1.503477e-02
#> 24: s1_faces.vhdr HEOG ICA8 -0.0033160702 1.503477e-02
#> 25: s1_faces.vhdr VEOG ICA22 -0.0729992398 1.342687e-02
#> 26: s1_faces.vhdr HEOG ICA22 -0.0075857794 1.342687e-02
#> 27: s1_faces.vhdr VEOG ICA19 -0.0113037545 1.276729e-02
#> 28: s1_faces.vhdr HEOG ICA19 0.0040990970 1.276729e-02
#> 29: s1_faces.vhdr HEOG ICA29 -0.0098960226 1.268464e-02
#> 30: s1_faces.vhdr VEOG ICA29 0.0082034132 1.268464e-02
#> 31: s1_faces.vhdr VEOG ICA10 0.0029772146 1.075145e-02
#> 32: s1_faces.vhdr HEOG ICA10 -0.0012705350 1.075145e-02
#> 33: s1_faces.vhdr VEOG ICA23 -0.0139531104 1.074652e-02
#> 34: s1_faces.vhdr HEOG ICA23 0.0080386572 1.074652e-02
#> 35: s1_faces.vhdr VEOG ICA1 0.0590125612 9.976900e-03
#> 36: s1_faces.vhdr HEOG ICA1 -0.0556782278 9.976900e-03
#> 37: s1_faces.vhdr VEOG ICA2 -0.0126428949 7.926431e-03
#> 38: s1_faces.vhdr HEOG ICA2 -0.0078820112 7.926431e-03
#> 39: s1_faces.vhdr VEOG ICA24 -0.0121801971 4.817643e-03
#> 40: s1_faces.vhdr HEOG ICA24 -0.0016042248 4.817643e-03
#> 41: s1_faces.vhdr HEOG ICA7 -0.0020891334 3.213751e-03
#> 42: s1_faces.vhdr VEOG ICA7 -0.0003910153 3.213751e-03
#> 43: s1_faces.vhdr VEOG ICA12 0.0144307012 1.612854e-03
#> 44: s1_faces.vhdr HEOG ICA12 0.0128083915 1.612854e-03
#> 45: s1_faces.vhdr HEOG ICA17 -0.0091874914 1.241979e-03
#> 46: s1_faces.vhdr VEOG ICA17 0.0017426840 1.241979e-03
#> 47: s1_faces.vhdr VEOG ICA18 -0.0131764353 8.322458e-04
#> 48: s1_faces.vhdr HEOG ICA18 -0.0117696565 8.322458e-04
#> 49: s1_faces.vhdr HEOG ICA4 -0.0190202998 5.031384e-04
#> 50: s1_faces.vhdr VEOG ICA4 -0.0171506626 5.031384e-04
#> 51: s1_faces.vhdr VEOG ICA6 -0.0089105463 4.393376e-05
#> 52: s1_faces.vhdr HEOG ICA6 0.0049753303 4.393376e-05
#> 53: s1_faces.vhdr VEOG ICA20 -0.0062150226 1.952630e-06
#> 54: s1_faces.vhdr HEOG ICA20 -0.0011994425 1.952630e-06
#> 55: s1_faces.vhdr HEOG ICA9 -0.0019097967 -4.193651e-04
#> 56: s1_faces.vhdr VEOG ICA9 0.0015343447 -4.193651e-04
#> 57: s1_faces.vhdr VEOG ICA15 -0.0265727074 -1.012212e-03
#> 58: s1_faces.vhdr HEOG ICA15 -0.0041488999 -1.012212e-03
#> 59: s1_faces.vhdr VEOG ICA3 0.0866623859 -7.459639e-03
#> 60: s1_faces.vhdr HEOG ICA3 -0.0019740756 -7.459639e-03
#> .recording EOG .ICA cor varWe’ll look closer at the ICA that are more likely to be related to eye movements, and we’ll compare with signals that look like blinks and saccades:
faces_ica <- faces_ica %>%
eeg_artif_step(HEOG, VEOG, .threshold = 30, .window = 200, .unit = "ms", .freq = c(1, 10)) %>%
eeg_artif_peak(VEOG, .threshold = 100, .freq = c(1, 10))
#> Setting up band-pass filter from 1 - 10 Hz
#> Width of the transition band at the low cut-off frequency is 1 Hz
#> Width of the transition band at the high cut-off frequency is 2.5 Hz
#> # Number of intervals with artifacts: 407
#> Setting up band-pass filter from 1 - 10 Hz
#> Width of the transition band at the low cut-off frequency is 1 Hz
#> Width of the transition band at the high cut-off frequency is 2.5 Hz
#> # Number of intervals with artifacts: 158In order to investigate the ICA, we are using an experimental function not yet exported (that’s why we use the :::). This function will probably converted into a shiny app. It’s possible to get the same plots by using a combination of plot_topo, eeg_ica_show, eeg_ica_summary_tbl, and some ggplot and cowplot functions.
We’ll see how the components explaining the most variance and most correlated to the EOG channels behave:
#samples with blinks
s_peaks <- filter(events_tbl(faces_ica), str_starts(.description, "peak")) %>%
pull(.initial)
eeguana:::plot_ica.eeg_ica_lst(faces_ica, samples = seq(s_peaks[1]-2000,s_peaks[1]+2000))
#> Warning in eeguana:::plot_ica.eeg_ica_lst(faces_ica, samples = seq(s_peaks[1]
#> - : This is an experimental function, and it might change or disappear in the
#> future. (Or it might be transformed into a shinyapp)
#> Using recording: s1_faces.vhdr
#> EOG channels detected as: HEOG, VEOG
#> Calculating the correlation of ICA components with filtered EOG channels...
#> Setting up band-pass filter from 0.1 - 30 Hz
#> Width of the transition band at the low cut-off frequency is 0.1 Hz
#> Width of the transition band at the high cut-off frequency is 7.5 Hz
#> EOG channels detected as: HEOG, VEOG
#> # Downsampling from 512Hz to 256Hz.
#> # Object size in memory 53.9 Mb
#> The following columns of signal_tbl are not channels (or ICA components): ICA5ICA14ICA25ICA28ICA13ICA16ICA27ICA26ICA21ICA30ICA11ICA8ICA22ICA19ICA29ICA10
#> * To build a channel use `channel_dbl()` function, e.g. channel_dbl(0) to populate the table with a channel containing 0 microvolts.
#> * To copy the structure of an existing channel one can do `new_ch = existing_channel * 0 + ...`
The clearest ones is ICA5, and so we’ll just remove that one.
faces_icaed <- faces_ica %>%
eeg_ica_keep(-ICA5) %>%
as_eeg_lst()Now we’ll segment the data to appropriate check if there are still artifacts.
events_tbl(faces_icaed) <- events_tbl(faces_icaed) %>%
filter(!.type %in% "artifact")
faces_seg <- faces_icaed %>%
eeg_select(-description, -type) %>%
eeg_segment(.description %in% c("s70", "s71"), .lim = c(-.1, .5))
#> # Total of 200 segments found.
#> # Object size in memory 16.5 Mb after segmentation.
faces_seg_artif <- faces_seg %>%
eeg_artif_minmax(-HEOG, -VEOG, .threshold = 100, .window = 150, .unit = "ms") %>%
eeg_artif_step(-HEOG, -VEOG, .threshold = 50, .window = 200, .unit = "ms")
#> # Number of intervals with artifacts: 0
#> # Number of intervals with artifacts: 1
## extracts the ids of the segments with artifacts
bad <- filter(events_tbl(faces_seg_artif), .type == "artifact") %>%
pull(.id) %>%
unique()
## Show the segment with artifact and one before and after:
faces_seg_artif %>%
eeg_filter(.id %in% c(bad-1, bad, bad+1)) %>%
eeg_select(-VEOG, -HEOG) %>%
plot() +
annotate_events() +
theme(legend.position = "bottom")
faces_seg <- faces_seg_artif %>%
eeg_events_to_NA(.type == "artifact", .entire_seg = TRUE, .drop_events = TRUE)
summary(faces_seg)
#> # EEG data:
#> # Sampling rate: 512 Hz.
#> # Size in memory: 16.5 Mb.
#> # Total duration: 00:02:00.
#> # Summary of segments
#> .recording n_segments n_incomplete
#> 1: s1_faces.vhdr 200 1
#> # Summary of events
#> .type .description n
#> 1: Stimulus s70 100
#> 2: Stimulus s71 100Finally, we can baseline the segments:
faces_seg <- faces_seg %>% eeg_baseline()Visualization
We edit the segmentation information and add more descriptive labels.
faces_seg <- faces_seg %>%
eeg_mutate(
condition =
ifelse(description == "s70", "faces", "non-faces")
) %>%
eeg_select(-type)
faces_seg
#> # EEG data:
#>
#> # Signal table:
#> .id .sample Fp1 Fpz Fp2 F7 F3
#> 1: 1 -50 -0.2239737 2.926696 0.1088197 -0.1107286 3.793132
#> 2: 1 -49 0.9975430 3.910457 0.7585352 0.7172289 4.038747
#> 3: 1 -48 2.2477072 4.853125 1.5129425 1.7686192 4.340948
#> 4: 1 -47 3.3361201 5.596977 2.2208750 2.9135916 4.642705
#> 5: 1 -46 4.0950011 6.007830 2.7411833 4.0297893 4.899710
#> ---
#> 61596: 200 253 -7.4272958 -4.030794 -3.3414518 2.5704919 4.188729
#> 61597: 200 254 -7.0028675 -4.076854 -3.7139526 1.3456003 3.356811
#> 61598: 200 255 -6.3948136 -4.124318 -4.0795149 0.2178476 2.644884
#> 61599: 200 256 -5.7744059 -4.236721 -4.4670174 -0.8343465 1.996052
#> 61600: 200 257 -5.3225204 -4.474952 -4.9176401 -1.8464510 1.323714
#> Fz F4 F8 FC5 FC1 FC2
#> 1: 4.228278 2.840110 -3.15472635 1.6113910 1.592771 -2.364200
#> 2: 4.270218 2.382479 -2.35472265 2.0405680 1.860704 -2.861841
#> 3: 4.421612 2.013301 -1.24644410 2.4513040 2.219617 -3.123105
#> 4: 4.591250 1.709526 0.05225002 2.8087541 2.615841 -3.167049
#> 5: 4.700407 1.452931 1.40978255 3.1051410 3.005053 -3.015121
#> ---
#> 61596: -8.730353 -12.167556 -1.36556543 4.6994268 -3.708406 -19.005953
#> 61597: -9.453131 -11.842863 -2.17206941 3.4988400 -4.958880 -19.459708
#> 61598: -9.985577 -11.460650 -2.78238153 2.2688876 -6.309382 -19.836152
#> 61599: -10.366157 -11.047480 -3.20156158 1.0080235 -7.730844 -20.137723
#> 61600: -10.661645 -10.644149 -3.46204537 -0.3005228 -9.204715 -20.380105
#> FC6 M1 T7 C3 Cz C4
#> 1: -1.8391090 -2.030141 8.787014 -2.0086889 -5.711339 -0.795271
#> 2: -1.1455898 -1.922150 9.114149 -2.1197829 -5.858017 -1.282006
#> 3: -0.1976263 -1.785010 8.977623 -2.2059909 -5.721134 -1.762303
#> 4: 0.8983126 -1.654058 8.398226 -2.2445165 -5.330143 -2.204568
#> 5: 2.0401272 -1.545839 7.452935 -2.1993214 -4.723109 -2.557211
#> ---
#> 61596: -12.4167239 4.310870 4.908478 -0.5437559 -14.854063 -34.752050
#> 61597: -13.3369315 4.454409 4.680512 -2.1197706 -16.258766 -35.506742
#> 61598: -14.2119912 4.482182 4.315460 -3.8504649 -17.808367 -36.200622
#> 61599: -14.9925296 4.394789 3.730447 -5.6655686 -19.413279 -36.808126
#> 61600: -15.6416067 4.203761 2.864558 -7.5086186 -20.985826 -37.301315
#> T8 M2 CP5 CP1 CP2 CP6
#> 1: 3.910053 2.030141 2.1730493 -5.558311 -3.912675 2.0442154
#> 2: 5.477720 1.922150 1.6874719 -5.681516 -4.454747 1.4513845
#> 3: 7.069972 1.785010 1.0244766 -5.664134 -4.852158 0.8377522
#> 4: 8.401729 1.654058 0.2760940 -5.506764 -5.089553 0.2248167
#> 5: 9.193991 1.545839 -0.4454964 -5.204630 -5.150715 -0.3509800
#> ---
#> 61596: -25.822482 -4.310870 3.4352795 -5.102525 -31.824672 -29.5344073
#> 61597: -27.019577 -4.454409 2.4188361 -6.646814 -33.109751 -30.5440118
#> 61598: -27.809796 -4.482182 1.1739532 -8.403855 -34.341884 -31.5585170
#> 61599: -28.198220 -4.394789 -0.2539754 -10.263130 -35.433438 -32.5685325
#> 61600: -28.211914 -4.203761 -1.8208137 -12.116925 -36.299732 -33.5345508
#> P7 P3 Pz P4 P8 POz
#> 1: -1.8486407 -2.6437673 -3.191924 -0.6720984 0.5889623 -1.233383
#> 2: -2.6139953 -3.1920017 -3.790144 -1.4646696 -1.8056977 -2.397317
#> 3: -3.1638517 -3.6552420 -4.342437 -2.1770814 -4.0661883 -3.500260
#> 4: -3.4395628 -3.9730269 -4.813028 -2.7552677 -5.9136674 -4.440224
#> 5: -3.4148548 -4.0892851 -5.160562 -3.1471226 -7.1129100 -5.124399
#> ---
#> 61596: 2.7192767 2.7365646 -12.488676 -27.5485924 -23.0238529 -5.742807
#> 61597: 2.3539444 1.6113645 -13.716290 -28.1890858 -22.9764754 -6.553364
#> 61598: 1.6063736 0.2901116 -14.977399 -28.6721202 -22.1172213 -7.394155
#> 61599: 0.5699239 -1.1512356 -16.183334 -28.9837949 -20.4604781 -8.207270
#> 61600: -0.6454750 -2.6374742 -17.250028 -29.1134647 -18.1089704 -8.943773
#> O1 Oz O2 HEOG VEOG
#> 1: -3.5386535 -0.7969323 0.8721568 -8.172153 -4.070868
#> 2: -4.7377453 -2.2782486 -0.7106985 -15.399203 -10.286688
#> 3: -5.8496010 -3.7413040 -2.2901256 -17.091586 -4.218329
#> 4: -6.7658792 -5.0345596 -3.6833116 -22.001254 -2.728094
#> 5: -7.3928468 -6.0157379 -4.7178515 -14.368930 3.892999
#> ---
#> 61596: 2.5152642 -3.0701651 -5.1484950 40.897595 104.275391
#> 61597: 1.8963234 -3.5788828 -3.8985745 19.307751 67.568359
#> 61598: 1.0904936 -4.1790541 -2.9052267 24.585583 66.189453
#> 61599: 0.2077071 -4.8230480 -2.2435981 4.779675 43.698242
#> 61600: -0.6593847 -5.4755758 -1.9725447 -2.815540 36.453125
#>
#> # Events table:
#> .id .type .description .initial .final .channel
#> 1: 1 Stimulus s71 1 1 <NA>
#> 2: 2 Stimulus s71 1 1 <NA>
#> 3: 3 Stimulus s71 1 1 <NA>
#> 4: 4 Stimulus s71 1 1 <NA>
#> 5: 5 Stimulus s71 1 1 <NA>
#> ---
#> 196: 196 Stimulus s71 1 1 <NA>
#> 197: 197 Stimulus s70 1 1 <NA>
#> 198: 198 Stimulus s71 1 1 <NA>
#> 199: 199 Stimulus s71 1 1 <NA>
#> 200: 200 Stimulus s71 1 1 <NA>
#>
#> # Segments table:
#> .id .recording segment description condition
#> 1: 1 s1_faces.vhdr 1 s71 non-faces
#> 2: 2 s1_faces.vhdr 2 s71 non-faces
#> 3: 3 s1_faces.vhdr 3 s71 non-faces
#> 4: 4 s1_faces.vhdr 4 s71 non-faces
#> 5: 5 s1_faces.vhdr 5 s71 non-faces
#> ---
#> 196: 196 s1_faces.vhdr 196 s71 non-faces
#> 197: 197 s1_faces.vhdr 197 s70 faces
#> 198: 198 s1_faces.vhdr 198 s71 non-faces
#> 199: 199 s1_faces.vhdr 199 s71 non-faces
#> 200: 200 s1_faces.vhdr 200 s71 non-facesWith some ggplot skills, we can create customized plots. ggplot is overloaded to work on an eeg_lst by first downsampling the signal when necessary, and converting it to a long-format data frame that is feed into ggplot. This object can then be customized. (Notice that the channels, or component names are in a .key column and their amplitude in .value column, and instead of samples there are now .time in seconds).
## ggplot uses internally a table that looks like this:
faces_seg %>% eeg_filter(.id %in% 1:3) %>%
as_tibble()
#> # A tibble: 31,416 × 8
#> .time .id .key .value .recording segment description condition
#> <dbl> <int> <chr> <dbl> <chr> <int> <chr> <chr>
#> 1 -0.0996 1 Fp1 -0.224 s1_faces.vhdr 1 s71 non-faces
#> 2 -0.0977 1 Fp1 0.998 s1_faces.vhdr 1 s71 non-faces
#> 3 -0.0957 1 Fp1 2.25 s1_faces.vhdr 1 s71 non-faces
#> 4 -0.0938 1 Fp1 3.34 s1_faces.vhdr 1 s71 non-faces
#> 5 -0.0918 1 Fp1 4.10 s1_faces.vhdr 1 s71 non-faces
#> 6 -0.0898 1 Fp1 4.40 s1_faces.vhdr 1 s71 non-faces
#> 7 -0.0879 1 Fp1 4.19 s1_faces.vhdr 1 s71 non-faces
#> 8 -0.0859 1 Fp1 3.49 s1_faces.vhdr 1 s71 non-faces
#> 9 -0.0840 1 Fp1 2.36 s1_faces.vhdr 1 s71 non-faces
#> 10 -0.0820 1 Fp1 0.972 s1_faces.vhdr 1 s71 non-faces
#> # … with 31,406 more rows
faces_seg %>%
eeg_select(O1, O2, P7, P8) %>%
ggplot(aes(x = .time, y = .value)) +
geom_line(alpha = .1, aes(group = .id, color = condition)) +
stat_summary(
fun = "mean", geom = "line", alpha = 1, size = 1.5,
aes(color = condition)
) +
facet_wrap(~.key) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_vline(xintercept = .17, linetype = "dotted") +
theme(legend.position = "bottom")
We can see here the N170 component in the faces condition for the average trial on a single participant. We can investigate the signal by averaging the channels of the occipital and parietal lobes using transmute(), and the special function, chs_mean(), a wrapper for rowMeans(), which takes as arguments the relevant channels and whether missing values should be omitted from the calculations.
faces_seg %>%
eeg_transmute(
Occipital = chs_mean(across(starts_with("O")), na.rm = TRUE), #O1, O2, Oz
Parietal = chs_mean(across(starts_with("O")), na.rm = TRUE) # P3, P4, P7, P8, Pz
) %>%
ggplot(aes(x = .time, y = .value)) +
geom_line(alpha = .1, aes(group = .id, color = condition)) +
stat_summary(
fun = "mean", geom = "line", alpha = 1, size = 1,
aes(color = condition)
) +
facet_wrap(~.key) +
theme(legend.position = "bottom")
We can also calculate the ERPs and then plot them in their layout (and we’ll add the same theme that plot uses for eeg_lsts):
ERP_faces <- faces_seg %>%
eeg_group_by(.sample, condition) %>%
eeg_summarize(across_ch(mean, na.rm = TRUE))
ERP_plot <- ERP_faces %>%
ggplot(aes(x = .time, y = .value)) +
geom_line(aes(color = condition)) +
facet_wrap(~.key) +
theme(legend.position = "bottom") +
ggtitle("ERPs for faces vs non-faces") +
theme_eeguana()
ERP_plot %>% plot_in_layout()
Another possibility is to create a topographic plot of the two conditions, by first making segments that include only the interval .1-.2 s after the onset of the stimuli, creating a table with interpolated amplitudes and using the ggplot wrapper plot_topo().
faces_seg %>%
eeg_filter(between(as_time(.sample, .unit = "s"), .1, .2)) %>%
eeg_group_by(condition) %>%
eeg_summarize(across_ch(mean, na.rm = TRUE)) %>%
plot_topo() +
annotate_head() +
geom_contour() +
geom_text(colour = "black") +
facet_grid(~condition)
For more specialized plots or analyses, it might be necessary to extract the data in long data frame format first. Even though, we can do this without transforming the object, here, we transform the data first and then we visualize independent t-tests at every electrode and time point.
df <- faces_seg %>%
eeg_select(O1, O2, P7, P8) %>%
as_tibble() %>%
# We can use regular dplyr functions now
group_by(.key, .time) %>%
summarize(
`t-value` = t.test(
.value[condition == "faces"],
.value[condition == "non-faces"]
)$statistic
)
#> `summarise()` has grouped output by '.key'. You can override using the
#> `.groups` argument.
df
#> # A tibble: 1,232 × 3
#> # Groups: .key [4]
#> .key .time `t-value`
#> <chr> <dbl> <dbl>
#> 1 O1 -0.0996 0.204
#> 2 O1 -0.0977 0.227
#> 3 O1 -0.0957 0.234
#> 4 O1 -0.0938 0.219
#> 5 O1 -0.0918 0.173
#> 6 O1 -0.0898 0.0880
#> 7 O1 -0.0879 -0.0415
#> 8 O1 -0.0859 -0.221
#> 9 O1 -0.0840 -0.448
#> 10 O1 -0.0820 -0.713
#> # … with 1,222 more rowsThen we just load the data frame into ggplot.
ggplot(df, aes(x = .time, y = `t-value`)) + geom_line() +
facet_wrap(~.key)
However, this can be also done in the without transforming the eeg_lst object:
faces_seg_t <-
faces_seg %>%
eeg_select(O1, O2, P7, P8) %>%
eeg_group_by(.sample) %>%
eeg_summarize(across_ch(list(t = ~t.test(
.[condition == "faces"],
.[condition == "non-faces"]
)$statistic)))
faces_seg_t %>%
ggplot(aes(x = .time, y = .value)) +
geom_line(alpha = .1, aes(group = .id)) +
stat_summary(fun = "mean", geom = "line", alpha = 1, size = 1) +
facet_wrap(~.key) +
theme(legend.position = "bottom")