# パッケージ読み込み・環境設定
require(conflicted)
conflict_prefer("lag", "dplyr", quiet = T)
conflict_prefer("filter", "dplyr", quiet = T)
conflict_prefer("select", "dplyr", quiet = T)
conflict_prefer("slide", "slider", quiet = T)
conflict_prefer("slide_dbl", "slider", quiet = T)
conflict_prefer("slide2_dbl", "slider", quiet = T)
require(profvis)
require(microbenchmark)
require(compareDF)
require(skimr)
require(tidyverse)
require(ggthemes)
require(lubridate)
require(zoo)
require(purrr)
require(RcppRoll)
require(tsibble)
require(slider)
require(tibbletime)
require(runner)
require(rsample)
theme_set(
  theme_bw(base_size = 14, base_family = "Noto Sans CJK JP") +
    theme(axis.text.x = element_text(angle = - 45), axis.title.x = element_blank(),
          axis.title.y = element_blank(), axis.text.y = element_blank(),
          legend.title = element_blank(), legend.position = "bottom"
          )
  )

この記事を読んでわかること (概要)

時系列データの集計で使う移動平均など, どうも用語が統一されてないがわりとよく使う「移動集計」をdplyr時代の作法でスマートにやる方法を紹介する. 単なるサンプルプログラムの紹介で終わらないよう, 現実にありそうな用例を用意したので, 時系列データ特有の前処理についてもいくらか触れている.

結論から言うと, RcppRoll が一番速いが複雑な処理ならまだあまり知られていない(?) slider か runner パッケージを使うのがオススメ.

初めに/能書き

最近 R-wakalang でローリング集計のやり方について質問が出た. 決してマイナーな集計処理ではないが, 意外と細かく説明している資料は少ない. 時系列データを扱うパッケージには zoo が昔からあるが, 今や dplyr ひいては tidyverse が広く普及しており, データフレームへの処理をパイプ演算子でつないでワンライナーで書いていくイマドキのやり方に対して zoo の構文は古臭く感じてしまう¹. 加えて, 最近 dplyr ver. 1.0.0 が登場し, 新たな使い方の可能性が見えてきた². そこで, dplyr/tidyverseと親和性の高いイマドキのローリング集計パッケージについて調べた結果をまとめる.

補足: そもそもdplyrはじめtidyverseを知らないと言う人は, 以下を参考に

自分が本当に初心者だと思う人向け: 『データハンドリング入門』
dplyr 『dplyr — 高速data.frame処理』
tidyr 『htidyr — シンプルなデータ変形ツール』『tidyr 1.0.0の新機能 pivot_*() / tidyr-pivot』
purrr 『purrr — ループ処理やapply系関数の決定版』『{purrr} mapを導入しよう。 - Qiita』
tibble 『Hadley Ecosystem 2016』

dplyr/tidyr/purrr の3つを組み合わせた総合的な解説はあまりない.

ただし上記は dplyr 1.0.0 の大きな更新を反映していない, または詳しく書いてないので, 以下も参照すると良い

時系列データ処理の基本的なこと

このセクションでわかること

dplyr にはラグ関数や順位を付ける関数が充実している
ウィンドウ関数でデータを取り出す範囲をスライドさせる
ローリング集計はウィンドウ関数で取り出したデータの部分を集計すること
dplyr の関数だけでローリング集計するのは大変

`dplyr` の時系列データ処理

移動集計の前に, dplyr で用意されている時系列データ処理に使える関数をいくつか紹介しておく.

そもそもdplyr自体に時系列データの処理を意識した関数がいくつかある. 公式マニュアルでも “Window functions • dplyr” というページが設けられている. 後述するようにこれはlag()/lead()以外windows関数とはちょっと違うのではないか, という気がするが, 時系列データ処理にも使える話なのでいちおう言及しておく. 同じ列のいくつか前/後の値を参照するには lead(), lag() が役に立つ. 売上の前年比とか*成長率とか, もラグを利用した集計処理であり, 複雑高度なデータ分析をしなくても必要な場面が存在することがわかる. 実際, グループ化とラグ集計の合わせ技は R-wakalang でもわりとよく訊かれる質問の1つだと思う.

例えば原系列 origin に対して, 前日比 day_over_day と前日からの成長率 growth は以下のように簡単に計算できる.

set.seed(42)
d <- tibble(date = seq(ymd("2019-12-26"), ymd("2020-01-10"), by = "1 day")) %>%
  mutate(origin = rnorm(n = n()))
d %>% mutate(
  day_over_day = origin / lag(origin),
  growth = day_over_day - 1
  ) %>% head()

date	origin	day_over_day	growth
2019-12-26	1.3709584	NA	NA
2019-12-27	-0.5646982	-0.4119003	-1.4119003
2019-12-28	0.3631284	-0.6430487	-1.6430487
2019-12-29	0.6328626	1.7428066	0.7428066
2019-12-30	0.4042683	0.6387932	-0.3612068
2019-12-31	-0.1061245	-0.2625101	-1.2625101

また, 累積計算もいくつか用意されている. 累積和の cumsum(), 累積平均 (チェザロ平均) の cummean(), 総乗 cumprod(), 累積最小値 cummin(), 累積最大値 cummax() はRに最初から用意されている基本関数だが, dplyr にはさらにいくつか便利なものがある. 冒頭のページで紹介されている関数を挙げると,

cumall(): 累積的な all 条件演算, つまり一度でも条件を満たさない値が現れると以降を全て FALSE とする
cumany(): 累積的な any 条件演算, つまり一度でも条件を満たした値が現れると以降を全て TRUE とする

また, 順位計算に使う関数も, 時系列データ処理で役に立つことが多い. dplyr には以下のような関数が用意されている

percent_rank(): 対応するパーセンタイルを出力する
cum_dist(): 対応する累積パーセンタイルを出力する
row_number(): 行番号を振る. base::rank() でタイの処理を"first"にした場合と同じ
min_rank(): base::rank() でタイの処理を "min" にした場合と同じで, 順位を下に詰める.
dense_rank(): min_rank() の順位のギャップを詰めるバージョン
ntile(): データをいくつかに分割し, 分割したタイル単位で順位を決める

がある(私は教科書やネット上に書いてあることをコピペしてまたネット上に公開するのが嫌いなので, 詳しい解説は上記ページやヘルプを見て欲しい. 最近は機械翻訳でも大意をつかめることが多いし.)

移動集計 (ローリング集計) とは

情報が少ないのはそもそも「移動集計」という言葉が定着していないからではないか, とふと思ったので書いておく. この手の処理は英語でも rolling aggregation という他に running — とか, moving — とか呼ぶのでではっきりしない. 日本語でもそれを反映して「ローリング集計」「ウィンドウ集計」などとも呼ばれる例があるが, これといって決まったものがないように思える. しかし, 時系列データの集計ではそこそこ使われる処理である. 今回は以降から「ローリング集計」で統一することにする.

まず, 集計ではないローリング集計処理とは, データの取得期間を少しづつづらして取り出し, 処理することである. これはスライシング(slicing)またはスライディング(sliding)とも呼ばれるし, データの一部だけが見えるように切り取る様子から「窓かけ」(windowing)と表現したり, そのような関数のことを 窓(ウィンドウ)関数と呼んだり³ ⁴, スライディングウィンドウという表現も見られる. つまり, ローリング集計とはウィンドウ関数が次々と取り出したデータの部分に対してそれぞれ集計をすることである.

例えば2020/1/1 から 2020/1/10 のデータに対して, 1週間のローリング集計をするならば, ある日付に対して取り出す期間は図のようになる⁵.

1週間ウィンドウ幅のイメージ

このように互いに重なった期間でグループ集計するのがローリング集計で, 集計とはいうものの, 元のデータの1件1件に対して計算するため, 集計後のデータの件数は必ずしも減らない. 集計期間が少しづつスライドするという点では, lead()/lag() も広義のローリング集計のためのウィンドウ関数である. 一方で, 累積計算は過去の全ての値を参照するため, ローリング集計を表現できるとは限らない.

典型的なローリング集計は, 各時点の過去数日間平均をとる移動平均(moveing average) で, 7時点前の移動平均を数式で表すなら以下のようになる.

\[ y_t = \frac{1}{7}\sum_{s = t-6}^t y_s \]

