Bayesian Modeling of Event-Related Potentials 🧠

Dr. Cheng-Han Yu
Department of Mathematical and Statistical Sciences
Marquette University

2024-05-03

Why Bayesian and Event-Related Potentials (ERP)?

• Why Bayesian?
  • Uncertainty quantification
  • Unified hierarchical framework
  • Incorporate information not provided by the data

Source: Schoot et al. (2021)

• I don’t know how to do research using frequentist approaches. XD
• Quantify the uncertainty about the parameters easily and naturally once we obtain the posterior distributions.
• Also, it’s quite straightforward to build a hierarchical model using a Bayesian approach. And so for this project, we are doing the inference at the individual subject level and the population group level for different groups in an unified hierarchical framework.
• And as you all may know, Bayesian models can incorporate useful prior information into the model.

Why Bayesian and Event-Related Potentials (ERP)?

• Why Bayesian?
  • Uncertainty quantification
  • Unified hierarchical framework
  • Incorporate information not provided by the data
• Why ERP?
  • noninvasive
  • high temporal resolution
  • less costly
  • lack of sophisticated statistical methods

Source: https://www.eegget-it.nl/erp.html

• The second question is Why ERP?
• Well, ERP signals in fact come from EEG, which is an noninvasive neuroimaging modality with pretty high time resolution.
• It is relatively less costly comparing to other neuroimaging modalities, such as fMRI.
• More importantly, while there are lots of statistical research on EEG signals, not many sophistocated statistical methods out there for analyzing ERP data.

• Develop Bayesian hierarchical models for inference about the (1) underlying ERP waveform, (2) the latency and (3) amplitude of ERP components.

• So our goal is to develop Bayesian hierarchical models for inference about latency, amplitude as well as the underlying ERP waveform.
• We in particular focus on the latency of some ERP components because it’s related to cognitive processes, but there are not many good approaches that help us identify or quantify it.
• So let’s talk about ERP signals in more detail.

Electroencephalography (EEG)

Source: https://en.wikipedia.org/wiki/Electroencephalography

• Let’s first talk about EEG signals.
• Basically when a EEG experiment is conducted, a person’s brain activity is recorded from the electrodes placed on the person’s head.
• When you look at the raw EEG, it’s a pretty complicated signal. It’s basically the sum of everything the brain is doing.
• It can indicate if someone is awake or asleep, but in its raw form, the EEG cannot reveal the specific neural processes that underlie perception, cognition or emotion.

The raw EEG can indicate general state, e.g., awake or asleep, but it can’t tell us much about specific neural processes that underlie perception, cognition or emotion.

EEG to Event-Related Potentials (ERP)

Source: ERP Boot Camp (Luck and Kappenman, 2020)

• However, if we apply a simple set of signal processing operations to the raw EEG, we can pull out portions of the neural activity that are related to specific sensory, cognitive, affective and motor processes.
• So to do this, we basically record the EEG while a person is performing some task in response to some external stimuli, for example look at an image, picture, listen to some words, and a typical experiment would involve many stimuli and many responses.

Event-Related Potentials (ERP)

Source: Luck (2012)

• To extract the so-called event related potentials, we first grab the period of EEG following each stimulus
• Then we take all the EEG epochs, and line them up in time.
• Finally we simply average across the epoch from some event or stimuli responses to get the ERP which is an electrical potential or voltage generated by the brain that is related to some event.

What ERP Can Do

Source: ERP Boot Camp (Luck and Kappenman, 2020)

• So what can we do with ERPs?
• Read words…

Estimate ERP Components

• Interested in estimating the amplitude and latency of some ERP component.

Source: Rossion & Jacques (2012)

• In a typical ERP study, we focus on estimating the amplitude and latency of some ERP components or the important peaks or dips that are related to human sensory, cognitive or motor processes where amplitude is the magnitude of the component and the latency is the time when the ERP component occurs.

Inference Methods

Current methods

• Grand ERP by averaging across subjects: lost subject-level information.

• No uncertainty quantification.

