1. How HMP works: practical, methodological and theoretical grounds#

Data#

Simulation#

We will start by simulating data to demonstrate the HMP analysis on simulated data. First, we need to load several libraries, after which we simulate a dataset with 4 events and the resulting 5 inter-event times (stimulus -> event 1, event 1 -> event 2, …, event 4 -> response.

[!NOTE] Most of this code is specific to simulation cases, to apply to your own data, use the code of the next tutorials

[1]:

## Importing these packages is specific for this simulation case
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gamma

## Importing HMP
import hmp
from hmp import simulations

[2]:

cpus = 2 # For multiprocessing, usually a good idea to use multiple CPUs as long as you have enough RAM

n_trials = 100 #Number of trials to simulate
sfreq = 500 #Sampling rate of the simulated data

##### Here we define the sources of the brain activity (event) for each trial
n_events = 4 #How many events to simulate
frequency = 10. #Frequency of the event defining its duration, half-sine of 10Hz = 50ms
amplitude = .2e-7 #Amplitude of the event in nAm, defining signal to noise ratio
shape = 2 #shape of the gamma distribution, with the defined means below it dictates the onset time probabilities
means = np.array([60, 150, 200, 100, 80])/shape #Mean duration of the between event times in ms
#Which source to activate for each event (see atlas when calling simulations.available_sources())
names = ['inferiortemporal-lh','caudalanteriorcingulate-rh','bankssts-lh','superiorparietal-lh','superiorparietal-lh']

sources = []
for source in zip(names, means): #One source = one frequency, one amplitude and a given by-trial variability distribution
    sources.append([source[0], frequency, amplitude, gamma(shape, scale=source[1])])

# Function used to generate the data
file = simulations.simulate(sources, n_trials, cpus, 'dataset', overwrite=False, sfreq=sfreq, seed=0)
#load electrode position, specific to the simulations
info = simulations.positions()

/home/gabriel/ownCloud/projects/RUGUU/hmp/hmp/simulations.py:219: UserWarning: /home/gabriel/ownCloud/projects/RUGUU/hmp/docs/source/notebooks/./dataset_raw.fif exists no new simulation performed
  warn(f"{subj_file} exists no new simulation performed", UserWarning)

Next, we load the simulated data, reflecting a realistic EEG dataset of a single participant. As you can see, we use standard MNE functions here.

[3]:

#Recovering the events to epoch the data (in the number of trials defined above)
events = np.load(file[1])
resp_trigger = int(np.max(np.unique(events[:,2])))#Resp trigger is the last source in each trial
event_id = {'stimulus':1}#trigger 1 = stimulus
resp_id = {'response':resp_trigger}
#Keeping only stimulus and response triggers

#Visualising the raw simulated EEG data
import mne

raw = mne.io.read_raw_fif(file[0], preload=False, verbose=False)
raw.pick_types(eeg=True).plot(scalings=dict(eeg=1e-5), events=events, block=True);

../_images/notebooks_1-How_HMP_works_7_0.png

Data Format#

We now prepare the data for the HMP package. First, we read the data into an xarray.

[4]:

### Reading the data
from hmp import io

epoch_data = io.read_mne_data(file[0], data_format='raw', event_id=event_id, resp_id=resp_id, sfreq=sfreq,
            events_provided=events, tmin =-.5, tmax=2, verbose=False,)#We add an offset to ensure a smooth estimation at the end of the RT
print(epoch_data)

Processing participant /home/gabriel/ownCloud/projects/RUGUU/hmp/docs/source/notebooks/./dataset_raw.fif's raw eeg
<xarray.Dataset> Size: 30MB
Dimensions:      (participant: 1, epoch: 100, channel: 59, sample: 1251)
Coordinates:
  * epoch        (epoch) int64 800B 0 1 2 3 4 5 6 7 ... 92 93 94 95 96 97 98 99
  * channel      (channel) <U7 2kB 'EEG 001' 'EEG 002' ... 'EEG 059' 'EEG 060'
  * sample       (sample) int64 10kB -250 -249 -248 -247 ... 997 998 999 1000
    event_name   (epoch) object 800B 'stimulus' 'stimulus' ... 'stimulus'
    rt           (epoch) float64 800B 0.35 0.554 0.624 0.364 ... 0.48 0.468 0.36
  * participant  (participant) <U5 20B 'sub-0'
Data variables:
    data         (participant, epoch, channel, sample) float32 30MB -3.249e-0...
Attributes:
    sfreq:      500.0
    lowpass:    40.0
    highpass:   0.10000000149011612
    reference:  None
    n_trials:   100
    tmin:       -0.5
    tmax:       2

The previous function epochs the data around our (simulated) stimuli and stores any additional infomation along (here ‘event_name’).

xarray is a handy way of using matrices with named dimensions. For example we can now very easily draw an ERP:

[5]:

#example usage of xarray
epoch_data.sel(channel=['EEG 001','EEG 002','EEG 003'], sample=range(400))\
    .data.groupby('sample').mean(['participant','epoch']).plot.line(hue='channel')
plt.show()

../_images/notebooks_1-How_HMP_works_11_0.png

Data preprocessing#

Next, we transform the data for HMP:

We cut the data at the interval of interest (here the time between simulated stimuli and simulated response denoted as rt in epoch_data)
Apply PCA and select top components to lower computational load
zscore the data per trial

For the PCA we select 3 components for simplicity (more on this with real data in Tutorial 3).

[6]:

from hmp import transformers

transformed = transformers.ProjPCA(epoch_data, n_comp=3)

100 positively defined intervalsbetween 0.002 and inf seconds.
Rejection summary:
 0 trials rejected based on threshold of None
 0 trials rejected based on interval limit of (0.002, inf)
 0 trials detected with no interval (e.g. preprocessing or interval exceeding epoch))