適当に作った2019/12/26-2020/1/10までの日次時系列データyとそれに対して1週間移動平均を計算して, 先ほどの図と同じ期間でグラフを描き比較すると, 以下のようになる. 1/1にならないと過去の値が7時点揃わないので, それ以前の期間は省略している.

典型的なローリング集計: 1週間移動平均

この「期間が重複しているグループ集計」というのが厄介で, 例えばSQLでやろうとすると, 単純にGROUP BYだけでは集計できない. 実装ごとにローリング集計のために機能が追加されていることも多いが, 一方でどのSQLでも動くように書くとなると(そうせざるをえない状況があるのかは知らないが…), 期間が増えるたびに条件をずらしたjoinを追加するとか, そういうとても嫌なクエリになってしまう. dplyrでも同様で, dplyr の関数だけで書こうとするなら同様にいくつもjoinするか, 以下のように lag()を使うことになるだろう⁶.

d %>% mutate(
  lag1 = lag(origin),
  lag2 = lag(origin, 2),
  lag3 = lag(origin, 3),
  lag4 = lag(origin, 4),
  lag5 = lag(origin, 5),
  lag6 = lag(origin, 6)
  ) %>% 
  rowwise %>% mutate(ma7 = mean(c(origin, lag1, lag2, lag3, lag4, lag5, lag6))) %>% ungroup %>%
  filter(date >= ymd("2020-01-01"))

date	origin	lag1	lag2	lag3	lag4	lag5	lag6	ma7
2020-01-01	1.5115220	-0.1061245	0.4042683	0.6328626	0.3631284	-0.5646982	1.3709584	0.5159882
2020-01-02	-0.0946590	1.5115220	-0.1061245	0.4042683	0.6328626	0.3631284	-0.5646982	0.3066142
2020-01-03	2.0184237	-0.0946590	1.5115220	-0.1061245	0.4042683	0.6328626	0.3631284	0.6756316
2020-01-04	-0.0627141	2.0184237	-0.0946590	1.5115220	-0.1061245	0.4042683	0.6328626	0.6147970
2020-01-05	1.3048697	-0.0627141	2.0184237	-0.0946590	1.5115220	-0.1061245	0.4042683	0.7107980
2020-01-06	2.2866454	1.3048697	-0.0627141	2.0184237	-0.0946590	1.5115220	-0.1061245	0.9797090
2020-01-07	-1.3888607	2.2866454	1.3048697	-0.0627141	2.0184237	-0.0946590	1.5115220	0.7964610
2020-01-08	-0.2787888	-1.3888607	2.2866454	1.3048697	-0.0627141	2.0184237	-0.0946590	0.5407023
2020-01-09	-0.1333213	-0.2787888	-1.3888607	2.2866454	1.3048697	-0.0627141	2.0184237	0.5351791
2020-01-10	0.6359504	-0.1333213	-0.2787888	-1.3888607	2.2866454	1.3048697	-0.0627141	0.3376829

本題: ローリング集計をするためのパッケージ

このセクションでわかること

zoo と RcppRoll はtidyverse より前からあるパッケージ
tidyverse 前提の tibbletime::rollify() がある
2019年に runner が, 今年に入ってから slider が登場した
- 正確には slider は tsibble の発展版
slider::slide_vec() を使えばpurrr:map()に似た構文でローリング集計を書け, 複雑な条件に柔軟に対応できる
- ウィンドウ幅を「何件前(後)まで」だけでなく「何日前まで」「何ヶ月前まで」でも指定できる
- 1日1件のデータに整形しなくとも集計できる

新しく登場したパッケージたち

前節に書いたように, dplyrにも時系列処理を意識した関数はいくらか用意されているものの, dplyr 単体ではローリング集計はやや面倒である. そこで外部パッケージに頼ることになる. RcppRoll や zoo パッケージは dplyr 登場以前から存在し, どちらもローリング集計の機能を提供する. mutate() 内でこれらの関数を呼び出せばローリング集計できるが, これらは dplyr 以前に作られたため, ちょっと挙動にくせがある (慣れてしまえば大して気にならない? 全くそのとおり).

そこで, ポストdplyr時代のローリング集計パッケージがないか探してみた. dplyr と親和性のあるパッケージとして, 今回は slider をメインに以下7つを紹介する(最後の2つはおまけだが). 結論から言うと slider に要注目なのだが, そもそも他のパッケージもあまり知られてなさそうなので機能比較をしつつ紹介する. 以降は tidyverse と lubridate パッケージと組み合わせる前提で話を進める.

slider (tsibble)
tibbletime
RcppRoll
zoo
runner
purrr
rsample

各パッケージのバージョン
package	version
tibbletime	0.1.5
slider	0.1.4
tsibble	0.9.1
runner	0.3.7
RcppRoll	0.3.0
zoo	1.8.8
rsample	0.0.7

なお私は参照透過性の弱い書き方を強いる言語・構文が嫌いなので(~~例えば一部のSQLとか~~), そのへんも評価の対象に入っている

`slider` (+ `tsibble`)

tsibble公式マニュアル

sliderパッケージが登場したのはつい最近である(CRANの登録は2020/2/6). tsibbleのほうは以前『CausalImpact でできること, できないこと』でも少し言及したように, 時系列データをdplyrなどと同じように扱うことを目指したパッケージである.tsibbleの使い方全般に関しては開発者による公式チュートリアルを読めばよいだろう. といっても通常のデータフレーム(あるいは tibble)と tsibble の提供するデータフレームとの相互変換くらいしか気にするところはない. 現時点(ver. 0.9.1)でもtsibbleにもローリング集計用の関数は残っているが, 非推奨扱いで, 代わりに sliderパッケージを使うことを開発者は奨めている (なので将来廃止されるかもしれない).

sliderパッケージの関数は, tsibbleの古い関数とよく似た構文だが後で紹介するRcppRollやtibbletimeより後に作られた. これらの機能で開発者が足りないと感じていたものが補われており, 後ろ向きだけでなく前向き集計も簡単にできたり, 期間の集計単位を調整できたりと, ローリング集計に柔軟に対応できるようになっている. 詳しくは以下の公式ドキュメントと, 開発者の書いた紹介記事を参考に.

“Sliding Window Functions • slider”
“slider 0.1.0 - tidyverse”

sliderのこれらの関数はpurrr::map*()シリーズに似て, slide_<返り値の型> という名前で何種類か用意され, さらに変数を2つ取るなら slide2*(), 多くの変数を取るなら pslide*() がある. 移動平均を取りたい場合は, slide_dbl() でできるが, sliderの特徴として変数の型で関数を変えるのが面倒な人向けにslide_vec()が用意されている. 実際今回は mutate() 内で使うこと前提なので, これは便利である⁷.

d %>% mutate(ma7 = slide_vec(.x = origin, .f = mean, .before = 6))  %>%
  filter(date >= ymd("2020-01-01"))

date	origin	ma7
2020-01-01	1.5115220	0.5159882
2020-01-02	-0.0946590	0.3066142
2020-01-03	2.0184237	0.6756316
2020-01-04	-0.0627141	0.6147970
2020-01-05	1.3048697	0.7107980
2020-01-06	2.2866454	0.9797090
2020-01-07	-1.3888607	0.7964610
2020-01-08	-0.2787888	0.5407023
2020-01-09	-0.1333213	0.5351791
2020-01-10	0.6359504	0.3376829

なお, 前向き集計したい場合は, マイナスで指定するか, .before= の代わりに .after= で指定する. 両方指定すれば前後から取れる.

そしてsliderにはウィンドウ幅の設定をもっと簡単にすべくslide_index()やslide_period()が用意されている. slide()関数のウィンドウはあくまで「何件前のレコードか」で設定するが, これら2つの関数は時間インデックスを基準にウィンドウ幅を決める. データに常にレコードが1日1件あるとは限らず, 抜けている日がある, 1日に複数件のレコードがある場合, slide() ならそのような事情に関係なく6件前まで遡るが, slide_index() なら6日前まで遡るし, 時刻でソートされていないとエラーを返してくれる⁸.

d  %>% filter(date != ymd("2019-12-31"), date != ymd("2020-01-05")) %>%
  slide(.x = .$date, .f = ~.x, .before = 6) %>% .[6:14]

