【2021-06-15 更新】本記事の内容は時代遅れです。局所解なので、アクセスログを見ながら優先度をつけてアップデートします。

実時間データの分析には色んなやっかいごとがあるように思えます。例えば以下のデータは瞳孔経を300Hzで記録したもので、時間と瞳孔経(mm)がそれぞれstamp列とpupil列にあります。

subject stamp pupil order item condition
P05 1 1.54 3 6 d
P05 5 1.70 3 6 d
P05 8 2.14 3 6 d
P05 11 1.59 3 6 d
P05 15 1.53 3 6 d
P05 18 2.05 3 6 d

ただデータが多すぎるせいで、可視化する時に重くなったり偽陽性がでてしまったりします1。あとはほとんど変化がないためデータとして冗長だったり、機材によってサンプリング周波数が違ったり…といったところです。

もっと前向きなモチベとしては、特定の区間にあるデータを同じ区間とみなし、グループ全体で平均や和を計算したい、なんてこともあるかもしれません。しかしこのままだとどこが一つの区間かがわからないため不可能です。

そういう時に連続値を任意の間隔で区切ってカテゴリーを与える処理をビン分割なんて言ったりします。そこで今回は使いまわせそうなビン分割関数を組んでみます2。

データの取得

まずはライブラリをimportします。ただ以下のライブラリは書きやすくなる、読みやすくなるだけで実際は一つも要りません。

library(dplyr)
library(magrittr)

使うデータは大学のサーバーに置いてあるのでurlで取ってきてloadします3。summaryで内容を確認できます。最初に画像刺激を呈示してから10秒近いデータがとれてますね。この被験者さんの瞳孔系は2mm前後というところみたいです。

# データ読み込み
load(url("http://phiz.c.u-tokyo.ac.jp/~kisiyama/gca2018/pupil.sample.Rdata"))
summary(pupil.sample)
#     subject         stamp           pupil          order   item     condition
#  P05    :2909   Min.   :    1   Min.   :1.28   Min.   :3   6:2909   a:   0   
#  P07    :   0   1st Qu.: 2584   1st Qu.:2.13   1st Qu.:3            b:   0   
#  P08    :   0   Median : 5154   Median :2.19   Median :3            c:   0   
#  P09    :   0   Mean   : 5276   Mean   :2.19   Mean   :3            d:2909   
#  P10    :   0   3rd Qu.: 8017   3rd Qu.:2.27   3rd Qu.:3                     
#  P11    :   0   Max.   :10440   Max.   :3.45   Max.   :3                     
#  (Other):   0

今回のミソはstamp とpupilです。もう少し細かく見るためにheadして最初の20行をみます。

pupil.sample %>% head(20)
subject stamp pupil order item condition
P05 1 1.54 3 6 d
P05 5 1.70 3 6 d
P05 8 2.14 3 6 d
P05 11 1.59 3 6 d
P05 15 1.53 3 6 d
P05 18 2.05 3 6 d
P05 21 1.53 3 6 d
P05 25 1.85 3 6 d
P05 28 1.66 3 6 d
P05 31 1.69 3 6 d
P05 35 1.76 3 6 d
P05 48 1.28 3 6 d
P05 61 2.03 3 6 d

さて本題です。今回の目標はstampを区切ってカテゴリを与えることです。

例えば20msごとに区切りたい場合は以下のようにbin(三列目)を追加したいわけです。20,40,60,80と20ms区切りのカテゴリがそれぞれの行に与えられています。そして以下のようにビン分割するには hist関数が役立ちます。

subject stamp bin pupil order item condition
P05 1 20 1.54 3 6 d
P05 5 20 1.70 3 6 d
P05 8 20 2.14 3 6 d
P05 11 20 1.59 3 6 d
P05 15 20 1.53 3 6 d
P05 18 20 2.05 3 6 d
P05 21 40 1.53 3 6 d
P05 25 40 1.85 3 6 d
P05 28 40 1.66 3 6 d
P05 31 40 1.69 3 6 d
P05 35 40 1.76 3 6 d
P05 48 60 1.28 3 6 d
P05 61 80 2.03 3 6 d

hist関数でビン分割

histはたぶんRを勉強し始める時にもっとも最初のころに触れる関数でヒストグラムをだしてくれます。上のテーブルのstampが0msから60msまでの範囲を20msでビン分割すると下のような図ができます。

これはhistです。

ちなみに上の図は以下で描けます。

bin = 20
pupil.sample %>% filter(0<=stamp,stamp<=60) %$%
    hist(stamp, breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin))

ここで大切なのは一つのビンの中にある個数がそのまま列の数になるという点です。図の0–20は6つ、21–40は5、41–60は1であることが分かりますが、これが上の表のbinと対応しています。

このヒストグラムから上のビン分割した列を作るためにはステップが3つ必要です。

  1. 対応するビンの最大値(20,40,60,80…)を取得
  2. 0–20の個数、21–40の個数とそれぞれのビンにある要素の数を取得
  3. ビンの個数分、対応する最大値を持つ行を作成

