---
title: "Plotting a grid of waveforms. This is adapted from some code pieces that Matt Craddock posted on his blog."
output: html_notebook
fontsize: 11pt
geometry: margin=1in
---

##Loading and formatting the data

Load some libraries.

```{r}
library(tidyverse)
library(gridExtra)
library(Hmisc)
```

* Set paths and load data

```{r set_path_to_data_and_load, results = 'hide', warning=FALSE, message=FALSE, error=FALSE}
data.base.path <- "/Users/jsprouse/Desktop/"
erp.path <- "exported ERPs"
erp.files <- file.path(data.base.path, erp.path) %>% 
  dir(full.names = TRUE)
erp.labels <- file.path(data.base.path, erp.path) %>% 
  dir() %>%
  stringr::str_replace(".txt", "")

## Load the data
erp.data.list <- map(erp.files, read_tsv)
```

Preprocessing
=============

* load the data and preprocess it (to long format)

```{r preprocess_data}
## -------------------
## Preprocess the data
## -------------------
## 1. Bind the individual subject data
## 2. Create Subject and Condition labels
## 3. Restructure the data to long format
## 4. Filter out time points that won't be used

erp.data <- bind_rows(erp.data.list) %>%
  dplyr::mutate(Subject_Cond = rep(erp.labels, 
                                   each = 32)) %>% #this is the number of channels
  tidyr::separate(col = Subject_Cond, 
           into = c("subject", "condition"), sep = "_") %>%
  tidyr::gather(key = timet, value = mV, 
         -time, -subject, -condition) %>%
  dplyr::filter(!is.na(mV)) %>%
  dplyr::rename(channel = time, time.ms = timet) %>%
  dplyr::mutate(mV = as.double(mV)) %>%
  dplyr::mutate(time.ms = as.double(time.ms)) %>%
  arrange(condition) %>%
  dplyr::mutate(channel = as.factor(channel), subject = as.factor(subject), condition = as.factor(condition))
```



##Electrode locations

We may need electrode coordinates (Cartesian) too. Here is a giant list, formed into a data frame, taken from Matt's code. These are not the same as the coordinates we use in the topographic waveform plot (taken from the erpR package).

