Analysis.Rmd

---
title: "Antelope Project Analysis"
author: "Andrea Stocco"
date: "April 29, 2020"
output:
  html_document:
    code_folding: hide
    theme: yeti
    toc: yes
    toc_depth: 3
    toc_float: yes
  pdf_document:
    toc: yes
    toc_depth: '3'
  word_document:
    toc: yes
    toc_depth: '3'
---

```{r setup, include=FALSE, warning=FALSE}
library(tidyverse)
library(kableExtra)
library(xtable)
library(data.table)
library(ggplot2)
library(ggthemes)
library(ggExtra)
library(colorspace)
library(RColorBrewer)
library(gridExtra)
library(ggdendro)
library(viridis)

source("./code/QEEG_Emotiv_Analysis_Rscript.R")
knitr::opts_chunk$set(echo = TRUE)
```

# Behavioral Data

## Demographics

First, let's load the data and transform it into a Wide format table

```{r}
behav <- read_tsv("data/behav/behav_data.txt",
                  col_types=cols())

behav$material[behav$material=="Swahili"] <- "Vocabulary"
behav$material[behav$material=="US Maps"] <- "Maps"

behav <- behav %>% rename(Gender = gender, Age = age)

wbehav <- behav %>% pivot_wider(values_from = alpha,
                                id_cols = c(Subject, Gender, Age),
                                names_from = material)


```

At this point, we can look at the participant demographocs

```{r}
ggplot(data=wbehav, aes(x=Age, col=Gender)) +
  geom_histogram(aes(col="white", fill=Gender),
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 1) +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  ggtitle("Age Distribution by Gender") + 
  theme_pander()
```

## Rate of Forgetting (Alpha parameter) for the two tasks

We can now examine the distribution and correlation of the rate of forgetting 
for the verbal (Vocabulary) and Maps: 

```{r}
mu <- behav %>% 
  group_by(material) %>% 
  summarize(alpha=mean(alpha))

p1 <- ggplot(behav, aes(x=alpha, fill=material)) +
  geom_histogram(col="white", binwidth = 0.025, 
                 alpha=0.5, position="identity") +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  geom_vline(data=mu, aes(xintercept=alpha, color=material),
             linetype="dashed") +
  xlab(expression(paste("Estimated value of ", alpha))) +
  theme_pander() +
  ylab("Count") +
  ggtitle("(A) Distribution By Materials") +
  theme(legend.position = c(0.2, 0.8)) +
  theme(plot.title = element_text(hjust = 0.5))

p2 <- ggplot(wbehav, aes(x=Vocabulary, y=Maps)) +
  geom_point(size=4, alpha=0.5, col="black") + #col=K[7]) +
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() + 
  scale_y_continuous() + 
  ggtitle("(B) Correlation Across Materials") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab(expression(paste(alpha, " Maps"))) +
  geom_rug(col="black", lwd=1, alpha=.5) +
  annotate("segment", x=0.1, y=0.1, xend=0.5, 
           yend=0.5, col="grey", lwd=1, lty=2) +
  theme(plot.title = element_text(hjust = 0.5))
  
grid.arrange(p1, p2, ncol=2)
```

```{r fig.width=5, fig.height=5 }
p3 <- ggplot(wbehav, aes(x=Vocabulary, y=Maps)) +
  geom_point(size=4, alpha=0.5, col="black") + #col=K[7]) +
  geom_smooth(method = "lm", formula = y ~ x, col="red", 
              fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() + 
  scale_y_continuous() + 
  coord_fixed() +
  ggtitle("Correlation of Alpha Parameters\nAcross Materials") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab(expression(paste(alpha, " Maps"))) +
  annotate("segment", x=0.1, y=0.1, xend=0.5, 
           yend=0.5, col="grey", lwd=1, lty=2) +
  annotate("text", label = paste("r =", round(cor(wbehav$Vocabulary,
                                                  wbehav$Maps),
                                              2)), 
                                 x=0.2, y=0.45, col="red") +
  
  theme(plot.title = element_text(hjust = 0.5))

ggMarginal(p3, col="white", fill="black", alpha=0.5,
           type="hist", bins=13)
```

