i Armitage Q. Bonferroni 1 SAS ver9.1.3 version up 2 *1 *2 FWE *3 2.1 vs vs vs 2.2 5µg 10µg 20µg 5µg 10µg 20µg vs 5µg vs 10µg vs 20µg *1 *2 *3 FWE 1

多重比較について

（その２−１：具体的手法（

i

）

）

第３回

Armitage

勉強会：資料２

土居正明

Q.多重比較の必要性については分かりました。では、具体的手法としてはBonferroni以外にはどのような方法がありま すか。

1

はじめに

本稿におけるSASのプログラムはver9.1.3の段階のものです。version up後に動く保証はありませんので、予めご了承 ください。 また、検定は基本的に両側検定のみを扱います。ご了承ください。

2

多重比較の必要なとき

多重比較が必要なシチェーションは様々ありますが、今回は３群以上を１つの評価項目に対して比較することを前提とし てご説明します*1_。 多重比較をする際に、最も大事なのは「どの群とどの群を比較するべきなのか」ということです*2_{。「多重比較について} （その１）」から、明らかに 検定の数が増えれば増えるほど、FWEは大きくなる という性質がありますので、できるだけ検 定の回数（比較する組合せ）を少なくしたい、という風に考えます*3_。 そのような視点を持ちつつ、まずは対比較・対照との比較・対比についてご説明しましょう。具体的な手法はあとまわし にして、多重比較が必要とされる典型的なシチェーションを整理します。

2.1

対比較

対比較とは、簡単に言いますと「全群総当りの検定」です。 例えば、「プラセボ・実対照薬・新薬」の３種類を比較しようとする際に、「プラセボ vs 新薬」「実対照薬vs 新薬」「プ ラセボ vs実対照薬」の全てのペアを比較したい、というようなシチェーションを考えます。このようなシチェーションを 対比較 と言います。

2.2

対照との比較

次に、対照との比較 です。 これは例えば、「プラセボ・5µg・10µg・20µg」の４群を比較しようとするときで、「5µg群・10µg群・20µg群のどれ か１群でもプラセボに勝っていたらそれでよしとする」という場合です。「プラセボ vs 5µg」「プラセボvs 10µg」「プラ セボ vs 20µg」の比較のみ行います。 *1_{他の例としては、２群間で２つ以上の評価項目を同時に比較する「多重エンドポイント」もありますし、有効中止を狙う中間解析も「１つの評価項} 目に対する２群比較を複数時点で行う」多重比較となります。 *2_{逆の言い方をすれば「どの群とどの群は比較する必要がないか」ということでもあります。} *3_{逆から見ますと、FWE が適切に制御されているとき、検定の回数が減れば減るほど検出力が上がるのです。}

2.3

対比

最後は、対比 です。 これにはいくつかの種類があります。例えば、薬A・B・C・D・Eの５種類があったときに、 (i)「A・B・C・Dの各群の効き目の平均値よりも、Eの効き目が大きい」 (ii)「A・B・Cの各群の効き目の平均値よりも、D・Eの各群の効き目の平均値が大きい」 (iii)「Aの効き目よりも、B・C・D・Eの各群の効き目の平均値が大きい」 のように「いくつかの投与群をグループ化して比較したい」という場合や、 (iv)「A・B・C・D」において、効き目が「A < B = C = D」や「A < B < C < D」のような傾向性を知りたい というような場合に用いられます。

3

準備

3.1

最初に注意

3.1.1 棄却限界について 一つ最初に大事な注意をしておきます。多重比較と言っても要は検定をしたいので、何か統計量を作ります。その統計量 の従う分布が分かれば、それによって棄却限界も分かる、ということになります。従って我々のするべきことは、しかるべ き統計量を作って、その統計量の従う分布においてFWEをα以下に抑える棄却限界を手に入れる ことです。 そこで確率密度関数を考えて、上側何％点かを求めることになるのですが、厳密に数式で確率密度関数を導出する、とい うことにそれほど興味のない方は、シミュレーションをもとにして棄却限界を作っている、と考えてもよいかと思います*4_。 つまり、 (i)まず２群比較の検定統計量を作る。

(ii)次に、(i)を参考にしてFWEを調整できる統計量を作る。

(iii)乱数を使って帰無仮説のもと*5_{での(ii)の統計量をたくさん（例えば1,000,000個）作ってヒストグラムを作る。} (iv) (iii)のヒストグラムの上から(100· α)％目の値を上側(100· α)％点として、その値を棄却限界とする。 の順番です。 3.1.2 棄却限界の作り方 では次に、棄却限界の作り方について、大体のイメージをご説明します。たとえば、ある２つの検定統計量Z1, Z2が、帰 無仮説のもとで独立にN (0, 1)に従うとしましょう。このとき、この両方を有意水準0.05で検定すると |Z1| > 1.96, |Z2| > 1.96 で棄却すればよくなります。しかし、（その１）で見てきましたように、このように「それぞれ5％ずつ」ではFWEは0.05 を超えてしまいます。計算すると F W E = P (|Z1| > 1.96 or |Z2| > 1.96)≒0.10 となります。そこで、棄却限界1.96をもう少し大きくしてやります。たとえば棄却限界を2.24とすれば、 F W E = P (|Z1| > 2.24 or |Z2| > 2.24)≒0.049 *4_{厳密には、棄却限界は理論的な導出が行われている場合もたくさんありますので、こういう表現は多重比較の専門家の先生には叱られるかもしれま} せん。ですがシミュレーションは式が簡単に表現できない場合でも棄却限界が作れますので、大変便利です。ですので、こういう発想ができるよう になることにはメリットがあると思います。 *5_{ここでいう「帰無仮説」とは、「全ての帰無仮説が全てなりたつ」ことを指します。そうする意味はもう少しあとにご説明します。}

