02: Reading the results

Get the numbers out of a fit, one by one or all in a table.

Introduction

In the first tutorial, report did everything at once: fit, summary and figure.

Here we look at what a fit gives back, and how to use it: read one number, check that the fit is good, keep all the numbers in a table, and choose the part of the curve to describe.

import numpy as np

from evoparam import fit_response, gen_response

Fitting without drawing

fit_response does only the fit: no figure, no summary. It is the one to use when you have many curves, or when you only want the numbers.

We use the same curve as in the first tutorial: a sharp bump at 25 ms and a broad one at 150 ms.

times, signal = gen_response([[25, -150, 10], [150, -80, 80]], seed=3)

fit = fit_response(times, signal, preset="ccep")

fit now holds everything evoparam found. print(fit) still shows the summary at any time.

print(fit)

==================================================================
                    evoparam - RESPONSE MODEL

                     fit window: 9 to 500 ms
          background: exponential, time constant 35.4 ms

                      2 component(s) found:
            label        latency     amplitude     width
            N1           25.0 ms    -152.1 µV     9.5 ms
            N2          151.9 ms     -77.4 µV    79.2 ms

           fit quality: R2 = 0.94, noise level = 7.4 µV
  (latency: time of the peak; width: full width at half maximum)
==================================================================

One trial, or the average?

fit_response describes one curve: one row of numbers over time. That curve can be a single trial (the response to one stimulus) or the average of many trials. The tool does not mind which; it fits what it is given.

In practice you give it the average of the trials of one electrode. One trial is mostly noise: the brain’s own activity is often larger than the response. Averaging keeps what comes back at every trial (the response) and cancels what does not (the noise): with 20 trials, the noise is about 4.5 times smaller (the square root of 20).

Below we make 20 noisy trials of the same response, then fit one trial alone and the average of the 20.

trials = [
    gen_response([[25, -150, 10], [150, -80, 80]], noise=60, seed=trial)[1] for trial in range(20)
]

fit_one = fit_response(times, trials[0], preset="ccep")

average = np.mean(trials, axis=0)
fit_average = fit_response(times, average, preset="ccep")

print("one trial:", len(fit_one.components), "bump(s)")
print("average of 20:", len(fit_average.components), "bump(s)")

one trial: 0 bump(s)
average of 20: 2 bump(s)

On the single trial the tool finds nothing: the response is there, but hidden in the noise, and evoparam does not invent bumps it cannot tell from noise. On the average it finds both bumps again.

The average gives the clean picture. To follow how each trial differs (does the N1 come a little earlier or later from one trial to the next?), evoparam has a separate function, fit_trials: it starts from the bumps of the average, then measures them again in every trial, and says when the noise is too strong to tell the trials apart. It works on all the trials of a recording at once (a TrialSet), cut around each stimulus beforehand.

Reading the bumps

The bumps are in fit.components, a list in order of time. Each bump has the three numbers of the first tutorial:

• latency: when the bump peaks;

• amplitude: its height above the background, in µV;

• fwhm: its width at half its height (“full width at half maximum”).

Note

Inside evoparam, times are in seconds (latency = 0.025 means 25 ms). To read them in milliseconds, add _ms to the name: latency_ms, fwhm_ms, onset_ms…

for bump in fit.components:
    print(bump.label, round(bump.latency_ms, 1), "ms,", round(bump.amplitude, 1), "µV")

N1 25.0 ms, -152.1 µV
N2 151.9 ms, -77.4 µV

To keep one number, take the bump by its place in the list: [0] is the first one.

n1 = fit.components[0]
print(round(n1.latency_ms, 1))

25.0

Each bump also has a few more numbers, all measured on the bump alone:

• onset_ms and offset_ms: when the bump starts and ends;

• area: the surface of the bump, the coloured area of the first tutorial (µV × seconds);

• rise_time_ms and decay_time_ms: how long it takes to reach its peak, and to come back.

print(round(n1.onset_ms, 1), round(n1.offset_ms, 1))

16.3 33.6

Reading the background

The background is in fit.background. Here it is a curve that fades away, and the number to read is its time constant: how long it takes to lose about two thirds of its height. The smaller it is, the faster it fades.

print(round(fit.background_ms["tau"], 1))

35.4

Is the fit good?

Before using the numbers, check that the model really follows the curve. Two numbers in fit.gof (“goodness of fit”) say it:

• r2: how much of the curve the model explains, from 0 (nothing) to 1 (everything);

• rmse: the typical distance between the curve and the model, in µV. It should be close to the level of the noise.

print(round(fit.gof.r2, 2), round(fit.gof.rmse, 1))

0.94 7.1

Note

Here R² is about 0.94: the model explains 94 % of the curve, and what is left is the noise we put in. A low R² means the model missed something: look at the figure (report) before trusting the numbers.

All the numbers in a table

to_dataframe puts every number of every bump in a table, one row per bump. This table can be sorted, plotted or saved like any spreadsheet (it needs the pandas package, already there in most notebooks).

table = fit.to_dataframe()
table[["label", "latency", "amplitude", "fwhm", "area"]]

	label	latency	amplitude	fwhm	area
0	N1	0.024951	-152.057403	0.009520	-1.491724
1	N2	0.151916	-77.413645	0.079193	-6.317850

To keep it, save it as a file that any spreadsheet program opens:

table.to_csv("my_results.csv")

Choosing the part of the curve to describe

The preset decides which part of the curve is described: for "ccep", from 9 to 500 ms after the stimulus. To look at another part, give a window, in seconds.

Here we keep only the first 100 ms: the broad bump at 150 ms is left out, and only the sharp one is described.

early = fit_response(times, signal, preset="ccep", window=(0.009, 0.1))
print(early)

==================================================================
                    evoparam - RESPONSE MODEL

                     fit window: 9 to 100 ms
          background: exponential, time constant 50.0 ms

                      1 component(s) found:
            label        latency     amplitude     width
            N1           24.9 ms    -152.4 µV     9.5 ms

           fit quality: R2 = 0.97, noise level = 7.4 µV
  (latency: time of the peak; width: full width at half maximum)
==================================================================

Note

The preset matters too: it holds the settings made for one kind of recording. "ccep" is for responses to electrical pulses in the brain; "erp" is for responses recorded on the scalp (EEG), "sep" for touch, "tep" for magnetic pulses (TMS), and there are a few more.

Conclusion

A fit gives back its bumps (fit.components), its background (fit.background) and how good it is (fit.gof), and to_dataframe puts everything in one table. With window and preset, you choose what is described.