And here is a summary of the data:

```{r}
summary(wbehav) %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

# Eyes-Closed EEG Data

Now we look at the neural correlates of the rates of forgetting in EEG data. We
begin with _eyes closed_ EEG data, which is the most common type of 
resting-state EEG recording.

## (Optional) Preprocess the Data

Set the `process_raw` variable to `TRUE` to reprocess the raw EEG data.
That might take a __very__ long time. 

```{r}
process_raw = F
if (process_raw) {
  setwd("./data/eyes_closed/")
  for (sub in dir()[grep("A[1-9]", dir())]) {
    setwd(sub)
    analyze.logfile(sub, "closed")
    setwd("..")
  }
  setwd("../..")
}
```

## Load the Eyes-Closed Spectra and Summary Files

The preprocessing step produces a number of text files. Here, we are interested
in the _summary_ file, which contains many useful statistics, including the 
Individual Alpha Frequency (i.e., the Alpha Peak) for each channel, the 
_spectra_ file, which contains the estimated log power for each channel in
increments of 0.5 Hz.

```{r}
ec_spectra <- NULL
for (sub in dir("data/eyes_closed/")[grep("A[1-9]",
                                        dir("data/eyes_closed/"))]) {
  table <- read_tsv(paste("data/eyes_closed/", sub, "/", 
                          sub, "_closed_spectra.txt", sep=""),
                    col_types = cols())
  if (is.null(ec_spectra)) {
    ec_spectra <- table
  } else {
    ec_spectra <- ec_spectra %>% bind_rows(table)
  }
}

ec_summary <- NULL
for (sub in dir("data/eyes_closed/")[grep("A[1-9]", 
                                          dir("data/eyes_closed/"))]) {
  table <- read_tsv(paste("data/eyes_closed/", sub, "/", 
                          sub, "_closed_summary.txt", sep=""),
                    col_types = cols())
  if (is.null(ec_summary)) {
    ec_summary <- table
  } else {
    ec_summary <- ec_summary %>% bind_rows(table)
  }
}
```

### Individual Alpha frquency

The IAF is perhaps the most important subject level characteristic of the
resting state spectrum. For this reason, we want to get some basic statistics.
First, we calculate the subject-level IAF by computing the mode IAF across all
channels for a given participant.

```{r}
iaf_cols <- paste(c("AF3", "AF4", "F3", "F4", 
                    "F7", "F8", "FC5", "FC6", 
                    "T7", "T8", "P7", "P8", 
                    "O1", "O2"), 
                  "IAF", 
                  sep="_")

ec_iafs <- ec_summary %>%
  select(c("Subject", iaf_cols)) %>%
  filter(Subject %in% wbehav$Subject) %>%
  pivot_longer(cols = iaf_cols, names_sep="_", 
               names_to = c("Channel", "Discard"), 
               values_to = "IAF_EC") %>%
  select(Subject, Channel, IAF_EC) %>%
  group_by(Subject) %>%
  summarize(IAF_EC = Mode(IAF_EC))
```

The mean IAF thus calculated is ```r round(mean(ec_iafs$IAF_EC),2)```, which is
exactly in the middle of our alpha frequency band.

```{r, fig.width=5, fig.height=5}
ggplot(data=ec_iafs, aes(x = IAF_EC)) +
  geom_histogram(aes(col="white"), fill="black", 
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 0.5) +
  ggtitle("IAF Distribution, Eyes Closed") + 
  xlab("IAF (Eyes Closed)") +
  theme_pander()
