Preprocessing of ERP Data (Both averaged within-subject and single-trial data)

The first part of this script just takes the matlab (EEGLAB) output with subject-level averages for tot/no-tot conditions for each electrode 250-700ms and saves it as a .rda object in the data/ folder for analysis.

erpSubjectAveraged <- read.table("data-raw/ERPLAB_AverageAmpERP.txt", header = T)
devtools::use_data(erpSubjectAveraged)

This rest of this script takes inputs of a text file of all the single trial master EEG data (output from EEGLAB) and combines it with the behavioral data, converting it to long format such that each row is 1 timepoint, at 1 electrode, during 1 trial. Note, this script will not run with the files currently on github, please email pab2163@columbia.edu for access to large raw data files

require(lme4)
require(tidyverse)
require(grid)
require(gridExtra)
require(dplyr)


# loading file ------------------------------------------------------------

# FOR AREA AND AMPLITUDES
#NOTE -- this file was not uploaded because of size! Please email pab2163@columbia.edu if you would like access 
filename = "data-raw/singleTrials.txt"

d <- read.delim(filename, header = F) %>%
  as_tibble() %>%
  select(-V4) %>%
  mutate(V5 = ifelse(V5 == 1, "yes", "no"))
names(d)[1:5] <- c("id","epoch","chan","tot","rejects")
mean(d$rejects)
View(d)

Add channels from text file

chan <- read.csv("data-raw/channels.txt", header=T) %>%
  as_tibble()
names(chan)[2] <- "chan"

d <- left_join(d, chan) %>%
  select(-chan)

Add behavioral data file

b <- read.csv("data-raw/behavioral_data_cleaned.csv", header=T) %>%
  as_tibble() %>%
  select(-question_num, subj_answer, -question, -correct_answer, -subj_answer, -test_trialnum) %>%
  filter(tot != "N/A") %>%
  mutate(recall = ifelse(recall == "I", 0, 1))

View(b)
unique(b$id)

# remove excluded subjects
b <- subset(b, id!= 2)
b <- subset(b, id != 5)
b <- subset(b, id != 19)

# comparing # of epochs ---------------------------------------------------

b1 <- b %>%
  group_by(id, tot) %>%
  summarise(n=n())
View(b1)

d1 <- d %>%
  filter(chlabel == "Pz") %>%    # only 1 elec for counting otherwise repeats
  group_by(id, tot) %>%
  summarise(n=n())

View(d1)
################## DO NOT PROCEED UNLESS THEY MATCH UP!!! ################## 
b1$n == d1$n

Sort by id and check again

# sorting by  id to line up -----------------------------------------------
b1 <- b[order(b$id, b$trialnum),] %>%
  group_by(id) %>%
  mutate(epoch = seq_len(n()))
d1 <- d %>%
  filter(chlabel == "Pz")

################## DO NOT PROCEED UNLESS THEY MATCH UP!!! ################## 
#b1$tot == d1$tot
length((b1$tot == d1$tot) == TRUE) == dim(b1)[1] #for summary

Join behavior and ERP frames

b1 <- select(b1, -trialnum)

#removing reject files
beh_erp <- left_join(d, b1, by = c("id", "tot", "epoch")) %>%
  filter(rejects == 0) %>%
  select(-rejects)
names(beh_erp)

#remove mistake files
beh_erp <- filter(beh_erp, mistake == 0)

More Artifact Rejection - Remove any trials for which amplitude was more than 50uV from baseline at any timepoint

# make a data frame with only eeg amplitude values, no other columns
beh_erp_nums <- beh_erp[-c(1:3, 454:458)]
names(beh_erp_nums)

# calculate max and min amplitude for each trial, for each electrode
beh_erp$maxAmp <- apply(beh_erp_nums,1, max)
beh_erp$minAmp <- apply(beh_erp_nums,1, min)

# if for an electrode on a trial, the absolute value of the amplitude is above 50, remove that electrode for that trial
beh_erp_clip <- filter(beh_erp, maxAmp < 50 & minAmp > -50)

# Check to make sure filtering is working
numCheck <- beh_erp_clip[-c(1:3, 454:460)]
names(numCheck)

beh_erp_clip$maxAmp <- apply(numCheck,1, max)
beh_erp_clip$minAmp <- apply(numCheck,1, min)

# Should all be within 50
summary(beh_erp_clip$minAmp)
summary(beh_erp_clip$maxAmp)

Put data in long format

long_clip <- beh_erp_clip %>%
  gather(time, uv, V7:V456)

names(long_clip)

# Get Accurate Timestamps -------------------------------------------------

# Function for turning the current column names into timestamps
mysubfunction <- function(x){
  mynewx1 <- gsub("V", "", x)
  mynewx1 <- as.numeric(mynewx1)
  mynewx2 <- ((mynewx1 -7)*2)-200 
  return(mynewx2)
}


#Apply that function via sapply
long_clip$timestamp <- sapply(long_clip$time, mysubfunction)
unique(long_clip$timestamp)

# SAVE!! (not run here) -------------------------------------------------
# save(long_clip, file = '../data/long_clip.rda', compress = 'xz)


pab2163/TOT_ERP documentation built on May 16, 2019, 8:13 p.m.