ぶんさんぶんせき

分散分析 analysis of variance にまつわるノート。 入門者向けではなく,統計学にユーザとして接する院生向けです。

これまでに相談を受けた内容から集めました。多分に誤記の可能性があるので,発見次第ビシビシご指摘ください。 (ただ,厳密には誤りだと認識しつつ,説明の平易さのためにあえて大まかな記述をしている部分もあります。)

  1. since 2003-04-06
  2. updated on 2011-05-02

このノートでは,同じ意味でも様々な語を使っていることがあるので,その対応を挙げておきます。 統計基礎用語集(英和編)を参考に。

分散分析と因果

分散分析は実験計画法に密接に結びついた統計解析であるため, それが適用されるのは実験データであることがほとんどで,また, 分析結果から導かれた結論として因果関係に言及するケースが多い。

しかし他の解析手法,例えば重回帰分析などと同様に,分散分析であっても それは相関関係の分析,統計的結論妥当性 (statistical conclusion validity) の確保に過ぎない。 すなわち,この分析法自体の力によって因果関係についての知見がもたらされることはない。 分散分析で結論され得るのは「偶然とは言えないような変数間の関係性がある(要因水準間の違いがある)」 ということであり,それ以上ではない。

データを分散分析にかければ(そして有意な結果を得れば)因果が結論できると考えるのは大きな誤りである。 因果関係を結論できるかどうかは分析法ではなくデータの収集法に依存する。 これを常に心に留めておくこと。

ランダム サンプリング

これは分散分析だけに限らない話であるが,分散分析も含めて, 検定統計量に関する確率分布を利用した検定をするということは,母集団についての推論を行っている。 すなわち,その場にある被験体群(標本)のみに関するものではない

標本から母集団について統計的に正確な推測を行うためには,ランダムサンプリングが欠かせない。 ランダムサンプリングとは,母集団のどの要素もゼロでない既知の確率で標本に含まれうるような,サンプリング方法である。 たとえば単純ランダムサンプリングでは,母集団(実践的にはサンプリングフレーム)の可能なすべての標本(要素のすべての組み合わせ)が等確率で選ばれる。すなわち抽出確率は既知。 ランダムサンプリングをしていなければ母集団についての正確な統計的推測はできないから, 分散分析などを適用しても母集団について正しい結論を得ることを一定の範囲で保証することができない。

ただし,心理学や医学などでの通常の適用場面では,ランダムサンプリングはしていないが,ランダム割当て(random assignment)を行っていることがある。 このような場合,一番適切なのは randomization test (無作為化テスト) であるが, 分散分析などの方法を適用しても,randomization test による正確なp値の近似値を得ることができる。 よって,ランダムサンプリングをしていなくても,ランダム割当てされたデータへの近似的な方法として 分散分析の適用は認められるのである。 ただし,近似であるという点を忘れてはいけない。このような近似的な方法の適用が行われてきたのは, randomization test を実行とすると莫大な量の計算が必要で,実行が難しかったからである。 ところが,近年は計算機の性能も飛躍的に上昇してきたから, randomization test を実行して正確なp値を得ることは難しくない。 そうすると,そのような場合の分散分析はお払い箱なのである。

分散分析は,母集団からランダムサンプリングをした場合という正当なケースや, サンプルサイズが非常に(randomization test ができないほど)大きいときの近似的代用法として などで適用されるべき。

理論的前提と平方和の分解

分散分析の根底にある基本中の基本。

あるいくつかの仮定の下で,全体の平方和は各コンポーネント(成分)に分解できる (=全体変動を各変動成分の単純な和として表せる)。

その仮定の一つが,線型性(linearity)。 もうひとつ,基本的な仮定として注目しておくべきなのは,直交性(orthogonality)。 そして三つ目に,なぜだか教科書でこれだけが強調されている(それゆえ有名な),誤差に関する諸々の仮定

線型性

言わずもがな,分散分析モデルは線型モデルの一種である。 分散分析がこの線型性の仮定に乗っかっているということは, モデルの線型性が前提できない現象へ分散分析を根拠なく適用してはいけない ということ。これはとっても忘れられがち。 もう少し実践的に言うと,線型モデルを当てはめるのがマズイと,事前に,あるいはデータの分布を見て,わかる場合は, もっと適切な分析モデルを当てはめるべきである,ということだ。 これは「何でもかんでも分散分析するな」の根拠の1つである。

ここで言っている線型性とは,独立変数と従属変数の関係の話ではなく,母数(parameter)についての線型性であることに注意。 式で書けばこういうこと。

y_i = sum_{j=1}^{p}x_{ij}beta_j + varepsilon_i qquad qquad i = 1,...,n

行列表記にすると,

y = X beta + varepsilon

ここで,y は被説明変数ベクトル(n×1),X は説明変数行列(計画行列とも呼ばれる)(n×p), beta は母数ベクトル(p×1),varepsilon は誤差変数ベクトル(n×1)。 ただしこれは母数モデルの場合。変量モデルの場合は,beta は変量効果を表す確率変数ベクトル。

たとえば,2要因デザイン(それぞれ2水準,交互作用なし),繰り返し数 3 の場合,次のようになる。

((y_1),(y_2),(vdots),(y_12)) = 
((1,1,1),(1,1,1),(1,1,1),
 (1,1,0),(1,1,0),(1,1,0),
 (1,0,1),(1,0,1),(1,0,1),
 (1,0,0),(1,0,0),(1,0,0))
((beta_0),(beta_1),(beta_2)) + 
((varepsilon_1),(varepsilon_2),(vdots),(varepsilon_12))

このように dummy code で計画行列を作った場合は母数の解釈に注意する必要がある。 例えば, beta_0 は総平均ではない。 母数について,ふつうの分散分析に近い解釈をしたい場合は, effect coding あるいは orthogonal coding を使うべし。

誤差分散

説明変数以外に帰属される分散(=説明変数に帰属できない分散)は誤差分散と呼ばれる。 この誤差分散に関しても,分散分析モデルではいくつか仮定を置いている。

正規分布

誤差は平均 0 ,分散 sigma_e^2 の正規分布に従うと仮定される。これは実践としてはかなり強い仮定に思われる。 なぜなら,実際に扱われるのは誤差がきっちり正規分布しているようなデータばかりではないから。

ただし,分散分析においては,この仮定は有意確率を精密に計算するためにあるようなものなのである。 というのは,確かに,誤差が正規分布であれば検定統計量であるF値はF分布に従うが, そうでなければF値がどのような分布を取るか明確でない。 しかし,F値は比であるから,誤差がどのような分布であろうともF比の分子が分母に比べて非常に大きければF値が有意であることは間違いなく, 微妙な計算を要求される場合しか誤差の分布型は問題にならない。そういう意味で,この点に関しては分散分析は頑健である。

i.i.d.

i.i.d.は independent and identically distributed の略で,そのまま,独立同一分布のこと。 分散分析においては,誤差分布はすべて独立で,同一の分布に従うと仮定されている。

上記のとおり,誤差分布は平均 0 の正規分布だという仮定があるが, 「同一」の分布なのだから,さらに分散もすべて同一だと仮定されている。これがいわゆる「等分散性の仮定」である。

もうひとつの独立性の仮定は,誤差についてこれまでに挙げた仮定以上に重要であるにもかかわらず, 初歩の教科書で注目されていることが少ない。 ここで,誤差が独立であるとは,誤差間に何の関係もないということ。 例えば,調査一人目 A さんについて何らかの属性を測定したとしよう。この値は母集団の平均値とまず一致しない。 すなわち,平均値からの誤差(個人差,測定誤差,etc.)が含まれている。 プラスの誤差かもしれないし,マイナスかもしれない。大きいかもしれないし,小さいかもしれない。 その確率は,正規分布に従っていると分散分析では仮定する。 調査二人目 B さんについても測定した。B さんの値にも誤差が含まれているだろう。 ここで,B さんの値にどのくらいの誤差が含まれているかの確率は,A さんにおけるそれに関係なく(A さんの誤差がどのくらいであっても一切左右されずに), 同じ分布に従っていると仮定される。三人目 C さんの誤差も,A さんの誤差とも B さんの誤差とも関係ない。以下同様。これが独立な場合。 逆に,独立でない場合というのは,A さんの誤差の値によって B さんの誤差の分布が違ってくるということである。

誤差が互いに依存していると分散分析の結果に大きく影響するので非常にマズいのだが, しばしば行われる同一被験体の反復測定などでは相関することのほうが多い。 特に気をつけよう。

直交性

上記の前提の1つ,直交性について。直交した釣り合い型データ(balanced data; 各要因水準のデータ数がきっちり同数のデータ) であれば特に問題はないが, 非釣り合い型データ(unbalanced data; 各要因水準のデータ数が不揃いのデータ)の場合, 直交性が崩れることがしばしば。

ここでいう直交とは,説明変数ベクトル同士の内積がゼロ(ベクトル間の角度が90°)という意味である。 これはすなわち,説明変数が互いに相関していないということ。

この直交性の定義から考えると,データ数が非釣り合いであっても直交し得る場合があるということがわかる。 つまり,各セルのデータ数が水準によって違ったとしても周辺度数には比例する場合,内積はゼロであるから直交である。 クロス表における独立性の話を思い浮かべてもらうとわかりやすい。

平方和のタイプ

分散分析は因子が直交している場合に成り立つ平方和の加法性を利用して平方和を分解するものであり, そもそも釣り合い型データに適用するよう考案された分析である。 ただ,直交性が崩れていると即アウトというわけではない。 もしそうなら現実にはアウトのデータが無残にも山ほどできてしまう。 そこでアンバランスでもなんとか個々の変動要因の平方和を計算する方法を考える。 その案がタイプI~IVの平方和というやつ(たぶん名前はSAS起源)。 だからこの「平方和のタイプ」は幸運にも自分のデータが直交釣り合い型データのときは まったく関係ない(どれでも同じ)。

基本的に,タイプの違いというのは,「非直交なデータでは要因間に相関があることになるから, 平方和の加法性が単純には成立しないので,その相関の分を如何にコントロールするか」, というコントロール方法の違いに拠る。ここでいう「コントロール」という概念は 偏相関係数や偏回帰係数におけるそれと同様だと思ってもらってよい。

Type I SS (タイプ I 平方和)

複数の変動要因の効果(主効果,交互作用)を1つずつモデルに追加していく時, モデルの平方和の増加分をその効果の平方和と考える。このことから,「逐次平方和」とも呼ばれる。 言い換えれば,ある要因の効果は, それ以前にモデルに投入されたすべての要因によってコントロールされた状態で検討される,ということ。

タイプ I の平方和は,要因の投入順序により異なった値を取る。 ただし,上で述べたように,直交デザインの場合は投入順序には依存しない(コントロールする分がないから)。 低次の効果から順に投入する(主効果,二重交互作用,三重交互作用…)というルールは守るように。

Type II SS (タイプ II 平方和)

タイプIと異なり,タイプIIでは同レベルの効果はモデルに一斉投入される。

主効果に関しては,まず,すべての要因の主効果だけを投入したモデルにおいて全体の平方和を計算する。 次にそこからある要因の主効果を取り除いたモデルの平方和を計算する(当然ながらこちらのほうが小さくなる)。 そしてその減少分をその要因の主効果の平方和とする。 同じ手順で,全要因についての主効果が計算される。 これは,ある要因の主効果が,それ以外のすべての要因の主効果でコントロールされた状態で検討されるということ。

一次交互作用に関しては,まず,すべての主効果とすべての一次交互作用を投入して平方和を計算する。 次にそこからある交互作用を取り除いたモデルの平方和を計算する。 そしてその減少分をその交互作用の平方和とする。 これは,ある一次交互作用が,それ以外のすべての一次交互作用とすべての主効果で コントロールされた状態で検討されるということ。 より高次の交互作用についても同様。

タイプIIの平方和は,投入順序には依存しない。

Type III SS (タイプ III 平方和)

タイプIIIの平方和では,すべての効果がモデルに同時投入される。 主効果,交互作用に関わらず,ある効果の平方和は,全効果を投入したモデルからその効果を取り除いたときの, 平方和の減少分として計算される。 これは,ある効果が,それ以外のすべての効果(主効果,交互作用)でコントロールされた状態で検討されるということ。

タイプIII の平方和も,投入順序には依存しない。

タイプIIとタイプIIIは両方とも「偏平方和」と呼ばれるが,これらの違いは, 同レベルとそれより低次の効果だけでコントロールする(タイプII)か, 全効果でコントロールする(タイプIII)か,という点。 2要因デザインを例にすると,ある要因の主効果の平方和の算出において, もう一つの要因の主効果以外に交互作用でもコントロールするか(タイプIII)しないか(タイプII), という違い。

Type IV SS (タイプ IV 平方和)

これはタイプIIIの類似版で,データが1つもないセル(空のセル)が存在するときに使う(それ以外は使う必要がない)。 空のセルが存在すると,そのセルについて何もわからないのだから,ふつうの分散分析を適用するのはよくない。 とくにタイプIIIは交互作用でコントロールするから。 よって,そのような場合,データが存在するセルだけでできる範囲の比較を行い, そこから各要因の主効果や(可能ならば)交互作用を推定するという方法が(統計ソフトの内部で)取られる。 その際に計算される平方和である。 コントロールについてはタイプIIIと同様に考えているので,空のセルがなければタイプIIIと一致する。

