« 読了：岩崎 et al. (2015) ツイッター上の炎上を予知する | メイン | 読了: 紀要論文に垣間見るさまざまな人生 »
2015年2月 4日 (水)
調査だか実験だかの参加者に刺激 P と Q を提示し、それぞれに対してなんらかの反応を得た。反応の分布を比べたとき、PとQのあいだに統計的に意味のある差はあるか。
反応が量的ならば、分布の平均に注目し、いわゆる「対応のある二水準」のあいだの平均の差の検定を行うことが多いだろう。手法としていちばんポピュラーなのは t 検定である。反応が二値ならば、(仮に反応を{0,1}とコード化したとして)各刺激に対する1の割合に注目し、「対応のある二水準」の間の割合の差の検定を行うことが多いだろう。McNemar検定が広く用いられている。
では、刺激が3つ以上あったらどうか。反応が量的ならばANOVAを用いることが多い。反応が二値の場合はCochranのQ検定を行うことが多いだろう。しかしこれらはオムニバス検定、すなわち、K個の水準のあいだの「どこか」に差があるといえるかどうかを調べている検定である。「どこ」に差があるのかを調べるためには別の手法が必要だ。
K個の水準から2つを取り出す全てのペアに関心があるとしよう。反応が量的である場合については、多重比較法と呼ばれる手法がたくさんあり、良い解説もたくさんある。問題は反応が二値の場合だ。いちばん簡単に思いつくのは、McNemar検定をすべてのペアについて繰り返す方法だが、ファミリーワイズの第一種過誤(Type I FWE)が増大する。そこで思いつくのは、Bonferroni法かなんかでp値を調整しちゃう、という手である。ちょっといま手元にないのだけれど、この手法を推奨している参考書が多いんじゃないかと思う。しかし、ちょっと不思議ではありませんか？ Bonferroni型の調整手法は汎用的な手法であって、個別の課題においては必ずしも最適でない。反応が量的な場合なら多重比較法はよりどりみどりなのに、反応が二値になると急に品揃えが悪くなっちゃうわけ？
手元の本を探してみたところ、困ったときの助け神、森・吉田本には、McNemar検定とRyan法を併用しろ、と一行だけ書いてある。森・吉田本は他の箇所でもRyan法推しなのである。Ryan法ってなに？ わたくし、それってTukey-Welsch法の別名なんじゃないかと疑っているのだが、詳しいことはわからない。SASやSPSSのANOVAの機能には多重比較法としてREGW法というのが搭載されてるけど(REGWQとか)、あれとの違いもよくわからない。
泣く子も黙る権威、AgrestiのCDA本は、章末の注でWestfall, Troendle, & Pennello (2010) を参照せよと書いている。私はこの本のいうことは頭から信じることに決めている。だって高かったんだもん。私費で買ったんだもん。
正直言って、細かい話だと思う。もし私の好き勝手が許されるならば、検定しないで済ませる方法を全力で考える。いかんともしがたい場合は、McNemar検定を繰り返し、Type I FWEが増大してますけどそれがなにか？ と開き直る方法を考える (←冗談ではない。これがもっとも正しい態度である場合も少なくない)。それがだめでも、細かいことを考えるのはうんざりなので、反復測定データのロジスティック回帰だかGLMMだかGEEだかに持ち込み、推定したパラメータをTukey法かなにかでツルッと多重比較することを考える。
しかし、娑婆は細かい話であふれておりまして...
Westfall, P.H., Troendle, J.F., & Pennello, G. (2010) Multiple McNemar Tests. Biometrics, 66(4), 1185-1191.
そんなこんなで、仕事の都合で読んだ。Agresti先生ご紹介の、McNemar検定を多水準に拡張して多重比較する手法についての論文。単にBonferroni法やSheffe法で調整するよりも良い方法をご提案します、とのこと。具体的には、検定は正確McNemar検定、多重比較の理屈はHolm法だと思う。
多重比較には疎いので、よくわからない箇所も多いんだけど...
まずは二群のケース。
IIDな二変量ベルヌーイ試行系列$(Y_{i1}, Y_{i2})$について考える。平均ベクトルを$(\theta_1, \theta_2)$とする。帰無仮説 $H: \theta_1 = \theta_2$について正確検定しよう。
$Y_1$と$Y_2$の同時確率分布を$\theta_{00}, \theta_{01}, \theta_{10}, \theta_{11}$とする。以下、$H$の下での非対角セルの確率を$\theta_d = \theta_{01} = \theta_{10}$とする。
$H$の下で、2x2クロス表の非対角セルの観察頻度の和を $N_d = N_{01} + N_{10}$としよう。$N_d=n_d$のときの$N_{01}$の条件分布は、言うまでもなく$B(n_d, .5)$だ。$N_{01}$の実現値が$n_{01}$だったとしよう。$B_{n_d, .5} \sim B(n_d, .5)$として、p値は
$p(n_{01}, n_d; upper) = P(B_{n_d, .5} \gt n_{01})$
$p(n_{01}, n_d; lower) = P(B_{n_d, .5} \lt n_{01})$
$p(n_{01}, n_d; two) = 2 \min \{ p(n_{01}, n_d; upper), p(n_{01}, n_d; lower) \}$
どれでもいいけど、とにかくMcNemar検定において、p値は$n_{01}$と$n_d$の関数として決まるわけである。
多重比較に拡張します。
仮説を$H_1, H_2, \ldots, H_m$とする。仮説$l$に対応するデータを$\{ y_{i1}^{(\ell)}, y_{i2}^{(\ell)} \}, i=1, \ldots, n^{(\ell)}$とする。仮説のあいだでデータが重複しているかどうかは問わない。個々の仮説について正確McNemar検定をしたときの、$m$個を通したFWERを制御したい。
もし$m$個の仮説のなかに真の仮説がなかったら、タイプIエラーもないことになるから、まあひとつは真の仮説があるとしよう。真の仮説の集合を$I$として、
FWER($I$) = sup $P$(Reject $H_l$ for some $\ell \in I$)
ただしsupは上限を表す記号。
実際にはどの仮説が$I$に属しているのかわからないわけだけど、それがどれであってもFWER($I$) $\leq \alpha$ となるようにしたいものである。こういうのをstrong controlという。これに対して、すべての仮説が真であったときにFWER($I$) $\leq \alpha$ となることを weak controlという。以下ではstrong controlについて考える。
仮説の論理積を表す仮説(つまり「どこにも有意差がない」仮説) $H_I = \cap_{\ell \in I} H_\ell$の検定統計量について考える。ある仮説$\ell$のp値を$p(N_{01}^{(\ell)}, N_d^{(\ell)})$として、
$T_I (N^+_I) = \min_{\ell \in I} p_\ell (N_{01}^{(\ell)}, N_d^{(\ell)})$
$N_I = \{N^{(\ell)}_d | \ell \in I\}$が与えられれば、この検定統計量について、Type I Error率が$\alpha$以下になるような臨界値$c^\alpha_I(N_I)$を定めることができる由(詳細略)。
検定統計量の棄却域を決めるより、対応するp値をつくったほうが楽である。すなわち:
$\tilde{p}_I (n^+_I) = \sum_{\ell \in I} P_{H_\ell} \{ p_\ell (N_{01}^{(\ell)}, N_d^{(\ell)}) \leq T_I (n^+_I) | N_d^{(\ell)} = n_d^{(\ell)} \}$
このp値を使ってdiscrete Bonferroni-Holm法を用いる由。知識不足でちょっと理解できないのだが、すべてのp値を大きい順に並べてステップダウン、仮説集合の p 値として上の式を使う... ということだろうか。ほかにブートストラップを使った手法も提案しているが、そちらは力尽きたのでパス。
後半は数値例とシミュレーション。パス。
どうでもいいけど、大文字のアイと小文字のエルを併用するのはやめてくださいな、先生。しょうがないので後者はわざわざ$\ell$と書いた。
読んだはいいけど、実際の計算はなにでどうやればいいかしらん。SASでもRでも、まるごとの実装はみあたらなかった。WestfallさんってSASのMULTTESTプロシジャやRのmultcompパッケージの開発者だと思うんだけど。見落としてるのかな。
論文：データ解析(2015-) - 読了：Westfall, Troendle, & Pennello (2010) 多重マクニマー検定