読了：Frisen (2009) 時系列の監視

Frisen, M. (2009) Optimal Sequential Surveillance for Finance, Public Health, and Other Area. Sequential Analysis, 28, 310-337.

　仕事の都合で読んだやつ。時系列監視についての解説論文。
　著者のFrisenさんはこの分野の有名な人だと思う(前に仕事で翻訳をやったことがある)。google様いわく、本論文の被引用件数は82。招待論文という位置付けで、10人のコメントと返答がついている。っていうか、Sequential Analysisなんていうジャーナルがあるのね。掲載論文を読むのははじめてだと思う。

　時系列の監視の話って、いろんな分野に専門家がいて全然違う用語を使うのでどうもとっつきにくく、あまり関わりたくない話題のひとつである。とはいえ、それをいうならそもそも万物に関わり合いを持ちたくないし(枕と布団を除く)、仕事とあらば好き嫌いはいってられない。がんばりましょう、と気合をいれて…

1. イントロダクション
　系列的サーヴェイランスの目的は、データを生成するプロセスにおける重要な変化をタイムリーに検出することである。[…]

2. 統計的サーベイランス問題
2.1 特徴づけ
　観察された量だか統計量だか、はたまたスカラーだかベクトルだか知らないが、とにかくサーベイランスの対象である過程を\( X = \{X(t): t =1, 2, \ldots, \} \)とする。変化が一度だけ、時点\(\tau\)で起きるとしよう。分布は変化の前は族 \( f^D \)に従い変化の後は族 \( f^C\)に従うとする。
　私たちは決定時点\(s\)において、出来事\(C(s)\)と\(D(s)\)を区別したい。たいていは、\(C(s) = \{\tau \leq s\}, D(s) = \{\tau \gt s\}\)である。
　システムの状態を決めている(おそらくは確率的な)過程のことを\(\phi(t)\)とする。たとえば分布のパラメータね。たいていの研究は、\(\phi(t)\)がまず\(\phi^0\)で、時点\(\tau\)に\(\phi^1\)に変わる、と想定している。
　私たちは観察\(X_s = \{X(t); t \leq s\}\)を使って警告統計量\(p(X_s)\)と制御限界\(G(s)\)をつくる。警告が発せられる時点は$$ t_A = \mathrm{min} \{s; p(X_s) \gt G(s)\} $$ である。

　サーベイランスには、洪水の警告のようにパッシブなものと(警告したからといって洪水は止まらない)、アクティブなものがある(帝王切開手術中の胎児の監視とか)。両者のあいだで過誤率の意味が違う。Frisen & de Mare (1991 Biometrika)をみよ。
　\(\tau\)は確率変数としてみるのが適切な場合と、未知の決定論的な値とみるのが適切な場合がある。前者の場合、確率\(\pi(t) = P(\tau = t)\)を考えることができる(ベイジアンなら事前分布とみなせる)。[…]
　サーベイランスにおける統計的推論は仮説検定とちがうという点に注意。系列的仮説検定であっても仮説は固定されている。しかしサーベイランスには固定された仮説がない。分娩を監視しているとして、胎児の首にへその緒が巻き付いているという仮説が棄却されたからといって、監視をやめるわけにはいかないでしょ？ [このたとえ話、前に読んだな。著者お気に入りの事例なのであろう]

2.2 評価
　有意水準、検定力、特異性、敏感性、などなどのご存じの指標による評価はいっけん便利にみえるけど、サーベイランスの文脈では解釈が難しい。[…] ではどういう指標を使うのが良いかというと…

  • false alarm.
    • 一番よくつかわれるのは\( ARL^0 = E(t_A|D) \)である[average run lengthの略]。平均じゃなくてメディアンをつかうこともある。
    • 理論研究では \( PFA = P(t_A \lt \tau)\)を使うこともある[false alarm prob.]。
  • 警告の遅延。
    • 一番よくつかわれるのは\( ARL^1\)、すなわち(変化がサーベイランス開始時点に起きていた時に)変化を検出するまでの平均時間だが、その定義は手法によって異なるし、うまく定義できないこともある。
    • \(\tau\)から\(t_A\)までの遅延の期待値 $$ ED(t) = E[max(0, t_A – t) | \tau = t]$$を使えという話もある。
    • \(ED(t)\)は\(t\)が増えるとゼロに近づくので、条件づけた期待値$$ CED(t) = E[t_A – \tau | t_A \geq \tau = t] = \frac{ED(t)}{P(t_A \geq t)} $$ を使えという話もある。ちなみに\(ARL^1 = CED(1) + 1 = ED(1) + 1\)である。
    • たいていの手法では\(CED(t)\)は定数に収束する。収束先、つまり安定状態平均遅延をSADTという。\(\tau\)がすごく大きい場合についてのみ考えているわけで、つまり\(ARL^1\)の正反対である。[なるほど]
    • 警告がそれ以上遅れちゃうと元も子もない、という定数\(d\)を考えて$$ PSD(d, t) = P(t_A – \tau \leq d | t_a \geq \tau = t) $$ を使うこともある[prob. successful detection]。重要な変化を素早く見つけたい場合は\(d\)を小さく、長期的な検出能が大事な場合は大きくする。
  • 予測値。$$PV(t) = P(\tau \leq t_A | t_A = t)$$ つまり、警告が出た時に変化が起きている確率。逆に警告が出ていないときの予測値を考えることもできる。この指標には、false alarmsのリスクとか、期待される遅延とか、変化についての信念とかが含まれている。