特徴の比較と選択

繰り返しになるが,釣り合い型データではどれでも同じである。 非釣り合い型だが直交したデータの場合,タイプIとタイプIIは一致するが, タイプIIIは交互作用で調整するので一致しない。 非直交なデータであっても,交互作用がないモデルの場合は,タイプII,タイプIII,タイプIVは同じである。 これらは平方和がどうコントロールされるか考えれば明らか。

交互作用に関して言うと,タイプIIは主効果の計算において交互作用を無視する。 タイプIIIは交互作用の大きさを「バランスよく」算定し(dummy codeでなくeffect codeをイメージせよ), 主効果の検討において交互作用を考慮する。

逆に,セルサイズに関して言うと, タイプIIでは各セルのデータ数で重みをつけた加重平均が周辺平均として計算され,主効果が検討される。 それに対し,タイプIIIでは,セルサイズで重みづけをしない(=等しく重みづけた)平均が計算されることになる。 でも空のセルがあっても等しく重みづけるというのが変だということで, タイプIVなんてのが出てきたりするわけだ。

このようにタイプIIもタイプIIIもそれぞれの言い分があるから, どっちがいいかというのは一概に議論が決着しないところ。ユーザは各自の判断で使うしかない。

要因間に因果的な順序が想定される場合などは,タイプIを使うだろう。 それ以外はたいていタイプIIかタイプIIIだが,重みづけをするべきかどうかは, そのデータの背景や研究分野の理論などとも相談しなければならないだろう。 主効果に主たる関心がある場合はタイプII, 交互作用に関心がある場合はタイプIIIという考え方もあるらしい。

青木さんのサイトで解説されている計算方法 はタイプIIの平方和。

以下,参考にできるページ。

固定効果と変量効果

分散分析における線型な統計モデルに組み込まれ得る要因効果は, 固定効果と変量効果の2種類に分類される(2種類あることに注意!)。 固定効果だけのモデル,変量効果だけのモデルもあり得るが,両方入ることもある(「混合効果モデル」を参照)。 おそらく心理学者にとっては固定効果のほうが馴染み深いだろう。 なぜならユーザ向け統計学の教科書の中でこれらそれぞれについて解説している分量は, 固定効果についてのほうが変量効果についてよりも圧倒的に多いからである。 しかし,もちろんながら,どんな教科書にも変量効果について一切書いていないというわけではないので, これを知っているかどうかが,きちんと勉強しているかどうかの指標になり得る。(という現状が悲しい。)

固定効果(fixed effect; 母数効果とも呼ばれる) とは,扱われている要因の各水準は(同じ実験を繰り返しても)一定の固定された(=変わらない)効果を持つ と想定されるタイプの効果である。 例えば,新薬A,既存薬B,プラセボ薬Cをそれぞれ別の被験者に投与して, 新薬の効果を見るという実験はこれに当たる。 「与えた薬剤」という要因の各水準(薬A,薬B,薬C)の平均値やそれらの差の大きさそのものに興味があり, 水準としてその特定の薬が選ばれるのが必然(それを調べたいのだから)という状況である。 この場合,同じ実験を繰り返すことは,同じ薬剤を投与することに他ならない。

一方,変量効果(random effect) とは,扱われている要因の各水準が(その実験毎に) ランダムに選ばれたものであると想定されるタイプの効果である。 例えば,体操競技の採点に審査員による違いがあるかどうかを検討するという実験において, 審査員は現存する審査員全体からランダムに一定人数選ばれたとすれば,これに当たるだろう。 そこでは,ある特定の審査員の採点の平均値にも,別の特定の審査員との平均値の差にも興味がない。 興味があるのは,全体として審査員の採点には個人差があるか(甘い審査員,辛い審査員がいるか)ということや, もし個人差があるなら,その要因(審査員の個人差)による分散が競技得点の全体の分散に占める割合は どれくらいか,ということであろう。 この場合,実験を繰り返すときには,その実験ごとにランダムに審査員を選んでくる,すなわち, 要因の水準1,水準2,水準3…として誰が選ばれてくるかは実験ごとに変わることになる。

要するに,伝統的分散分析においては,変量因子では各水準の値は確率変数の値であり, その水準は平均 0,分散 sigma^2 の正規分布からランダムに選ばれたものと仮定される。 分散は固定値であるが未知であり,これが推定される。 これに対し,固定因子では,各水準の値は固定値(母数)であり, 各水準の値をすべて合計すると 0 であるという前提(というか制約)がついた上で それぞれの値が推定される。

固定効果と変量効果の区別は,主効果にのみ関する話ではない。交互作用にもこれらの区別は同様に適用される。 特に覚えておくべきは,交互作用を構成する各因子に1つでも変量因子が含まれている場合, その交互作用も必ず変量である。が, ある交互作用が変量であるとしても,各主効果や低次の交互作用に変量が含まれているとは限らない。 すべて固定効果の因子であっても,それらの交互作用が変量であるというモデルもあり得る。

さて,実践において問題となるのは,自分の実験で扱っている因子が固定因子なのか変量因子なのか という判別であろう。これは,その因子が何であるか(審査員,薬,etc.)では決まらない。 むしろ,研究者が何を調べようとしているか(そしてその結果どのような水準の決め方をしたか)で決まる。 上に述べたように,変量効果を持つ因子の水準はランダムに決定されたと見なされるものであり, それらの個々の水準には興味がない。対して,固定効果を持つ因子の水準は, 特定の水準をその効果を見るために選ぶものである。

例えば,ハングルに関して文字の違いによって覚えやすさに違いがあるかを調べる実験をするとしよう。 全ハングル文字の中からランダムにいくつか選んできて全体として違いがあるかを実験するならば, この「文字の違い」という因子は変量因子である。 しかし,特定の文字Aと文字Bと文字Cの間に違いがあるかを調べたい(その特定の文字でなければならない)ときは, この因子は固定因子である。 判断のヒントになるのは,同じ目的で同じ実験を反復するときに,水準が変わるかどうか。

心理,教育,医学の実験でよくある「被験者」という因子(ようするに乱塊法,反復測定,被験者内配置)は, たいていの実験においては変量因子である。これが固定因子である実験というのは, 特定のAさんとBさんの違いを調べたい,という場合である。 (私が知る限りこれらの領域でそんなケースは滅多にない。) だから,たいていの人はこれを変量因子として分析しているはずであり(そうでなければ誤用だろう), 変量因子は実は非常に身近なものである。

帰無仮説の違い

固定効果と変量効果では検定における帰無仮説が異なる。

固定効果における帰無仮説とは,要因 A の水準数を p とすると,

H_0: a_1 = a_2 = ... = a_p

というもので,各水準の効果(母数 i.e. 未知の固定値)がすべて等しい, もっと言えば,各水準の平均値がすべて等しい,という仮説である。

もちろん,分散分析においては,

sum_{i=0}^{p}a_i = 0

という制約がつくので,上の仮説は,

H_0: a_1 = a_2 = ... = a_p = 0

ということである。

これに対し,変量効果における帰無仮説とは,

sigma_A^2 = 0

であり,要因A(確率変数)の分散がゼロである,という仮説である。

F比を構成するもの

分散分析において効果の検定に用いられるF値は,分母と分子で構成される。すなわち,比が計算される。 よくある分散分析の利用場面では,分母には誤差項の平均平方,分子には主効果や交互作用の平均平方があてられることが多い。 しかし,これは常にそうであるわけではない。

上述の「よくある利用場面」,たとえば一元配置分散分析において,分母に誤差項の平均平方 MS_e , 分子にA因子の主効果の平均平方 MS_a で比が計算される理由は, その場合の MS_e の期待値が sigma_e^2MS_a の 期待値が sigma_e^2 + n sigma_a^2 だからである。 両者の違いは sigma_a^2 の含まれる第2項のみであり, sigma_a^2 が 0 ならば違いはなくなる。 この比を計算したとき,もし sigma_a^2 = 0 が真(Aの効果がない)ならば,比は 1 になるし, sigma_a^2 > 0(Aの効果がある)ならば,比は 1 より大きくなる。 これを利用して,F比が 1( sigma_a^2 = 0 のケース)よりも偶然といえないほど離れているかどうか, 確率を計算する。

この「よくある利用場面」とは母数モデルを指している。母数モデルにおいては,いくつ要因が入った計画であっても, 各効果の平均平方の期待値は 「 誤差の分散 sigma_e^2 + その効果の分散×K 」 となる(Kは定数,水準数やセル内のデータ数で決まる)。 よって,分母に誤差項の平均平方,分子に各効果の平均平方をもってきてF比を作ればよい。 しかし,変量効果が入るときにはそうはいかない

以下に,デザインの例と,そのときの各平均平方の期待値を挙げる。 表中の delta_A などは,その要因が変量効果の場合 1 ,固定効果の場合 0 である。 すなわち,delta の付いている項は delta の添え字の要因が固定効果であれば消える。 p,q,r はそれぞれ要因A,B,Cの水準数,nはセル内のデータサイズである(釣り合い型データの場合を考えている)。

一元配置
変動要因自由度平均平方の期待値
A p - 1 σ2e + n σ2A
誤差p(n - 1)σ2e
全体pn - 1
二元配置
変動要因自由度平均平方の期待値
A p - 1 σ2e + δB n σ2A×B + qn σ2A
B q - 1 σ2e + δA n σ2A×B + pn σ2B
A×B (p - 1)(q - 1)σ2e + n σ2A×B
誤差pq(n - 1) σ2e
全体pqn - 1
三元配置
変動要因自由度平均平方の期待値
A p - 1 σ2e + δBδC n σ2A×B×C + δB rn σ2A×B + δC qn σ2A×C + qrn σ2A
B q - 1 σ2e + δAδC n σ2A×B×C + δA rn σ2A×B + δC pn σ2B×C + prn σ2B
C r - 1 σ2e + δAδB n σ2A×B×C + δA qn σ2A×C + δB pn σ2B×C + pqn σ2C
A×B (p - 1)(q - 1) σ2e + δC n σ2A×B×C + rn σ2A×B
A×C (p - 1)(r - 1) σ2e + δB n σ2A×B×C + qn σ2A×C
B×C (q - 1)(r - 1) σ2e + δA n σ2A×B×C + pn σ2B×C
A×B×C(p - 1)(q - 1)(r - 1)σ2e + n σ2A×B×C
誤差pqr(n - 1) σ2e
全体pqrn - 1

このように,変量効果がモデルに入ってくると,母数モデルの場合とは平均平方の期待値が異なり(deltaの付いた項が加わる), その結果,F比を作るために使う平均平方が変わったり適切なF比を容易に作れなかったりする場合がある。 F比が作れないということは,古典的な分散分析のやり方で検定できない,ということ。 複雑な計算をしたり,後述の混合効果モデルを使ったりする必要がある。

次にいくつか具体例を挙げてどのような分散分析表になるか見てみる。

乱塊法

2因子乱塊法(randomized block design)では,次のようなモデルを設定する。

y_{ijk} = mu + A_i + B_j + (A xx B)_{ij} + C_k + varepsilon_{ijk}

説明変数として組み込まれているのは要因A,要因B(およびそれらの交互作用),ブロックC,である。 ここでは,ブロックと効果を調べたい要因(A,B)との交互作用はないものと仮定されている。 この式を見ればわかるように,ブロック(例えば「被験体」)というのは,1つの因子としてモデルに組み込まれる (研究者がその効果を調べたいものではないかもしれないが)。 すなわち,これは3要因計画の一種である。ブロックと他の要因との交互作用がないと仮定される特徴があるとしても, そのようなモデリング(交互作用項を外すとか)は乱塊法に限らず行われ得るので, 乱塊法も3要因計画の文脈でより一般的に考えることができる。

ただし,注意しなければならないのは,ここに組み込まれている要因(A,B,ブロック)が変量因子かどうか, ということ。 効果を調べたい要因A,Bは固定因子であることが多いと思われるが,ブロック因子は変量因子である場合もある。 例えば「被験者」というブロック因子は典型的に変量因子だと見なされる。 上記の分散分析表の一覧をもとに,このようなブロックが変量因子の場合の分散分析表を書いてみよう。 セル内のデータサイズ(n)が 1 の場合(繰り返しなし)が心理学では典型的なようなのでこの場合を考える。

上記の三元配置の表をもとに n = 1 を入れて書き直す。そうすると,誤差の自由度が 0 と計算される。 つまり,繰り返しがないので A×B×C のフルモデルを当てはめてしまうと誤差が同定できない。 しかし,この場合,ブロックC と他の効果との交互作用は無いものと仮定されているので,変動要因AC,BC,ABCの行は消され, これらの和が誤差の平方和と自由度となる。 さらに,平均平方の期待値について,ブロックC と他の効果との交互作用は無いという仮定だから,これらの入っている項を消す。 delta も C に関係するもの以外は 0 なので,これの付いた項は消す。 するとこのような表になる。