```

### Average Spectrum

First, let's examine the spectrograms for each channel to make sure it looks 
normal and the IAF (individual alpha frequency) is reasonable, i.e., the same 
for each channel and roughly in the middle of the band definition

```{r fig.width=6, fig.height=6}
l_ec_spectra <- pivot_longer(ec_spectra, cols=names(ec_spectra)[3:130], 
                         names_to="Freq")

l_ec_spectra <- l_ec_spectra %>% 
  rename(Power = value) %>%
  add_column(Recording = "Eyes Closed")

l_ec_spectra$Frequency <- as.numeric(substr(l_ec_spectra$Freq, 
                                            0, str_length(l_ec_spectra$Freq) -2))
l_ec_spectra <- l_ec_spectra %>% 
  add_column(Band="Delta", BandMin=0, BandMax=4) 
l_ec_spectra$Band[l_ec_spectra$Frequency <= 40] <- "Gamma"
l_ec_spectra$Band[l_ec_spectra$Frequency < 30] <- "High Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 18] <- "Upper Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 15] <- "Low Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 13] <- "Alpha"
l_ec_spectra$Band[l_ec_spectra$Frequency < 8] <- "Theta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 4] <- "Delta"

l_ec_spectra$BandMin[l_ec_spectra$Frequency <= 40] <- 30
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 30] <- 18
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 18] <- 15
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 15] <- 13
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 13] <- 8
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 8] <- 4
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 4] <- 0

l_ec_spectra$BandMax[l_ec_spectra$Frequency <= 40] <- 40
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 30] <- 30
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 18] <- 18
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 15] <- 15
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 13] <- 13
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 8] <- 8
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 4] <- 4

l_ec_spectra$Band <- factor(l_ec_spectra$Band, 
                        levels = c("Delta", "Theta", "Alpha", 
                                   "Low Beta", "Upper Beta",
                                   "High Beta", "Gamma"))
```

Now, we remove the three participants for which we have poor quality data 

```{r}
l_ec_spectra <- l_ec_spectra %>%
  filter(Subject %in% behav$Subject)
```

We can visualize the power spectra to ensure that our data looks normal

```{r}
gd <- l_ec_spectra %>% 
        group_by(Band) %>% 
        summarise(
          Min = mean(BandMin),
          Max = mean(BandMax),
          Power =mean(Power),
          Frequency = mean(BandMin)
          #Channel=
        )

ggplot(data=l_ec_spectra, aes(x=Frequency, y=Power, Channel)) +
  geom_rect(data = gd, aes(xmin = Min, xmax = Max, fill = Band), 
            ymin=0, ymax=Inf, colour=NA, alpha=0.5) +
  stat_summary(fun.data=mean_sdl, 
               geom = "ribbon", colour = "white", alpha = 0.4) +
  stat_summary(fun = mean, geom = "line", lwd = 1) +
  facet_wrap(~ Channel, ncol=4) +
  scale_alpha_manual(values = seq(0.1, 0.9, 0.1)) +
  ggtitle("Eyes-Closed Power Spectrum Across Channels") +
  ylab("Log Power") + 
  xlab("Frequency (Hz)") + 
  theme_pander() +
  coord_cartesian(xlim=c(1,40), ylim=c(5, 17)) + #ylim=c(5,17)) +
  scale_fill_brewer(palette = "Set3") +
  theme(plot.title = element_text(hjust = 0.5))
```

Now, let's create a mean spectrogram with annotations that mark the different
frequency bands. It serves as a visual reference that our data is not grossly 
out of whack.

```{r fig.width=6, fig.height=3}
K = brewer.pal(7, "Set3")

al_ec_spectra <- aggregate(l_ec_spectra[c("Power")],
                       list(Subject = l_ec_spectra$Subject, 
                            Frequency = l_ec_spectra$Frequency, 
                            Band = l_ec_spectra$Band), 
                       mean)