```{r}
elec.coord = data.frame(channel=c("Fp1",	"AF7",	"AF3",	"F1",	"F3",	"F5",	"F7",	"FT7",	"FC5",	"FC3",	"FC1",	"C1",	"C3",	"C5",	"T7",	"TP7",	"CP5",	"CP3",	"CP1",	"P1",	"P3",	"P5",	"P7",	"P9",	"PO7",	"PO3",	"O1",	"Iz",	"Oz",	"POz",	"Pz",	"CPz",	"Fpz",	"Fp2",	"AF8",	"AF4",	"AFz",	"Fz",	"F2",	"F4",	"F6",	"F8",	"FT8",	"FC6",	"FC4",	"FC2",	"FCz",	"Cz",	"C2",	"C4",	"C6",	"T8",	"TP8",	"CP6",	"CP4",	"CP2",	"P2",	"P4",	"P6",	"P8",	"P10",	"PO8",	"PO4",	"O2",	"LO1",	"LO2",	"SO1",	"IO1",	"M1",	"M2", "FT9", "FT10"), x=c(-0.158507951196727,	-0.306220792940307,	-0.160894602703784,	-0.111231689343979,	-0.221253907813213,	-0.327718412610652,	-0.426753882970127,	-0.505736594237533,	-0.381956848927333,	-0.25547983525988,	-0.127945796584781,	-0.13319,	-0.26669,	-0.3999,	-0.53318,	-0.505742347629983,	-0.381951819153741,	-0.255491477354236,	-0.127956990743876,	-0.111218290253606,	-0.221267739734655,	-0.327693903532661,	-0.426737596055039,	-0.529417172201933,	-0.306206021969086,	-0.160887810154621,	-0.158542158546356,	7.75813747259848E-17,	6.20516286659972E-17,	4.65292304868045E-17,	3.10301005967956E-17,	1.55064777708038E-17,	0,	0.158507951196727,	0.305957920131658,	0.160894602703784,	0,	0,	0.111132030087498,	0.220965559664675,	0.327483020901138,	0.426495945228168,	0.505736594237533,	0.381956848927333,	0.25547983525988,	0.127945796584781,	0,	0,	0.13348,	0.26667,	0.3999,	0.53318,	0.505559973212276,	0.381951819153741,	0.255491477354236,	0.127956990743876,	0.111118643002187,	0.220979398075575,	0.327458529427387,	0.426512248599182,	0.529417172201933,	0.305972696504555,	0.160887810154621,	0.158542158546356,	-0.438655260304613,	0.438655260304613,	-0.246525910621202,	-0.315269162639131,	-0.735005330945149,	0.735005330945149, -0.735005330945149,	0.735005330945149), y=c(0.489989723879405,	0.423151810667294,	0.389182995158822,	0.255900863979949,	0.264176032935205,	0.280839046672563,	0.311045993013766,	0.164831379442813,	0.144085733011642,	0.133419721838956,	0.128281197126804,	8.1555353589218E-18,	1.63300527432304E-17,	2.44868127489513E-17,	3.26478590184693E-17,	-0.164813725804961,	-0.144099065733086,	-0.133397426509471,	-0.128270031261287,	-0.255906687706798,	-0.264164447745561,	-0.280867644429768,	-0.311068337368441,	-0.39447167981708,	-0.423162499531643,	-0.389185803240111,	-0.489978656742783,	-0.6335,	-0.50669,	-0.37994,	-0.25338,	-0.12662,	0.50669,	0.489989723879405,	0.423317241804193,	0.389182995158822,	0.37994,	0.25338,	0.255671586783967,	0.264299964892311,	0.280637327206245,	0.311382616091412,	0.164831379442813,	0.144085733011642,	0.133419721838956,	0.128281197126804,	0.12662,	0,	8.17329273750943E-18,	1.63288280964312E-17,	2.44868127489513E-17,	3.26478590184693E-17,	-0.165340139063697,	-0.144099065733086,	-0.133397426509471,	-0.128270031261287,	-0.255677405292983,	-0.264288394800371,	-0.280665904422416,	-0.311360284389113,	-0.39447167981708,	-0.423306561600138,	-0.389185803240111,	-0.489978656742783,	0.487176021789961,	0.487176021789961,	0.528676521412089,	0.61875057065737,	-0.135163946310443,	-0.135163946310443, 0.135163946310443,	0.135163946310443))
            
```

Find the electrodes in our dataset:
```{r}
coords = elec.coord %>% 
  filter(channel %in% levels(erp.data$channel)) %>%
  droplevels()
```

Combine the electrode locations to the data using left_join():

```{r}

#first we make sure that the order of the levels of each dataframe are identical. This eliminates any warnings or errors.
erp.data$channel=factor(erp.data$channel, levels=sort(levels(erp.data$channel)))

coords$channel=factor(coords$channel, levels=sort(levels(coords$channel)))

#now we combine using left_join
full.locs <- erp.data %>% left_join(coords, by = "channel")
```

##Plot separate ERPs for each condition (we will do difference waves later)
Let's plot for a 9 channel grid first, we can worry about topography later

```{r}
channelList = c("Fz", "F3", "F4", "Cz", "C3", "C4", "Pz", "P3", "P4")

full = filter(full.locs, channel %in% channelList)

full = droplevels(full)
```

Make subsets of the conditions? Since we can calculate difference waves in the script, here I will just focus on the two conditions.
```{r}
all = filter(full, condition == "congruent" | condition == "incongruent")
all = droplevels(all)
```


I like to call the subset I am working on "plotting". If we had other subsets of conditions, this step would matter more.
```{r}
plotting = all
```

Reformat the levels as necessary.

```{r}
#re-order the levels of the factors the way we want them
#reorder the channel levels
plotting$channel = factor(plotting$channel, levels=c('F3', 'Fz', 'F4', 'C3', 'Cz', 'C4', 'P3', 'Pz', 'P4'))


#reorder the condition levels
#plotting$condition = factor(plotting$condition, levels=c("incongruent", "congruent"))

```

Calculate means
```{r}
plotting.means = plotting %>%
  #group_by(subdesign, channel, condition, time.ms) %>%
  group_by(channel, condition, time.ms) %>%
  summarize(meanVoltage = mean(mV, na.rm=TRUE))
```


Plot 1: Means of both conditions, no confidence intervals.

