Problems with Gaussian Random Walk Model

Questions version agnostic
,
11

Hi again. and thanks again for the feedback. I have managed to solve the problem, or sort of. My reasoning was going back to the basics for a moment. That is, pondering whether event related potentials (ERPs) can actually be modelled by a gaussian random walk (GRW). Although ERPs have several sources of noise that may be random, i.e. from unknown or uncertain sources (such as unrelated brain processes, unidentified electrical sources in the environment, etc.), they are not a stochastic process on their own (i.e. they are postsynaptic electrical activity elicited by a specific stimuli in a controlled experiment, etc.). So, maybe, a GRW wasn’t the appropriate prior after all. I changed the GRW prior for a Normal distribution, such as to have multivariate normal from the product with the Cholesky decomposition covariance prior (using the Cholesky prior directly produced convergence problems). The resulting model looks like this:

with pm.Model() as mod: #Multivariate Normal model with LKJ prior
    sd = pm.HalfNormal.dist(1.0)
    chol, corr, std = pm.LKJCholeskyCov("chol", n=E, eta=6.0, sd_dist=sd, compute_corr=True)
    cov = pm.Deterministic("cov", chol.dot(chol.T))
    m = pm.Normal('m', mu=0, sigma=1.0)
    g = pm.Normal("g", mu=m, shape=(S,E,G))
    M = pm.Deterministic("M", tt.dot(cov, g.T))
    e = pm.HalfNormal("e", 5.0)+1.0
    y = pm.Normal("y", mu=M, sigma=e, observed=amps)

with mod:
    trace = pm.sample(2000, tune=2000, chains=4, cores=8, init='adapt_diag')

I run some prior predictive checks to tune up the values assigned to each prior, and it works the best with a relatively high eta parameter for the Cholesky prior, and relatively high sigma for the likelihood’s HalfNormal. At a single electrode, the prior predictive look like this:

ppcs_pz
ppcs_pz1800×1200 323 KB

This is quite reasonable, as an amplitude variation between -10uV and 10uV is to be expected for ERPs such as P300. Naturally, the observations at the baseline should always be near zero (if baseline correction was applied), but observations should make these priors adapt. Indeed, the model samples very well, with all ESS over 2000 and all R_hats =~ 1. Here is an example of the approximated signal at electrode Pz:

pz_posterior
pz_posterior1800×1200 368 KB

An these are topographical maps of the scalp showing estimated activity averaged across a time window (350-750ms):

learners_topo
learners_topo1311×397 113 KB

non-learners_topo
non-learners_topo1313×397 125 KB

The only small issue I see is that the estimated peak of the Learners wave is slightly off mark, but not by much. The correlations from the Cholesky prior are way more sensible now, showing high correlations with adjacent electrodes, but negative correlations with frontal electrodes (as they tend towards the opposite polarity):

corrs_topo
corrs_topo1800×1200 219 KB

Sorry I have added so much unnecessary detail, at least for those already familiar with the stats used in this model, but I thought it may come handy for those who may run into this post trying to troubleshoot something or searching for alternative models for ERP data. I hope it helps, I’m sure this could use some extra improvements. I’ll keep working on those and if someone is interested, just give me a shout and I can post/share them (if I have any). Or if anyone has some advice, please just let me know. Thanks in advance.

3 Likes