# Analyzing EEG / MEG data with PyMC3

**URL:** <https://discourse.pymc.io/t/analyzing-eeg-meg-data-with-pymc3/7993>\
**Category:** Questions\
**Created:** [September 8, 2021, 8:37am UTC](https://discourse.pymc.io/t/analyzing-eeg-meg-data-with-pymc3/7993 "2021-09-08T08:37:43Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![thht](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/thht/32/4417_2.png) [@thht](https://discourse.pymc.io/u/thht)\
**Post date:** [September 8, 2021, 8:37am UTC](https://discourse.pymc.io/t/analyzing-eeg-meg-data-with-pymc3/7993/1 "2021-09-08T08:37:43Z")

</div>

Hi,

just to provide some background: I am a scientist working with EEG / MEG data. I recently came across PyMC3 and was wondering whether we can apply this to our data.

The probably most common use case is this: I collected data from N participants in a paradigm that has 2 conditions. The data is sampled from a number of sensors (102 magnetometers in my case). Because the data is very noisy, I average over something like 100 trials within each of the conditions for each participant separately.

So, what I end up with is a matrix with the shape of (n\_participants, n\_channels, n\_samples) for each condition. Because I am interested in the difference between the conditions, I use the difference of the two matrices for further analysis.

The hypothesis I want to test is that at some channels and some samples the mean of the differences over all participants is not 0.

I can quite safely assume that my the differences follow a student-t distribution. So, I tried to come up with a model and used the models defined in the [BEST toolbox](https://pypi.org/project/best/) as a template.

This is what I cam up with:

```python
# normalisation seems to be needed because the values are in the range of 1e-15 and that seems to be a problem for the sampling algorithm
norm_factor = np.std(data_for_model)
data_for_model = data_for_model / norm_factor

nu_min = 2.5
nu_mean = 30
nu_mean - nu_min

with pm.Model() as model:
    mu = pm.Normal('mu', mu=0, sigma=1, shape=(n_channels, n_samples)) # we want to estimate the mean for each channel and sample
    std = pm.Uniform('std', lower=1e-3, upper=1e3, shape=(n_channels, 1)) # but we can assume that the std deviation is the same within each channel
    nu = pm.Exponential('nu', 1 / (nu_mean - nu_min)) + nu_min
    modeled_data = pm.StudentT('modeled', mu=mu, nu=nu, sigma=std, observed=data_for_model, total_size=data_for_model.shape)

    trace = pm.sample(1000)

```

Here are my questions:

1. Does this make sense to you? Or did I miss something important? This is my first attempt…
2. It seems that it is important to scale the really tiny values because otherwise, I get lots of errors from the sampler. Does that make sense?
3. Any way to speed this up? It takes like 15 hours on my laptop. I could send it to a quite powerful cluster that has ~200 cores. But it seems that the `cores` argument to `pm.sample` then only increases the number of chains, so it would not be done quicker. I already looked into providing a `multiprocessing` context but I need to limit the amount of cores available to it and I could not find out, how to do this.

Any help is highly appreciated! Thanks in advance!

Thomas

---

<div class="post-metadata">

**Author:** ![vb690](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/vb690/32/3432_2.png) [@vb690](https://discourse.pymc.io/u/vb690)\
**Post date:** [September 10, 2021, 7:35am UTC](https://discourse.pymc.io/t/analyzing-eeg-meg-data-with-pymc3/7993/2 "2021-09-10T07:35:09Z")

</div>

Hi @thht

if you don’t mind the opinion from a PyMC3 novice here are some ideas

1. You might want to consider a more informative prior for the scale parameter, at the moment is **extremely** wide. I usually go for an HalfCauchy distribution. You can center it on a very small value if you are afraid that large values can cause problems.
2. Given the nature of EEG signal (if I recall correctly from my master) you might want to model the correlation between electrodes.
3. It might be an overkill here, but PyMC3 offers [approximate variational inference](https://docs.pymc.io/notebooks/variational_api_quickstart.html) which can provide consistent speed-up.

Finally,  
can you expand a bit on the nature of n\_samples? I struggle to depict the experimental design (I was expecting a (n\_participants, n\_channels) matrix).

---

<div class="post-metadata">

**Author:** ![thht](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/thht/32/4417_2.png) [@thht](https://discourse.pymc.io/u/thht)\
**Post date:** [September 10, 2021, 8:01am UTC](https://discourse.pymc.io/t/analyzing-eeg-meg-data-with-pymc3/7993/3 "2021-09-10T08:01:55Z")

</div>

hi,  
thanks for the input.

1. What do you mean by the “scale” parameter? The `std`?
2. Yes, you are right. and for the MEG, I actually have a noise covariance matrix that, so i know how they covary/correlate. But i have no idea, how to model that and/or use that strong prior information. do you have any idea?
3. Actually, using variational inference with minibatches is the way I do it now. It decreases the computation time from ~20 hours to 20 minutes.

Concerning your question regarding the shape of the matrix: EEG or MEG data is sampled at a regular interval like 1000Hz. It is the strength of these methods that they allow you to infer when an effect happens. So, let’s say, I want to analyze the cortical response to an auditory (e.g., a sound) compared to a visual (an image, for instance) stimulus. I would present the audio ~100 times and the visual ~100 times as well. I have triggers in my data, so I know exactly, when this stimulation happened. I would then cut out the data around these triggers (like 300ms before the stimulation to 700ms after). Because of the very low signal to noise ratio, I would then average the 100 visual and 100 auditory trials within each participant. Assuming a sampling rate of 1000Hz, I would have 1000 samples per channel for each condition and participant.

This leaves me with a (n\_channels, n\_samples) matrix per participant and condition, which reduces to just one of these matrices per participant before I apply the statistics, because I am interested at what channels (where?) and what samples (when?) there is a difference between the two conditions.

So, in the end, I concatenate all these matrices for all participants, which leaves me with said (n\_participants, n\_channels, n\_samples) matrix.

Does this make it clearer?

---

<div class="post-metadata">

**Author:** ![vb690](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/vb690/32/3432_2.png) [@vb690](https://discourse.pymc.io/u/vb690)\
**Post date:** [September 10, 2021, 9:21am UTC](https://discourse.pymc.io/t/analyzing-eeg-meg-data-with-pymc3/7993/4 "2021-09-10T09:21:44Z")

</div>

Hi

yes that is very clear now.

1. Yes, by scale parameter I meant the `std`, might be quite challenging to sample from an uniform prior in the range 0.001 - 1000.

2. I only played around with it a couple of times, but I believe you could use a Multivariate Normal with your noise matrix as covariance matrix. The lecture by [Richard McElreath on Gaussian Processes](https://www.youtube.com/watch?v=pwMRbt2CbSU) is the most accessible resource I found so far on the subject.

3. The same logic I suspect it also applies to samples, since they represent different points in time you might expect some form of temporal structure to be present. This would probably require a spatio-temporal model but this is something **way** out of my depth.

---

<div class="post-metadata">

**Author:** ![Simon](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/simon/32/5566_2.png) [@Simon](https://discourse.pymc.io/u/Simon)\
**Post date:** [December 7, 2022, 11:42am UTC](https://discourse.pymc.io/t/analyzing-eeg-meg-data-with-pymc3/7993/5 "2022-12-07T11:42:18Z")

</div>

Hi, I’m more than a year late, but just in case you’re still interested in this topic I’ll leave a comment here. If you’re interested in the difference between conditions, the best way to conduct inference in that respect, to avoid a big loss of information (acknowledged uncertainty), would be to include both conditions in your model and then take the difference between the posteriors of each condition (ex-post, rather than taking the difference between groups from the observed signal). In that case you could do something like:

```auto
with pm.Model() as model:
    mu = pm.Normal('mu', mu=0.0, sigma=1.0, shape=(n_channels, n_samples, n_conditions)) 
    std = pm.HalfNormal('std', 1.0)
    nu = pm.Exponential('nu', 1 / (nu_mean - nu_min)) + nu_m
    y = pm.StudentT('y', mu=mu, nu=nu, sigma=std, observed=data)

    trace = pm.sample(1000)

```

Ideally, however, you would like to have the mu parameter as a linear regression formula (i.e. with an intercept and a slope), so it can tell the change induced as a function of condition, that could allow you to add varying intercepts and/or slopes by stimuli and participant as well. I have taken an approach like that here: [https://onlinelibrary.wiley.com/doi/full/10.1111/ejn.15470](https://onlinelibrary.wiley.com/doi/full/10.1111/ejn.15470) (that still needs much improvement though). If you’re interested in accounting for the time domain across the epoch, then you could give a look to this: [Bayesian models for event-related potentials time-series and electrode correlations estimation | bioRxiv](https://www.biorxiv.org/content/10.1101/2022.08.02.502520v1.full) . That implements a method for getting covariance/correlation matrices via LKJ priors. These could also be implemented in more conventional linear regression models, I would also recommend a lecture by McElreath on this topic, though a different one: [Statistical Rethinking 2022 Lecture 14 - Correlated Varying Effects - YouTube](https://www.youtube.com/watch?v=XDoAglqd7ss&t=2011s).

Sorry if I’m too late and this is irrelevant now. But maybe someone else may find it useful as well.