histはヒストグラムを描くための情報を保持しているので、それらの情報にアクセスして1と2を実現させます。

bin取得

まずは横軸の情報、つまり何を基準にhistがビン分割したか、を取得できます。ただそのままやると毎度ヒストグラムが生成されるので、plot=FALSEと指定して図を出力しないようになります。

bin = 20
bins =
    pupil.sample %>% filter(0<=stamp,stamp<=60) %$%
    hist(stamp,
         breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin),
         plot=FALSE) %$%
    breaks[-1] 
bins
# [1] 20 40 60

結果は20, 40, 60と20msでビン分割した結果の数字が取れています。あとはそれぞれのビンの頻度を求めます。

bin個数取得

各ビン内のデータが幾つあるか、という情報はcountsに入っています。

bin = 20
n_bins =
    pupil.sample %>% filter(0<=stamp,stamp<=60) %$%
    hist(stamp,
         breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin),
         plot=FALSE) %$%
    counts
n_bins
# [1] 6 5 1

結果は6,5,1なので、0–20までが6つ、20–40までが5つ,そして40–60までが1つという結果が返ってきました。あとはこれを反映したベクトルがあれば先ほどの表を再現できます。

mapply

0–20までが6つ、20–40までが5つ, そして40–60までが1つという結果を反映したベクトルを生成するためにはビンに対してビン個数を対応させたベクトルをつくります。これに関しては一度見た方が分かりやすいと思います。

mapply(rep,bins,n_bins)
# [[1]]
# [1] 20 20 20 20 20 20
# 
# [[2]]
# [1] 40 40 40 40 40
# 
# [[3]]
# [1] 60

これはrepという関数をbinsとn_binsの各要素のそれぞれに適応させています。つまりrep(bins[1],n_bins[1]),rep(bins[2],n_bins[2]),...というのを全ての要素分おこないます。これをunlistするとベクトルになるので、それを列に加えます。

最終的な関数と利用方法

以上をまとめた関数が以下になります。

# take a list of timestamp and bin, return the same number of rows with binned stamps 
binner = function(stamp, bin){
    bins  = hist(stamp,
                 breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin),
                 plot=FALSE) %$%
            breaks[-1] 
    n_bins= hist(stamp,
                 breaks=seq(0,ceiling(max(stamp)/bin)*bin, bin),
                 plot=FALSE) %$%
            counts
    mapply(rep,bins,n_bins) %>% unlist %>%
    return
}

後はこの関数をどうやって適用するかですね。subject, itemにgroup_byしてあげると、続く操作をsubjectとitemごとに行えます。その操作というのがビン分割なので、stampとbinを引数に与えてビン分割したベクトルを返します。そのままcbindするともともとのデータフレームと順序が違うので(subject, itemをベースにした順序になっている)、もとのデータフレームにarrangeを当ててあげています。

pupil.sample %>% group_by(subject,item) %>%
    do(data.frame(bins=binner(.$stamp,bin=20))) %$%
    cbind(bins,arrange(pupil.sample,subject,item)) %>% head(13)

する以下のような結果になります。ちゃんとstampに対応したbinが作られていますね。

bin subject stamp pupil order item condition
20 P05 1 1.54 3 6 d
20 P05 5 1.70 3 6 d
20 P05 8 2.14 3 6 d
20 P05 11 1.59 3 6 d
20 P05 15 1.53 3 6 d
20 P05 18 2.05 3 6 d
40 P05 21 1.53 3 6 d
40 P05 25 1.85 3 6 d
40 P05 28 1.66 3 6 d
40 P05 31 1.69 3 6 d
40 P05 35 1.76 3 6 d
60 P05 48 1.28 3 6 d
80 P05 61 2.03 3 6 d

まとめ

あとはsubject,item,condition,binあたりでgroup_byしてsummariseすれば平均や標準偏差を求めるなど、好きな操作ができます。

参考

  • https://note.nkmk.me/python-pandas-cut-qcut-binning/
  • https://www.ibm.com/support/knowledgecenter/ja/SSLVMB_24.0.0/spss/base/idh_bander_gating.html#idh_bander_gating

  1. ここらへんの議論はMirman(2014)のp.19あたりを参照。 ↩

  2. 実際にビン分割する際はデータの性質（どれくらい時間センシティブなデータなのか）や研究分野（心理言語学なのか教育学なのか、はたまた金融工学なのか）、研究課題によります。あとpythonのpandasなんかだとcut関数があるのですが、Rで該当する関数が見つけられず、また先輩から頂いたスクリプトもreshapeとか駆使しており、初見の時難しいのに加えて一部の関数がdeprecateしており他の人に継承できそうになかった次第です。 ↩

  3. 余談ですがデータ共有はバイナリー形式がおすすめです。 tsvとかcsvとかだとOSによりBOMの問題があります。 ↩