こんにちは。データサイエンスチームの坂本と申します。
みなさん、データ処理における「外れ値」と聞くと、真っ先に除外すべきものというイメージをお持ちではないでしょうか。実は、必ずしもそうとはいえない場合もあるのです。
TechBlog第八回では、ウェアラブルデータ（※）解析の面白さを知っていただくために、この「外れ値」をめぐるひとつの事例をご紹介したいと思います。
※ウェアラブルデータ……スマートウォッチを代表とするウェアラブル端末で測定された生体情報等のデータ
ウェアラブルデータの分析手法
ウェアラブルデータの分析では、目的に応じてさまざまな手法を使います。
代表的なものには、２つの群の差を取り扱う手法（t検定、分散分析等）や、データとデータの関連性を捉える手法（相関分析、一般化線形モデル等）、反復測定を考慮した手法（マルチレベル分析等）、分類や予測に適した手法（機械学習、深層学習等）、グループ分けに適した手法（クラスター分析等）などがあります。
これらの手法を駆使してデータの切り口を工夫したり、的確な分析手法を選択したり、ときには分析手法自体をオーダーメイドすることで、データからいろいろな情報を抽出することができます。
ウェアラブルデータは非常に豊かな情報源です。同じデータに対しても、さまざまなアイデアと手法を駆使することで、多様な研究に使える・発展する可能性を秘めています。言いかえれば、解析者が自身のアイデアでデータの持つ可能性を広げられる点こそが解析の面白さでもあります。
未知なる可能性を秘めたウェアラブルデータの解析。なかでも今回注目したいのが、先ほども触れた「外れ値」です。
極端に高い心拍数の捉え方
研究やデータ解析において、「外れ値」は除外対象として捉えられるのが一般的です。外れ値が、たとえばセンサーの誤作動といったような「適切に測定されなかったことに由来するノイズ」である場合、それを含めて分析してしまうと、結論に誤りが生じてしまう可能性があるためです。
しかし、一見すると除外対象のように思える「外れ値」も、ノイズとは異なる重要な意味を持っている場合があります。それを本稿では「意味のある外れ値」と表現しています。
たとえば、解析対象のデータが「前方に道路を横切る歩行者が現れてから、運転者がブレーキを踏むまでの時間」であったとします。このようなケースでは、ブレーキを踏むまでの平均時間よりも、極端な踏み遅れ（＝意味のある外れ値）が重要な意味を持つことになります。事故につながるリスクが格段に高い現象だからです。意味のある外れ値がどのような条件の時に生じやすいのか。運転者が強い眠気を感じていたのか、歩行者を視認しにくい特別な交通状況があったのか……がわかれば、対策を立てることができます。
こう考えると、意味のある外れ値を除外せず、研究することの重要性がよくわかると思います。
同じくして、「極端に高い心拍数」にも重要な意味があると考えられます（もちろん、測定ミスの場合を除いて）。その背景には、急激な運動、動機や息切れの発生、心疾患などの病気や、手術後で弱った身体への負担、強度の緊張の発生など、さまざまな理由が想像できます。極端に高い心拍数の発生を定量的に評価することができれば、さまざまな研究への応用が期待できそうです。
では、実際に「定量的に評価する」方法にはどのようなものがあるでしょうか？
指数−ガウス分布
外れ値のような極端に大きな値を除外せず、解析的に考慮するための手法の一つに、「指数−ガウス分布（執筆者訳。一般的にはex-Gaussian Distributionや Exponentially Modified Gaussian Distribution等と呼ばれる）」を用いた分析があります。先ほど例に挙げた「ブレーキを踏むまでの時間」など、反応にかかる所要時間の分析に用いられることがあります。
指数−ガウス分布は、基本的なデータのばらつき部分を表現する正規分布と、外れ値に相当する成分を表現する指数分布とを合成（畳み込み積分により導出）した確率分布です。「山なりで右に裾が長い（重い）分布」という特徴的な形状を持ちます（図1）。
合成元となる正規分布の平均（μ）と分散（または標準偏差σ）、指数分布の期待値（τ、場合により、その逆数のλ）がパラメータであり、特に「外れ値の大きさや出現頻度に関する情報」をパラメータτで捉えることができます。データが指数−ガウス分布にしたがう（近似できる）と仮定し、パラメータ推定を実施することで、これらを定量化することができます。
睡眠中の心拍数データの分析
図2は、特定の期間に測定された、睡眠中の心拍数の度数分布です。2名分の実際の測定データ（※）を可視化しています。
※ 社内メンバーの提供データであり、お預かりした研究データではありません。
※最小値と最大値を破線、中央値を実線で示しています。
横軸は１分あたりの心拍数です。その心拍数が睡眠中の何分間、観測できたかを表すのが縦軸です。グラフからは、Ｂさんのみ1分あたり100以上の心拍数が観測されていることがわかります（赤矢印の箇所）。
Ａさんは、睡眠中の最小心拍数が約50、最大心拍数が約80です。測定データはこの区間において「山なりで、右に裾が長い（重い）形状」で分布していることがわかります。
他方、Ｂさんも「山なりで、右に裾が長い（重い）形状」であることに変わりはありませんが、その裾がAさんよりも長い（重い）です。最大値も180を超えており、ときおり極端に高い心拍数が観測されていることがわかります。
「睡眠中の1分あたりの心拍数が80を超えたものは極端な値（＝外れ値）」として、その割合を求めるなどの方法も考えられます。ただ、個人によって心拍数の平均値などが違うことを考えると、一律に「心拍数が80を超えたもの」のようなしきい値を決めるのではなく、個人差を考慮できる何らかの統計分析手法を適用した方がよいでしょう。
そこで今回は、ＡさんとＢさんの２名について、それぞれのデータを用いて指数−ガウス分布のパラメータを推定します。統計解析環境R（version4.3.3）と、ハミルトニアン・モンテカルロ法（NUTSアルゴリズム）を搭載したソフトウェアSTAN（version 2.32.2）を使用します（iteration = 2000, warmup = 1000, thinning = 4, chains = 4, total samples = 1000）。分析モデルの詳細は、本稿末尾をご覧ください。
分析結果
パラメータ推定の結果、全てのパラメータについて適切に推定できていることが確認できました（Rhat統計量<1.01、実効サンプルサイズ>10%、モンテカルロ標準誤差/標準偏差 < 10%）。
実際の度数分布と事後予測分布を重ね、視覚的にチェックしてみましょう（図3）。
※事後予測分布…既に観測されたデータに基づいて、将来のデータや新しい観測がどのような値をとるかを予測するために使う分布。ここでは先ほどパラメータ推定した指数−ガウス分布のこと。
※破線は測定データの最小値と最大値を示す。
図3では、実際の測定データの度数分布（ヒストグラム）と、曲線で示した事後予測分布（5つのiterationを無作為に抜粋）を重ねて表示しています。指数−ガウス分布によって、睡眠中の心拍数のばらつきかた（分布）を的確に捉えられていると判断してよさそうです。
次に、パラメータの推定結果も見ていきましょう。推定の結果、Ａさんより、Ｂさんの方が、外れ値の頻度や大きさを表すパラメータτの値が高いことが確認されました（δ（τA - τB） = -0.508[-0.887, -0.1234]）。結果を図4に示します。
※点は事後分布の期待値（Expected A Posteriori）を示し、エラーバーは95%信用区間（Credible Interval）を示します。
技術の応用と研究への発展
このように、指数−ガウス分布を用いた解析は、「意味のある外れ値」を定量的に評価する際に有用だと考えられます。
さらにここから、さまざまな分析モデルへと発展させることもできそうです。「年齢が高くなるほど、パラメータτの値が高くなる」ことを仮定した回帰モデルや、「対照群より疾患群の方がパラメータτの値が高い」という個人と集団の階層性を仮定した分析モデル、そして、「手術直後からの3ヶ月間でパラメータτの推定値が徐々に低下した」などを取り扱う時系列モデルなど、やり方次第で多くの研究に貢献できる可能性があります。
もちろん、適用範囲は「極端に高い心拍数」だけではありません。極端に長い睡眠時間の発生や、極端に多い運動量など、さまざまな測定データに対して活用できる可能性があります。
「極端に高い」や「極端に長い」がキーワードになる現象を取り扱う際には、指数−ガウス分布の活用を検討してみてください。ウェアラブルデータを用いた詳細な研究等のご相談は、ぜひ弊社まで。
分析モデルとスクリプト（.stanファイル）
分析モデルを記載したStanスクリプト
// ExGaussian.stan // 入力データ data { int<lower=0> N; // 対象者数(今回はN=2) int<lower=0> r; // 観測したすべての心拍数の数(レコード数) int<lower=1> ID[r]; // 観測レコードに対応した対象者のID番号(A=1, B=2) vector[r] Y; // 観測した1分ごとの心拍数の値 int Const; // 心拍数のスケールを調整する定数(今回はConst = 60) } // データの前処理 transformed data { // 0付近の正の実数で安定した推定をするために、スケールを秒単位に。 vector[r] rescaled_Y; rescaled_Y = Y / Const; } // 推定パラメータ parameters { // 個人ごとのパラメータ vector[N] mu; vector<lower=0>[N] sigma; vector<lower=0>[N] tau; } // モデル model { // 事前分布 mu ~ normal(0, 10); sigma ~ student_t(3, 0, 1); tau ~ gamma(0.01, 0.01); // 尤度関数 for(i in 1:r) { // stanでは指数-ガウス分布の引数にλを使用するため、1/τを所与する rescaled_Y[i] ~ exp_mod_normal(mu[ID[i]], sigma[ID[i]], inv(tau[ID[i]])); } } // 事後予測分布の生成 generated quantities { vector[r] predY; for(i in 1:r) { // 生成した予測分布を元のスケールに戻す predY[i] = Const * exp_mod_normal_rng(mu[ID[i]], sigma[ID[i]], inv(tau[ID[i]])); } }
書いた人：坂本
