日本経営システム学会

Size: px

Start display at page:

Download "日本経営システム学会"

もりよりわにべ
5 years ago
Views:

統計的検定における多重比較に関する一考察名古屋大学古橋武 A Study on Multiple Comparisons in Statistical Test Nagoya University Takeshi Furuhashi 1 はじめに多重比較は研究者にとって間違いやすく, いささかならずややこしい問題である.

1 統計的検定における多重比較に関する一考察名古屋大学古橋武 A Study on Multiple Comparisons in Statistical Test Nagoya University Takeshi Furuhashi 1 はじめに多重比較は研究者にとって間違いやすく, いささかならずややこしい問題である. 本稿は統計的検定における多重性の問題について, 全て Excel のシミュレーションによる具体例を紹介しながら, 多重比較法の基本的な考え方を解説する. シミュレーションは Excel2007 を用いている. 本稿中の Excel ファイルは全て次の URL よりダウンロードできる. 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)

2 第 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 = である. 図の例ではという値が得られている.1000 組のシミュレーションでは確率的ゆらぎが大きいので, この 1000 組のシミュレーションを 100 回繰り返すマクロを作成して, 第 1 種の過誤の確率の平均値を求め, その 95% 信頼区間を計算した. マクロは当該 Excel ファイルのメニューバーの開発マクロ実行とクリックすることで実行できる. マクロの中身は実行の代わりに編集をクリックすることで見ることができる. 結果を第 26 行目に示す. 第 1 種の過誤の確率の真値が 95% の確率で ~ の範囲内にあることが分かる. 図 2.2 は, 独立な事象 X 1, X 2 の場合の t 分布の結合確率密度関数の鳥瞰図 ( 同図 (a)) と等高線を描いた平面図 ( 同図 (b)) である. ステューデントの t 検定において, 両側検定の場合,p 値 = 0.05 に対応する検定統計量 t = である. 図 2.3 は図 2.2 において, t > の領域を示す. 鳥瞰図では t < の領域の値を 0 としている. t > の領域は小さな山として残っている. 平面図 (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 の領域を黒く塗ってある.

シダックの方法は, 各検定の p 値が名義水準 α' を下回った場合に帰無仮説を棄却する. この方法は, FWER = α を保証する. 図 3.1 はシダックの方法により 2 重の検定を 1000 組実施するシミュレーション画面の抜粋である. 図 2.1 との違いは, 各検定における p 値の閾値をセル G4 にて名義水準 α' = 0.0253 としている点だけである.

3 (a) 鳥瞰図 t 2 t α 図 3.1 シダックの方法による検定を 1000 組実施するシミュレーションの画面 ( シタックの方法 (2 事象 ) のシミュレーション.xlsm) -t α -t α (b) 平面図図 2.3 t 検定を 2 重に繰り返す場合の第 1 種の過誤の領域 (t α = 2.306) 図 2_2 では t > の領域を黒く塗ってある. 第 1 種の過誤の確率は t > の領域の体積で与えられる. この体積がである. このように 1 つの組の中で検定を複数回繰り返すことで, 第 1 種の過誤が有意水準 α を超えてしまう問題が多重性の問題と呼ばれる. t α t 1 α = 1 (1 FWER) 1 2 (3.3) = と求められる.α' は名義水準と呼ばれる. シダックの方法は, 各検定の p 値が名義水準 α' を下回った場合に帰無仮説を棄却する. この方法は, FWER = α を保証する. 図 3.1 はシダックの方法により 2 重の検定を 1000 組実施するシミュレーション画面の抜粋である. 図 2.1 との違いは, 各検定における p 値の閾値をセル G4 にて名義水準 α' = としている点だけである. 第 20 行目の各セルでは,p 値が α' を下回った場合に 1 が出力されている. 第 26 行目より,FWER = 0.05 となる結果が得られてい 3 多重比較法 [1, 2, 3] る. 図 3.2 は図 2.3 と比較すると, 閾値 t a' = と多重比較法は, 有意水準を α とすると 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 回検定を繰り返す場合に拡張した場合と較べて

4 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% 信頼区間は ~ と得られている. 図 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 = (3.10) と与えられる. 図 3.5 はこの閾値により第 1 種の過誤が起きる領域を黒く塗りつぶして示す. 実際には q 23 値の軸を加えた 4 次元空間内の黒色部分の

7) 式の q ij の絶対値が q max0 の値を超えた場合に, その平均値の差に有意差有りとしている.27 行目に求められているように,FWER 0.05 が得られている. q max0 図 3.5 テューキーの方法 (q max0 = 2.

5 (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 組実施するシミュレーションの画面の抜粋である.

6 図 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 = との比較がなされている. 閾値 q max0 はステューデント化された範囲の数表より, データ群数 N' = N - 1= 3, 自由度 ν = N(n - 1) = 32 のとき q max0 =3.475/ 2 = である.FWER 0.05 が達成されている. もし, 上位検定を介さずに,4 群に対してテューキーの方法を適用する場合は, 数表より N = 4, ν = 32 のとき q max0 =3.832/ 2 = であるので, ヘイターフィッシャーの方法により閾値の緩和が実現されている.

図 4.3 ホルムの方法 (5 事象 ) を 1000 組実施するシミュレーション ( ホルム (5 事象 ) によるシミュレーション.xlsm) 進む. (3) t t -1 とする.t = 0 であれば,(4) へ進む.(4.8) 図 4.

3 ホルムの方法 [9] ホルムの方法は下位検定法として位置づけられているわけではないが, 前節までの下位検定における閾値緩和の考え方を発展させた方法として位置づけることができる.

1 節に記したボンフェローニの方法に対して, 以下の手順で閾値 ( 名義水準 )α th を緩和する. (1) 帰無仮説 H i を p 値の小さな順に並べる. pp jj pp kk pp ll 1 jj, kk,, ll N (4.8) であったとする.i = j, t = N とする.

7 図 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 が出力されている. 結果として帰無仮説 H 1, H 2, H 4, H 5 は棄却されず, 検定が終了している. セル B25 では B24~F24 にて 1 回以上 1 があるときに 1 が出力される. セル B27 では 1000 組の検定において第 1 種の過誤が起きた頻度が求められている. 図の例では 43 組で第 1 種の過誤が起きていた.

8 閾値マクロにより, この 1000 組の検定を 100 回繰り返して, 第 1 種の過誤の確率の 95% 信頼区間を求めた結果が第 29 行目に表示されている.FWER 0.05 である. FDR_TST 法新手法? ボンフェローニの方法ホルムの方法検定の回数図 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 法の手順は以下の通りである.

(TST2) EE (VV + SS) = V + S として (TST3) へ進む. (TST3) 次式により αα 2 を決定し,BH 法を適用する. αα 2 = αα 1 NN/(NN EE (VV + SS)) (5.7) 終了したら (TST4) へ進む. (TST4) 検定を終了する., 図 5.1 は TST 法によるシミュレーション画面の抜粋である.

9 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 組について検定を実施するシミュレーション画面の抜粋である.

?? この手法は致命的な間違いを犯している. 理由は, 本稿を前章まで読まれた方にはお分かりいただけると思いますが, いかがでしょうか? 7 結論本稿では Excel を用いたシミュレーションにより多重性の問題を例示し, 基本的な多重比較法の解説を, 閾値の緩和の観点から展開した. もとより, 本稿ではカバーできなかった多重比較法がたくさんある. 別の機会に報告したい.

10 参考文献図 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 (2007) [5] H. I. Braun(Ed.): The Collected Works of John W. Tukey, Vol. III Multiple Comparisons: , 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 (1986) [8] N. Salkind(Ed.): Fisher's LSD, Encyclopedia of Measurement and Statistics, Vol.1, pp (2007) [9] S. Holm: A Simple Sequentially Rejective Multiple Test procedure, Scandinavian J. of Statistics, Vol. 6, No. 2, pp (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 (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 (2006)

ビジネス統計統計基礎とエクセル分析　正誤表

ビジネス統計統計基礎とエクセル分析　正誤表ビジネス統計統計基礎とエクセル分析ビジネス統計スペシャリストエクセル分析スペシャリスト公式テキスト正誤表と学習用データ更新履歴平成 30 年 5 月 14 日現在公式テキスト正誤表頁場所誤正修正 6 知識編第章 -3-3 最頻値の解説内容たとえば, 表.1 のデータであれば, 最頻値は 167.5cm というたとえば, 表.1 のデータであれば, 最頻値は 165.0cm ということになります