[예제] Fitbit 심박 데이터로 Detrended flucatuation analysis (DFA) 해보기 with python

 

 

 데이터에 트렌드가 올라타버리는 non-stationary data에서 Long-memory 효과, 그러니까 multi-scale에서의 작동 메커니즘을 분석하기에 유용한 방법으로써 Detrended fluctuation analysis(DFA)가 있다. 관련해서 개념을 정리하는 글을 몇일 전에 썼는데, 마침 1년 가까이 모아온 심박 데이터가 있어 DFA를 직접 해봤다.

 

심박과 관련해서 DFA를 적용한 분석은 다음 논문이 원류인 것으로 보인다. 

Peng, C. K., Havlin, S., Stanley, H. E., & Goldberger, A. L. (1995). Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos: an interdisciplinary journal of nonlinear science, 5(1), 82-87.

 

문제는, 심박수(beat per minute)가 아닌, Heart rate variability(HRV)를 분석했다는 것. HRV는 한 심박에서 다음 심박까지 걸린 시간을 말하는데, 대략 bpm과 반비례 관계에 있다. 게다가 multiscale을 정의하기 위해 trend를 제거하는 window의 기준인 n도 시간이 아니라 beat 수 기준이다. 그러니까 x축이 beat, y축이 HRV인 데이터에 대한 DFA인 것.

 

Heart rate를 bpm으로 측정한 결과로도 DFA를 분석한 연구가 있는지 찾아봤는데 이건 따로 없었다. 그래서 먼저, DFA의 scale을 이미 알고 있는 시그널을 대상으로 bpm값으로도 DFA가 비슷한 결과를 낼 지 확인해 봤다. (DFA에서 white noise은 0.5, pink noise는 1, brown noise는 1.5의 scale을 보임)

실제 DFA를 위한 코드만 보려면 1은 건너 뛰고 2로.

 

 

1. HRV 데이터 만들기 - Brown, Pink, White noise

 

각 Noise를 생성하기 위해 사용한 코드는 stackoverflow에서 JeffreyLai의 답변을 참조했다.

import numpy as np

def noise_psd(N, psd = lambda f: 1):
        X_white = np.fft.rfft(np.random.randn(N));
        S = psd(np.fft.rfftfreq(N))
        # Normalize S
        S = S / np.sqrt(np.mean(S**2))
        X_shaped = X_white * S;
        return np.fft.irfft(X_shaped);

def PSDGenerator(f):
    return lambda N: noise_psd(N, f)

@PSDGenerator
def white_noise(f):
    return 1;

@PSDGenerator
def brownian_noise(f):
    return 1/np.where(f == 0, float('inf'), f)

@PSDGenerator
def pink_noise(f):
    return 1/np.where(f == 0, float('inf'), np.sqrt(f))

각 noise의 특징인 frequency 프로파일을 이용해 data를 생성하는 방법이다. 이 함수를 이용해서 대략 1분에 60번의 맥박을 상정, 1000ms을 중심으로 오르락 내리락 하는 가짜 심박 데이터를 만들었다.

Brown_HRV = brownian_noise(86400*5) * 300 + 1000
Pink_HRV = pink_noise(86400*5) * 300 + 1000
White_HRV = white_noise(86400*5) * 300 + 1000

데이터 길이는 대략 5일이 되도록 잡았다. 샘플을 300개만 추려보면 대략 이런 느낌.

 

2. Detrended fluctuation analysis - fathon

DFA는 fathon이라는 package를 이용했다. 코드는 다음과 같다. 계산에 시간이 좀 걸리므로 주의 (2~3분 정도). 

 

import fathon
from fathon import fathonUtils as fu

def CAL_DFS(Data):
    # zero-mean cumulative sum
    CS = fu.toAggregated(Data)

    # initialize dfa object
    pydfa = fathon.DFA(CS)
    wins = fu.linRangeByStep(5, 600)
    n, F = pydfa.computeFlucVec(wins, revSeg=True, polOrd=2)
    H, H_intercept = pydfa.fitFlucVec()
    return H, n, F


