--- title: "Plotting a isovoltage topoplots. This is adapted from some code 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(akima) #for interpolation library(scales) library(mgcv) #another interpolation option library(gridExtra) library(png) library(grid) ``` Load data from individual ERPs. One could load directly from a grandaveraged ERP, but since we used individual ERPs for the waveforms, we can use them again here. * 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 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") ``` This assumes that we have three bins: the two conditions, and a bin that is the difference between them. If you don't have a difference bin, you need to calcluate a difference. There is code a bit further down to do that. But we can skip it here. ```{r} all = filter(full.locs, condition == "congruent" | condition == "incongruent" | condition == "N400") all = droplevels(all) ``` Choose a subset of time that includes the words we want ```{r} plottable = all %>% filter(time.ms>= -200 & time.ms<=1000) ``` Do you need to re-baseline? If so, here is some code to do it right here in R. ```{r} #calculate the baseline value # baseline = plottable %>% # filter(time.ms >= -200 & time.ms<=0) %>% # group_by(subject, channel, condition) %>% # summarise(microvolts = mean(mV)) %>% # arrange(subject,channel, condition) %>% # rename(baseline = microvolts) # # #left_join the baseline value to the original dataset # plottable = left_join(plottable, baseline) # # #create the re-baselined dataset # plottable = plottable %>% # mutate(mV.original = mV, mV = mV.original - baseline) ``` Choose a subset to work on. In this case, we only have one subset, so this seems like a waste. But if we had more than one subset, we could choose then one we want here. ```{r} plotting = plottable ``` Here is some code to calculate difference effects. If you already have this in your data set, you don't need this. ```{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), condition=rep("difference", times=length(N400.difference)), mV = N400.difference) # # #now bind it to the original data set # plotting2 = rbind(plotting, select(N400.effect, channel, subject, condition, time.ms, mV, x, y)) ``` Calculate mean values in the windows that we want to plot. This also labels the windows. ```{r} startList = c(100, 200, 300) endList = c(150, 250, 500) for(i in 1:length(startList)){ #select a time temp = plotting %>% filter(time.ms >= startList[i] & time.ms<=endList[i]) %>% group_by(condition, channel, y, x) %>% summarise(microvolts = mean(mV)) if(startList[i] == 100 | startList[i] == 600 | startList[i] == 1100){ temp2 = data.frame(window = factor(rep("N1", times=nrow(temp))), temp)} if(startList[i] == 200 | startList[i] == 700 | startList[i] == 1200){ temp2 = data.frame(window = factor(rep("P2", times=nrow(temp))), temp)} if(startList[i] == 300 | startList[i] == 800 | startList[i] == 1300){ temp2 = data.frame(window = factor(rep("300-500", times=nrow(temp))), temp)} if(i==1){ plottable=temp2 }else{ plottable = rbind(plottable, temp2) } } ``` ##Draw the head, and choose some colors Do we want to define our own colors? ```{r} #this is matlab's jet palette jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")) #this is the 9 color diverging palette from colorbrewer #it comes out a little pale jon1=colorRampPalette(c("#2166ac", "#4393c3", "#92c5de", "#d1e5f0", "#f7f7f7", "#fddbc7", "#f4a582", "#d6604d", "#b2182b")) #This is 13 color sequence that I made up by taking colors from two sequential palettes jon2=colorRampPalette(c("#08306b","#08519c","#2171b5","#4292c6","#6baed6","#9ecae1", "#f7f7f7","#fc9272", "#fb6a4a", "#ef3b2c", "#cb181d","#a50f15","#67000d")) #This is 7 of the 13 colors above jon3=colorRampPalette(c("#08306b","#2171b5","#6baed6", "#f7f7f7", "#fb6a4a", "#cb181d","#67000d")) #We can also use colorbrewer palettes directly ``` Draw the head (taken directly from Matt's code): ```{r} theme_topo <- function(base_size = 12) { theme_bw(base_size = base_size) %+replace% theme( rect = element_blank(), line = element_blank(), axis.text = element_blank(), axis.title = element_blank() ) } circleFun <- function(center = c(0,0),diameter = 1, npoints = 100) { r = diameter / 2 tt <- seq(0,2*pi,length.out = npoints) xx <- center[1] + r * cos(tt) yy <- center[2] + r * sin(tt) return(data.frame(x = xx, y = yy)) } headShape <- circleFun(c(0, 0), round(max(coords$x)), npoints = 100) # 0 nose <- data.frame(x = c(-0.075,0,.075),y=c(.495,.575,.495)) ggplot(headShape,aes(x,y))+ geom_path()+ geom_text(data = coords, aes(x, y, label = channel))+ geom_line(data = nose, aes(x, y))+ theme_topo()+ coord_equal() ``` ##Interpolation The trick here is that we want to interpolate based on each grouping factor (effect and time). I can't get dplyr pipelines to work here. I don't know why. So I am just going to loop through the dataframe and re-assemble it ```{r} gridRes <- 150 # Specify the number of points for each grid dimension i.e. the resolution/smoothness of the interpolation #loop through effects, and through times #first, define a results dataframe. We are going to rbind() into this. results = data.frame() for(i in 1:length(levels(plottable$condition))){ for(j in 1:length(levels(plottable$window))){ tempData = filter(plottable, condition == levels(plottable$condition)[i] & window == levels(plottable$window)[j]) #need to check if there is data here, because the final word in the sentence doesn't have a 500-800 window if(nrow(tempData)>1){ tmpTopo <- with(tempData, interp(x = x, y = y, z = microvolts, xo = seq(min(x)*2, max(x)*2, length = gridRes), yo = seq(min(y)*2, max(y)*2, length = gridRes), linear = FALSE, extrap = TRUE) ) interpTopo <- data.frame(x = tmpTopo$x, tmpTopo$z) names(interpTopo)[1:length(tmpTopo$y)+1] <- tmpTopo$y interpTopo <- gather(interpTopo, key = y, value = microvolts, -x, convert = TRUE) interpTopo$incircle <- sqrt(interpTopo$x^2 + interpTopo$y^2) < .7 # mark grid elements that are outside of the plotting circle interpTopo <- interpTopo[interpTopo$incircle,] #remove the elements outside the circle tempInterpTopo = data.frame(condition = rep(as.character(levels(plottable$condition)[i]), times=nrow(interpTopo)), window=rep(as.character(levels(plottable$window)[j]), times=nrow(interpTopo)), interpTopo) #add to results file results = rbind(results, tempInterpTopo, row.names=NULL) } else {next} #move to next value of window loop if there is no data to work with } } ``` Here we change some things for aesthetic reasons. ```{r} #reorder the levels of the effect factor results$condition = factor(results$condition, levels=c("congruent", "incongruent", "N400")) results$window = factor(results$window, levels=c("N1", "P2", "300-500")) #define a list of points for the channel locations in the plot dots = filter(plottable, channel !="FT9" & channel !="FT10" & channel !="M1" & channel !="M2") ``` ```{r} maskRing <- circleFun(diameter = 1.42) #create a circle round the outside of the plotting area to mask the jagged edges of the interpolation akimaPlot <- ggplot(results, aes(x = x, y = y, fill = microvolts) ) + facet_grid(condition ~ window) + geom_raster() + stat_contour(aes(z = microvolts), colour = "black", binwidth = 0.5, alpha=0.2) + #this sets the opacity/transparency of the lines in the plot #theme_topo()+ theme_minimal()+ theme(line=element_blank(), axis.text = element_blank(),axis.title = element_blank())+ scale_fill_gradientn(colours = jon1(10), limits = c(-2,2), guide = "colourbar", oob = squish) + #theme(legend.position = "bottom") + #guides(fill = guide_colorbar(direction = "horizontal")) + geom_path(data = maskRing, aes(x, y, z = NULL, fill =NULL), colour = "white", size = 6)+ geom_point(data = dots, aes(x, y), size = 1)+ geom_path(data = headShape, aes(x, y, z = NULL, fill = NULL), size = 1.5)+ geom_path(data = nose, aes(x, y, z = NULL, fill = NULL), size = 1.5)+ coord_equal() akimaPlot ``` plot it and save it ```{r} quartz(width = 9, height = 9) akimaPlot ggsave("isovoltage.N400.pdf") ```