A metric for event detection quality

[GWeindel]

Context: It's sometimes hard to judge how much an event is supported by the EEG data as opposed to supported by the distributions. E.g. some trials might have very low probability over long times samples without a clear peak.

Proposition: create a measure of probability update due to the EEG data, a bit like a prior/posterior contraction index in Bayesian analysis.

Implementation: we can create, for each event, a null model where all magnitudes and parameters are used except the one from the studied event. We don't re-estimate the model but just compare the event probabilities of each trial in the estimated model vs the null model.
See for example:

fake_mags = np.zeros(estimates.magnitudes.shape)
fake_mags[:2] = estimates.magnitudes[:2]
null = init.fit_single(4, parameters=estimates.parameters, magnitudes=fake_mags, maximization=False)

randidx = list(np.random.randint(0, 20000, 5))
fig, ax = plt.subplots(5, sharex=True, sharey=False, dpi=200)
for i, idx in enumerate(randidx): 
    for ev in range(4):
        ax[i].plot(estimates.samples*1000/256, estimates.eventprobs.isel(trial_x_participant=idx)[:,ev], color=colors_ev[ev])
        ax[i].plot(null.samples*1000/256, null.eventprobs.isel(trial_x_participant=idx)[:,ev], ls='--', color=colors_ev[ev])
plt.xlim(0,1000)

In this case we're looking at how the event 3 improved with the estimation of magnitudes compared to just using the distribution and the convolution with the probabilties from the previous events

We could thus compute something like the ratio of the probabilities given the real estimates compared to the null model. In pseudo-code this would give something like:

x = argmax(estimates.eventprobs.sel(trial=idx))

update = estimates.eventprobs.sel(x)/null.eventprobs.sel(x)

What do you think? @jelmerborst @lvmaanen @rickdott amongst others

[GWeindel]

Using the following code on a real data set:

estimates_maxes = estimates.eventprobs.argmax('samples')
update_factor = np.zeros(estimates.eventprobs.sel(samples=0).shape)*np.nan
for event in estimates.event.values:
    fake_mags = np.zeros(estimates.magnitudes.shape)
    fake_mags[:event] = estimates.magnitudes[:event]
    estimates_test = init.fit_single(4, parameters=estimates.parameters, magnitudes=fake_mags, maximization=False)
    update_factor[:, event] = estimates.eventprobs.sel(samples=estimates_maxes.sel(event=event), event=event)-\
    estimates_test.eventprobs.sel(samples=estimates_maxes.sel(event=event), event=event)

Gives me an interesting plot for the distributions of this update factor for each event:

We can then even inspect participants that might be misfitted? (x = event, y=update_factor, color = participants)

[lvmaanen]

This seems like a similar problem as @rickdott worked on in the context of SoCOM where we wanted to quantify the probability that an event occurred. However, it seems that SoCOM also applies in cases where the time series under consideration are not only HMP events +/- 1, but can be quite different.

[rickdott]

Cool idea, difficult for me to interpret. How would a misfitted participant be reflected in their update factor? Seems like if its generally higher/lower than the mean the pattern found at the peak of that participant's HMP probabilities is different from the generally estimated pattern for that event?

[GWeindel]

I guess this calls for simulations, but basically events/trials/participant with low support from the EEG data will be closer to 0.

I changed the formula as a ratio or a subtraction still depends mostly on the speed of the event. This is an example where I compute (probs - null_probs) / null_probs

In this case it's quite obvious that the last event (x-axis) is very weakly informed by the actual EEG data and more by the gamma distributions compared to the other ones