ggplot(data=al_ec_spectra, aes(x=Frequency, y=Power, col=Power)) +
  annotate("rect", xmin = 0, xmax = 4, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[1]) +
  annotate("rect", xmin = 4, xmax = 8, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[2]) +
  annotate("rect", xmin = 8, xmax = 13, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[3]) +      
  annotate("rect", xmin = 13, xmax = 15, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[4]) +
  annotate("rect", xmin = 15, xmax = 18, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill = K[5]) +
  annotate("rect", xmin = 18, xmax = 30, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[6]) +
  annotate("rect", xmin = 30, xmax = 40, ymin = 0, ymax = Inf,
           alpha = 0.5, fill=K[7]) +
  stat_summary(fun.data = mean_sdl, geom = "ribbon", 
               alpha = 0.5, fill = "grey50", col = "white") +
  stat_summary(fun=mean, geom="line", lwd=2) +
  coord_cartesian(xlim=c(0,40), ylim=c(5,18)) +
  xlab("Frequency (Hz)") +
  theme_pander() +
  annotate("text", x=2, y=5, label="Delta", angle=90, hjust=0) +
  annotate("text", x=6, y=5, label="Theta", angle=90, hjust=0) +
  annotate("text", x=10.5, y=5, label="Alpha", angle=90, hjust=0) +
  annotate("text", x=14, y=5, label="Low Beta", angle=90, hjust=0) +
  annotate("text", x=16.5, y=5, label="Upper Beta", angle=90, hjust=0) +
  annotate("text", x=24, y=5, label="High Beta", angle=90, hjust=0) +
  annotate("text", x=35, y=5, label="Gamma", angle=90, hjust=0) 
```

## Correlations with Rate of Forgetting

Having ensured that our EEG recordings look normal, reflect known 
neurophysiological characteristics, and that our frequency bands boundaries
are reasonable, we can proceed with correlating the rates of forgetting with
resting state QEEG characteristics.

First, we calculate the mean power for every band and channel, remove 
participants who were not included in the behavioral tests, and correct for 
multiple comparisons by frequency band: 

```{r}
wbehav <- wbehav %>%
  mutate(GlobalAlpha=(Vocabulary + Maps)/2)

Adata_ec <- l_ec_spectra %>% 
  group_by(Subject, Channel, Band, Recording) %>% 
  summarize(Power=mean(Power))

Fdata_ec <- inner_join(Adata_ec, behav, by="Subject")

Rdata_ec <- Fdata_ec %>%
  group_by(material, Channel, Band, Recording) %>%
  summarise(r = cor(Power, alpha),
            p = cor.test(Power, alpha)$p.value)

Rdata_ec <- Rdata_ec %>%
  group_by(material, Band, Recording) %>%
  mutate(q = p.adjust(p, method="BH")) 
```

Here is the full table of statistics:

```{r}
write_csv(Rdata_ec, "correlations_channel_band_eyes_closed.csv", col_names = T)

Rdata_ec %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

And here is the corresponding distributions of _r_ and _p_ values.
The dashed lines correspond to a significant threshold of _p_  < .05 on 
either a _r_ or a _p_-value scale. 

```{r}
ggplot(Rdata_ec, aes(x = Band, y = r, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  facet_wrap(~ material) +
  ggtitle("Correlation by Frequency Band") +
  ylab("r value") +
  annotate("segment", x=-Inf, xend=Inf, y=0.28, yend = 0.28, lty=2) +
  annotate("segment", x=-Inf, xend=Inf, y=-0.28, yend = -0.28, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme_pander()

ggplot(Rdata_ec, aes(x = Band, y = p, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", alpha=0.5, geom = "errorbar") +
  ggtitle("p-value, by Band") +
  ylab("p value") +
  facet_wrap(~ material) +
  annotate("segment", x=-Inf, xend=Inf, y=0.05, yend = 0.05, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  scale_y_log10() +
  theme_pander()
```

From the original list, we can extract only those channels that survive 
the FDR correction:

```{r}
survivors <- Rdata_ec %>%
  filter(q < 0.05)
  
survivors %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

Only one channel P8, in one frequency band (low beta) survives FDR correction.
We can plot the correlation between power in P8 and rate of forgetting to see
the relationship:

```{r fig.width=5, fig.height=5 }

focus <- Fdata_ec %>%
  filter(material == "Vocabulary", 
         Channel %in% survivors$Channel,
         Band == "Low Beta") %>%
  rename(Alpha = alpha)

p <- ggplot(focus, aes(x = Alpha, y = Power)) +
  geom_point(size = 4, alpha = 0.5, col = "black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() +   
  scale_y_continuous() + 
  ggtitle("Rate of Forgetting\nand Eyes-Closed Beta Power") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab("Low Beta (13-15 Hz) Power Over P8") +
  geom_text(data=survivors, col="red",
            mapping=aes(x=0.25, y=10.5, 
                        label= paste("r =", round(r, 2)))) +
  theme(plot.title = element_text(hjust = 0.5))

ggMarginal(p, col = "white", fill = "black", alpha = 0.5,
           type = "hist", bins = 13)
```

# Eyes-Open EEG Data

## (Optional) Preprocess the Data

Set the `process_raw` variable to `TRUE` to reprocess the raw EEG data.
That might take a __very__ long time. 

```{r}
process_raw = F
if (process_raw) {
  setwd("./data/eyes_open/")
  for (sub in dir()[grep("A[1-9]", dir())]) {
    setwd(sub)
    analyze.logfile(sub, "open")
    setwd("..")
  }
  setwd("../..")
}
```


## Load the Eyes-Open Spectra and Summary Files

As in the case of eyes-closed data, we are going to examine the spectral and
IAF characteristics of eyes-open recordings though the _summary_ files produced 
by preprocessing.

```{r}
eo_spectra <- NULL
for (sub in dir("data/eyes_open/")[grep("A[1-9]", 
                                        dir("data/eyes_open/"))]) {
  table <- read_tsv(paste("data/eyes_open/", sub, "/", 
                          sub, "_open_spectra.txt", sep=""),
                    col_types = cols())
  if (is.null(eo_spectra)) {
    eo_spectra <- table
  } else {
    eo_spectra <- eo_spectra %>% bind_rows(table)
  }
}

eo_summary <- NULL
for (sub in dir("data/eyes_open/")[grep("A[1-9]", 
                                        dir("data/eyes_open/"))]) {
  table <- read_tsv(paste("data/eyes_open/", sub, "/", 
                          sub, "_open_summary.txt", sep=""),
                    col_types = cols())
  if (is.null(eo_summary)) {
    eo_summary <- table
  } else {
    eo_summary <- eo_summary %>% bind_rows(table)
  }
}
```

### Individual Alpha frquency

Since eyes-open recordings are somewhat less common, we need to run a few sanity
checks on them. The first and most obsvious concerns the IAF, and whether the 
IAFs during eyes-open recordings are similar and correlated to those observed
during eyes-closed sessions. 

```{r}
eo_iafs <- eo_summary %>%
  select(c("Subject", iaf_cols)) %>%
  filter(Subject %in% wbehav$Subject) %>%
  pivot_longer(cols = iaf_cols, names_sep="_", 
               names_to = c("Channel", "Discard"), 
               values_to = "IAF_EO") %>%
  select(Subject, Channel, IAF_EO) %>%
  group_by(Subject) %>%
  summarize(IAF_EO = Mode(IAF_EO))

iafs <- inner_join(eo_iafs, ec_iafs, by="Subject")
```

The mean IAF thus calculated is ```r round(mean(eo_iafs$IAF_EO),2)```, which is,
once more, exactly in the middle of our alpha frequency band.

```{r, fig.width=5, fig.height=5}
ggplot(data=eo_iafs, aes(x = IAF_EO)) +
  geom_histogram(aes(col="white"), fill="black", 
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 0.5) +
  ggtitle("IAF Distribution, Eyes-Open") + 
  xlab("IAF (Eyes Open)") +
  theme_pander()
```

As it can be seen, there is a significant correlation in the IAFs between the
two recordings (_p_ = ```r round(cor.test(iafs$IAF_EO, iafs$IAF_EC)$p.value, 6)```).
Between recordings, the individual IAF remains the 
same +/- ```r round(sd(iafs$IAF_EO - iafs$IAF_EC), 2)```.


```{r fig.width=6, fig.height=5}
ggplot(iafs, aes(x=IAF_EO, y=IAF_EC)) +
  geom_count(alpha=0.5, col="black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  xlab("IAF Eyes Open (Hz)") +
  ylab("IAF Eyes Closed (Hz") +
  geom_text(data = summarize(iafs, Mean=mean(IAF_EC)),  # Avoid label overlaps
            x= 9.5, y=11.5, 
            col="red", 
            label=paste("r =", 
                        round(cor(iafs$IAF_EC, iafs$IAF_EO), 2))) +
  ggtitle("Correlation Between Eyes-Closed and Eyes-Open IAF") +
  theme_pander()
```

### Average Spectrum

We can also investigate whether the average spectrum is comparable to that of 
eyes-closed recordings. To do so, we will calculate the average spectrum across
channels, much in the same ways as it was done for eyes-closed data:

```{r}
l_eo_spectra <- pivot_longer(eo_spectra, cols=names(ec_spectra)[3:130], 
                         names_to="Freq")

l_eo_spectra <- l_eo_spectra %>% 
  rename(Power = value) %>%
  add_column(Recording = "Eyes Open")

l_eo_spectra$Frequency <- as.numeric(substr(l_eo_spectra$Freq, 
                                            0, str_length(l_eo_spectra$Freq) -2))
l_eo_spectra <- l_eo_spectra %>% 
  add_column(Band="Delta", BandMin=0, BandMax=4) 
l_eo_spectra$Band[l_eo_spectra$Frequency <= 40] <- "Gamma"
l_eo_spectra$Band[l_eo_spectra$Frequency < 30] <- "High Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 18] <- "Upper Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 15] <- "Low Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 13] <- "Alpha"
l_eo_spectra$Band[l_eo_spectra$Frequency < 8] <- "Theta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 4] <- "Delta"

l_eo_spectra$BandMin[l_eo_spectra$Frequency <= 40] <- 30
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 30] <- 18
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 18] <- 15
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 15] <- 13
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 13] <- 8
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 8] <- 4
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 4] <- 0

