ERP package
Outline
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).
preprocessing
1. list2df
2. AddExpCondSub
3. ind_select
4. ele_select
5. conj_select
6. aggraw
EDA
1.edaplot
2.bootplot
tests
1.testplot
2.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(simerp);data(erpcz);data(ERPsets)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){
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){
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:2,427:430)],3) # look at the data Channel value.1 value.426 Experiment Condition Subject
1 Fp1 0.36995740 9.971283 Exp1 word subj1
2 Fp2 -0.02759907 9.678596 Exp1 word subj1
3 F3 0.42971300 4.413941 Exp1 word subj13. 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)
}
S <- ind_select(data = ERP_df,
frames = 1:426,
datacol = 2:427,
subcol = 430,
chcol = 1,
othvarcol = 428:429,
select_sub=c("subj1","subj2"))
dim(S) # 2 conditions *34 electrodes * 2 subjects [1] 136 4304. 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 = ERP_df,
frames = 1:426,
datacol = 2:427,
subcol = 430,
chcol = 1,
othvarcol = 428:429,
select_ch=c("Fp1","Fp2"))
dim(Fp) # 2 conditions * 20 subjects * 2 electrode[1] 80 4305. conj_select
select channels and subjects
conj_select <- function(data,frames,datacol,subcol,chcol,othvarcol=NULL,
select_sub,
select_ch,...){
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 = ERP_df,
frames = 1:426,
datacol = 2:427,
subcol = 430,
chcol = 1,
othvarcol = 428:429,
select_sub=c("subj1"),
select_ch=c("Fp1","Fp2"))
dim(SE) # 1 sub * 2 Channel * 2 Condition[1] 4 4306. 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 = ERP_df,
frames = 1:426,
datacol = 2:427, #you could put all elements you want to aggregate in datacol argument
subcol = 430,
chcol = 1,
othvarcol = 428:429,
aggvarcol= c(1,429),
fun=median) # mean median and so on
tail(A[,1:6]) 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.68893100EDA
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(ERP_df,
frames = 1:426,
datacol=2:427,
subcol=430,
chcol=1,
othvarcol=428:429,
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
#head(erpcz)
erpcz$Score <- rep(round(rnorm(20,50,5)),each=2) # for demostration
edaplot(erpcz,
frames=seq(0,1001,4),
erpcol = 1:251,subcol=252,chcol=NULL,othvarcol=c(253,254))+
geom_line(aes(col=Score))+ # you could specify the color argument
facet_wrap(~Subject)+
theme_bw()+
theme(legend.position="none")+
xlim(-100,1100)+
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.
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])
}
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
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)
}
#head(ERP_df)
#multiple electrode
bootplot(ERP_df, # input the data
frames=1:426,
datacol=2:427,
subcol=430,
chcol=1,
othvarcol=428,
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"))
# single channel
erpcz$Eletrode <- "Cz" # create channel variable
bootplot(erpcz,
frames=1:251,
datacol = 1:251,
subcol=252,
chcol=255,
othvarcol=254,
cpvarcol=253,
fun=median,
bootnum=300,
bootintval=c(.025,.975),
alpha=0.5)
tests
1.testplot
2.signalplot
1. testplot
plot the difference curve and highlight the significant value
testplot <- function(tests_rst,frames,
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.
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)
}
return(plot)
}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)
testplot(tests, # input the test results
frames,
cor=F, # make sure you model are not in a correlation setting
siglinecol = "blue", # indicate the color of the significant line
sigp=T, #Significant points are added.
sigpcol = "yellow", #The color of Significant points (when sigp=T)
sigp_yloc = -10)+ #The y location of significant points (when sigp=T)
xlab("Time (ms)")+
ylab("Difference ERP curves")+
theme_bw()
tests = erpfatest(simerp[,1:251],design=model.matrix(~y,data=simerp),
nbf=3,min.err=1e-01,maxiter=10)
frames = seq(0,1001,4)
testplot(tests,frames,cor=T, #Correlation should be TRUE (cor=T).
sigp=T,sigpcol="yellow",19,sigp_yloc=-1)+
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)
# there will be some warning messages if you cover the original color value.
ylim(-2,2)+
xlim(-100,1100)+
theme_bw()+
theme(legend.position="none")+
labs(list(x="Time (ms)",y="Correlation",title="Simulation"))
##### 2. signalplot
##### plot original signal curve and highlight the significant location
##### and it is possible to add some bootstrap interval
signalplot <- function(tests_rst,
data,
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) {
dta <- data
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)
return(plot)
}signalplot(tests,erpcz,
datacol=1:251,subcol=252,chcol=255,othvarcol=254,
cpvarcol=253,
sigp_yloc= -10,
sigpcol="pink",
sigpshape=19,
wantbootplot=T,
fun=mean,
bootnum=10,
bootintval=c(.025,.975),
alpha=0.3)
####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
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)
}
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",...){
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"){
stop("cpvar should be factor")
}
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,...)
}ERP_df_agg <- aggraw(ERP_df,frames=1:426,2:427,430,1,othvarcol=c(428,429),
aggvarcol=c(428,429,1),
fun=mean)
#head(ERP_df_agg)
ERP_df_agg$Condition<-as.factor(ERP_df_agg$Condition)
scalpplot(ERP_df_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,...)
}##Example
ERP_df_agg <- aggraw(ERP_df,frames=1:426,2:427,430,1,othvarcol=c(428,429),
aggvarcol=c(428,429,1),
fun=mean)
ERP_word <- subset(ERP_df_agg,ERP_df_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 