neuropsychology / NeuroKit

NeuroKit2: The Python Toolbox for Neurophysiological Signal Processing
https://neuropsychology.github.io/NeuroKit
MIT License
1.56k stars 419 forks source link

Confusion between exhalation peaks/troughs and inhale/exhale onsets [RSP] #528

Closed shanelindsay closed 3 years ago

shanelindsay commented 3 years ago

Neurokit labels peaks in the respiration signal as inhalation peaks and troughs as exhalation troughs (e.g. the rsp_plot graph). The example https://neurokit2.readthedocs.io/en/latest/examples/rrv.html doesn't specify the nature of the respiration measurment (but perhaps should), though here refers to inspiration peaks and exhalation troughs. This had potential for confusion, as peaks and troughs in a respiration signal have different meaning depending on the type of recording and the measure (i.e. volume, spirometer or breathing belt). The code which calculates whether a sample is inspiration or expiration assumes that a trough is the start of inhalation, which would the case for breathing belt data (and I assume spirometer users might not be using neurokit for analysis).

However, some of the code (if I am reading it right) seems to assume that "peaks" are inhalation onsets. For example, https://neurokit2.readthedocs.io/en/latest/_modules/neurokit2/rsp/rsp_rrv.html#rsp_rrv bbi = np.peaks(peaks) / sampling rate * 1000 Here bbi is breath to interval, which is used to calculate a number of measures. Traditionally these measures would use the onset of the inhalation as the start of a breath cycle not the start of the exhale (i.e. the "inspiration peak" in neurokit terminology), so it seems the code should be np.peaks(troughs). So this suggests a few possible changes - fixing the code to use troughs not peaks (and I think a few other functions are affected), being more explicit about the respiration recording type assumed in the docs, and perhaps shifting more to labelling peaks and troughs as inhale onsets and exhale onsets instead of inspiration/exhalation troughs and inspiration peaks (which would reduce ambiguity of the meaning of a "inspiration/exhalation peak/trough"). [note, some functions might work very similar regardless of using peaks and troughs, but in terms of reporting of measures inhale onsets would be better to fit with standard methods)

welcome[bot] commented 3 years ago

Hi 👋 Thanks for reaching out and opening your first issue here! We'll try to come back to you as soon as possible. ❤️ kenobi

DominiqueMakowski commented 3 years ago

Hi @shanelindsay, thanks for chiming in. It's true overall that our RSP processing hasn't received the same amount of love as others signals. So we're really glad to receive suggestions and feedback from experts such as yourself so that together we can make the RSP features top-notch!

The example https://neurokit2.readthedocs.io/en/latest/examples/rrv.html doesn't specify the nature of the respiration measurment (but perhaps should)

by nature of measurement you mean like whether it comes from a respiration belt etc.?

The code which calculates whether a sample is inspiration or expiration assumes that a trough is the start of inhalation, which would the case for breathing belt data

It is true that currently the code implicitly assumes respiration belt data, and we don't have options or methods more adapted to other techniques. To be honest, because we don't use other techniques in our research, our experience and understanding of the specificities are fairly limited (hence the importance of your suggestions). That said, I'm not entirely sure which processing steps in particular could be differentiated based on the recording technique... the peak detection? The filtering/cleaning? the features extraction 🤔 Perhaps we could list like all the techniques measuring "RSP" and for each, the specificities that we need either to specifically support, or that would make them incompatible with some of the existing NK algos.

Traditionally these measures would use the onset of the inhalation as the start of a breath cycle not the start of the exhale

Good catch, we need to address that.

zen-juen commented 3 years ago

Great call! Interesting that the shading of our rsp plot is actually in line with what was raised here, in that the trough is the start of the breath cycle (and hence inhalation onset), although our operationalization of bbi as np.diff(peaks) implies that each peak is the inhalation onset. image

And so if I'm understanding correctly, here the labels of exhalation troughs should be inhalation onsets and inhalation peaks should be exhalation onsets. Adding this to our general roadmap!

shanelindsay commented 3 years ago

Thanks for the response guys. Yeah @zen-juen I went through the code to double check and saw that the graph and other code that does assume troughs as inhale onsets. The reason why I looked into it is that I ran some analysis that showed expiration to shorter than inspiration, which is biologically implausible in normal populations, and just wanted to double check the code. Which checked out for this. I found the probable reason for my issue that was pauses being coded as inhales thus inflating mean inhale length, suggesting to me a need to code pauses (which I will return to below).

One suggestion would be to add a default argument to some of the functions of type = "bb" or "RIP" for breathing belt data, which would make things more explicit, and open the possibility in future of other types of signals.

by nature of measurement you mean like whether it comes from a respiration belt etc.?

Yes. It could also perhaps mention that it assumes a single belt recording and the user would need to combine dual recordings (i.e. thorax and ab) before inputting into neurokit.

That said, I'm not entirely sure which processing steps in particular could be differentiated based on the recording technique... the peak detection? The filtering/cleaning? the features extraction 🤔 Perhaps we could list like all the techniques measuring "RSP" and for each, the specificities that we need either to specifically support, or that would make them incompatible with some of the existing NK algos.

There are a lot of methods for measuring RSP! See https://iopscience.iop.org/article/10.1088/1361-6579/ab299e

For example, khodadad2018 algorithm is based on electrical impedance tomography. While the algorithm seems to work ok, there are some features that are very different compared with breathing belt data, e.g. apneas here:

https://thorax.bmj.com/content/thoraxjnl/72/1/83/F2.large.jpg

I am guessing the main demand for neurokit users is breathing belt data i.e. respiratory inductance plethysmography (RIP). Potentially spirometry and estimation from RR (such as from PPG) might also be of interest (see link above for developments in that area).

