分散分析の視覚的理解ー多重比較禁止令

分散分析の視覚的理解ー多重比較禁止令
Page content

他のメンバーみたいに真面目っぽい記事を書いてみた Jan.23.2021


警告:これは由緒正しき正規分布帝国の布教活動です。ノンパラメトリック自治領の方は(以下略)

※わかりやすさ重視なので厳密性はゆるゆるです。素人がふわっと書いてるだけなので絶対に間違いが入っています。決して教科書ではないのでご了承ください。

※一応t検定の基礎あたりまでを勉強している前提で書いてますが、筆者も数学は苦手なので、えっと、、、これ以降の内容が電波系の戯言だと感じなかったら大丈夫です。


世界の始まり、由緒正しき正規分布皇帝は”平均あれ”と言われた。

すると人々の目は開かれ、あらゆる事物の平均が目に映った。人々は物事の本質が見えたと喜んだ。

しかし、なぜ平均値がモノゴトの本質を表すのか?それはここが中心極限定理に守られた由緒正しき正規分布帝国だからである。

時は流れ、我々は今でもこの言葉に憑りつかれている忠実である。 これを我々はこう表しました。

Y = μ + ε

(日本語訳:)こいつは平均(μ)くらいなはずですが、平均からズレてる分については誤差(ε)としましょう。

まず、世界全体の平均があり、世の全てのモノゴトは平均からの距離、すなわち誤差の大きさで説明される。

智に働けば角が立つ。情に棹させば流される。

そりゃあうまくいかねえわけです。

そこで宗教革命!

効果なるものを唱えたヤツがおるわけです。

Y = μ + ε

Y = μ + A + ε

(日本語訳:)こいつは平均(μ)からAの影響分ズレています。それで説明できない部分は誤差(ε)としましょう。

人々はAを基準にモノゴトをグループに分け始めました。そしてグループごとの平均でグループどうしを比べ始めました。こうやって世界は分断されてしまいました。

でもグループの平均の違いをどう考えましょうか?

どれくらい違いがあれば私はあなたと違うんですか?絶対に違うって?あなたはいっつもそうよね…

平均ってなんでしょうか?

例えば任意の2つの実数を足し合わせて2で割ると、2つの値のバラツキは相殺されますよね。もし元の値がわからなくなったら平均値だけ見ても元の2つの値はわかりません。当たり前の話ですね。

平均=4 → (2, 6), (3, 5), (1, 9)・・・無限に存在

例えば平均=4を考えたとき、もし元の値が(3, 5)ならだいたい4ぐらいだって言っても良さそうです。でも、もし元の値が(-2006, 2014)とかだと、だいたい4ぐらいって言うのは間違ってはいないものの、あんまり実用的じゃないですよね。

値のバラツキが大きいほど、総平均だけで語る現実性はなくなっていきます。そこで世界をグループに分けてみたのです。

1つのグループの平均を取ると、グループとは無関係なバラツキは全部相殺されます。平均で見ると、”グループ”が”全体”からどれくらいズレているかがわかります。当たり前の話ですね。

各グループの平均をあつめてくると、グループとは無関係の影響を全部相殺して、純粋に”グループだけ”のバラツキがわかります。

この”グループのバラツキ”が大きければグループによって違いがありすぎる、すなわち世界を総平均だけでは語れないと言うことができます。

逆にグループ分けなんて無意味なら”グループのバラツキ”は大きくないはずですよね?そんなの世界に不要な分断と対立を生むだけです。

問:グループ分けには本当に意味があるのか?

そこでこの”バラツキ”をどう比べようか、って話なんですが、、、バラツキを指標化した値は皆さんご存じ分散ですね。分散でモノゴトを評価する分析法なのでANOVA(分散分析、analysis of variance)と言います。


one-way between ANOVA

リアルな数字を登場させてみましょう。

例えば、Itamae・ぽ・Johsonの3人にそれぞれチュパカブラの卵を3つずつ渡して育ててもらって、1年後に体長を測って成長を見たとしましょう。

誰が育てたかによってチュパカブラの成長には差があるのでしょうか? one-way between ANOVA(対応の無い一元配置分散分析)で分析してみましょう。

IV(independent variable, 独立変数/説明変数):飼い主(水準はItamae, ぽ, Johsonの3条件)