変動要因自由度平均平方の期待値F比
A p - 1 σ2e + qr σ2AMSA/MSe
B q - 1 σ2e + pr σ2BMSB/MSe
C r - 1 σ2e + pq σ2CMSC/MSe
A×B(p - 1)(q - 1)σ2e + r σ2A×BMSA×B/MSe
誤差(pq - 1)(r - 1)σ2e
全体pqr - 1

表を見ればわかるとおり,このデザインではどの効果も誤差の平均平方で割ってF値を出せばよい。

次に,ブロックと要因A,Bとの間に交互作用があると考えられる場合にどうなるかを見てみる。 モデルはこのようになる。

y_{ijk} = mu + A_i + B_j + (A xx B)_{ij} + C_k + (A xx C)_{ik} + (B xx C)_{jk} + (A xx B xx C)_{ijk} + varepsilon_{ijk}

ここでは,繰り返しなしの場合を考えているのでA×B×Cの二次交互作用と誤差を分離できないが, 繰り返しがあれば分離することができる。

先ほどと同じく,上記の三元配置の表をもとに n = 1 を入れて書き直す。 繰り返しが無いので誤差の平方和を取り出せない。 さらに,平均平方の期待値について,delta_C 以外のdeltaの付いた項は消す。 するとこのような表になる。

変動要因自由度平均平方の期待値F比
A p - 1 σ2e + q σ2A×C + qr σ2AMSA/MSA×C
B q - 1 σ2e + p σ2B×C + pr σ2BMSB/MSB×C
C r - 1 σ2e + pq σ2C
A×B(p - 1)(q - 1) σ2e + σ2A×B×C + r σ2A×BMSA×B/MSA×B×C
A×C(p - 1)(r - 1) σ2e + q σ2A×C
B×C(q - 1)(r - 1) σ2e + p σ2B×C
A×B×C(p - 1)(q - 1)(r - 1)σ2e + σ2A×B×C
全体pqr - 1

ここで注目すべきは,変動要因AとBの主効果やA×Bの交互作用を検討するときに作るF比。 この表にある平均平方の期待値を見ればわかるように, sigma_A^2sigma_B^2 が 0 かどうかを検定するには, 分母にそれぞれA×C交互作用の平均平方,B×C交互作用の平均平方,を用いる。 A×B交互作用の検定には分母にA×B×C交互作用の平均平方を用いる。これら以外の効果については, どれか2つの平均平方の組み合わせで適切なF比をつくることができないもちろんこれはブロック (例えば「被験者」) を変量因子と考えているからで,固定因子であればこのようなことにはならない。

SPSSなどではデフォルトでこの「ブロックとの交互作用あり」のモデルによって分析される。

分割区画デザイン

分割区画デザイン (split-plot design) の場合も乱塊法と同様,変量因子に注意しなければならない。 例として,ブロック間1要因,ブロック内1要因の繰り返しなしデザイン(ブロックは変量因子)を考える。モデルは次のとおり。

y_{ijk} = mu + A_i + B_j + (A xx B)_{ij} + C_{k(i)} + (B xx C)_{jk(i)} + varepsilon_{ijk}

Aがブロック間要因,Bがブロック内要因,Cがブロックである。 C の添え字に()が使われているのは,ブロックCが要因Aにネストしている,という意味 (次の「分散分析表と平均平方の期待値」の節も参照)。 ブロックCは要因Aとクロスしていない。 すなわち,ブロックCの各水準と要因Aの各水準の全組み合わせについてデータが取られているわけではない。 要因Aの各水準の下にブロックCの水準のいくつかが割り当てられているのである。だからこそ,要因Aはブロック間要因,と呼ばれる。

具体例を見るとよりわかりやすいだろう。 要因Aは2水準でそれぞれブロック2つ,要因Bも2水準,繰り返しなし,の場合の水準の割当てを示してみる。

{:
(y_1:, A_1, C_1, B_1),
(y_2:, A_1, C_1, B_2),
(y_3:, A_1, C_2, B_1),
(y_4:, A_1, C_2, B_2),
(y_5:, A_2, C_3, B_1),
(y_6:, A_2, C_3, B_2),
(y_7:, A_2, C_4, B_1),
(y_8:, A_2, C_4, B_2)
:}

これを見るとたしかに A はブロック間, Bはブロック内の割当てになっていることが確認できるだろう。 ここではあえてブロックC に 1~4 の添え字を付けたが, 要因A の水準ごとに 1, 2, 1, 2 と付けることもできる(クロスの場合みたく)。 しかし,そうした場合,1つ目の C_1C_2 と 2つ目の C_1C_2 とは同じブロックではない(同じ効果が期待されない)。 よってこの場合の C には,C_{k(i)} という風に C の添え字だけでなく A の添え字の条件付けが必要となる。

上記の一覧にある三元配置分散分析表から,このようなデザインの場合の分散分析表を作ると,下のようになる。 C(A) すなわち A にネストした C の効果は,変動要因C と AC の和になる。C の主効果が A×C の交互作用に吸収されたという イメージで考えるとよい。実際,C 単独の主効果は無い(取り出せない)ものと前提され, そうして膨らんだ A×C 交互作用を C(A) の効果として扱う。B×C(A) の効果も同様に,B×C と A×B×C が合わさる。 また,ここでも n = 1 の場合を考えているので,誤差は分離できない。

変動要因自由度平均平方の期待値F比
A p - 1 σ2e + q σ2C(A) + qr σ2AMSA / MSC(A)
B q - 1 σ2e + σ2B×C(A) + pr σ2BMSB / MSB×C(A)
A×B(p - 1)(q - 1) σ2e + σ2B×C(A) + r σ2A×BMSA×B / MSB×C(A)
C(A)p(r - 1) σ2e + q σ2C(A)
B×C(A)p(q - 1)(r - 1)σ2e + σ2B×C(A)
全体pqr - 1

見てわかるとおり,要因A の主効果を検定するには,誤差の平均平方でなく C(A) の平均平方で割り, 要因Bの主効果,A×B の交互作用の検定には B×C(A) の平均平方で割る必要がある。 その他の効果を検定するための適切なF比は作れない(繰り返しがあれば可能)。

分散分析表と平均平方の期待値

もっと一般的に,任意数の要因の組み合わせにおいて平均平方の期待値がどのようになるかを見てみたい人は, 下の可変テーブルが使える。いろいろ変えて試してみてください。

要因数,各要因が固定効果か変量効果か,ある要因が別の要因にネストしているかどうか,を選択できる。 固定効果/変量効果の選択における delta というのは,上の説明における delta と同じである。 固定,変量を選べばこの記号がなくなるので見やすいかも。

階層のチェックは,「(行)が(列)にネストしている」(入れ子になっている)という意味である。 例えば「B が A にネストしている」とは,B の各水準はAの各水準ごとに異なっているということ。 個体を4つ選んで,それぞれについて5つ体細胞を取り,薬品検査をしたなら,検査結果の変動に対する体細胞の要因は個体の要因にネストされている。 こういうデザインの実験は,「枝分かれ実験」とか「階層計画 (hierarchical design) の実験」とも呼ばれている。

上述の乱塊法や分割区画の場合を自分で作ってみてほしい。

(JavaScript で書いているのでこれを利用するならスクリプトを有効にしてください。)

混合効果モデル

混合効果モデル(混合モデル)とは,単純には,固定効果と変量効果を両方含むモデルのことである。 これは理論的には古くから考えられてきたが,実践的な計算では変量効果を扱う上で限界があった。

だが,近年の統計学とコンピュータの発展によって(ようするに最尤法が実用できるようになった),両方の効果が入ったモデルを適切に扱うことができるようになった。 それだけでなく,相関のあるデータ(欠損値がある場合や,変量因子間や誤差間の共分散)も扱うことができる。 このような拡張があることから,分散分析モデルの一種としての範囲を超えて, 別な手法として「(線型)混合モデル」と呼ばれ確立されているようだ。

詳しくは下節「線型混合効果モデル」を参照のこと。

交互作用の検討

交互作用効果が有意であった場合,全体的 (over all) 効果だけで満足せずに, 交互作用の詳細なパターンを調べようとすることがしばしばある。

単純効果

心理学・医学・教育学などでは,交互作用が有意だった場合に,事後的に単純効果 (simple effects) の検定が行われることが多い。 これは交互作用効果の下位検定と呼ばれる。 しかしこれは,交互作用効果の詳細なパターンを確認する手段の1つに過ぎないし, その前の分散分析での交互作用効果の検定と整合しない。その意味で,厳密に「下位」ではない。 例えば,交互作用が有意でないのに単純主効果は有意であったり,交互作用が有意なのにどの単純主効果も有意でなかったり, などという場合があり得る。

「下位検定」は英語にすると subeffect tests だろうか。しかし私はこの語を自分の読んだ論文/書籍で見かけた覚えがない。 それに比べて日本語の文献では「下位検定」をよく見かける。日本人が作った用語なのだろうか? post hoc analysis, post hoc test という言葉ならしばしば聞くが,これを訳すなら「事後分析」「事後検定」だろうし。

ある要因の水準を限定して,それに当てはまる部分のデータだけでの効果というものを考えたのが,単純効果という概念である。 そして,これを検定するのが単純効果の検定である。 例えば,2要因(それぞれ3水準)の完全無作為計画で実験をした場合, 有意な交互作用効果が得られ,その後,要因Aの水準a1のみに限って(水準a2,a3のデータは除いて), 要因Bの効果を検討しようとするかもしれない。 これは要因Aの水準a1における要因Bの単純主効果 (simple main effect) と呼ばれる。 合わせて,要因Aの水準a2やa3においても要因Bの効果を検討する。 これらの結果に基づいて,交互作用がどのようなパターンであるかの分析を試みることが多い。 繰り返すが,これは交互作用の詳細なパターンを調べる(これだけで完璧ではない)手段のうちの1つである。 このような単純主効果の検討は,簡単に言えば,主効果と交互作用効果を合わせた効果を調べている。 その検定結果がどうなるかはケースバイケースであり(すなわち,交互作用の有意性だけでは決まらず), 前述したように「交互作用効果の検定と整合しない」。

単純効果という概念には,単純主効果だけでなく,単純交互作用効果 (simple interaction effects) も含まれる。 これも前述と同様に,ある要因の水準を限定して,それに当てはまる部分のデータだけでの交互作用効果,というものを考える。

単純効果検定の誤差項

単純効果の検定を行う場合,基本的には,水準を限定した一部のデータ(サブセット) に対してもう一度分散分析を行うことになるのだが,その際の誤差項の扱いに注意。

チェックポイントは,その前の全データでの分散分析において, (いま単純効果を調べようとしている)変動要因の効果の検定に用いた誤差項と, それより高次の交互作用効果の検定(おそらく単純効果の検討を動機づけただろうソレ)に用いた誤差項とが, 同一であったかどうか。

これが同一ならば,迷うことなく,その誤差項の平均平方を単純効果の検定に用いるF値の分母として使えばよい (分子はもちろんサブセットにおける関心のある変動要因の平均平方)。 例えば,上述のとおり,すべて固定効果の完全無作為計画ではどの効果を検討するにも同じ誤差項を使うので,この場合に該当する。

同一でない場合は論争に巻き込まれる。例えば,上述のように分割区画デザインの場合は変動要因によって 用いる誤差項が変わるため,この場合に該当する。こちらの場合に,誤差項をどのように計算するかについて, 心理学界や生物学界で受け入れられている方法が2通りある。 ひとつはプールされた誤差項を使う方法,もうひとつは水準別誤差項を使う方法。

プールされた誤差項

検討しようとしている単純効果に関係する誤差項が複数ある場合, それらを全部まとめて1つの誤差項を計算するという手が考えられる。 これが「プールされた誤差項」 (pooled error terms) と呼ばれる方法。 「まとめる」というのは,単純に,複数の誤差項の平方和と自由度をそれぞれ合計した後で割り算して 平均平方を求めるというもの(要するに加重平均)である。 しかし,異質なものを合併していることからこれだけでは済まず, 自由度の調整 (Winer, Brown, & Michels, 1991) あるいは修正F分布表の使用 (Kirk, 1982) などの対処が合わせて必要だとされている。

水準別誤差項

もっとシンプルで理解しやすいのは,サブセットについて分散分析を適用し, その結果をそのまま採用するという方法。この方法だと, 誤差項は対象となる一部のデータのみをもとに新たに(検定毎に)計算される (すなわち誤差項の平均平方がサブセットごとに変わり得る)ので,水準別誤差項(separate error terms)と呼ばれる。

こちらの方法を適用する場合は事後分析という前提は必要ない。

水準別誤差項のほうがわかりやすいが,わざわざプールされた誤差項という込み入ったことをする理由の1つは 「できるだけたくさんのデータから推測した方が精度が高くなる」ということのようだ。

しかし,Maxwell & Delaney (1990) によれば,水準別誤差項での対比較は MANOVA と整合性がある,すなわち, 有意差の有無が一致する,らしい。(ふつうの ANOVA では整合しない。)

検定の多重性