## [[1]]
## [1] "2019-12-26" "2019-12-27" "2019-12-28" "2019-12-29" "2019-12-30"
## [6] "2020-01-01"
## 
## [[2]]
## [1] "2019-12-26" "2019-12-27" "2019-12-28" "2019-12-29" "2019-12-30"
## [6] "2020-01-01" "2020-01-02"
## 
## [[3]]
## [1] "2019-12-27" "2019-12-28" "2019-12-29" "2019-12-30" "2020-01-01"
## [6] "2020-01-02" "2020-01-03"
## 
## [[4]]
## [1] "2019-12-28" "2019-12-29" "2019-12-30" "2020-01-01" "2020-01-02"
## [6] "2020-01-03" "2020-01-04"
## 
## [[5]]
## [1] "2019-12-29" "2019-12-30" "2020-01-01" "2020-01-02" "2020-01-03"
## [6] "2020-01-04" "2020-01-06"
## 
## [[6]]
## [1] "2019-12-30" "2020-01-01" "2020-01-02" "2020-01-03" "2020-01-04"
## [6] "2020-01-06" "2020-01-07"
## 
## [[7]]
## [1] "2020-01-01" "2020-01-02" "2020-01-03" "2020-01-04" "2020-01-06"
## [6] "2020-01-07" "2020-01-08"
## 
## [[8]]
## [1] "2020-01-02" "2020-01-03" "2020-01-04" "2020-01-06" "2020-01-07"
## [6] "2020-01-08" "2020-01-09"
## 
## [[9]]
## [1] "2020-01-03" "2020-01-04" "2020-01-06" "2020-01-07" "2020-01-08"
## [6] "2020-01-09" "2020-01-10"

d %>%  filter(date != ymd("2019-12-31"), date != ymd("2020-01-05")) %>%
  slide_index(.x = .$date, .i = .$date, .f = ~.x, .before = 6) %>% .[6:14]

## [[1]]
## [1] "2019-12-26" "2019-12-27" "2019-12-28" "2019-12-29" "2019-12-30"
## [6] "2020-01-01"
## 
## [[2]]
## [1] "2019-12-27" "2019-12-28" "2019-12-29" "2019-12-30" "2020-01-01"
## [6] "2020-01-02"
## 
## [[3]]
## [1] "2019-12-28" "2019-12-29" "2019-12-30" "2020-01-01" "2020-01-02"
## [6] "2020-01-03"
## 
## [[4]]
## [1] "2019-12-29" "2019-12-30" "2020-01-01" "2020-01-02" "2020-01-03"
## [6] "2020-01-04"
## 
## [[5]]
## [1] "2020-01-01" "2020-01-02" "2020-01-03" "2020-01-04" "2020-01-06"
## 
## [[6]]
## [1] "2020-01-01" "2020-01-02" "2020-01-03" "2020-01-04" "2020-01-06"
## [6] "2020-01-07"
## 
## [[7]]
## [1] "2020-01-02" "2020-01-03" "2020-01-04" "2020-01-06" "2020-01-07"
## [6] "2020-01-08"
## 
## [[8]]
## [1] "2020-01-03" "2020-01-04" "2020-01-06" "2020-01-07" "2020-01-08"
## [6] "2020-01-09"
## 
## [[9]]
## [1] "2020-01-04" "2020-01-06" "2020-01-07" "2020-01-08" "2020-01-09"
## [6] "2020-01-10"

さらなる応用として, slide_index() は時間インデックスを指定する必要がある. よって, 日次データを月単位のインデックスに変換して与えることで, 元のデータを変更せずに月単位でローリング集計することも可能になる.

d %>%  filter(date != ymd("2019-12-31"), date != ymd("2020-01-05")) %>%
  mutate(ma_month = slide_index_dbl(.x = origin, .i = as.yearmon(date), .f = mean, .before = 1))

slide_period() もよく似た機能を持つが, こちらはさらに, dplyr::summarise() と組み合わせることでより広く応用できる集計関数としても使える. こちらは例が複雑になるので後の実践パートで具体的な使い方を紹介する.

書かなくても類推できると思うが現時点で日本語情報が全く無いので slide*() 系関数の一覧リストを書いておく. 公式リファレンスはここ.

slide(): 基本. 返り値はリスト
slide_vec(): slide()の結果をベクトル型にして返す. 以降の slide_dbl(), slide_lgl() などを書きわけるのが面倒な人向け.
slide_dbl(): slide()の結果をnumeric型で返す.
slide_int(): slide()の結果をinteger型で返す.
slide_lgl(): slide()の結果をlogical型で返す.
slide_chr(): slide()の結果をcharacter型で返す.
slide_dfr(): 1行のdata.frameとして返す.
slide_dfc(): 1列のdata.frameとして返す.
slide2(): slide() の集計対象の変数を2つ指定できるバージョン, slide() 同様に slide2_vec(), slide2_dbl(), … がある.
pslide(): slide() の集計対象の変数を複数指定できるバージョン. slide() 同様に slide2_vec(), slide2_dbl(), … がある.
slide_index(): slide() とほぼ同じだが, 行数ではなく, 指定された時刻インデックスによってウィンドウ幅のを評価する. slide() 同様に slide_index_vec(), slide_index_dbl(), … がある.
slide_index2(): slider_index() の集計対象の変数を2つ指定できるバージョン. slide() 同様に slide_index2_vec(). slide_index2_dbl(), … がある.
pslide_index(): slide_index() の集計対象の変数を複数指定できるバージョン. slide() 同様に pslide_index_vec(), pslide_index_dbl(), … がある.
slide_period(): slide_index() に週, 月, 四半期, 年単位の集計機能を追加したもの. 例えば period = "month, .every = 1" で1ヶ月単位の集計とする slide() 同様に slide_period_vec(). slide_period_dbl(), … がある.
slide_period2(): slide_period() の集計対象の変数を2つ指定できるバージョン. slide() 同様に slide_period2_vec(). slide_period2_dbl(), … がある.
pslide_period(): slide_period() の集計対象の変数を複数指定できるバージョン. slide() 同様に pslide_period_vec(), pslide_period_dbl(), … がある.
block(): slide() 系関数で指定されたウィンドウ幅に分割した時間インデックスだけを出力する関数. 返り値はリスト.
hop(): slide() の低水準関数バージョン. ウィンドウ幅の境界を細かく指定したいとき, 返す値の長さが入力と違うとき, などより複雑な処理をしたい時に使う. hop2(), hop_index(), hop2_index() も同様.

補足: `tsibble` のウィンドウ関数

なお, tsibble 側に残してある機能として, partial_slider()/partial_pslider() がある. これは指定のウィンドウ幅で切り取って, 対応する値をリストで返すと言う処理(slider::block() はウィンドウ幅に対応する時刻インデックスを返す). これを使って移動平均を出したいなら,

as_tsibble(d, index = date) %>% mutate(
  window7 = partial_slider(.x = origin, .size = 7),
  ma = map_dbl(window7, mean)
  )  %>%
  filter(date >= ymd("2020-01-01")) %>% head

## Warning: `slide()` is deprecated as of tsibble 0.9.0.
## Please use `slider::slide()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

date	origin	window7	ma
2020-01-01	1.5115220	1.3709584, -0.5646982, 0.3631284, 0.6328626, 0.4042683, -0.1061245, 1.5115220	0.5159882
2020-01-02	-0.0946590	-0.56469817, 0.36312841, 0.63286260, 0.40426832, -0.10612452, 1.51152200, -0.09465904	0.3066142
2020-01-03	2.0184237	0.36312841, 0.63286260, 0.40426832, -0.10612452, 1.51152200, -0.09465904, 2.01842371	0.6756316
2020-01-04	-0.0627141	0.63286260, 0.40426832, -0.10612452, 1.51152200, -0.09465904, 2.01842371, -0.06271410	0.6147970
2020-01-05	1.3048697	0.40426832, -0.10612452, 1.51152200, -0.09465904, 2.01842371, -0.06271410, 1.30486965	0.7107980
2020-01-06	2.2866454	-0.10612452, 1.51152200, -0.09465904, 2.01842371, -0.06271410, 1.30486965, 2.28664539	0.9797090

のようになる. partial_slider() の返り値が入った window7 がリストなので, purrr::map_dbl() で平均をとれば移動平均になる. つまり, 集計するには purrrか lapply()が必須である. ローリング集計をしたいだけならば slider パッケージの関数だけで事足りるので, 現状ではあまり使う機会がなさそうだ