DV(dependent variable, 従属変数/目的変数):チュパカブラの体長(長さ㎝、比尺度)

実際のANOVAに入る前にY=μ+A+εのモデルの意味を図示して確認してみましょう。Y=μ+A+εは行列式なので9匹のチュパカブラについて9本の関係式が成り立ちます。

Y = μ + A + ε

113 = 120 - 2 - 5

125 = 120 - 2 + 7

116 = 120 - 2 - 2

106 = 120 - 6 - 8

120 = 120 - 6 + 6

116 = 120 - 6 + 2

133 = 120 + 8 + 5

124 = 120 + 8 - 4

127 = 120 + 8 - 1

下の図では例として、Johsonの124㎝のチュパカブラに注目しています。でも各点について一つずつ上の式をあてはめてみるとY=μ+A+εの全体像が見えてきます。

※図の中の正規分布は模式的に置いたもので、形や大きさは正確ではない

これがGLM(general linear model, 一般線形モデル)としてのデータの捉え方です。GLMは正規分布帝国の最強兵器ですから、前提条件としてバラツキが全て正規分布に従います(※)。

Aは3つの条件だ!って思いこみがちですが、見方を考えれば正規分布から抽出したN=3のサンプルだと見ることもできるわけです。後で出てきますがこれ超重要ポイントです。

※同じGLMでもgeneralized linear model(一般化線型モデル)は違います。こいつは正規分布帝国すらを支配下に置いてしまった終末兵器です。実はさらにそれを超える兵器が…

以降、上の図のイメージを念頭に置きながら読んでいくと少しは理解しやすくなると思います。

まずは分散で評価するということなので、まずはデータの各値を下のように変形します。

Y = μ + A + ε

Y - μ = A + ε

普通の分散の定義はY-μ(偏差)の二乗の合計をN数で割ることです。

しかし今回はt検定で使った不偏分散と同じで、データ一つあたりのバラツキは単純にNで割るだけでは出せません。

なのでデータ一つあたりで考えるのは後にして、まずはバラツキの合計量で見ていきしょう。

モデルの各項のバラツキは、各値を二乗して合計した二乗和で考えます。この二乗和のことを別名SS(偏差平方和、sum of squared deviations/squared sum)とも言います。

Σ{(Y-μ)²} = Σ(A²) + Σ(ε²)

同じ式を書き換えると

SS(total) = SS(A) + SS(ε)

まず主効果Aを条件ごとにまとめて(これをプールすると言います)、それぞれで平均を出します(水準平均)。

総平均μ = 120

水準平均(Itamae) = 118

水準平均(ぽ) = 114

水準平均(Johson) = 128

次に総平均μと水準平均の差分をとってAを出します(偏差)。これは言い換えればチュパカブラの個体差を無視した時、飼い主によってμからどれくらい変動するか、という値です。

A₁ = 118 - 120 = -2

A₂ = 114 - 120 = -6

A₃ = 128 -120 = 8

残差εはε=Y-μ-Aで定義されます。モデルのε以外の項が全部埋まったあと、最後に丁度良く等式が成立するように都合が良い数が入ります。だから”残された差”で残差と呼ばれていて一匹ずつ違う値になります。

ここで、モデルと対比させながら考えます。

Y = μ + A + ε →→ SS(total) = SS(A) + SS(ε)

①もし条件による差が全く影響していなければ(=誰がチュパカブラを育てても関係なく同じくらいの大きさになるなら)、条件間の差が出ずに条件平均は総平均と全く同じになるはずです。残る要因はチュパカブラの個体差だけです。 Y = μ+ε → SS(A)=0, SS(total)=SS(ε)

②もしデータ内のバラツキが全て条件による影響Aだけで説明される(=チュパカブラがどれくらい大きくなるかは誰が育てたかだけで決まる)としたら、各条件内では全く同じ値を取っているはずです。チュパカブラを育てた人が同じなら必ず同じ大きさに成長します。 Y = μ+A → SS(total)=SS(A), SS(ε)=0

よって、SS(total)を構成するSS(A)とSS(ε)を比べて主効果Aの強さを評価しようという考え方をします。(2つの比→効果量:後述)

でも、SS(A)とSS(ε)を直接比較するはできません。上の図にもあるように主効果の強さによってSS(total)の大きさが変化しますよね?