l_eo_spectra$BandMax[l_eo_spectra$Frequency <= 40] <- 40
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 30] <- 30
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 18] <- 18
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 15] <- 15
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 13] <- 13
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 8] <- 8
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 4] <- 4

l_eo_spectra$Band <- factor(l_eo_spectra$Band, 
                        levels = c("Delta", "Theta", "Alpha", 
                                   "Low Beta", "Upper Beta",
                                   "High Beta", "Gamma"))
```

Now, we remove the three participants for which we have poor quality data 

```{r}
l_eo_spectra <- l_eo_spectra %>%
  filter(Subject %in% behav$Subject)
```

And we can visualize the eyes-open power spectra by channel:

```{r}
gd <- l_eo_spectra %>% 
        group_by(Band) %>% 
        summarise(
          Min = mean(BandMin),
          Max = mean(BandMax),
          Power =mean(Power),
          Frequency = mean(BandMin)
        )

ggplot(data=l_eo_spectra, aes(x=Frequency, y=Power, Channel)) +
  geom_rect(data = gd, aes(xmin = Min, xmax = Max, fill = Band), 
            ymin=0, ymax=Inf, colour=NA, alpha=0.5) +
  stat_summary(fun.data=mean_sdl, 
               geom = "ribbon", colour = "white", alpha = 0.4) +
  stat_summary(fun = mean, geom = "line", lwd = 1) +
  facet_wrap(~ Channel, ncol=4) +
  scale_alpha_manual(values = seq(0.1, 0.9, 0.1)) +
  ggtitle("Eyes-Open Power Spectrum Across Channels") +
  ylab("Log Power") + 
  xlab("Frequency (Hz)") + 
  theme_pander() +
  coord_cartesian(xlim=c(1,40), ylim=c(5, 17)) + #ylim=c(5,17)) +
  scale_fill_brewer(palette = "Set3") +
  theme(plot.title = element_text(hjust = 0.5))