tsibbleのウィンドウ関数は以下の4つである. slider::block()と違ってウィンドウ幅に含まれる時刻インデックスではなく, 値を返す.

slider(): ベクトルから指定のウィンドウ幅の要素のリストを返す. 返り値の長さは N - .size なので mutate() で使えない
psluder(): slider() の複数変数バージョン
partial_slider(): ベクトルから指定のウィンドウ幅の要素のリストを返す. .bind = T で元の長さに合わせるためにNAで埋めることができる
partial_pslider(): partial_slider() の複数変数バージョン

`tibbletime`

tibbletimeはtsibbleと同様にdplyrの時系列データ対応版として作られたパッケージである. 今はローリング集計の話に限定するが, 2つのパッケージのどっちが良いかという問題に関しては, Tokyo.R 69回の発表資料『tsibbleとtibbletimeの紹介』の9ページ目の比較表が端的で参考になるだろう.

tibbletimeではrollify() という関数が用意されている⁹. これはローリング集計用関数を生成する, いわゆる生成関数の機能を持つ.

用法を一般化して説明すると, rolify()は適用する関数とウィンドウ幅を指定できるので, slider_vec() と雰囲気は同じである. ただし generating functionなので, 一旦関数を定義してから書く必要がある.

get_ma7 <- rollify(mean, window = 7, na_value = 0)
d %>% mutate(ma7 = get_ma7(origin)) %>%
  filter(date >= ymd("2020-01-01")) %>% head

date	origin	ma7
2020-01-01	1.5115220	0.5159882
2020-01-02	-0.0946590	0.3066142
2020-01-03	2.0184237	0.6756316
2020-01-04	-0.0627141	0.6147970
2020-01-05	1.3048697	0.7107980
2020-01-06	2.2866454	0.9797090

ワンライナーに収めるならこうなる.

d %>% mutate(ma7 = (rollify(mean, window = 7))(origin)) %>%
  filter(date >= ymd("2020-01-01"))

slider と比べると

slide_index() のように時刻インデックス基準でウィンドウ幅を適用する機能がない.
集計関数をその都度定義する必要がある. 特にウィンドウ幅を変えるだけでも別の関数になる.
期間が揃わない間の期間は全てNA扱いになる
ウィンドウ幅を後ろ向きにしか指定できない
必ずベクトルで返さなければならない (リストをネストできない, データフレームを返せない, などの不便さがある)

という仕様が少し使いにくい. (2) は括弧が増えるだけと言えばそうで, (3) 以降も不便になるの状況は限定的ではあるが.

`RcppRoll`

RcppRollのgithubリポジトリ

名前の通り Rcpp (C++) を使って書かれているため, 実行速度に優れているというのが開発者のアピールポイント. 2013年に公開されたパッケージだが, これもあまり情報がない, というか開発者も用例集やマニュアルの類を用意してなかったりする. まあわざわざ説明が必要なほど複雑なものではないので関数ヘルプで十分なのだが.

Rcppは, 速度を重視してか任意の集計関数を指定できない. roll_*() という関数群の中から選ぶことになる. 選べるのは和, 積, 最大・最小値, 平均値, 中央値, 分散, 標準偏差といった基本的な要約統計量だけで, それぞれに左寄せ(=前向き), 右寄せ(=後ろ向き)版関数がある. ただし align= でも指定できる. また, weight= を指定できるので, やろうと思えばガウシアンカーネルによる平滑化などもできることだろう.

d %>% mutate(ma7 = roll_meanr(origin, 7, )) %>%
  filter(date >= ymd("2020-01-01")) %>% head

date	origin	ma7
2020-01-01	1.5115220	0.5159882
2020-01-02	-0.0946590	0.3066142
2020-01-03	2.0184237	0.6756316
2020-01-04	-0.0627141	0.6147970
2020-01-05	1.3048697	0.7107980
2020-01-06	2.2866454	0.9797090

しかし一方で,

ウィンドウ幅の足りない期間は計算できない¹⁰
デフォルトで後ろ向きではなく前後で取り出してしまう
時間インデックス基準のウィンドウ機能がない

という弱点がある. またウィンドウ幅も後ろ向きではなく中央揃えがデフォルト, であることに注意する. 後ろ向き集計する時は, 末尾に r のついた関数を使うと良いだろう.

ヘルプに書いてあることそのままだが, RcppRollの関数一覧もいちおう書いておこう.

roll_sum(): 和. roll_suml(), roll_sumr() はそれぞれ前向き, 後ろ向き集計版だが, align= で指定することもできる.
roll_prod(): 積. roll_prodl(), roll_prodr() はそれぞれ前向き, 後ろ向き集計版だが, align= で指定することもできる.
roll_min(): 最小値. roll_minl(), roll_minr() はそれぞれ前向き, 後ろ向き集計版だが, align= で指定することもできる.
roll_max(): 最大値. roll_maxl(), roll_maxr() はそれぞれ前向き, 後ろ向き集計版だが, align= で指定することもできる.
roll_mean(): 平均値. roll_meanl(), roll_meanr() はそれぞれ前向き, 後ろ向き集計版だが, align= で指定することもできる.
roll_median(): 中央値. roll_medianl(), roll_medianr() はそれぞれ前向き, 後ろ向き集計版だが, align= で指定することもできる.
roll_var(): 分散. roll_varl(), roll_varr() はそれぞれ前向き, 後ろ向き集計版だが, align= で指定することもできる.
roll_sd(): 標準偏差 roll_sdl(), roll_sdr() はそれぞれ前向き, 後ろ向き集計だが, align= で指定することもできる.

`runner`

runner公式マニュアル

runnerはおそらく今回紹介するパッケージのうち, 一番知名度がない. このパッケージは簡単な操作でできることを意識して作られているようだ(個人的にはヘルプに図解を入れてるのが好感度高い). ほとんどの操作は runner::runner() 関数に集約されている.

d %>% mutate(ma7 = runner(x = origin, f = mean, k = 7)) %>%
  filter(date >= ymd("2020-01-01")) %>% head

date	origin	ma7
2020-01-01	1.5115220	0.5159882
2020-01-02	-0.0946590	0.3066142
2020-01-03	2.0184237	0.6756316
2020-01-04	-0.0627141	0.6147970
2020-01-05	1.3048697	0.7107980
2020-01-06	2.2866454	0.9797090

runner は slider に比べてオプションがシンプルにまとまっており, なおかつそれなりに柔軟なことができる. 例えば,

ウィンドウ幅だけデータが揃わない最初の期間を NA に置き換えるか, 除去するかの na_pad. zoo と違い, F の場合はウィンドウ幅に足りなくても計算する
ウィンドウをシフトさせる lag (これで後ろ向き・前向きにとどまらないより柔軟なアラインメント指定ができる).
k を指定しない場合ウィンドウ幅が過去全ての値を対象となる
idx に Date や POSIX 型の時刻インデックス別を指定すれば, ラグやウィンドウ幅の指定を "1 month" というふうに書ける
入力データにはベクトルだけでなくデータフレームも与えられるため, 複数列を使った集計 (例えば相関係数) も可能

これらは他のパッケージでは見られない使い勝手を工夫した仕様だと思う. 例えば, 今回の例は以下のように書くこともできる.

d %>% mutate(ma7 = runner(x = origin, f = mean, idx = date, k = "1 week")) %>%
  filter(date >= ymd("2020-01-01")) %>% head

date	origin	ma7
2020-01-01	1.5115220	0.5159882
2020-01-02	-0.0946590	0.3066142
2020-01-03	2.0184237	0.6756316
2020-01-04	-0.0627141	0.6147970
2020-01-05	1.3048697	0.7107980
2020-01-06	2.2866454	0.9797090

一方で, 集計関数には tidyverse 系でよくある ~sum(.x) のような formula 式には対応していない.

`zoo`

既に書いたようにzooはかなり古いパッケージで, dplyr より前から存在する. 機能も豊富で, ローリング集計に関しては rollmean(), rollmax など RcppRoll の関数に対応する(というよりRcppRollの関数名がzooを参考にしている)ものや, 任意の集計関数を適用できる rollaplly(), rollapplyr() があるほか, 時系列データの処理に役に立つ関数がいくつも用意されている. しかし, data.frame や tibble ではなく matrix を拡張した zoo クラスでデータを保持する. データフレームで扱えないわけではないが, ものによってはクラス変換に気をつける必要がある.