単純効果の検定の場合も,検定の多重性が問題になるはずだが,これを考慮に入れて調整しているような研究はほとんどない (=editorにそれをしなくてもOKだと容認されている)。 多くの研究者は3つ以上の水準を取る要因の水準ペアの多重比較には敏感なくせに, こちらに敏感にならないのはなんとも奇妙で,ダブルスタンダードにすら聞こえる(しかしそれは大多数の実情である)。

そもそも,複数の要因効果の検討にも多重性の問題はあるはず。 しかしそれを明示的に調整しないのは,ANOVA全体で調整されているからではない。 下記シミュレーションを実行してみよう。

x1 <- gl(2, 12)
x2 <- gl(2, 6, 24)
ps <- replicate(1000, {
  y <- rnorm(24)
  anova(lm(y ~ x1 * x2))[[5]][-4]
})
apply(ps < 0.05, 1, sum)
sum(apply(ps < 0.05, 2, any))

調整しないことは,要因効果の検討がそれぞれで 1 family だという言い分で正当化できると思うかもしれない。 しかし,現場での実際としては,1つの実験計画やその集まりに単位があることが多い。 つまり,効果を発見をしたと主張するかどうか(≒公表するかどうか. 残念ながら)は,行った実験や調査にて1つ以上の有意差が得られたかどうかで決められている。 すると帰無仮説は「すべての効果がない」である。 多重対比較での躍起と両立して正当化するのは難しいだろう。

対比とその交互作用

Kirk (1982) によれば,単純効果の検定だけでなく,個々の対比の交互作用を調べるのがよい。

ここでいう対比とは,ある要因にp個の水準があり, 水準iにおける平均値を mu_i としたとき,

psi = sum_{i=1}^p c_i mu_i

で定義されるpsiである。ただし,

sum_{i=1}^p c_i = 0

を満たす。要するに,各水準の平均値の重みづけ和である。 重み c_i は上の条件を満たす限り任意に決めてよい(意味のあるものは限られるが)。 この psi の推定値

hat{psi} = sum_{i=1}^p c_i bar{y}_i

はサンプルから計算できるが,母集団において psi が 0 かどうか,つまり帰無仮説

H_0: psi = 0

を検定したいとする。例えば,帰無仮説の下では

{:
( E[hat{psi}] ,=, sum_{i=1}^p c_i E[bar{y}_i] ,=, 0 ),
( V[hat{psi}] ,=, sum_{i=1}^p c_i^2 V[bar{y}_i] ,=, sigma_e^2 sum_{i=1}^p frac{c_i^2}{n_i} )
:}

であるから, t 統計量を求めて

t = frac{sum_{i=1}^p c_i bar{y}_i}{sqrt{hat{sigma}_e^2 sum_{i=1}^p c_i^2 / n_i}} qquad sim t(n - p)

という具合に自由度 n - p の t分布を使った検定を行うことができる。ここで hat{sigma}_e^2 は分散分析でも得られる誤差分散の推定値

hat{sigma}_e^2 = frac{sum_{i=1}^p sum_{k=1}^{n_i} (y_{ik} - bar{y}_{i.})^2}{n - 1}

である。F分布を使うなら

F = frac{(sum_{i=1}^p c_i bar{y}_i)^2}{hat{sigma}_e^2 sum_{i=1}^p c_i^2 / n_i} qquad sim F(1, n - p)

さて,例えば c_i = 1, c_2 = -1 (残りはすべて 0)とすると, psi = mu_1 - mu_2 であり, 検定の帰無仮説は mu_1 - mu_2 = 0 ,すなわち mu_1 = mu_2 であるから, 水準1と水準2の平均値に差があるかどうかの検定と同じ意味になる。 従って明らかに,対比がカバーする範囲は各水準の対比較 (pairwise comparison) の範囲を完全に含む

また,例えば c_1 = 1/2, c_2 = 1/2, c_3 = -1 (残りはすべて 0)とすると, 水準1,2の平均と水準3との比較,といったことも可能になる。 対比のバリエーションはかなり豊富にある(というか,ストレートに解釈できないものも含めれば無限にある)。

ちなみに,計画的対比 (planned contrasts) でない場合に,多重検定を考慮して, 上の F 統計量に対してもっと保守的な臨界値

F ge (p - 1) F(p - 1, n - p, alpha)

を使って同時検定を行う(あるいは同時信頼区間を求める)やり方は,Scheffé の方法と呼ばれている。

処理-対比交互作用

処理-対比 交互作用 (treatment-contrast interactions) とは,例えば2要因計画の場合,ある要因Bの処理水準によって, もう一つの要因Aについての対比が一定でない,という交互作用である。 これを検討することで,3水準以上ある場合の細かな交互作用パターンがわかる。

要因Aに関する対比 psi_A を要因Bの水準 B_j 毎に考えるとき,これを psi_{Aj} として,一般に処理-対比 交互作用の平方和は

"SS" _{hat{psi}_A xx B}
 = n { sum_{j=1}^q (hat{psi}_{Aj} - {sum_{j=1}^q hat{psi}_{Aj}} / q)^2 } / {sum_{i=1}^p c_i^2}
 = n { sum_{j=1}^q hat{psi}_{Aj}^2 - 1/q (sum_{j=1}^q hat{psi}_{Aj})^2 } / {sum_{i=1}^p c_i^2}

ここで p は要因Aの水準数,q は要因Bの水準数, n はセル内の繰り返しの数である。自由度は q - 1 である。

例えば,3水準ある要因Aについて psi_A = 1 xx mu_1. + (-1) xx mu_2. + 0 xx mu_3. という対比を考えるとしよう。そして,要因Bの各水準 B_j (j = 1, 2, 3) の間で,この要因Aに関する対比 psi_{Aj} に違いがあるかを検討する。 平方和は

{:
( "SS" _{hat{psi}_A xx B} ,=, n { sum_{j=1}^3 hat{psi}_{Aj}^2 - frac{1}{3}(sum_{j=1}^3 hat{psi}_{Aj})^2 } / {sum_{i=1}^3 c_i^2} ),
( ,=, n/2 (sum_{j=1}^3 (bar{y}_{1j} - bar{y}_{2j}) ^2 - 1/3 (sum_{j=1}^3 bar{y}_{1j} - bar{y}_{2j})^2) )
:}

となる。これと誤差の平方和とでF比をつくって検定すればよい。 ただしここでも,上記の単純効果の場合と同様,要因Bの一部の水準について検討する場合は誤差項をプールするかどうかの選択がある。

対比-対比交互作用

処理-対比 交互作用の検定で有意な結果を得たなら,さらに詳細な検討をする。 対比-対比 交互作用 (contrast-contrast interactions) とは,例えば2要因計画の場合, ある要因Aについての対比によってもう一つの要因Bについての対比が一定でない,という交互作用である。 つまり,処理-対比 交互作用では片方の要因の水準毎に調べていたのを,両方とも対比にして検討するということ。 Kirk (1982) の3つの中では,これのみが純粋に交互作用を検討する方法である (単純主効果や処理-対比交互作用は主効果と confound している)。

要因Bの水準 B_j (j=1,...,q) 毎の要因Aに関する対比を psi_{Aj} , 要因Aの水準 A_i (i=1,...,p) 毎の要因Bに関する対比を psi_{Bi} とすると 平方和は

"SS" _{hat{psi}_A xx hat{psi}_B}
 = n ( sum_{j=1}^q c_j hat{psi}_{Aj} )^2 / {sum_{i=1}^p c_i^2 sum_{j=1}^q c_j^2}
 = n ( sum_{i=1}^p c_i hat{psi}_{Bi} )^2 / {sum_{i=1}^p c_i^2 sum_{j=1}^q c_j^2}

検定の仕方は処理-対比交互作用と同様。自由度は 1 である。

一般線型仮説

固定効果のみの分散分析を含め,(単変量)線型モデル

y = X beta + varepsilon

における beta に関する帰無仮説は,一般的に

H_0: C beta = h

と表せる。これを線型仮説 ( linear hypothesis ; 一般線型仮説 general linear hypothesis , 線形仮説 ) と呼ぶ。 ここで,C は k × q 行列(仮説行列と呼ばれることもある), beta は q × 1 の母数ベクトル, h は k × 1 の定数ベクトルである。 例えば,

X = ((1,1,0),(,vdots,),(1,0,1),(,vdots,),(1,-1,-1),(,vdots,)), quad beta = ((beta_1),(beta_2),(beta_3))

のとき,

H_0: mu_1 = mu_2 = mu_3

という,三群の平均値が等しいという仮説は

H_0: beta_2 = beta_3 = 0

と同じであり,これは

H_0: ((0,1,0),(0,0,1))((beta_1),(beta_2),(beta_3)) = ((0),(0))

と書ける。beta の推定値 b は, W = (X^TX)^-1 とおくと

b = W X^T y qquad sim N(beta, sigma^2 W)

帰無仮説のもとでの制約付き推定値 b_H は, V = ( C W C^T )^-1U = W - W C^T V C WD = C b - h とおくと

{:
( b_H ,=, U X^T y + W C^T V h ),
(     ,=, b - W C^T V D )
:}

線型仮説を検定するための統計量は一般に