• Filtering for removing drifts or noises: no rule of thumb for which frequency/bandwidth should be used.

• Pre-specify component search window for quantifying amplitude: no systematic ways of selecting it.

• Separate models for control-experimental comparison.

Source: Fisher et al. (2010)

• Current methods have several issues or drawbacks.
• In order to obtain precise point estimates of latency and amplitude, or the entire ERP waveform, people tend to get the grand ERP waveform by averaging across subjects.
• By doing so, subject-level information is lost.
• NO UQ
• To get a smooth ERP curve, some kinds of filters are necessary to remove signal drifts and noises. But no rule of thumb for which frequency/bandwidth should be used.
• To analyze ERP components, a time search window needs to be specified to frame the range of the components and an inappropriate time window may result in false analysis.
• But so far there are no systematic or consistent ways of selecting such windows, and the window is chosen based on previous scientific findings.

Solve these issues by

• Bayesian hierarchical model with a derivative Gaussian process prior: Yu et al., Biometrics (2022), Li et al., JASA T & M (2023), Yu et al., Data Science in Science (2024)

Gaussian Variables

Multivariate Gaussian

• \mathbf{Y}\in \mathbb{R}^d \sim N_d\left(\boldsymbol \mu, \boldsymbol \Sigma\right)

• Bivariate Gaussian (d=2):

\mathbf{Y}= \begin{pmatrix} Y_1 \\ Y_2 \end{pmatrix} \boldsymbol \mu= \begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix} \boldsymbol \Sigma= \begin{pmatrix} \sigma_1^2 & \rho_{12}\sigma_1\sigma_2\\ \rho_{12}\sigma_1\sigma_2 & \sigma_2^2 \end{pmatrix}

Suppose we draw a bunch of Gaussian variables, and we want to visualize them.

Two Independent Standard Normals

\boldsymbol \mu= \begin{pmatrix} 0 \\ 0 \end{pmatrix} \boldsymbol \Sigma= \begin{pmatrix}1 & 0\\ 0 & 1 \end{pmatrix}

Suppose we draw a bunch of Gaussian variables, and we want to visualize them.

Two Different Independent Normals

\boldsymbol \mu= \begin{pmatrix} 0 \\ 1 \end{pmatrix} \boldsymbol \Sigma= \begin{pmatrix}1 & 0\\ 0 & 0.2 \end{pmatrix}

Two Correlated Standard Normals

\boldsymbol \mu= \begin{pmatrix} 0 \\ 0 \end{pmatrix} \boldsymbol \Sigma= \begin{pmatrix}1 & 0.9\\ 0.9 & 1 \end{pmatrix}

Three or More Variables

• Hard to visualize in dimensions > 2, so stack points next to each other.

\boldsymbol \mu= \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \boldsymbol \Sigma= \begin{pmatrix}1 & 0.99\\ 0.99 & 1 \end{pmatrix}

y1 and y2 are highly correlated, y1 and y2 have similar values, and so when we put them next to each other, the lines do not go up or down that much. Instead, the lines are quite horizontal.

d = 5

\boldsymbol \mu= \begin{pmatrix} 0 \\ 0 \\ 0\\0\\0\end{pmatrix} \quad \quad \boldsymbol \Sigma= \begin{pmatrix}1 & 0.99 & 0.98 & 0.97 & 0.96\\ 0.99 & 1 & 0.99 & 0.98 & 0.97 \\ 0.98 & 0.99 & 1 & 0.99 & 0.98 \\ 0.97 & 0.98 & 0.97 & 1 & 0.99\\0.96 & 0.97 & 0.98 & 0.99 & 1 \end{pmatrix}

• Each line is one sample (path).

d = 50

\boldsymbol \mu= \begin{pmatrix} 0 \\ 0 \\ \vdots \\0\\0\end{pmatrix} \quad \quad \boldsymbol \Sigma= \begin{pmatrix} 1 & 0.99 & 0.98 & 0.97 & 0.96 & \cdots \\ 0.99 & 1 & 0.99 & 0.98 & 0.97 & \cdots\\ 0.98 & 0.99 & 1 & 0.99 & 0.98 & \cdots\\ 0.97 & 0.98 & 0.97 & 1 & 0.99 & \cdots\\ 0.96 & 0.97 & 0.98 & 0.99 & 1 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \end{pmatrix}