fig = plt.figure()
H, n, F = CAL_DFS(Brown_HRV)
plt.scatter(n, F, label='Brown ' + str(round(H, 2)))
H, n, F = CAL_DFS(Pink_HRV)
plt.scatter(n, F, label='Pink ' + str(round(H, 2)))
H, n, F = CAL_DFS(White_HRV)
plt.scatter(n, F, label='White ' + str(round(H, 2)))
plt.xscale('log')
plt.yscale('log')
plt.grid(True, which='both')
plt.xlabel('n')
plt.ylabel('F')
plt.legend()
plt.show()

utility인 'toAggregated' 함수는 데이터를 Cumulative sum 해주고, 평균값을 빼주는 함수고, 'linRangebyStep'은 분석할 window size를 배열로 뽑아주는 함수다. 여기서는 조사할 window size(n)를 5부터 600까지로 정했다.

 

window size n에 따른 fluctuation을 계산해주는 함수인 computeFlucVec의 인자로는 사용할 window size의 배열인 wins, window를 데이터 앞쪽+뒤쪽에서부터 잡아서 계산할지 정하는 revSeg, trend를 제거하기 위해 사용할 다항회귀의 차수를 결정하는 polOrd가 있다. 연구에 따라 다르긴 한데 대체로 2차 다항식 회귀를 이용하는 것으로 보인다.

 

위를 실행하면 아래 figure가 출력된다.

각 시그널의 DFA scale이 1.5, 0.99, 0.51로 예상했던 값인 1.5, 1.0, 0.5와 근접하게 나온 것을 확인할 수 있다.

 

 

3. HRV → bpm 변환값의 DFA scale

이제, HRV를 bpm(beat per second)로 바꾼 데이터에서도 유효한지 확인해 보기 위한 코드를 짜보았다.

 

def CAL_BPM_from_HRV(data):
    BT = np.cumsum(data) / 1000
    BPM = []
    while len(BT) != 0:
        target = BT[BT < 60]
        BPM += [len(target)]
        BT = BT - 60
        BT = BT[BT >= 0]
    return np.array(BPM[:-1])


Brown_BPM = CAL_BPM_from_HRV(Brown_HRV)
Pink_BPM = CAL_BPM_from_HRV(Pink_HRV)
White_BPM = CAL_BPM_from_HRV(White_HRV)

 

HRV를 Cumulative sum한 값이 심박이 찍힌 시계열이고, HRV단위가 ms니까 1000으로 나누어 초단위로 바꿔줬다. 그래서 60초 내에 있는 데이터 수를 카운트 해서 분당 bpm으로 바꿨다. (마지막은 bpm은 분단위로 딱 안떨어지는 경우가 발생하니까 버렸다)

 

이렇게 구한 bpm을 위와 같은 방식으로 플롯해보면 

Brown noise의 감도는 좀 떨어지지만, 대략 비슷한 scale로 나온다. Brown이 대략 80 부근에서 scale이 꺾이는데, 이를 기준으로 나눠서 분석해봤다.

fig = plt.figure()
data = Brown_BPM
# zero-mean cumulative sum
CS = fu.toAggregated(data)
# initialize dfa object
pydfa = fathon.DFA(CS)
wins = fu.linRangeByStep(5, 80)
n, F = pydfa.computeFlucVec(wins, revSeg=True, polOrd=2)
H, H_intercept = pydfa.fitFlucVec()
plt.scatter(n, F, label='Brown low ' + str(round(H, 2)))
# initialize dfa object
pydfa = fathon.DFA(CS)
wins = fu.linRangeByStep(80, 600)
n, F = pydfa.computeFlucVec(wins, revSeg=True, polOrd=2)
H, H_intercept = pydfa.fitFlucVec()
plt.scatter(n, F, label='Brown high ' + str(round(H, 2)))
plt.xscale('log')
plt.yscale('log')
plt.grid(True, which='both')
plt.xlabel('n')
plt.ylabel('F')
plt.legend()
plt.show()

