ERP package
Outline
data
1. ERPdata : maybe a formal data for package
2. NIRSdata
preprocessing
1. list2df
2. AddExpCondSub
3. ind_select
4. ele_select
5. conj_select
6. aggraw
EDA
1.edaplot
2.bootplot
tests
1.multichtest
2.testplot
3.signalplot
ERP only
1. scalpplot
2. topograph
Data
1. erpcz : a single electrode, 2 within subject conditions, 20 subjects data.
2. simerp: a single electrode, 1 condition, 20 subjects, a score between subject (y).
3. ERPsets : a list dataset from erpR package: 20 subjects , 2 conditions (40 lists).
library(ERP);require(mnormt);require(fdrtool);library(reshape);library(ggplot2);library(dplyr)
library(erpR);require(akima);library(reshape2)
data(ERPsets)
load("ERPdata.RData")preprocessing
1. list2df
2. AddExpCondSub
3. ind_select
4. ele_select
5. conj_select
6. aggraw
1. list2df : It is a function that could transformed a list to a dataframe for our package.
###list2data
list2data <- function(list_data,frames){
# some check function
for (i in 1:length(list_data)){
list_data[[i]]$frames <-frames
list_data[[i]]$list.name <- names(list_data)[i]
}
Data_list <- lapply(list_data,melt,id=c("frames","list.name"))
Data_list <- lapply(Data_list,reshape,
timevar = "frames",idvar = c("variable","list.name"),
direction = "wide")
Ana_data <- Data_list[[1]]
for (i in 2:length(ERPsets)){
Ana_data <- rbind(Ana_data,Data_list[[i]])
}
colnames(Ana_data)[2] <- "Channel"
rownames(Ana_data) <- 1:dim(Ana_data)[1]
return(Ana_data)
}2. AddExpCondSub :
After transforming,
the user have to append Subject, Condition, and other variables by themself
(depend on the design).
But if your “List Name” have a common grammar : Exp_Condition_Subject,
AddCondSub function could help.
AddExpCondSub <- function(data,list.name_col){
# some check functions
dta <- data
num <- dim(dta)[2]
for (i in 1 : dim(dta)[1]){
pos_list = gregexpr("_", data[i,list.name_col])
pos1 <- pos_list[[1]][1]
pos2 <- pos_list[[1]][2]
dta[i,num+1] <- substr(dta[i,list.name_col],
1,pos1-1)
dta[i,num+2] <- substr(dta[i,list.name_col],
pos1+1,pos2-1)
dta[i,num+3] <- substr(dta[i,list.name_col],
pos2+1,nchar(dta[i,list.name_col]))
}
colnames(dta)[(num+1):(num+3)] <- c("Experiment","Condition","Subject")
return(dta)
}ERP_df <- list2data(list_data=ERPsets, # input list_data
frames= 1:426) # input frames
ERP_df <- AddExpCondSub(data = ERP_df, # input the transformed data
list.name_col = 1) # list.name column
#ERP_df <- ERP_df[,-1] # remove list name column, Depend on user's need
head(ERP_df[,c(1:3,427:430)],3) # look at the data list.name Channel value.1 value.425 value.426 Experiment
1 Exp1_word_subj1 Fp1 0.36995740 14.160980 9.971283 Exp1
2 Exp1_word_subj1 Fp2 -0.02759907 13.750350 9.678596 Exp1
3 Exp1_word_subj1 F3 0.42971300 6.269233 4.413941 Exp1
Condition
1 word
2 word
3 word3. ind_select
select a single subject or several subjects.
ind_select <- function(data,frames,datacol,subcol,chcol=NULL,othvarcol=NULL,
select_sub,...){
# need some check function
dta <- data
num <- length(select_sub)
data_list <- list()
for (i in 1:num){
data_new <- subset(dta,dta[,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]])
}
}
rownames(data_select) <- 1:dim(data_select)[1]
return (data_select)
}
S1S2 <- ind_select(data = ERPdata,
frames = 1:426,
datacol = 2:427,
subcol = 430,
chcol = 1,
othvarcol = c(428:429,431:432),
select_sub=c("subj1","subj2"))
dim(S1S2) # 2 conditions *34 electrodes * 2 subjects [1] 136 4324. ch_select
select a single channel or several channels
ch_select <- function(data,frames,datacol,subcol=NULL,chcol,othvarcol=NULL,
select_ch,...){
# some check function
dta <- data
num <- length(select_ch)
data_list <- list()
for (i in 1:num){
data_new <- subset(dta,dta[,chcol]==select_ch[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]])
}
}
rownames(data_select) <- 1:dim(data_select)[1]
return (data_select)
}Fp <- ch_select(data = ERPdata,
frames = 1:426,
datacol = 2:427,
subcol = 430,
chcol = 1,
othvarcol = c(428:429,431:432),
select_ch=c("Fp1","Fp2"))
dim(Fp) # 2 conditions * 20 subjects * 2 electrode[1] 80 4325. conj_select
select channels and subjects
conj_select <- function(data,frames,datacol,subcol,chcol,othvarcol=NULL,
select_sub,
select_ch,...){
# some check function
dta <- data
ind_data <- ind_select(dta,frames,datacol,subcol,chcol,othvarcol,
select_sub = select_sub)
ind_ele_data <- ch_select(ind_data,frames,datacol,subcol,chcol,othvarcol,
select_ch = select_ch)
rownames(ind_ele_data) <- 1 :dim(ind_ele_data)[1]
return(ind_ele_data)
}SE <- conj_select(data = ERPdata,
frames = 1:426,
datacol = 2:427,
subcol = 430,
chcol = 1,
othvarcol = c(428:429,431:432),
select_sub=c("subj1"),
select_ch=c("Fp1","Fp2"))
dim(SE) # 1 sub * 2 Channel * 2 Condition[1] 4 4326. Aggregate Data:
We could aggregate data by the column like conditions, subjects or channel.
aggraw <- function(data,frames,datacol,subcol=NULL,chcol=NULL,othvarcol=NULL,
aggvarcol,
fun=mean,...){
#some check function
options(warn=-1)
# should close the warnings?
dta <- data
agglength <- length(aggvarcol)
aggvar_list <- list(dta[,aggvarcol[1]])
if (agglength > 1){
for (i in 2:agglength ){
aggvar_list <- append(aggvar_list,list(dta[,aggvarcol[i]]))
}
}
aggdata <- aggregate(dta[,datacol],by=aggvar_list,
fun,...)
aggdata <- aggdata[,1:(agglength+length(datacol))]
for (i in 1: agglength){
colnames(aggdata)[i] <- colnames(dta)[aggvarcol[i]]
}
rownames(aggdata) <- 1:dim(aggdata)[1]
return(aggdata)
}A <-aggraw(data = ERPdata,
frames = 1:426,
datacol = 2:427, #you could put all elements you want to aggregate in datacol argument
subcol = 430,
chcol = 1,
othvarcol = c(428:429,431:432),
aggvarcol= c(1,429),
fun=median) # mean median and so on
tail(A[,1:6]);dim(A) Channel Condition value.1 value.2 value.3 value.4
63 OZ word 0.3915177 0.3250047 -0.01918097 -0.00771116
64 FT7 word 0.4533806 0.3932142 -0.23910702 -0.41419665
65 FT8 word 0.9787850 1.3420155 1.02323850 1.11230300
66 A2 word 0.3218264 0.8343187 0.93161350 0.87334650
67 VEOG1 word 15.2228400 22.3858300 20.08051500 21.02715500
68 HEOG1 word 3.4621365 6.0697715 5.82959150 6.68893100[1] 68 428EDA
1.edaplot
2.bootplot
1.edaplot
The function “edaplot”" uses ggplot2 grammar, so it has several flexible options.
# group comparison please put in data after aggregate_raw
edaplot <- function(data,frames=NULL,datacol,subcol=NULL,chcol=NULL,othvarcol=NULL,
outlinesub=NULL,outcolor="red",
select_sub=c(NULL),
select_ch=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_ch)==FALSE){
dta <- conj_select(data = data,frames = frames,
datacol,subcol,chcol,othvarcol,
select_sub=select_sub,
select_ch=select_ch)
} else if (is.null(select_sub)==FALSE ) {
dta <- ind_select(data=data,frames=frames,
datacol,subcol,chcol,othvarcol,
select_sub=select_sub)
} else if (is.null(select_ch)==FALSE) {
dta <- ch_select(data,frames,
datacol,subcol,chcol,othvarcol,
select_ch=select_ch)
} else {
dta <- data
}
# plot
subvar <- variable.names(dta)[subcol]
dta$groupvar <- rownames(dta)
datalong <- melt(dta,
id=c(variable.names(dta)[c(subcol,chcol,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)
}#Although Full data (all trials) work fine, I recommened that the argument erpdata should be a single subject data or an aggregate data.
edaplot(ERPdata,
frames = 1:426,
datacol=2:427,
subcol=430,
chcol=1,
othvarcol=c(428:429,431:432),
outlinesub="subj9",outcolor = "blue",# highlight a single subject with the color you want
select_sub = c("subj5","subj8","subj9"), # choose the subject you want to show(optional)
select_ch = c("F3","F4"))+# choose the channel you want to show(optional)
facet_grid(Channel~Condition)+ # other flexible setting
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
edaplot(ERPdata,
frames=1:426,
datacol = 2:427,subcol=430,chcol=1,othvarcol=c(428:429,431:432),
select_ch = c("F3","F4","Cz"))+
geom_line(aes(col=IQ))+
facet_grid(Condition~Channel)+
theme_bw()+
theme(legend.position="none")+
xlim(-100,500)+
ylim(-15,15)+
labs(list(title = "Flexible",x="time",y="signal"))
2.bootplot
bootplot use a bootstaping method to draw the confidence interval on two (or more) compared
variable.
bootplot <- function(data,frames,datacol,subcol=NULL,chcol=NULL,othvarcol=NULL,
cpvarcol,
fun=mean,
bootnum=10,
bootintval=c(.025,.975),
alpha=0.3){
# need some check function
bootstrap <- function(x,bootnum,bootfun,bootintval=c(0.025,0.975),quantilenum=1){
mat <- matrix(NA,bootnum,1)
for (m in 1 : bootnum){
new <- sample(x,replace=T)
mat[m,1] <- bootfun(new)
}
return(quantile(mat,bootintval,na.rm = T)[quantilenum])
}
dta <- data
funname <- as.character(substitute(fun))
dta[,cpvarcol] <- as.factor(dta[,cpvarcol])
dta[,chcol] <- as.factor(dta[,chcol])
data_fun <- aggraw(dta,frames,
datacol,subcol,chcol,othvarcol,
aggvarcol = c(chcol,cpvarcol),
fun=fun)
data_Q1 <- aggraw(dta,frames,
datacol,subcol,chcol,othvarcol,
aggvarcol = c(chcol,cpvarcol),
fun=bootstrap,
bootnum=bootnum,
bootfun=fun,
bootintval=bootintval,
quantilenum=1)
data_Q2 <- aggraw(dta,frames,
datacol,subcol,chcol,othvarcol,
aggvarcol = c(chcol,cpvarcol),
fun=bootstrap,
bootnum=bootnum,
bootfun=fun,
bootintval=bootintval,
quantilenum=2)
data_fun_long <- melt(data_fun,id=c(colnames(dta)[chcol],colnames(dta)[cpvarcol]))
data_fun_long <- data_fun_long[order(data_fun_long[,1],
data_fun_long[,2],
data_fun_long[,3]),]
colnames(data_fun_long)[4] <- "FUN"
data_Q1_long <- melt(data_Q1,id=c(colnames(dta)[chcol],colnames(dta)[cpvarcol]))
data_Q1_long <- data_Q1_long[order(data_Q1_long[,1],
data_Q1_long[,2],
data_Q1_long[,3]),]
data_Q2_long <- melt(data_Q2,id=c(colnames(dta)[chcol],colnames(dta)[cpvarcol]))
data_Q2_long <- data_Q2_long[order(data_Q2_long[,1],
data_Q2_long[,2],
data_Q2_long[,3]),]
data_fun_long$Q1 <- data_Q1_long[,4]
data_fun_long$Q2 <- data_Q2_long[,4]
data_fun_long$frames <- c(rep(frames,(dim(data_fun_long)[1]/length(frames)))) #change
colnames(data_fun_long)[1:2] <- c("Channel","Condition")
plot <- ggplot(data_fun_long,aes(x=frames,group=Condition))+
geom_ribbon(aes(x=frames, ymax=Q2, ymin=Q1,fill=Condition), alpha=alpha)+ # set alpha
#geom_errorbar(aes(x=frames, ymax=Q2, ymin=Q1,fill=group,col=group,alpha=alpha))+
geom_line(aes(y = FUN,col=Condition))+
labs(y="Signal")+# Need some changes ?
facet_wrap(~Channel)
return(plot)
}
bootplot(ERPdata, # input the data
frames=1:426,
datacol=2:427,
subcol=430,
chcol=1,
othvarcol=c(428:429,431:432),
cpvarcol=429, # Important : the column of that single variable you want to compare
fun=mean, # the function use to draw boot interval and line
bootnum=10, # bootsraping number
bootintval=c(.025,.975), # interval
alpha=0.5)+ # the value of alpha on the plot
ylim(-5,5)+
scale_fill_manual(values=c("red","blue"),
name="Word or Non Word",label=c("No","Yes"))+
scale_colour_manual(values=c("red","blue"),
name="Word or Non Word",label=c("No","Yes"))
tests
1.multichtest
2.testplot
3.signalplot
1.multichtest
multichtest <- function(data,datacol,chcol,
testtype="erpfatest",
# do not specify model.matrix like original test function###
design_model,
# do not specify model.matrix like original test function
design0_model=NULL,...){
dta <- data
dta[,chcol]=as.factor(dta[,chcol])
levelnum <- length(levels(dta[,chcol]))
if (levelnum == 1) {
design <- model.matrix(design_model,data=data)
if (is.null(design0_model)==F){
design0 <- model.matrix(design0_model,data=data)
}
if ( testtype == "erpavetest" ){
test_list <- erpavetest(dta[,datacol],design,design0,...)
}
if (testtype == "erpfatest") {
test_list <- erpfatest(dta[,datacol],design,design0,...)
}
if (testtype == "erptest") {
test_list <- erptest(dta[,datacol], design, design0,...)
}
if (testtype == "gbtest") {
test_list <- gbtest(dta[,datacol], design, design0,...)
}
} else {
test_list=list()
dta_list=list()
for (i in 1:levelnum) {
dta_list[[i]] <- subset(dta,dta[,chcol]==(levels(dta[,chcol])[i]))
design <- model.matrix(design_model,data=dta_list[[i]])
if (is.null(design0_model)==F){
design0 <- model.matrix(design0_model,data=dta_list[[i]])
}
if ( testtype == "erpavetest" ){
test_list[[i]] <- erpavetest(dta_list[[i]][,datacol],
design,design0,...)
}
if (testtype == "erpfatest") {
test_list[[i]] <- erpfatest(dta_list[[i]][,datacol],
design,design0,...)
}
if (testtype == "erptest") {
test_list[[i]] <- erptest(dta_list[[i]][,datacol],
design, design0,...)
}
if (testtype == "gbtest") {
test_list[[i]] <- gbtest(dta_list[[i]][,datacol],
design, design0,...)
}
names(test_list)[i] <- levels(dta[,chcol])[i]
}
}
return(test_list)
}
ERP_test <- multichtest(ERPdata,2:427,chcol=1,testtype="erpfatest",
design_model=(~Subject+Condition),
design0_model=(~Subject))2. testplot
plot the difference curve and highlight the significant value
testplot <- function(tests_rst,frames,multi=FALSE,
siglinecol="red",
cor=FALSE,
sigp=FALSE,sigpcol="red",sigpshape=8,sigp_yloc=0,...){
# add some check function
options(warn=-1)
# Whether check off the warn is good or not is an open question.
if (multi==FALSE){
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",siglinecol))
} else {
plot <- ggplot(data,aes(x=frames,y=signal,group=group))+
geom_line(aes(col=significant))+
scale_colour_manual(values=c("black",siglinecol))
}
if (sigp == TRUE){
plot <- plot+
geom_point(aes(x=sign_frames,y = group),shape=sigpshape,col=sigpcol)
}
} else {
listlen <- length(tests_rst)
data_list <- list()
for (k in 1:listlen){
data_list[[k]] <- data.frame(signal=as.numeric(tests_rst[[k]]$signal))
data_list[[k]]$frames = frames
data_list[[k]]$significant <- ifelse(data_list[[k]]$frames %in% data_list[[k]]$frames[tests_rst[[k]]$significant],"sig","non-sig")
data_list[[k]]$sign_frames <- ifelse(data_list[[k]]$frames %in% data_list[[k]]$frames[tests_rst[[k]]$significant],data_list[[k]]$frames,NA)
data_list[[k]]$group <- rep(sigp_yloc,length(data_list[[k]]$signal))
data_list[[k]]$r2 <- tests_rst[[k]]$r2
data_list[[k]]$Channel <- names(tests_rst)[k]
}
data_plot <- data_list[[1]]
for (j in 2 : listlen){
data_plot <- rbind(data_plot,data_list[[j]])
}
if (cor == TRUE){
plot <- ggplot(data_plot,aes(x=frames,y=sign(signal)*sqrt(r2),group=group))+
geom_line(aes(col=significant))+
scale_colour_manual(values=c("black",siglinecol))+
facet_wrap(~Channel)
} else {
plot <- ggplot(data_plot,aes(x=frames,y=signal,group=group))+
geom_line(aes(col=significant))+
scale_colour_manual(values=c("black",siglinecol))+
facet_wrap(~Channel)
}
if (sigp == TRUE){
plot <- plot+
geom_point(aes(x=sign_frames,y = group),shape=sigpshape,col=sigpcol)
}
}
return(plot)
}testplot(ERP_test,
frames=1:426,
multi=T, # Multiple Plot = TRUE (several Channels)
siglinecol="red",
cor=FALSE,
sigp=T,sigpcol="red",sigpshape=8,sigp_yloc=-5)+
ylim(-5,5)+
xlim(-100,500)+
theme_bw()+
theme(legend.position="none")+
labs(list(x="Time (ms)",y="Correlation",title="Simulation"))
##### 3. signalplot
##### plot original signal curve and highlight the significant location
##### and it is possible to add some bootstrap interval
signalplot <- function(tests_rst,
multi=F,
data,
frames,
datacol,subcol=NULL,chcol=NULL,othvarcol=NULL,
cpvarcol,
sigp_yloc=0,
sigpcol="red",
sigpshape=8,
wantbootplot=FALSE,
fun=mean,
bootnum=10,
bootintval=c(.025,.975),
alpha=0.3) {
# some check functions
bootstrap <- function(x,bootnum,bootfun,bootintval=c(0.025,0.975),quantilenum=1){
mat <- matrix(NA,bootnum,1)
for (m in 1 : bootnum){
new <- sample(x,replace=T)
mat[m,1] <- bootfun(new)
}
return(quantile(mat,bootintval,na.rm = T)[quantilenum])
}
dta <- data
if (multi == F ){
data_sign <- data.frame(frames=frames)
data_sign$significant <-
ifelse(data_sign$frames %in% data_sign$frames[tests_rst$significant],
"sig","non-sig")
data_sign$sign_frames <-
ifelse(data_sign$frames %in% data_sign$frames[tests_rst$significant],
data_sign$frames,NA)
data_sign$group <- rep(sigp_yloc,length(frames))
plot <- bootplot(dta,frames,
datacol,subcol,chcol,othvarcol=othvarcol,
cpvarcol,
fun,bootnum,bootintval=bootintval,
alpha=ifelse(wantbootplot==T,alpha,0))+
geom_point(data=data_sign,
mapping = aes(x=sign_frames,y = group,group=NULL),
shape=sigpshape,col=sigpcol)
}
if (multi == T){
listlen <- length(tests_rst)
test_list = list()
for (i in 1: listlen){
test_list[[i]] <- data.frame(frames=frames)
test_list[[i]]$significant <-
ifelse(test_list[[i]]$frames %in% test_list[[i]]$frames[tests_rst[[i]]$significant],"sig","non-sig")
test_list[[i]]$sign_frames <-
ifelse(test_list[[i]]$frames %in% test_list[[i]]$frames[tests_rst[[i]]$significant],test_list[[i]]$frames,NA)
test_list[[i]]$group <- rep(sigp_yloc,length(frames))
test_list[[i]]$Channel <- names(tests_rst)[i]
}
test_plot <- test_list[[1]]
for (k in 2:listlen){
test_plot <- rbind(test_plot,test_list[[k]])
}
funname <- as.character(substitute(fun))
dta[,cpvarcol] <- as.factor(dta[,cpvarcol])
dta[,chcol] <- as.factor(dta[,chcol])
data_fun <- aggraw(dta,frames,
datacol,subcol,chcol,othvarcol,
aggvarcol = c(chcol,cpvarcol),
fun)
data_Q1 <- aggraw(dta,frames,
datacol,subcol,chcol,othvarcol,
aggvarcol = c(chcol,cpvarcol),
fun=bootstrap,
bootnum=bootnum,
bootfun=fun,
bootintval=bootintval,
quantilenum=1)
data_Q2 <- aggraw(dta,frames,
datacol,subcol,chcol,othvarcol,
aggvarcol = c(chcol,cpvarcol),
fun=bootstrap,
bootnum=bootnum,
bootfun=fun,
bootintval=bootintval,
quantilenum=2)
data_fun_long <- melt(data_fun,id=c(colnames(dta)[chcol],colnames(dta)[cpvarcol]))
data_fun_long <- data_fun_long[order(data_fun_long[,1],
data_fun_long[,2],
data_fun_long[,3]),]
colnames(data_fun_long)[4] <- "FUN"
data_Q1_long <- melt(data_Q1,id=c(colnames(dta)[chcol],colnames(dta)[cpvarcol]))
data_Q1_long <- data_Q1_long[order(data_Q1_long[,1],
data_Q1_long[,2],
data_Q1_long[,3]),]
data_Q2_long <- melt(data_Q2,id=c(colnames(dta)[chcol],colnames(dta)[cpvarcol]))
data_Q2_long <- data_Q2_long[order(data_Q2_long[,1],
data_Q2_long[,2],
data_Q2_long[,3]),]
data_fun_long$Q1 <- data_Q1_long[,4]
data_fun_long$Q2 <- data_Q2_long[,4]
data_fun_long$frames <- c(rep(frames,(dim(data_fun_long)[1]/length(frames))))#change
colnames(data_fun_long)[1:2] <- c("Channel","Condition")
data_plot <- merge(data_fun_long,test_plot,by =c("Channel","frames"))
plot <- ggplot(data_plot,aes(x=frames,group=Condition))+
geom_ribbon(aes(x=frames, ymax=Q2, ymin=Q1,fill=Condition)
,alpha=ifelse(wantbootplot==T,alpha,0))+
geom_line(aes(y = FUN,col=Condition))+
facet_wrap(~Channel)+
geom_point(aes(x=sign_frames,y = group,group=NULL),
shape=sigpshape,col=sigpcol)+
labs(y="Signal")
}
return(plot)
}signalplot(tests_rst = ERP_test,
multi = T,
data = ERPdata,
frames = 1:426,
datacol = 2:427 , subcol = 430 , chcol = 1 , othvarcol = c(428:429,431:432),
cpvarcol = 429,
sigp_yloc = -7,
sigpcol = "red",
sigpshape = 8,
wantbootplot=T, # do you want the bootinterval
fun=mean,
bootnum=10,
bootintval=c(.025,.975),
alpha=0.3)+
ylim(-10,5)
ERP only
1. scalpplot
2. topograph
1. scalpplot
revised from erpR package
plot the scalp :
BUT I think the bootplot is a better choice !!!!!
# revised scalp in erpR package
scalpplot <- function(elect_erpdata,frames,erpcol,subcol=NULL,elecol,othvarcol=NULL,
cpvarcol=NULL,
startmsec=frames[1],endmsec=frames[length(frames)],
layout=1,ylim=10, legend=TRUE, legend.lab="default",...){
scalp <- function (categ, smo = NULL, layout = 1, ylims = "auto", yrev = TRUE,
startmsec = -200, endmsec = 1200, lwd = 1, lty = 1,
color.list = c("blue", "red", "darkgreen"), legend = F, legend.lab = "default",
t.axis = seq(-100, endmsec, 200), scalp.array = NULL) {
if (class(categ) != "list") {
stop("input object must be a list!")
}
if (length(legend.lab) == 1 & legend.lab[1] == "default") {
legend.lab <- names(categ)
#legend.lab = deparse(substitute(categ))
#legend.lab = gsub("\\(", "", legend.lab)
#legend.lab = gsub("\\)", "", legend.lab)
#legend.lab = gsub("^list", "", legend.lab)
#legend.lab = gsub(" ", "", legend.lab)
#legend.lab = unlist(strsplit(legend.lab, ","))
}
if (length(lwd) == 1) {
lwd = rep(lwd, length(categ))
}
if (length(lty) == 1) {
lty = rep(lty, length(categ))
}
if (layout[1] == 1) {
electrodes = c("yaxis", "Fp1", "blank", "Fp2", "legend",
"F7", "F3", "FZ", "F4", "F8", "FT7", "FC3", "FCZ",
"FC4", "FT8", "T3", "C3", "CZ", "C4", "T4", "TP7",
"CP3", "CPZ", "CP4", "TP8", "T5", "P3", "PZ", "P4",
"T6", "xaxis", "O1", "OZ", "O2", "blank")
}
if (layout[1] == 2) {
electrodes = c("yaxis", "Fp1", "FPZ", "Fp2", "legend",
"F7", "F3", "FZ", "F4", "F8", "FT7", "FC3", "FCZ",
"FC4", "FT8", "T7", "C3", "CZ", "C4", "T8", "TP7",
"CP3", "CPZ", "CP4", "TP8", "P7", "P3", "PZ", "P4",
"P8", "xaxis", "O1", "OZ", "O2", "blank")
}
if (layout[1] == 3) {
electrodes = c("yaxis", "Fp1", "Fpz", "Fp2", "legend",
"F7", "F3", "FZ", "F4", "F8", "FT7", "FC3", "FCz",
"FC4", "FT8", "T3", "C3", "Cz", "C4", "T4", "TP7",
"CP3", "CPz", "CP4", "TP8", "T5", "P3", "PZ", "P4",
"T6", "xaxis", "O1", "blank", "O2", "blank")
}
if (layout[1] == 4) {
electrodes = c("yaxis", "Fp1", "blank", "Fp2", "legend",
"blank", "AF3", "blank", "AF4", "blank", "F7", "F3",
"Fz", "F4", "F8", "FC5", "FC1", "FCz", "FC2", "FC6",
"T7", "C3", "Cz", "C4", "T8", "blank", "CP1", "CPz",
"CP2", "blank", "P7", "P3", "Pz", "P4", "P8", "blank",
"O1", "blank", "O2", "blank")
}
if (layout[1] == 5) {
electrodes = c("yaxis", "Fp1", "Fpz", "Fp2", "legend",
"blank", "AF3", "blank", "AF4", "blank", "F7", "F3",
"Fz", "F4", "F8", "FC5", "FC1", "blank", "FC2", "FC6",
"T7", "C3", "Cz", "C4", "T8", "CP5", "CP1", "blank",
"CP2", "CP6", "P7", "P3", "Pz", "P4", "P8", "PO7",
"PO3", "POz", "PO4", "PO8", "blank", "O1", "Oz",
"O2", "blank")
}
if (length(layout) > 1) {
electrodes = layout
}
if (ylims == "auto") {
alldata = NULL
for (i in 1:length(categ)) {
alldata = rbind(alldata, categ[[i]])
}
ymax = max(alldata)
ymin = min(alldata)
yedge = max(c(ymax, abs(ymin)))
yedge = c(-yedge, yedge)
}
if (ylims != "auto") {
yedge = ylims
yedge = c(-ylims, ylims)
}
if (yrev == TRUE) {
yedge = sort(yedge, decreasing = T)
}
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(7, 5), mai = c(0, 0, 0, 0))
if (layout[1] == 5) {
par(mfrow = c(10, 5), mai = c(0, 0, 0, 0))
}
if (layout[1] == 4) {
par(mfrow = c(8, 5), mai = c(0, 0, 0, 0))
}
if (!is.null(scalp.array)) {
par(mfrow = scalp.array, mai = c(0, 0, 0, 0))
}
for (i in 1:(length(electrodes))) {
if (electrodes[i] == "yaxis") {
plot(1, type = "n", frame.plot = FALSE,
xlim = c(1,dim(categ[[1]])[1]),
xaxt = "n", yaxt = "n",
ylim = c(yedge[1] + yedge[1]/3, yedge[2] + (yedge[2]/3)))
axis(side = 2, pos = dim(categ[[1]])[1]/2,
at = c(round(ceiling(yedge[1]),0),
round(ceiling(yedge[1])/2, 0), 0,
round(floor(yedge[2])/2,0), round(floor(yedge[2]), 0)),
cex.axis = 0.8,
las = 2)
text((dim(categ[[1]])[1]/2) + (dim(categ[[1]])[1]/8),
0, labels = expression(paste(mu, "V")), cex = 1.4)
}
if (electrodes[i] == "blank") {
plot.new()
}
if (electrodes[i] == "legend") {
plot.new()
if (legend == "TRUE") {
legend("center", legend = legend.lab, col = color.list,
cex = 1.2, lty = lty, lwd = lwd)
}
}
if (electrodes[i] == "xaxis") {
plot(1, type = "n", frame.plot = FALSE, xlim = c(1,
dim(categ[[1]])[1]), xaxt = "n", yaxt = "n",
ylim = c(yedge[1] + yedge[1]/3, yedge[2] + (yedge[2]/3)))
axis(1, pos = 0, at = msectopoints(t.axis, dim(categ[[1]])[1],
startmsec, endmsec),
labels = paste(t.axis))
}
if (!electrodes[i] %in% c("xaxis", "yaxis", "legend",
"blank")) {
if (!is.null(smo)) {
el=smooth.spline(categ[[1]][[electrodes[i]]][1:dim(categ[[1]])[1]],
spar = smo)
}
else {
el = categ[[1]][[electrodes[i]]][1:dim(categ[[1]])[1]]
}
plot(el, type = "l", ylim = c(yedge[1] + yedge[1]/3,
yedge[2] + (yedge[2]/3)), col = color.list[1],
main = "", ylab = "", xlab = "", cex.main = 0.85,
xlim = c(1, dim(categ[[1]])[1]), xaxt = "n",
yaxt = "n", frame.plot = FALSE, lwd = lwd[1],
lty = lty[1])
totalendmsec = endmsec + abs(startmsec)
zeropoint = msectopoints(0, dim(categ[[1]])[1], startmsec = startmsec,
endmsec = endmsec)
segments(x0 = zeropoint, y0 = -0.8, x1 = zeropoint,
y1 = 0.5, lwd = 1.5)
abline(h = 0, lty = "longdash")
mtext(electrodes[i], side = 3, line = -2)
if (length(categ) > 1 & electrodes[i] != "blank") {
for (k in 2:length(categ)) {
if (!is.null(smo)) {
el = smooth.spline(categ[[k]][[electrodes[i]]][1:dim(categ[[1]])[1]],
spar = smo)
}
else {
el = categ[[k]][[electrodes[i]]][1:dim(categ[[1]])[1]]
}
lines(el, col = color.list[k], lwd = lwd[k],
lty = lty[k])
}
}
}
}
par(oldpar)
}
data <- elect_erpdata
if (is.null(cpvarcol)==TRUE){
data_list <- list()
data_list[[1]] <- data
rownames(data_list[[1]]) <- data_list[[1]][,elecol]
data_list[[1]] <- data_list[[1]][,erpcol]
data_list[[1]] <- data.frame(t(data_list[[1]]))
names(data_list) <- "ERP Signal"
} else {
if ( class(data[,cpvarcol]) != "factor"){
data[,cpvarcol] <- as.factor(data[,cpvarcol])
}
factor_num <- length(levels(data[,cpvarcol]))
data_list <- list()
for ( i in 1 : factor_num ) {
data_list[[i]] <- subset(data,
data[,cpvarcol]==(levels(data[,cpvarcol])[i]))
}
names(data_list) <- levels(data[,cpvarcol])
for (i in 1:factor_num){
rownames(data_list[[i]]) <- data_list[[i]][,elecol]
data_list[[i]] <- data_list[[i]][,erpcol]
data_list[[i]] <- data.frame(t(data_list[[i]]))
}
}
scalp(data_list,layout=layout,ylim=ylim,
startmsec=startmsec,endmsec=endmsec,
legend=legend,legend.lab=legend.lab,...)
}ERPdata_agg <- aggraw(ERPdata,frames=1:426,2:427,430,1,othvarcol=c(428,429,431,432),
aggvarcol=c(428,429,1),
fun=mean)
scalpplot(ERPdata_agg,frames=1:426,erpcol=4:429,subcol=NULL,elecol=3,othvarcol = 1,
cpvarcol=2)
2. topograph
revised from erpR
but I have not yet really understand the meaning of their setting.
topograph <- function(elect_erpdata,frames,chcol, # only
startmsec=0, endmsec=1000, win.ini=0,win.end=1000,...){
# need some change
data <- elect_erpdata
rownames(data) <- data[,chcol]
data <- data[,-chcol]
topo_data <- data.frame(t(data))
notfound <- topoplot(topo_data, return.notfound=TRUE)
topoplot(topo_data,startmsec,endmsec,win.ini,win.end,exclude=notfound,...)
}ERP_word <- subset(ERPdata_agg,ERPdata_agg$Condition=="word")
topograph(ERP_word[,-c(1:2)], # you should unput ONLY the erpdata with electrode
frames = 1:426,
chcol = 1,
startmsec=1, endmsec=426,
win.ini=1,win.end=426,
draw.nose=T,projection="equalarea")WARNING: your data (after interpolation) are out of range as compared to the zlims specified.
Your data range is: -8.25, 5.63
Your zlims are: -8, 8 