• Each line is one sample (path).

• Think of Gaussian processes as an infinite dimensional distribution over functions

  • all we need to do is change notation

Gaussian Process (GP)

A Gaussian process (GP) is a process f(\cdot): \mathbb{R}^p \rightarrow \mathbb{R} where any finite n-dimensional distributions are Gaussian, for any x_1\dots, x_n \in \mathbb{R}^p:

f(x_1), \dots, f(x_n) \sim N_n\left(\boldsymbol \mu, \boldsymbol \Sigma\right)

• Write f(\cdot) \sim GP to denote that the function f is a GP.

• The underlying true ERP waveform is assumed a GP.

• Before getting into our model, let me briefly talk about GP.

• Here shows a one-dimensional white noises and a realization or sample path of a GP.

• Here shows a two-dimensional white noises and a realization or sample path of a GP.

Mean and Covariance Function

• To fully specify a GP, need the mean and covariance function, f(\cdot) \sim GP\left(m(\cdot), k(\cdot, \cdot)\right) where

\text{E}(f(x)) = m(x) \text{Cov}(f(x), f(x')) = k(x, x')

• Popular choices of m(\cdot) are m(x) = 0 or \text{const} for all x, or m(x) = \beta'x

• Care more about the covariance function or kernel as it governs the how the process looks like by defining the similarity between data points.

Squared Exponential (SE) Kernel

k(x, x' \mid \tau, h) = \tau^2 \exp\left(-\frac{(x - x')^2}{2h^2} \right)

k(x, x') = \exp\left(-\frac{1}{2}(x - x')^2 \right)

• The parameter h is the characteristic length-scale that controls the number of level-zero upcrossings.

Squared Exponential (SE) Kernel

k(x, x') = \exp\left(-\frac{1}{2}\frac{(x - x')^2}{ {\color{red}{0.25}}^2} \right)

Squared Exponential (SE) Kernel

k(x, x') = \exp\left(-\frac{1}{2}\frac{(x - x')^2}{ {\color{red}{4}}^2} \right)

Squared Exponential (SE) Kernel

k(x, x') = {\color{red}{10^2}}\exp\left(-\frac{1}{2}(x - x')^2 \right)

• The parameter \tau^2 is the variance of f(x) that controls the vertical variation of the process.

Derivative Gaussian Process (DGP)

• Force f' = 0 at x = -4, 0, and 4, and all paths follow the constraint.
• The stationary points can be a maximum, minimum or inflection point.
• Add 3 constraints, but a path will have at least 3 stationary points.

• When there are some observations, these data points are assumed the realized values of a Gaussian process, and so when we simulate a GP path, the sample path must go through the data points if there is no noises or random errors.
• DGP is basically a GP with derivative constraints or information.
• For example here on the left, we have GP sample paths of SE kernel without derivative constraints.
• But we could add derivative information to the GP.
• Like the one on the right that (…)

Shaped-Constrained DGP Regression

\begin{aligned} y_i &= f(x_i) + \epsilon _i, ~~ \epsilon _i \sim N(0, \sigma^2), ~~ i = 1, \dots, n, \\ f(\cdot) &\sim GP\\ f'(t_m) &= 0, ~~m = 1, \dots, M. \end{aligned}

• In ERP applications,
  • \{y_i\}_{i=1}^n: the ERP time series observations.
  • f(\cdot): the true underlying smooth ERP waveform.
  • t_1, \dots, t_M: the latency of the ERP components.

• Data/Likelihood level is a GP regression with the random error term that is assumed N(0, \sigma^2)
• Then the regression function follows a DGP prior
• The idea is that the regression function can be quite flexible with a GP prior, but the function cannot be too flexible because we think there are M stationary points at t1, …, tM. Note that t_1, …, t_M are parameters we are interested and want to learn.