普通の分散ならSS(=全体のバラツキ量)をNで割れば一匹あたりのバラツキがわかって他のサンプルの分散と比較できます。これが標準化です。

でもSS(A)を標準化するには何で割ればいいか、、、と考えていると普通ならAの水準の数で割れば良さそうです。

ここで最初の方で触れた不偏分散を思い出してください。t検定で出てきた推定値の分散はNではなく(N-1)で割ってましたよね?

そうです、Nのかわりにdf(自由度)で割るのです。dfはそれぞれの項におけるN数的な存在なのです(これまた誤解を生みそうな表現ですが…)。

まとめ:SSをdfで割ると不偏分散が出ます。これをMS(mean of squares, 平均平方)と呼んでます。

SS/df = MS

MS(A)=SS(A)/df(A)=312/2=156

MS(ε)=SS(ε)/df(ε)=224/6=37.33

SS(A)とSS(ε)を標準化してMS(A)とMS(ε)とすることでこの2つを比べることができるのです。さて具体的にどう比べましょうか。

ここで思い出してください。ここは由緒正しき正規分布帝国です。全ての項について正規分布だと言えるのです(拍手)。

正規分布するサンプルの分散(自由度df₁)と、同じく正規分布のサンプルの分散(自由度df₂)の比をとったらなんか決まった確率分布になりそうじゃないですか?それがF分布です。

source:wikipedia

はい、長くなるのでかなり説明を端折ってます(※)。端折ってしまったのでt検定は独立変数が2水準しかない特殊形のANOVAなんですよ~の話もできません泣。興味がある方は今までの説明を二人の場合で考えてみてください。最後にt(df)=√F(1,df)でつながるはずです。

※統計を学ぶ時、教科書上ではだいたい正規分布→χ²分布→F分布→T分布の順番で載っていますが、分析方法を学ぶ順番はだいたいT検定→ANOVA、と逆なので皆チンプンカンプンになります。正規分布帝国を歩む最初の感動はやはり、ANOVAを通じてこれらの分布を教科書通りに一本の繋がった流れで理解することです。興味ある方は是非勉強してみてください。

F分布はF値(=分散の比)のpdf(確率密度関数)です。だからMS(A)/MS(ε)=Fを出してF分布に照らし合わせれば分布確率がわかります。t検定の時と全く同じ考え方ですが、ANOVAの場合の対立仮設はMS(A)»MS(ε)なので上側だけ見る片側検定です。

source:BellCurve ※今回の例は上図と自由度が違う

F = MS(A)/MS(ε) = 156/37.33 = 4.18 ~ F(2, 6)

有意水準をα=5%とするとF分布表より

Fα(2, 6)=4.76 > 4.18

飼い主の違いによる体長には有意な差は認められませんでした(※)。

以上、one-way between ANOVA おしまい!

、、、と安心する前に

統計的仮説検定は考え方からして根本的に不完全な分析法です。特にN数に対して敏感です。t検定でCohen’sdを見たようにANOVAでも効果量を必ず確認しましょう。

ANOVAでまず出てくる効果量は分散説明率、η²です。

η² = SS(A) / SS(total)

式の形を見れば意味は一目瞭然です。ちなみにチュパカブラの例で計算してみるとη²=.58、かなり強いです。

しかし、η²は完璧ではありません。

もし仮にIVが増えてしまったらどうなるでしょうか?モデルによってη²の解釈が相対的に変わってしまうのです。

そこで偏η²参上。

偏η² = SS(A) / {SS(A)+SS(ε)}

これでモデルに依存しない効果量になりました。ちなみにone-way between ANOVAのモデルにはそもそも主効果が一つしかないので、η²=偏η²になります。

結論:効果量も大きいので主効果は無視できない。

晴れて分析おしまい!

※当たり前ですが今回の分析は非現実的です。普通に考えてチュパカブラ3匹ずつじゃあんまり意味無いです。現実的には3人にはもっとたくさんのチュパカブラを育ててもらわないと統計的に意味のある結果は出ないでしょう。(具体的に何匹が妥当?という話は長くなるので省略します)


さてここで、もし仮に、チュパカブラの例でANOVA有意となった場合を考えてみましょう。この時、飼い主によって違いが出るとまでは言えますが、飼い主の中で具体的に誰のチュパカブラが他の誰かより大きく成長するのかはわかりません。

