Spectral Analysis of Sleep EEG Signals
In rsleep: Analysis of Sleep Data

options(scipen=999)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>")

This vignette provides a basic introduction to spectral analysis of Electroencephalography (EEG) signal from a sleep record file.

EEG refers to all the methods of recording, analysis and interpretation of the electrical activity of the brain. In clinical EEG, multiple electrodes are usually placed on the scalp, measuring its superficial activity over time. Electrodes are typically arranged using the standardized International 10-20 system @niedermeyerElectroencephalographyBasicPrinciples2005, in order to enhance analysis reproducibility.

EEG is a major component in sleep analysis. Sleep stages, such as slow wave sleep or paradoxical sleep are partly defined over visual EEG characteristics @AASMScoringManual. Many sleep related disorders can be identified in EEG data. Polysomnography (PSG), the gold standard exam in sleep medicine, includes EEG along many other physiological signals @ibanezSurveySleepAssessment2018.

Sleep data

Sleep records are usually stored using European Data Format (EDF) @Kemp1992. The R package edfReader reads .edf files. Reading an .edf file takes two steps: First reading the headers of the file, then reading the selected signals. The following spectral analysis will be performed on a single channel of the EEG, the C3-M2 central derivation.

library(edfReader)

fname <- "15012016HD.edf"
if(!file.exists(fname)){
  download.file(
  url = "https://rsleep.org/data/15012016HD.edf",
  destfile = fname,
  method = "curl")}

h <- readEdfHeader(fname)

s <- readEdfSignals(h, signals = "C3-M2")

library(edfReader)

download.file(
  url = "https://rsleep.org/data/15012016HD.edf",
  destfile = "15012016HD.edf",
  method = 'curl')

h <- readEdfHeader("15012016HD.edf")

s <- readEdfSignals(h, signals = "C3-M2")

The rsleep function spectrogram() plots the spectrogram of the signal. But first, install rsleep latest version from Github.

devtools::install_github("boupetch/rsleep")

library(rsleep)

spectrogram(
  signal = s$signal,
  sRate = s$sRate,
  startTime = s$startTime)

library(rsleep)

spectrogram(
  signal = s$signal,
  sRate = s$sRate,
  startTime = s$startTime)

rsleep reads this particular format with the read_events_noxturnal function. The plot_hypnogram function then plots the hypnogram.

if(!file.exists("15012016HD.csv")){
  download.file(
  url = "https://rsleep.org/data/15012016HD.csv",
  destfile = "15012016HD.csv")}

events <- rsleep::read_events_noxturnal(
  "15012016HD.csv", h$startTime)

rsleep::plot_hypnogram(events)

if(!file.exists("15012016HD.csv")){
  download.file(
  url = "https://rsleep.org/data/15012016HD.csv",
  destfile = "15012016HD.csv",
  method="curl")}

events <- read_events_noxturnal("15012016HD.csv")

# Remove last epoch as signal stops before.
events <- head(events, -1)

# Remove other events
events <- events[events$event %in% c("AWA", "REM", "N1", "N2", "N3"),]

plot_hypnogram(events)

The hypnogram show sleep stages over time using consecutive 30 seconds epochs. This record was manually scored by well-trained technicians according to the American Academy of Sleep Medicine manual @AASMScoringManual. Five sleep stages can be observed:

AWA: Awakeness.
REM: Rapid-Eye-Movement (REM) sleep, or paradoxical sleep, highlighted in red.
N1: Light sleep transitional stage.
N2: Light sleep.
N3: Slow-wave sleep.

Visual scoring is an empirical science requiring a considerable amount of knowledge and training. Alternative methods like spectral estimations techniques such as the Fourier transform must be used to quantify information carried in the physiological signals @tongQuantitativeEEGAnalysis2009.

Epoching

Epoching is the first step of sleep analysis. Physiological signal, such as EEG, is splitted into consecutive epochs of a discrete duration. Epochs usually start at lights turnoff, when the patient or subject starts the night.

As the example record already has a hynogram, the EEG signal can be splitted using these scored epochs. The epochs function from the rsleep package split the signal according to these parameters. It returns a list of signal vectors.

epochs <- epochs(
  signals = s$signal,
  sRates = s$sRate,
  epoch = events,
  startTime = as.numeric(as.POSIXct(h$startTime)))