[7]:

from hmp.visu import plot_components_sensor

plot_components_sensor(transformed.weights, info)
plt.show()

../_images/notebooks_1-How_HMP_works_15_0.png

From now on the estimation performed by HMP will be done on these 3 virtual channels

HMP and implementation assumptions#

Expected pattern#

HMP assumes that a transition or cognitive event is represented in the data as a specific pattern. Therefore, before estimating a model one needs to build such an expected event characteristics.

The standard 50ms half-sine can be created by invoking the HalfSine class. The class will then store the expected pattern given the sampling frequency of the data.

[8]:

from hmp.patterns import HalfSine

event_properties = HalfSine.create_expected(sfreq=epoch_data.sfreq)
event_properties

[8]:

HalfSine(sfreq=500.0, width=25, location=25, template=array([0.00502324, 0.01499051, 0.02472136, 0.03406234, 0.04286614,
       0.05099392, 0.05831749, 0.06472136, 0.07010453, 0.07438212,
       0.07748665, 0.07936918, 0.08      , 0.07936918, 0.07748665,
       0.07438212, 0.07010453, 0.06472136, 0.05831749, 0.05099392,
       0.04286614, 0.03406234, 0.02472136, 0.01499051, 0.00502324]))

event_properties.template contains the default shape that we expect in the data, but we can replace that with whatever template we want if we have specific hypothesis:

[9]:

plt.figure(figsize=(2,3))
plt.plot(event_properties.template, 'x', label='default')
plt.ylabel('Normalized value')
plt.xlabel('Samples NOT time')

# Alternative template, just for fun
from scipy import signal

fancy = signal.windows.flattop(len(event_properties.template))
fancy = fancy/np.sum(fancy**2)#Normalizing as the original one
plt.plot(fancy, 'x', label='fancy')
plt.legend(bbox_to_anchor=(1,1))
plt.show();

../_images/notebooks_1-How_HMP_works_22_0.png

If we wanted the fancy one we could simply do:

from hmp.patterns import Arbitrary
event_properties = Arbitrary.create_expected(sfreq=epoch_data.sfreq, template=fancy)

We then create a TrialData class that combines the expected pattern with the data. The class will automatically cross-correlate the data with the provided template.

Below, we plot those cross-correlated values; the vertical lines indicate the start of the events in the simulated data.

[10]:

from hmp import transformers
from hmp.trialdata import TrialData

trial_data = TrialData.from_transformer(transformed.data, pattern=event_properties.template)
plt.plot(trial_data.cross_corr[trial_data.starts[0]:trial_data.ends[0]])

#Recover the actual time of the simulated events
sim_source_times = np.reshape(np.ediff1d(events[:,0],to_begin=0)[events[:,2] > 1], (n_trials, n_events+1))
plt.vlines(sim_source_times[0,:-1].cumsum(), -2.5, 2.5, 'k') #overlaying the simulated event onset times
plt.show()

../_images/notebooks_1-How_HMP_works_25_0.png