そこで多くの場合、誰と誰の間に差があるのかを突き止めるためにpost-hoc(事後検定, 多重比較)というものを行います。でも、ただt検定で比較を繰り返していくと、比較を重ねるごとにαエラーがかけられてしまいます。

例えば、チュパカブラの例で飼い主の総当たり比較を行うために₃C₂=3回繰り返しt-testを行うとします。

比較①:Itamae - ぽ

比較②:ぽ - Johson

比較③:Johson - Itamae

t検定の有意水準をα=5%とすると、3回繰り返した後のαエラーの確率は

1 - α = (1 - 0.05)³ = 0.86

∴ α = 0.14

と10%以上になってしまいます。

なのでこれを補正するために様々な手法が考案されてきました。多重比較についてはまた後で言及します。


two-way/within ANOVA

one-way between ANOVA がわかればtwo-wayだろうがwithinだろうが基本はどれも同じです。なので細かいことは教科書に投げて今回は大まかなことだけをかる~く触れます。

どんなANOVAも本質的には全く同じ操作をします。

①モデルを考える

②各項のSS→MSを求める

③効果のMSをMS(ε)で割って効果ごとにF値を求める

④求めたF値でF検定をする。(+偏η²も見る)

two-way between ANOVA(対応の無い二元配置分散分析)は2つのIV、AとBそれぞれの効果、それに加えてAとBの交互作用というものを考える必要があります。

Y - μ = A + B + AB + ε

交互作用ABというのはAとBの相乗効果のことです。プロフェッサーの説明↓

風邪をひいたとき、風邪薬を飲みますね。市販の風邪薬にはエフェドリンとリン酸ジヒドロコデインって麻薬成分が含まれてるんだけど、それぞれ単独で飲んでも大したことはないのが、あわせて呑むとハイになれるんですね~。

OTC薬の量で実際に効果出るのかは筆者も気になるところではあります(笑)。

ちなみに交互作用は強く出るだけ以外にも拮抗する形のこともあります。例えば、テトロドトキシンとアコニチンを同時に摂取すると効果が拮抗してすぐには死なないよ、みたいな話です。

SSを求めるのも効果の数が増えただけでやることはone-wayと同じです。ひとつひとつ水準ごとにプールして平均を出す、ってのを効果ごとに行います。

A → A₁ A₂ A₃ ・・

B → B₁ B₂ B₃ ・・

A×B → A₁⋂B₁ A₁⋂B₂ ・・A₂⋂B₁ A₂⋂B₂・・

先ほどの例をちょっと手を加えて、3人にチュパカブラを3匹ずつ育ててもらうかわりに1人につき、メキシコ産のチュパカブラ2匹とチリ産のチュパカブラ2匹の計4匹を育ててもらうことにしました。IVは飼い主(効果A)と産地(効果B)です。

モデルのうち、A + B + AB の部分だけを引っ張り出すと下図のようになります。ABのdfはAとBのdfをかけると出てきます。

Y - μ = A + B + AB + ε

SS(total) = SS(A) + SS(B) + SS(AB) + SS(ε)

μ=120

Aの水準平均:(itamae, ぽ, Johson)=(118, 114, 128)

Bの水準平均:(メキシコ, チリ)=(119, 121)

ABの水準平均:

itamae(メキシコ, チリ)=(115.5, 120.5)

ぽ(メキシコ, チリ)=(113, 115)

Johson(メキシコ, チリ)=(128.5, 112.5)

これであとの操作は全く同じです。

②A, B, AB, εのSS→MSを求める

③A, B, ABのMSをMS(ε)で割って3つのF値を求める

④求めたF値で1つずつF検定をする。(+偏η²も見る)

以上、two-wayおしまい!


within ANOVA(対応のある分散分析)については、実は普通のbetween ANOVAのモデルに効果をもうひとつ追加しただけです。被験者をIVとして追加して、一人ひとりの被験者を条件と捉えているだけなのです。

Y = μ + A + S + ε

one-way within ANOVAでは上のように被験者効果Sというものが加わります。被験者効果と言っていますが、この先の取り扱いはtwo-way between ANOVAと全く同じです。

(期待してる突っ込み):主効果Aと被験者効果Sの交互作用ASは考えなくても良いの?

