Atypical brain responses to auditory spatial cues in adults with autism spectrum disorder
Jon Brock
31/07/2017
========================================================
1 Background
This code plots and analyzes data from a study of auditory scene analysis in adults with autism. Participants listened to dichotic pitch stimuli in which inter-aural timing (ITD) and amplitude (IAD) differences inserted into broadband noise result in the illusory perception of a tone or pitch in a different spatial location to the noise.
Behavioural data include performance on a task in which participants perceptually matched ITD and IAD stimuli and a detection task in which they had to indicate whether or not they could hear a pitch sound separate from the noise.
During the detection task, we recorded brain responses using EEG. We use the openxlsx package to import the pre-processed EEG data from a folder of excel workbooks and ggplot2 to plot event-related potentials for groups and for individuals. Topographic plots are made using the erpR package (other R packages such as eegkit do not currently support 128 channel EEG).
We then calculate the area under the curve for two ERP components, the Object-Related Negativity (ORN) and the P400 and perform analysis of variance of behavioural and EEG data using the ezANOVA function from the ez package.
1.1 References
Giorgio Arcara and Anna Petrova (2014). erpR: Event-related potentials (ERP) analysis, graphics and utility functions. R package version 0.2.0. https://CRAN.R-project.org/package=erpR
Michael A. Lawrence (2016). ez: Easy Analysis and Visualization of Factorial Experiments. R package version 4.4-0. https://CRAN.R-project.org/package=ez
Veema Lodhia, Michael Hautus, Blake W. Johnson and Jon Brock (in press). Atypical brain responses to auditory spatial cues in adults with autism spectrum disorder. European Journal of Neuroscience.
Alexander Walker (2017). openxlsx: Read, Write and Edit XLSX Files. R package version 4.0.0. https://CRAN.R-project.org/package=openxlsx
H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2009.
2 Libraries
library(reshape2) # for melt and cast
library(stringr) # for str_splt_fixed
library(ez) # for ezANOVA
library(plyr) # for ddply
library(ggplot2) # for ggplot
library(knitr) # for kable
library(cowplot) # for plot_grid
library(openxlsx) # for importing data from excel
library(erpR) # for topoplot
library(dplyr) # for data wrangling, filtering, selecting
library(akima) # required by erpR for extrapolation between electrodes
library(magick) # image manipulation
library(PerformanceAnalytics) # for correlation grids
sessionInfo()## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12.6
##
## locale:
## [1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
##
## attached base packages:
## [1] tcltk stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] PerformanceAnalytics_1.4.3541 xts_0.9-7
## [3] zoo_1.8-0 magick_0.3
## [5] akima_0.6-2 dplyr_0.5.0
## [7] erpR_0.2.0 rpanel_1.1-3
## [9] openxlsx_4.0.0 cowplot_0.7.0
## [11] knitr_1.15.1 ggplot2_2.2.1
## [13] plyr_1.8.4 ez_4.4-0
## [15] stringr_1.1.0 reshape2_1.4.2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.8 nloptr_1.0.4 tools_3.3.2
## [4] digest_0.6.11 lme4_1.1-12 evaluate_0.10
## [7] tibble_1.2 nlme_3.1-128 gtable_0.2.0
## [10] lattice_0.20-34 mgcv_1.8-15 Matrix_1.2-7.1
## [13] DBI_0.6-1 yaml_2.1.14 parallel_3.3.2
## [16] SparseM_1.74 MatrixModels_0.4-1 rprojroot_1.1
## [19] grid_3.3.2 nnet_7.3-12 R6_2.2.0
## [22] rmarkdown_1.3 sp_1.2-4 minqa_1.2.4
## [25] car_2.1-4 magrittr_1.5 backports_1.0.4
## [28] scales_0.4.1 htmltools_0.3.5 MASS_7.3-45
## [31] splines_3.3.2 assertthat_0.1 pbkrtest_0.4-6
## [34] colorspace_1.3-2 quantreg_5.29 stringi_1.1.2
## [37] lazyeval_0.2.0 munsell_0.4.33 Behavioural data
3.1 Matching task
Import data for matching task. Parse the condition label using str_splt_fixed from the stringr library.
MatchingDataWide <- read.csv("data/BehaviouralData/matchingTask.csv", stringsAsFactors=TRUE)
MatchingData <- melt(MatchingDataWide, id=c("SubNum","Group"), variable.name="Condition", value.name="Score")
MatchingData$Pitch <- str_split_fixed(MatchingData$Condition, "_", 2)[,1]
MatchingData$Pitch <- mapvalues(MatchingData$Pitch, from = c("Control", "Tone"), to = c("No Pitch", "Pitch"))
MatchingData$Laterality <- str_split_fixed(MatchingData$Condition, "_", 2)[,2]
cols <- c("Pitch", "Laterality")
MatchingData[cols] <- lapply(MatchingData[cols], factor)
str(MatchingData)## 'data.frame': 120 obs. of 6 variables:
## $ SubNum : Factor w/ 30 levels "ASDP1","ASDP10",..: 16 23 24 25 26 27 28 29 30 17 ...
## $ Group : Factor w/ 2 levels "ASD","TD": 2 2 2 2 2 2 2 2 2 2 ...
## $ Condition : Factor w/ 4 levels "Control_Left",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Score : num -0.3 -0.5 -0.4 -0.9 -0.5 -0.5 -0.4 -0.3 -0.7 -0.4 ...
## $ Pitch : Factor w/ 2 levels "No Pitch","Pitch": 1 1 1 1 1 1 1 1 1 1 ...
## $ Laterality: Factor w/ 2 levels "Left","Right": 1 1 1 1 1 1 1 1 1 1 ...3.1.1 Boxplot
Plot horizontal boxplot using ggplot2.
MatchingBoxplot <- ggplot(MatchingData, aes(factor(Laterality), y=Score, fill=Group))+
facet_grid(Pitch~.)+
geom_boxplot(outlier.colour=NA)+
scale_fill_brewer(palette="Oranges") +
theme_bw() +
theme(text=element_text(size=16), plot.title = element_text(hjust = 0.5), axis.line.x = element_line(color="black", size = 2), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), strip.background = element_blank()) + scale_x_discrete("") +
geom_point(position=position_jitterdodge(dodge.width=0.75, jitter.width=0.3), alpha=0.2, size=3)+
scale_y_continuous("Lateralization of IAD stimulus", limits=c(-1.1,1.1), breaks=seq(-1.0, 1.0, 0.5))+
theme(axis.title.y = element_text(vjust=1.0))+
coord_flip()
MatchingBoxplot
3.1.2 ANOVA
ANOVA on ArcSin-transformed absolute values of responses using ezANOVA from the ez package.
MatchingData$AbsScore <- ifelse(MatchingData$Laterality=="Left", -1*MatchingData$Score, MatchingData$Score)
MatchingData$AsinScore <- asin(MatchingData$AbsScore)
kable(ezANOVA(MatchingData, dv=AsinScore, wid=SubNum, within=.(Pitch,Laterality), between=Group, type = 3), row.names=FALSE, caption="ANOVA: Matching Task Performance", digits=3)
|
3.2 EEG behavioural performance
Import data for performance during the EEG recording.
EEG_TaskData <- read.csv("data/BehaviouralData/EEG_Task.csv", stringsAsFactors=FALSE)
EEG_TaskData_Long <- melt(EEG_TaskData, id=c("SubNum","Group"), variable.name="Condition", value.name="Score")
EEG_TaskData_Long$Cue <- str_split_fixed(EEG_TaskData_Long$Condition, "_", 2)[,1]
EEG_Data_Acc <- ddply(EEG_TaskData_Long, .(SubNum, Group, Cue), summarize, Score=mean(Score))
cols <- c("SubNum", "Group", "Cue")
EEG_Data_Acc[cols] <- lapply(EEG_Data_Acc[cols], factor)
str(EEG_Data_Acc)## 'data.frame': 60 obs. of 4 variables:
## $ SubNum: Factor w/ 30 levels "ASD1","ASD10",..: 1 1 2 2 3 3 4 4 5 5 ...
## $ Group : Factor w/ 2 levels "ASD","TD": 1 1 1 1 1 1 1 1 1 1 ...
## $ Cue : Factor w/ 2 levels "IAD","ITD": 1 2 1 2 1 2 1 2 1 2 ...
## $ Score : num 95.9 91.8 85.2 72.1 77.5 ...3.2.1 Boxplot
Plot boxplot using ggplot2.
EEG_Data_Acc$Cue <- factor(EEG_Data_Acc$Cue, levels = c("ITD", "IAD"))
AccuracyBoxplot <- ggplot(EEG_Data_Acc, aes(factor(Cue), y=Score, fill=Group))+
geom_boxplot(outlier.colour=NA)+
scale_fill_brewer(palette="Oranges") +
theme_bw() +
theme(text=element_text(size=16), plot.title = element_text(hjust = 0.5), axis.line.x = element_line(color="black", size = 2), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), strip.background = element_blank()) + scale_x_discrete("")+
geom_point(position=position_jitterdodge(dodge.width=0.75, jitter.width=0.3), alpha=0.2, size=3)+
scale_y_continuous("Proportion correct", limits=c(0,100), breaks=seq(0, 100, 20))+
theme(axis.title.y = element_text(vjust=1.0))
AccuracyBoxplot
3.2.2 ANOVA
ANOVA on ArcSin-transformed proportion correct.
EEG_Data_Acc$AsinScore <- asin(EEG_Data_Acc$Score/100)
kable(ezANOVA(EEG_Data_Acc, dv=AsinScore, wid=SubNum, within=.(Cue), between=Group, type = 3), row.names=FALSE, caption="ANOVA: Task Performance During EEG Recording", digits=3)
|
3.3 Multi-panel figure
Combine the two boxplots using plot_grid from the cowplot package. This is figure 2 in the paper.
Plots <- plot_grid(MatchingBoxplot + theme(legend.position="none"), AccuracyBoxplot + theme(legend.position="none"), labels = c('A', 'B'), label_size = 14, rel_widths = c(1.5, 1), vjust = -1, scale = 0.85)
legend <- get_legend(MatchingBoxplot)
tiff(filename="output/figures/Figure2.png", width=1000, height=500)
plot_grid(Plots, legend, rel_widths = c(2, .3))
dev.off()## quartz_off_screen
## 3knitr::include_graphics("output/figures/Figure2.png")
4 EEG data
4.1 Data processing
4.1.1 Import data from Excel
Offline analysis of EEG data included epoching, artefact rejection, averaging, and filtering. The resultant ERPs were exported as Excel workbooks, one per participant. Each workbook contains four sheets, one for each condition. Each sheet contains the event-related potential for each EEG channel.
The script below reads each worksheet in each workbook in turn using functions from the openxlsx library. It parses the name of each worksheet to determine the Condition, Group, Subject Number and adds these as columns to the dataframe. It then adds each new worksheet and each new workbook to a dataframe, AllChannelsData, which contains the data for each participant, condition, and channel. This is exported as a .csv file.
DatafileNum <- 0
fnames = dir("data/ERPs", full.names=TRUE, all.files=FALSE) # Check no hidden files in folder because that messes things up!!!
for(j in 1:length(fnames)) {
DatafileNum <- DatafileNum + 1
CurrentWorkBook <- loadWorkbook(fnames[j], xlsxFile = NULL)
SheetNames <- sheets(CurrentWorkBook)
for(i in 1:length(SheetNames)) {
SheetName <- SheetNames[i]
CurrentSheet <- readWorkbook(CurrentWorkBook, sheet = SheetName)
SheetNameLength <- nchar(SheetName)
SheetID <- ifelse (substr(SheetName,SheetNameLength,SheetNameLength)==" ", substr(SheetName,1,SheetNameLength-1), SheetName)
if(substr(SheetID,nchar(SheetID)-1,nchar(SheetID))=="_P"){
CurrentSheet$Pitch <- "Pitch"
CurrentSheet$Cue <- substr(SheetID,nchar(SheetID)-4,nchar(SheetID)-2)
CurrentSheet$SubjectID <- substr(SheetID,1,nchar(SheetID)-6)
} else {
CurrentSheet$Pitch <- "NoPitch"
CurrentSheet$Cue <- substr(SheetID,nchar(SheetID)-5,nchar(SheetID)-3)
CurrentSheet$SubjectID <- substr(SheetID,1,nchar(SheetID)-7)
}
CurrentSheet$Group <- ifelse(substr(SheetID,1,3)=="ASD","ASD","TD") # Extract group from first 3 characters of worksheet ID
CurrentSheet$SheetID <- SheetID
CurrentSheet$TimeMS <- CurrentSheet[,1]
CurrentSheet$Time <- as.numeric(substr(CurrentSheet$TimeMS, 1, nchar(CurrentSheet$TimeMS)-2))
CurrentSheet[,1] <- NULL
CurrentSheet$TimeMS <- NULL
if(i==1){
Waveforms <- CurrentSheet
} else {
Waveforms <- rbind(Waveforms, CurrentSheet)
}
}
if(DatafileNum==1) {
AllChannelsData <- Waveforms
} else {
AllChannelsData <- rbind(AllChannelsData, Waveforms)
}
}
AllChannelsData$SheetID <- NULL
AllChannelsDataLong <- melt(AllChannelsData, id.vars = c("Time","Group","SubjectID","Cue", "Pitch"), value.name = "Amplitude", variable.name = "Electrode")
write.csv(AllChannelsData, "output/data/AllChannels.csv")4.1.2 Calculate grand average difference waveform
Average across subjects and conditions and then calculate the difference at each time point and channel between Pitch and NoPitch condition.
GrandAveDataLong <- ddply(AllChannelsDataLong, c("Time", "Electrode", "Pitch"), summarise, Amplitude = mean(Amplitude))
GrandAveData <- dcast(GrandAveDataLong, Time + Electrode ~ Pitch, value.var="Amplitude")
GrandAveData <- mutate(GrandAveData, Amplitude = Pitch - NoPitch) 4.1.3 Butterfly plot
The butterfly plot shows the grand mean difference waveform for each channel. Note the peaks at around 230ms (ORN) and 450ms (P400).
ggplot(data=GrandAveData, aes(x = Time, y = Amplitude, group = Electrode)) +
geom_freqpoly(stat = "identity", size = I(.3)) +
ylab(expression(Amplitude~'/'~mu~V)) +
scale_x_continuous("Time / ms",limits=c(-100,750),breaks=seq(-100, 700, 100)) +
theme_bw() +
theme(text=element_text(size=14), plot.title = element_text(hjust = 0.5), legend.position="none", axis.line.x = element_line(color="black", size = 2), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()) +
geom_hline(yintercept = 0, size = 0.2) +
theme(strip.background = element_blank())
4.1.4 GFP plot
Plot the global field power by calculating the standard deviation across channels at each time point.
StDev_Diff <- ddply(GrandAveData, c("Time"), summarise, StDev = sd(Amplitude))
ggplot(data=StDev_Diff, aes(x = Time, y = StDev)) +
geom_freqpoly(stat = "identity", size = I(.3)) +
scale_y_continuous("Standard Deviation") +
scale_x_continuous("Time") +
theme_bw() +
theme(text=element_text(size=14), plot.title = element_text(hjust = 0.5), legend.position="none", axis.line.x = element_line(color="black", size = 2), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()) +
geom_hline(yintercept = 0, size = 0.2) +
theme(strip.background = element_blank())
4.1.5 Extract peak times
Extract peak time for ORN and P400. These times are used for the initial topographic plots, which are used to select the electrodes used in further analyses.
ORN_Window <- filter(StDev_Diff, Time>100, Time<300)
ORN_PeakTime <- ORN_Window$Time[which.max(ORN_Window$StDev)]
P400_Window <- filter(StDev_Diff, Time>300, Time<500)
P400_PeakTime <- P400_Window$Time[which.max(P400_Window$StDev)]4.1.6 Topographic plots for grand mean ORN and P400
For topographic plots, we first need to read in sensor position data.
SensorPos <- read.csv("data/SensorPos/SensorPos.csv", row.names=NULL)
SensorPos$x <- SensorPos$x / max(SensorPos$x)
SensorPos$y <- SensorPos$y / max(SensorPos$y)
SensorPos$z <- SensorPos$z / max(SensorPos$z)We use the topoplot function in the erpR package. This requires the data in wide format with a column for each channel and a row for each timepoint. We also need to remove the Time column from the dataframe.
GrandAveWide <- select(dcast(GrandAveData, Time ~ Electrode, value.var = "Amplitude"), -Time)The topographic plot for the ORN shows a circular cluster of electrodes with a clear negativity.
topoplot(GrandAveWide, startmsec=-100, endmsec=748, win.ini=ORN_PeakTime, win.end=ORN_PeakTime, elec.coord = SensorPos,
projection="orthographic", zlim = c(-0.7, 0.7),
palette.steps = 100, palette.col = "jet", interpolation = "linear", draw.elec.lab = TRUE, contour=FALSE,
draw.nose=TRUE, draw.elec.pos = TRUE)
For the P400, the topographic plot shows a positivity in a similar cluster of electrodes.
topoplot(GrandAveWide, startmsec=-100, endmsec=748, win.ini=P400_PeakTime, win.end=P400_PeakTime, elec.coord = SensorPos,
projection="orthographic", zlim = c(-0.5, 0.5),
palette.steps = 100, palette.col = "jet", interpolation = "linear", draw.elec.lab = TRUE, contour=FALSE,
draw.nose=TRUE, draw.elec.pos = TRUE)
4.1.7 Average across selected channels
Based on these topographic plots, we average across a symmetrical cluster of channels showing robust ORN and P400 effects. Note that this biases the results towards findings a main effect of Pitch but should be orthogonal to any interactions involving Group or Cue.
AllChannelsData$Amplitude <- rowMeans(subset(AllChannelsData, select = c(e13, e6, e112, e30, e7, e106, e105, e37, e31, e129, e80, e87, e54, e55, e79)), na.rm = TRUE)
WaveformData <- AllChannelsData[,c("Time", "Amplitude", "Pitch", "Cue", "SubjectID", "Group")]
cols <- c("Pitch", "Cue", "SubjectID", "Group")
WaveformData[cols] <- lapply(WaveformData[cols], factor)
knitr::include_graphics("images/image003.png") 
Here we show the waveforms for each participant averaged across conditions. This shows that all participants have clear ERPs to the auditory stimuli.
IndWFData <- ddply(WaveformData, c("Group", "SubjectID", "Time"), summarise, Amplitude = mean(Amplitude))
GroupWFData <- ddply(WaveformData, c("Group", "Time"), summarise, Amplitude = mean(Amplitude))
GroupWFData$SubjectID <- "XXX"
GroupWFData$Colour <- "a"
IndWFData$Colour <- "b"
PlotData <- rbind(GroupWFData, IndWFData)
ggplot(data=PlotData, aes(x = Time, y = Amplitude, group = SubjectID, colour = Colour)) +
geom_freqpoly(stat = "identity", size = I(.3)) +
facet_grid(. ~ Group)+
scale_y_continuous(expression(Amplitude~'/'~mu~V), limits=c(-4.5,4.5), breaks=seq(-5, 5, 1)) +
scale_x_continuous("Time / ms",limits=c(-100,750),breaks=seq(-100, 700, 100)) +
theme_bw() +
theme(text=element_text(size=14), plot.title = element_text(hjust = 0.5), legend.position="none", axis.line.x = element_line(color="black", size = 2), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()) +
geom_hline(yintercept = 0, size = 0.2) +
theme(strip.background = element_blank()) +
scale_colour_manual(values=c("black", "grey"))
4.2 Determine ORN and P400 time windows
Next, we calculate the difference waveforms.
WaveformData_Wide <- dcast(WaveformData, SubjectID + Group + Cue + Time ~ Pitch, value.var="Amplitude")
WaveformData_Wide <- mutate(WaveformData_Wide, Difference=Pitch - NoPitch)Average across cues and subjects to create the Grand Mean Difference Waveform. Determine the full width half max windows for the ORN and P400 based on the Grand Mean Difference Waveform.
GMDifferenceData <- ddply(WaveformData_Wide, .(Time), summarise, Difference = mean(Difference))
ORN_Peak <- min(GMDifferenceData$Difference)
GM_ORN_Window <- GMDifferenceData[(GMDifferenceData$Difference < ORN_Peak / 2),]
ORN_Start <- min(GM_ORN_Window$Time)
ORN_End <- max(GM_ORN_Window$Time)
P400_Peak <- max(GMDifferenceData$Difference)
GM_P400_Window <- GMDifferenceData[(GMDifferenceData$Difference > P400_Peak / 2),]
P400_Start <- min(GM_P400_Window$Time)
P400_End <- max(GM_P400_Window$Time)
Windows.DF <- data.frame(Component = c("ORN", "P400"),
StartTime = c(ORN_Start, P400_Start),
EndTime = c(ORN_End, P400_End))
kable(Windows.DF, row.names=FALSE, caption="Start and end times for ORN and P400", digits=0)| Component | StartTime | EndTime |
|---|---|---|
| ORN | 196 | 268 |
| P400 | 360 | 564 |
Plot group averaged waveforms for each condition. Superimpose the ORN and P400 time windows (geom_rect). Export as PNG. This is figure 3 in the paper.
WaveformPlotData <- ddply(WaveformData, c("Time","Pitch","Cue","Group"), summarise, Amplitude = mean(Amplitude))
WaveformPlotData$Pitch <- mapvalues(WaveformPlotData$Pitch, from = "NoPitch", to = "No Pitch")
WaveformPlotData$Cue <- factor(WaveformPlotData$Cue, levels = c("ITD", "IAD"))
Fig3 <- ggplot(data=WaveformPlotData, aes(x = Time, y = Amplitude, colour = Pitch)) +
geom_rect(aes(xmin=ORN_Start, xmax=ORN_End, ymin=-Inf, ymax=Inf), colour="gray97", fill="gray97") +
geom_rect(aes(xmin=P400_Start, xmax=P400_End, ymin=-Inf, ymax=Inf), colour="gray97", fill="gray97") +
geom_hline(yintercept = 0, size = 0.2) +
geom_freqpoly(stat = "identity", size = I(1.0)) +
facet_grid(Cue ~ Group)+
scale_y_continuous(expression(Amplitude~'/'~mu~V), limits=c(-2.5,2.5), breaks=seq(-5, 5, 1)) +
scale_x_continuous("Time / ms",limits=c(-100,750),breaks=seq(-100, 700, 100)) +
theme_bw() +
theme(text=element_text(size=20), plot.title = element_text(hjust = 0.5), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), strip.background = element_blank(), legend.title=element_blank(), legend.position="none") +
scale_colour_manual(values=c("dark grey", "dark orange", "grey"))
tiff(filename="output/figures/Figure3.png", width=800, height=800)
Fig3
dev.off()## quartz_off_screen
## 2knitr::include_graphics("output/figures/Figure3.png")
4.3 Topographic plots
Plot topographies for the ORN and P400 by group and cue. First, we average the data within each group.
GroupData <- ddply(AllChannelsDataLong, c("Time","Pitch","Cue","Group", "Electrode"), summarise, Amplitude = mean(Amplitude))Then calculate the difference between Pitch and NoPitch conditions
GroupDataWide <- dcast(GroupData, Group + Cue + Electrode + Time ~ Pitch, value.var="Amplitude", fun.aggregate=mean)
GroupDataWide <- mutate(GroupDataWide, Difference = Pitch - NoPitch)
GroupDataWide <- select(GroupDataWide, -Pitch, -NoPitch)Then arrange the data so that it can be read by the topoplot function (below). For each group and cue, there is one column for each channel and no other columns.
AllDiff <- dcast(GroupDataWide, Group + Cue + Time ~ Electrode, value.var = "Difference")
ASD_ITD <- select(filter(AllDiff, Group=="ASD", Cue=="ITD"), -Group, -Cue, -Time)
ASD_IAD <- select(filter(AllDiff, Group=="ASD", Cue=="IAD"), -Group, -Cue, -Time)
TD_ITD <- select(filter(AllDiff, Group=="TD", Cue=="ITD"), -Group, -Cue, -Time)
TD_IAD <- select(filter(AllDiff, Group=="TD", Cue=="IAD"), -Group, -Cue, -Time)The DrawTopo function calls the the topoplot function from the erpR package.
DrawTopo <- function (Data, WindowStart, WindowEnd, FigName="Topology", width=500, height=500, Title="", DrawLabels=FALSE) {
ImageFileName <- paste0("output/figures/", FigName, ".png")
png(filename=ImageFileName, width=width, height=height)
topoplot(Data, startmsec=-100, endmsec=748,
win.ini=WindowStart, win.end=WindowEnd,
elec.coord = SensorPos,
projection="orthographic", interpolation = "linear", extrap = FALSE, interp.points = 200,
palette.steps = 50, palette.col = "jet", contour=FALSE,
draw.elec.lab = DrawLabels, draw.elec.pos = TRUE, draw.nose=TRUE,
zlim = c(-0.7, 0.7),
mask=TRUE)
dev.off()
image_read(ImageFileName)
}Plot ORN topology
ASD_ITD_ORN <- DrawTopo (ASD_ITD, ORN_Start, ORN_End, "ASD_ITD_ORN")
ASD_ITD_ORN <- image_crop(ASD_ITD_ORN, "360x390+80+0")
ASD_ITD_ORN <- image_annotate(ASD_ITD_ORN, "ASD", size = 40, gravity = "northwest", location = "+220+20", color = "black")
TD_ITD_ORN <- DrawTopo (TD_ITD, ORN_Start, ORN_End, "TD_ITD_ORN")
TD_ITD_ORN <- image_crop(TD_ITD_ORN, "420x390+80+0")
TD_ITD_ORN <- image_annotate(TD_ITD_ORN, "TD", size = 40, gravity = "northwest", location = "+240+20", color = "black")
ASD_IAD_ORN <- DrawTopo (ASD_IAD, ORN_Start, ORN_End, "ASD_IAD_ORN")
ASD_IAD_ORN <- image_crop(ASD_IAD_ORN, "360x360+80+40")
TD_IAD_ORN <- DrawTopo (TD_IAD, ORN_Start, ORN_End, "TD_IAD_ORN")
TD_IAD_ORN <- image_crop(TD_IAD_ORN, "420x360+80+40")
ITD_ORN <- image_append (c(ASD_ITD_ORN, TD_ITD_ORN))
IAD_ORN <- image_append (c(ASD_IAD_ORN, TD_IAD_ORN))
ORN_Topo <- image_append (c(ITD_ORN, IAD_ORN), stack = TRUE)Plot P400 topology
ASD_ITD_P400 <- DrawTopo (ASD_ITD, P400_Start, P400_End, "ASD_ITD_P400")
ASD_ITD_P400 <- image_crop(ASD_ITD_P400, "360x390+80+0")
ASD_ITD_P400 <- image_annotate(ASD_ITD_P400, "ASD", size = 40, gravity = "northwest", location = "+220+20", color = "black")
TD_ITD_P400 <- DrawTopo (TD_ITD, P400_Start, P400_End, "TD_ITD_P400")
TD_ITD_P400 <- image_crop(TD_ITD_P400, "420x390+80+0")
TD_ITD_P400 <- image_annotate(TD_ITD_P400, "TD", size = 40, gravity = "northwest", location = "+240+20", color = "black")
TD_ITD_P400 <- image_annotate(TD_ITD_P400, "ITD", size = 40, gravity = "northwest", location = "+490+210", color = "black", degrees = 90)
ASD_IAD_P400 <- DrawTopo (ASD_IAD, P400_Start, P400_End, "ASD_IAD_P400")
ASD_IAD_P400 <- image_crop(ASD_IAD_P400, "360x360+80+40")
TD_IAD_P400 <- DrawTopo (TD_IAD, P400_Start, P400_End, "TD_IAD_P400")
TD_IAD_P400 <- image_crop(TD_IAD_P400, "420x360+80+40")
TD_IAD_P400 <- image_annotate(TD_IAD_P400, "IAD", size = 40, gravity = "northwest", location = "+490+220", color = "black", degrees = 90)
ITD_P400 <- image_append (c(ASD_ITD_P400, TD_ITD_P400))
IAD_P400 <- image_append (c(ASD_IAD_P400, TD_IAD_P400))
P400_Topo <- image_append (c(ITD_P400, IAD_P400), stack = TRUE)Combine the ORN and P400 plots.
Topographies <- image_append (c(ORN_Topo, P400_Topo))Add a legend (palette). To ensure the colours in the legend match the figure, we have to create a dummy plot with the same properties
PalettePlot <- topoplot(ASD_ITD, startmsec=-100, endmsec=748,
win.ini=0, win.end=100,
elec.coord = SensorPos,
projection="orthographic", interpolation = "linear", extrap = FALSE, interp.points = 200,
palette.steps = 50, palette.col = "jet", contour=FALSE,
draw.elec.lab = FALSE, draw.elec.pos = FALSE, draw.nose=TRUE,
zlim = c(-0.7, 0.7),
mask=TRUE)
tiff(filename="output/figures/TopoPal.png", width=image_info(Topographies)$width, height=image_info(Topographies)$height/4)
plot.new()
topoplot.palette(cols=PalettePlot$palette, pos = c(0.5,0.5), p.width=0.7, p.height=0.3, horizontal = TRUE, rev= FALSE, palette.lwd=1, palette.bord="black", palette.lim=PalettePlot$zlim, labels=TRUE, lab.dist = 1, lab.cex = 1, lab.font = 1, lab.family = "")
dev.off()## quartz_off_screen
## 2Palette <- image_read("output/figures/TopoPal.png")
Topographies <- image_append (c(Topographies, Palette), stack = TRUE)Export the topography image
image_write(Topographies, "output/figures/Fig5.png")
knitr::include_graphics("output/figures/Fig5.png")
4.4 Difference waveforms
Plot the difference waveforms, including individual subjects and group averages. Export as PNG.
IndDiffWFData <- subset(WaveformData_Wide, select = -c(NoPitch, Pitch))
GroupDiffWFData <- ddply(WaveformData_Wide, c("Group", "Time", "Cue"), summarise, Difference = mean(Difference))
GroupDiffWFData$SubjectID <- "XXX"
GroupDiffWFData$Colour <- "a"
IndDiffWFData$Colour <- "b"
PlotData <- rbind(GroupDiffWFData, IndDiffWFData)
PlotData$Cue <- factor(PlotData$Cue, levels = c("ITD", "IAD"))
DiffWFs <- ggplot(data=PlotData, aes(x = Time, y = Difference, group = SubjectID, colour = Colour)) +
geom_rect(aes(xmin=ORN_Start, xmax=ORN_End, ymin=-Inf, ymax=Inf), colour="gray97", fill="gray97") +
geom_rect(aes(xmin=P400_Start, xmax=P400_End, ymin=-Inf, ymax=Inf), colour="gray97", fill="gray97") +
geom_hline(yintercept = 0, size = 0.2) +
geom_freqpoly(stat = "identity", size = I(.3)) +
facet_grid(Cue ~ Group)+
scale_y_continuous(expression(Amplitude~'/'~mu~V),limits=c(-2.5,2.5), breaks=seq(-5, 5, 0.5)) +
scale_x_continuous("Time / ms",limits=c(-100,750),breaks=seq(-100, 700, 100)) +
theme_bw() +
theme(text=element_text(size=20), plot.title = element_text(hjust = 0.5), legend.position="none", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), strip.background = element_blank()) +
scale_colour_manual(values=c("black", "orange"))
tiff(filename="output/figures/Fig4.png", width=800, height=800)
DiffWFs
dev.off()## quartz_off_screen
## 2knitr::include_graphics("output/figures/Fig4.png")
4.5 Area under the curve analyses
4.5.1 Data processing
Calculate the under the curve for the ORN and P400. Combine into one dataframe.
ORN_Data <- ddply(WaveformData_Wide[(WaveformData_Wide$Time>=ORN_Start & WaveformData_Wide$Time<=ORN_End),], c("Cue","Group","SubjectID"), summarise, Amplitude = mean(Difference))
P400_Data <- ddply(WaveformData_Wide[(WaveformData_Wide$Time>=P400_Start & WaveformData_Wide$Time<=P400_End),], c("Cue","Group","SubjectID"), summarise, Amplitude = mean(Difference))
ORN_Data$Component <- factor("ORN")
P400_Data$Component <- factor("P400")
AUC.DF <- rbind(ORN_Data, P400_Data)
str(AUC.DF)## 'data.frame': 120 obs. of 5 variables:
## $ Cue : Factor w/ 2 levels "IAD","ITD": 1 1 1 1 1 1 1 1 1 1 ...
## $ Group : Factor w/ 2 levels "ASD","TD": 1 1 1 1 1 1 1 1 1 1 ...
## $ SubjectID: Factor w/ 30 levels "ASD_P1","ASD_P10",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Amplitude: num -0.46 -0.207 -0.543 -0.452 -0.712 ...
## $ Component: Factor w/ 2 levels "ORN","P400": 1 1 1 1 1 1 1 1 1 1 ...4.5.2 T-tests
Run one-sample t-tests on each component and print as table.
TTest.DF <- data.frame(Group = c("TD", "ASD", "TD", "ASD", "TD", "ASD", "TD", "ASD"),
Component = c("ORN", "ORN", "ORN", "ORN", "P400", "P400", "P400", "P400"),
Cue = c("IAD", "IAD", "ITD", "ITD", "IAD", "IAD", "ITD", "ITD"),
df = numeric(length(8)),
t = numeric(length(8)),
p = numeric(length(8)))
for(i in 1:8) {
TTestData <- AUC.DF[((AUC.DF$Group==TTest.DF$Group[i]) & (AUC.DF$Component==TTest.DF$Component[i]) & (AUC.DF$Cue==TTest.DF$Cue[i]) ),]
TTest <- t.test (TTestData$Amplitude, mu=0)
TTest.DF$df[i] <- TTest$parameter
TTest.DF$t[i] <- TTest$stat
TTest.DF$p[i] <- TTest$p.value
}
kable(TTest.DF, row.names=FALSE, caption="One-sample t-test for each component", digits=3)| Group | Component | Cue | df | t | p |
|---|---|---|---|---|---|
| TD | ORN | IAD | 14 | -3.947 | 0.001 |
| ASD | ORN | IAD | 14 | -3.740 | 0.002 |
| TD | ORN | ITD | 14 | -3.157 | 0.007 |
| ASD | ORN | ITD | 14 | -0.743 | 0.470 |
| TD | P400 | IAD | 14 | 3.767 | 0.002 |
| ASD | P400 | IAD | 14 | -0.607 | 0.554 |
| TD | P400 | ITD | 14 | 5.796 | 0.000 |
| ASD | P400 | ITD | 14 | 1.891 | 0.080 |
4.5.3 Boxplots
Generate boxplot in ggplot. Save the legend for later.
AUC.DF$Cue <- factor(AUC.DF$Cue, levels = c("ITD", "IAD"))
AUC_Boxplot <- ggplot(AUC.DF, aes(factor(Cue), y=Amplitude, fill=Group))+
geom_hline(yintercept = 0, size = 0.2) +
facet_grid(.~Component)+
geom_boxplot(outlier.colour=NA)+
scale_fill_brewer(palette="Oranges") +
scale_x_discrete("")+
geom_point(position=position_jitterdodge(dodge.width=0.75, jitter.width=0.3), alpha=0.2, size=3)+
scale_y_continuous(expression(Amplitude~'/'~mu~V),limits=c(-1.5,1.5), breaks=seq(-1.6, 1.6, 0.2)) +
theme_bw() +
theme(text=element_text(size=16), plot.title = element_text(hjust = 0.5), axis.line.x = element_line(color="black", size = 2), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), strip.background = element_blank()) +
theme(axis.title.y = element_text(vjust=1.0))
legend <- get_legend(AUC_Boxplot)Add asterisks to boxplot to indicate significant effects.
TTest.DF$Sig <- as.factor(ifelse(TTest.DF$p<0.05,1,0))
TTest.DF$Amplitude = rep(1.3, 8)
Asterisks <- filter(TTest.DF, Sig == 1)
AUC_Boxplot <- AUC_Boxplot +
geom_point(data=TTest.DF, aes(factor(Cue), colour=Sig, fill=Group), shape=42, size = 8, position=position_dodge(width=0.8)) +
scale_color_manual(values=c("white", "black"))Export figure (using original legend).
Fig6 <- plot_grid(AUC_Boxplot + theme(legend.position="none"), legend, rel_widths = c(2, .5))
png(filename="output/figures/Fig6.png", width=600, height=500)
Fig6
dev.off()## quartz_off_screen
## 2knitr::include_graphics("output/figures/Fig6.png")
4.6 Analysis of variance
4.6.2 P400
Two-way ANOVA between Group and Cue on the P400 data.
kable(ezANOVA(P400_Data, dv = Amplitude, wid = SubjectID, within = .(Cue), between = .(Group), type = 3), row.names=FALSE, caption="ANOVA: P400", digits=3)
|
4.7 Correlations
AUC_Wide <- dcast(AUC.DF, Group + SubjectID ~ Component + Cue, value.var="Amplitude")
AUC_Wide$SubNum <- str_replace(AUC_Wide$SubjectID, "_P", "")
Acc_Wide <- dcast(EEG_Data_Acc, SubNum ~ Cue, value.var="AsinScore")
AUC_Wide$ScoreIAD <- Acc_Wide$IAD[match(AUC_Wide$SubNum, Acc_Wide$SubNum)]
AUC_Wide$ScoreITD <- Acc_Wide$ITD[match(AUC_Wide$SubNum, Acc_Wide$SubNum)]
Demographics <- read.csv("data/Demographics/Demographics.csv", stringsAsFactors=FALSE)
Demographics$SubjectID <- paste0(Demographics$Group, "_", Demographics$SubNum)
AUC_Wide$SCQ <- Demographics$SCQ_Lifetime[match(AUC_Wide$SubjectID, Demographics$SubjectID)]
CorrDF_ASD <- select(filter(AUC_Wide, Group=="ASD"), -SubNum, -SubjectID, -Group)
chart.Correlation(CorrDF_ASD, histogram=TRUE, pch=19)
CorrDF_TD <- select(filter(AUC_Wide, Group=="TD"), -SubNum, -SubjectID, -Group, -SCQ)
chart.Correlation(CorrDF_TD, histogram=TRUE, pch=19)