The specific pattern to be expected is an open discussion, given the transient and oscillatory nature of E/MEG signals, a half-sine is a good starting point. The specific width (or frequency) of the half-sine is more crucial (see tutorial on assumption and the simulations in the paper), overall the shorter the pattern the better up to the point were high frequency noise gets in the way of the estimation.

Based on experience, a pattern width between 30 and 50ms generally provides the most robust solutions (but caution regarding the number of points in the pattern given the sampling rate, e.g. no use to take a 100Hz signal with a 30ms half-sine as the pattern would only be defined on three points)

Expected time distribution#

Next HMP assumes that the time of the event peak follow a specific one-parameter time distribution. By default we expect event peak times to vary per-trial according to a Gamma distribution with a shape of 2 and a scale free to vary:

[11]:

from hmp.distributions import Gamma

gamma_2 = Gamma()
gamma_2

[11]:

<hmp.distributions.Gamma at 0x72fd8b589c70>

For example an event happening on average 100ms following stimulus would have the blue distribution given a Gamma with a shape of 2 or the orange one with a Weibull and a shape of 5.

[12]:

example_mean = 100
x = np.linspace(0, 500, 1000)
plt.plot(x, gamma_2.pdf(x, gamma_2.shape,
        scale=gamma_2.mean_to_scale(example_mean)),
        label=r'Default Gamma$_2$')

## Other distribution to illustrate
from hmp.distributions import Weibull

other_time_distribution = Weibull(shape=5)
plt.plot(x, other_time_distribution.pdf(x, other_time_distribution.shape,
        scale=other_time_distribution.mean_to_scale(example_mean)),
        label=r'Example Weibull$_5$')
plt.legend()
plt.show()

../_images/notebooks_1-How_HMP_works_31_0.png

Overall the specific distribution doesn’t really matter (see tutorials on assumptions and the appendix of the paper). The most important thing is to not exclude likely times (e.g. the Weibull doesn’t allow event to happen directly after the stimulus for example).

[!NOTE] Declaring a time distribution is not a requirement, if you don’t provide one, HMP will assume a Gamma distribution with a shape of 2 and a scale free to vary. This is the default setting and works well in most cases.

Estimating a model without the ground truth#

in real data, to estimate this unkown set of parameters (time and channel parameters) we use the expectation maximization algorithm on the event probability as obtained through the Baum-Welsh algorithm. We start with neutral guess on the time and channel parameters (0’s for channel_pars, and mean RT / n_events for the time_pars).

[24]:

#r Defining a fresh model as we would do for real data
model = EventModel(event_properties, gamma_2, n_events=n_events)
loglikelihood, estimates = model.fit_transform(trial_data)#This function simultaneously estimates the paraneters and transform the data to event probability space

hmp.visu.plot_topo_timecourse(epoch_data, estimates, info)
plt.show()

Estimating 4 events model with 1 starting point(s)

../_images/notebooks_1-How_HMP_works_68_1.png

Also without providing the true magnitudes and parameters, the algorithm manages to recover the events, on average, in the simulated data.

EM algorithm#

Under the hood, HMP uses an expectation maximization algorithm to find the solution that best explains the data. The way it works is:

make a first proposal on the time location and channel contribution towards each event
Compute the most likely location of the events given this proposition
Update the parameters according to this most likely position
Repeat steps 2-3 up to the moment where the loglikelihood of the data stops increasing more than the tolerance criterion (default 1e-4)

More details about the algorithm can be found in the accompanying 2024 paper

We can inspect the traces of this EM algorithm to see how the log-likelihood of the fit developed with iterations. This shows how many iterations it needed and at what Log-likelihood (LL) it stopped according to the default tolerance on the increase in LL (default 1e-4).

[25]:

plt.plot(model.traces)
plt.ylabel('Log-likelihood')
plt.xlabel('EM iteration')
plt.show()

../_images/notebooks_1-How_HMP_works_72_0.png

Along this we can take a look at how the position of the events were updated, i.e. how the time parameters changed according to the evolution of the EM algorithm:

[26]:

plt.plot(model.time_pars_dev[:,0,:,1], label=['Interval '+str(x) for x in [1,2,3,4,5]])
plt.ylabel(r'Gamma$_2$ scale')
plt.xlabel('EM iteration')
plt.legend(loc=1)
plt.show()

../_images/notebooks_1-How_HMP_works_74_0.png

