Analysis_RR.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)
library(robcor)

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")
d3pal <- pal_jco()
K = d3pal(7)


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}
r2p <- function(r, n) {
  df = n - 2
  t = (r * sqrt(df)) / sqrt(1 - r**2)
  2 * pt(-abs(t), df = df)
}

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),
            rob_r =robcor(Power, alpha, method="ssd"),
            p = cor.test(Power, alpha)$p.value)


Rdata_ec <- Rdata_ec %>%
  group_by(material, Band, Recording) %>%
  mutate(q = p.adjust(p, method="BH"),
         rob_p = r2p(rob_r, n=50))  %>%
  mutate(rob_q = p.adjust(rob_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(rob_q < 0.05 | 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("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}
iafplot <- ggplot(iafs, aes(x=IAF_EO, y=IAF_EC)) +
  geom_count(alpha=0.5, col="black") + 
  scale_size_area() +
  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("IAF Across Recordings") +
  theme_pander()

iafplot
```

### 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 Across Recordings

```{r, fig.width=10, fig.height=4}
test <- al_spectra %>%
  group_by(Subject, Band, Recording) %>%
  summarize(Power = mean(Power)) %>%
  ungroup()

test$Recording <- recode(test$Recording, 
                         `Eyes Open` = "EyesOpen", 
                         `Eyes Closed` = "EyesClosed")

by_subj <- test %>%
  pivot_wider(id_cols = c("Subject", "Band"), 
              values_from = c("Power"), 
              names_from = c("Recording")) 

global <- by_subj %>%
  group_by(Band) %>%
  summarize(r=cor(EyesOpen, EyesClosed), 
            peo=mean(EyesOpen),
            pec=mean(EyesClosed))

global$y <- seq(8,11, 0.5)

bandplot <- ggplot(by_subj, aes(y=EyesClosed, x=EyesOpen, fill=Band, color=Band)) +
  geom_point(stat="identity") +
  geom_smooth(method="lm") +
  #scale_fill_brewer(palette="Set3") +
  #scale_color_brewer(palette="Set3") +
  scale_color_jco() +
  scale_fill_jco() +
  xlab("Power Eyes Open") +
  ylab("Power Eyes Closed") +
  geom_text(data = global, 
            x = global$peo + 1.5, y=global$pec -.5,
            label=paste("r =", 
                        round(global$r, 2)),
            show.legend = F) +
  ggtitle("(B) Spectral Power Across Recordings") +
  theme_pander() +
  theme(legend.position="bottom") 


ggplot(global, aes(x=Band, y=r, fill=Band, color=Band)) +
  geom_bar(stat="identity", alpha=0.75) +
  scale_fill_brewer(palette="Set3") +
  scale_color_brewer(palette="Set3") +
  ylim(0,1) +
  theme_pander() 

iafplot <- iafplot + theme(legend.position="bottom") +   
  ggtitle("(A) IAF Across Recordings") 

grid.arrange(iafplot, bandplot, ncol=2)
g <- arrangeGrob(iafplot, bandplot, ncol=2) #generates g
ggsave(file="figure5.png", g, dpi=300, width = 9, height=5, units = "in")
```

## 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,
            rob_r = robcor(Power, alpha))

Rdata_eo <- Cdata_eo %>%
  group_by(material, Band, Recording) %>%
  mutate(q = p.adjust(p, method="BH")) %>%
  mutate(rob_p = r2p(rob_r, n=50)) %>%
  mutate(rob_q = p.adjust(rob_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 | rob_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))
```

# Summary

In general, posterior (P8) and frontal channels (AF3/4, F7/8) are the most commonly associated with rate of forgetting. The correlations mostly happen in the beta range.

```{r}
Rdata_g <- full_join(Rdata_ec, Rdata_eo)
ggplot(filter(Rdata_g, p < 0.05), aes(x=Band, fill=Channel)) + 
  geom_bar(col="white", alpha=0.75) +
  theme_pander() + 
  scale_fill_jco()

Rdata_g$Recording <- recode(Rdata_g$Recording, 
                            `Eyes Open` = "EyesOpen", 
                            `Eyes Closed` = "EyesClosed")


Rdata_g %>% 
  pivot_wider(id_cols = ("Channel", "Band"), 
              values_from=c("Recording")) %>%
  group_by(Recording, material)
```

# Predicting Rate of Forgetting with Group Lasso

So far, we have shown that neurophysiological signals exists that are 
correlated with the rate of forgetting. This, however, does not directly 
answer wheter we could "read" an individual's rate of forgetting from their
EEG data directy. 

## Group LASSO

To do so, we need to use a ML approach. Specifically, here we will use 
_group LASSO_ regression. Group LASSO is a way to combine regression, validation,
and feature selection into a single approach. 

LASSO works like a normal regression, with the additional constraint that the 
$\beta$ regressors need to minimize the sum of the residual sum of squares 
(RSS) and the quantity $- \lambda||\beta||_{1}$, where $||\beta||_{1}$ is the 
first-order norm (the sum of absolute values of each regressor weight). For
$\lambda = 0$, LASSO reduces to normal regression but, as the value of 
$\lambda$ grows, more and more regressors are "pushed" to zero and removed from
the pool. In _group_ LASSO, regressors are combined
and selected by groups, which is a natural choice when natural groupings exists,
as in the cases of frequency bands.

As a demo, we will pilot this approach to eyes-closed recordings and the 
combined rates of forgetting.an 

## Data Preparation

First, we need to transform the data appropriately, converting from the long 
to the wide format, which is better suited for regression: 

```{r}
library(gglasso)
wfdata <- Fdata_ec %>% 
  pivot_wider(id_cols = c("Subject", "Channel", 
                          "Band", "material", "alpha"), 
              names_from = c("Channel", "Band"), 
              values_from = c("Power"))

wide_data_ec <- wfdata %>%
  pivot_wider(names_from = c("material"), 
              values_from = "alpha")
```

From the wide format data, we are going to define the data that will enter
the regression, the matrix of observations _X_ and the vector of variables 
_Y_. The matrix of observations _X_ is 50 x 98 matrix, corresponding
to 50 participants _x_  98 possible power values (14 channels _x_ 7 frequency
bands). The _Y_ vector is made of 50 observed rates of forgetting. To make
the prediction reliable, we will average the rates of forgetting across 
materials (Maps + Vocabulary)

```{r}

channel_names <- c("AF3", "AF4", "F3", "F4", 
                   "F7", "F8", "FC5", "FC6", 
                   "T7", "T8", "P7", "P8", 
                   "O1", "O2")


band_names <- c("Delta", "Theta", "Alpha", "Low Beta", 
                "Upper Beta", "High Beta", "Gamma")

by_band_order <- as.vector(
  outer(channel_names, 
        band_names, 
        function(x, y){
          paste(x, y, sep="_")
          }))

# Predictors
X <- wide_data_ec[,c(-1, -100, -101)] %>%
  #select(ends_with("Beta")) %>%
  as.matrix()

X <- X[,by_band_order]
# We need to reorder X by frequency bands


# Maps
Ymaps <- wide_data_ec[,c(101)] %>%
  as.matrix()

# Vocab
Yvoc <- wide_data_ec[,c(101)] %>%
  as.matrix()

# Average across materials to obrtain stable estimates
Y <- (Ymaps + Yvoc)/2

# Let's save the mean (will be halpful later)
uY <- mean(Y)

# Center Y to remove dominant effect of intercept
Y <- Y - uY
```

Finally, we need to define our set of grouping variables, defined as a vector
of numbers that identify each group. We are going to create two such vectors, 
one for the freqency _bands_ and one for the _channels_. In this example
we are going to use _bands_ as the grouping variables.

```{r}
bands <- rep(1:7, each=14)
channels <- sort(rep(1:14, 7))
```

## Fitting the LASSO

Now, with that in place, we can fit the LASSO model. The LASSO model is fit
much like a normal regression, with an additional partameter $\lambda$. 

Because we do not know _a priori_ which value of $lambda$ to pick, the canonical
solution is to fit fit different models for different values of $\lambda$.

```{r}
fit <- gglasso(x = X, y = Y, 
               group = bands,
               loss = 'ls', 
               nlambda=100,
               lambda.factor = 0.001)
coef.mat=fit$beta
```

## Finding the optimal $\lambda$ with Cross-Validation

Our goal, of course, is generalizability across participants. Although the fit 
model is important, what we really care is a different measure, that of 
_cross-validation_. For cross-validation, we are going to use the _n_-fold 
approach, in which the data is randomly split so that 1/_n_ of it is saved for
test. (if _n_ is identical to the number of participants, it becomes the
"Leave-One-Out" validation scheme). Here, we used the Leave-One-Out method, 
which guarantees that the results are the same every time this script is ruin
(because there only one possible permutation of the data). So, _n_ = 50.

```{r}
fit.cv <- cv.gglasso(x = X,
                     y = Y,
                     group = bands,
                     pred.loss = "L1",
                     nlambda = 100,
                     lambda.factor = 0.001,
                     nfolds = 50)
```

### Visualizing Validation Results

Here are the fit and cross-validation plots, with the vertical lines marking
the points corresponding to $\lambda_{min}$ (optimal for cross-validation) and  
$\lambda_{1SE}$ (more conservative model). Because we have used the 
Leave-One-Out validation scheme, we are can use the the more aggressive 
$\lambda_{min}$ estimate, since our sample of the folds was exhaustive.

```{r}
# Best lambda

lmin <- fit.cv$lambda.min  # Optimal lambda, maximizes generalization
l1se <- fit.cv$lambda.1se  # Conservative model, min + 1SE

plot(fit.cv, 
     pch=20,
     lwd=2)

plot(fit, group=T, 
     col = inferno(length(unique(fit$group)) + 1), 
     lwd = 2)

legend(x="topright", 
       legend=c("Delta", "Theta", "Alpha",
                "Low Beta", "Upper Beta", "High Beta",
                "Gamma"),
       bty="n",
       ncol=2,
       col=inferno(length(unique(fit$group)) + 1),
       lty=1, lwd=2)
abline(v=log(lmin), lty=2)
abline(v=log(l1se), lty=2)
```

Using the Lambda min value, we can now compare predicted vs. observed values.

```{r}
prediction <- predict(object = fit,
                    newx = X,
                    s = lmin,
                    type='link')

observed <- data.frame(Subject = wide_data_ec$Subject,
                       RoF = as_vector(Y) + uY,
                       Condition = "Observed") %>%
  rename(RoF = Maps)

predicted <- data.frame(Subject = wide_data_ec$Subject,
                       RoF = as_vector(prediction) + uY,
                       Condition = "Predicted") %>%
  rename(RoF = X1)

comparison <- as_tibble(rbind(observed, predicted))

plt = cbind(Y,
            predict(object = fit,
                    newx = X,
                    s = lmin,
                    type='link'))

# matplot(plt,
#         main="Predicted vs Actual",
#         type='b',
#         lwd=2,
#         pch="o",
#         col=viridis(3),
#         xlab="Subject",
#         ylab="Rate of Forgetting (alpha)")

ggplot(acomparison, aes(x=reorder(Subject, RoF), y=RoF), col=Condition) +
  geom_point(aes(col=Condition)) +
  geom_line(alpha=0.5, lty=2) +
  xlab("Subject") +
  annotate("text", x= 40, y = 0.20, 
           label=(paste("r =", 
                        round(cor(Y, prediction),
                              2)))) +
  theme_pander() +
  ylab("Rate of Forgetting")+
  ggtitle("Cross-Validation:\nPredicted vs. Observed Rates of Forgetting") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) 
```

```{r}
wcomparison <- comparison %>%
  pivot_wider(names_from=Condition, values_from = RoF)

p <- ggplot(wcomparison, aes(x = Predicted, y = Observed)) +
  geom_point(alpha=0.5, size=4) +
  geom_abline(slope = 1, intercept = 0,
              lwd=1, col="blue", lty=2) +
  annotate("text", x= 0.40, y = 0.25, 
           label=(paste("r =", 
                        round(cor(Y, prediction),
                              2)))) +
  theme_pander() +
  #coord_fixed() +
  coord_equal(ratio=1) +
  ylab("Rate of Forgetting (Observed)") +
  xlab("Rate of Forgetting (Predicted)") +
  ggtitle("Cross-Validation") 

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

In terms of overall error, the root mean squared error between predicted and 
observed Rate of Forgetting is _RMSE_ 
= ```r round(sqrt(mean((prediction - Y)**2)), 3)```


And now, let's print the full table of beta regression weights (for
inspectio and analysis) at the desiored value of $\lambda_{min}$. We
can also save the weights for future analysis.

```{r}

get_channel <- function(s) {
  str_split(s, "_", 2)[[1]][1]
}

get_band <- function(s) {
  str_split(s, "_", 2)[[1]][2]
}

vget_channel <- Vectorize(get_channel)
vget_band <- Vectorize(get_band)


coef(fit.cv, 
     s=lmin) -> betas
  
bdf <- data.frame("Predictor" = rownames(betas),
                  "Beta" = as.numeric(betas)) %>%
  mutate(Channel = vget_channel(Predictor),
         Band = vget_band(Predictor)) 

bdf %>%
  filter(Beta != 0) %>%
  xtable() %>%
  kable(digits = 5) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))


write_csv(bdf, "group_lasso_weights.csv", col_names = T)
```

And here is the visualization of the corresponding weights as a topomap:

![Topomap of Beta weights](images/topo_group_lasso.png){width=9in}