統計的検定における多重比較に関する一考察 名古屋大学古橋武 A Study on Multiple Comparisons in Statistical Test Nagoya University Takeshi Furuhashi 1 はじめに多重比較は研究者にとって間違いやすく, いささかならずややこしい問題である. 本稿は統計的検定における多重性の問題について, 全て Excel のシミュレーションによる具体例を紹介しながら, 多重比較法の基本的な考え方を解説する. シミュレーションは Excel2007 を用いている. 本稿中の Excel ファイルは全て次の URL よりダウンロードできる. http://www.echo.nuee.nagoya-u.ac.jp/~furuhasi/educat ion/multiple_comparison/index.html 多重比較法の進展は, 各手法が提案された順序を無視すれば,FWER = αを確保しながら, いかにして閾値を緩めるか, という観点から見ることができる. シダックの方法は事象間が独立であるとして閾値を決定している. テューキーの方法は, 事象間に相関がある検定において, 閾値を緩めている. 下位検定法 ( ヘイター フィッシャー方法 ) は上位検定法 ( 分散分析 ) よりもデータ群数を 1 つ減らして閾値を決定でき, 閾値を緩められる. ホルムの方法では, この考えを推し進めて, データ群数を 1 つずつ減らしながら下位検定を繰り返すことで, 閾値のさらなる緩和を実現している.FDR は FWER とは異なる新しい発想で, 閾値の大幅緩和に成功している. 本稿では, これらの手法をシミュレーション例を付して紹介する. 本稿の最後では, 多重比較法を深く理解するための課題を提示する. かったのに帰無仮説を棄却しない確率は第 2 種の 過誤と呼ばれる. 通常, 統計的検定は第 1 種の過 誤の確率を有意水準以下とする. 図 2.1 は 1 つの組の中でステューデントの t 検 定を 2 回繰り返すことを 1000 組について実施する シミュレーション画面の抜粋である. 図説の xxx.xlsm はこのシミュレーション用の Excel のフ ァイル名である. 多重比較とは 1 つの組の中で検 定を複数回繰り返すことである. 図の例では, 独 立な事象 X 1,X 2 の母分散をそれぞれ µ 1,µ 2 とする と, 帰無仮説 H 1, H 2 を HH 1 : μμ 1 = μμ 0 (2.1) HH 2 : μμ 2 = μμ 0 としている. ここで, 帰無仮説 H 1, H 2 はファミリ ーと呼ばれる. 多重比較における Family-wise Error Rate (FWER) は FWER = PPPP {VV 1} (2.2) と定義される.V は帰無仮説のファミリーの中で 2 多重性の問題 帰無仮説が正しかったのに棄却してしまう過誤 は第 1 種の過誤と呼ばれる. 逆に対立仮説が正し 図 2.1 ステューデントの t 検定を 2 回繰り返すことを 1000 組実施するシミュレーションの画面 (t 検定 (2 事象 ) のシミュレーション.xlsm)
第 1 種の過誤が起きる回数である.FWER はファ ミリーの中で第 1 種の過誤が少なくとも 1 回起き る確率である. 各事象のサンプルサイズ n = 9 である. セル B7 ~B15 に事象 X 1 の第 1 組のサンプル値が計算され ている. 各サンプル値は NORMDIST(RAND(), µ 0, σ 0) 関数により平均 µ 0, 標準偏差 σ 0 の正規乱数が生 成されている.µ 0 と σ 0 はセル B4, C4 にてそれぞれ 1 に指定されている. セル B16 で平均値 xx 1 が, セ ル B17 で不偏分散 vv ee 2 が計算されている. セル B18 では tt 1 = xx 1 μμ 0 vv ee 2 /nn (2.3) により検定統計量が計算されている. この t 値は 自由度 8 (= n - 1) の t 分布に従うことが知られてい る. また, セル B19 では TDIST(t, n-1, 2) 関数によ り両側検定の場合の p 値が求められている. ここ で,3 番目の引数の 2 は両側検定の指定である. セル B20 では p 値がセル F4 の有意水準を下回っ た場合に 1 が出力される. この例では, セル B20 が 1 であるので, 第 1 組の事象 X 1 の帰無仮説 H 1 は棄却されている. サンプルは µ 0 と σ 0 に従った正 規乱数であるので, 正しい帰無仮説が棄却された ので, ここで第 1 種の過誤が起きている. C 列のセル C6~C20 では事象 X 2 の 9 個のサン プルに対して t 検定が実施されている. セル B21 では B20,C20 の少なくともいずれかで 1 となっ た場合に 1 が出力される. これが 1 であれば, 帰 無仮説のファミリーの中で第 1 種の過誤が起きた ことを意味する. この Excel ファイルではセル B5~C21 を横に 1000 組コピーしてある. セル B23 では第 21 行目 の中の 1 の数を数えて, 全体の組数 1000 で割り, 第 1 種の過誤の起きた頻度が求められている. そ れぞれの検定における閾値を α とすると 2 重の検 定の場合の FWER は で,FWER = 0.0975 である. 図の例では 0.092 という値が得られている.1000 組のシミュレーションでは確率的ゆらぎが大きいので, この 1000 組のシミュレーションを 100 回繰り返すマクロを作成して, 第 1 種の過誤の確率の平均値を求め, その 95% 信頼区間を計算した. マクロは当該 Excel ファイルのメニューバーの 開発 マクロ 実行 とクリックすることで実行できる. マクロの中身は 実行 の代わりに 編集 をクリックすることで見ることができる. 結果を第 26 行目に示す. 第 1 種の過誤の確率の真値が 95% の確率で 0.0959~0.1003 の範囲内にあることが分かる. 図 2.2 は, 独立な事象 X 1, X 2 の場合の t 分布の結合確率密度関数の鳥瞰図 ( 同図 (a)) と等高線を描いた平面図 ( 同図 (b)) である. ステューデントの t 検定において, 両側検定の場合,p 値 = 0.05 に対応する検定統計量 t = 2.306 である. 図 2.3 は図 2.2 において, t > 2.306 の領域を示す. 鳥瞰図では t < 2.306 の領域の値を 0 としている. t > 2.306 の領域は小さな山として残っている. 平面図 (a) 鳥瞰図 t 2 t 1 FWER = 1 (1 α) 2 (2.4) となる. セル F4 にて α = 0.05 と設定されているの (b) 平面図 図 2.2 独立な事象 X 1, X 2 の場合の t 分布の結合確率密度関数
(a) 鳥瞰図 t 2 t α 図 3.1 シダックの方法による検定を 1000 組実施するシミュレーションの画面 ( シタ ックの方法 (2 事象 ) のシミュレーション.xlsm) -t α -t α (b) 平面図 図 2.3 t 検定を 2 重に繰り返す場合の第 1 種の過誤の領域 (t α = 2.306) 図 2_2 では t > 2.306 の領域を黒く塗ってある. 第 1 種の過誤の確率は t > 2.306 の領域の体積で与え られる. この体積が 0.0975 である. このように 1 つの組の中で検定を複数回繰り返 すことで, 第 1 種の過誤が有意水準 α を超えてし まう問題が多重性の問題と呼ばれる. t α t 1 α = 1 (1 FWER) 1 2 (3.3) = 0.0253 と求められる.α' は名義水準と呼ばれる. シダッ クの方法は, 各検定の p 値が名義水準 α' を下回 った場合に帰無仮説を棄却する. この方法は, FWER = α を保証する. 図 3.1 はシダックの方法により 2 重の検定を 1000 組実施するシミュレーション画面の抜粋で ある. 図 2.1 との違いは, 各検定における p 値の 閾値をセル G4 にて名義水準 α' = 0.0253 としてい る点だけである. 第 20 行目の各セルでは,p 値が α' を下回った場合に 1 が出力されている. 第 26 行目より,FWER = 0.05 となる結果が得られてい 3 多重比較法 [1, 2, 3] る. 図 3.2 は図 2.3 と比較すると, 閾値 t a' = 2.743 と 多重比較法は, 有意水準を α とすると FWER α (3.1) とする手法である. 3.1 シダックの方法 [4] 同じ組の中で検定を 2 回繰り返す例を取りあげ る. 事象 X 1, X 2 は独立とする. 帰無仮説のファミ リーは HH ii :μμ ii = μμ 0 ii = 1, 2 (3.2) である.(2.4) 式において,FWER = 0.05 とする閾 値を α' とすると して, 黒塗りの部分を狭くしてある. この t > t a' の領域が, 事象 X 1, X 2 の少なくとも一方において 検定統計量の t 値が閾値を超えた場合に, 対応す る帰無仮説を棄却する領域である. この領域の体 積を 0.05 (= FWER) としている. なお, ボンフェローニの方法はシダックの方法 の簡略法である. 同じ組の中で N 回検定を繰り返 す場合は α = FWER/NN (3.4) とする. ボンフェローニの方法は (3.3) 式を N 回検 定を繰り返す場合に拡張した場合と較べて
FWER/NN < 1 (1 FFFFFFFF) 1/NN (3.5) が常に成立する. 計算機が身近になかった時代に は便利な方法であった. 3.2 テューキーの方法 [5] 前節の例では, 検定対象である事象 X 1, X 2 の平 均値は独立であった. 検定の課題には各事象の母 平均の差の検定がある. 帰無仮説のファミリーを HH 12 :μμ 1 = μμ 2 HH 13 :μμ 1 = μμ 3 (3.6) HH 23 :μμ 2 = μμ 3 とする. 事象 X 1, X 2, X 3 が互いに独立であっても, 平均値の差 XX 1 XX 2, XX 1 XX 3, XX 2 XX 3 には従属関係 がある. t α 図 3.2 シダックの方法 (t α = 2.743) t α - t α - t α t 1 図 3.3 平均値の差の検定 (3 事象, 多重比較を考慮せず ) を 1000 組実施するシミュレーション ( 平均値の差の検定 (3 事象 ) シミュレーション.xlxm) t 2 まず, 多重比較を考慮しないシミュレーション を実施する. 図 3.3 は, 互いに独立な事象 X 1, X 2, X 3 の平均値の差について t 検定を 3 回繰り返すこと を,1000 組実施した場合のシミュレーション結果 である. セル B19, C19, D19 でそれぞれ xx 1 xx 2, xx 1 xx 3, xx 2 xx 3 が計算されている. セル B20, C20, D20 の検定統計量 q ij (i, j = 1, 2, 3, ii jj) は次 式による. qq iiii = xx ıı xx ȷȷ / 2vv 2 EE /nn (3.7) ただし,vv EE 2 は各事象の不偏分散の平均値であり, 平均平方と呼ばれる. 検定統計量 q ij は, データ群 数 ( 事象の数 ) を N, 各事象のサンプルサイズを n とすると, 自由度 ν = NN(nn 1) (3.8) の t 分布に従う. 図 3.3 の例では N = 3, n = 9 であ る. セル F4 の閾値 q 0 は,TINV(α, N(n-1)) ( ただし, α = 0.05) によって求められた値である. もし, 平 均値の差 XX 1 XX 2, XX 1 XX 3, XX 2 XX 3 が互いに独立で あれば,(2.4) 式を参考に,FWER=0.143 となると 期待されるが, これらには従属関係があるため, FWER の 95% 信頼区間は 0.1152~0.1198 と得られ ている. 図 3.4 は N = 3, n = 9 のデータを 1 千万組生成し て,q 12 値と q 13 値の分布を求めた結果である.q 12 値 と q 13 値の結合確率密度関数の鳥瞰図と平面図であ る. 図中の縞模様は等高線である.q 12 と q 13 は独立 ではないため結合確率分布は楕円形状となっている. テューキーの方法は, 検定統計量 q ij の最大値 q max の確率分布から, FWER = PP{qq mmmmmm qq mmmmmm0 } = αα (3.9) とする閾値 q max0 を与える.α = 0.05 のとき q max0 はステューデント化された範囲の 5% 点と呼ばれ る.q max0 は解析的には求められないため [1][3], 数 値計算により求められた数表が作成されている. このステューデント化された範囲の 5% 点の数表 より N = 3, ν =24 のとき qq mmmmmm0 = 3.532/ 2 = 2.498 (3.10) と与えられる. 図 3.5 はこの閾値により第 1 種の 過誤が起きる領域を黒く塗りつぶして示す. 実際 には q 23 値の軸を加えた 4 次元空間内の黒色部分の
(a) 鳥瞰図 q 13 (b) 平面図 図 3.4 平均値の差の q 値の結合確率密度関数 q q 12 図 3.6 テューキーの方法 (3 事象 ) を 1000 組実施するシミュレーション ( テューキーの方法 (3 事象 ) のシミュレーション.xlsm) 4 上位検定と下位検定 超体積が 0.05 (=FWER) となる. なお, 数表の値を 1/ 2 倍して用いるのは, 数表が作られた歴史的経緯によるようである. 図 3.6 はテューキーの方法による検定を 1000 組実施するシミュレーション結果である. 図 3.3 との違いは, セル F4 に q max0 の値が与えられた点だけである. 第 21 行目にて (3.7) 式の q ij の絶対値が q max0 の値を超えた場合に, その平均値の差に有意差有りとしている.27 行目に求められているように,FWER 0.05 が得られている. q max0 図 3.5 テューキーの方法 (q max0 = 2.498) q 13 q max0 - q max0 - q max0 q 12 N 重比較の FWER = α とする閾値は,N が大き くなるとともに厳しくなり, 仮説の有意差が出に くくなる. 上位検定 下位検定と呼ばれる手法が あるが, これを閾値緩和の視点で見ていく. 4.1 分散分析 ( 上位検定 ) [6] 一元配置の分散分析について紹介する. 帰無仮 説のファミリーを HH 123 :μμ 1 = μμ 2 = μμ 3 (4.1) とする. このとき対立仮説は 対立仮説 :μμ 1 μμ 2 or μμ 1 μμ 3 or μμ 2 μμ 3 (4.2) である. 帰無仮説 H 123 を棄却できる場合, 対立仮 説のいずれかが有意であることが分かるが, どれ であるかは下位検定を実施しなければ分からない. 分散分析の検定統計量 F は, データ群数が N, 各群のサンプルサイズが等しく n である場合 FF = SS AA/(NN 1) (4.3) SS EE /(nn(nn 1)) により与えられる. この F 値は自由度 N -1, n (N - 1) の F 分布に従う. ここで, NN SS AA = nn(xx ii xx ) 2 ii=1 NN nn SS EE = xx iiii xx ii 2 ii=1 jj=1 (4.4)
図 4.1 分散分析 (3 事象 ) を 1000 組実施するシミュレーション ( 分散分析 (3 事象 ) のシミュレーション.slxm) であり,S A, S E はそれぞれ水準間変動, 水準内変動 と呼ばれる.xx は全体の平均値,xx ii はデータ群 i の平 均値である. 図 4.1 は分散分析による検定を 1000 組実施する シミュレーションの画面の抜粋である.N = 3, 帰 無仮説が (4.1) 式の場合である. セル B18, B19 で (4.4) 式の S A, S E がそれぞれ計算され, セル B20 で (4.3) 式の F 値が求められている. セル B21 では, この F 値から FDIST() 関数により求められる p 値 がセル F4 に入力された有意水準を下回った場合 に 1 が出力される. このとき帰無仮説 H 123 が棄却 される. 第 26 行目の 95% 信頼区間より,FWER 0.05 であることが分かる. 4.2 ヘイター フィッシャーの方法 ( 下位検定 ) [7,8] 分散分析の結果, 帰無仮説が棄却された場合, どの対立仮説が有意であるかを調べるのが下位検 定である. 下位検定法にヘイター フィッシャー の方法 [7] がある. これは, データ群数が N で, 各 群のサンプルサイズが等しい場合, データ群数 N' = N - 1 と 1 つ減らしてテューキーの方法を適用す る方法である. 元のデータ群数 N = 3 の場合にの み適用できる下位検定法にフィッシャーの PLSD 法 [8] があるが, ヘイター フィッシャーの方法は N = 3 のときフィッシャーの PLSD 法と一致する. データ群数を 1 つ減らして N' = N - 1 として良い理由は次のように説明される. 例えば 4 群の場合で各群の母平均が µ 1 >> µ 2, µ 3, µ 4 であったとする. これらの平均値の差により, 上位検定では有意差有りとなる. このとき注意すべきは,µ 2 - µ 3, µ 2 - µ 4, µ 3 - µ 4 については 検定されずに 下位検定に渡されることである. そこで, 下位検定ではデータ群数 N' = 4 ではなく,N' = 3 として,µ 2 - µ 3, µ 2 - µ 4, µ 3 - µ 4 について多重性を考慮した検定を行えば十分である. なお,µ 1 - µ 2, µ 1 - µ 3, µ 1 - µ 4 に関しては検定により有意差有りと出るが, これは正しい結論であり, 第 1 種の過誤ではない. N = 4 とする. 分散分析の帰無仮説は HH 1234 :μμ 1 = μμ 2 = μμ 3 = μμ 4 (4.5) である. 下位検定の帰無仮説のファミリーは HH 12 :μμ 1 = μμ 2 HH 13 :μμ 1 = μμ 3 HH 14 :μμ 1 = μμ 4 HH 23 :μμ 2 = μμ 3 (4.6) HH 24 :μμ 2 = μμ 4 HH 34 :μμ 3 = μμ 4 である. 図 4.2 は 4 群の場合の上位検定を分散分析, 下位検定をヘイター フィッシャーの方法として, 1000 組の検定を実施するシミュレーション画面の抜粋である. 第 22 行目までで図 4.1 と同様の分散分析がなされている. セル B26~E26, B27, C27 にて (3.7) 式の q 12, q 13, q 14, q 23, q 24, q 34 の値 ( サンプルサイズ n = 9, 自由度 ν = N(n - 1)= 32) が求められ, セル B28~E28, B29, C29 にてそれぞれセル E24 に入力されている閾値 q max0 = 2.457 との比較がなされている. 閾値 q max0 はステューデント化された範囲の数表より, データ群数 N' = N - 1= 3, 自由度 ν = N(n - 1) = 32 のとき q max0 =3.475/ 2 = 2.457 である.FWER 0.05 が達成されている. もし, 上位検定を介さずに,4 群に対してテューキーの方法を適用する場合は, 数表より N = 4, ν = 32 のとき q max0 =3.832/ 2 = 2.710 であるので, ヘイター フィッシャーの方法により閾値の緩和が実現されている.
図 4.3 ホルムの方法 (5 事象 ) を 1000 組実施するシミュレーション ( ホルム (5 事象 ) によるシミュレーション.xlsm) 進む. (3) t t -1 とする.t = 0 であれば,(4) へ進む.(4.8) 図 4.2 ヘイターフィッシャーの方法 (4 事象 ) を 1000 組実施するシミュレーション ( ヘイター フィッシャーの方法 (4 事象 ) のシミュレーション.xlsm) 式の順に, 次に大きい p 値の添え字を i に代入 し,(2) へもどる. (4) 検定を終了する. 4.3 ホルムの方法 [9] ホルムの方法は下位検定法として位置づけられ ているわけではないが, 前節までの下位検定にお ける閾値緩和の考え方を発展させた方法として位 置づけることができる. すなわち, 最も検定統計 量の絶対値の大きいものに多重比較法を適用し, 対応する帰無仮説を棄却できれば, データ群を 1 つ削減して閾値を緩和し, 以降同様にして, 多重 比較法の適用, 対応する帰無仮説の棄却, 閾値緩 和を繰り返す方法である. 簡単のため, 帰無仮説のファミリーを HH ii :μμ ii = μμ 0 ii = 1, 2,, NN (4.7) とする.3.1 節に記したボンフェローニの方法に対 して, 以下の手順で閾値 ( 名義水準 )α th を緩和 する. (1) 帰無仮説 H i を p 値の小さな順に並べる. pp jj pp kk pp ll 1 jj, kk,, ll N (4.8) であったとする.i = j, t = N とする. (2) 名義水準 α th = α/t とし,pp ii < αα tth であれば,H i を棄却し,(3) へ進む.pp ii α tth であれば,(4) へ 図 4.3 はデータ群数 N = 5, サンプルサイズ n = 9 の場合にホルムの方法を適用して,1000 組の検定を繰り返すシミュレーション画面の抜粋である. 帰無仮説の数 N = 5 である. セル D4 にて,FWER の有意水準 α = 0.05 が入力されている. 第 23 行目にて第 21 行目の p 値が小さい順に並べ替えられている. セル B24 では, セル B23 の p 値がセル H4 の有意水準 α 5 = α/5 = 0.05/5 = 0.01 を下回ったので, 有意差有りとして 1 が出力されている. 対応する帰無仮説 H 3 が棄却されている. もともと帰無仮説が正しい設定であるので, ここで第 1 種の過誤が起きている. セル C24 ではセル C23 の p 値がセル G4 の有意水準 α 4 = α/4 = 0.0125 より大きかったので,0 が出力されている. 結果として帰無仮説 H 1, H 2, H 4, H 5 は棄却されず, 検定が終了している. セル B25 では B24~F24 にて 1 回以上 1 があるときに 1 が出力される. セル B27 では 1000 組の検定において第 1 種の過誤が起きた頻度が求められている. 図の例では 43 組で第 1 種の過誤が起きていた.
0.06 0.04 0.02 0 閾値 マクロにより, この 1000 組の検定を 100 回繰り 返して, 第 1 種の過誤の確率の 95% 信頼区間を求 めた結果が第 29 行目に表示されている.FWER 0.05 である. FDR_TST 法 新手法? ボンフェローニの方法 ホルムの方法 1 2 3 4 5 検定の回数 図 4.4 ホルムの方法, ボンフェローニの方法, FDR_TST 法, 新手法 (?) の閾値 (N = 5, E(V+S) = 1) 図 4.4 はホルムの方法とボンフェローニの方法 の閾値 α th と検定の回数の関係を示す. ホルムの方 法により, 検定の回数とともに閾値が緩和されて いく様子が分かる. 他の方法の閾値は後述する. 5.1 FDR の定義 FDR は多重比較に対する新しい考え方である. 表 5.1 は帰無仮説の棄却 / 保留と第 1 種 / 第 2 種 の過誤の関係を表す. 帰無仮説の数 N の検定にお いて, 帰無仮説が真である場合に帰無仮説を保留 した数が U, 棄却した数が V, 帰無仮説が偽であ る場合に帰無仮説を保留した数が T, 棄却した数 が S である.V は第 1 種の過誤が起きた数であり, T は第 2 種の過誤が起きた数である.N = U+V+T+S の関係がある. また, 帰無仮説を棄却 した総数が V+S である.FDR は FDR = EE(VV/(VV + SS)) (5.1) と定義される.E(x) は x の期待値である.FDR は 採択した対立仮説の中に含まれる間違った対立仮 説の比率の期待値である.FDR に基づく検定では FDR αα (5.2) となるように, 検定の閾値を緩和し, FWER > α (5.3) となることは許容する.FWER が上式のようであ っても, 採択された対立仮説に含まれる誤った対 立仮説の比率が一定値以下であればそれでよしと する, 新しい考えに基づく多重検定法である. 5 FDR_TST 法 [10, 11] 帰無仮説の数が増えると, 検定の閾値が厳しく なり, 有意差は出にくくなる.FWER を抑えると いう観点では問題ないのであるが, 対立仮説の中 に真のものが含まれていても有意差なしとされる 第 2 種の過誤の確率が高くなってしまう. ゲノム 解析などでは仮説の数に対してサンプルサイズが 小さく, 有意差が出にくいことが問題となってい た. そこで,FWER そのものを見直す考え方とし て FDR(False Discovery Rate) が提唱されている. 表 5.1 帰無仮説を保留帰無仮説を棄却 計 U 第 1 種の過誤 V U + V 第 2 種の過誤 T S T + S 計 U + T V + S N 帰無仮説真 帰無仮説偽 帰無仮説の保留 / 棄却と第 1/2 種の過誤 5.2 BH 法 [10] FDR に基づく検定法の基本的方法は BH 法と呼 ばれる. 簡単のため, 帰無仮説のファミリーを HH ii :μμ ii = μμ 0 ii = 1, 2,, NN (5.4) とする. 以下の手順で検定を行う. なお,α' は次 節の TST 法では αα 1, αα 2 として与えられる. (BH1) 帰無仮説 H i を p 値の大きな順に並べる. pp jj pp kk pp ll 1 jj, kk,, ll NN (5.5) であったとする.i = j, t = N とする. (BH2) 名義水準 α th を αα tth = αα tt/nn (5.6) とし,pp ii < αα tth であれば,H i および H i より p 値 の小さな帰無仮説を全て棄却し,(BH4) へ進む. pp ii α tth であれば,H i を保留して (BH3) へ進む. (BH3) t t - 1 とする.t = 0 であれば,(BH4) へ 進む.(5.5) 式の順に, 次に小さい p 値の添え字 を i に代入し,(BH2) へもどる. (BH4) 検定を終了する.
5.3 TST (Two-stage linear step-up procedure) 法 [11] BH 法を 2 段階 (Two-stage) で適用する方法 (TST 法 ) が提案されている.TST 法は,1 段階目では BH 法を E(V + S) の推定に適用し,2 段階目で FDR α とする検定を行う.TST 法の手順は以下の通り である. (TST1) αα 1 = α /(1 + α), として,BH 法を適用する. 棄却された帰無仮説の数を V + S とする.V + S = 0 であれば,(TST4) へ進む.V + S = N であれ ば, 全ての帰無仮説を棄却して (TST4) へ進む. 0 < V + S < N であれば,(TST2) へ進む. (TST2) EE (VV + SS) = V + S として (TST3) へ進む. (TST3) 次式により αα 2 を決定し,BH 法を適用する. αα 2 = αα 1 NN/(NN EE (VV + SS)) (5.7) 終了したら (TST4) へ進む. (TST4) 検定を終了する., 図 5.1 は TST 法によるシミュレーション画面の 抜粋である. 帰無仮説のファミリーは (5.4) 式で, N = 10 である. ただし,10 のデータ群の内, 事象 X 8, X 9, X 10 はセル D4 の平均値 µ 1 を用いて正規乱数 が生成されている. 表 5.1 において S = 3 となる設 定である. 第 24 行目までが (TST1) の計算である. 第 23 行目で (BH1) の並べ替えがなされ, セル F4 図 5.1 FDR_TST 法 (N = 10, 1000 組 ) のシミュレーション (FDR_TST(10 群 ) によるシミュレーション.xlsm) にて α = 0.05 と与えられているので, セル G4~K6 にて (5.6) 式により α th が計算されている. 第 24 行 目で (BH2), (BH3) の比較がなされ, セル B26 に V + S が得られている.(TST2) にてこの V + S を EE (VV + SS) としている. セル D26 と第 28, 29 行目で (TST3) が実行されている. セル D26 にて (5.7) 式の 計算がなされ, 第 28 行目で (5.6) 式の閾値が求めら れている. 第 29 行目で各帰無仮説の p 値と閾値 α th の比較がなされている. 第 30 行目では, 偽の 帰無仮説の数 S = 3 が分かっているとして,V/(V+S) が計算されている. 図の例では 1/4 = 0.25 である. 第 35 行目の FDR の 95% 信頼区間より, 少し小さ めではあるが FDR αα の値が得られている. 図 4.4 に帰無仮説のファミリーにおいて N = 5, EE (VV + SS) = 1 の場合の閾値と検定の回数の関係を 示してある.FDR_TST 法では, ホルムの方法より もさらに閾値の緩和が実現されている. 6 課題 以下のアルゴリズムを考えてみた. 帰無仮説のフ ァミリーを とする. HH ii :μμ ii = μμ 0 ii = 1, 2,, NN (6.1) (1) 帰無仮説 H i を p 値の大きな順に並べる. pp jj pp kk pp ll 1 jj, kk,, ll NN (6.2) であったとする.i = j, t = N とする. (2) 名義水準 α th を αα tth = αα (6.3) とし,pp ii < αα tth であれば,H i および H i より p 値 の小さな帰無仮説を全て棄却し,(5) へ進む. pp ii α tth であれば,H i を保留して (3) へ進む. (3) t t - 1 とする.t = 1 であれば,(4) へ進む. (6.2) 式の順に, 次に小さい p 値の添え字を i に 代入し,(2) へもどる. (4) 名義水準 α th1 を αα tth1 = NN(1 αα)nn 1 αα + (1 αα) NN (1 αα) NN(1 αα) NN 1 (6.4) とする.pp ll < αα tth1 であれば,H l を棄却し,pp ll α tth1 であれば,H l を保留して (5) へ進む. (5) 検定を終了する.
参考文献 図 6.1 新手法 (?)(N = 5, 1000 組 ) のシミュレーション ( 新手法 (?)(5 群 ) によるシミュレーション.xlsm) ステップ (4) は FWER = α となるように最後の段 階で調整している. 図 6.1 はこの新手法 (?) により,N = 5 の帰無 仮説のファミリーに対して,1000 組について検定 を実施するシミュレーション画面の抜粋である. セル F6 にて α th1 の値が (6.4) 式により計算されてい る. ただし,α = 0.05 である. このシミュレーシ ョンでは帰無仮説は全て正しい設定である. 第 30 行目において FWER 0.05 である. 図 4.4 にこの新手法 (?) の閾値の変化の様子を 書き加えて示す. この手法は,FWER 0.05 を保 ちながら, 閾値の大幅な緩和を実現している??? この手法は致命的な間違いを犯している. 理由 は, 本稿を前章まで読まれた方にはお分かりいた だけると思いますが, いかがでしょうか? 7 結論 本稿では Excel を用いたシミュレーションによ り多重性の問題を例示し, 基本的な多重比較法の 解説を, 閾値の緩和の観点から展開した. もとより, 本稿ではカバーできなかった多重比 較法がたくさんある. 別の機会に報告したい. [1] 広津千尋 : 実験データの解析- 分散分析を超えて- 共立出版(1992) [2] 永田 吉田 : 統計的多重比較法の基礎 サイエンティスト社 (1997) [3] 古橋武編 : 統計 多変量解析とソフトコンピューティング ( 改訂版 ) 共立出版(2014) [4] N. Salkind(Ed.): BonferroniTest, Encyclopedia of Measurement and Statistics, Vol.1, pp.103-107 (2007) [5] H. I. Braun(Ed.): The Collected Works of John W. Tukey, Vol. III Multiple Comparisons: 1948--1983, Chapman & Hall (1994) [6] 石村貞夫 : 分散分析のはなし 東京図書 (1992) [7] A. J. Hayter: The Maximum Familywise Error Rate of Fisher's least Significant Difference Test, J. of American Statistical Association, Vol. 81, No. 396, pp.1000-1004 (1986) [8] N. Salkind(Ed.): Fisher's LSD, Encyclopedia of Measurement and Statistics, Vol.1, pp.359-360 (2007) [9] S. Holm: A Simple Sequentially Rejective Multiple Test procedure, Scandinavian J. of Statistics, Vol. 6, No. 2, pp.65-70 (1979) [10] Y. Benjamini and Y. Hochberg: Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing, J. of the Royal Statistical Society. Series B(Methodological), Vol. 57, No. 1, pp.289-300 (1995) [11] Y. Benjamini, A. M. Krieger and D. Yekutieli: Adaptive Linear Step-up Procedures that Control the False Discovery Rate, Biometrika, Vol. 93, No. 3, pp. 491-507 (2006)