ERP package : erp_group_ggplot
data
preprocessing
EDA
test plot
group_bootstrap
library(ERP);require(mnormt);require(fdrtool);library(reshape);library(ggplot2);library(dplyr)
data(simerp);data(erpcz)Data Generate for Example
set.seed(123)
erp <- erpcz[,1:251]
m <- apply(erp,2,mean)
s <- apply(erp,2,sd)
for (j in 41:300){
for (i in 1:251){
erp[j:300,i] <- rnorm(1,m[i],s[i])
}
}
erp$Subject <- sort(rep(paste("S",1:10,sep=""),30))
erp$Electrode <- rep(c(rep("cz",10),rep("fz",10),rep("pz",10)),10)
erp$Condition <- rep(c(rep("Exp",5),rep("Con",5)),30)
erp$Gender <- c(rep("Male",30),rep("Female",30))
erp$Scores <- sort(rep(round(rnorm(10,70,5),2),30))
erp$RT <- abs(round(rnorm(300,1,5),2))
head(erp,100)[,251:257]## X1000 Subject Electrode Condition Gender Scores RT
## 1 -1.533214808 S1 cz Exp Male 65.23 3.05
## 2 3.632015705 S1 cz Exp Male 65.23 4.08
## 3 1.431424856 S1 cz Exp Male 65.23 3.17
## 4 4.562942982 S1 cz Exp Male 65.23 1.53
## 5 -1.563839793 S1 cz Exp Male 65.23 9.93
## 6 -1.327963471 S1 cz Con Male 65.23 0.88
## 7 1.621153355 S1 cz Con Male 65.23 0.43
## 8 1.150131226 S1 cz Con Male 65.23 4.60
## 9 2.264639139 S1 cz Con Male 65.23 5.15
## 10 -1.985386491 S1 cz Con Male 65.23 6.16
## 11 5.289675713 S1 fz Exp Male 65.23 5.91
## 12 -2.535944462 S1 fz Exp Male 65.23 6.44
## 13 -6.433563232 S1 fz Exp Male 65.23 3.53
## 14 -0.977460086 S1 fz Exp Male 65.23 4.69
## 15 -2.001248837 S1 fz Exp Male 65.23 2.78
## 16 1.214420080 S1 fz Con Male 65.23 4.61
## 17 -1.661339641 S1 fz Con Male 65.23 6.76
## 18 -4.182179928 S1 fz Con Male 65.23 0.03
## 19 3.180048943 S1 fz Con Male 65.23 2.43
## 20 -7.429838181 S1 fz Con Male 65.23 3.36
## 21 -0.004267288 S1 pz Exp Male 65.23 8.03
## 22 -1.869916916 S1 pz Exp Male 65.23 3.24
## 23 -1.038703322 S1 pz Exp Male 65.23 1.83
## 24 -3.446271420 S1 pz Exp Male 65.23 2.08
## 25 -2.430668592 S1 pz Exp Male 65.23 5.55
## 26 -4.899024963 S1 pz Con Male 65.23 2.06
## 27 -1.065935373 S1 pz Con Male 65.23 4.68
## 28 -7.618601799 S1 pz Con Male 65.23 5.96
## 29 0.392958701 S1 pz Con Male 65.23 0.50
## 30 -6.202139378 S1 pz Con Male 65.23 0.12
## 31 -1.812034965 S10 cz Exp Female 65.63 9.68
## 32 -0.265757233 S10 cz Exp Female 65.63 1.81
## 33 -3.333134413 S10 cz Exp Female 65.63 9.70
## 34 -7.367009640 S10 cz Exp Female 65.63 2.26
## 35 -0.580906868 S10 cz Exp Female 65.63 4.72
## 36 -4.479408264 S10 cz Con Female 65.63 2.67
## 37 -1.075744629 S10 cz Con Female 65.63 8.36
## 38 -20.677419660 S10 cz Con Female 65.63 1.10
## 39 -1.764858484 S10 cz Con Female 65.63 1.26
## 40 6.576538086 S10 cz Con Female 65.63 2.77
## 41 -3.466308342 S10 fz Exp Female 65.63 0.63
## 42 -6.280486029 S10 fz Exp Female 65.63 0.99
## 43 0.572520481 S10 fz Exp Female 65.63 2.45
## 44 -2.357986621 S10 fz Exp Female 65.63 9.47
## 45 4.123173844 S10 fz Exp Female 65.63 3.96
## 46 7.928965749 S10 fz Con Female 65.63 9.39
## 47 0.221929806 S10 fz Con Female 65.63 0.40
## 48 -4.686863284 S10 fz Con Female 65.63 3.52
## 49 -6.919152230 S10 fz Con Female 65.63 4.40
## 50 5.507145856 S10 fz Con Female 65.63 1.46
## 51 -4.693926960 S10 pz Exp Female 65.63 3.92
## 52 -4.049228904 S10 pz Exp Female 65.63 1.00
## 53 -0.554156951 S10 pz Exp Female 65.63 0.95
## 54 -1.501209337 S10 pz Exp Female 65.63 1.92
## 55 -7.886915485 S10 pz Exp Female 65.63 3.22
## 56 1.416506659 S10 pz Con Female 65.63 9.27
## 57 0.943498066 S10 pz Con Female 65.63 4.79
## 58 -2.238611795 S10 pz Con Female 65.63 3.84
## 59 3.210088883 S10 pz Con Female 65.63 6.68
## 60 -8.112075518 S10 pz Con Female 65.63 6.97
## 61 -4.709261640 S2 cz Exp Male 66.15 2.68
## 62 1.958666009 S2 cz Exp Male 66.15 2.84
## 63 -7.828519498 S2 cz Exp Male 66.15 2.27
## 64 0.572771526 S2 cz Exp Male 66.15 3.66
## 65 -7.922105245 S2 cz Exp Male 66.15 1.68
## 66 3.328339040 S2 cz Con Male 66.15 8.08
## 67 -2.105812141 S2 cz Con Male 66.15 1.98
## 68 2.156455510 S2 cz Con Male 66.15 4.01
## 69 -7.385025706 S2 cz Con Male 66.15 3.80
## 70 1.808731984 S2 cz Con Male 66.15 3.73
## 71 -1.703246031 S2 fz Exp Male 66.15 8.14
## 72 -2.402564814 S2 fz Exp Male 66.15 3.25
## 73 -3.638497377 S2 fz Exp Male 66.15 0.38
## 74 -17.543610435 S2 fz Exp Male 66.15 4.61
## 75 -1.923190766 S2 fz Exp Male 66.15 7.51
## 76 2.536401658 S2 fz Con Male 66.15 0.51
## 77 2.773909693 S2 fz Con Male 66.15 7.00
## 78 0.124280198 S2 fz Con Male 66.15 7.46
## 79 2.028445710 S2 fz Con Male 66.15 0.67
## 80 0.358292427 S2 fz Con Male 66.15 0.16
## 81 -2.694107124 S2 pz Exp Male 66.15 0.72
## 82 3.434254922 S2 pz Exp Male 66.15 1.16
## 83 -4.743837746 S2 pz Exp Male 66.15 1.10
## 84 -5.626897686 S2 pz Exp Male 66.15 6.78
## 85 2.825942466 S2 pz Exp Male 66.15 0.29
## 86 -6.101178421 S2 pz Con Male 66.15 2.23
## 87 0.152494237 S2 pz Con Male 66.15 1.47
## 88 -5.240117610 S2 pz Con Male 66.15 3.95
## 89 -4.880834905 S2 pz Con Male 66.15 7.54
## 90 0.622705433 S2 pz Con Male 66.15 12.78
## 91 3.955105303 S3 cz Exp Female 67.01 0.81
## 92 -0.040470940 S3 cz Exp Female 67.01 2.85
## 93 -2.173208264 S3 cz Exp Female 67.01 1.16
## 94 -11.813141223 S3 cz Exp Female 67.01 1.75
## 95 -3.584641609 S3 cz Exp Female 67.01 7.20
## 96 4.115430314 S3 cz Con Female 67.01 2.17
## 97 -7.233071381 S3 cz Con Female 67.01 1.89
## 98 1.035752736 S3 cz Con Female 67.01 6.56
## 99 -2.600796293 S3 cz Con Female 67.01 9.55
## 100 -1.891068715 S3 cz Con Female 67.01 5.39#dim(erp)
frames <- seq(0,1001,4)preprocessing
1.Subject selection
2.electrode selection
3.conjuction selection (individual and electrode)
4.aggregate
- Subject selection
ind_select <- function(erpdata,frames,erpcol,subcol,elecol,othvarcol=NULL,
select_sub,...){
# need some check function
data <- erpdata
num <- length(select_sub)
data_list <- list()
for (i in 1:num){
data_new <- subset(data,data[,subcol]==select_sub[i])
data_list[[i]] <- data_new
}
data_select <- data_list[[1]]
if (num != 1) {
for (i in 2:num) {
data_select <- rbind(data_select,data_list[[i]])
}
}
return (data_select)
}Subject selection_example
S <- ind_select(erp,frames,1:251,252,253,254:257,
select_sub=c("S1","S10","S8"))
dim(S)## [1] 90 2572.Electrode Selection
elect_select <- function(erpdata,frames,erpcol,subcol,elecol,othvarcol=NULL,
select_ele,...){
# some check function
data <- erpdata
num <- length(select_ele)
data_list <- list()
for (i in 1:num){
data_new <- subset(data,data[,elecol]==select_ele[i])
data_list[[i]] <- data_new
}
data_select <- data_list[[1]]
if (num != 1){
for (i in 2:num){
data_select <- rbind(data_select,data_list[[i]])
}
}
return (data_select)
}Electrode Selection _example
el <- elect_select(erp,frames,1:251,252,253,254:257,
select_ele=c("cz","fz"))
dim(el)## [1] 200 257- Conjunction_Selection
conj_select <- function(erpdata,frames,erpcol,subcol,elecol,othvarcol=NULL,
select_sub,
select_ele,...){
data <- erpdata
ind_data <- ind_select(data,frames,erpcol,subcol,elecol,othvarcol,
select_sub = select_sub)
ind_ele_data <- elect_select(ind_data,frames,erpcol,subcol,elecol,othvarcol,
select_ele = select_ele)
return(ind_ele_data)
}Conjunction_Selection_Example
SE <- conj_select(erp,frames,1:251,252,253,254:257,
select_sub=c("S1","S2"),
select_ele=c("cz","fz"))
dim(SE)## [1] 40 257- Aggregate Raw data
aggraw <- function(erpdata,frames,erpcol,subcol,elecol,othvarcol=NULL,
aggvarcol=c(subcol),
fun=mean...){
#some check function
data <- erpdata
agglength <- length(aggvarcol)
aggvar_list <- list(data[,aggvarcol[1]])
if (agglength > 1){
for (i in 2:agglength ){
aggvar_list <- append(aggvar_list,list(data[,aggvarcol[i]]))
}
}
aggdata <- aggregate(data[,c(erpcol,aggvarcol)],by=aggvar_list,fun)
aggdata <- aggdata[,1:(agglength+length(erpcol))]
for (i in 1: agglength){
colnames(aggdata)[i] <- colnames(data)[aggvarcol[i]]
}
return(aggdata)
}Aggregate Raw data _Example
A <-aggraw(erp,frames,1:251,252,253,othvarcol=254:257,
aggvarcol=c(253,252),
fun=median)
# you could put all element you want to aggregate in the argument erpcol
head(A)[,1:7]## Electrode Subject X0 X4 X8 X12
## 1 cz S1 -0.044902303 0.04303762 -0.04744685 -0.08189084
## 2 fz S1 -0.004841731 0.06401845 0.10148365 -0.17240361
## 3 pz S1 0.741473704 1.19190216 1.33562088 0.84747371
## 4 cz S10 -0.333857536 -0.44890246 -0.37539005 -0.24336576
## 5 fz S10 0.344259185 -0.30369284 0.47529451 0.74365250
## 6 pz S10 -0.070352518 -0.25235034 0.19857943 0.09318284
## X16
## 1 -0.2842869
## 2 -0.3055670
## 3 0.8820590
## 4 -0.6097555
## 5 -0.2303083
## 6 0.4208192EDA
The function “erp_ggplot”" uses ggplot2 grammar, so it has several flexible options. Users have to enter…
1. erpdata = whole data
2. frames = frames
3. erpcol : the column of erp signal 4. subcol = the number of “subject”" column, ex: subject is in column 252, enter 252
5. elecol = the numer of eletrode column
- othvarcol = the number of other variables (optional)
- outlinesub = which subject the user want to outline (optional)
outcolor = the outline color (optional)
- select_sub = the subject you would like to put in the plot (optional)
select_ele = the eletrode you would like to put in the plot (optional)
# group comparison please put in data after aggregate_raw
erp_ggplot <- function(erpdata,frames=NULL,erpcol,subcol,elecol,othvarcol=NULL,
outlinesub=NULL,outcolor="red",
select_sub=c(NULL),
select_ele=c(NULL),...){
#some check functions
#if (is.null(frames) == FALSE)
# if (length(frames) != (ncol(data)-1-length(othvarcol)))
# stop(paste("frames should be either null or of length",
# (ncol(data)-1-length(othvarcol))))
#if (is.null(frames) == FALSE) {
# if (any(frames != sort(frames)))
# stop("frames should be an ascending sequence of integers")
# }
#if (is.null(frames))
# frames = 1:(ncol(data)-1-length(othvarcol))
#selection
if (is.null(select_sub)==FALSE & is.null(select_ele)==FALSE){
data <- conj_select(erpdata,frames,erpcol,subcol,elecol,othvarcol,
select_sub=select_sub,
select_ele=select_ele)
} else if (is.null(select_sub)==FALSE ) {
data <- ind_select(erpdata,frames,erpcol,subcol,elecol,othvarcol,
select_sub=select_sub)
} else if (is.null(select_ele)==FALSE) {
data <- elect_select(erpdata,frames,erpcol,subcol,elecol,othvarcol,
select_ele=select_ele)
} else {
data <- erpdata
}
# plot
subvar <- variable.names(data)[subcol]
data$groupvar <- rownames(data)
datalong <- melt(data,
id=c(variable.names(data)[c(subcol,elecol,othvarcol)],
"groupvar"))
datalong <- datalong[order(datalong$groupvar),]
datalongorder <- datalong
datalongorder$frames <- rep(frames,length(datalongorder[,1])/length(frames))
if (is.null(outlinesub) == FALSE){ # how to outline several subjects (and color)
data2 <- subset(datalongorder,datalongorder[,1]==outlinesub)
plot <- ggplot(datalongorder,
aes(x=frames,y=value,group=groupvar,...))+
geom_line()+
geom_line(data=data2,aes(x=frames,y=value),col=outcolor)
# need warning for covering geom_line()
} else {
plot <- ggplot(datalongorder,
aes(x=frames,y=value,group=groupvar))+
geom_line()
}
return(plot)
}Example
#Although Full data (all trials) work fine, I recommened that the argument erpdata should be a single subject data or aggregate data.
# Aggregate Data
subjectmean <-aggraw(erp,frames,c(1:251,256,257),252,253,othvarcol=254:257,
aggvarcol=c(253,252,254,255),
fun=mean)
#head(subjectmean[order(subjectmean$Subject),],50)[,c(1:5,255:257)]
erp_ggplot(subjectmean,frames,erpcol=5:256,subcol=2,elecol=1,othvarcol=c(3,4,256:257),
outlinesub="S9",outcolor = "blue",
# you could highlight a single subject
select_sub = c("S5","S8","S9"),
# you can choose which subject you want to show in plot
select_ele = c("cz","fz"))+
# you can choose which electrode you want to show in plot
facet_grid(Electrode~Condition)+
theme_bw()+
stat_summary(aes(group=NULL),fun.y = "mean",
colour = "red", size = 0.5, geom = "line")
# put on the summary line NOTE: need a group=NULL argument
erp_ggplot(subjectmean,frames,erpcol = 5:256,subcol=2,elecol=1,othvarcol=c(3:4,256:257))+
geom_line(aes(col=Scores))+ # you could specify the color argument
facet_grid(Electrode~Condition)+
theme_bw()
# Single Subject : S1
S1 <- ind_select(erp,frames,1:251,252,253,254:257,
select_sub=c("S1"))
#head(S1)
erp_ggplot(S1,frames,erpcol=1:251,subcol = 252,elecol=253,othvar=254:257)+
geom_line(aes(col=RT))+
theme(legend.position="none")+
xlim(-500,1500)+
ylim(-30,30)+
labs(list(title = "Flexible",x="time",y="signal"))+
facet_grid(Condition~.)+
stat_summary(aes(group=NULL),fun.data="mean_cl_boot", colour="red",
geom="crossbar", width=0.2)
# iI am not really sure whether the stat_summary function (mean_cl_boot) works right(or fine) or not, but I still provide for this demo.tests plot
- warning on or off
- some check function
erp_tests_ggplot
User have to enter
1. test_rst = test results of the test function (erpfatest, erpavetest, erptest, gbtest)
2. frames = frames
3. cor = show correlation or not (TRUE or FALSE)
4. sigp = Do you want to show the additional points indicate the significant points ? (TRUE or FALSE)
5. sigpcol = the color of additional points (when sigp = TRUE, default = “red”)
6. sigp_yloc = the location of additional points (when sigp = TRUE, default = 0)
erp_test_ggplot <- function(tests_rst,frames,cor=FALSE,
sigp=FALSE,sigpcol="red",sigp_yloc=0,...){
# add some check function
options(warn=-1)
# Whether check off the warn is good or not is an open question.
data <- data.frame(signal=as.numeric(tests_rst$signal))
data$frames = frames
data$significant <- ifelse(data$frames %in% data$frames[tests_rst$significant],
"sig","non-sig")
data$sign_frames <- ifelse(data$frames %in% data$frames[tests_rst$significant],
data$frames,NA)
data$group <- rep(sigp_yloc,length(data$signal))
data$r2 <- tests_rst$r2
if (cor == TRUE){
plot <- ggplot(data,aes(x=frames,y=sign(signal)*sqrt(r2),group=group))+
geom_line(aes(col=significant))+
scale_colour_manual(values=c("black","red"))
} else {
plot <- ggplot(data,aes(x=frames,y=signal,group=group))+
geom_line(aes(col=significant))+
scale_colour_manual(values=c("black","red"))
}
if (sigp == TRUE){
plot <- plot+
geom_point(aes(x=sign_frames,y = group),shape=8,col=sigpcol)
}
return(plot)
}Example 1
tests = erpfatest(erpcz[,1:251],design=model.matrix(~Subject+Instruction,data=erpcz),
design0=model.matrix(~Subject,data=erpcz),nbf=3,
s0=c(1:50,226:251),min.err=1e-01)
frames = seq(0,1001,4)
#Original Plot
#plot(frames,tests$signal,type="l",xlab="Time (ms)",ylab="Difference ERP curves")
#points(frames[tests$significant],rep(0,length(tests$significant)),pch=16,col="blue")
#title("Paired comparison at electrode CZ")#New Plot
erp_test_ggplot(tests,frames,cor=F) # cor = FALSE 
# the simplest example
erp_test_ggplot(tests,frames,cor=F,
sigp=T,#Significant points are added.
sigpcol = "yellow",#The color of Significant points
sigp_yloc = -10) #The y location of significant points 
Example 2 : Correlation Setting
tests = erpfatest(simerp[,1:251],design=model.matrix(~y,data=simerp),
nbf=3,min.err=1e-01,maxiter=10)
frames = seq(0,1001,4)
#plot(frames,sign(tests$signal)*sqrt(tests$r2),type="l",
# xlab="Time (ms)",ylab="Correlation",ylim=c(-1,1))
#points(frames[tests$significant],rep(-1,length(tests$significant)),
# pch=16,col="blue")
#title("Simulation")#New Plot
erp_test_ggplot(tests,frames,cor=T,sigp=T,"yellow",-1)+ #Correlation should be TRUE (cor=T).
scale_colour_manual(values=c("blue","green"),
# change the color (cover the original one)
name="Sig",label=c("Y","N"))+
# change the label (cover the original one)
ylim(-2,2)+
xlim(-100,1100)+
theme_bw()+
theme(legend.position="none")+
labs(list(x="Time (ms)",y="Correlation",title="Simulation"))
group comparison
- bootinterval : some check function, change (mean)
erp_grpcp_bootintval
erp_bootintval <- function(erpdata,frames,erpcol,subcol,elecol,othvarcol=NULL,
cpvarcol,
fun=mean,
bootnum=300,
bootintval=c(.025,.975)){
# need some check function
data <- erpdata
if ( class(data[,cpvarcol]) != "factor"){
stop("cpvar should be factor")
}
factor_num <- length(levels(data[,cpvarcol]))
data_list <- list()
#funname <- as.character(fun) # some problem
for ( i in 1 : factor_num ) {
data_list[[i]] <- subset(data,data[,cpvarcol]==(levels(data[,cpvarcol])[i]))
data_list[[i]] <- data_list[[i]][,erpcol]
}
quan_list <- list()
for ( j in 1:factor_num) {
quan_list[[j]] <- as.data.frame(matrix(NA,3,length(erpcol)))
rownames(quan_list[[j]]) <- c("FUN","Q1","Q2") #how to change,# some problem
colnames(quan_list[[j]]) <- colnames(data[,erpcol])
quan_list[[j]][1,] <- apply(data_list[[j]][,1:length(erpcol)],2,fun)
}
for (n in 1 : factor_num){
for (k in 1:length(erpcol)) {
mat <- matrix(NA,bootnum,1)
for (m in 1 : bootnum){
new <- sample(data_list[[n]][,k],replace=T)
mat[m,1] <- fun(new)
}
quan_list[[n]][2,k] <- quantile(mat[,1],bootintval)[1]
quan_list[[n]][3,k] <- quantile(mat[,1],bootintval)[2]
}
}
names(quan_list) <- levels(data[,cpvarcol])
return(quan_list)
}erp_boot_ggplot <- function(erpdata,frames,erpcol=NULL,subcol,elecol,othvarcol=NULL,
cpvarcol,
quan_list=NULL,
fun=mean,bootnum=300,bootintval=c(.025,.975),alpha=0.3) {
#need some check function
data <- erpdata
if (is.null(quan_list)){
quan_list <- erp_bootintval(data,frames,erpcol,subcol,elecol,othvarcol,
cpvarcol,
fun,bootnum,bootintval)
}
factor_num <- length(levels(data[,cpvarcol]))
cpvar_names <- colnames(data)[cpvarcol]
data_plot <- as.data.frame(t(quan_list[[1]]))
for (i in 2:factor_num){
data_plot <- rbind(data_plot,as.data.frame(t(quan_list[[i]])))
}
data_plot$frames <- (rep(frames,factor_num))
data_plot[,5] <- sort(rep(names(quan_list),length(frames)))
colnames(data_plot)[5] <- cpvar_names
data_plot$group <- data_plot[,5]
plot <- ggplot(data_plot,aes(x=frames,group=group))+
geom_ribbon(aes(x=frames, ymax=Q2, ymin=Q1,fill=group), alpha=alpha)+ # set alpha
geom_line(aes(y = FUN,col=group))+
labs(y="Signal") # Need some changes ?
return(plot)
}frames <- seq(0,1001,4)
subjectmean <-aggraw(erp,frames,c(1:251,256,257),252,253,othvarcol=254:257,
aggvarcol=c(253,252,254,255),
fun=mean)
#head(subjectmean)[,5:256]
submeancz <- elect_select(subjectmean,frames,
erpcol=5:255,subcol=2,elecol=1,othvarcol=c(3:4,256,257),
select_ele = "cz")
#submeancz[,c(1:6,255:257)]
submeancz$Condition <- as.factor(submeancz$Condition)
quan <- erp_bootintval(submeancz,frames,
erpcol=5:255,subcol=2,elecol=1,othvarcol=c(3:4,256,257),
cpvarcol = 3,
fun=median,bootnum=10,bootintval = c(0.01,0.99))#report a quantile list
erp_boot_ggplot(submeancz,frames,
erpcol=5:255,subcol=2,elecol=1,othvarcol=c(3:4,256,257),
cpvarcol = 3,
fun=mean,bootnum=10,bootintval = c(0.05,0.95),
alpha=0.5)
#erpcz
erpcz$Ele <- "cz"
erp_boot_ggplot(erpcz,frames,
erpcol=1:251,subcol=252,elecol=254,othvarcol = 253,
cpvarcol = 253,
fun=median,bootnum=10,bootintval = c(0.025,0.975),
alpha=0.3)
quan
plot(x=frames,y=quan[[1]][1,],type="l",ylim=c(-10,10))
lines(x=frames,y=quan[[1]][2,],type="l",col="red")
lines(x=frames,y=quan[[1]][3,],type="l",col="red")
lines(x=frames,y=quan[[2]][1,],type="l",col="blue")
lines(x=frames,y=quan[[2]][2,],type="l",col="green")
lines(x=frames,y=quan[[2]][3,],type="l",col="green")