Here we see that the the first EM iteration starts with the assumption that all events are equally spaced in the RT. Then for example the first interval progressively became shorter while the second one became longer. As we progress in the EM iteration the time distribution parameters get more stable. We can increase the tolerance but this, most of the time, wouldn’t have much consequences except making the fit slower:

[27]:

#re-defining the model to remove previous run
lowtol_model = EventModel(event_properties, gamma_2, n_events=n_events, tolerance=1e-10)
lowtol_loglikelohood, lowtol_estimates = lowtol_model.fit_transform(trial_data)
plt.plot(lowtol_model.time_pars_dev[:,0,:,1], label=['Interval '+str(x) for x in [1,2,3,4,5]])
plt.ylabel(r'Gamma$_2$ scale')
plt.xlabel('EM iteration')
plt.legend(loc=1)
plt.show()

Estimating 4 events model with 1 starting point(s)

../_images/notebooks_1-How_HMP_works_76_1.png

Here the profile is exactly the same except that the EM algorithm went through a few more iterations.

Using the single trial estimation#

If we select for all trials the most likely sample to contain a transition (e.g. the maximum probability in the single trial plot with p(events) on the y-axis) we can take a look at the recovery of our simulated times on all trials:

[28]:

# Estimate the duration between events based on the most time for each event at each trial
estimated_times = hmp.utils.event_times(estimates, duration=True, mean=False, add_rt=True).T

fig, ax = plt.subplots(n_events+1,2, figsize=(10,2.5*n_events+1), dpi=100, sharex=True)
i = 0
ax[0,0].set_title('point-wise')
ax[0,1].set_title('Histograms')
ax[-1,0].set_xlabel('Real by-trial peak')

for event in estimated_times:
    ax[i,0].plot([np.min(event), np.max(event)], [np.min(event), np.max(event)],'--', color='k')
    ax[i,0].plot(sim_source_times[:,i], event, 'o', color=colors[i])
    ax[i,0].set_ylabel(f'estimated peak {i+1} time')
    ax[i,1].hist(event, color='red', bins=20, histtype='step', label='Estimated')
    ax[i,1].hist(sim_source_times[:,i].T, color='k', bins=20, histtype='step', label='True')
    i+= 1
ax[0,-1].legend()
plt.xlim(0,250)
plt.tight_layout()
plt.show()

../_images/notebooks_1-How_HMP_works_80_0.png

Overall the exact time of the simulated event is detected as the most likely time by HMP! Of course this will depend on the signal to noise ratio and the non-violation of HMP’s assumptions (see the simulations in the paper)

Saving and loading a model#

Depending on the size of the data (with main factor the sampling rate) the HMP model can take a while to estimate. Therefore, it is often useful to save the model and the estimates for later use.

[29]:

import pickle

with open("model.pkl", 'wb') as output:
    pickle.dump(model, output)


# Same goes fpr the estimates
with open("estimates.pkl", 'wb') as output:
    pickle.dump(estimates, output)

For those who are not comortable in continuing in python you also have the option to just export the estimated event probabilities to a csv for use in other languages/softwares

[30]:

io.save_eventprobs_csv(estimates, 'estimates.csv')
import pandas as pd

df_estimates = pd.read_csv('estimates.csv')
print(df_estimates)

Saved at estimates.csv
        sample  event participant  epoch  group  eventprobs
0            0      0       sub-0      0      0         0.0
1            0      0       sub-0      1      0         0.0
2            0      0       sub-0      2      0         0.0
3            0      0       sub-0      3      0         0.0
4            0      0       sub-0      4      0         0.0
...        ...    ...         ...    ...    ...         ...
253995     634      3       sub-0     95      0         0.0
253996     634      3       sub-0     96      0         0.0
253997     634      3       sub-0     97      0         0.0
253998     634      3       sub-0     98      0         0.0
253999     634      3       sub-0     99      0         0.0

[254000 rows x 6 columns]

1. How HMP works: practical, methodological and theoretical grounds#

Data#

Simulation#

Data Format#

Data preprocessing#

HMP and implementation assumptions#

Expected pattern#

Expected time distribution#

An HMP model#

HMP parameters#

Time parameters#

Channel parameters#

Transforming the data to event probability space using these parameters#

Visualizing the average time location and topographies associated with the estimated events#

Estimating a model without the ground truth#

EM algorithm#

Using the single trial estimation#

Saving and loading a model#