대략 80분을 기준으로 그보다 작은 window에서는 scale이 1.5부근으로 나오다가, 그보다 큰 window에서 1.17로 떨어진다. 이건 넓은 window에서 fluctuation이 HRV에 비해 과소평가된다는 뜻인데, 시간 기준으로 window가 넓어지다보면 빠른 심박 구간에서는 더 많이 beat가 찍히고, 느린 심박 구간에서는 더 적은 beat가 찍히기 때문인 것으로 보인다. 그러니까 시간 기준의 fluctuation에서는 상대적으로 빠른 심박(짧은 HRV)에 의한 fluctuation이 과소평가 되고, 느린 심박(긴 HRV)에 의한 fluctuation이 과대평가 되는데, 이게 window가 넓을수록 + 더 long-memory를 지닌 signal일수록 왜곡이 심해지는 것으로 보인다.

 

종합해보면, BPM값에서 포함되는 beat 수의 변동폭이 커지는 높은 window구간에서 multi-scale이 나타난다면 꺾이기 전인 아래 구간을 보는 편이 좋을 것 같다.

 

 

4. Fitbit으로 측정한 심박 데이터의 DFA scale

 

예전 사용하던 Fitbit에서 뽑아놓은 2017년도 심박 데이터를 사용했다. 샤워할 때를 빼고 거의 300일간 계속 차고 있었으니, 꽤 쓸만한 데이터다. 데이터는 일주일 씩 끊어서 각각 DFA scale을 핏팅했다.

 

* fitbit data 받는 방법
(일본어 사이트인데 구글 사이트 번역으로도 볼만하다 + 예전에 참고한 방법이라 지금은 조금 다를지도)

https://oku.edu.mie-u.ac.jp/~okumura/python/fitbit.html

https://qiita.com/osapiii/items/1e53531193ab73f6388b
https://qiita.com/fujit33/items/2af7c4afdb4e07601def

이렇게 구한 Scale 값을 평균내보면 1.16이 나온다. HRV의 경우 건강한 사람에게서 대략 1.0이 나온다고 하는데, 조금 높은 값이다. 어떤 주간은 높은 n 영역에서 scale이 꺾이는 것처럼 보이기도 하는데, 눈대중으로 60분 이상/이하로 나뉘는 것으로 보인다.

 

그래서 60을 기준으로 Low, High 영역의 alpha를 값을 나눠서 시간에 따른 변화를 살펴봤는데, 변동의 정도는 조금씩 달랐지만, 내려갔다 올라갔다 하는 트렌드는 비슷한 경향을 보였다.

흥미로운 건, 동일인물임에도 alpha값에 변동이 꽤 있다는 점이다. alpha값은 건강한 사람에게서 robust하게 1이 나오고(fractality), 심장질환이 있는 사람에게서 깨진 패턴이 발견된다는 게 선행연구들의 baseline인 걸 보면 생각해볼 수 있는 이유로는 내가 건강하지 않거나, 데이터가 역시나 더러웠다, bpm이 아니라 HRV를 봐야 제대로다, 이 중 하나일 것 같다. 또, 운동에 의해 alpha값이 변하는 경우가 있다거나, 원래 단일 scale이 아니라는 연구들도 후속연구로 나오고 있기에, 그러한 가능성도 고려할 필요도 있는 것 같다.

 

그럼에도 불구하고, alpha값 변동이 완전히 noisy하지 않고 일종의 trend가 보인다는 점은 정황적으로 흥미롭다. 2017년 4월까지는 회사를 다니다가 5월에 관뒀다. 이 때 alpha값이 떨어지는데 아마도 퇴사가 심장에 더 좋았을지도 모르겠다. 그러다가 7월 8월에는 이상한 패턴을 보이는데, 이 때 대학원 들어가겠다고 논문 읽고, 연구계획 짜고 밤도 새고 그랬던 기억이 난다. 9월 이후로는 대학원 생활이 시작되는데, 역시나 밤을 새는 일도 있었지만 스케쥴에 치인다기 보다는 자발적인 밤샘이 많았다. 그러니까 몰입해서 밤새고 다음날 푹 잔다거나 하는 패턴이었다는 점에서 7-8월과는 조금 다르다.

 

뭐 이것만 가지고 좋다 나쁘다는 따질 수 없겠지만, 꽤 흥미로운 index인 것 같다. fitbit이 살아있었다면 매일 체킹해보는 것도 재미있었을 것 같다. 여튼.. 가볍게 해본다는 게 꽤 길어졌다.