threshold_filter can generate negative RR values

[noah10]

The cubic spline interpolation used in the threshold_filter function can generate negative values. These then cause the interpolated data set to fail validation. Here's some sample code:

from hrv.rri import RRi
from hrv.filters import threshold_filter

rri_vals = [ 795, 832, 827, 897, 916, 902, 870, 827, 819, 757, 749, 725,
  727, 743, 725, 722, 719, 773, 784, 832, 908, 897, 859, 835,
  859, 824, 805, 754, 760, 765, 711, 849, 528, 762, 671, 725,
 1309, 560, 962, 892, 876, 900, 838, 873, 784, 830, 584, 579,
  803, 611, 620, 703, 768, 795, 765, 644, 714, 722, 717, 706,
  754, 714, 695, 730, 708, 722, 760, 849, 867, 878, 859, 857,
  881, 892, 819, 789, 835, 876, 770, 838, 862, 884, 876, 840,
  760, 822, 892, 935, 932, 824, 830, 795, 792, 873, 916, 983,
 1013, 989, 962, 959, 937, 959, 908, 897, 889, 816, 768, 752,
  741, 692, 620, 789, 611, 671, 663, 501, 797, 1517, 708, 601,
  894, 708, 730, 711, 733, 717, 719, 711, 671, 636, 641, 649,
  595, 636, 620, 628, 703, 725, 808, 792, 851, 867, 881, 927,
  894, 884, 873, 892, 916, 862, 870, 908, 967, 964, 943, 851,
  822, 843, 916, 989, 1010, 1010, 878, 797, 913, 784, 884, 951,
 1018, 1088, 1037, 1029, 929, 876, 805, 738, 706, 652, 655, 692,
  665, 606, 676, 673, 679, 757, 725, 1983, 2207, 671, 687, 932,
  590, 776, 838, 784, 808, 795, 808, 690, 819, 746, 795, 835,
  981, 929, 956, 916, 865, 851, 770, 735, 768, 706, 762, 698,
  727, 773, 795, 854, 994, 935, 889, 814, 733, 717, 725, 1072,
 1018, 784, 660, 784, 620, 770, 743, 768, 811, 897, 892, 840,
  811, 827, 816, 884, 886, 876, 752, 851, 859, 789, 768, 760,
  800, 811, 921, 851, 819, 859, 865, 787, 849, 876, 827, 773,
  717, 760, 762, 854, 897, 997, 991, 919, 886, 843, 784, 757,
  760, 735, 730, 1355, 620, 614, 665, 676, 668, 644, 668, 733,
  770, 700, 749, 867, 897, 805, 789, 749, 703, 719, 770, 628,
  768, 412, 932, 652, 749, 760, 773, 905, 865, 876, 711, 752,
  722, 795, 803, 781, 765, 730, 711, 738, 725, 746, 811, 859,
  881, 800, 822, 584, 956, 595, 800, 749, 1506]

rri = RRi(rri_vals)
corrected_rr_vals = threshold_filter(rri, threshold='low', local_median_size=5)

And here's the stack trace that it generates:

Traceback (most recent call last):
  File "hrv-bug2.py", line 34, in <module>
    corrected_rr_vals = threshold_filter(rri, threshold='low', local_median_size=5)
  File "/Users/home/.virtualenvs/pyhrv-kubios-A4er0XLO/lib/python3.8/site-packages/hrv/filters.py", line 228, in threshold_filter
    return RRi(cubic_spline(rri_time), rri_time)
  File "/Users/home/.virtualenvs/pyhrv-kubios-A4er0XLO/lib/python3.8/site-packages/hrv/rri.py", line 57, in __init__
    self.__rri = _validate_rri(rri)
  File "/Users/home/.virtualenvs/pyhrv-kubios-A4er0XLO/lib/python3.8/site-packages/hrv/rri.py", line 464, in _validate_rri
    raise ValueError("rri series can only have positive values")
ValueError: rri series can only have positive values

[rhenanbartels]

HI @noah10, how are you doing? Sorry for the late reply, I'm returning from vacation.

Thank you very much for reporting the bugs. I will check the code for both issues (issue #22 and issue #23) in the following days.

Cheers!

[noah10]

No worries about the delay and thanks for looking into them. Let me know if you need any more info.

[rhenanbartels]

Hi @noah10, I dug into the code and figured out that the CubicSpline function is generating the negative values. I've tried other interpolation methods, but they also produce negative values.

I think because the value that needs to be interpolated is in the edge of the signal (last value), the interpolation method is not able to properly identify a reasonable surrogate for it.

I am still thinking about the best solution for this situation. I am considering creating a class RRiWithNegativeValues which allows negative values and prints a warning message about them. The problem with this solution is: accepting values that were mistakenly calculated will affect all results. So maybe just keeping throwing the exception would be a better approach. In this scenario, the user will need to check the signal and re-prepare it.

I've noticed that removing the last value of rri_vals array will do the trick and stop the exception.

Any thoughts?

Cheers,
Rhenan

[noah10]

Yes, I had noticed that when I first found this issue. I did a little bit of digging to learn about positive-only interpolation and came up with the following techniques:

Monotone cubic interpolation
Interpolating the log of the data and returning the exp of the result
pchip interpolation

Unfortunately I don't have enough signal processing experience to really know the pros and cons of each of these approaches in the context of HRV analysis.

I don't think that allowing negative RRi values is a good approach, though. They're fundamentally wrong - nobody can have -32 ms between heartbeats. And suggesting that people massage their data in order to work around a problem created by trying to massage their data also doesn't seem right. :-) (Removing the last value isn't a universal fix, of course.) In the absence of any other solution I'd suggest picking one of the algorithms above and using that.