Rで心拍数からストレスを可視化する
経緯
以前俺のドキドキを可視化するという話をした。
俺のドキドキは大体心理的ストレスからくるのであり、日々の重圧にさらされながら俺の心臓はそれでも負けじと自律神経と手を取り合って俺の恒常性を保つべく頑張り続けている。
ということで心拍数からストレスを可視化できるんだろうなーとか思って調べたら当然のごとく商用サービスがある。
商用サービスはたとえばこれ。
商用になってるということはなにかしらインデックスを算出する方法があるということでありさらに調べるとまとまった文献が見つかる。
ここなんて仕組みを解説した上にmatlabのコードまで掲載されている。最高。
ということでRでストレスを可視化する。
概観
心拍変動を周波数解析することで、各周波数の波に分解できる。
どの周波数の波が多いかでストレス状況を可視化できるという理屈。
もう少し具体的にいうと、まずこの波の分解は「パワースペクトル密度の算出」と呼ばれている。
分解された結果出てくる波のうち高周波の波(0.15-0.40Hz)をHF、低周波の波(0.05-0.15Hz)をLFという。
で、LFの積算値がHFの積算値に比べて大きい、つまり交感神経優位だと、ストレスがかかっている、という対応関係がある。
LF積算値のHF積算値に対する相対的な大きさはLH比と呼ばれており、これは単に前者を後者で割って算出する。
詳しい話は上で紹介したサイトを参照いただきたい。
ということで実質Rでやることはパワースペクトル密度の算出だけで、あとは既定の周波数ごとに積算してLFをHFで割ればよい。
心拍数データを読み込む
github上にデータを置いているのでそちらを読み込む。
httpsはRの組み込みではサポートされていないのでreadrパッケージを用いた。ただ、readrで読み込むとtbl形式になってしまい、これはdata.frameと互換性があるとうたっているものの時々変な挙動を示すのでdata.frameのみの状態に戻している。
library("readr") f <- "https://raw.githubusercontent.com/dichika/mydata/master/heartrate.csv" d <- read_csv(f,col_types = "cd") d <- as.data.frame(d) d$time <- as.POSIXct(d$time) library("ggplot2") ggplot(d, aes(x=time, y=heart_rate)) + geom_line()
可視化するとこのようになる。心拍数は５秒単位で計測している。
12時過ぎたあたりから心拍数がぐいぐい上がっているがこれは昼ご飯を食べた後移動したからである。
また18時付近で大きめのスパイクがあるが、これは部下に対して腹を立てていた時間帯に対応している。あれは大人気なかった。
LH比を算出する。
5分単位でLH比を算出する。
ローリング処理なのでもっと速い書き方/パッケージはあるはず。
パワースペクトル密度の算出についてはいくつか方法があり、ここではRの組み込み関数 spectrum で可能なAR/ピリオドグラムによる算出法を用いた。
# LH比を抽出する関数 getLHratio <- function(data, m="ar"){ res_spec <- spectrum(data, method = m) lf <- sum(res_spec$spec[res_spec$freq>=0.05 & res_spec$freq<=0.15]) hf <- sum(res_spec$spec[res_spec$freq>=0.15 & res_spec$freq<=0.40]) res <- lf/hf } # 5秒間隔で計測しているので5分単位で算出する場合は300秒/5=60個ずつの枠をずらしながらLH比を算出する window <- 60 ix <- nrow(d)-window # ARを用いて算出する場合 res_lh_ar <- numeric(ix) for(i in seq(ix)){ res_lh_ar[i] <- try(getLHratio(d$heart_rate[seq(i,i+window-1)]), TRUE) } d$lh_ar <- c(res_lh_ar, rep(NA, window)) # ピリオドグラムを用いて算出する場合 res_lh_per <- numeric(ix) for(i in seq(ix)){ res_lh_per[i] <- try(getLHratio(d$heart_rate[seq(i,i+window-1)], m="pgram"), TRUE) } d$lh_per <- c(res_lh_per, rep(NA, window))
心拍数とLH比を併せて可視化
上記で求めたLH比を心拍数と併せて可視化すると以下のようになる。
ARよりピリオドグラムを用いた方が外れ値（この場合は持続時間の短い心拍数のスパイク）に強いような印象がある。
AR、ピリオドグラムともにLH比が５ぐらいまでは正常範囲でストレスがかかるとそれ以上になる様子がわかる。
# 心拍数と併せて可視化 library("tidyr") d_g <- gather(d, var, val, -time) d_g$val <- as.numeric(d_g$val) ggplot(d_g, aes(x=time, y=val)) + geom_line() + facet_grid(var~., scales="free_y")