3. 最適性
3.1 最適性基準

  • ARL最適性。よく用いられる最適性指標は、\(ARL^0\)を固定したときの最小\(ARL^1\)である。これはつまり、全時点が\(f^C\)に属しているという仮説と全時点が\(f^D\)に属しているという仮説を弁別しようとしているわけで、全観察に等しいウェイトをつけるのが正解だということになる。しかしサーベイランスというのは直近の観察に最も重みを付けるべきなのであって、つまりこの最適性指標を使っていると初期の変化に関心を集中させていることになる。記述としてはかまわないけど、最適性基準としてはおかしい。[…]
  • SADT最適性。これは逆に、直近の変化を強調する指標となる。
  • 最小期待遅延。
    • 期待遅延は\(\tau\)に依存するので、平均した\(\sum w(t) CED(t)\)を用いる。結局\(\tau\)と\(t_A\)の分布を通じた\(ED(i)\)の期待値\(ED = E[ED(i)]\)となる[なぜ急に\(i\)になったんだろうか]。PFAを固定してこれを最小化することをED基準という。
    • Shiryaev(1963)は以下を提案している。効用関数\(u(\tau, t_A)\)を考える。\(t_A \lt \tau\)ならば\(h(t_A-\tau)\)とし(たいていは定数\(b\)とする。false alarmのコストはふつう時点を問わないから)、そうでなければ\(a_1(t_A – \tau) + a_2\)とする。で、その期待値\(U = E[u(\tau, t_A)]\)を考える。結局$$ U = b P(t_A \leq \tau) + a_1 ED + a_2$$ となる。[なるほどね、これが一番腑に落ちるなあ]
  • ミニマックス最適性。\(\tau\)の分布を使わないというのがポイント。多くの理論研究がこの指標をつかっている。[… 簡単に説明してあるんだけど、私の理解の範囲を超えるのでメモ省略。sup ess sup とかって言われても困るよ…]

3.2 最適手法
　[…] サーベイランスの最適性についての多くの結果は、変化前と変化後のデータがiidだと仮定している。本節でもそうする。[…]
　サーベイランス手法の多くは部分尤度比の結合として表現できる。\(\tau\)を固定したときの尤度比は$$ L(s, t) = \frac{f_{X_s}(x_s | \tau = t)}{f_{X_s}(x_s | D) } $$ である。式の要素をどうするかは状況によるけれど、もっとも多いのは、独立な正規分布でシフトは1ステップだけだとするやりかたである。

  • 完全尤度比法(LR)。最小期待遅延基準と、効用関数の幅広いクラスに対して最適である。完全尤度比が基準を超えたら警告する。\(K\)を定数として$$ t_A = \mathrm{min} \left\{ s; \frac{f_{X_s}(x_s | C(s))}{f_{X_s}(x_s | D(s))} \gt \frac{P(\tau \gt s)}{P(\tau \leq s)} \times \frac{K}{1-K} \right\} $$ 部分尤度比をつかって書き直すと、密度\(\pi(t)\)に比例した\(w(s,t)\)をもってきて$$ t_A = \mathrm{min} \left\{ s; \sum_{t=1}^s w(s,t) L(s,t) \gt G(s) \right\} $$ となる。\(G(s)\)は警告限界だが時間に依存する。事後分布を使って$$ t_A = \mathrm{min} \left\{ s; P(C(s)|X_s = x_s) \gt K\right\}$$ と書いてもよい。なお、この手法のことをベイズ法ということがあるが、\(\tau\)の分布が事前分布といいうるようなものかどうかは場合による。
  • Shiryaev-Roberts法。単純に、\(\sum_{t=1}^s L(s,t) \gt G\)となったら警告する。\(\tau\)の事前分布が無情報な時のLR法だといってもよい。
  • Shewhart法[「シューハート管理図」のシューハートだと思う]。最後の観察だけ見て、\(L(s,s) \gt G\)のときに警告する。
  • CUSUM法。これを尤度比法と呼ぶこともあるんだけど、LRと混同しないこと。ミニマックス最適である。$$ t_A = \mathrm{min} \left\{ s; \mathrm{max}(L(s,t); t = 1,\ldots,s) \geq G \right\} $$
  • EWMA法(exponentially weighted moving average)。ターゲットの値を\(Z_0\)として(0にしてもよい)、$$ Z_s = (1-\lambda) Z_{s-1} + \lambda Y(s)$$ を警告統計量にする。\(\lambda\)は0より大, 1より小の定数ね。1にしたらShewhart法に帰着する。[よくわからん。\(Y(s)\)ってなによ？ \(X(s)\)のこと？ つまり、尤度比だとか一切考えず、観察の指数平滑移動平均が閾値を超えたら警告しますってこと？]
    　\(t_A\)は、漸近的に$$ t_A = \mathrm{min} \{ s; Z_s \gt L \sigma_Z\} $$ とする手と、そうではなくて正確なSDを使う手がある。
    　ここで大事なのは\(\lambda\)の設定である。ARL最適なのは\(\lambda = 0\)だけど、それはさすがにナンセンスである。
    　なおEWMA法は、部分尤度で直接には表現できないけれど、完全LR法における部分尤度の結合を線形近似したものだと捉えることができる[ふーん？？？]。