となり、FWE≤ 0.05となります*6_{。このように、「通常の棄却限界だとFWEが大きくなるので、棄却限界を大きくする} ことでFWEの増大を抑える」という方法をとることにします*7_。 さて、その際に F W E = P (|Z1| > 2.24 or |Z2| > 2.24) を考えるのですが、「どちらかが2.24を超える確率」というのは「大きい方が2.24を超える確率」という風に言い換える こともできます。すなわち F W E = P µ max i=1,2|Zi| > 2.24 ¶ です。このようにmax i=1,2|Zi|、すなわち 複数の統計量の最大値 に注目して棄却限界を調べていきます。これが、これからご 説明しますTukey-Kramerの方法やDunnettの方法の棄却限界の作り方の基本です。 3.1.3 多重比較の回数について 仮に「プラセボ・5µg・10µg・20µg」の４群の比較をしたいと思ったとします。そこで「対比較（総当たり：比較６回） と対照との比較（プラセボとの比較：比較３回）のどちらを使おう？」と考えたとしましょう。一つの案としては、「対比較 は総当たりなので、対照との比較も含んでいるから、対比較をやっていれば問題ない」という風に考えるかたもいらっしゃ るかも知れません。ところがこれは結構問題なのです。 そもそもデータ解析の一番の基本は データは、より正しい手法で解析すればより信頼できる答えが出てくる、というもの です。ですので、「対照との比較が知りたい」（比較３回でよい）ときに、わざわざ対比較で６回も比較してしまうと、無駄 が大きくなってしまうのです。 具体的に、前回みましたBonferroniの方法で見てみましょう。これは全体の有意水準をαで抑えたいときに、各検定の 有意水準を「α÷（検定の回数）」で求める方法でした。今、全体の有意水準を0.05としますと、 （対比較） 0.05 6 = 0.0083, （対照との比較) 0.05 3 = 0.0167 と、１回あたりの有意水準が倍違ってしまいます。これは要は「検定を６回すると、３回のときに比べて有意差が出にくく なっている」ことになります。つまり、対照との比較で済むところを対比較してしまったら、検出力が大きく低下してしま うのです。繰り返しますが、データは、より正しい手法で解析することが重要 です。計画段階できちんと「どの比較が必要 でどの比較は必要ないのか」について吟味することが必要です。 3.1.4 データについての仮定 では、対比較・対照との比較・対比のそれぞれについて、具体的手法を見ていくことにしましょう。今回は、(i) データ はそれぞれ独立に正規分布していることと、(ii) 各群の分散が等しいことの ２つを仮定 します。つまり、データは 第１群 x11 x12 · · · x1n1 ∼N (µ1, σ 2₎ 第２群 x21 x22 · · · x2n2 ∼N (µ2, σ 2₎ .. . ... ... · · · ... 第I群 xI1 xI2 · · · xInI ∼N (µI, σ 2₎ *6_{少し面倒ですが計算していきますと、大体 P (|Z1| > 2.29) = 0.025 ですので、} F W E = P (|Z1| > 2.24 or |Z2| > 2.24) = 1− P (|Z1| ≤ 2.24 and |Z2| ≤ 2.24) = 1− P (|Z1| ≤ 2.24) · P (|Z2| ≤ 2.24) (* Z1, Z2は独立) = 1− 0.975 · 0.975 ≒ 0.049 となります。 *7_{棄却限界が大きくなっていますから、通常より「差がつきにくく」なり検出力が下がる、という流れになります。}

で、全て独立と仮定します*8_{。そして、各群の標本平均値をx¯1, ¯x2,}· · · , ¯_{xI、標本不偏分散をそれぞれ} c σ2 1= 1 n1− 1 n1 X j=1 (x1j− ¯x1)2 c σ2 2= 1 n2− 1 n2 X j=1 (x2j− ¯x2)2 .. . ... c σ2 I = 1 nI− 1 nI X j=1 (xIj− ¯xI)2 とおくこととします。 では、個別の手法を見ていくことにしましょう。

4

全ての対比較：

Tukey-Kramer

の方法

最初に全ての対比較をする方法を考えましょう。

4.1

仮説

まずは仮説です。上の通り薬剤群がI群あり、全ての比較を考えますので帰無仮説と対立仮説は*9 H12: µ1= µ2 H12A : µ16= µ2 H13: µ1= µ3 H13A : µ16= µ3 .. . ... H1I : µ1= µI H1IA : µ16= µI H23: µ2= µ3 H23A : µ26= µ3 .. . ... H(I_−1)I : µI−1= µI H_(IA_−1)I : µI−16= µI となります。なお、帰無仮説の数（=対立仮説の数）はそれぞれ全体I群の中から2群を取り出す取り出し方の総数に等し いので、 IC2= I(I− 1) 2 個になります。

4.2

２群比較の検定統計量

では次に、２群比較の検定統計量です。２群ごとの比較を考えればよいので、検定統計量としては通常のt統計量を考え ればよいような気がします。しかし、せっかく全群で分散が同じということを仮定していますので、それを利用してやりま しょう。つまり、「全群の分散をプールする」ということを考えるのです*10_。 通常の両群の分散が等しいことを仮定した両側t検定では、第i群と第j群の比較をする際の検定統計量は、 tij= | ¯ xi− ¯xj| r³ 1 ni + 1 nj ´ c σ2 ij *8_{各群の例数は異なっていてもよいものとします。} *9_{添え字が多くて面倒なので、Hijを帰無仮説、H}A ijを対立仮説とします。 *10_{この点は、分散分析の F 検定でも同じように考えます。詳しくは「正規分布・t 分布・χ}2_{分布・F 分布とは何か？」の資料をご覧ください。}

とおきます*11_{。通常のt検定では、σ}2_{の推定値σ}c2 ij を考える際にi群とj群のσ2の推定値σc2i とσc2j を併合してやって c σ2 ij= 1 (ni− 1) + (nj− 1) ³ (ni− 1)cσ2i + (nj− 1)cσj2 ´ = 1 ni+ nj− 2 Ã_n i X k=1 (xik− ¯xi)2+ nj X k=1 (xjk− ¯xj)2 ! を考えます。しかし、このσc2 ij の役割は「本当は真の分散σ 2_{を使いたいのだけれど、分散は試験をしてみるまで分からな} いので、仕方なく推定値で代用している」というものです。そのため、「σ2の推定はできるだけ精度よくやりたい」と思う のは自然でしょう。今回、せっかく 全群通してデータの分散σ2_{が等しい ことを仮定していますので、σ}2_{の推定にはI群} のデータを全て使ってやるのです。つまり、 c σ2₌ 1 (n1− 1) + (n2− 1) + · · · + (nI− 1) ³ (n1− 1)cσ12+ (n2− 1)cσ22+· · · + (nI− 1)cσI2 ´ = 1 n− I  Xn1 j=1 (x1j− ¯x1)2+ n2 X j=1 (x2j− ¯x2)2+· · · + nI X j=1 (xIj− ¯xI)2   (n = n1+ n2+· · · + nIとおく) = 1 n− I I X i=1 ni X j=1 (xij− ¯xi)2 を用いて、統計量を Tij = | ¯ xi− ¯xj| r³ 1 ni + 1 nj ´ c σ2 とした検定をTukey-Kramerの方法 と言います。そのなかで特に各群の例数が等しいときをTukeyの方法 と呼びます*12_。 ここで、「正規分布・t分布・χ2_{分布・F分布とは何か？」の資料のt分布の導出のところで少し触れましたが、T} ijから 分子の絶対値を除いたものは自由度(n− I)のt分布に従います*13_。

4.3

FWE

を制御する統計量

FWEは「試験全体としてのαエラー」つまり「正しい帰無仮説が、１つでも棄却される確率」でした。これを制御する ためには、「3.1.2 棄却限界の作り方」で見ました通り、I(I₂−1)個の統計量{Tij}のうち、最も棄却されやすいもの、つま り最大のものを制御すればよい、ということになります。つまり T = max i<j Tij = max i<j | ¯xi− ¯xj| r³ 1 ni + 1 nj ´ c σ2 の値を制御することでFWEが制御できます*14_{。ここで、たとえばT12とT13は、どちらも第１群のデータを用いていま} すので、独立ではありません。この点については、補足で少しご説明します。 *11_{通常は分子に絶対値はつけませんが、両側検定を考える場合本質的には同じことなので、今回はつけることとします。なお、絶対値がついている場} 合、有意水準 0.05 の両側検定の棄却域は t > t(n1+ n2− 2, 0.975) で構成されます (t(n1+ n2− 2, 0.975) は自由度 (n1+ n2− 2) の t 分布の 下側 97.5 ％点)。 *12_{歴史的にはおそらく話が逆で、先に例数が等しい Tukey の方法ができて、それを拡張したのが Tukey-Kramer の方法のはずです。} *13_σ2_{の推定に全データを用いることで推定精度がよくなり、t 分布の自由度が増加して裾が正規分布の方に降りて行くため、より小さい値で有意差} が出るようになります。つまり、検出力が増加するのです。 *14_max

i<j の「i < j」は、「i, j 全部」としても構いません。絶対値がついているので| ¯xi− ¯xj| = | ¯xj− ¯xi| となるため片方だけとればよい、ということ で半分に減らしただけです。また、i = j のときは分子が０になり、最大になることはまずありえませんので不要です。

また、大変重要なこととして、対比較・対照との比較の場合FWEを制御する棄却限界は「全ての帰無仮説が正しい」状 況のもとでのT の分布から求めます*15_{。詳しくは「6帰無仮説とFWEに関する注意」でご説明しますが、大体のイメー} ジとしては「FWEは『正しい帰無仮説が１つでも棄却される確率』なので、正しい帰無仮説が多い方がFWEは大きい。 そのため、最も正しい帰無仮説が多い状況（＝全ての帰無仮説が正しい状況）でFWEを制御しておけば、真の状況が他の 場合でも大丈夫」という風に思っていただければ結構です。

4.4

棄却限界

では、棄却限界はどのように考えたらよいでしょうか？ 4.4.1 （方法１）シミュレーションを用いた構成法 例えばI = 3で、n1= 10, n2= 12, n3= 11のとき、帰無仮説・対立仮説は H12: µ1= µ2 H12A : µ16= µ2 H13: µ1= µ3 H13A : µ16= µ3 H23: µ2= µ3 H23A : µ26= µ3 となります。棄却域は帰無仮説のもとで構成されますので、「すべての帰無仮説が正しい」という状況のもとで上のT がど のような分布をするか、ということをシミュレーションにより確認すればよいのです*16_{。具体的な方法は以下の通りです。}   仮にいまµ1= µ2= µ3= 0, σ2= 1だとして、T の分布を考えますと、シミュレーションとして (i)N (0, 1)から独立なデータを(10 + 12 + 11) = 33個取り出します。 (ii)最初の10個からx¯1を、次の12個からx¯2を、最後の11個からx¯3を作ります。 (iii)全データを用いてcσ2_{を作ります。}

(iv) (i)∼(iii)よりT12, T13, T23を作り、最大値をT とします。

(v) (i)∼(iv)を例えば1,000,000回*17_{繰り返して、T のヒストグラムを作ります。} (vi) Tのヒストグラムの上側5％点を選んで、そこを両側5％棄却限界とします*18_。 とすればよくなります。  このように考えると、シミュレーションをもとにして棄却限界を作ることができます。ところが、今の例では µ1 = µ2 = µ3 = 0, σ2 = 1という特殊な状況を考えたのでした。一般の場合*19はどうなっているのでしょうか？ 実は 上と全く同じように、N (0, 1)の元で乱数を発生させてよい のです。このことの詳細なご説明は式の計算が結構大変なので 補足に回しましょう。 4.4.2 (方法２）SASを用いた計算法 では次にSASでどのように実行するかを見ていきましょう。さらに２つの方法があります。(i)ProbMC関数 を用いて 棄却限界を求める方法、(ii)proc glmを使う方法、の2種類です*20_。 (i)ProbMC関数を使う方法

棄却限界(ca)を求めるのに必要なパラメータは「FWE (alpha)」「自由度(df)」「群の数(k)」です。ここで、「自由度」 は（全群の例数）−（群の数）で求まります。データステップで *15_{対比較・対照との比較以外では一般にこうとは限りませんので、状況に応じて考えてください。} *16_{この「すべての帰無仮説が正しい」状況で考えることの妥当性は「6 帰無仮説と FWE に対する注意」で確認します。} *17_{回数は適当ですが 1,000,000 回やれば十分でしょう。} *18_{絶対値をとっていますので、「＋方向に大きい」「−方向に大きい」の両方が「絶対値が大きい」となります。従って、上側だけの棄却域を考えるこ} とで両側検定となります。 *19_{とはいえ全ての帰無仮説が成り立っている場合を考えますので、µ1= µ2= µ3は仮定しています。} *20_{proc mixed でもできますが、今回は省略します。}

data d1;  c = ProbMC(”Range”, . , 1-alpha, df, k);  ca = c /sqrt(2); run; とするとcaが棄却限界となります。ここでそれぞれ ”Range” : 「Tukey-Kramerの方法」を意味します。 . : 2番目にあると「下側％点が欲しい」を意味します。棄却限界が欲しい場合はこうします。 1-alpha : 「Student化した範囲の分布（補足参照）の下側確率」です。分かりにくい方は1-(有意水準）と覚えてください。 df : 「自由度」です。先にも書きましたが「(全例数)−(群の数)」で計算できます。 k : 「群の数」です。

という意味があります。cは「Student化した範囲の分布における下側100(1-alpha)％点」です。Student化した範囲 とは何か、なぜcaが直接求まらずcを求めてから√2で割るのか、については補足でご説明します。ただ、ProbMCで 求めた値cを√2で割らないと棄却限界が求まらないという事実だけは覚えておいてください。 (ii)proc glmを使う方法  たとえばデータセットをd2、投与群を示す変数をdose、被説明変数をbpとしますと、 proc glm data=d2;  class dose;  model bp = dose;

 means dose /tukey; run; で「どの群とどの群で有意差があるか」などが求まります*21_。

4.5

検定の方法

棄却限界を求めると、検定ができます。帰無仮説Hij : µi= µjの有意水準alphaの両側検定は Tij > ca のときに棄却されます(caは上のProbMCで計算したものを√2で割ったもの、もしくはシミュレーションにより得られ た棄却限界)。このようにして、全てのi, jの組み合わせについて、Hijの検定を行えばよいのです。

4.6

balanced

と

unbalanced

balancedとunbalancedの場合について整理しておきましょう。balancedのときをTukeyの方法、unbalancedのとき

をTukey-Kramerの方法と言います。どちらでも同じ方法をとりますが、proc glmやprobMC関数で棄却限界を求める場

合、Tukeyの方法は厳密な方法、Tukey-Kramerの方法は近似的な方法です*22_。

*21_{本当は SAS アウトプットも載せたいのですが、SAS アウトプットを L}A_{TEX に読み込む方法を知りませんのでまだできません。ご存知の方、教え}

て下さい。

5

対照との比較：両側

Dunnett

の方法

では次に対照との比較です。これについては、実は手法としてはTukey-Kramerの方法とほとんど変わりません。最 大の違いは、比較する仮説の数がTukey-Kramerに比べて少なくなりますので、その分差が検出しやすくなる ということ です。あとは、片側検定に意味が出てくる ことも違います。ただ基本的に臨床試験などでは両側検定がメインですので、今 回は両側検定のみ扱うこととします。

5.1

仮説

では仮説です。今度は例えば「第１群との比較」に注目するとすると、「第１群vsその他の群」という比較を考えますの で、全部で(I− 1)回の比較で済みます。したがって仮説は H12: µ1= µ2 H12A : µ16= µ2 H13: µ1= µ3 H13A : µ16= µ3 .. . ... H1I : µ1= µI H1IA : µ16= µI となります。

5.2

２群比較の検定統計量

２群比較の検定統計量は、Tukey-Kramerの方法とほぼ全く同じです。t統計量の分散σ2_{の推定値にはTukey-Kramer} と同様、全データを用いたcσ2_{を利用した} T1j= | ¯ x1− ¯xj| r³ 1 n1 + 1 nj ´ c σ2 とします。

5.3

FWE

を制御できる統計量

Tukey-Kramerの方法と同じく、これもT1j全体の最大値となります。つまり、 T = max 2≦j≦I T1j = max 2≦j≦I | ¯x1− ¯xj| r³ 1 n1 + 1 nj ´ c σ2 を用います。

5.4

棄却限界

棄却限界についてもほとんど同じです。 5.4.1 (方法１)シミュレーションを用いた方法 シミュレーションの手順も Tukey-Kramerと同じように考えることができます。例えばI = 3 で、n1 = 10, n2 = 12, n3= 11のとき、帰無仮説と対立仮説は H12: µ1= µ2 H12A : µ16= µ2 H13: µ1= µ3 H13A : µ16= µ3

となります。

 Tukey-Kramerのときと同じくµ1= µ2= µ3= 0, σ2= 1だとしてよいので、 (i) N (0, 1)から独立なデータを(10 + 12 + 11) = 33個取り出します。

(ii) 最初の10個からx¯1を、次の12個からx¯2を、最後の11個からx¯3を作ります。

(iii) 全データを用いてσc2_{を作ります。}

(iv) (i)∼(iii)よりT12, T13,を作り、最大値（大きい方）をT とします。

(v) (i)∼(iv)を例えば1,000,000回繰り返して、Tのヒストグラムを作ります。

(vi) Tのヒストグラムの上側5％点を選んで、そこを両側5％棄却限界とします。の順に行います。

5.4.2 (方法２）SASを用いた計算法

Tukey-Kramerと同じく(i)ProbMC関数 を使う方法、(ii)proc glmを使う方法、の２つがあります。

(i)ProbMC関数を使う方法 これについては、Tukey-Kramerとよく似ています。違うのは４点で、 (a)balancedとunbalancedで指定が違うこと (b)最初の”Range”の部分を”Dunnett2”とすること*23 (c)√2で割らなくてもよいこと (d)全体k群の中で「１群が対照」なので、「群の数」の代わりに「(群の数)-1」（対照の群を除いた数）とすること*24 だけです。データセットで棄却限界を求めるには (balanced) balancedの場合、 data d1;  ca = ProbMC(”Dunnett2”, . , 1-alpha, df, k-1); run; だけで十分です。Tukey-Kramerの方法と同様に ”Dunnett2” : 「両側Dunnettの方法」を意味します。 . : 2番目にあると「下側％点が欲しい」を意味します。棄却限界が欲しい場合はこうします。 1-alpha : 「下側確率」です。分かりにくい方は1-(有意水準）と覚えてください。 df : 「自由度」です。Tukey-Kramerと同じく「(全例数)-(対照群を含めた群の数)」です。 k-1 : 「(対照群を含めた群の数)-1」です。 なお、Tukey-Kramerにおける「Student化した範囲の分布」に対応するものについては、補足で少しだけ触れます。 *23_{”Dunnett2”の「2」は両側検定を意味します。上で述べましたとおり Dunnett の検定には片側検定がありますので、そちらを使いたい場合} は”Dunnett1”としてください。 *24_{この数は「対照との比較」においては「比較の回数」にも一致します。}

(unbalanced)  では次にunbalancedな方です。３群の比較で１群目（対照群）の例数が１０人、２群目が１３人、３群目が１４人とし て、具体的数値を代入しますと data d1;  n1=10; n2=13; n3=14;  k=3;  df = n1 + n2 + n3 - k;  lambda1 = sqrt(n2/(n1+n2)); lambda2 = sqrt(n3/(n1+n3));

 ca = ProbMC(”Dunnett2”, . , 1-alpha, df, k-1, lambda1, lambda2); run; で求まります。 (ii)proc glmを使う方法 Tukey-Kramerと同様に、データセット：d2、投与群：dose、被説明変数：bp、とおくと (balanced) proc glm data=d2;  class dose;  model bp = dose;

 means dose /dunnett; run;

(unbalanced) proc glm data=d2;

 class dose;

 model bp = dose;

 lsmeans dose / pdiff cl adjust = dunnett; run; で同様に棄却限界・有意差のある組み合わせ、などが求まります。

5.5

検定の方法

Tukey-Kramerと同じく、帰無仮説H1j : µ1= µjの有意水準alphaの両側検定は T1j> ca のときに棄却すればよくなります。

5.6

balanced

と

unbalanced

今回は、balancedとunbalancedの場合の両方をDunnettの方法と呼びます。しかし、Tukey-Kramerとは違い、proc

6

帰無仮説と

FWE

に対する注意

対比・対照との比較のどちらも棄却限界を「すべての帰無仮説が正しい」という状況のもとでのTの分布を考えました。 これによって、きちんとFWEが制御されている ことをみていきます。そのためには、「真の状況がどのようであっても FWEがα以下に抑えられている」ことを確認しないといけません。 例えば、４群の対照との比較で、 H12: µ1= µ2 H13: µ1= µ3 H14: µ1= µ4 の場合を考えましょう。全ての帰無仮説が正しい場合 は F W Eall= P ¡ H₁₂A ∪ H₁₃A ∪ H₁₄A | µ1= µ2= µ3= µ4 ¢ となります。次に、H12, H13が正しく、H14は間違っている場合 は、H14はFWEつまり「正しい帰無仮説が誤って棄却 される確率」には寄与しませんので、 F W E12,13 = P (H12A ∪ H A 13| µ1= µ2= µ3= µ4) となります。これより、明らかに F W E12,13 ≤ F W Eall であることはお分かりいただけると思います。同様に他の場合も考えていきますと、「すべての帰無仮説が正しい」場合 のF W Eallを制御しておけば、それ以外のどの状況に対してもFWEが制御されているということが分かります*25。対比 較についても同様に考えていただけば同じことが成り立ちます。  なお、今確認したのはあくまで 対比較・対照との比較の場合だけ であり、それ以外の場合に「すべての帰無仮説が正しい もとでTの分布を考える」ので十分かどうかは分かりません。それぞれの場合ごとに検討することを忘れないでください。 *25_{厳密には、「真の状況」に正しい対立仮説が多い場合、FWE は α よりだいぶ小さくなってしまいますが、先にも書きました通り「真の状況」は我々} には分かりませんので、保守的にしておくより仕様がありません。なお、今まで FWE を「α以下」に制御する、という表現をしてきたのは、この ような正しい対立仮説がある場合を考慮してのことです。

7

補足１：

FWE

を制御する統計量に関する注意

例えばTukey-Kramerの方法において、FWEを制御する統計量は T = max i<j Tij でした。ここでたとえば T12= | ¯ x1− ¯x2| r³ 1 n1 + 1 n2 ´ c σ2 , T13= | ¯ x1− ¯x3| r³ 1 n1 + 1 n3 ´ c σ2 の２つは共に第１群のデータを用いているため、独立ではありません。 もし仮にこの２つが独立なら、たとえばある定数cに対して P (max{T12, T13} > c) = P (T12> c or T13> c) = 1− P (T12≤ c and T13≤ c) = 1− P (T12≤ c) · P (T13≤ c) となり、FWEを制御する棄却限界を計算するために、「T12 の分布」「T13の分布」のみを考慮すればよくなり、さらに T12, T13はt検定の統計量とほぼ同じであることから、計算はおそらくそれほど難しくありません*26。しかし、実際は独立 でなく、統計量同士の相関が生じてしまいますので、計算は相当大変なはずです。 なお、シミュレーションを用いて棄却限界を計算するときは、そのようなことは一切考慮しなくてもよいですので大変楽 ですが、その分その後ろにある数理的構造などがさっぱり分かりません。どちらも一長一短だと思います。

8

補足２：

Tukey-Kramer

の方法のシミュレーションによる棄却限界について

Tukey-Kramerの方法のシミュレーションによる棄却限界を考えるときに、現実のµ1= µ2 = µ3, σ2がどのような値で あっても、xij をN (0, 1)から 独立にとって考えればよい、という風に述べました。このことをご説明しましょう。実は、 単純な２群比較のt統計量を帰無仮説のもとで考えれば十分なので、その場合にご説明します。 言いたいことは、「N (µ, σ2_{)からのサンプリングによる自由度nのt統計量の分布の上側5％棄却限界が、N (0, 1)から} のサンプリングによる自由度nのt統計量の分布の上側5％棄却限界と等しいこと」です。 N (µ, σ2)のもとでのt統計量とN (0, 1)のもとでのt統計量では、大きさが等しい、ということが言いたいのです。 データとしては、独立な x11,· · · , x1n1 ∼N (µ, σ 2₎ x21,· · · , x2n2 ∼N (µ, σ 2₎ が欲しいのですが、これらのデータを、「N (0, 1)からのデータを加工したもの」とみることにしましょう。つまり、直接 x11,· · · , x1n1, x21,· · · , x2n2を発せさせるのではなく、まず最初に y11,· · · , y1n1 ∼N (0, 1) y21,· · · , y2n2 ∼N (0, 1) を独立に発生させます。そして、 xij = √ σ2_y ij+ µ (i = 1, 2, j = 1,· · · , ni) *26_{著者はきちんと計算を追いかけたことがないので正確には分かりませんが。もしかしたら、群ごとに例数が異なっているときは意外と大変かもしれ} ません。

として{xij}を作ります。このときx11,· · · , x2n2 が全て独立なこととy11,· · · , y2n2が独立なことは同値になります。こう してy11,· · · , y2n2をもとにして、独立なx11,· · · , x2n2が発生したことになります。 では、x11,· · · , x2n2 におけるt統計量を{yij}の言葉で書いてみましょう*27。 ¯ x1= 1 n1 n1 X j=1 x1j = 1 n1 n1 X j=1 ³√ σ2_y ij+ µ ´ =√σ2_y¯ 1+ µ   ¯y1= n1 X j=1 y1jとおく   同様に、 ¯ x2= √ σ2_y¯ 2+ µ とおけます。また、 c σ2 1 = 1 n1− 1 n1 X j=1 ³ x1j− ¯x1 ´2 = 1 n1− 1 n1 X j=1 ³ (√σ2_y 1j+ µ)− ( √ σ2_y¯ 1+ µ) ´2 = σ 2 n1− 1 n1 X j=1 ³ y1j− ¯y1 ´2 同様に c σ2 2 = σ2 n2− 1 n1 X j=1 ³ y2j− ¯y2 ´2 ですので、統合すると c σ2₌ 1 (n1− 1) + (n2− 1) ³ (n1− 1)cσ12+ (n2− 1)cσ22 ´ = 1 (n1− 1) + (n2− 1)  (n1− 1) · σ2 n1− 1 n1 X j=1 ³ y1j− ¯y1 ´2 + (n2− 1) · σ2 n2− 1 n1 X j=1 ³ y2j− ¯y2 ´2   = σ 2 n1+ n2− 2  Xn1 j=1 ³ y1j− ¯y1 ´2 + n1 X j=1 ³ y2j− ¯y2 ´2   = σ 2 n1+ n2− 2 2 X i=1 ni X j=1 (yij− ¯yi)2 *27_{xijを全て消して yijや µ, σ}2_{で表現する、ということです。}

となります。これより、２群比較のt統計量は t12= | ¯ x1− ¯x2| r³ 1 n1 + 1 n2 ´ c σ2 = s³ | ¯x1− ¯x2| 1 n1 + 1 n2 ´ 1 n1+n2−2 2 P i=1 ni P j=1 (xij− ¯xi)2 = |( √ σ2_y¯ 1+ µ)− ( √ σ2_y¯ 2+ µ)| s³ 1 n1 + 1 n2 ´ σ2 n1+n2−2 2 P i=1 ni P j=1 (yij− ¯yi)2 = √ σ2| ¯_y 1− ¯y2| s³ 1 n1 + 1 n2 ´ σ2 n1+n2−2 2 P i=1 ni P j=1 (yij− ¯yi)2 = √ σ2| ¯_y 1− ¯y2| √ σ2s³1 n1 + 1 n2 ´ 1 n1+n2−2 2 P i=1 ni P j=1 (yij− ¯yi)2 = s³ | ¯y1− ¯y2| 1 n1 + 1 n2 ´ 1 n1+n2−2 2 P i=1 ni P j=1 (yij− ¯yi)2 となります。µ, σ2がすべて消えていて、これはN (0, 1)からのサンプル{yij}におけるt統計量と全く同じ形です。  このように、N (µ, σ2_{)からの乱数でt統計量を作っても、N (0, 1)からのt統計量を作っても、全く同じ形になることが} 分かりました。 つまりどういうことかと言いますと、N (µ, σ2_{)をもとにしたt分布の棄却限界が知りたいとき、N (0, 1)をもとにしたシ} ミュレーションをすればよい、ということです。より具体的に言いますと、このようなt12をN (0, 1)のもとで1,000,000 個発生させて上側５％点をとれば、それが自由度(n1+ n2− 2)のt分布の有意水準５％の片側検定の棄却限界になる、と いうことになります。 Tukey-Kramerの方法とt検定の統計量の違いは、分散をpoolする部分と絶対値の有無だけで、上で議論した部分につ いては全く違いはないですのでTukey-Kramerの方法*28_{において、シミュレーションによって棄却限界を計算しようと} するとき、N (0, 1)からの乱数をもとに統計量を計算してもよいということが分かりました。

9

補足３：

Student

化した範囲の分布について

9.1

分布の導出と検定統計量との関係

Tukey-Kramerのところで出てきましたStudent化した範囲の分布についてみていきます。まずは定義を、それから Tukey-Kramerの統計量との比較を考えていきます。 *28_{Dunnett なども同じです。}

「定義：Student化した範囲の分布」 u1,· · · , uk, χ2 : 互いに独立 u1,· · · , uk∼N (0, 1), χ2∼χ2(ν)のときに、 QRk,ν = max i<j |u_qi− uj| χ2 ν のことをStudent化した範囲 とよび、このQR_k,ν の従う分布のことをStudent化した範囲の分布と呼ぶ*29_。  では、先のTukey-Kramerの統計量 Tij = | ¯ xi− ¯xj| r³ 1 ni + 1 nj ´ c σ2 とStudent化した範囲の分布の関係を見ていきましょう。実は、厳密な関係が成り立つのは例数が同じTukeyの方法の時 であり、Tukey-Kramerは近似的な方法なのです。このとき、ni= njをnとおくと、上のTij は Tij = | ¯xi− ¯xj| q¡₁ n+ 1 n¢ cσ 2 = | ¯xqi− ¯xj| 2 nσc 2 (1) となります。ここで今、xi1,· · · , xin∼N (µi, σ2)とxj1,· · · , xjn∼N (µj, σ2)（全て独立）が成り立つことから、 ¯ xi∼N µ µi, σ2 n ¶ , ¯xj ∼N µ µj, σ2 n ¶ です。ここで、帰無仮説Hij: µi= µj (= µとおく)のもとでは ¯ xi, ¯xj∼N µ µ,σ 2 n ¶ となります。これよりQR k,νと対応をつけるためには、x¯i, ¯xjを変形してN (0, 1)に従うようにしてやりたくなります。今、 ¯ x_qi− µ σ2 n ,x¯qj− µ σ2 n ∼N (0, 1) となることから、(1)に変形を加えてやりますと Tij = | ¯ xi− ¯xj| q 2 nσc 2 = |( ¯xi− µ) − ( ¯q xj− µ)| 2 nσc 2 さてここで、分子分母をともに q σ2 n で割ってやります。すると Tij = ¯¯ ¯¯ ¯ ( ¯xi−µ) q σ2 n −( ¯xj−µ) q σ2 n ¯¯ ¯¯ ¯ r 2 nσc2 σ2 n = ¯¯ ¯¯ ¯ ( ¯x_qi−µ) σ2 n −( ¯x_qj−µ) σ2 n ¯¯ ¯¯ ¯ q 2cσ_σ22

*29_{SAS の ProbMC 関数 (”Range”) によって求まるのは、この Q}R

となります。ここでいくつかチェックしましょう。 (i)分子の ( ¯xi−µ) q σ2 n と( ¯xj−µ) q σ2 n は、独立にN (0, 1)に従います。 (ii)分母のうち σ_σc22 の部分については、今回はν = (n1+ n2+ n3)− 3とおくと、νσc 2 σ2 ∼χ2(ν)になります*30*31。 (iii)分子と分母はi, jがいくつであっても互いに独立になります*32_。  これで分子は対応がつきました。分母について見ていきましょう。今 c σ2 σ2 = νcσ_σ22 ν より、(ii)と見比べていただくと、σ_σc22 はχ2(ν)に従う分布をνで割ったもの となりますので、 q c σ2 σ2 がQRk,ν の分母 q χ2 ν に対応します。

これより√2とmax以外の部分では、分子・分母ともに対応がつきました。そこで、T = maxi<jTij を考えることで、

T = max i<j Tij = max i<j ¯¯ ¯¯ ¯ ( ¯x_qi−µ) σ2 n −( ¯xj−µ) q σ2 n ¯¯ ¯¯ ¯ q 2cσ_σ22 =√1 2maxi<j ¯¯ ¯¯ ¯ ( ¯xi−µ) q σ2 n −( ¯xj−µ) q σ2 n ¯¯ ¯¯ ¯ q c σ2 σ2 =√1 2Q R 3,ν となります*33_。 従って、たとえばT の下側95％点が欲しい場合、QR 3,νの下側95％点を求めてそれを √ 2で割ってやればよい、という ことになります。

9.2

注意

最初にも述べましたとおり、今の議論は全群の例数が等しいとき、つまりTukeyの方法の場合です。例数が異なる

Tukey-CramerはTukeyの方法を例数が等しくない場合にも拡張した方法です。Tukey-Cramerでもこの手法でFWE

が保たれていることが知られています。 *30_{証明は略します。分散分析講義「第９回：多変量正規分布」の資料中に書いてありますので、興味ある方はそちらをご覧ください。} *31_{今回は３群でしたのでこのように書けます。一般には、ν = (全例数)}− (群の数) としたときに、νσc2 σ2 ∼ χ2(ν) となります。 *32_{これも証明は略します。これも「多変量正規分布」にほぼ同じことの証明が書いてあります。} *33_QR 3,νの「3」は薬剤群の数、ν は (全例数)-(薬剤群の数) です。

10

補足４：両側

Dunnett

の方法

(balanced)

の棄却限界について

Tukeyの方法では、わざわざ√2で割るのが面倒でした。そこで最初から√2で割った分布を基本に持ってくれば面倒が ない、ということでDunnettの方法においては、以下のQD2 k−1,νにおいて分母に √ 2を入れてあります。対照との比較であ ることを考慮して、以下の分布を定義します。なお、対照とするのは第１群とします。 「定義：Dunnettの両側範囲の分布*34_」 u1,· · · , uk, χ2 : 互いに独立 u1,· · · , uk∼N (0, 1), χ2∼χ2(ν)のときに、 QD2_k−1,ν = max 2≤i≤k |u1− ui| q 2χ_ν2 の従う分布をDunnettの両側範囲の分布と呼ぶ。  このようにおけば、先のStudent化した範囲の分布の議論とほとんど同じ議論で、QD2_k−1,ν の下側95％点が有意水準5 ％の両側検定の棄却限界となります。

*34_{「Dunnett’s two-sided range distribution」を適当に訳しました。もしかしたら他のよく知られた日本語名があるかもしれませんので、一応英}