This just sets up the plot
```{r}
p = ggplot(data=plotting.means, aes(x=time.ms, y=meanVoltage, group=condition)) +
	ylab("microvolts") + 
	xlab("time (in ms)") +
	geom_hline(aes(yintercept=0)) +
  geom_vline(aes(xintercept=0)) +
  theme_minimal() +
  theme(legend.position="bottom", legend.background = element_rect(fill="transparent", color="transparent"), legend.title=element_blank())
```

This creates the edges of boxes to draw on the plot. I like to use them to indicate when the word was on the screen (versus when the word was off the screen). You can also use them to indicate analysis windows if you are focusing on a specific time window. Here we only have one box, so the start of the first word and the start of the last word are the same. 
```{r}
startOfFirstBox = 300
startOfLastBox = 300
widthOfBox = 200
widthBetweenBoxes = 200

words = data.frame(xmin = c(seq(startOfFirstBox, startOfLastBox, by=widthOfBox+widthBetweenBoxes)), xmax=c(seq(startOfFirstBox+widthOfBox, startOfLastBox+widthOfBox, by=widthOfBox+widthBetweenBoxes)), ymin=-Inf, ymax=Inf)
```

This plots the means, with word boxes:
```{r}
m=p + 
  geom_line(aes(color = condition)) +
  facet_wrap(~channel, ncol=3) +
  scale_color_manual(values = c('#000000', '#e41a1c')) +
  #word boxes below here
  scale_fill_manual(values = c('#000000', '#e41a1c')) + 
  guides(fill="none") +
  geom_rect(data=words, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), color="grey20", alpha=.25, inherit.aes=F) + 
	scale_x_continuous(limits=c(-200,1000), breaks=c(seq(-200, 1000, by=100)), labels=c(seq(-200, 1000, by=100)))

#This command would plot negative up
#scale_y_continuous(trans = "reverse") +

#This command would remove the x-axis and grid lines
#	theme(axis.title.x = element_blank()), axis.ticks.x = element_blank(), panel.grid.major.x=element_blank(), panel.grid.minor.x=element_blank())

#This command would add words as labels instead of time on the x-axis
#scale_x_continuous(limits=c(0,3000), breaks=c(seq(150, 2650, by=500)), labels=c("the\nthe","package\nsinger", "arrived\ncoughed", "after\nafter", "the\nthe", "deadline\nsong"))  +
```



To make it at a specific size, we need the quartz() command:
```{r}
quartz(height=9, width=15)
m
ggsave("waveform.means.pdf")
```



Plot 2: means and confidence intervals 

Calculate means and confidence intervals for both conditions
```{r message=FALSE, warning=FALSE, include=FALSE, paged.print=FALSE}
plotting.ci = plotting %>%
  group_by(channel, condition, time.ms) %>%
  do(data.frame(rbind(mean_cl_boot(.$mV))))
```

Save the results of confidence interval bootstrap so we don't have to run it again.
```{r}
write.csv(plotting.ci, file="confidence.intervals.csv", row.names=F)
```


Now we can make the same plot with bootstrapped confidence intervals

This just sets up the plot
```{r}
ci = ggplot(data=plotting.ci, aes(x=time.ms, y=y, group=condition)) +
	ylab("microvolts") + 
	xlab("time (in ms)") +
	geom_hline(aes(yintercept=0)) +
  geom_vline(aes(xintercept=0)) +
  theme_minimal() +
  theme(legend.position="bottom", legend.background = element_rect(fill="transparent", color="transparent"), legend.title=element_blank())
```

We can create boxes again:
```{r}
startOfFirstBox = 300
startOfLastBox = 300
widthOfBox = 200
widthBetweenBoxes = 200

words = data.frame(xmin = c(seq(startOfFirstBox, startOfLastBox, by=widthOfBox+widthBetweenBoxes)), xmax=c(seq(startOfFirstBox+widthOfBox, startOfLastBox+widthOfBox, by=widthOfBox+widthBetweenBoxes)), ymin=-Inf, ymax=Inf)
```