交互作用ASは存在します。実はone-way within ANOVAの場合は交互作用ASは残差εに一致します。だからY=μ+A+S+ASと書いてε=0とすることもできます。

ちなみにone-way betweenも同様に考えると残差εはチュパカブラの個体差だけが影響しているバラツキなので実は誤差ε=被験者効果Sなのです。なのでtwo-way within ANOVAも残差の意味が違います。

これは各項の自由度をひとつずつ考えてみると納得できます。実は自由度について考えるというのは、”SSを切り分ける”行為についてすご~く本質的な議論に繋がるのですが、その話は教科書に投げます。

以上、within おしまい!



Contrast

さて、ANOVAで有意となった場合、普通はpost-hoc(事後検定、多重比較)をします。UMINのサイトにわかりやすくまとめてあるのでここでは説明しません。

http://plaza.umin.ac.jp/~beehappy/stat/com-ph.html#post-hoc

それにしてもさすが臨床試験やりまくってるだけありますね~。

でも今回は違うアプローチを紹介します。

ここから先の内容は残念ながら日本国内ではあまり浸透していません。

ですが!

この手法は実験統計の世界ではメジャーな存在です。 そして大学の実験計画法の授業では必ず学びます(少なくとも欧米では)。

それがコントラストコーディング。

それでは、まず根本的なところから考えてみましょう。

ANOVAはGLM(一般線型モデル)の特殊形です。post-hoc法が論理的に誤り、というのはANOVAの回帰モデルを全く無視した解析を無理くり追加で組み合わせているからです。

チュパカブラの例を考えると、もしpost-hocで有意水準を厳しめに調整したt検定を3回行うと(Bonferroni法)、3つの比較のSSの合計値は元のモデルのSS(A)の1.5倍に膨れ上がって一致しません。

・このIVは明らかにDVに強い影響がある

・この群とこの群には明らかな差がある

post-hocを使うとこれらは全く違う分析になってしまうということです。

そもそも論、ANOVAというのはSSを我々に都合が良いように切り分けているだけです。ならば、主効果のSSをさらに分解すれば回帰式を崩さずに条件間比較ができるのではないか?という考え方が第一選択であるべきです。

one-way between ANOVAを例に見てみましょう。

今まではこう書いてました。

Y = μ + A + ε

ここで、中身は全く変えないまま、主効果をこう書き直してみます。

Y = μ + βX + ε

なんか、GLMの基本形(Y=βX+ε)っぽくなった!って思った方は鋭い!

そうなんです。μ = β₀X₀と考えるとGLM基本形そのものになって、βX=(バラメータ行列×デザイン行列)なのです。

話を戻してβXを以下のように分解してみましょう。

A = βX = β₁X₁ + β₂X₂ + ・・・ + β(df)X(df)

主効果Aをdfの数の項に分解すれば条件間の比較ができるってことです!(あとで図示します)

→しかし、うまく条件間比較をするためにはデザインマトリックスをうまい具合に設定する必要があります。この行列Xをどう設定するか(コーディング)が非常に重要なポイントなのです。

コーディングには、いくつか種類があります(例えば、名義尺度変換で使うダミーコーディング)。ここで使うのはANOVAのモデルを崩さないコントラストコーディングです。コントラストコーディングには以下の2つの特徴があるます。

①全ての条件が比較が比較に含まれていること

②各比較が独立であること

βXは自由度の数しか項がないのでコントラストを組む際はデータの背景にある論理的文脈を考慮する必要があります。

※逆に総当たりの比較をしなければならない場合はpost-hocで多重比較が用いられます。

説明だけじゃピンとこないと思うので、まずは、2条件の比較(independent t-test)について考えてみましょう。


t検定とは、A1、A2の2条件について以下のようなコントラストでANOVAを行うということです。Xは(1,-1)で足し合わせると0なので独立な比較です。

コントラストはX = (1, -1)となっていますが、逆にしても符号が逆転するだけで解析は本質的に変わりません。

これを以下の式にあてはめます。

Y = μ+βX+ε

ここで、Xは列ベクトルですので実際には条件別にこういうことになります。

A1 : Y = μ +β +ε

A2 : Y = μ -β +ε