`purrr`

purrr 公式マニュアル

紹介はするが一応これでもできますよという程度で, 私としてはこれで書くのはおすすめしない. purrr と言えば並列処理を簡単に書くためのパッケージである. しかし並列処理というのは1つ1つを独立して処理するため, 1つ前の値に応じて処理を変える再帰的な処理はできない. しかし時系列データの処理の多くは移動平均やラグのように, 直前の値を使って計算するものが多い. そのような処理をするために, accumulate(), reduce() という関数が用意されている. 日本語で言及しているものに Heavy Watal 氏の紹介ページがある. しかし, 基本的に1件前のレコードを参照できるだけなので, 2時点以上を参照するローリング集計をしようとすると一気に複雑になる, というか冒頭のラグを手動で書くのと大差ない.

`rsample`

rsample 公式マニュアル

rsampleは機械学習や統計モデリングのために使う再標本処理のためのパッケージで, その中にウィンドウ関数に基づいて入力データを分割する機能がある. よって移動平均などの計算もできることはできるのだが, そういった集計をするには構文がものものしくなりすぎている.

まず, rolling_origin() は指定されたウィンドウ幅に従ってデータを分割し, さらにそれぞれを analysis と assessment と呼ばれる2つに分割する¹¹. 分割したサブデータセットは analysis()/assessment() でアクセスできる. これをパイプでつなげて移動平均の集計に使う場合, 以下のようになる.

rolling_origin(d,  initial = 7, assess = 0, cumulative = F) %>%
  mutate(
    date = map(splits, ~max(analysis(.x)$date)),
    ma7 = map_dbl(splits, ~mean(analysis(.x)$origin))
    ) %>% unnest(date) %>% select(-splits, -id) %>% head

date	ma7
2020-01-01	0.5159882
2020-01-02	0.3066142
2020-01-03	0.6756316
2020-01-04	0.6147970
2020-01-05	0.7107980
2020-01-06	0.9797090

このように, ウィンドウ関数による分割までが rsample の提供する機能なので, さらに集計関数を挟む必要があり, 結局 tsibble パッケージの slider() と大差ない. よって今回は参考として挙げるだけに留めておく.

実例で各パッケージを比較

既に言及した例や公式マニュアルのは簡単な使用法だけなので, ここでは実用性を意識してもう少し凝った例で slider, tibbletime, RcppRoll, zoo, runner の5つのパッケージの使い勝手を検証する. なお, lubridate パッケージの日付処理関数群や最近出た dplyr ver. 1.0.0 の構文も多用しているので注意.

Q1 (簡単): 毎日の売上集計

df1 は, 2020/1/1-2020/12/31 の期間の4種類の商品の販売記録である(という設定). sale_1-sale_4はそれぞれ商品1-4の売上額である. 商品別に, 1日あたりの売上額について1週間(7日間)の移動平均を取れ. また, それをグラフに表せ.

set.seed(42)
df1 <- expand_grid(date = seq(ymd("2020-01-01"), ymd("2020-12-31"), by = "1 day"), item = 1:4) %>%
  mutate(sale = round(abs(rnorm(n = n())) * 100, 0)) %>%
  pivot_wider(id_cols = date, names_from = item, values_from = sale, names_prefix = "sale_")
head(df1)

date	sale_1	sale_2	sale_3	sale_4
2020-01-01	137	56	36	63
2020-01-02	40	11	151	9
2020-01-03	202	6	130	229
2020-01-04	139	28	13	64
2020-01-05	28	266	244	132
2020-01-06	31	178	17	121

これはかなり簡単で, 5つのパッケージどれでもほぼおなじ書き方になる. dplyr 1.0.0 の across() を早速使って4つの列に同じローリング集計を適用してもいいし, 一旦 sale_* をロングにして, group_by() でグループ別にローリング集計することもできる. グラフに表示することを考えるとロングの方が ggplot2 で扱いやすいため, 後者でやってみる.

