Salmon, M., Schumacher, D., Hohle, M. (2016) Monitoring Count Time Series in R: Aberration Detection in Public Health Surveillance. Journal of Statistical Software. 70(10).

　カウント時系列監視のためのRパッケージsurveillanceの解説。実戦投入しようかな？ と思ってめくってみた。サーベイランスと言っても広うございますが、これは疫学の文脈での、アウトブレイク検出を意図したパッケージである。
　個人的な好みの問題だと思うけど、品質管理系の論文より10倍くらいわかりやすいような気がする…

1. イントロダクション
　Rパッケージsurveillanceは、カウント・割合の時系列における異常(aberration)検出のいろんな手法を実装している。手法のレビューとしてBuckeridge(2007 J.BiomedicalInformatics), Unkel et al.(2012 JRSS-A)をみよ。[…]

2. surveillanceパッケージ入門
2.1 時系列と関連情報の保存方法
　以下、時系列の長さを\(n\)、エンティティ(地域とか)の数を\(m\)とし、時点\(t\)、エンティティ\(i\)でのカウントを\(y_{it}\)とする。
　時系列はS4クラスstsに格納する。多変量でもよい。[…スロットの説明…　オブジェクトの作り方… メソッド…]

2.2 異常検出アルゴリズムの使い方
　単変量時系列\(\{y_t\}\)について考える。未知の時点\(\tau\)における重要な変化をとらえるというのがお題である。
　任意の時点 \(t_0 \geq 1\)において利用可能な過去のカウントを\(\mathbf{y}_{t_0} = \{y_t; t \leq t_0\}\)とする。そこから統計量\(r(\mathbf{y}_{t_0})\)をつくる。既知の閾値を\(g\)として、警告時点は$$ T_A = \min(t_0 \geq 1: r(\mathbf{y}_{t_0}) \gt g) $$ となる。[なんなのこのコロン… \(r(\mathbf{y}_{t_0}) \gt g\)を満たす最初の\(t_0\)ってことだろうけど]

　まずは、本パッケージが提供する異常検出アルゴリズムをひとつ紹介しよう。
　CDC[米疾病予防管理センター]のEARS (Early Aberration Detection System) という方法が、earsC() 関数で実装してある。手法の詳細についてはFricker, Hegler, Dunfee(2008 Stat.Med.)をみよ。C1, C2, C3という3つのバージョンがある。C1, ベースライン7時点で説明する。
　評価時点\(t_0\)の前の7時点の平均を期待値とする。きちんと書くと$$ \bar{y}_{t_0} = \frac{1}{7} \sum_{i=t_0 – 7}^{t_0 – 1} y_i $$ $$ s^2_{t_0} = \frac{1}{7-1} \sum_{i=t_0 – 7}^{t_0-1} (y_i – \bar{y}_{t_0})^2 $$ ね。で、$$ C_{t_0} = \frac{y_{t_0} – \bar{y}_{t_0}}{s_{t_0}} $$ を求める。
　アウトブレイクなしという\(H_0\)の下で、\(C_{t_0} \sim N(0, 1)\)と仮定できる。[えええええ！？正規近似すんの！？ そうか、疫学の世界では標本サイズが大きいからね…]
　上界\(U_{t_0}\)を、標準正規分布の\(1-\alpha\)分位点 \(z_{1-\alpha}\)を使って$$ U_{t_0} = \bar{y}_{t_0} + z_{1-\alpha} s_{t_0} $$ として、\(y_{t_0} \gt U_{t_0}\)になったら警告する。
　[ま・じ・か。要はシューハート管理図じゃんか]

　EARS法C1, C2, C3は直近の情報しか使わない。短期間のデータしかないとか、カウントが一定だと期待できる場合にはこれでいいんだけど、トレンドや季節性がある場合位は困る。

3. いくつかの文脈におけるサーベイランス手法
　では、アルゴリズムを紹介しよう。現状では以下の方法を実装している: earsC, farrington, farringtonFlexible, bayes, cdc, boda, bodaDelay[本文中に説明がない], glrnb, glrpois, cusum, rogerson, outbreakP, categoricalCUSUM, pairedbinCUSUM.

