読了：Sonesson & Bock (2003) 公衆衛生における時系列監視

Sonesson, C., Bock, D. (2003) A review and discussion of prospective statistical surveillance in public health. Journal of Royal Statitical Society, A. 166, 5-21.

　公衆衛生分野での時系列監視のレビュー。やれやれ、なぜこんな面白くもない話を勉強しているのか… (すいません、頭の悪い奴のひがみです)
　とはいえ、品質管理の話よりはとっつきやすい。あとファイナンスね！ いったいなんなのあれ。疫学や品質管理の専門家は人類にとって必要だけど、ファイナンスの専門家なんて全員拉致して離島に閉じ込めてしまえば世界はかえって平和になるのではないだろうか。(ごめんなさい、貧乏人のそねみです)

1. イントロダクション
　[…]
　罹患数などの時系列データの使い道として後ろ向きretrospectiveなのと前向きprospectiveなのがある。後者の場合、インシデンスが増えているかどうかをシーケンシャルに判断しなければならない。統計的観点からいうと実は後ろ向き分析より難しい。本論文は公衆衛生分野における前向きな統計的サーベイランス手法に焦点を当てる。
　[品質管理におけるサーベイランスとのちがい…]
　[先行するレビューとのちがい…]

2. 統計的サーベイランスの一般的概念
　[Frisen(2009)と定式化がびっくりするくらい似ている。よく見たら同じ同じ研究チームなのね。スウェーデン・イェーテボリ大学経済学部だ]
　統計的サーベイランスとは、確率過程\(X = \{X(t); t = 1,2,\ldots\} \)をオンライン監視して、未知の時点\(\tau\)における過程の重要な変化をできるだけ早く正確に検出するという問題である。決定時間\(s\)のそれぞれにおいて、システムがin-control状態 \(D(s)\)なのかout-of-control状態\(C(s)\)なのかを弁別したい。そのために、\(X_s = \{ X(t); t \leq s\}\)を使って警告集合\(A(s)\)をつくる。もし\(X_s ^in A(s)\)だったらシステムが\(C(s)\)ですという警告を発するわけである。これは通常、警告関数\(p(X_s)\)と制御限界\(g(s)\)を用いて実現される。警告時間を\(t_A\)として$$ t_A = \mathrm{min} \{ s; p(X_s) \gt g(s)\} $$ と書ける。

　\(D(s), C(s)\)は場面によっていろいろだが、もっとも多いのは\(D(s) = \{ \tau \gt s\}, C(s) = \{ \tau \leq s\} \)である。
　検出すべき変化も場面によっていろいろである。\(X\)の分布のパラメータの変化に関心が持たれることが多い。多くの研究はある定数レベルから別の定数レベルへの変化を扱っているが、徐々に起きる変化とか、線形な変化とか、指数的増大について扱っている研究もある。

3. サーベイランス手法の評価
　まず、誤警告の分布をin-control runの長さの平均$$ ARL^0 = E[t_A | \tau = \infty] $$ で要約することがある。また、誤警告確率$$ P(t_A \lt \tau) = \sum_{t=1}^\infty P(\tau = t) P(t_A \lt \tau | \tau = t ) $$ を使うこともある。しかしこれは\(\tau\)の分布についての仮定が必要になる。幾何分布だと仮定されることが多いが、そうすると各時点のシフト強度が一定だと仮定していることになる。
　さて、サーベイランス手法の評価においては誤警告と警告遅延の短さとのトレードオフに直面することになる。仮説検定と同じように考えれば、\(ARL^0\)なり誤警告確率なりを第一種過誤だと思って固定して比べることになる。
　真の変化に関する評価の方法はたくさんある。一番よく用いられているのはout-of-control runの長さの平均$$ ARL^1 = E[t_A | \tau = 1] $$ である。しかしこれは、製造プロセスでは自然だが[あっ、そういう発想なのか、ようやく腑に落ちたよ]、公衆衛生ではそうでもない。代わりに用いられている指標として条件付きの期待遅延 $$ CED(t) = E[t_A – \tau | t_A \geq \tau = t] $$ がある。\(\tau\)の分布を仮定すれば、期遅延 $$ ED_\tau = \sum_{t=1}^\infty P(\tau = t) P(t_A \geq t) CED(t) $$ を用いることもできる。
　遅延期間\(d\)が許容できる場合には、成功検出確率 $$ PSD(d, t) = P(t_A – \tau \leq d | t_A \geq \tau = t) $$ を使うこともできる。
　サーベイ手法によって誤警告の分布が異なる。従って特定の時点での警告に占める誤警告の割合も異なる。言い換えると警告の信頼性が異なるわけである。その指標として予測値$$ PV(t) = P \{ C(t) | t_A = t\}$$ が使える。願わくば時間を通じて一定であってほしい。
　[以下、referenceがたくさん載っている]

