---
title: Gammaband responses to equiluminant colors in early human visual cortex
output:
  html_document:
    toc: yes
    df_print: paged
  html_notebook:
    toc: yes
    number_sections: true
  pdf_document:
    toc: yes
editor_options:
  chunk_output_type: inline
---

# Introduction
This notebook contains statistical analyses reported in Stauch, Peter, Ehrlich, Nolte, & Fries (2021), _Gammaband responses to equiluminant colors in early human visual cortex._

To run it, the dataset color_dataset.csv [(available from Zenodo)](https://doi.org/ADDLINK) needs to be stored in the same folder as the notebook.

First, we'll load in required libraries as well as the dataset, remove invalid trials and prepare data formats.

```{r,warning=FALSE,message=FALSE}
library(magrittr)
library(lme4)
library(lmerTest)
library(afex)
library(dplyr)
library(data.table)
library(tidyverse)
library(Hmisc)
library(emmeans)
library(tidymodels)
library(MuMIn)
library(BayesFactor)
data <- fread('color_dataset.csv', stringsAsFactors=TRUE)
data$ID <- factor(data$ID)
data$stimulusVerbal <- relevel(data$stimulusVerbal, ref = "grating")
data %<>% filter(answer !='prepress' & answer !='pause')

peaks <- fread('color_peaks.csv', stringsAsFactors=TRUE)

```


# Results
We'll now go through the subsections of our Results and provide code for statistical analyses.

## Behavior

### Get average per-subject accuracy and reaction times.
```{r,warning=FALSE,message=FALSE}
data %>% 
    filter(correctReport==1|correctReport==0) %>%
    filter(wasCatch==0 & wasFixated==1) -> behavData
behavData %>% 
    group_by(ID, stimType) %>%
    dplyr::summarize(meanAcc = mean(correctReport, na.rm=TRUE),
                     meanRT = mean(rt, na.rm=TRUE)) -> meanBehav
behavData %>% 
    group_by(ID, stimulusVerbal) %>%
    dplyr::summarize(meanAcc = mean(correctReport, na.rm=TRUE),
                     meanRT = mean(rt, na.rm=TRUE)) -> meanBehavAllStims
behavData %>% 
    group_by(ID, stimulusVerbal) %>%
    filter(repNr > 50) %>% # repetitions during which staircases were stable
    dplyr::summarize(meanContrast = mean(tContrast, na.rm=TRUE)) ->
                       meanBehavAllStimsStable
```

### Get average RT.
```{r,warning=FALSE,message=FALSE}
meanBehav %>%
    group_by(stimType) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanRT*1000))))
```

### Get average accuracy.
```{r,warning=FALSE,message=FALSE}
meanBehav %>%
    group_by(stimType) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanAcc))))
```

### Compare reaction times over stimuli.
```{r,warning=FALSE,message=FALSE}
meanBehavAllStims %>%
    group_by(stimulusVerbal) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanRT*1000))))
meanBehavAllStims %>% ungroup() %>%
    filter(stimulusVerbal != "grating") -> dat

m0 <- aov_ez(id = "ID", dv = "meanRT", within = "stimulusVerbal",
              data = dat)
summary(m0)
emmeans(m0, specs = pairwise~stimulusVerbal, adjust ='holm')
```

### Compare staircased target contrast over stimuli.
```{r,warning=FALSE,message=FALSE}
meanBehavAllStimsStable %>%
    group_by(stimulusVerbal) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanContrast*100))))
meanBehavAllStimsStable %>% ungroup() %>%
    filter(stimulusVerbal != "grating") -> dat

m1 <- aov_ez(id = "ID", dv = "meanContrast", within = "stimulusVerbal",
              data = dat)
summary(m1)
m1
emmeans(m1, specs = pairwise~stimulusVerbal, adjust ='holm')
```

### Test blue target contrast against all other colors.
```{r,warning=FALSE,message=FALSE}
meanBehavAllStimsStable %>%
  filter(stimulusVerbal != 'grating') %>%
  mutate(isBlue = factor(stimulusVerbal=='blue')) %>%
  group_by(isBlue, ID) %>%
  dplyr::summarize(meanContBN = mean(meanContrast)) %>%
  group_split(isBlue) -> s1

t.test(s1[[1]]$meanContBN, s1[[2]]$meanContBN, paired = TRUE)
```

## ERFs
### Get averages over stimulus types
```{r,warning=FALSE,message=FALSE}
data %>% 
    filter(wasFixated==1) -> erfData
erfData %>% 
    group_by(ID, stimType) %>%
    dplyr::summarize(meanC1 = mean(C1, na.rm=TRUE)) %>%
    ungroup() -> meanERFType
erfData %>% 
    group_by(ID, stimulusVerbal) %>%
    dplyr::summarize(meanC1 = mean(C1, na.rm=TRUE)) %>%
    ungroup() -> meanERFStim
```

### Compare between stimuli
```{r,warning=FALSE,message=FALSE}
meanERFStim %>%
    group_by(stimulusVerbal) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanC1))))
meanERFStim %>% ungroup() %>%
    filter(stimulusVerbal != "grating") -> dat

m2 <- aov_ez(id = "ID", dv = "meanC1", within = "stimulusVerbal",
              data = dat)
summary(m2)
m2
emmeans(m2, specs = pairwise~stimulusVerbal, adjust ='holm')

```

### Test grating ERF vs color ERF.
```{r,warning=FALSE,message=FALSE}
meanERFType %>%
  group_split(stimType) -> s1

t.test(s1[[1]]$meanC1, s1[[2]]$meanC1, paired = TRUE)
```
## Gamma power and frequency
### Get average per-subject gamma power and frequency per stimulus
```{r,warning=FALSE,message=FALSE}
data %>% 
    filter(wasFixated==1) -> gammaData
gammaData %>% 
    group_by(ID, stimType) %>%
    dplyr::summarize(meanGPower = mean(percGammapower, na.rm=TRUE),
                     meanGFreq = mean(realGammapeak, na.rm=TRUE)) %>%
    ungroup() -> meanGP
gammaData %>% 
    group_by(ID, stimulusVerbal) %>%
    dplyr::summarize(meanGPower = mean(percGammapower, na.rm=TRUE),
                     meanGFreq = mean(realGammapeak, na.rm=TRUE)) %>%
    ungroup() -> meanGPStims
peaks %>%
  group_by(ID, stimType) %>%
  dplyr::summarize(meanPeak1 = mean(peak1, na.rm=TRUE),
                   meanPeak2 = mean(peak2, na.rm=TRUE)) %>%
  ungroup() -> peaksStimtype
```

### Compare gamma power between stimulus types
```{r,warning=FALSE,message=FALSE}
meanGP %>%
    group_by(stimType) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanGPower))))

meanGP %>%
  group_split(stimType) -> s1

t.test(s1[[1]]$meanGPower, s1[[2]]$meanGPower, paired = TRUE)
```

### Compare gamma frequency between stimulus types
```{r,warning=FALSE,message=FALSE}
peaksStimtype %>%
    group_by(stimType) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanPeak1, na.rm=TRUE))))

peaksStimtype %>%
  group_split(stimType) -> s1

t.test(s1[[1]]$meanPeak1, s1[[2]]$meanPeak1, paired = TRUE)
cor.test(s1[[1]]$meanPeak1, s1[[2]]$meanPeak1)

peaksStimtype %>%
    group_by(stimType) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanPeak2, na.rm=TRUE)))) 
t.test(s1[[1]]$meanPeak2, s1[[2]]$meanPeak2, paired = TRUE)
```
### Compare gamma power between stimuli
```{r,warning=FALSE,message=FALSE}
meanGPStims %>%
    group_by(stimulusVerbal) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$meanGPower))))
meanGPStims %>% ungroup() %>%
    filter(stimulusVerbal != "grating") -> dat

m3 <- aov_ez(id = "ID", dv = "meanGPower", within = "stimulusVerbal",
              data = dat)
summary(m3)
m3
emmeans(m3, specs = pairwise~stimulusVerbal, adjust ='none')
```
### Compare gamma power between red and green - Bayes factor
```{r,warning=FALSE,message=FALSE}
meanGPStims %>% ungroup() %>%
    filter(stimulusVerbal == "red" | stimulusVerbal == "green") %>%
    group_split(stimulusVerbal)-> s

1/ttestBF(s[[1]]$meanGPower, s[[2]]$meanGPower, paired = TRUE, rscale = 1)
```
### Check how many peaks were fitted for all stimuli
```{r,warning=FALSE,message=FALSE}
peaks %>% ungroup() %>%
  group_by(stimType,ID) %>%
  dplyr::summarize(onePeak1Found = paste("Peak1", 
                                         as.character(sum(peak1found)>0)),
                   onePeak2Found = paste("Peak2", 
                                         as.character(sum(peak2found)>0))) -> d
prop.table(table(d$stimType, d$onePeak1Found),1)
infer::prop_test(d, response = onePeak1Found, explanatory = stimType, 
                 success = "Peak1 TRUE")

prop.table(table(d$stimType, d$onePeak2Found),1)
infer::prop_test(d, response = onePeak2Found, explanatory = stimType, 
                 success = "Peak2 TRUE")
```

### Compare the lower gamma peak frequency between stimuli
```{r,warning=FALSE,message=FALSE}
peaks %>%
    group_by(stimulusVerbal) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$peak1))))
peaks %>% ungroup() %>%
    filter(stimulusVerbal != "grating") -> dat

m4 <- aov_ez(id = "ID", dv = "peak1", within = "stimulusVerbal",
              data = dat)
summary(m4)
m4
emmeans(m4, specs = pairwise~stimulusVerbal, adjust ='holm')
```

### Compare the higher gamma peak frequency between stimuli
As the second peak is missing for many participant-stimulus combinations, we're 
using a linear mixed-model instead of a RMANOVA here.
```{r,warning=FALSE,message=FALSE}
peaks %>%
    group_by(stimulusVerbal) %>%
    do(data.frame(rbind(Hmisc::smean.cl.boot(.$peak2))))
peaks %>% ungroup() %>%
    filter(stimulusVerbal != "grating") -> dat

m5 <- lmer(peak2 ~ (1|ID) + stimulusVerbal, data = dat)
anova(m5)
emmeans(m5, specs = pairwise~stimulusVerbal, adjust ='holm')
```

## Methods
### Number of excluded trials
```{r,warning=FALSE,message=FALSE}
nanCounts <- data %>% 
    group_by(ID) %>%
    dplyr::summarize(non_na_count = sum(!is.na(gammapower)))
sprintf('Gamma power reported in %.0f%% of trials.', 
        100*(mean(nanCounts$non_na_count)/540))
```