van Doorn, J., Ly, A., Marsman, M., Wagenmakers, E.J. (2020) Bayesian rank-based hypothesis testing for the rank sum test, the signed rank test and Spearman’s rho. Journal of Applied Statistics.
仕事の都合で慌てて読んだ奴。順位和検定、符号順位検定、順位相関の検定をベイズ・ファクターでやるにはどうすればよいかという解説論文。
1. イントロダクション
本論文は、パラメトリックなベイズ・ファクター(BF)をランク・ベースの手法に拡張する。
ランク・ベースの手法は、外れ値とか分布の仮定からの逸脱とかに強く、単調変換に対して不変である。また世の中には本質的に順序尺度であるデータも多いので便利である(リッカート尺度とか)。メジャーなのは、Mann-Whitney-Wilcoxon rank sum test (Wilcoxon rank sum test, Mann-Whitney-Wilcoxon U test. これは2標本t検定のランク・ベース版), Wilcoxon sighed rank test (対応のある t検定のランク・ベース版)、そしてSpearmanの\(\rho\) (ピアソンの相関係数のランク・ベース版)。
しかしこれらは頻度主義のパラダイムで開発されたので、BFの求め方はまだ提案されてない。ここで難しいのは、尤度関数がなにかがはっきりしないということだが、なんらかの潜在的なデータがあってそれは正規分布に従っていると考えればどうにかなる。
2. 一般的な方法論
関心あるパラメータを\(\theta\)とする。その事後分布は、尤度と事前分布の積に比例する: $$ \pi(\theta|data) \propto f(data|\theta) \times \pi(\theta) $$ さて、
問題は、ランク・ベース手法の場合尤度\(f(data|\theta)\)がわからんということだ。そこで潜在正規フレームワークを使う。
[ここで歴史のレビュー。メモ省略]
2標本の場合について考える。順序データのベクトルを\( (r^x, r^y) \)とし、それに対応する潜在正規スコアを \( (z^x, z^y) \)とする。潜在正規スコアの事後分布は $$ \pi(z^x, z^y, | x, y) \propto f(r^x, r^r | z^x, z^y) \times f(z^x, z^y| \theta) \times \pi(\theta) $$ 右辺第2項は潜在正規構造を表す。第1項は指標関数のセットで、観察された順位を観察されていない潜在正規スコアに結びつけている。
この式に基づき、観察された順序情報を入力にとって、\(\theta, z^x, z^y\)を生成するMCMCサンプラーを構築できる。
\( f(r^x, r^y | z^x, z^y) \) は指標関数で、潜在正規尤度 \( f(z^x, z^y| \theta) \) を切断することによって、潜在スコア\( z^x, z^y \)がデータにおける順序情報を保持することを保証する。[つまり、\( (z^x, z^y) \)の推定値の大小が、観察されている\( (r^x, r^y) \) の大小と矛盾する場合には、尤度関数\( f(z^x, z^y| \theta) \)がなにを抜かそうが問答無用で尤度をゼロにしてやるぜ、ということであろう]
これはどういうことかというと、いま潜在変数\( z^x_i \)があるとして、その範囲は、下の下閾値\(a^x_i\), 上閾値\(b^x_i\)
によって切断されるということである。$$ a^x_i = max_{j: r^x_j < r^x_i} (z^x_j) $$ $$ b^x_i = min_{j: r^x_j > r^x_i} (z^x_j) $$ [一本目の式はこういうことだ。いま正規分布に従うスコア\(z^x_1, z^x_2, \ldots\)があり、これに順位\(r^x_1, r^x_2, \ldots\)を振って(小さい順に順位を振る)、元のスコアは消してしまったとする。\(r^x_i\)の背後にあったスコア\(z^x_i\)はどんな値か。それはわからないんだけど、すくなくともこういうことはいえる。俺の順位が\(r^x_i\)だとしよう。全員に振られた順位\(r^x_1, \ldots, r^x_{n_1}\)から、俺の順位\(r^x_i\)よか若い奴らを抜き出してくる。そいつらが持っていたスコアのなかで一番大きい奴を\(a^x_i\)とする。それは俺が持っていたスコアより必ず小さい、つまり俺のスコアの下限だ]
このように、潜在変数のある値の上限と下限はその潜在変数の他の値で決まるので、MCMCにおいては強い自己相関が生じる[そうそう、気持ち悪いよね…]。これに対処するため、自己相関を消すステップを追加する。Liu & Sabatti(2000 Biometrika), Morey, Rouder & Speckman(2008 J.Math.Psych.)をみよ。[ふうん…]
[このくだり、なにもこんな風に考えなくても、たとえば2(名義)x3(順序)のクロス表のU検定で順序変数の背後に潜在正規変数を考えるのならば、順序変数の2つの閾値をパラメータにしたほうが早いんじゃないの？ と不思議に思ったんだけど、ここではクロス表の検定だけじゃなくて、\( (r^x, r^y) \)の水準数がすごく多い場合も考慮しているから、きっと閾値をパラメータにしたくないのだろう]
この同時事後分布が手に入れば、\(\theta\)の周辺事後分布も手に入るし、\(H_0: \theta = \theta_0\)の予測性能をなんらかの\(H_1\)の予測性能と比べるBFを求めることもできる。BFは、関心ある仮説のオッズが事前から事後へと変化した程度であると解釈できる。$$ \frac{p(H_1)}{p(H_0)} \times \frac{p(data|H_1)}{p(data|H_0)} = \frac{p(H_1|data)}{p(H_0|data)} $$ 左から順に、事前オッズ、\(BF_{10}\), 事後オッズである。\(BF_{10}=7\)だったら、\(H_1\)の下でのデータのもっともらしさは\(H_0\)の下でのもっともらしさの7倍だということである。
ネストされたモデルであれば、BFはSavege-Dickey密度比を使って簡単に手に入る。$$ BF_{10} = \frac{p(\theta_0|H_0)}{p(\theta_0|data, H_1)}$$ [Savege-Dickey密度比については、Wagenmakers, Lodewyckx, Kuriyal, Grasman (2010 CogPsy)というのがreferされている。CogPsyか…]
3. ケース1. Wilcoxon rank sum 検定
データを\( x = (x_1, \ldots, x_{n_1}), y = (y_1, \ldots, y_{n_2}) \)とする。ふたつを積み上げて順位に変換して切り離す。\(x_i\)の順位を\(r^x_i\), \(y_i\)の順位を\(r^y_i\)とする。
Wilcoxon順位和検定の検定統計量\(U\)は、\(r^x_i\)を合計して\(n_1(n_1+1)\)を引いたもの、ないし\(r^y_i\)を合計して\(n_2(n_2+1)\)を引いたものである[添字が\(n_x, n_y\)になっていたけど\(n_1, n_2\)に直した]。これを差がないという帰無仮説の下での値\(n_1 n_2/2\)と比べる。
\(U\)はサンプルサイズに依存するので、かわりにその標準化された効果量である点双列相関(2値変数と順序変数の連関の指標)を使おう。$$\rho_{rb} = 1 – \frac{2U}{n_1 n_2} $$ これはですね、\(Q(\cdots)\)を符号を返す関数として、$$ \rho_{rb} = \frac{\sum_{i=1}^{n_1} \sum_{j=1}^{n_2} Q(x_i – y_i)}{n_1 n_2} $$と書いても良い。
さて、ふたつの検定手法があるとき、同じレベルの検定力に到達するために必要な観察数の比をAREという。データが正規分布しているとして、2標本 t 検定に対するrank sum検定のAREはおよそ0.955。ちょっとしか下がらないことがわかっている。分布が裾を引くとAREは1より大きくなる。
さあ！これをベイズ化します。
モデルは次のとおり。$$ \delta \sim Normal(0, g)$$ $$ g \sim InverseGamma \left( \frac{1}{2}, \frac{\gamma^2}{2} \right) $$ $$ Z^x_i \sim Normal \left( -\frac{1}{2} \delta, 1 \right)$$ $$Z^y_i \sim Normal \left( \frac{1}{2} \delta, 1 \right)$$ $$ r^x_i \leftarrow Rank(Z^x_i)$$ $$r^y_i \leftarrow Rank(Z^y_i)$$ [面倒くさいので略記したけど、\(Rank(Z^x_i), Rank(Z^y_i)\)はいずれも\((Z^x_1, \ldots, Z^x_n, Z^y_1, \ldots, Z^y_n)\)のなかでの順序である。あれ？ 観察数は\(n_1, n_2\)だったのにいつのまにか両方\(n\)になっておる… 誤記だろうな…]
\(\delta\)の事前分布は尺度\(\gamma\)のコーシー分布にしたいのだが、計算を簡単にするため上のようにする。[このくだりRouder, et al.(2009 PsychonomicBull.Rev.)というのがreferされている。これ読んでおいたほうがよさそうだ…]
で、次のようにGibbsサンプリングする。[以下では\(n_1, n_2\)じゃなくて\(n_x, n_y\)になっている… ちゃんと校正してほしい…]
- 個々の\(x_i\)について切断正規分布をつくって、そこから潜在変数スコアをサンプリングする。$$ (Z^x_i | z^x_i, z^y_i, \delta) \sim N_{(a^x_i, b^x_i)} \left( -\frac{1}{2} \delta, 1 \right) $$
- 同様に、個々の\(y_i\)について、$$ (Z^y_i | z^y_i, z^y_i, \delta) \sim N_{(a^y_i, b^y_i)} \left( \frac{1}{2} \delta, 1 \right) $$
- \(\delta\)のサンプリング。$$(\delta | z^x, z^y, g) \sim N(\mu_\delta, \sigma_\delta)$$ $$ \mu_\delta = \frac{2g(n_y \bar{z}^y – n_x \bar{z}^x)}{g(n_x + n_y)+4}$$ $$\sigma^2_\delta = \frac{4g}{g(n_x+n_y)+4} $$
- \(g\)のサンプリング。$$ (G|\delta) \sim InverseGamma \left(1, \frac{\delta^2+\gamma^2}{2} \right)$$
[途中から写経になっちゃったんだけど、1と2の式がどうにも納得いかない。
まず式の右辺の上下閾値だけど、この書き方だと、\(x\)の\(i\)番目の値の順位\(r^x_i\)の背後にある潜在変数\(Z^x_i\)について両側を切断した正規分布を考えるとき、その閾値はそのイテレーションの時点での\( (z^x_1, \ldots, z^x_{n_1}) \)のみに依存して決まることになるのではなかろうか。でも、たとえば\(r^y_j < r^x_i\)だったら\( z^y_j > r^x_i \)でないとおかしいから、閾値は\( (z^y_1, \ldots, z^y_{n_2}) \)にも依存するのではないかと思う。ほんとは\(r = (r^x_1, \ldots, r^x_{n_1}, r^y_1,\ldots, r^y_{n_2}), z = (z^x_1, \ldots, z^x_{n_1}, z^y_1 \ldots, z^y_{n_2})\)という風に定義しておいて$$ a^x_i = max_{j: r_j < r^x_i} (z_j) $$ $$ b^x_i = max_{j: r_j > r^x_i} (z_j) $$と書かないといかんのではなかろうか？ (追記: 著者らが配っているRコードを確認したけど、やはりそうなっていた)
それから式の左辺、\(Z^x_i\) は新しく得るサンプル、\(z^x_i\)は前のイテレーションで得たサンプルということではないかと想像しているんだけど(きちんと説明されてない)、\(Z^x_i\)が\(z^x_i\)に条件付けられるっていうのはおかしい(実のところ、\(Z^x_i\)はいろんなものに条件付けらているけど\(z^x_i\)だけは関係ないのではないかしらん)。これって、\( ( Z^x_i | z^x, z^y, \delta ) \)の間違いじゃありませんかね???
…もっとも、こういうときって、たいていは著者の誤記じゃなくて俺の理解不足なんだよな… 辛いなあ…]
シミュレーション。
データを生成します。動かす要因は:
- \( \Delta \) [これって\( \delta \)の間違いじゃないの?] : (0, 0,5, 1,5)
- \(n\): (10, 20, 50)
- 分布: skew-normal, Cauchy, logistic, uniform
- 2群が同じ分布に従うか、片方の群は正規分布か
それぞれ1000個データを生成した。それぞれについてUとかBFとかを求めた。
結果。横軸に(Wilcoxon順位和検定の検定統計量の代用として)\( \rho_{rb} \)をとり、縦軸に\(BF_{01}\)をとると、左右対称の山形となり、\(\rho_{rb}=0\)でBFが最大になる。\(n\)が大きいと山が真ん中に集中する。そりゃまあそうなってほしいですよね。
今度はパラメトリックなBFと提案手法のBFを比べてみる。散布図を書くと大体一致するが、提案手法のほうが安定している。また、提案手法は分布形の影響をあまりうけない。
データ例。[実データ分析例の短い紹介。パス]
4. ケース2. Wilcoxon signed rank test [いま関心ないのでパス]
5. ケース3. Spearmanの\(\rho_s\) [いま関心ないのでパス]
6. 結論
このフレームワークでは潜在正規分布を仮定したけど、別の他の分布でも良い。
… 読んでみたはいいけれど… 勉強になったし、ちゃんとRのコードも配られていて助かるんだけど… 嗚呼、助けて神様… たかがクロス表のベイズ・ファクターを求めるために、いちいちサンプリングなんてしたくない… BICでいい感じに近似したい…
今抱えているのは、2(名義) x K(順序)のクロス表があって、Mann-Whitney-Wilcoxon 順位和検定で2群を比較するかわりにBFを求めたい、それもできるかぎり手っ取り早く求めたい、という問題なのである。なにかの一般化線型モデルでBICを求めてBFを近似するとして、その一般化線形モデルってなんだろう。順序反応を目的変数にした比例オッズモデルに、説明変数として名義変数を入れたときのBICと入れないときのBICを求めればよいのだろうか。いや、そうするといつのまにか比例オッズ性という強い仮定が入っちゃうよね… うううう、わからん…