Prior on Stationary Points

• t_1, \dots, t_M are parameters to be estimated.

• In ERP studies, M is usually known.

• Use an univariate prior t \sim \pi(t) that corresponds to assuming that f has at least one critical point, and utilize the posterior of t to infer all critical points.

An univariate prior of t leads to

• Efficient computation:

  • The covariance matrix only adds one more dimension to from n by n to (n+1) by (n+1).
  • Simple one-dimensional sampling on t.

• Avoids possible misspecification of M that distorts the fitted result.

• Theorem: The posterior of t converges to a mixture of Gaussians as n \rightarrow \infty, and the estimated number of mixture components converges to the number of stationary points.

Mixture of Gaussians

Posterior Inference via Monte Carlo EM

\begin{aligned} y_i &= f(x_i) + \epsilon_i, ~~ \epsilon_i \sim N(0, \sigma^2), ~~i = 1, \dots, n, \\ f &\sim GP(0, k(\cdot, \cdot; \boldsymbol \theta)) , ~~ f'(t) = 0,\\ t &\sim \pi(t), ~~ \sigma^2 \sim \pi(\sigma^2), ~~ \boldsymbol \theta\sim \pi(\boldsymbol \theta) \end{aligned}

For iteration i = 1, 2, \dots until convergence:

• Monte Carlo E-step: Metropolis-Hastings on t and Gibbs sampling on \sigma^2.

• M-step: Given samples of t and \sigma^2, update \boldsymbol \theta to \hat{\boldsymbol \theta}^{(i+1)} by setting \begin{aligned} \hat{\boldsymbol \theta}^{(i+1)} &= \arg \max_{\boldsymbol \theta} \frac{1}{J} \sum_{j=1}^J \log p(\bm{y} \mid t_j^{(i)}, (\sigma_j^2)^{(i)}, \boldsymbol \theta^{(i)}) \end{aligned}

• Data/Likelihood level is a GP regression with the random error term that is assumed N(0 \sigma^2)
• Then the regression function follows a DGP prior where t is treated as a random variable.
• After t is estimated, we learn when those stationary points happen.
• And after \theta is estimated, we are able to estimate the entire function, or the ERP waveform.

Simulation: Uncertainty Quantification of t and f

Two stationary points t_1 = 0.44 and t_2 = 1.46

\pi(t) = \text{Unif}(0, 2)

• f(x) = 0.3 + 0.4x+0.5\sin(3.2x) + 1.1 /(1 + x^2)
• x\in [0, 2]

Oracle DGP

• Better curve fitting if true critical points were known.

One thing I’d like to mention is that if the true s.p were known to us, which we called oracle DGP, its curve fitting would be better than the basic GPR. The reason is that now we have more information about how the regression function looks like, and the information is 100% correct. So in general, the uncertainty band of oracle DGP will be narrower than GPR or DGP when s.p is uncertain.

Model Extension

Current ERP Studies

• Multi-subject
  • Averaging across subjects to obtain the grand ERP waveform.
• Difference in latency and/or amplitude of ERP components for factors such as age, gender, etc.
  • Basic ANOVA is conducted.
• A time search window needs to be specified to frame the range of the components.
  • No systematic and consistent ways of selecting such windows.

Source: Rossion & Jacques (2012)

The model we just talked about is the basic skeleton for ERP analysis. But we need a more sophisticated one for addressing current ERP study issues.

Semiparametric Latent ANOVA Model (SLAM)

• Estimate all parameters in a unified framework.

• Supports a multi-factor ANOVA setting that is capable of examining how multiple factors affect the latency and amplitude.

• Bayesian probabilistic statements about the amplitude and latency

• The hierarchical model provides estimates at both group and subject levels.

• Simultaneously estimates the amplitude and latency as well as the ERP waveform in a unified framework.

• Supports a multi-factor ANOVA setting that is capable of examining how multiple factors or covariates affect ERP components’ latency and amplitude.

• Bayesian approach having probabilistic statements about the amplitude and latency of ERP components.

• The hierarchical nature of the model provides estimates at the grand group level as well as the individual subject level that is masked by grand averaging.