```

At a first sight, the eyes-open power spectrum looks remarkably similar, with
the obvious (and expected) difference of a less prominent pean in the alpha 
band. To better compare the two types of recordings, we can visualize the 
avearage spectra on top of ech other:

```{r fig.width=6, fig.height=3}
l_spectra <- rbind(l_ec_spectra, l_eo_spectra)

al_spectra <- l_spectra %>%
  group_by(Subject, Frequency, Band, Recording) %>%
  summarise(Power = mean(Power))

al_spectra <- al_spectra %>% ungroup()

ggplot(data=al_spectra, aes(x = Frequency, y = Power, 
                            col = Recording, fill = Recording)) +
  annotate("rect", xmin = 0, xmax = 4, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[1]) +
  annotate("rect", xmin = 4, xmax = 8, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[2]) +
  annotate("rect", xmin = 8, xmax = 13, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[3]) +      
  annotate("rect", xmin = 13, xmax = 15, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[4]) +
  annotate("rect", xmin = 15, xmax = 18, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill = K[5]) +
  annotate("rect", xmin = 18, xmax = 30, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[6]) +
  annotate("rect", xmin = 30, xmax = 40, ymin = 0, ymax = Inf,
           alpha = 0.5, fill=K[7]) +
  stat_summary(fun.data = mean_sdl, geom = "ribbon", 
               alpha = 0.5, col = "white") +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  stat_summary(fun = mean, geom = "line", lwd = 2) +
  coord_cartesian(xlim=c(0,40), ylim=c(5,18)) +
  xlab("Frequency (Hz)") +
  ylab("Log Power") +
  theme_pander() +
  annotate("text", x=2, y=18, label="Delta", angle=90, hjust=1) +
  annotate("text", x=6, y=18, label="Theta", angle=90, hjust=1) +
  annotate("text", x=10.5, y=18, label="Alpha", angle=90, hjust=1) +
  annotate("text", x=14, y=18, label="Low Beta", angle=90, hjust=1) +
  annotate("text", x=16.5, y=18, label="Upper Beta", angle=90, hjust=1) +
  annotate("text", x=24, y=18, label="High Beta", angle=90, hjust=1) +
  annotate("text", x=35, y=18, label="Gamma", angle=90, hjust=1) 
```

## Correlations with Rate of Forgetting

Like we did for eyes-closed data, we can now calculate the correlations between
eyes-open EEG power spectra and the behavioral rate fo forgetting for verbal and
visual materials. 

```{r}
wbehav <- wbehav %>%
  mutate(GlobalAlpha=(Vocabulary + Maps) / 2)

Adata_eo <- l_eo_spectra %>% 
  group_by(Subject, Channel, Band, Recording) %>% 
  summarize(Power=mean(Power))

Fdata_eo <- inner_join(Adata_eo, behav, by = "Subject")

Cdata_eo <- Fdata_eo %>%
  group_by(material, Channel, Band, Recording) %>%
  summarise(r = cor(Power, alpha),
            p = cor.test(Power, alpha)$p.value)

Rdata_eo <- Cdata_eo %>%
  group_by(material, Band, Recording) %>%
  mutate(q = p.adjust(p, method="BH")) 