3.1 カウントデータのための万能手法
　farrington, farringtonFlexibleを紹介する。前者はFarrington et al.(1996 JRSS-A)。後者はその改良版でNoufaily et al.(2012 Stat.Med.)。
　まず、\(\mathbf{y}_{t_0}\)にoversisperseつきのポアソンGLMモデル(対数リンク)をあてはめる。\(E(y_t) = \mu_t, \log \mu_t = \alpha + \beta t, Var(y_t) = \phi \mu_t \) ただし\(\phi \geq 1\)、というモデルである。[びびってはいかん。モデルの右辺にトレンドをいれているだけだ。季節性を入れる、参照期間を少し削る、などいくつかの話題が載っているが、メモ省略]
　次に、カウントの期待値\(\mu_{t_0}\)を予測する。で、予測区間の上界\(U_{t_0}\)を求める。正規近似する方法などいくつかご用意している[メモ省略]。
　最後に、\(y_{t_0}\)と\(U_{t_0}\)を比較し、警告の有無を判断する。
　なお、GLMのフィッティングにはいろいろ工夫があって…[メモ省略]
　この方法と似た方法として以下がある。トレンドとか、過去のアウトブレイクの割引とか、GLMのフィッティング時の工夫とかがない。あまりお勧めしません。

• bayes: 共役事前分布を用い、負の二項分布のパラメータを推定する
• cdc: 週次データを4週ごとにまとめて正規分布を使う

3.2 ベイジアンによる改良
　bodaを紹介する。これはfarrington, farringtonFlexibleのベイズ版で、ベイジアンGAMを使う。なかではINLAパッケージを呼んでいる。
　[面白そうなんだけど、時間の都合でパス]

3.3 1時点検出を超えて
　いろいろ工夫はあるものの、Farrington法は結局は1時点検出である。アウトブレイク初期に徐々にシフトが起きたときには検出できない。そういうのは統計的プロセス管理でいうところの管理図で見つけられるんだけど、今度はoverdispersionや季節性に対応できない。
　Hohle & Paul(2008 Comp.Stat.DataAnal.)は伝統的なサーベイランス手法と統計的プロセス管理の合成手法を提案している。glrnb, glrPoisがそれである。

　管理図では、ICとOCの2状態を考えて、各時点で尤度を比較する。
　共変量を\(\mathbf{z}_t\)とし、IC分布を\(f_{\theta_0} (y_t | \mathbf{z}_t) \)とする。これを二項分布(glrnb)ないしポアソン分布(glrPois)、対数リンクのGLMで推定する。IC平均は[… \(\log \mu_{0,t}\)の式が示されている。調和回帰で季節効果を表現していてめんどくさいのでメモは省略するけれど、要するに、切片項、トレンド、季節項のみのモデルである]
　[…残りのメモも省略するけど、要するにこういう話だ。フェイズ1でパラメータを学習し、フェイズ2でサーベイランスする。OC平均は\(\mu_{1,t} = \mu_{0,t} \exp(\kappa)\)ないし\(\mu_{1,t} = \mu_{0,t} + \lambda y_{t-1}\)とモデル化する。で、各時点で尤度比ないし一般化尤度比を求める。式を見た感じだと、どうやらCUSUMみたいに累積するようだ。で、警告を判断する]

　ほかに、CUSUM法としてcusum, rogersonを用意している。Hohle & Paul(2008)が比較しているのでみよ。
　outbreakPというのも用意している。これは定数レベルから単調に増大するような変化(たとえばインフルエンザ流行の始まり)を検出するセミパラメトリックな方法である。

3.4 カテゴリカルデータの監視方法　
　[時間の都合でパス。categoricalCUSUMという方法の紹介である。末尾に数行だけpariedbinCUSUMが出てくる]

3.5 その他のアルゴリズム
　表をつくったのでみろ。

　アルゴリズムの選び方について。以下の点に注目するとよい。履歴データをどのくらい持っているか(EARSは直近だけで大丈夫だがFarringtonはたくさん必要)。そのアルゴリズムが履歴データのうちどこを使っているか。一時点以上の検出か、もっと長いシフトの検出か。ファインチューニングはどのくらい大変か。

4. サーベイランスを通常の監視に組み込む
　[面白そうなんだけどな… いまちょっと読んでる時間がない]

5. 結論
　このパッケージはほかにも、モデリングとか、ナウキャスティングとか、バックプロジェクションとかの機能を持っている[それぞれreferenceを挙げている]。
　今後の課題: 階層構造のあるデータ、報告遅延の調整。
　このパッケージはもともと公衆衛生の文脈で作られたものだが、ファイナンス、職場の安全性監視、環境サーベイランスなど多方面でご利用いただける。

　関連パッケージの紹介。統計的プロセス管理ならspc, qccがある。構造変化ならstrucchange, カウント時系列のGLMならtscountがある。疫学モデリングはEpiEstim, outbreaker, Outbreak-Toolsがある。Rじゃないけど、異常検出のSaTScanというソフトもある。
　云々。