q1 <- df1 %>% pivot_longer(cols = sale_1:sale_4, names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  group_by(item) %>% mutate(ma6_sale = RcppRoll::roll_meanr(sale, n = 7))

上記は RcppRoll の場合だが, 他のパッケージの場合はこの部分をそれぞれ以下のように置き換える. zoo::rollmeanr() はウィンドウ幅が足りない最初の6日分を切り捨ててしまうので, na.pad = T の指定が必要である. slider::slide_vec() はデフォルトでウィンドウ幅に満たなくとも計算するので, 他の関数と結果を一致させるためにここでは complete = T を指定している. 一方で, ウィンドウ幅に満たなくとも計算するという処理は他の3つの関数ではできない. また, slider の関数はどれもウィンドウ幅の指定に現時点を含まないので, .before = 6 を指定する.

(tibbletime::rollify(.f = mean, window = 7))(sale)
zoo::rollmeanr(sale, k = 7, na.pad = T)
slider::slide_vec(sale, .f = mean, .before = 6, .complete = T)
runner::runner(sale, f = mean, k = 7, na_pad = T)

ただし, 現状 runner のみ, ベクトルを与えた場合, na_pad = T しても NA ではなく NaN になることが多いので注意 (長さゼロのベクトルを返してしまうため. 多分開発者のミス).

グラフにするとこうなる

各商品の売上額1週間移動平均

Q2 (まあ簡単): 複数店舗の売上集計

df2 は, 2020/1/1-2020/12/31 の期間の4種類の商品の, 3つの店舗での販売記録である sale_1-sale_4はそれぞれ商品1-4の売上額である. 商品別に, 1日あたりの全店舗の合計売上額について1週間(7日間)の移動平均を取れ. また, それをグラフに表せ.

set.seed(42)
df2 <- expand_grid(date = seq(ymd("2020-01-01"), ymd("2020-12-31"), by = "1 day"), shop = LETTERS[1:3], item = 1:4) %>%
  mutate(sale = round(abs(rnorm(n = n())) * 100, 0)) %>%
  pivot_wider(id_cols = c(date, shop), names_from = item, values_from = sale, names_prefix = "sale_")
head(df2)

date	shop	sale_1	sale_2	sale_3	sale_4
2020-01-01	A	137	56	36	63
2020-01-01	B	40	11	151	9
2020-01-01	C	202	6	130	229
2020-01-02	A	139	28	13	64
2020-01-02	B	28	266	244	132
2020-01-02	C	31	178	17	121

まずRcppRollを使う例を紹介する. 既に説明したように1日1件のデータでないと集計できないため, 一旦店舗をグループ集計で縮約する必要がある. 今回紹介する関数はいずれも mutate() 内で呼び出せるので, 売上額をlongにしてから group_by して, roll_meanr() を適用する.

q2_rcpproll <- df2 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(item, date) %>% group_by(item, date) %>% summarise(sale = sum(sale), .groups = "drop") %>%
  group_by(item) %>% mutate(ma6_sale = RcppRoll::roll_meanr(sale, n = 7)) %>% ungroup

tibbletime, zoo も引き続き同じようにやるしかない

(tibbletime::rollify(.f = mean, window = 7))(sale)
zoo::rollmeanr(sale, k = 7, na.pad = T)

RcppRollでの集計結果

slider の場合も slide_vec() に置き換えることで同様にできるが, slide_period() を使うことで1日1件のデータでなくとも集計できるので試してみる. ここでは slide_period_dfr() を使ってデータフレームで返す. 開発者がtidyverse.orgに書いた記事では集計関数を先に定義しているが, もちろんインラインで書くこともできる. ただし, 今回は商品 (item) 別に集計する必要があるため, 一旦 tidyr::nest() で折りたたんでから mutate(map(...)) 内で並列処理している. また, summarise() 内に日付を自分で書く必要があることに注意

dplyr 1.0.0 以降で, summarise() が複数の値を返す集計関数を許容するようになった¹²ため, 以前よりも簡単に書ける. 一方で, 集計関数にそのまま mean を与えると間違った結果になる. 日数ではなく7日間に含まれる件数で割ってしまうからだ. また, summarise() 内で読んでいるせいでせっかく NA でパディングしてくれたものが消され, 日付が 1/7 からになってしまうので集計関数をいじって NA を出している (日付リストを完全にすることにこだわらないなら ~sum(.x) / 7 で十分)

q2_slider <- df2 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(item, date) %>% group_by(item) %>%
  summarise(ma6_sale = slider::slide_period_dbl(sale, .i = date, .period = "day", .f = ~if (is.null(.x)) NA else sum(.x)/7, .before = 6, .complete = T),
            date = unique(date), .groups = "drop")

というわけで思ったより簡単にはならない…

そして dplyr のバージョンが古い場合は以下のようにさらに複雑になり, 実行速度も低下してしまう.

q2_slider <- df2 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(item, date) %>% group_by(item)  %>% nest %>% mutate(
    data = map(data, ~slider::slide_period_dfr(.x, .i = .x$date,  .period = "day", .f = ~summarise(.x, date = max(date), ma6_sale = sum(sale)/7), .before = 6, .complete = T))
  ) %>% unnest(data) %>% ungroup

runner も内部で集計ができるので試してみる. 集計関数のところが少し厄介で, そのまま mean を与えると間違った結果になるのは slider と同じだ.

q2_runner <- df2 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(item, date) %>% group_by(item) %>% summarise(
    ma6_sale = runner::runner(sale, f = function(x) if(length(x) == 0) NA else sum(x) / 7, k = "7 days", idx = date, at = unique(date), na_pad = T) %>% if_else(is.nan(.), NA_real_, .),
    date = unique(date),
    .groups = "drop")

Q3 (ちょっと厄介): 売上が計上されないときもある

Q2では毎日A, B, C の3店舗の売上が計上されていた. しかし今, df3は売上がない時はレコードが計上されないことがある. df3をQ2と同様に集計せよ.

set.seed(42)
df3 <- group_by(df2, date) %>% nest %>% ungroup %>% sample_frac(size = .8, replace = F) %>% arrange(date) %>% unnest(cols = data)

date	shop	sale_1	sale_2	sale_3	sale_4
2020-01-07	A	62	95	54	58
2020-01-07	B	77	46	89	110
2020-01-07	C	151	26	9	12
2020-01-09	A	113	146	8	65
2020-01-09	B	120	104	100	185
2020-01-09	C	67	11	42	12
2020-01-10	A	19	12	3	11

まずRcppRollを使う例を紹介する. Q2同様に1日1件のデータに整形してからローリング集計関数を適用する必要があるので, まず1日単位で集計し, それから脱落している日を補い, 欠損値として新たに行を作る必要がある. つまり, データの「脱落」と「欠損」の違いは以下のようになる.

date.x	sale.x	date.y	sale.y
2020-01-01	-0.234365277305921	2020-01-01	-0.234365277305921
		NA	NA
		NA	NA
2020-01-04	-0.271763715111397	2020-01-04	-0.271763715111397
2020-01-05	0.947951995875196	2020-01-05	0.947951995875196

366日分の日付データ tibble(date = seq(ymd("2020-01-01"), ymd("2020-12-31"), by = "1 day")) を結合する手もあるが, 今回は tidyr::complete() を使う¹³. そして日付順にソートした上で(忘れるとおかしな結果になる), 前回同様に移動平均を集計できる. また, 今度は NA を含む可能性があるため忘れずに集計関数に na.rm = T を指定する.

q3_rcpproll <- df3 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(item, date) %>% group_by(item, date) %>%
  summarise(sale = sum(sale, na.rm = T), .groups="drop") %>%
  complete(date = full_seq(date, 1), item) %>% arrange(item, date) %>% group_by(item) %>%
  mutate(ma6_sale = RcppRoll::roll_meanr(sale, n = 7, na.rm = T))

この例では 366日分あるが, 営業日のみにしたいというのなら complete() の直後に祝休日を除外する. 例: filter(., !wday(date) %in% c(0, 6), !date %in% unlist(zipangu::jholiday(2020))) のようにする. zipangu::jholiday() はある年の日本の祝日を計算する関数である.

Q2, RcppRoll使用

tibbletime の場合も関数以外全く同じになる. NA が存在するので, ~sum(.x, na.rm = T) という無名関数を与える

q3_tibbletime <- df3 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(item, date) %>% group_by(item,date) %>%
  summarise(sale = sum(sale, na.rm = T), .groups="drop") %>%
  complete(date = full_seq(date, 1), item) %>% arrange(item, date) %>% group_by(item) %>%
  mutate(ma6_sale = tibbletime::rollify(~mean(.x, na.rm = T), window = 7)(sale))

zoo も欠けた日付を補う関数としてna.locf() があるが, これはついでに値も勝手に補間してしまう. さらに言うなら, 現時点ではデータフレームを zooやxtsに適切に変えられる関数が存在しない¹⁴. 結局 complete() や完全な日付の列を結合したほうがはるかに簡単である.

q3_zoo <- df3 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(item, date) %>% group_by(item, date) %>%
  summarise(sale = sum(sale, na.rm = T), .groups="drop") %>%
  complete(date = full_seq(date, 1), item) %>% arrange(item, date) %>% group_by(item) %>%
  mutate(ma6_sale = zoo::rollmeanr(sale, k = 7, na.pad = T, na.rm = T))

slider::slide_period() は日付基準でウィンドウ幅を適用できるが, 日付そのものに欠けがあるとその日の移動平均が集計できないので, sliderを使う場合も日付を補う必要がある. tidyr::complete() の別解としてtsibble::fill_gaps() を使う例を紹介する.

runner も同様に, 以下のように書く.

q3_runner <- df3 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(item, date) %>% group_by(item, date) %>%
  summarise(sale = sum(sale, na.rm = T), .groups="drop") %>%
  as_tsibble(index = date, key = item) %>% fill_gaps %>% as_tibble %>%
  arrange(item, date) %>% group_by(item) %>%
  mutate(ma6_sale = runner::runner(sale, f = mean, k = 7, na_pad = T, na.rm = T) %>% if_else(is.nan(.), NA_real_, .))

Q4 (だいぶ厄介): 各商品の過去1週間の売上シェアを毎日計算する

df4 は, 2020/1/1-2020/12/31 の期間の4種類の商品の, 3つの店舗での販売記録である sale_1-sale_4はそれぞれ商品1-4の売上額である. 商品別に, 毎日, 過去1週間の売上シェアを計算せよ

これはR-wakalangに投稿された質問を少し改変したもの. 売上シェアの移動平均, つまり1日毎にシェアを計算してその移動平均を取るわけではないことに注意. RcppRoll などは複数の変数を同時に扱えないため, 分子と分母をそれぞれ移動合計を取ってから計算する必要がある.

df4 <- expand_grid(date = seq(ymd("2020-01-01"), ymd("2020-12-31"), by = "1 day"), item = 1:4, shop = LETTERS[1:4]) %>%
  mutate(sale = round(abs(rnorm(n = n())) * 100, 0)) %>%
  pivot_wider(id_cols = c(date, shop), names_from = item, values_from = sale, names_prefix = "sale_")

date	shop	sale_1	sale_2	sale_3	sale_4
2020-01-01	A	120	135	73	138
2020-01-01	B	47	2	100	205
2020-01-01	C	27	24	126	102
2020-01-01	D	39	94	125	3
2020-01-02	A	70	120	74	87
2020-01-02	B	97	19	5	97

よって, rcpproll ではこれまで同様にまず集計しやすいlongな形式に変換してから計算する.

q4_rcpproll <- df4 %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
  arrange(date, item) %>% group_by(date, item) %>% summarise(sale = sum(sale), .groups = "drop_last") %>%
  mutate(total = sum(sale)) %>% ungroup %>% arrange(item, date) %>% group_by(item) %>%
  mutate(across(c(sale, total), ~RcppRoll::roll_sumr(.x, n = 7)), share = sale / total)

zoo, tibbletime もほぼ同じ書き方しかないので省略するとして, slider, runner は多変数の集計関数を扱えるため¹⁵再び工夫の機会がある. 多変数を扱えることを利用して, ワイドな状態から直接シェアの移動集計をするように書いてみる.

q4_slider <- df4 %>% rowwise %>% arrange(date) %>%
  mutate(total = sum(c_across(starts_with("sale")))) %>% ungroup %>%
  summarise(across(starts_with("sale"),
                   ~slide_period2_vec(.x, total, .i = date, .period = "day", .f = ~sum(.x) / sum(.y), .before = 6, .complete = T)),
            date = unique(date)) %>%
  pivot_longer(starts_with("sale"), names_pattern = "sale_(.+)", names_to = "item", values_to = "share")

同様に, runner でも以下のように書く.

q4_runner <- df4 %>% rowwise %>% arrange(date) %>% mutate(total = sum(c_across(starts_with("sale")))) %>% ungroup %>%
  summarise(
    across(starts_with("sale"), ~runner::runner(tibble(n = .x, d = total), idx = date, at = unique(date), f = function(x) sum(x$n) / sum(x$d), k = 7, na_pad = T)),
    date = unique(date)) %>%
  pivot_longer(cols = starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "share")

runner() に summarise() を与えるやり方もあるがこちらのほうが遅い & 複雑になる. 例えば以下のように(ただしこれは体裁を完全に整えてはいない).

q4_runner <- df4 %>% rowwise %>% arrange(date) %>% mutate(total = sum(c_across(starts_with("sale")))) %>% ungroup %>%
  runner::runner(idx = .$date, at = unique(.$date), f = function(x) summarise(x, across(starts_with("sale"), ~sum(.x) / sum(total)), date = max(date)), k = 7, na_pad = T)

ここからさらに「1ヶ月ごと」「四半期ごと」など条件を変えていけば slider や runner のほうが便利になる場合も増えていくのだが, 例題を考えるのがめんどくさくなってきたのでこの辺にしておく.

実行時間の比較

最後に, それぞれの書き方で実行速度がどれくらい変わるかを確認する. 1つはローリング集計関数だけの呼び出し速度, 処理が複雑なQ4の実行速度を調べた.

set.seed(42)
res <- list()
for(n in 2:5){
  print(10^n)
  res[[as.character(10^n)]] <- microbenchmark(
    roll_meanr(x, n = 5),
    (rollify(.f = mean, window = 5))(x),
    tibbletime_ma(x),
    rollmeanr(x, k = 5, na.pad = T),
    rollapplyr(x, width = 5, FUN = mean),
    slide_vec(x, .f = mean, .before = 4, .complete = T),
    slide_dbl(x, .f = mean, .before = 4, .complete = T),
    runner(x, f = mean, k = 5, na_pad = T),
    setup = {
      tibbletime_ma <- rollify(mean, 5)
      x <- rnorm(n = 10^n)
      }, times = 100
    ) %>% as_tibble %>%
    mutate(name = c("RcppRoll", "tibbletime-inline", "tibbletime-predefined", "zoo-rollmeanr", "zoo-rollapplyr", "slider-vec", "slider-dbl", "runner")[as.integer(expr)]) %>%
    mutate(size = 10^n)
}
res <- bind_rows(res) %>% mutate(time = time / 1e8)

各パッケージの所要時間

1万件程度まではそれぞれの実行速度は体感でほとんど変わらないが, 10万件になるとかなり差が付いてくる. その中ではやはり RcppRoll が最も早かった. が, zoo も割と追いついている. tibbletime::rollify() は群を抜いて遅く¹⁶, 特別な機能もないため現状では使うメリットはなさそうだ. slider と runner は前者のほうが少し速く, slide_vec より slide_dbl で型指定をした場合の差はわずかに早くなる程度だった.

当初は slider や runner が便利だと喧伝したかったが, 思ったよりメリットが少ない. RcppRoll より高速だとは思っていなかったが, 古いzooにもそこそこ差をつけられているのはいただけない.

そこで, 極端に遅い rollify() を除外した上でさらにQ4の回答例のプログラムの実行速度を比較する.

set.seed(42)
res2 <- list()
for(n in 2:4){
  print(10^n)
  s <- Sys.time()
  df <- expand_grid(date = seq(ymd("1900-01-01"), ymd("1900-01-01") + days(10^n), by = "1 day"), item = 1:4, shop = LETTERS[1:4]) %>%
        mutate(sale = round(abs(rnorm(n = n())) * 100, 0)) %>%
        pivot_wider(id_cols = c(date, shop), names_from = item, values_from = sale, names_prefix = "sale_")
  res2 <- c(res2, microbenchmark(
    df %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
      arrange(date, item) %>% group_by(date, item) %>% summarise(sale = sum(sale), .groups = "drop_last") %>%
      mutate(total = sum(sale)) %>% ungroup %>% arrange(item, date) %>% group_by(item) %>%
      mutate(across(c(sale, total), ~roll_sumr(.x, n = 7)), share = sale / total),
    df %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
      arrange(date, item) %>% group_by(date, item) %>% summarise(sale = sum(sale), .groups = "drop_last") %>%
      mutate(total = sum(sale)) %>% ungroup %>% arrange(item, date) %>% group_by(item) %>%
      mutate(across(c(sale, total), ~rollsumr(.x, k = 7, na.pad = T)), share = sale / total),
    df %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
      arrange(date, item) %>% group_by(date, item) %>% summarise(sale = sum(sale), .groups = "drop_last") %>%
      mutate(total = sum(sale)) %>% ungroup %>% arrange(item, date) %>% group_by(item) %>%
      mutate(across(c(sale, total), ~slide_dbl(.x, .f = mean, .before = 6, .complete = T)), share = sale / total),
    df %>% rowwise %>% arrange(date) %>%
      mutate(total = sum(c_across(starts_with("sale")))) %>% ungroup %>%
      summarise(across(starts_with("sale"),
                       ~slide_period2_dbl(.x, total, .i = date, .period = "day", .f = ~sum(.x) / sum(.y), .before = 6, .complete = T)),
                date = unique(date)) %>%
      pivot_longer(starts_with("sale"), names_pattern = "sale_(.+)", names_to = "item", values_to = "share"),
    df %>% rowwise %>% arrange(date) %>% mutate(total = sum(c_across(starts_with("sale")))) %>% ungroup %>%
    summarise(
      across(starts_with("sale"), ~runner(tibble(n = .x, d = total), idx = date, at = unique(date), f = function(x) sum(x$n) / sum(x$d), k = 7, na_pad = T)),
      date = unique(date)) %>%
    pivot_longer(cols = starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "share"),
    times = 100
    ) %>% as_tibble %>%
    mutate(name = c("RcppRoll", "zoo", "slider-dbl-long", "slider-dbl-wide", "runner")[as.integer(expr)]) %>%
    mutate(size = 10^n) %>% list)
  print(Sys.time() - s)
}
n <- 5
df <- expand_grid(date = seq(ymd("1900-01-01"), ymd("1900-01-01") + days(10^n), by = "1 day"), item = 1:4, shop = LETTERS[1:4]) %>%
        mutate(sale = round(abs(rnorm(n = n())) * 100, 0)) %>%
        pivot_wider(id_cols = c(date, shop), names_from = item, values_from = sale, names_prefix = "sale_")
s <- Sys.time()
print(10^n)
res2 <- c(res2, microbenchmark(
    df %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
      arrange(date, item) %>% group_by(date, item) %>% summarise(sale = sum(sale), .groups = "drop_last") %>%
      mutate(total = sum(sale)) %>% ungroup %>% arrange(item, date) %>% group_by(item) %>%
      mutate(across(c(sale, total), ~roll_sumr(.x, n = 7)), share = sale / total),
    df %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
      arrange(date, item) %>% group_by(date, item) %>% summarise(sale = sum(sale), .groups = "drop_last") %>%
      mutate(total = sum(sale)) %>% ungroup %>% arrange(item, date) %>% group_by(item) %>%
      mutate(across(c(sale, total), ~rollsumr(.x, k = 7, na.pad = T)), share = sale / total),
    df %>% pivot_longer(starts_with("sale"), names_to = "item", names_pattern = "sale_(.+)", values_to = "sale") %>%
      arrange(date, item) %>% group_by(date, item) %>% summarise(sale = sum(sale), .groups = "drop_last") %>%
      mutate(total = sum(sale)) %>% ungroup %>% arrange(item, date) %>% group_by(item) %>%
      mutate(across(c(sale, total), ~slide_dbl(.x, .f = mean, .before = 6, .complete = T)), share = sale / total),
    df %>% rowwise %>% arrange(date) %>%
      mutate(total = sum(c_across(starts_with("sale")))) %>% ungroup %>%
      summarise(across(starts_with("sale"),
                       ~slide_period2_dbl(.x, total, .i = date, .period = "day", .f = ~sum(.x) / sum(.y), .before = 6, .complete = T)),
                date = unique(date)) %>%
      pivot_longer(starts_with("sale"), names_pattern = "sale_(.+)", names_to = "item", values_to = "share"),
    times = 100
) %>%
  as_tibble %>%
  mutate(name = c("RcppRoll", "zoo", "slider-dbl-long", "slider-dbl-wide")[as.integer(expr)]) %>%
  mutate(size = 10^n) %>% list)
print(Sys.time() - s)

今度は runner の遅さが目立ったため, 10万件の場合のみ runner の結果を省略した.

Q4の実行時間比較

どうやら slider でまとめてやるよりも RcppRoll で個別に集計したほうが速いらしい. slider開発者はRcppRoll に対してsliderのオーバーヘッドは集計関数が組み込まれていて変更できないか, 自由に決めて呼び出せるかの違いで発生すると述べている¹⁷が, それ以外にも高速化の余地はありそうだ.

slider はもっと複雑な移動集計にも活用できるはずだが, 実際には tidyverse でかなり対処できてしまうため, tidyverse でデータを tidy な構造にしてから RcppRoll で集計するというアプローチのほうがやりやすいことが多い. 単に移動平均を取りたいだけならば RcppRoll::roll_meanr() が最速なので, 殆どの場合で RcppRoll, やや複雑な条件でやりたい時(そしてデータが極端に多くない時)には slider や runner の出番, ということになる.

tibbletime::rollify() が現状では使える場面がないとは既に言った通りだが, そうなると zoo パッケージも rollapply()/rollapplyr()以外はRcppRollの関数が上位互換になっている¹⁸. よって, ほとんどの場面で RcppRoll でやるのがよく, 基本的な要約統計量以外の集計, とくに複数の変数で計算する必要がある場合は sliderかrunnerを使うというのが良さそうだ¹⁹.

この記事を読んでわかったこと (まとめ)

RcppRoll パッケージはなんだかんだで早いができることが少ない
slider/runner は柔軟なローリング集計ができるがまだ速さでRcppRollにかなわない
dplyr 1.0.0 の新機能に across() が取り沙汰されるが summarise() の仕様変更も侮れない
tsibble::fill_gaps() で時系列データの行を補間
zoo は使わなくてもなんとかなる

matrix 形式のデータは ggplot2 とも相性が悪い zoo は組み込みデバイスに対応したグラフ描画関数を用意しているが, これも ggplot2 に慣れてしまった今となっては億劫になる原因だ.↩︎
先日の Tokyo-R でも, メインテーマではないが最新の dplyr 1.0.0 を使ったプログラムを紹介する予定だったのだが, アクシデントにより叶わなかった.↩︎
広義の窓関数は信号処理に使われ, 今回紹介するような集計処理ではあまり使わないものが多い. しかし Python のデータ集計モジュールの pandas では, ローリング集計用のクラスだけでなくなぜか信号処理を想定してるかのように窓関数が充実している↩︎
集計ではなく単にラグを取るだけの処理も広義のローリング集計なので, dplyr の公式マニュアルでは “window function” という見出しで lag() 関数など時系列データ操作に使う関数全般を紹介している↩︎
これは各時点より後ろの期間を取る, 「後ろ向き」集計である「前向き」の1週間を取りたいこともあるが, 単純な集計では後ろ向きを知りたいことのほうが多いはずだし, 本質的な違いはないので以降は後ろ向きの集計に限定して説明する. また, 各時点を中心として前後で取り出したい場合もある. これは「平滑化」と呼ばれる処理である.)↩︎
dplyr 1.0.0 以降ならば, 平均を取る箇所を rowwise %>% mutate(ma = rowMeans(across(c(origin, starts_with("lag"))))) %>% と書くこともできるが, 結局大量にラグ変数を作る必要があるため面倒さはあまり変わらない↩︎
purrr でも map_vec() を用意すると言う話はあるが未実装である https://github.com/tidyverse/purrr/issues/435 ↩︎
slide_index()のこの機能は data.table で利用する際にも役立つと開発者が述べている↩︎
日本語では, 誰が作ったのか不明で, またほとんどがヘルプの用例そのままだが, http://delta0726.web.fc2.com/packages/finance/tibbletime.html で使用例が紹介されている.↩︎
partial という引数が用意されているが, “currently unimplemented” とのこと. なお更新は2年前から止まっている.↩︎
なぜさらに2分割するかというと, 統計モデルの評価で必要になる場合があるのだが, 意外なことにこのことをちゃんと説明している資料は日本語はおろか英語でもほとんどない. なのでそのうち詳しい説明をどこかで書く.↩︎
dplyr 1.0.0 において, 個人的にはこれが一番影響の大きな変更だと思うのだが, 既にいくつか出ている紹介記事の中ではここくらいでしか触れられていない.↩︎
詳しい使い方は https://notchained.hatenablog.com/entry/2017/04/04/220236 などを参考に↩︎
もし知っている人がいたら教えてください↩︎
zoo::rollaply() でも, 各系列に並行して同じ処理を適用するだけなので2系列の相関係数などを計算するのは不可能である. その他 zoo をベースにした特定の集計のみを行える関数がいくつかみられる.↩︎
ソースコードを確認したら, 他のパッケージでは部分的に Rcpp で高速化を図っているものが多いのに対しこのウィンドウ関数は apply() を呼び出して計算しているだけだった.↩︎
https://davisvaughan.github.io/slider/index.html ↩︎
zooパッケージには他にも時系列データの補間などよく使われる機能が用意されているが, 今はそれらの処理も zoo を使わず tidyverse 系のパッケージだけでできるようになっていることが多い. しかし, ファイナンス系のパッケージの多くはいまだ zoo に依存しているようなのでそっちは乗り換えにまだ時間がかかりそうだ. tidyquant というパッケージはあるが.↩︎
いくつか未確認な点が残されている. 例えばzooは系列を並列処理してくれるので, itemの数が多い時により有利になるかもしれない. また, ウィンドウ幅の増減でパフォーマンスにどう差が出るかも確認していない. しかしそういう話は細かすぎるので今回は省略した.↩︎