これを図で表すとこうです↓。※ただし、A1>A2→β>0、A1<A2→β<0。残差εはA1またはA2から各点へのバラツキ。

各条件は平均からβズレているということです(A1とA2でNが違うとそのままでは成り立たない)。この時の帰無仮説はβ = 0、すなわちA1 = A2 = μです。

でも実際にこれでANOVAをしようと思うと困りますよね。SS(A)はSS(βX)ということになりますが、βがわかりません。

SS(total) = SS(βX) + SS(ε)

※回帰分析として捉えるとSS(total)→SS(Y)、SS(A)→SS(reg.)などと表記を変えることがありますが、式の中身は全く同じです。ここでは混乱を避けるために表記は変えません。

でも心配しなくてもβの導出ができなくてもANOVAでは困りません(βの一般解はYにXの一般化逆行列をかけます)。

導出は省略しますが、二変数の積の合計SP(sum of products)を使ってSS(βX)を出します。

SP(XY) = ΣXY

SS(βX) = (SP(XY))² / SS(X)

これで今までと同じようにdfで割ってMSを出して、そこからF値を出すことができます。


ここまでは数字や文字が違うだけで普段のindependent t-testと全く同じ形です。ここでさらに、3条件(one-way between ANOVA)についてみてみましょう。

以下のようなコントラストを用いたANOVAを考えてみましょう。

X1ではA1とA2&A3を比較して、X2ではA2とA3を比較します。X1もX2も縦に見るとそれぞれ総和が0なので独立な比較です。

ちなみにコントラストの組み方としてはこれ以外にも、、、

①A1⇔A2②A3⇔A1A2、とか、①A1A3⇔A2②A1⇔A3、みたいなのも考えられます。

X1、X2をどう組んでそれぞれ何を比較するかは実験や解析の論理的な文脈で決めます。(※詳しく話すとトレンド分析などの話もありますが、今回は省略)

上のコントラストに戻って、これも先ほど同様にGLMの式に入れます。

Y = μ+β₁X₁+β₂X₂+ε

Xを各条件ごとにバラして並べてみると、、

A1:Y = μ+2β₁+ε

A2:Y = μ – β₁+β₂+ε

A3:Y = μ-β₁-β₂+ε

各条件が平均からβ2つ分ズレているという形で表記されていることがわかります。df=2個に分解できたということです。

これも図に表してみると以下のようになります。※必ずしも図のようにA1>A2>A3(すなわちβ₁,β₂>0)だとは限らない

これはすなわちX₁について見ると、

帰無仮説β₁=0を棄却→A1とA2&A3には差がある

と言うことが言えますし、X₂も同様に、

帰無仮説β₂=0を棄却→A2とA3には差がある

と言える、というロジックです。

あとは先ほどと一緒で、SSをだしてdfで割ってMS出してF値出すの流れです。

SS(total) = SS(β₁X₁) + SS(β₂X₂) + SS(ε)

SS(β₁X₁) = (SP(X₁Y))² / SS(X₁)

SS(β₂X₂) = (SP(X₂Y))² / SS(X₂)

2つのSS(βX)はコントラストに対応しているのでそれぞれF検定を行うと群間比較ができます。また、通常のANOVAはSS(A)についてF検定を行うことと同じです。

SS(A) = SS(β₁X₁) + SS(β₂X₂)

これ、post-hocより遥かに美しいですよね??

結論:統計に強い現場で理由もなく”とりあえず”で多重比較を持ち出すと人権を失います。

以上、おしまい。


まだまだ2-wayだのwithinだのありますが、基本は全部同じ考え方です。そしてその先には混合モデルや階層モデルといった応用、MANOVA、ANCOVA、ANACONDA(ちがう)、、、ANOVA系のバリエーションがあなたを待っています泣。

フワッと理解が目的だったので前提条件(同じN数、等分散性、正規性)が崩れた場合の対処や統計ソフトの使用法やSCFによるSS計算などなど実運用に関することは相当端折りました。なので実際に使用する方は是非教科書を読んでください。

で、素晴らしいことにめっちゃわかりやすい名著がネット上で公開されているんです!良い時代~

https://www.researchgate.net/profile/Barbara_Tabachnick/publication/259465542_Experimental_Designs_Using_ANOVA/links/5e6bb05f92851c6ba70085db/Experimental-Designs-Using-ANOVA.pdf