```

We save the data so that it can be visualized as a topomap with Python MNE. 
Here is the full table of statistics:

```{r}
write_csv(Rdata_eo, "correlations_channel_band_eyes_open.csv", col_names = T)

Rdata_eo %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

And here is the corresponding distributions of _r_ and _p_ values. 
The dashed lines correspond to a significant threshold of _p_  < .05 on either 
a _r_ or a _p_-value scale. 

```{r}
ggplot(Rdata_eo, aes(x = Band, y = r, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  facet_wrap(~ material) +
  ggtitle("Correlation, by Band") +
  ylab("r value") +
  annotate("segment", x=-Inf, xend=Inf, y=0.28, yend = 0.28, lty=2) +
  annotate("segment", x=-Inf, xend=Inf, y=-0.28, yend = -0.28, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme_pander()

ggplot(Rdata_eo, aes(x = Band, y = p, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  ggtitle("p-value, by Band") +
  ylab("p value") +
  facet_wrap(~ material) +
  annotate("segment", x=-Inf, xend=Inf, y=0.05, yend = 0.05, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  scale_y_log10() +
  theme_pander()
```

Now, we might wonder how much the correlations change from eyes closed to
eyes open recordings. The best way to chack is to visualize then as a 
scatterplot:

```{r}
Rdata <- bind_rows(Rdata_ec, Rdata_eo)
wrdata <- pivot_wider(Rdata, id_cols = c("material", "Channel", "Band"),
                      names_from = Recording, 
                      values_from = r)

ggplot(wrdata, aes(x = `Eyes Open`, y = `Eyes Closed`, 
                   col = material)) +
  annotate("segment", x = 0, y=-Inf, xend=0, 
           yend=Inf, col="grey", lwd=1, lty=1) +
  annotate("segment", x = -Inf, y=0, xend=Inf, 
           yend=0, col="grey", lwd=1, lty=1) +
  
  geom_point(alpha=0.5, size=3) +
  annotate("segment", x = -0.25, y=-0.25, xend=0.5, 
           yend=0.5, col="black", lwd=1, lty=2) +
  theme(plot.title = element_text(hjust = 0.5)) +
  #geom_smooth(method = "lm", formula = y ~ x, 
  #            fullrange = T, lwd = 1) +
  scale_color_brewer(palette="Dark2") +
  ggtitle("Correlation Coefficients Across Recordings") +
  theme_pander()
```

Let's tabulate only those channels that survive the FDR correction:

```{r}
survivors <- Rdata_eo %>%
  filter(q < 0.05)
  
survivors %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```


In the case of eyes-open data, we have a much larger set of channels that
survive correction for multiple comparison. Although these correlations 
generally reflect the same trends we have seen the eyes-closed data (i.e., 
significant positive correlations for prefrontal and right parietal sites), they
are now stroger for _visual_ (not _verbal_) materials and centered at at 
slightly higher frequency (upper beta instead of low beta). 

Here is an overview of the scatterplots and correlations:

```{r fig.width=9, fig.height=9}
focus <- Fdata_eo %>%
  filter(material == "Maps", 
         Channel %in% survivors$Channel,
         Band == "Upper Beta") %>%
  rename(Alpha = alpha)

ggplot(focus, aes(x=Alpha, y=Power)) +
  geom_point(size=4, alpha=0.5, col="black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() +   
  scale_y_continuous() + 
  ggtitle("Rate of Forgetting and Eyes-Open Beta Power") +
  xlab(expression(paste(alpha, " Maps"))) +
  ylab("Upper Beta Power (15-18 Hz)") +
  facet_wrap(~ Channel, scales="free_y") +
  geom_text(data=survivors, col="red",
            mapping=aes(x=0.25, y=9.5, 
                        label= paste("r =", round(r, 2)))) +
  theme(panel.spacing = unit(1.5, "lines")) +
  theme(plot.title = element_text(hjust = 0.5))
```