[R] dplyr 1.0.0 時代の時系列データ処理 — 特に移動集計 —

ill-identified

この記事を読んでわかること (概要)

初めに/能書き

時系列データ処理の基本的なこと

このセクションでわかること

`dplyr` の時系列データ処理

移動集計 (ローリング集計) とは

本題: ローリング集計をするためのパッケージ

このセクションでわかること

新しく登場したパッケージたち

`slider` (+ `tsibble`)

補足: `tsibble` のウィンドウ関数

`tibbletime`

`RcppRoll`

`runner`

`zoo`

`purrr`

`rsample`

実例で各パッケージを比較

Q1 (簡単): 毎日の売上集計

Q2 (まあ簡単): 複数店舗の売上集計

Q3 (ちょっと厄介): 売上が計上されないときもある

Q4 (だいぶ厄介): 各商品の過去1週間の売上シェアを毎日計算する

実行時間の比較

この記事を読んでわかったこと (まとめ)

[R] dplyr 1.0.0 時代の時系列データ処理 — 特に移動集計 —

ill-identified

この記事を読んでわかること (概要)

初めに/能書き

時系列データ処理の基本的なこと

このセクションでわかること

dplyr の時系列データ処理

移動集計 (ローリング集計) とは

本題: ローリング集計をするためのパッケージ

このセクションでわかること

新しく登場したパッケージたち

slider (+ tsibble)

補足: tsibble のウィンドウ関数

tibbletime

RcppRoll

runner

zoo

purrr

rsample

実例で各パッケージを比較

Q1 (簡単): 毎日の売上集計

Q2 (まあ簡単): 複数店舗の売上集計

Q3 (ちょっと厄介): 売上が計上されないときもある

Q4 (だいぶ厄介): 各商品の過去1週間の売上シェアを毎日計算する

実行時間の比較

この記事を読んでわかったこと (まとめ)

`dplyr` の時系列データ処理

`slider` (+ `tsibble`)

補足: `tsibble` のウィンドウ関数

`tibbletime`

`RcppRoll`

`runner`

`zoo`

`purrr`

`rsample`