Periodogram

The Fourier transform (FT) may be the most important function in signal analysis. It decomposes the signal into its constituent frequencies and computes its power spectral densities (PSD). However, EEG signals also carry a lot of noise. This noise is easily intereprted by the FT and can jam the results. To solve this problem, Welch's method split the signal into overlapping segments to average power spectral densities from the Fourier transform.

Single epoch

The pwelch function rsleep computes a periodogram using Welch's method. The following example computes and plot the periodogram of the 120th epoch. This epoch has been scored N2, or light sleep, by the expert.

p <- pwelch(epochs[[120]], sRate = s$sRate)

summary(p)

This epoch periodogram shows high PSD in lower frequencies of the spectrum. As values are normalized using log, PSD are negative.

Stages profiles

To compute average periodograms by stage, hypnogram and epochs can be iterated simultaneously using the mapply function. Periodograms can be filtered at this step to discard values over 30 Hertz.

periodograms <- mapply(
  x = epochs, 
  y = events$event,
  FUN = function(x,y){
    p <- pwelch(x, sRate = s$sRate, show = FALSE)
    p <- as.data.frame(p[p$hz <= 30,])
    p$stage <- y
    p
}, SIMPLIFY = F)

mapply returns a list that can be coerced to a dataframe using rbind combined to do.call.

periodograms_df <- do.call("rbind", periodograms)

Once coerced to a dataframe, raw periodogram values can be averaged by stage.

avg_periodograms <- aggregate(psd ~ hz+stage, periodograms_df, mean)

Aggregated periodograms can then be plotted using ggplot2.

library(ggplot2)

palette <- c("#5BBCD6","#00A08A","#F2AD00","#F98400","#FF0000")

ggplot(avg_periodograms, aes(x=hz,y=psd,color=stage)) +
  geom_line() + theme_bw() +
  theme(legend.title = element_blank()) + 
  scale_colour_manual(name = "stage",
                      values = palette) +
  xlab("Frequency (Hertz)") + ylab("PSD")

Each sleep stage show a distinct average periodogram. If the N3 stage averages higher PSD values in the lower spectrum, it show way lower PSD in the upper frequencies compared to other stages.

Bands

The traditional way to simplify the EEG periodogram is to cut the frequencies of the spectrum into bands or ranges @niedermeyerElectroencephalographyBasicPrinciples2005:

Delta: Below 3.5 Hertz, the Delta band is associated with slow-wave sleep in adults subjects.
Theta: Between 3.5 and 7.5 Hertz, the Theta band is associated with drowsiness in adults and teens subjects.
Alpha: Between 7.5 and 13 Hertz, the Alpha band is associated with a relaxed state and eyes closed.
Beta: Between 13 and 30 Hertz the Beta band is associated with active thinking, focus, high alert or anxiousness.
Gamma: Over 30 Hertz the Gamma band is mostly used for animal EEG analysis.

Bands can be computed using the bands_psd of the rsleep package. Those bands can be normalized by the spectrum range covered by the bands.

bands <- lapply(epochs,function(x){
    bands_psd(
      bands = list(c(0.5,3.5), # Delta
                   c(3.5,7.5), # Theta
                   c(7.5,13), # Alpha
                   c(13,30)), # Beta
      signal = x, sRate = s$sRate)
})

As lapply returns a list, results must be reshaped in order to obtain a dataframe object.

bands_df <- data.frame(matrix(unlist(bands), nrow=length(bands), byrow=TRUE))

colnames(bands_df) <- c("Delta","Theta","Alpha","Beta")

Stages can be retreived from the hypnogram.

bands_df$stage <- rsleep::hypnogram(events)$event

Now that the epochs bands PSD and their corresponding stages are stored in a dataframe, they can easily be plotted using boxplots from ggplot2.

bands_df_long <- reshape2::melt(bands_df, "stage")

palette <-c("#F98400", "#F2AD00", "#00A08A", "#FF0000", "#5BBCD6")

ggplot(bands_df_long,
       aes(x=stage,y=value,color=stage)) +
  geom_boxplot() +
  facet_grid(rows = vars(variable),scales = "free") +
  scale_colour_manual(name = "stage",
                      values = palette) +
  theme_bw() + xlab("") + ylab("PSD") + 
  theme(legend.position = "none")