SLAM One-way ANOVA - Likelihood with DGP

y_{igs} = f_{gs}(x_i) + \epsilon_{igs}, ~~ \epsilon_{igs} \stackrel{\rm iid}{\sim}N(0, \sigma^2).

• Time point i = 1, \dots, n

• Factor level g = 1, \dots, G, (g = \color{blue}{\text{young}, \text{old}} for Age factor)

• Subject s = 1, \dots, S_{g} (\color{blue}{S_{young} = 18; S_{old} = 11})

f_{gs} \mid \bm{t}_{gs} \sim \text{DGP}(0, k(\cdot, \cdot; \boldsymbol \theta), \bm{t}_{gs}), where \bm{t}_{gs} = \{t_{gs}^m\}_{m = 1}^M.

• M is pre-specified.

SLAM One-way ANOVA - Prior

\begin{aligned} t_{gs}^m \mid r_{g}^m, \eta_{g}^m &\stackrel{\rm iid}{\sim}\text{gbeta}\left(r_{g}^m\eta_{g}^m, \, (1 - r_{g}^m)\eta_{g}^m, \, a^{m}, \, b^{m}\right), \, s = 1, \dots, S_{g}. \end{aligned}

• Subject-level t_{gs}^m \in [a^{m}, b^{m}] (\color{blue}{t_{young, 1}^{N100} \in [60, 140]; t_{old, 2}^{P200} \in [140, 220]})

• Group-level r_{g}^m \in (0, 1): prior mean location in [a^{m}, b^{m}]

• {\mathbb E}\left(t_{gs}^m\right) = (1 - r_{g}^m)a^{m} + r_{g}^mb^{m}

• For the subject-level parameters, we assume it follows general beta distribution that has support a^m, b^m interval, which can be viewed as the most coarse search window for a latency.
• The shape parameters are parameterized using r and \eta, so that the group-level parameter r indicates the location of the prior mean in the search window, and tells us a or b is closer to the mean location.
• It represents the relative location of the interval, or the weight between the two end points of the search window.

SLAM One-way ANOVA - Link Function

With the link function \phi(\cdot) \in (-\infty, \infty), and z_j \in \{0, 1\}, j = 1, \dots, G-1,

\phi\left(r_{g}^m\right) = \beta_0^m + \beta_{1}^mz_{1} + \cdots + \beta_{G-1}^mz_{G-1}.

\color{blue}{\phi\left(r_{young}^{N100}\right) = \beta_0^{N100}}.

\color{blue}{\phi\left(r_{old}^{N100}\right) = \beta_0^{N100} + \beta_{1}^{N100}}.

• Accommodate any linear model with categorical and numerical variables

• Being able to describe interactions

• Any link function with domain (0, 1) can be used, logit, probit, complementary log-log, etc.

• To do the inference about r, we introduce an ANOVA structure to model how the factor the location of ERP component.
• Since r is bounded between 0 and 1, we use a link function \phi() followed by an ANOVA which is written as a regression or linear model using dummy variables.

SLAM One-way ANOVA - Other Priors

\begin{aligned} \beta_{j}^m &\sim N(\mu_j^m, (\sigma_j^m)^2), ~~ j = 0, 1, \dots, G-1\\ \eta_{g}^m &\sim Ga(\alpha_{\eta}, \beta_{\eta}),\\ \pi(\sigma^2, \tau^2, h) &= \pi(\sigma^2) \pi(\tau^2) \pi(h) = IG(\alpha_{\sigma}, \beta_{\sigma})IG(\alpha_{\tau}, \beta_{\tau})Ga(\alpha_{h}, \beta_{h}). \end{aligned}

SLAM One-way ANOVA - Inference

• Still use Monte Carlo EM for posterior inference.

• Everybody’s tired. Check algorithm details in the paper. Let’s see some analysis results!

For iteration i = 1, 2, \dots until convergence:

• Monte Carlo E-step:

  • Metropolis-Hastings on \mathbf{t}, \beta_{j}^m, \eta_g^m, and Gibbs sampling on \sigma^2.