F = frac{ "SS" _H // nu_H}{ "SS" _E // nu_E} qquad sim F_(nu_H,nu_E)

が使える。ここで,"SS" _E"SS" _Hnu_Hnu_E はそれぞれ,

{:
( "SS" _E ,=, || y - hat y || ^2 ),
( ,=, (y - hat y)^T (y - hat y) ),
( ,=, (y - hat y)^T y ),
( ,=, (y^T - b^T X^T) y ),
( ,=, y^T ( I - X W X^T ) y qquad sim sigma^2 chi_{nu_E}^2 )
:}
{:
( "SS" _H ,=, || y - hat y_H ||^2 - || y - hat y ||^2 ),
( ,=, || y - X (b - W C^T V D) ||^2 - || y - hat y ||^2 ),
( ,=, ( y^T - y^T X W X^T + D^T V C W X^T )( y - X W X^T y + X W C^T V D ) - || y - hat y ||^2 ),
( ,=, y^T ( I - X W X^T ) y + D^T V D - || y - hat y ||^2 ),
( ,=, D^T V D  qquad sim sigma^2 chi_{nu_H}^2 ),
( ,=, || X W C^T V D ||^2 = || hat y - hat y_H ||^2 )
:}
nu_H = "rank" (C) = k
nu_E = n - q

反復測定

球形仮定

上記の通り,古典的な分散分析では,要因が被験者間だとか被験者内だとか,集団で実験しているとかそうでないとか,そんなことにはお構いなく, すべての観測に関して誤差分布は独立であることを前提として理論が組み立てられている。 例えば,上節の分割区画デザインなどで被験者内要因を組み込んだモデルでは,被験者間要因の各水準における分散共分散行列 Sigma_m を考えると,独立な誤差の他に「被験者」因子による共分散が入って,

Sigma_m = (
(sigma_e^2+sigma_s^2, sigma_s^2, ..., sigma_s^2),
(sigma_s^2, sigma_e^2+sigma_s^2, , sigma_s^2),
(vdots, , ddots, vdots),
(sigma_s^2, sigma_s^2, ..., sigma_e^2+sigma_s^2)
)

と仮定される。分散共分散行列に対するこのような条件は複合対称性 (compound symmetry) と呼ばれる。 Huynh & Feldt (1970) はこの条件を満たす行列を S 型行列 (type S matrix) と呼んでいる。

ところが,実際に被験者内要因の入った実験や調査を行ったときにこのように整った分散共分散行列が得られるとは限らない― 例えば共分散が抽出誤差程度でないほどバラバラになる。 するとそのような場合,前提が満たされていないのだから,分散分析という方法の適用はまったく誤りであるかのように思われるかもしれない。 しかし,そうではない。

古典的な分散分析は最終的に統計量 F を用いて検定を行うのであるが,そのロジックは, 上の節で述べたような諸仮定(誤差の等分散など)が満たされているならば所定の計算で得られる統計量 F が帰無仮説のもとで F 分布と呼ばれる分布に従いますよ,ということである。 この逆は述べられていない。つまり,諸仮定が満たされていることは十分条件であり, 諸仮定を満たしていなくても F 分布に従う余地があるのだ。

では,反復測定を入れている場合,帰無仮説のもとで F 分布に従うための必要十分条件とは何か。 それは,循環性仮定 (circularity assumption)と呼ばれており,共分散行列 Sigma_m が次の形をとることである:

Sigma_m = (
(sigma_1^2, frac{sigma_1^2+sigma_2^2}{2}-lambda, ..., frac{sigma_1^2+sigma_p^2}{2}-lambda),
(frac{sigma_1^2+sigma_2^2}{2}-lambda, sigma_2^2, , frac{sigma_2^2+sigma_p^2}{2}-lambda),
(vdots, , ddots, vdots),
(frac{sigma_1^2+sigma_p^2}{2}-lambda, frac{sigma_2^2+sigma_p^2}{2}-lambda, ..., sigma_p^2)
) = (
(2 alpha_1 + lambda, alpha_1 + alpha_2, ..., alpha_1 + alpha_p),
(alpha_1 + alpha_2, 2 alpha_2 + lambda, , alpha_2 + alpha_p),
(vdots, , ddots, vdots),
(alpha_1 + alpha_p, alpha_2 + alpha_p, ..., 2 alpha_p + lambda)
) = lambda I_p + alpha 1_p^T + 1_p alpha^T

ここで alpha = ((alpha_1, alpha_2, ..., alpha_p))^T である。 この形の行列は H 型行列 (type H matrix) と呼ばれたり,Huynh & Feldt の名をとって H-F 共分散構造と呼ばれたりする。

この仮定は,共分散行列 Sigma_m を考えている p 個の反復測度 m_i のうち2つの差の分散が,どの組をとっても一定の値 2 lambda であるということを意味する:

V[m_i - m_j] = sigma_i^2 + sigma_j^2 - 2 sigma_{ij} = 2 lambda qquad qquad i != j

これからさらに次の事が導かれる: k × p ( k < p ) の正規直交対比行列 (orthonormal contrast matrix) C (すなわち,C C^T = I_kC 1_p = 0_k) と正定数 lambda でもって

V[Cm] = C Sigma_m C^T = C (lambda I_p + alpha 1_p^T + 1_p alpha^T) C^T = lambda I_k

これは,共分散行列の直交分解が球形 (spherical),すなわち,どの方向にも同じ大きさになる,ということを意味している。 よってこの仮定は,球形仮定 (sphericity assumption) と呼ばれている。 あるいは,直交変換によって等分散になる,と言った方がよいだろうか。 球形仮定は,球面性仮定とか球状性仮定とも呼ばれる。等方性とも呼ばれる。 以下では千野(1993)に習って球形仮定と呼ぶことにする。

上の式からわかるように,反復測度が2水準しかなければ,球形仮定は常に成り立つ。 よってその場合,球形の確認のための検定や自由度の修正などの複雑な議論はまったく関係ない。

また,上節での乱塊法の説明にて,交互作用を含めたモデルとそうでないモデルを紹介しているが, 球形仮定が成り立つならばブロックとの交互作用はない。 つまり,(大局的)球形仮定が満たされていると考えられる場合ならば,交互作用なしモデルでよい,ということだ。

ところで,下節で紹介する一般化多変量分散分析は,共分散行列 Sigma_m の構造にまったく仮定を置かない。 つまり,分散共分散行列を p×p 個の自由パラメタとして扱う。 よって,このような球形仮定の話とは無縁になれる。サンプルサイズが十分な場合には,そちらがお勧め。

さらには,共分散行列が球形仮定を満たしていなくても,その他の構造を満たしているなら,それを適切にモデリングした方がよい。 線型混合モデルを使おう。

自由度の修正

球形仮定が成り立っていれば,F 比は F 分布に従う。成り立っていない場合は一般に,F 検定の第一種過誤の確率はインフレする。 インフレを修正するために Box (1956) は次のように定義される自由度の修正係数 epsi を導入した:

epsi = frac{["tr" (C Sigma_m C^T)]^2}{k "tr" [(C Sigma_m C^T)^2]} = frac{(sum_{i=1}^{k} lambda_i)^2}{k sum_{i=1}^{k} lambda_i^2}

ここで lambda_iC Sigma_m C^T の固有値である。 この定義から,epsi1//k le epsi le 1 の範囲の値となる。「epsi の下限」と呼ばれるのは,この左辺の 1//k のことである。 この epsi を F 値の分母,分子の自由度にかけて p 値を出すことで,第一種過誤のインフレを防げる。

上の定義式に母共分散行列 Sigma_m が入っていることから, 実践的には epsi の値は未知のため,推定しなければならない。 そこで Greenhouse & Geisser (1959) が推定値を提案した。 この G-G の推定値はサンプルサイズが小さいときに保守的すぎるとして,Huynh & Feldt (1976) が G-G の改善案を提案した。 しかし,H-F は甘すぎるとの指摘もある(Maxwell & Delaney, 1990)。 現在,心理学界などで広く用いられているのは,G-G の推定値のようである。 Girden (1992) は G-G の推定値が 0.75 より大きいときは H-F を使用し,そうでないときは G-G を使用するよう勧めている。

球形仮定が成り立っているかどうかを判断するため統計ソフトでは Mauchly の検定が自動的に出力されたりするが, そのような球形仮説の検定の頑健性が低いとの観点から, 球形仮説の検定をせずに Greenhouse & Geisser (1959) の3ステップ手続き(つまり,最初に通常の,次に下限,最後にG-GまたはH-F)を勧める論者もいるようだ。

局所的球形仮定

上記の球形仮定の定式化にて,対比行列の行数を k < p としているように,球形仮定は k = p - 1 以外の場合にも問われる。 k = p - 1 の場合,すなわち,反復測定として計画される全要因全水準を含めて考える場合,大局的球形仮定(global sphericity assumption)と呼ばれる。 一方,k < p - 1 の場合,すなわち一部の対比に限定して考える場合,局所的球形仮定(local sphericity assumption)と呼ばれる。 これらはそれぞれ,全体的循環性(overall circularity)と局所的循環性(local circularity)と呼ばれることもある。

わざわざ一部の対比を考えるということは,つまり,全体的な球形仮定が満たされていない場合でも, 主効果や交互作用効果のそれぞれの検討において局所的な球形仮定が満たされていたりいなかったりが併存し得るということだ。 よって,自由度の調整で球形仮定の違反をしのごうとするならば,どのように調整するか(あるいはしないか)は個々の要因効果ごとに検討すべきということになる。

局所的球形仮定は,特に,上節で説明している単純効果を調べる場合に重要になる。 水準別誤差項の考え方では,誤差分布を別々に推定するのだから,それぞれにおいて球形仮定が満たされている共分散構造なのかどうかをチェックすべきということになる。

多標本球形仮定

デザインが個体内の反復測定要因だけでなく個体間要因も含んでいる場合,個体間要因の各条件を仮に別々の母集団と考え, それぞれの標本を入手しているということで,多標本だと見なす。 つまりは,上節の分割区画デザインのような場合のことである。

その場合,F 比が F 分布に従うのに必要な仮定としては,個体間要因の各条件において球形仮定が満たされていることだけでなく, 条件間でその球形な共分散行列が同一でなければならない。 これを多標本球形仮定(multisample sphericity assumption)あるいは多標本循環性(multisample circularity)と呼ぶ。

この条件間の等質性の仮定は,多変量分散分析でも同様にあるので,多変量分散分析を使ったとしても避けられない。 避けたければ,ノンパラか,他の heteroscedasticity を許すよう一般化されたパラメトリックモデルに持ち込もう。

事後分析

独自の統計量を使うような多重比較用に開発された手法のほとんどは,i.i.d.を仮定している。 よって,球形仮定が成り立たない場合に使ってはならないことはもちろん,複合対称性の下でしか使ってはならないものもある。(Boik, 1981)

有意水準を調整するタイプのうち Bonferroni 法が頑健でまだマシだと言われているので (Keselman & Keselman, 1988) , どうしてもというならこれを使って t 検定などする。 しかしこれでも,異分散やサイズが不釣り合いなら場合によっては第一種過誤の確率が水準をオーバーする。

第一種過誤にうるさい人は,かなり保守的であるが多変量の同時信頼区間を使うのがよい。

もっと一般化されたモデル

「一般化された」は「制約を外された」と換言できる。 制約(必要な前提)がより少ない分析方法が使えるのはすばらしいことだ。 古い方法はますます出番が無くなる。分散分析もこの運命をたどるように見える。 実践からは姿を消して学習教材としてしか残らないかもしれない。

「分散分析モデルは一般線型モデルに含まれる」と述べる書籍やウェブページがあるが, 一般線型モデルという呼び名は通常,固定効果のみを扱うモデルを指しているようだ。 分散分析モデルは元来,固定効果,変量効果の両方を扱うモデルであるから,従って,そのような意味ならば 「一般線型モデルに含まれる」というのは正しくないと思われる。 歴史上複数のルーツで開発され別々に扱われていた統計解析手法を統合した,という表現もしばしば一般線型モデルに対して使われているのを見かけるが, 回帰分析,分散分析,共分散分析の完全な統合を行ったのは混合効果モデルだと言うほうが整合性があるだろう。

すでにさらなる一般化がなされた解析法(一般化線型混合モデル,構造方程式モデルなど) が実践で活躍しているが(R でも実行可能), これらもパラメタに関する線型性を含意していることを忘れてはいけない。 そして,パラメタに関して非線型なモデルでしか適切に捉えられない現象も自然には多々ある。 心理学者たるあなたが今調べようとしている現象がまさにその種の現象かもしれないのである。 そのような場合は非線型回帰モデルなどにも目を向けるべし。 一般化加法モデル (generalized additive model) もすでに実践で活躍するようになっている(これも R で実行可能)。

多変量分散分析

多変量分散分析 (multivariate analysis of variance; MANOVA; multiple analysis of varianceと書いている文献もあるが紛らわしい) とは, 複数の従属変数を同時にモデルに組み込む分散分析であり,単変量の分散分析の拡張であると見なせる。 単変量の分散分析では各条件の平均値を比較するが,多変量分散分析では条件間で平均値ベクトルを比較する。

多変量分散分析を行う目的は大きく3つ挙げられる(時にはこれらのうちの複数を目的としているかもしれない)。

  1. 実験毎 (experimentwise) での第一種の過誤 (type I error) の確率を一定の水準以下に保つ
  2. 従属変数間の相関も考慮に入れる
  3. 単変量では検出できないかもしれない独立変数の効果を,複数の従属変数を組み合わせることで検出する

1つ目の目的は,いわゆる protected の理屈である。 ただし,多変量分散分析の後に従属変数ごとに単変量の分散分析を行うというやり方には批判がある。 これはフォローアップの節を参照。

2つ目の目的に関して,単変量の分散分析では変動要因ごとに平方和(分散)を分解するが, 多変量分散分析では平方和積和行列(分散共分散行列)を分解する。 つまり分散だけでなく共分散(すなわち相関情報)も分析対象にしているということである。

3つ目の目的に関しては,多変量分散分析は複数の従属変数の線型合成(重み付け和)を作成し,それに対する独立変数の効果を分析する。 しかも独立変数の効果(群間の違い)が最大になるように係数を決める。

多変量分散分析のモデルと仮定

多変量分散分析モデルは,

Y = X ccB + ccE

ここで,Y はデータ行列(n×p; 被説明変数), X は被験者間の計画行列(n×k), ccB は母数行列(k×p), ccE は誤差行列(n×p)。

さらに,一般化多変量分散分析モデル (generalized multivariate analysis of variance; GMANOVA) は,

Y = X ccB W + ccE

ここで,W は被験者内の計画行列(q×p; q≦p), また ccB は母数行列(k×q)である。

たとえば,被験者間2×2要因計画(交互作用なし)で2つの従属変数を持つ場合,次のようになる。

((y_11,y_12),(y_21,y_22),(vdots,),(y_(n1),y_(n2))) = 
((1,1,1),(,vdots,),
 (1,1,0),(,vdots,),
 (1,0,1),(,vdots,),
 (1,0,0),(,vdots,))
((beta_11,beta_12),(beta_21,beta_22),(beta_31,beta_32)) + 
((varepsilon_11,varepsilon_12),(varepsilon_21,varepsilon_22),(vdots,),(varepsilon_(n1),varepsilon_(n2)))

被験者内のみの2×2要因計画(交互作用あり)ならば,GMANOVAモデルを用いて,

((y_11,y_12,y_13,y_14),(y_21,y_22,y_23,y_24),(,vdots,,),(y_(n1),y_(n2),y_(n3),y_(n4))) = 
((1),(vdots),(1))
((beta_11,beta_12,beta_13,beta_14))
((1,1,1,1),(0,1,0,1),(0,0,1,1),(0,0,0,1)) + 
((varepsilon_11,varepsilon_12,varepsilon_13,varepsilon_14),
 (varepsilon_21,varepsilon_22,varepsilon_23,varepsilon_24),
 (,vdots,,),
 (varepsilon_(n1),varepsilon_(n2),varepsilon_(n3),varepsilon_(n4)))

GMANOVAモデルは,成長曲線モデル (growth curve model) とも呼ばれる。というのも,計画行列の作り方次第で,

((y_11,y_12,y_13),(y_21,y_22,y_23),(,vdots,),(y_(n1),y_(n2),y_(n3))) = 
((1,1),(vdots,),
 (1,0),(vdots,))
((beta_11,beta_12,beta_13),(beta_21,beta_22,beta_23))
((1,1,1),(1,2,3),(1,4,9)) + 
((varepsilon_11,varepsilon_12,varepsilon_13),(varepsilon_21,varepsilon_22,varepsilon_23),(,vdots,),(varepsilon_(n1),varepsilon_(n2),varepsilon_(n3)))

こんなこともできるからだ。これは3つの従属変数を3時点での反復測定と考えたときに2条件それぞれにおいて2次曲線をあてはめるモデルになっている。

多変量分散分析においては誤差成分 ccE に関して,多変量正規性,共分散行列の等質性,独立性の3つの仮定があり, これらはそれぞれ,上のほうで述べている単変量の分散分析における正規性,等分散性,独立性の仮定の拡張版だと考えればよい。

多変量正規性

上記の誤差行列 ccE については,各行ベクトルが期待値 0 の多変量正規分布に従うと仮定される:

varepsilon_(i.) = (varepsilon_(i1),varepsilon_(i2),...,varepsilon_(ip)) sim N_p(0, itSigma)

これは,従属変数ごとに誤差がそれぞれ正規分布している,ということではない。もっと大変。 多変量正規分布の密度関数は,変量の数を p としたとき,

frac{1}{sqrt{(2pi)^p|itSigma|}}exp{-frac{1}{2}(x - mu)^T itSigma^{-1} (x - mu)}

である。ここで,mu は p 個の平均値ベクトル,itSigma は p × p の分散共分散行列である。

多変量正規分布に従う確率変数ベクトルは,どんな(0でない)線型変換を施しても多変量正規分布する。 これは,p個の確率変数のうち任意の一部を取り出しても,多変量正規分布することを含意する。 つまり,一部というのが「1つ」の場合として,すべての確率変数がそれぞれ個別に見ると正規分布することも含意する。 しかも周辺分布と一致する。 さらに,そのベクトルの要素のいかなる線型結合も正規分布する,条件付き分布も正規分布する,無相関なら独立である,など,色々な性質が伴う。 まさに理想的な分布である。

単変量の分散分析と同様に多変量分散分析も多変量正規性の崩れには頑健であると一般に信じられている。 しかし,多変量分散分析の頑健性についての研究はまだ十分積み重ねられているわけではないようだ。 あまり盲信しない方がいい。

共分散行列の等質性

単変量の分散分析で正規性に加えて分散が等しいことを仮定したように,多変量分散分析でも上の多変量正規性に加えて, 条件内で等しいのはもちろんのこと条件間で共分散行列が等しいことを仮定する。 たまにある間違い:これは共分散行列 itSigma の成分である複数の分散や共分散が等しいという意味ではない。 上の表記で言えば,各 varepsilon_(i.) がすべて同一の多変量正規分布に従うということである。

この仮定をテストするために,通常,Box の M や Bartlettの chi^2 という指標が用いられる。 が,これらの検定を行う前に,まず多変量正規性を確認すべし。

共分散行列の等質性の仮定を満たさない場合,各条件の n がほぼ等しいなら,多少検出力が落ちるくらいで大きな影響はないと言われている(Stevens, 1986)。 もし n が大きく異なるなら,第一種過誤の確率が大きく変わる(上にも下にも)ので気をつけよ。

独立性

単変量の場合と同じく,この仮定が一番重要。 従属変数の個々の値は他の値とは独立に(互いに影響されずに)得られたものであると仮定される。 上の表記で言えば,各 varepsilon_(i.) が独立に分布するということである。

しかし例えば,集団で実験をやっていて,Aさんの反応はBさんの反応にひきずられるかもしれない。 こういうのは級内相関を用いてテストできる。だが,ちょっと相関があっただけで第一種の過誤の確率は何倍にもインフレする。 独立性の仮定が満たされないのは,まずもって実験手法の問題なことがほとんどなので,くれぐれも気をつけるべし。

多変量一般線型仮説

上で説明している一般線型仮説は,多変量の場合にも一般化できる。

H_0: C ccB A = G

ここで,C は k × q 行列 (k ≦ q), ccB は q × p の母数行列, A は p × r 行列 (r ≦ p), G は k × r の定数行列である。

ccB の推定値 B は, S = Y^T ( I - X (X^T X)^-1 X^T ) Y とおくと,

B = (X^T X)^-1 X^T Y S^-1 W^T (W S^-1 W^T)^-1

誤差と仮説の積和行列にあたる EH は,

{:
( E ,=, A^T ( W (Y^T ( I - X (X^T X)^-1 X^T ) Y)^-1 W^T )^-1 A ),
( H ,=, (C B A - G)^T ( C ( (X^T X)^-1 + (X^T X)^-1 X^T Y ( S^-1 - S^-1 W^T (W S^-1 W^T)^-1 W S^-1 ) Y^T X (X^T X)^-1 ) C^T )^-1 (C B A - G) )
:}

となる。W = I の場合,すなわち,被験者内要因のない MANOVA の場合は,

{:
( B ,=, (X^T X)^-1 X^T Y ),
( E ,=, A^T Y^T ( I - X (X^T X)^-1 X^T ) Y A ),
( H ,=, (C B A - G)^T ( C (X^T X)^-1 C^T )^-1 (C B A - G) )
:}

と単変量の場合に似た単純な形になる。

多変量検定

上記のとおり,多変量分散分析では平均値ベクトルについて仮説をテストするために 要因効果の平方和積和行列と誤差の平方和積和行列を比較したい。 しかし,複数の行や列を持つ行列だと,単変量の場合のようにそのまま比を計算できない。 このため,平方和積和行列の何らかの「大きさ」の比較をする統計量がほしい。 そして,どのように相対的「大きさ」を測るかで複数のやり方が提案されている。 大抵の統計ソフトでは,下の4つが出力できるようになっている。 以下 lambda_i は行列 H E^-1 の固有値である。

Wilks' lambda

文献でもっともよく見かけるもの。 Lambda = frac{|E|}{|H+E|} = frac{1}{|I+HE^-1|} = prod (1 + lambda_i)^-1 。 そのままの分布は複雑なので,変換して,特別な場合に F 分布に従うことや,chi^2 分布に近似することを利用して検定する。 意味としては単変量での 1 - R^2 っぽい。

Pillai's trace

V = "tr" ( H (H + E)^-1 ) = "tr" (( I + (H E^-1)^-1)^-1) = sum frac{lambda_i}{1+lambda_i} 。 Wilks' lambda より保守的であるが,頑健であるため,こちらの使用が推奨されている (Hand and Taylor, 1987; Olson, 1976)。R のデフォルト。 R^2 っぽい。

Hotelling-Lawley trace

U = "tr" (H E^-1) = sum lambda_i 。Wilks' lambda よりもリベラル。 たとえるなら F 比。

Roy's greatest root

最大の固有値 lambda_1 で定義。リベラルすぎて実践的には使わない。p値の下限値を知るのに利用。

上から3つの統計量はすべての正準相関(の二乗である固有値)を加味するが,最後のは正準相関のうち絶対値が最大のものだけを使う。

被験者間要因が2水準の時には,これらは違いはないのでどれでもよい。

多変量分散分析のフォローアップ

多変量効果が有意だった場合に,さてどうするか。多変量分散分析の難しいところでもある。

従属変数毎に単変量の分散分析

これは,よく見かけるが,あまりおすすめしない。

最初に多変量分散分析をやれば experimentwise の過誤の確率を適切に抑えられると期待するかもしれないが, protected の理屈が通るのは帰無仮説が真のときである。 研究上のたいていのケースでは,いくつかの従属変数では差があり他の従属変数では差がない,というようなことになっている。 つまり,多変量有意性検定における帰無仮説は偽であり,十分な検出力があれば多変量分散分析が帰無仮説を棄却する状況である。 このような場合,多変量分散分析はその後に普通にα=.05で単変量分散分析を繰り返すことへの protection になっておらず, 1つ以上の差がない従属変数について有意差があるという結果を得る可能性は,単変量分散分析を繰り返せば繰り返すほど増える。 従って,Bonferroni の修正などを行わずに単変量の分散分析を従属変数ごとに適用する仕方は,やってはいけない

また,個々の単変量の分散分析は,従属変数間の相関を考慮しないので,もったいない。 以下に紹介する他の方法を使うことをおすすめする。

重判別分析

多変量分散分析の独立変数を目的変数,従属変数を説明変数として(つまり逆にして),重判別分析を行う。 どの「従属」変数が寄与しているか(あるいは,条件間を異ならせている基本的な次元)を探す。

ステップダウン分析

多変量分散分析では従属変数間の相関が許されるが,この相関について,従属変数間の因果関係が予め想定される場合, その因果順序に従って分析に組み込んでいく仕方。

まず,因果連鎖において最初に登場する変数 A だけで単変量の分散分析を行う。 次に,2番目に登場する変数 B を分析するときには,最初の変数 A を共変量とした共分散分析を行う。 その次に,因果連鎖の3番目の変数 C を分析するときには,その前の変数 A,B を共変量とした共分散分析を行う。 以下同様。

Wilkinson (1975) の方法

従属変数のうちのどれか1つを除いた多変量分散分析をそれぞれ行い,F値の変化を調べ, 多変量での有意性にどの変数が寄与しているかを探るというやり方。

Huberty & Smith (1982) も似た考え方を提唱している。

多変量対比

上の手法たちは従属変数を見分ける目的のものだが,多変量対比は独立変数の水準を見分けるもの。 単変量のANOVAと同様,対比によってどの群(の組み合わせ)の間に違いがあるのかを検討する。 2群を比較する際の多変量検定統計量として有名なのは Hotelling の T^2 である。 一般には,上記の多変量線型仮説の検討である。

有意な T^2 を得たら,これは多変量の比較であるから, さらにそのフォローアップを行いたいことが多いだろう。それには,上記の各手法を用いればよい。

一般化線型モデル

上節のとおり,線型モデル(固定効果の分散分析を含む)は一般に

y = X beta + varepsilon

と定式化される。期待値の観点からは

E[y] = X beta

と表せる。もちろん分布は正規分布である。これを

g(E[y]) = X beta

と拡張する。つまり,右辺の線型な式と y の期待値の間に関数 g を一段かませることを許す。 さらに,分布についても正規分布に限らず指数分布族(正規分布,二項分布,ポアソン分布など)であればよい,と拡張する。 これが 一般化線型モデル ( generalized linear model; GLM ) と呼ばれる統計モデルである。

関数 g はリンク関数 (link function; 連結関数) と呼ばれる。関数であれば何でもよいというわけには行かず,逆関数があることだけでなく, 統計学的によい性質を持っていることが望ましい。 よって実践では多くの場合,正準リンク関数 (canonical link function) を使う。よく使われる分布と正準リンク関数の組は次のとおり:

分布リンク関数逆関数
正規分布 ( gaussian )identityg(mu) = mug^-1(x) = x
二項分布 ( binomial )logitg(mu) = log(frac{mu}{1 - mu})g^-1(x) = frac{1}{1 + e^{-x}}
ポアソン分布 ( poisson )logg(mu) = log(mu)g^-1(x) = e^x
ガンマ分布 ( Gamma )inverseg(mu) = frac{1}{mu}g^-1(x) = frac{1}{x}
ワルド分布 ( inverse.gaussian )1/mu^2g(mu) = frac{1}{mu^2}g^-1(x) = frac{1}{sqrt(x)}

ここに挙げた名前は R での分布とリンク関数の指定に使える名前である。用意されている分布とリンク関数は他にもたくさんある上に,R ではリンク関数をユーザが作ることもできるのですばらしい。

二項分布でリンク関数に logit を使うモデルは伝統的にロジスティック回帰と呼ばれてきた。 ポアソン分布に対してリンク関数 log を使うモデルは伝統的にポアソン回帰あるいは対数線型モデルと呼ばれてきた。 つまり GLM はこれらを統合している。

古典的な線型モデルは,一般化線型モデルで gaussian(link="identity") を使った場合に相当する。よって,固定効果の分散分析は glm() でできる。 さらには,他の分布とリンク関数を試してどちらが当てはまりがよいか比較検討することもできるから,一般化線型モデルのほうが視野が広がる(=真実に近づく)。 よってこれからの時代は,最初から分散分析と言わず GLM なスタンスでいたほうがよい。 (が,実際に R で常時 glm() を使え,ということではない。lm() で足りるモデルなら lm() を使う方がよい。)

データ変換

分散分析を適用する前にデータ変換 (data transformation; あるいは変数変換 transformation of variables) を行って,その結果を分散分析する,という手続きがとられることがある。 この操作は,古典的な教科書には大抵載っている。 例えば,反応時間のデータが対数正規分布すると考え,収集したデータに対して対数変換を行い,その結果を分散分析にかける,など。

実践的には,データ変換は大きく分けて次の3つの目的のために行われる:

  1. 適切な解釈をしやすくする
    1. もとのデータとは異なるスケール(対数,逆数など)での単純な関係を見出す
    2. データをあるスケールへ変換し,平均値について線型モデルが当てはまるようにする
  2. 統計的推測の改善
    1. 等分散性 homoscedasticity を得る
    2. 正規性 normality を得る
  3. 視覚化のための便宜
    1. プロットが特定の誌面サイズに収まるようにする

ふつう,これらの目的は下へ行くほど重要性が低くなる。

例えば,二項分布する比率データに対して逆正弦変換(角変換)を適用しようとするかもしれない。それによって,分散を一定にし,正規分布に近づける。 しかしこの操作は,ロジスティック回帰を用いた場合のようなパラメタの解釈の容易さを無くする。 あるいは,例えばポアソン分布するカウントデータに平方根変換を適用して分散を一定にしようとするかもしれない。 しかしこれも,ポアソン回帰のような率直な解釈ができなくなる。 つまりは,2. の目的のために 1. を失っていて本末転倒なのだ。

一般化線型モデルは,言わば,1. の目的と 2. の目的を分けて達成する。よって,2. のために 1. を犠牲にしない。 これは実践的には非常に重要である。特に,もとの(変換前の)スケールで結果を解釈することが求められる研究ではなおさらだ。 これが,一般化線型モデルが使えるデータでは積極的に使うべき,な理由の1つである。

解釈を難しくすることの極端な例として,データ変換がイリュージョンを見せる場合すらある。 例えば,データ変換によって,本当は存在しない交互作用が有意だとされ (粕谷, 2004),交互作用があると勘違いしてしまう。 信じられない人は,ポアソン分布のデータに対して平方根変換や対数変換をした場合の R でのデモを伊東さんが書いてらっしゃるので試してみよう。 しばしば見かける「変数変換+分散分析を使い,あとは解釈時の注意でなんとかせよ」という指南は,この目の前に現れた交互作用の星印をユーザの解釈の力だけで打ち消してしまえ,というかなりハードな課題を与えているのである。

このようなアーティファクトが発生するのは,まさに,分析手法の前提条件を満たすためのデータ変換によって目的 2. はよろしくとも目的 1. が犠牲になる,すなわち平均構造が歪められることに起因する。

例えば,よく言われるように,対数変換で得られたデータの算術平均は,もとのスケールでは幾何平均であり,標準偏差も標準偏差でなくなって,解釈において意味がかなり違ってくる。 もちろん,研究上,幾何平均について物を言いたい場合は,対数変換は良き友である。一概に対数変換がダメなわけではない。 つまるところ変換の良し悪しは研究者の意図次第なのだが,その意図と実際がズレないようにするのが,何にもまして重要な仕事である。 下手をすれば査読者も気づかないで,要因効果が存在する態で研究が進んでしまう。

R で言えば,対数変換+分散分析に相当する

glm(log(y) ~ x, gaussian)

と,リンク関数に log を用いた GLM

glm(y ~ x, gaussian(link="log"))

が,分布について違う(前者は対数変換したものが正規分布,後者はもとのまま正規分布)だけでなく, 期待値についても違った構造を指定するモデル(前者は対数変換後の期待値が x の一次式,後者はもとの期待値の対数が x の一次式)であることを理解する必要がある, たとえ GLM を使わずに分散分析だけで生活する人であっても。 なぜならその分散分析のモデルが何を言わんとするものかは正確に把握していなければならないからだ。

GLM にほれ込むあまり,逆に,データ変換があらゆる場面で要らなくなると思ってもいけない。 上記の通り,解釈の面でも場合によってはデータ変換は欠かせない操作であり, また,それ以外の面でもデータ変換は有用であるから,データ変換の勉強はしっかりしておくべし。

線型混合効果モデル

混合効果モデル(mixed effects model; 混合モデル mixed model)とは,固定効果とランダム効果(=変量効果)を両方含むモデルのことである(上記 固定効果と変量効果 の節を参照)。 混合分布モデル (mixture model) と混同しないように注意。 特に線型であることを強調する場合は線型混合効果モデル(linear mixed effects model; 線型混合モデル linear mixed model) とも呼ぶ。 というのは,裏返せば,非線型混合効果モデル (nonlinear mixed effects model) というのも使える,ということだ。

混合効果モデルは,

という,これまた現場としては非常に嬉しい特徴がある。

線型混合効果モデルの一般式は

y = X beta + Z gamma + varepsilon

ここで,X は固定効果の計画行列(n×p),beta は固定効果ベクトル(p×1), Z はランダム効果の計画行列(n×q),gamma はランダム効果ベクトル(q×1)。 ランダム効果,誤差ともに正規分布に従うと仮定されており,また,期待値と分散を

{:
( E[(gamma), (varepsilon)] ,=, [(0),(0)] ),
( V[(gamma), (varepsilon)] ,=, [(G, 0), (0, R)] )
:}

と仮定している。 Gはランダム効果間の分散共分散行列(q×q),Rは誤差間の分散共分散行列(n×n)である。 つまり,ランダム効果間や誤差間の共分散が 0 でないモデリングをすることができる。 一方,固定効果との間はもちろんのこと,ランダム効果と誤差の間は無相関だと仮定している。 また,Rの対角要素はすべて同一である必要が無い。 すなわち,古典的な分散分析(を含む一般線型モデル)にあった,誤差の等分散性の仮定,独立性の仮定,を外すことができる。 被説明変数 y の分散は

V[y] = Z G Z^{T} + R

このようになり,誤差の分散共分散 R だけでなく ランダム効果の分散共分散 G にも規定される。 分析の実践においては,この RG の構造を指定することができる。

混合効果モデルから見れば,(古典的分散分析を含む)一般線型モデルは,R = sigma^2I , かつ, Z = 0 という,混合効果モデルの特殊ケースである。

この混合効果モデルを用いれば,上節の変量効果のところで書いた「被験者」因子のみならず, これまでランダム効果であるにもかかわらず固定効果として分析されていたものを適切にランダム効果と見なして分析することができる。 また欠損がある場合,通常の分散分析では扱い切れなかった (セルサイズが異なるときになんとか対策を練ってきたが) ものを listwise-deletion せずに正面切って扱うことができる(MCAR, MARの場合)。さらに,同一被験体のデータは相関することがよくあるが, それを無視して(ある程度ずれ過ぎていなればよしとして)分析してきた。 混合効果モデルではそのような誤差が独立でない場合にも適切に対応できる。 よって,今後は積極的に混合効果モデルを使用するのがよい。 R では nlme パッケージの lme 関数など(たくさんある)で線型混合モデルによる分析を行うことができる。

一般化線型混合モデル

一般化線型混合効果モデル(generalized linear mixed effects model; 一般化線型混合モデル generalized linear mixed model)は, 上記の,一般化線型モデルと線型混合効果モデルの合体である。つまりは,この両者の特徴をそれぞれ受け継いだ,ハイブリッドなスーパーマンだ。 例えば,ランダム効果を組み込みながら,分布を二項分布としてモデリングできるのだ。夢のよう。

ところが,まさに夢のような話で,世の中,都合の良いことばかりというわけには行かない。 このようにどんどん拡張していく上で壁として立ちふさがるのは,推定が難しくなることである。

ここまで説明してきた現代的な統計モデルの実践適用では,基本的にパラメタを最尤法(あるいは制限付き最尤法,残差最尤法)で推定する。 しかし,モデルが複雑になると,尤度関数の最大化問題は解析的に解けなくなる。 よって,この問題を数値的に解く様々な工夫がなされることになる。 最尤推定値を求める工夫には,今のところ,これ1つで万能という決定的なものはない。 すると,ユーザには,それらの数値計算手法の中から1つ選択するというこれまた平易でない課題が降りかかるのだ。

最近では,最尤推定ではなくベイズ推定も現場で用いられるようになってきた。 経験ベイズ (EB) とマルコフ連鎖モンテカルロ (MCMC) のおかげ。 現場は,背後の哲学なんかは無視して,統計学者から使えると聞けば何でもござれな感じだ。 もしこのサークルに参加されるなら,あなたを悩ませる選択肢をますます増やすことになる。

実のところ,ここまでくるとケースバイケースだから ―私の経験不足もあるが― ,一般的なことが言いづらい。 生じた問題に応じて個別に解決策を探っていくしかなかろう。 何年かすればもっとスマートなヒューリスティックを提示できるかも。

構造方程式モデリング

構造方程式モデリング (structural equation modeling, SEM; 共分散構造分析 covariance structure analysis) とは, 観測変数間の分散共分散(あるいはもっと一般にモーメント)に対して, 特に観測変数だけでなく潜在変数 (latent variables) を含んだ背後の構造を想定し, そこから予測される共分散と実際に収集されたデータの共分散が近づくようモデルのパラメタを推定する方法である。

簡単に言えば,SEM は潜在変数と顕在変数(manifest variables =観測された変数)に係るパラメタを含んだ連立方程式を解く,いわゆる同時方程式モデルである。 同時方程式モデルは経済学などで長年扱われてきたが,潜在変数の psychometric な測定を含んだ点に現代的な SEM の特徴がある。 確認してもらえばわかるように,上で説明している各種の統計モデルではすべて,1つの方程式でモデルを表現している。 SEM はそれを複数の方程式に拡張したものと見ることができる。 そのように連立方程式なモデルは,数式では複雑に見えてしまうかも知れないが, 実践ではそれをパス図 (path diagram) として表現することによりユーザインターフェイスを格段に易しくし,利用者の増大につながっている。

ふつう SEM のための統計ソフトは潜在変数を含まないモデルも扱えるようになっているので, 分散分析モデルもその中で表現でき,パラメタの推定や検定をすることができる。試しに,二元配置分散分析を R の sem でやってみる:

path.diagram(result1)
y <- c(rnorm(16, 1), rnorm(16, 0), rnorm(16, 1), rnorm(16, 2))
x1 <- c(rep(0, 32), rep(1, 32))
x2 <- rep(c(rep(0, 16), rep(1, 16)), 2)

mycov1 <- cov(cbind(y, x1, x2, "x1:x2"=x1*x2))

model1 <- matrix(c(
 '   x1  -> y    ', 'b1', NA,
 '   x2  -> y    ', 'b2', NA,
 'x1:x2  -> y    ', 'b3', NA,
 '    y <-> y    ', 'e1', NA,
 '   x1 <-> x1   ', 'v1', NA,
 '   x2 <-> x2   ', 'v2', NA,
 'x1:x2 <-> x1:x2', 'v3', NA,
 '   x1 <-> x2   ', 'c1', NA,
 '   x1 <-> x1:x2', 'c2', NA,
 '   x2 <-> x1:x2', 'c3', NA
), ncol=3, byrow=TRUE)

library(sem)
result1 <- sem(model1, mycov1, 64)
summary(result1)

R の sem ではこれは次のようにも書ける:

model2 <- matrix(c(
 '   x1  -> y', 'b1', NA,
 '   x2  -> y', 'b2', NA,
 'x1:x2  -> y', 'b3', NA,
 '    y <-> y', 'e1', NA
), ncol=3, byrow=TRUE)

result2 <- sem(model2, mycov1, 64, fixed.x=c('x1','x2','x1:x2'))
summary(result2)

この結果は次の線型モデルでの結果と同じになるはず:

summary(lm(y ~ x1 * x2))

ここまでならわざわざ SEM を使う利点は無いが,SEM で表現できることが理解できれば,その応用範囲が非常に広いことがわかるだろう。

例えば,測定の信頼性を考慮して,潜在変数を被説明変数とした分散分析ができる。

path.diagram(result3)
f1 <- c(rnorm(16, 1), rnorm(16, 2), rnorm(16, 2), rnorm(16, 3))
y1 <- f1 * 0.5 + rnorm(64) # 4回の測定 y1, y2, y3, y4
y2 <- f1 * 0.5 + rnorm(64) # のきなみ信頼性が低い
y3 <- f1 * 0.6 + rnorm(64) #
y4 <- f1 * 0.6 + rnorm(64) #
x1 <- c(rep(0, 32), rep(1, 32))
x2 <- rep(c(rep(0, 16), rep(1, 16)), 2)

mycov2 <- cov(cbind(y1, y2, y3, y4, x1, x2, "x1:x2"=x1*x2))

model3 <- matrix(c(
 '   f1  -> y1', 'a1', NA,
 '   f1  -> y2', 'a2', NA,
 '   f1  -> y3', 'a3', NA,
 '   f1  -> y4',   NA,  1,
 '   x1  -> f1', 'b1', NA,
 '   x2  -> f1', 'b2', NA,
 'x1:x2  -> f1', 'b3', NA,
 '   y1 <-> y1', 'e1', NA,
 '   y2 <-> y2', 'e2', NA,
 '   y3 <-> y3', 'e3', NA,
 '   y4 <-> y4', 'e4', NA,
 '   f1 <-> f1', 'd1', NA
), ncol=3, byrow=TRUE)

result3 <- sem(model3, mycov2, 64, fixed.x=c('x1','x2','x1:x2'))
summary(result3)

MANOVA の結果と比べてみる:

summary(manova(cbind(y1, y2, y3, y4) ~ x1 * x2))

他にも,説明変数を潜在変数にしたり (潜在クラス分析 latent class analysis),多層モデルを作ったり (階層線型モデル hierarchical linear models),成長曲線との交互作用を調べたり (潜在曲線モデル latent curve models),線型の範囲でもやれることは大量にある。

説明変数や被説明変数に潜在変数が使えることは,しばしば心理学,生物学,教育学などでの理論をより適切に反映することができ,望ましい。 例えば,誤差のない測定など無いのだから,指標から測定誤差の影響を除くことはほぼ普遍的に有意義である。ぜひやるべし。 古典的な分散分析,線型モデルは,測定誤差もその他の誤差も全部いっしょくたにして誤差項に押し込んでいるのだ。

一方で,SEM の結果は,ユーザが非常にフレキシブルに指定できるモデルに大きく依存しているので,モデルが適切かどうかが決定的に重要である。 SEM によって得られる恩恵は言わばこの「仮定の強さ」のおかげなので,その仮定が間違っていたら元も子もない。 モデルが複雑になればなるほど,その分野での理論のしっかりした検討を欠かさないようにしよう。

閉検定手順

業界では,仮説族あたりの過誤の確率 ( familywise error rate; FWER ) が重視されている。 要するに,仮説1つ1つではなく,仮説の集まり1つに対して,第一種過誤の確率を一定水準以下に抑えよ,という指令だ。 何を family だと設定するかの問題は横に置いておいて, とりあえず汎用的な対策として閉検定手順 (closed testing procedure) と呼ばれる手続きが提案され (Marcus, Peritz, & Gabriel, 1976) よく使われていて,分散分析とも併せて使えるので,ここに説明しておく。

  1. 関心のある帰無仮説(が表す母数部分空間)の集まり(仮説族) F を定める。
  2. F から部分空間の積 nn について閉じた仮説族 bar F を作る。

    すなわち,AA H_i in F (H_i in bar F) かつ AA H_i, H_j in bar F (H_i nn H_j in bar F) となる最小の bar F をとる。

  3. ある帰無仮説 H_i in F について, H_j sube H_i であるすべての H_j in bar F において水準 alpha で有意な結果を得たとき, H_i を棄却する。

これを守ることで,FWER を alpha 以下に抑えられる。 注目すべき特徴は,各 H_j in bar F の有意性を調べる際に有意水準は alpha でよい,ということ。 それによって,単純な Bonferroni の不等式にもとづく制御(有意水準を検定の個数 k で割る方法)を行うよりは検出力を高く保つことができる。 また,もう一つ注目すべき特徴として,各 H_j in bar F のテストは第一種過誤の確率 alpha を守るならどんな検定法でもよい。 もちろん分散分析を使うこともできる。この点も Dunn-Bonferroni 法と共通である。 したがって一般に,Dunn-Bonferroni 法は出番がない。

FWER 制御の理屈

閉検定手順によって FWER ≦ alpha となる理屈は次の通りである:

まず,n 個の事象 E_1, E_2, ..., E_n について一般に

P( E_1 ^^ E_2 ^^ ... ^^ E_n ) le P( E_k ) qquad "for all" quad 1 le k le n

が成り立つ [point(1)] 。すなわち,複数の事象が「かつ」で結ばれた場合,その生起確率はそれらの内の1つの事象の生起確率を超えることはない。 本来は事象の積には ^^ ではなく nn を使うべきだが, 上記の母数空間の nn と区別しやすくするために,ここでは ^^ を使うことにする。

閉検定手順では,ある仮説 H_i in F を(閉検定手順の意味で)棄却するには, それが部分空間として含む H_j in bar F がすべて(通常の意味で)棄却されなければならない。 この H_j in bar F の集まりを S(H_i) と書くことにする。 つまり,S(H_i) = { H_j in bar F quad | quad H_j sube H_i } である。 仮説の集まり A に含まれるすべての仮説が(通常の意味で)棄却されるという事象を C(A) ,仮説 H_i が(閉検定手順の意味で)棄却されるという事象を R(H_i) と書くことにすると,仮説 H_i が(閉検定手順の意味で)棄却される確率は

P( R(H_i) ) = P( C( S(H_i) ) )

である [point(2)] 。また,S(H_i) に含まれる仮説を G_1, G_2, ..., G_p とすると,

P( C( S(H_i) ) ) = P( C({G_1}) ^^ C({G_2}) ^^ ... ^^ C({G_p}) )

である [point(3)] 。

いま,同時に検討したい関心ある帰無仮説を H_1, H_2, ..., H_n とする。 もしこれらが全部正しくて独立に(通常の意味で)水準 alpha で検定されるなら, "FWER" = 1 - (1 - alpha)^n と批判されるところだ。 bar Fnn について閉じているので, H_1, ..., H_n のすべての共通部分をとった仮説 nnn H_ibar F の要素である。 そしてこの nnn H_iH_1, ..., H_n のどの仮説にも部分として含まれるから, H_1, ..., H_n のどの仮説を(閉検定手順の意味で)棄却するにもその前に nnn H_i が(通常の意味で)棄却されなければならない。 閉検定手順では nnn H_i の(通常の意味での)検定は有意水準 alpha で行われる。 ここでもし, nnn H_i が正しくない帰無仮説であったなら,おそらく alpha よりもだいぶ大きな確率で(通常の意味で)棄却されるであろう。 一方,nnn H_i が正しい帰無仮説であったなら,(通常の意味で)棄却される確率は alpha 以下になる [point(4)] 。

さて,ここまでの話では H_1, ..., H_n について特に断りはなく,正しい仮説でも正しくない仮説でもよかった。 しかし,FWER は,帰無仮説が正しい場合にしか関係がない(そういう概念だから)。 当然ながら,関心のある n 個の帰無仮説のうち真に正しいのはどれとどれかなんて研究者には不明である。 とりあえず,n 個の仮説のうち m 個が正しい,と仮定しよう ( 1 le m le n ) 。 この場合でも,正しい仮説 H_1, ..., H_m について共通部分 nnn^m H_i は必ず S(H_1), ..., S(H_m) のすべてに要素として含まれる [point(5)] 。 すると,H_1, ..., H_m の少なくとも1つ以上が(閉検定手順の意味で)棄却される確率,すなわち FWER について,次の不等式が成り立つ:

{:
( "FWER" ,=, P( R(H_1) vv ... vv R(H_m) ), ),
(  ,=, P( C(S(H_1)) vv ... vv C(S(H_m)) ), cdots "point" (2) ),
(  ,=, P( (C(S(H_1) setminus { nnn^m H_i } ) vv ... vv C(S(H_m) setminus { nnn^m H_i } )) ^^ C( { nnn^m H_i } ) ), cdots "point" (3) "and" (5) ),
(  ,le, P( C({nnn^m H_i}) ), cdots "point" (1) )
:}

このとき nnn^m H_i は正しい仮説なので,nnn^m H_i が(通常の意味で)棄却される確率 P( C({nnn^m H_i}) ) は point(4) のとおり alpha 以下である。 従って FWER は alpha 以下となる。

このことは,m 個の正しい仮説をどのように選んでも成り立つ。つまり,関心のある帰無仮説の集まりのうち本当に正しいものがいくつあろうとどれであろうと, すべての場合においてそれら正しい帰無仮説が1つ以上棄却される確率は alpha 以下となる。

Holm–Bonferroni 法

閉検定手順の応用の1つとして,Holm-Bonferroni 法と呼ばれる手法がしばしば利用されている。 これは,各仮説 H_j in bar F の(通常の意味での)検定にて,Bonferroni 修正を用いるという手法である。

例として,4群のうち1群を対照群として,これと残りの3群との平均値を比較したい場合を考える。関心のある仮説族 F は,

{:
( H_1:, mu_1 = mu_4 ),
( H_2:, mu_2 = mu_4 ),
( H_3:, mu_3 = mu_4 )
:}

である。ここから nn に関する閉包 (closure) bar F

bar F = { H_1, H_2, H_3, H_1 nn H_2, H_1 nn H_3, H_2 nn H_3, H_1 nn H_2 nn H_3}

となる。例えば,H_1 を(閉検定手順の意味で)棄却するには,H_1 が部分として含む,

H_1, H_1 nn H_2, H_1 nn H_3, H_1 nn H_2 nn H_3

をすべて(通常の意味で)棄却しなければならない。H_2 については

H_2, H_1 nn H_2, H_2 nn H_3, H_1 nn H_2 nn H_3

の(通常の意味での)棄却が必要である。ここで,複数の極小仮説 (minimal hypotheses) が nn でつながっているもの,例えば H_1 nn H_2H_1 nn H_2 nn H_3 などを(通常の意味で)第一種過誤の確率 alpha を守って検定するために,Bonferroni 不等式による調整を伴って対比較(例えば t 検定)しようと考える。 よって,H_1 nn H_2 nn H_3 を(通常の意味で)棄却するには,1回の t 検定の有意水準は alpha // 3 を用いることになる。

H_1 nn H_2 nn H_3 を(通常の意味で)棄却するためには,H_1, H_2, H_3 のどれかで alpha // 3 を下回る p 値を得ればよい。 よって,実践的には,H_1, H_2, H_3 のすべてで t 検定の p 値を算出して, それらを比べて一番小さな p 値を出す仮説を特定し(仮に H_1 としよう), その p 値が alpha // 3 を下回っているか確認する。 もし下回っていれなければ,H_1 nn H_2 nn H_3 は棄却できない。 この例の場合,H_1 nn H_2 nn H_3 の(通常の意味で)棄却は H_1, H_2, H_3 のすべてにおいて(閉検定手順の意味での)棄却に必要だから,H_1 nn H_2 nn H_3 が棄却できないなら H_1, H_2, H_3 はすべて(閉検定手順の意味で)棄却できないことになり,閉検定手順は終了する。

もし H_1p 値が alpha // 3 を下回っていれば,H_1 の t 検定の結果はその他に H_1 の(閉検定手順の意味での)棄却に必要なすべての仮説(H_1, H_1 nn H_2, H_1 nn H_3 )を棄却するだろう。なぜなら,Bonferroni 法に従えば,他の仮説の検定では alpha // 3 よりも大きな有意水準( alpha // 2 または alpha )を用いるからである。 このことから,H_1 の t 検定によって H_1 nn H_2 nn H_3 が(通常の意味で)棄却されたなら,即刻,H_1 は(閉検定手順の意味で)棄却されることがわかる。

この時点では H_2H_3 についてはまだ結論が出ない。 次に H_2H_3 の t 検定の p 値を比べて,小さい方を特定して(仮に H_2 としよう),残りの仮説の棄却に挑むことになる。 残りは H_2, H_3, H_2 nn H_3 だが,H_2 に関係し一番狭いのは H_2 nn H_3 なので,これに取り組む。H_2p 値が有意水準 alpha // 2 を下回っていなければ,そこで終了。 下回っていれば同様に,H_2 は自動的に(通常の意味で)棄却されるので,(閉検定手順の意味で)H_2 が棄却される。 最後に H_3 について有意水準 alpha で検討すれば終わりである。

このようにして,狭い( nn でたくさん連結された)仮説を検定する方法として Bonferroni 法(他にも方法はあるのに)を選択したとき, 上記の閉検定手順の一般論よりももう少し単純な規則でこの手順を表現することができる:

  1. 原子的な( nn の無い)帰無仮説の数が m のとき,その中で一番小さな p 値を得た仮説 H_1 について alpha // m を下回っているかどうかを調べる。下回っていなければ終了。下回っていれば H_1 を棄却して次へ。
  2. 二番目に小さな p 値を得た仮説 H_2 について alpha // (m - 1) を下回っているかどうかを調べる。下回っていなければ終了。下回っていれば H_2 を棄却して次へ。
  3. 同様に最後の1つの仮説まで続ける。

これが広く知られている Holm-Bonferroni 法の手順である。R の p.adjust.methods のデフォルト。

Shaffer の方法

別の例として,4 群の総当たり対比較をしたい場合を考える。関心のある仮説族 F は,

{:
( H_12:, mu_1 = mu_2 ),
( H_13:, mu_1 = mu_3 ),
( H_14:, mu_1 = mu_4 ),
( H_23:, mu_2 = mu_3 ),
( H_24:, mu_2 = mu_4 ),
( H_34:, mu_3 = mu_4 )
:}

よって,nn に関する閉包 bar F

bar F = { H_12 , H_13 , H_14 , H_23 , H_24 , H_34 ,
 H_12 nn H_23 , H_12 nn H_24 , H_12 nn H_34 , H_13 nn H_24 , H_13 nn H_34 ,
 H_14 nn H_23 , H_23 nn H_34 , H_12 nn H_23 nn H_34 }

例えば H_12 を(閉検定手順の意味で)棄却するには,

H_12, H_12 nn H_23, H_12 nn H_24, H_12 nn H_34, H_12 nn H_23 nn H_34

をすべて(通常の意味で)棄却する。

上節のように3つ以上の群が関わる仮説では Bonferroni 法を使うとしよう。H_12, ..., H_34 全部の p 値を計算する。仮に,H_12 が一番 p 値が小さく,alpha // 6 を下回るなら,H_12 は(閉検定手順の意味で)棄却される。

次に p 値が小さいのは H_13 だとしよう。残りの仮説のうち,H_13 の(閉検定手順の意味での)棄却に必要なのは

H_13, H_13 nn H_24, H_13 nn H_34

である。すると,有意水準 alpha // 3H_13 nn H_34 ( なぜならこれは H_13 nn H_34 nn H_14 と同じである)を検定するのが次のステップである。

上記の Holm–Bonferroni 法(最後に書き直した手順)に従えば,問答無用で alpha // m の分母を 1 ずつ減らしていくので,二番目に p 値が小さい仮説を検討する際は,有意水準は alpha // 5 になる。 これが Shaffer の方法と Holm–Bonferroni 法との違いである。Shaffer の方法もたんに閉検定手順の一般的な理屈に従っているのであるが,Holm–Bonferroni 法はそれを上記のような憶えやすい手順に直したため,関心のある帰無仮説のセットによっては閉検定手順から(保守的過ぎる方向へ)ズレてしまうのだ。

よって Shaffer の方法のほうが Holm–Bonferroni 法よりも検出力が高くなる。