Now let's plot:
```{r}
ci.complete=ci + 
  geom_ribbon( aes(ymin=ymin, ymax=ymax, fill = condition), alpha=.3) +
  geom_line(aes(color = condition)) +
  facet_wrap(~channel, ncol=3) +
  scale_color_manual(values = c('#000000', '#e41a1c')) +
  #word boxes below here
  scale_fill_manual(values = c('#000000', '#e41a1c')) + 
  guides(fill="none") +
  geom_rect(data=words, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), color="grey20", alpha=.25, inherit.aes=F) + 
	scale_x_continuous(limits=c(-200,1000), breaks=c(seq(-200, 1000, by=100)), labels=c(seq(-200, 1000, by=100)))
```

To make it at a specific size, we need the quartz() command:
```{r}
quartz(height=9, width=15)
ci.complete
ggsave("waveform.means.ci.pdf")
```



##Difference waves

Plot 3: Difference waves for all conditions

Here we can work with difference waves instead of each condition in isolation.

```{r}
#split by conditions
split.conditions = plotting %>%
  split(f=plotting$condition)

#calculate differences in voltage
N400.difference = split.conditions$incongruent$mV - split.conditions$congruent$mV

#recreate the dataset
N400.effect = data.frame(select(split.conditions$incongruent, -condition, -mV), mV = N400.difference)

```

Let's plot for a 9 channel grid first, we can worry about topography later

```{r}
channelList = c("Fz", "F3", "F4", "Cz", "C3", "C4", "Pz", "P3", "P4")

N400.effect = filter(N400.effect, channel %in% channelList)
N400.effect = droplevels(N400.effect)
```

Reformat the levels as necessary.

```{r}
#re-order the levels of the factors the way we want them
N400.effect$channel = factor(N400.effect$channel, levels=c('F3', 'Fz', 'F4', 'C3', 'Cz', 'C4', 'P3', 'Pz', 'P4'))
```

Calculate means and confidence intervals (better to do this once here than in the ggplot calls... saves time)

```{r}
N400.effect.ci = N400.effect %>%
  group_by(channel, time.ms) %>%
  do(data.frame(rbind(mean_cl_boot(.$mV))))

write.csv(N400.effect.ci, file="difference.confidence.intervals.csv", row.names=F)
```


This just sets up the plot (no colors this time)
```{r}
q = ggplot(data=N400.effect.ci, aes(x=time.ms, y=y)) +
	ylab("microvolts") + 
	xlab("time (in ms)") +
	geom_hline(aes(yintercept=0)) +
  geom_vline(aes(xintercept=0)) +
  theme_minimal() +
  theme(legend.position="bottom", legend.background = element_rect(fill="transparent", color="transparent"), legend.title=element_blank())
```

We can create boxes again:
```{r}
startOfFirstBox = 300
startOfLastBox = 300
widthOfBox = 200
widthBetweenBoxes = 200

words = data.frame(xmin = c(seq(startOfFirstBox, startOfLastBox, by=widthOfBox+widthBetweenBoxes)), xmax=c(seq(startOfFirstBox+widthOfBox, startOfLastBox+widthOfBox, by=widthOfBox+widthBetweenBoxes)), ymin=-Inf, ymax=Inf)
```


This plots means and confidence intervals
```{r}
q.ci=q + 
  geom_ribbon( aes(ymin=ymin, ymax=ymax), alpha=.3, fill = '#e41a1c') +
  geom_line(color = '#e41a1c') +
  facet_wrap(~channel, ncol=3) +
  guides(fill="none") +
  #word boxes
geom_rect(data=words, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), color="grey20", alpha=.25, inherit.aes=F) + 
		scale_x_continuous(limits=c(-200,1000), breaks=c(seq(-200, 1000, by=100)), labels=c(seq(-200, 1000, by=100)))

```

Now we can plot it in a new window, and save it:
```{r}
quartz(height=9, width=15)
q.ci
ggsave("difference.wave.N400.pdf")
```

Plot 4: Difference waves for each participant

We use the "q" plot from before, and add the N400.effect dataset:
```{r}
i.difference=q + 
  facet_wrap(~channel, ncol=3) +
  guides(fill="none", color="none", alpha="none")+
  geom_line(data=N400.effect, aes(x=time.ms, y=mV, group = subject, alpha=.2), color=gray(.5))
```

Now we can plot it in a new window, and save it:
```{r}
quartz(height=9, width=15)
i.difference
ggsave("difference.wave.individuals.N400.pdf")
```