• M-step: Given samples of \mathbf{t} and \sigma^2, update \boldsymbol \theta to \hat{\boldsymbol \theta}^{(i+1)} by setting \begin{aligned} \hat{\boldsymbol \theta}^{(i+1)} &= \arg \max_{\boldsymbol \theta} \frac{1}{L} \sum_{l=1}^L \log p(\bm{y} \mid \mathbf{t}_l^{(i)}, (\sigma_l^2)^{(i)}, \boldsymbol \theta^{(i)}) \end{aligned}

where

\bm{y}_{gs} \mid \mathbf{t}_{gs}, \sigma, \boldsymbol \theta\stackrel{\rm ind}{\sim}N\left(\boldsymbol \mu_{gs}, \boldsymbol \Sigma_{gs}(\boldsymbol \theta) +\sigma^2\mathbf{I}\right).

Simulation - Data

• 10 subjects (regression functions) in each group

• Data of size n=100 are generated with \sigma = 0.25

Point and Interval Estimation of Subject-level Latency

Point and Interval Estimation of Subject-level Latency

Uncertainty Quantification of Group-level Latency

Latency root mean square error (RMSE)

Method	Sine	Cosine
SLAM-posterior mean	0.004 (0.0000)	0.007 (0.0002)
LOESS1	0.012 (0.0021)	0.025 (0.0036)

• Smaller RMSE with less variation.

Method	Sine	Cosine
SLAM-posterior mean	0.0037 (0.0000)	0.0068 (0.0002)
SLAM-median	0.0036 (0.0001)	0.0068 (0.0002)
LOESS	0.0119 (0.0021)	0.0251 (0.0036)

1. Jeste SS, Kirkham N, Senturk D, Hasenstab K, Sugar C, Kupelian C, Baker E, Sanders A, Shimizu C, Norona A, et al. 2014. Electrophysiological evidence of heterogeneity in visual statistical learning in young children with ASD. Developmental Science. 18(1):90–105

Curve Fitting

Amplitude Estimates

Amplitude root mean square error

Method	Sine	Cosine
SLAM-posterior mean	0.04 (0.0004)	0.06 (0.0016)
LOESS	0.07 (0.0021)	0.07 (0.0066)

• The amplitude estimation relies on a search window that is specified by the posterior of t or r.

Method	Sine	Cosine
SLAM-posterior mean	0.0421 (0.0004)	0.0582 (0.0016)
SLAM-median	0.0421 (0.0005)	0.0581 (0.0016)
LOESS	0.0661 (0.0021)	0.0654 (0.0066)

ERP Data

• Experiment on speech recognition. Choose the word you hear from

  • date-tate vs. dape-tape
  • gate-cate vs. gake-cake.

• N100 (60 ~ 140 msec) component captures phonological (syllable) representation.

• P200 (140 ~ 220 msec) component represents higher-order perceptual processing.

• The experiment involved 18 college-age students and 11 older controls, as aging is known to cause difficulties in perceiving speech, especially in noisy environments.
• The older subjects have ages ranging from 47 to 91 years old with a mean age of 67.78 years.
• The experiment resulted in a total of 2304 trials.
• The six electrodes Fp1, F3, F7, C3, FC5, and T7

Latency Difference Between Subjects

Latency Difference Between Groups

• The older group has longer N100 and P200 latency.

• The time lag between N100 and P200 in the young group is shorter.

Amplitude Difference Between Groups

• The older group has larger N100 amplitude.

Conclusion

• Novel Bayesian SLAM estimates the latency and amplitude of ERP components with uncertainty quantification.

• Solves the search window specification problem by generating posterior distribution of latency.

• Incorporates the ANOVA and possibly a latent generalized linear model structure in a single unified framework.

• Examines both individual and group differences that facilitate comparing different characteristics or factors, such as age/gender.

• R package slam at https://github.com/chenghanyustats/slam

• Working in progress:

  • trial variability
  • autoregressive correlated noises
  • spatial modeling and mapping on electrodes.