4. 発生率増大の検出
　罹患数についてはポアソン過程が仮定されることが多い(4.1節で扱う)。もっと複雑な時間依存仮定が仮定されることもある(4.2節で扱う)。
　もうひとつの重要な仮定として、ベースライン率が既知かどうかがある。ベースライン率に季節効果があるという場合もある。要するに、\(C(s), D(s)\)の特徴づけが必要なわけである。

4.1 ポアソン過程の強度の変化の検出
4.1.1 イベント間の時間間隔を使う場合

　統計的プロセス管理の標準的手法としてCUSUMとEWMAがある。

  • CUSUMでは、警告統計量は観察と期待の差のCUSUMに基づく[cumulative sumのことね]。\(\tau = t\)のときの時点\(s\)の尤度比を\(L(s,t)\)として、$$ t_A = \mathrm{min}[s; \mathrm{max}_t \{L(s,t)\} \gt K] $$ となる。これはミニマックス最適であることが知られている。つまり、変化が起きる前の観察結果が最悪であったときの条件付き期待遅延の最大値を最小化している[ああ、ここでいうミニマックスってそういう意味だったのか…]。
  • EWMAでは、警告統計量は観察のEWMAに基づく。$$ Z_s = \lambda \sum_{t=1}^s (1-\lambda)^{s-t} X(t) $$ 警告限界はふつう定数とされる。

　指数分布する変数にたいしてCUSUMやEWMAを使った例は公衆衛生の文脈ではみあたらないが、他の分野ではみられる。

　先天性奇形のサーベイランスの分野では、奇形児の出産の間隔をその間に生まれた健常児の人数で測る。いろんな手法があって…[パス]

4.1.2 固定期間におけるイベント数を使う場合
[ここ、関心があるので細かくメモする]
　固定期間におけるイベント数が記録されている場合、過程についての情報は失われ、サーベイランス手法は過程における変化をできるだけ早く検出するという点に関して最適でなくなる。従って、固定期間をつかう理由は報告システムの実務的制約だけである。

  • よく用いられているのはポアソンCUSUM法である。各期間のイベント数を期待値と比べ、偏差の累積和で警告統計量をつくる。また、条件付き尤度比の最大値を警告統計量にすることもある。
  • ポアソンCUSUM法についての一般的レビューはLucas(1985 Technometrics)をみよ。ポアソンCUSUM法は先天性奇形のサーベイランスにもつかわれていて、前項の手法と比較した研究がたくさんある[…中略…]。ポアソンCUSUMのほうがよいといわれている。
  • Lie et al.(1993): 乳児がダウン症である確率の系列的二項尤度比検定を提案。警告限界は(alphaレベルではなくて)\(ARL^0\)を固定するように決められている。ポアソンCUSUMと比べると、\(ARL^1\)の観点からは系列的二項尤度比検定のほうが良いと主張されている。
  • Radaelli(1992): ポアソンCUSUMを累積スコア法と比べている[なにそれ]。
  • Rossi et al.(1999 Statist.Med.): ポアソンCUSUMをいろんな正規近似とくらべている。
  • Praus et al.(1993 Statist.Med.): ポアソンCUSUMを薬の副作用の上市後サーベイランスに使っている。
  • Hutwagner et al.(1997): ポアソンCUSUMをサルモネラのアウトブレイクに適用している。
  • Bjerkesal & Bakketeig(1975): 先天性奇形の事例数に対してポアソン・シューハート法を使った初期の例。シューハート型の方法とは、直近の観察のみを使い、$$ t_A = \mathrm{min} \{s; L(s,s) \gt K\} $$ とする方法である。一般的に効率的でなく、\(C(s) = \{\tau = s\} \)と定式化できるときだけに使える。つまり、即時検知だけに関心があるか、検出すべき変化が大きいときだけである。[そうか… シューハート管理図ってのも基本的には即時検知の方法なんでしょうね。運用面ではいろいろ工夫があるかもだけど]

4.1.3 移動窓のなかの過程を観察する方法
　後ろ向き分析の設定での使用についてはStroup et al.(1989, 1993)を参照。前向きの場合として[…reference…]。
　この方法は最適でない。たとえば、固定長の2つの連続した移動窓を比べて、徐々に起きている変化を検知するのは難しい。移動窓をつかうと過程についての情報が失われてしまう。[…] ベースライン率が全く未知の場合には使わざるを得ないと思うかもしれないけど、連続時点の差をサーベイランスするというような他の方法がある。Frisen(1992 Statist.Med.)をみよ。
　[事例。メモ省略]

