ERP analysis (example)¶
Running a model¶
An example of a model would look something like this:
# Load the mgcv library
suppressMessages(library(mgcv))
# Prevent the code below from running
run = FALSE
# Code to run a model
if(run) {
# Load the data
load("eegdata.rda")
# Create a few new columns
eegdata$Start = erps$Time == 0
eegdata$signal = "O1"
# Restrict temporal range
eegdata = eegdata[eegdata$Time >= 0 & eegdata$Time <= 600,]
# Run the model
model.gam = bam(signal ~ s(Word, bs = "re") + s(Text, bs = "re") + s(Phrase, bs = "re") + s(Time, k = 20) +
s(Trial, Subject, bs = "fs", m = 1) + s(Time, Subject, bs = "fs", m = 1) +
s(LogComplexity) + ti(Time, LogComplexity, k = c(20, 5)) +
s(logFrequencyWord) + ti(Time, logFrequencyWord, k = c(20, 5)) +
s(logFrequencyPhrase) + ti(Time, logFrequencyPhrase, k = c(20, 5)) +
s(logFrequencyPrep) + ti(Time, logFrequencyPrep, k = c(20, 5)) +
s(re) + ti(Time, re, k = c(20, 5)), data = eegdata,
rho = 0.75, AR.start = eegdata$Start)
}
As pointed out by Harald, you could furthermore include by-item factor smooths for Time, of the form s(Time, Item, bs="fs", m=1) if you so desire. Here, we decided not to do that, because the resulting model seemed overspecified upon inspection.
The argument AR.start takes as input a variable in your dataframe that indicates whether the current data point is the start of a new time series (in this case: a new trial).
Rather than running the model from scratch, we simply load the results here:
# Load the model
load("model.gam.rda")
The summary of this model looks like this:
# Show the summary of the model
summary(model.gam)
Plotting results¶
For each predictor, the model includes a main effect smooth, as well as a tinsor product interaction of Time by predictor. As an example, we plot the main effect smooth, as well as the tinsor product interaction for the predictor word frequency (logFrequencyWord) below.
Plot for the main effect smooth:
# Set plot parameters
options(repr.plot.width=5, repr.plot.height=2.8)
options(jupyter.plot_mimetypes = "image/png")
par(mar=c(4.5,4.5,3.5,0.5))
# Plot
plot.gam(model.gam,select=9,rug=FALSE,ylim=c(-2,2),xlab="Word Frequency",
ylab=expression(paste("Voltage (", mu, "V)", sep="")),
main="Main effect of Word Frequency (time domain)",
cex.lab="0.7",cex.axis="0.6",cex.main="0.75")
Plot for tinsor product interaction with Time:
# Set plot parameters
options(repr.plot.width=5, repr.plot.height=2.8)
options(jupyter.plot_mimetypes = "image/png")
par(mar=c(4.5,4.5,3.5,0.5))
# Source helpers.R to use an edited version of pvisgam from the itsadug package
source("helpers.R")
# Plot
pvisgam(model.gam, view = c("Time", "logFrequencyWord"), color="topo", plot.type="contour",
zlim = c(-1,1),xlab="time (ms)", ylab="Word Frequency",
main="Time by Word Frequency interaction (time domain)",
cex.lab="0.7",cex.axis="0.6",cex.main="0.75")
Timing information¶
Estimates of the temporal onset of an effect are calculated on the basis of predicted values for the relevant model term. Using an edited version of pvisgam (see helpers.R), we generate a matrix with predicted values, as well as the standard errors for each cell in this matrix.
We return the results by setting returnMatOnly to TRUE.
To be able to get accurate estimates, we need a high resolution. We therefore set n.grid to 600 (i.e., one predicted value per ms).
res = pvisgam(model.gam, view = c("Time", "logFrequencyWord"), plot.type="contour", n.grid=600,
returnMatOnly=TRUE)
The temporal onset of an effect is defined as the first point in time at which 0 is within 3 standard errors of the predicted value for a given value of a predictor. To make our life a bit easier, we define a function that finds this point in time for us:
# Define function that retrieves the temporal onset of an effect
getTiming.fnc = function(res, pred_range = NA, se = 3,
pred_vals = seq(13.60304, 18.34987, length.out = 600)) {
# Generate a matrix with minimum values
mat_min = res$fit - se * res$se
# Generate a matrix with maximum values
mat_max = res$fit + se * res$se
# Limit predictor values for which the temporal onset is retrieved
cols = 1:600
if(!is.na(pred_range[1])) {
cols = which(pred_vals>=pred_range[1] & pred_vals<=pred_range[2])
}
# Find time points at which a significant effect is observed
sign = sapply(1:600,FUN=function(x) { max(mat_min[x,cols] > 0) | min(mat_max[x,cols]) < 0 })
# Define the onset of an effect as the first point in time at which a significant effect is
# observed
onset = which(sign)[1]
# Return onset
return(onset)
}
The function getTiming.fnc() can now be used to establish the temporal onset of the effect of word frequency:
# Get the overall onset of the effect of word frequency
getTiming.fnc(res)
Or, if we want the onset of the effect of word frequency for a subset of the word frequency range:
# Get the onset of the effect of word frequency for low frequency words
getTiming.fnc(res,pred_range = c(13,15))
Frequency domain¶
To get information from the frequency domain, we did not carry out a separate analysis in the frequency domain. Instead, we took the results of a model fit to the data in the time domain and converted these results to the frequency domain. I would therefore refer to the procedure we followed as a " description of the results in the frequency domain", rather than a "frequency domain analysis of the data".
Below, I give an example of the conversion of results to the frequency domain for the Time by logWordFrequency interaction introduced above.
First, we generate a matrix with predicted values in the time domain:
# Create data matrix in the time domain
mat = pvisgam(model.gam, view = c("Time", "logFrequencyWord"), plot.type="contour", n.grid=600,
returnMatOnly=TRUE)$fit
# Set row names
rownames(mat) = 1:600
# Set column names of the matrix to match the range of logFrequencyWord
colnames(mat) = seq(13.60304, 18.34987, length.out = 600)
Next, we convert the matrix to a time-series object:
# Create a time series object
timeseries = ts(mat,start=0,end=1,frequency=1000)
We calculate the spectral density of the time-series with the spectrum() function. The spectral density provides a description of the time-series in the frequency domain:
# Calculate spectral density
spectrum = spectrum(timeseries,method="pgram",plot=FALSE)
Based on the spectrum, we create another matrix. This time, the matrix contains the results in the frequency domain.
# Create empty matrix
specmat = matrix(nrow=600,ncol=30,0)
# Enter data into the matrix for each point in time (0-600 ms) and
# for each of the 30 lowest frequencies in the spectrum (~ 0-29 Hz)
for(i in 1:600) {
specmat[i,] = spectrum$spec[,i][1:30]
}
# Set appropriate column names
colnames(specmat) = spectrum$freq[1:30]
We can now plot the representation of the results in the frequency domain, using an edited version of the filled.contour function.
# Set plot parameters
options(repr.plot.width=5, repr.plot.height=2.8)
options(jupyter.plot_mimetypes = "image/png")
par(mar=c(4.5,4.5,3.5,0.5))
# Plot
filled.contour.flex(x = seq(13.60304, 18.34987, length.out = 600), y = spectrum$freq[1:30],
z = specmat, ylim=c(spectrum$freq[1],10), color.palette=topo.colors,
xlab="Word Frequency",ylab="frequency (Hz)",
main="Time by Word Frequency interaction (frequency domain)",
cex.lab="0.7",cex.axis="0.6",cex.main="0.75")
Ignoring the blip at the lower right hand side of the plot, maximum spectral density is reached around 3 Hz.