4. サーベイランスを必要とする分野
[ファイナンス、公衆衛生の2つについて詳しく紹介し、あとは名前といくつかのreferenceを挙げる程度]

  • ファイナンスにおけるサーベイランス:
    • ファイナンスのデータに基づく取引方略には多大な関心がもたれる。特に取引のタイムリーさが重要である。[…]
    • 方略の評価の際には投資のリターンに注目することが多い。ボラタリティに注目することもある。
    • 多くの場合とても複雑で…[非線形だとか、期待値だけじゃなくて分散に関心があるとか、多変量だとか、なんとか。超めんどくさいのでメモ省略]
  • 公衆衛生におけるサーベイランス:
    • [インフルエンザとかの流行検出について…]
    • 悪い治療のサーベイランスという問題もある。これはHarold Shipman事件がきっかけで注目された[イギリスの医者が患者を大量に殺していた事件だそうな。へええええ]。
    • 公衆衛生のサーベイランスのひとつの特徴は、ポアソン過程に関心がもたれるという点である。時間の代わりにポジティブ事例数が用いられるという点も特徴的である[健康な赤ちゃんが何人生まれて、病気の赤ちゃんが生まれて、健康な赤ちゃんが何人生まれて…というようなデータのこと]。変化も複雑なことが多く、未知のベースラインから徐々に変化することが多い。[Woodall(2006, J.QualityTech.), Woodall et al.(2008, JRSS-A)というのを読むとよさそう]
    • 空間的クラスタリングという問題もある。Kulldorff & Nagarwalla (1995)のソフトの登場によって注目されるようになった。
  • 患者の監視
  • 産業的品質管理
  • 経済学におけるサーベイランス。マクロ経済指標についてAnderson et al.(2005 J.Forecasting; 2006 J.App.Stat.)、顧客離脱についてPettersson(2004 Qualty&ReliabilityEng.Int)がある。
  • 環境のサーベイランス

5. 複雑な状況のための方法
5.1 自己相関を持つデータ
　自己相関がある場合、ひとつのアプローチは警告限界を調整してfalse alarmのリスクをコントロールすることで、もうひとつのアプローチは残差の過程を監視することだ。こうしゃはGARCHのような、尤度を明示的に描けない複雑な時系列モデルで有益である。[…]

5.2 離散分布
　公衆衛生では、ある特性(罹患)を持っている対象の割合とか数とかを日次ないし週次で観察することが多い。レアな疾患の場合は前回の発生からの時間が問題になることもある。
　ポアソン過程に関心を持たれることが多い。時間間隔ごとの件数だけが利用可能なこともある。件数が多いときには正規近似でいけるけど、そうでないことも多い。

  • Andersson et al.(2008 Stat.Methods.Med.Res.): インフルエンザのモデル
  • Woodall(1997 J.QualityTech.): いろんなデータにShewhart法に適用する
  • Lucas(1985 Technometrics): ポアソンCUSUM法。
  • Vardeman & Ray (1985 Technometrics): イベントの間隔が指数分布しているときの指数CUSUM法。
  • Chen(1978), Gallus et al.(1986), Radaelli(1992 J.AppliedStat.): ある数のイベントが生じたときにのみ記録されるようなデータや、時間間隔が閾値を超えたかどうかだけが記録されるようなデータの監視。
  • Gan(1991 J.Stat.Comp.Sim.), Borror et al.(1998 J.Qual.Tech.): ポアソンEMWA法。
  • Kennet & Pollak (1996 J.AppliedStat.): イベント間隔が指数分布するポアソン過程における強度の変化をShiryaev-Roberts法で検出する。
  • Sonesson & Bock(2003 JRSS-A): 固定期間のイベント数についてのLR法。