For RIP data, the ideal would be to have have algorithms better adapted or validated for this method, compared with khodadad2018 and the other one. I will add a feature suggestion on a separate post about this.

But briefly, the Breathmetrics algorithm is designed for spirometry but easily adapted for RIP, so if that was implemented in neurokit, along with some additional metrics for volume measurement/estimation, then it would have support then for RIP and spirometry (two for the price of one). But supporting pauses could be a headache for existing code.

RIP data also can be quite prone to movement artefacts (say, compared to spirometry) and improved outlier detection would be beneficial here. I will see if I can add another post on that.

shanelindsay commented 3 years ago

Also to add a separate but related point:

In this code https://neurokit2.readthedocs.io/en/latest/_modules/neurokit2/rsp/rsp_findpeaks.html :

_rsp_findpeaks_extrema(rsp_cleaned):
    # Detect zero crossings (note that these are zero crossings in the raw
    # signal, not in its gradient).
    greater = rsp_cleaned > 0

I think it might be based on khodadad2018 ?

For RIP measurements are often in arbitrary units where zero has no meaning. So the raw signal might not have a zero, unless it is z-transformed. This differs from other methods like spirometry and electrical impedance tomography where 0 crossing has a specific meaning e.g. in spirometry downwards zero crossing marks exhale onset. So I haven't thought through the consequences of this code but raising it as a potential issue.

zen-juen commented 3 years ago

Hi @shanelindsay, in the previous PR (ref above) I've changed all instances of "inhalation peaks" and "exhalation troughs", to "peaks (exhalation onsets)" and "troughs (inhalation onsets)" respectively in the docstrings. The explicit labelling of onsets is done only in rsp_plot and I left the namings in the generated signal columns/dictionary info as they are (i.e., RSP_Peaks and RSP_Troughs) since it seems more consistent with our API like that. rsp_rrv() also now correctly uses troughs rather than peaks.

However, I have one question regarding the computing of respiratory rate - as you can see in our code, there are currently two methods to do so: https://github.com/neuropsychology/NeuroKit/blob/c1104386655724a5c624e740abf651c939fc5e48/neurokit2/rsp/rsp_rate.py#L24-L28 The latter is independent of the indices but for the former, does it matter whether periods are computed from peaks or troughs?

zen-juen commented 3 years ago

Regarding the above-mentioned issue about respiratory rate, see https://www.mdpi.com/1424-8220/20/6/1607 who concurs with the use of peaks from BIOPAC belt measurements, but with an additional threshold implemented:

The respiration rate was calculated by counting the peaks from the respiration belt, which were used as the ground truth. The threshold for detecting peaks in the respiration belt data was calculated using Equation (9). The term RGnd(n) is the signal measured from the respiration belt. If the peak was larger than the threshold, it was detected, and the peak was counted only when the interval between the peaks was 2.5 s or more. Threshold = 0.8 × (1/(max(RGnd(n)) − min(RGnd(n)))) (9)

as well as in https://www.mdpi.com/795524

Then, the average breath frequency has been extracted from the respiratory signal acquired through the wearable belt and compared with the one recorded by the spirometer. The average frequency has been extracted from the average breath-to-breath interval, estimated on the basis of the peak-to-peak distance.

and https://www.sciencedirect.com/science/article/pii/S2352648320300015 (also with BIOPAC belt sensor) image

shanelindsay commented 3 years ago

Great to see those changes.

In answer to your question, in one sense it might not matter, since you are going to get similar results either way. But I think it does matter in a couple of ways. In that first paper you cite, they only care about respiration rate, nothing else. And it is primarily about using radar to detect respiration rate, and maybe peaks are better here for that method. The others seem similar in that they are engineering solutions, whereas I would see as neurokit more aimed at scientific solutions, for psychophysiology research.

In neurokit, the respiration analysis emcompasses lots of aspects of the respiration signal.

So it seems a bit weird (and a reviewer might think it weird) in a paper to say we calculated the start of the breath as inspiration onset (for various metrics), but for calculation of respiratory rate we counted exhalation onset.

The code seems to suggest it would be straightforward to change peaks to troughs, so this is what I would advocate for.

On Tue, 14 Sept 2021 at 03:45, Zen @.***> wrote:

Regarding the above-mentioned issue about respiratory rate, see https://www.mdpi.com/1424-8220/20/6/1607 who concurs with the use of peaks from belt measurements, but with an additional threshold implemented:

The respiration rate was calculated by counting the peaks from the respiration belt, which were used as the ground truth. The threshold for detecting peaks in the respiration belt data was calculated using Equation (9). The term RGnd(n) is the signal measured from the respiration belt. If the peak was larger than the threshold, it was detected, and the peak was counted only when the interval between the peaks was 2.5 s or more. Threshold = 0.8 × (1/(max(RGnd(n)) − min(RGnd(n)))) (9)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/neuropsychology/NeuroKit/issues/528#issuecomment-918748791, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAVKQFEIRCUBVRJN4EN4ANDUB2ZNBANCNFSM5DG2X23A . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

zen-juen commented 3 years ago

Yes it's indeed true that many studies pit more modern/rare respiration measurement techniques (like radar) against the respiratory signal as measured by a belt as a reference. Based on a breadth-focused search, I haven't yet found any papers that are purely focused on evaluating the latter.

I also agree that it is concerning if our operationalization of BBI is inconsistent with the way we compute respiratory rate, as you correctly highlighted above. Given that respiratory rate is akin to the number of respiration events in time, it makes sense to compute it based on periods between successive inspirations (inspiration onset-to-inspiration onset) - which this paper also does. What do you think? @DominiqueMakowski

DominiqueMakowski commented 3 years ago

i.e. based on troughs?

zen-juen commented 3 years ago

i.e. based on troughs?

Yes 😄