4.1.4 尤度比法
　ポアソン過程では固定期間のイベント数よりイベント間の時間を使ったほうがよい。しかし、サーベイランス手法の性質は警告統計量と警告限界によっても決まる。要はシステムの望ましい特性に応じて決めるべきである。
　望ましい性質を表現する際には最適性基準が用いられることが多い。では最適性とはなにかというと、これがなかなか難しいのだが、誤警告の確率を固定した時の期待遅延の最小化というがひとつの自然な考え方であろう。[…]
　\(D(s) = \{\tau \gt s\}, C(s) = \{\tau \leq s\}\)で、過程の強度が\(\lambda_0\)から\(\lambda_1\)に変化するという場合について考えよう。期待遅延を最小化するのは尤度比法である。$$ t_A = \mathrm{min} \{ s; P(\tau \leq s | X_s = x_s) \gt K\} $$ $$ = \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)} \frac{K}{1-K} \right] $$ $$ = \mathrm{min} \left[ s; \sum_{t=1}^s \frac{P(\tau = t)}{P(\tau \leq s)} L(s,t) \gt \frac{P(\tau \gt s)}{P(\tau \leq s)} \frac{K}{1-K} \right] $$ と書ける。

　尤度比法の弱点は\(\tau\)の分布が必要だという点である。幾何分布が用いられることが多い。シフトの強度が低ければShiryaev-Roberts法で近似してもよい。$$ t_A = \mathrm{min} \left\{ s; \sum_{t=1}^s L(s,t) \gt K \right\}$$ 変化時点について無情報事前分布を仮定したと考えることもできる。

　尤度比法なりShiryaev-Roberts法なりをポアソン過程に拡張するのは容易である[Appendixで導出している]。これらの方法はイベント間間隔が時間で表現されていても件数で表現されていても使える。誤警報の確率を固定したとき、他の方法より期待遅延が短い。とはいえ、ある疫学的な文脈でいったいなにがも望ましい特性なのかということについては十分に検討する必要がある。S-R法の場合、予測値が一定であるという特性が満たされるのは、正規分布の平均がシフトするときである。ポアソン過程のシフトに関しても成り立つだろうが、他の場合どうかはわからない。

　手法の違いを理解する方法のひとつは、警告統計量における条件付き尤度比の使われ方に注目することだ。CUSUM法では最大値、Shewharts法では直近の値をみる。尤度比法では変化点の分布で重みづける。S-R法では均等な重みを使う。

4.2 時間依存性のある過程
　罹患数時系列が季節性や自己相関などの時間依存性を持っていたらどうするか。正しくモデル化できればそこからの逸脱を変化の指標と捉えることができる。ふつう定常性の仮定はいらない。

  • Watier et al.(1991): ARIMAモデルで、予測区間の上側を警告閾値にする。
  • Nobre & Stroup (1994): 予測誤差を使って逸脱検出の確率指標関数を求める[?]
  • Farrington et al.(1996): 伝染病のアウトブレイクを見つけるのに回帰アルゴリズムを使う。ベースライン率の予測区間を閾値にするので、シューハート型の指標と言える。
  • Williamson & Hudson (1999): ARIMAモデルをいろんな病気についてつくり、1期先予測誤差をサーベイランスに使う。
  • Vanbrackle & Williamson (1999): いろんな指標のARLを比較。
  • Smith & West(1983 Biometrics): 状態空間モデルとカルマンフィルタで状態のオンライン事後確率を求める。[おおお、やっと出てきた。なんで状態空間モデルを使わないのかなと不思議だったのである。メモは省略するけど、他にも論文が挙げられている: Smith et al.(1983 Statistician), Gordon & Smith (1990 JASA), Whittaker & Fruhwirth-Schnatter (1994 App.Statist.), Schlain et al.(1992 Statist.Med.), Schlain et al.(1993 Statist.Med.), Stroup & Thacker (1993 Epidemiology). 読むんなら最後の3件が楽しそう]

5. 空間的サーベイランス
[2頁。いまちょっと関心ないのでパス]

6. 考察と結語
[略]
——————
　勉強になりましたですー。
　自分の抱えている問題に引き付けて考えると、注目している変数について正しくモデル化することもさることながら、何をもって逸脱とするかという点を明確にすることがポイントだなと思った。ユーザの立場からすると、なんでもいいからなんかあったら教えてね、という感じであって、たとえば急増したら警告してほしいし、過去12期にわたって徐々に増えていてもそれはそれでそう教えてほしいんだけど、こういう異なる課題をごっちゃにすると良い警告システムが作れないのである。