Analyzing EEG data using GAMs

• Introduction
  • ERPs to study gender violations
  • Research question
• Design
• Methods: R code
• Results
• Discussion

• A P600 (a positivity 'around' 600 ms. after stimulus onset) is sensitive to grammatical violations
• An N400 (a negativity 'around' 400 ms. after stimulus onset) is modulated by semantic context and lexical properties of a word
  • The P600/N400 are found by comparing the incorrect sentences to the correct sentences
• Native speakers appear to show a P600 for gender violations
  • But analyzed by averaging over items and over subjects!
• Here we are interested in how non-native speakers respond to gender violations (joint work with Nienke Meulman)
• Gender is very hard to learn for L2 learners
• Even though behaviorally L2 learners might show correct responses, the brain may reveal differences in processing gender

• Is the P600 for gender violations dependent on age of arrival for the L2 learners of German?

• Today: analysis of single region of interest (ROI 8)

• 67 L2 speakers of German (Slavic L1)
• Auditory presentation of correct sentences or sentences with a gender violation (incorrect determiner; no determiners in L1)
• 48 items in each condition: 96 trials per participant (minus artifacts)
• Example:
  Nach der Schlägerei ist das/*der Auge des Angestellten von der Krankenschwester versorgt worden.
  After the fight theneut/*themasc eye of the worker was treated by the nurse

(the signal has been downsampled to 100 Hz)

(note the arbitrary age splits, however)

(for all participants; to limit CPU time, the time window was set to [500,1300])

dat$SubjectCor = interaction(dat$Subject, dat$Correctness)
dat$WordCor = interaction(dat$Word, dat$Correctness)
# duration discrete=F: 9300 seconds, so using discrete with nthreads is much faster
system.time(m0 <- bam(uV ~ s(Time, by = Correctness) + Correctness + s(Time, SubjectCor, bs = "fs", m = 1) + 
    s(Time, WordCor, bs = "fs", m = 1), data = dat, rho = rhoval, AR.start = dat$start.event, discrete = T, 
    nthreads = 2))

#    user  system elapsed 
#  853.89    4.82  449.54

system.time(smry0 <- summary(m0))  # storing summary as it also takes some time to compute

#    user  system elapsed 
#    77.1     0.6    77.8

par(mfrow = c(1, 2))
plot_smooth(m0, view = "Time", rug = F, plot_all = "Correctness", main = "", rm.ranef = T, eegAxis = T)
plot_diff(m0, view = "Time", comp = list(Correctness = c("incor", "cor")), rm.ranef = T, eegAxis = T)

par(mfrow = c(1, 2))
plot(m0b, select = 2, shade = T, rug = F, main = "Binary difference curve", ylim = c(3, -3))
plot(m0o, select = 2, shade = T, rug = F, main = "Ordered factor difference curve", ylim = c(3, -3))

(te is used to model a non-linear interaction)

m1 = bam(uV ~ te(Time, AoArr, by = Correctness) + Correctness + s(Time, SubjectCor, bs = "fs", m = 1) + 
    s(Time, WordCor, bs = "fs", m = 1), data = dat, rho = rhoval, AR.start = dat$start.event, discrete = T, 
    nthreads = 2)
smry1 = summary(m1)
smry1$p.table

#                  Estimate Std. Error t value Pr(>|t|)
# (Intercept)        -0.408      0.474  -0.861    0.389
# Correctnessincor    0.444      0.672   0.661    0.509

smry1$s.table

#                                    edf  Ref.df     F   p-value
# te(Time,AoArr):Correctnesscor     3.10    3.20 1.764  1.52e-01
# te(Time,AoArr):Correctnessincor   5.88    6.95 4.917  1.57e-05
# s(Time,SubjectCor)              111.94 1203.00 0.658 2.24e-119
# s(Time,WordCor)                 134.38 1727.00 0.275  5.93e-55

Note the default maximum number of edf's per smooth: 24 (5\(^2\) - 1)

plot_diff2(m1, view = c("Time", "AoArr"), rm.ranef = T, comp = list(Correctness = c("incor", "cor")))
fadeRug(dat$Time, dat$AoArr)  # indicate points without data

(note that the \(y\)-axes are not inverted here)

• Go to: http://eolomea.let.rug.nl/GAM/InterpretingInteractions

m2 = bam(uV ~ s(Time, by = Correctness) + s(AoArr, by = Correctness) + ti(Time, AoArr, by = Correctness) + 
    Correctness + s(Time, SubjectCor, bs = "fs", m = 1) + s(Time, WordCor, bs = "fs", m = 1), data = dat, 
    rho = rhoval, AR.start = dat$start.event, discrete = T, nthreads = 2)
smry2 = summary(m2)
smry2$p.table

#                  Estimate Std. Error t value Pr(>|t|)
# (Intercept)        -0.418      0.475  -0.881    0.378
# Correctnessincor    0.480      0.673   0.714    0.475

smry2$s.table

#                                    edf  Ref.df      F   p-value
# s(Time):Correctnesscor            1.04    1.08 0.0682  8.37e-01
# s(Time):Correctnessincor          3.27    4.25 6.2711  3.54e-05
# s(AoArr):Correctnesscor           1.06    1.07 2.6019  1.12e-01
# s(AoArr):Correctnessincor         1.06    1.06 3.8290  4.32e-02
# ti(Time,AoArr):Correctnesscor     1.04    1.07 2.1989  1.27e-01
# ti(Time,AoArr):Correctnessincor   1.53    1.89 0.3922  6.17e-01
# s(Time,SubjectCor)              111.82 1204.00 0.6515 1.32e-118
# s(Time,WordCor)                 134.38 1727.00 0.2755  5.95e-55

m3 = bam(uV ~ s(Time, by = Correctness) + s(AoArr, by = Correctness) + Correctness + s(Time, SubjectCor, 
    bs = "fs", m = 1) + s(Time, WordCor, bs = "fs", m = 1), data = dat, rho = rhoval, AR.start = dat$start.event, 
    discrete = T, nthreads = 2)
smry3 = summary(m3)
smry3$p.table

#                  Estimate Std. Error t value Pr(>|t|)
# (Intercept)        -0.414      0.475  -0.872    0.383
# Correctnessincor    0.486      0.673   0.722    0.470

smry3$s.table

#                              edf  Ref.df     F   p-value
# s(Time):Correctnesscor      1.04    1.08 0.385  5.74e-01
# s(Time):Correctnessincor    3.32    4.32 6.765  1.25e-05
# s(AoArr):Correctnesscor     1.07    1.08 2.690  1.06e-01
# s(AoArr):Correctnessincor   1.10    1.11 3.386  5.21e-02
# s(Time,SubjectCor)        111.78 1204.00 0.650 2.32e-118
# s(Time,WordCor)           134.38 1727.00 0.275  5.97e-55

par(mfrow = c(1, 2))
plot_diff(m3, view = "Time", comp = list(Correctness = c("incor", "cor")), rm.ranef = T, eegAxis = T)
plot_diff(m3, view = "AoArr", comp = list(Correctness = c("incor", "cor")), rm.ranef = T, eegAxis = T)

library(car)
par(mfrow = c(1, 2))
qqp(resid_gam(m3))
hist(resid_gam(m3))

• This type of residual distribution is common for EEG data
• These deviations are problematic and may affect \(p\)-values
• Solution: model fitting using scaled-\(t\) distribution
  • Implemented in gam: family="scat"
  • Adding rho to gam possible via custom function (created by Natalya Pya)
  • Unfortunately, gam is too slow, so a bam approach is in development
  • Current state: implemented for fREML estimation with pre-specified scaled-\(t\) parameters (\(\theta\)) and no discrete=T (see Meulman, Wieling, Sprenger, Stowe and Schmid, PLOS ONE, 2015)
  • Thus gam still necessary to determine (\(\theta\)) (but applied to subset / simpler model specification)

1. Fit an intercept-only scaled-t gam model (including rho, method="REML") to all data and extract \(\theta\)
2. Fit the complete model without random effects as a scaled-\(t\) gam model (including rho, method="REML") to all data and extract \(\theta\)
3. Compare the models' \(\theta\) values:
   • Relatively small differences: use \(\theta\) of the complex model in the bam estimation
   • Relatively large differences: determine \(\theta\) by fittting the complete model specification (including random effects) to a small subset of the data (randomly select trials, compare \(\theta\) across a few runs)

library(car)
par(mfrow = c(1, 2))
qqp(resid_gam(m3), main = "m3")
qqp(resid(m4.scat), main = "m4.scat")  # resid of scaled-t model takes rho into account

• Still work to do: e.g., testing the significance of other possibly important variables (proficiency, age, etc.)
• But don't make it too complex!
  • There is much variation present in EEG data and adding very complex surfaces will almost certainly improve your model significantly
  • Keep it simple, otherwise you won't be able to compute/interpret the results
  • Also note that for hypothesis testing, it is essential not to conduct step-wise approaches
• To conclude:
  • GAMs are very useful to analyze EEG and other time-series data
  • The method is very suitable to detect non-linear patterns, while taking into account individual variation and correcting for autocorrelation

• We have applied GAMs to EEG data and learned how to:
  • model difference smooths directly using binary predictors and ordered factors
  • use te(Time,AoArr) to model a non-linear interaction
  • decompose te(Time,AoArr) using ti() and s()
  • use a scaled-\(t\) model to improve residuals
    • (fitting this type of model for a complete data set may take a few weeks of CPU time, however)
• While we have analyzed a single region of interest, GAMs allow for spatial distribution analyses
  • E.g., via te(x,y,Time)
  • See http://eolomea.let.rug.nl/GAM/ScalpEEG
• After the break:
  • http://www.let.rug.nl/wieling/statscourse/lecture4/lab
• Next lecture: generalized additive logistic modeling using dialect data