5.3 複雑な変化

  • 未知のレベルの間の変化。パラメータ(たとえば期待値)が未知の値に変化するという場合も多い。また、変化前についても十分な情報がないということもある。
    • ひとつの方法は、データをベースラインに対して不変な統計量に変換するというものである。それぞれの観察を、それまですべての平均からの偏差に変換するとか。
    • ベースラインからの距離をたとえば尤度比で測るという方法もある。尤度ベースの方法で最大尤度を測るのを一般化尤度比法(GLR)という[? ここよくわかんなかった…]
    • ベイジアンな方法もあって… [説明が短すぎてよくわからん]
  • 徐々に変化する場合。[集中力が途切れてきたのでほぼ逐語訳]
    　Yaschin(1993)は、CUSUM法とEWMA法を徐々に起きる変化の検出へと一般化している。Han & Tsung(2006)は動的変化の検出のためにCUSUM法のリファレンス・フリーな方法を提案している。
    　以下では単調な変化について論じる。変化点のサーベイランスについては4.1節, 4.2節で触れた。[…]
    　曲線の形状が未知である場合には、ノンパラメトリックな方法に関心が持たれる。前節で述べた、ベースラインとシフトが未知だという問題も回避できる。ここでもGLRアプローチが使える。Frisen(1994)は、順序制約の下での最大尤度推定を提案している。このアプローチは、曲線の形についてはノンパラメトリックだが、指数分布族に属する分布を使っているという点ではパラメトリックであり、よってセミパラメトリックであるといえる。Andersson(2002 J.AppliedStat.)をみよ。
    　\(\tau=i\)におけるある変化に関する部分最大尤度比は$$ \frac{\mathrm{max} f(x_s|C^i)}{\mathrm{max} f(x_s|D)} $$ とかける。\(\mathrm{max}\)は\(C^i\)ないし\(D\)における分布の族を通じての最大値である。もし時点\(i\)における変化の確率がウェイト\(w_i\)に比例するならば、完全最大尤度は$$ \frac{\mathrm{max} f(x_s|C)}{\mathrm{max} f(x_s|D)} = \sum_{i=1}^s w_i \times \frac{\mathrm{max} f(x_s|\mu = \hat{\mu}^{C_i})}{\mathrm{max} f(x_s|\mu = \hat{\mu}^D)} $$ と書ける。ただし、\(\hat{\mu}^{C_i}, \hat{\mu}^{D}\)は順序制約のもとでの最大尤度推定値である。これらのウェイトを使って、完全尤度比法の持つED最適性を近似するような性質を持つ完全最大尤度比法が可能になる。部分最大尤度比の組み合わせを変えて、3.2節で論じた他の最適性基準を近似することもできる。[… この後はreferenceが続く。うーむ、写経してみたけれどさっぱりわからない。徐々に変化するのを検出したいという話がなぜ上の式につながるのかが理解できていないのである。参ったなあ。上のAnderssonというのを読むのがいいのかな。ファイナンスの話は苦手だから、疫学分野のBock, Andersson, Frisen (2007, Biometrical J.)というのもよさそう]

5.4 多変量サーベイランス
[3ページ。疲れたのでメモは省略。とりあえず次元縮約するのが初手らしい。false alarmはFDRでコントロールしますとか、めんどくさい話が続く。後半は空間サーベイランスの簡単な紹介であった]

6. 将来の方向
　サーベイランスの研究はながらく品質管理に支配されていたんだけど、最近になって公衆衛生分野での研究が増えてきた(伝染病とかバイオテロとかのせいで)。今後はファイナンスも増えてくるだろう。
　[云々云々…]
——————
　もやもやしていた問題がいくつか解消した。レビューの徳といえよう。
　私が今考えている問題に近いのは公衆衛生のサーベイランスのようだ。手元のデータをあれこれいじるのも結構なことだが、やっぱり座学も大事だね。次はSonesson & Bock(2003 JRSS-A)を読むのがよさそうだ。