年末年始にやっていた実験データの分析の中で、分散分析を何回も行ったのですが、Rで分散分析をやるときに基本関数では平方和のタイプを選べないんですよね。
簡単な内容なのですが、平方和のタイプを選びながら分析する方法を、メモしておきます。
あわせて、平方和のタイプってそもそもなんだっけ？ってのも、自分のためのメモを兼ねて書いておきます。詳しくは統計の解説書を当たるべきだと思いますが。私は南風原教科書で復習しました（そういえばこれの続編が最近出版されて、いまパラパラと読んでいます）。
重回帰分析の場合
まず、重回帰分析を行って、回帰係数の有意性を検定する場合がありますよね。ある独立変数の係数がゼロ（つまりあってもなくても従属変数に影響なし）という帰無仮説を検定するわけです。で、t検定で行う場合が多いとは思いますが、F検定でやることもできます。
その場合、検定対象となる独立変数（複数の場合は変数群）を投入する前と後で、分散説明率（R^2）がどんだけ増えたかを考えます。Texを書くのがめんどうなので日本語で説明するとｗ、「（1-変数追加後の分散説明率）÷（サンプルサイズ-追加後の変数の個数-1）*1」を分母とし、「分散説明率の増分÷増えた変数の個数*2」を分子とする検定統計量が、F分布に従うことを利用するわけです。
ここでポイントは、「その変数を新たに投入するとモデルの分散説明率がどんだけ増えるか」を問題にする検定統計量ってのは、その変数を何番目に投入したかで変わるということです。なぜならこの値は、「その変数から他の全ての変数の影響を取り除いた成分と、従属変数との、部分相関」に基づいて算出されてるからですね。
つまり、その変数を投入する前にすでに投入してある変数群との相関の分だけ、分散説明率があらかじめ奪われた状態で、上記の検定統計量の中に登場してくるということです。だからその「すでに投入してある変数群」が何なのかによって、値が変わるわけですね。これは不公平ですｗ
研究目的や理論上の理由から、明らかに「優先的な変数」がある場合はべつに良いとされているのですが*3、とくに理由がないのであれば、この不公平は問題になります。Rで計算するときも、同じ変数をモデルの前の方に入れるとその係数が有意になりやすくなったりして、恣意的であると言えたり言えなかったり。
こういう問題があるので、「その変数を最後に投入した時の分散説明率の増分」をベースに算出するというルールを定めたのが、「タイプIIIの平方和」ってやつですね。（あとで検索とかで引っ掛けることを考えると、「タイプIII」と書くべきなのか「タイプⅢ」と書くべきなのか、はたまた「タイプ3」のほうがいいのか、よくわかりませんが。）
この方法だと、不公平さがなくなるので気楽に使えます。逆に、単純にモデルの先頭から順番に計算していくのは「タイプI（タイプ1）の平方和」となります。
分散分析の場合
２要因以上の分散分析の場合も同じ問題が生じます。分散分析はカテゴリ変数の水準間の平均値差の検定に使われるわけですが、統計モデルとしては、重回帰分析と同じ感じで、線形モデルに「水準数−１個」のダミー変数を投入したものとして表現できます。いわゆる分散分析は、この線形モデルの切片を無視して、回帰係数がゼロでないかどうかを検定してるわけですね。
水準が３つ以上ある場合は、ダミー変数は「ダミー変数群」となり、このダミー変数群を投入することによって分散説明率がどんだけ増えるかを考えることになります。
で、２要因以上の分散分析で、全てのセルのサンプルサイズ（被験者数、ケース数）が同じであれば、じつは要因間の相関がゼロになるので、「その変数をモデルに登場する順番」は考えなくていいことになります。重回帰のところで説明したように、相関があるから「先に投入した変数によって平方和が一部奪われてしまう」のが問題なのであって、相関がなければどういう順番で投入しても同じです。
ところが、往々にしてセル間のサンプルサイズは異なってくる（Unbalanced Desingと言われる）わけで、どうしてもタイプIではつらい。そこで、心理学の研究では、タイプIIIの平方和がよく使われることになっているようです。
Rで平方和のタイプを選択する
私はいつも、Rで分散分析をやる場合は、線形モデルのlm()関数の計算結果をanova()に渡す形で、
> x1 <- rep(c("A", "B", "C"), 100) > x2 <- c(rep("F", 100), rep("M", 200)) > set.seed(10) > y <- rnorm(300) > model1 <- lm(y ~ x1 * x2) > anova(model1) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x1 2 6.094 3.04698 3.3327 0.03705 * x2 1 0.684 0.68381 0.7479 0.38783 x1:x2 2 0.367 0.18369 0.2009 0.81809 Residuals 294 268.791 0.91425 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
みたいにやります。
ところがこのanova()関数（や、aov()関数）だと、平方和のタイプを選べないらしく、強制的にタイプIになっているようです。
anova()やaov()がタイプIを使っているということを確認するには何を見ればいいのか分かってないのですが、適当にググるとそれっぽいことが書いてあるので、たぶんそうなんでしょう。（参考1・参考2）
ためしにx1とx2の順番を入れ替えてみると、
> model2 <- lm(y ~ x2 * x1) > anova(model2) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) x2 1 0.713 0.71337 0.7803 0.37778 x1 2 6.064 3.03220 3.3166 0.03764 * x2:x1 2 0.367 0.18369 0.2009 0.81809 Residuals 294 268.791 0.91425 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 >
すこし値が変わりましたね。これはX2の水準「F」が100個であるのに対して「M」が200個というUnbalanced Design（セルの大きさが不釣り合い）になっていて、相関が発生してしまうからです。
それで、調べたら{car}というパッケージがあって、その中に大文字から始まるAnova()という関数があり、これを使えばタイプII平方和やタイプIII平方和が計算できるみたいです。
とりあえず以下のとおりタイプIII平方和でやってみます。
後日のエントリに書いたように、options(contrasts = c("contr.sum", "contr.sum"))ってのを設定しとかないといけないです。
> library(car) > options(contrasts = c("contr.sum", "contr.sum")) > Anova(model1, type=3) Anova Table (Type III tests) Response: y Sum Sq Df F value Pr(>F) (Intercept) 1.916 1 2.0952 0.14883 x1 5.480 2 2.9972 0.05145 . x2 0.688 1 0.7531 0.38621 x1:x2 0.367 2 0.2009 0.81809 Residuals 268.791 294 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
x1についてもx2についても、タイプI平方和の時と値が変わってますね。
Intercept（切片）項も出力されますが、無視すればいいです。分散分析は独立変数の回帰係数のところだけ見るものなので。
参考
↓の記事に、{car}パッケージの解説が少し載ってました。
Colorless Green Ideas:Rの分散分析表で逐次平方和でなく、調整平方和を出す方法
普通、調整平方和というと、Type-IIの平方和のことを指す。これは、他の主要因（一次）の影響を除いたときの平方和を求めるものである。Type-IIIの平方和は、他の主要因（一次）だけでなく二次以上の項の影響も除いて計算した平方和である。
と書かれている部分の「主要因」ってのは分散分析で言えば主効果のことで、「二次以上の項」ってのは交互作用のところですね。タイプIIは、「主効果を交互作用よりも優先投入する」という点ではタイプIの発想を用い、「複数の主効果の間では、最後に投入した時の平方和を用いる」という意味でタイプIIIの発想を用いるものです。
↓の井関龍太氏のサイトでは、「ANOVA君」という分散分析の関数が公開されていて、これは平方和のタイプを選ぶ機能がついているのですが、デフォルトを何にするかで迷ったという話が書かれており、平方和のタイプの長短についてもいろいろ論争があるということが分かり、参考になります。
ANOVA君/平方和のタイプ - 井関龍太のページ
ANOVA君は、データフレームの作り方から指定されていることもあり、私はまだ使ったことがないのですが、今度つかってみたいと思います。
なお、統計と関係ない余談なのですが、井関氏の心理学の論文（物語を読解するときの心理状態を研究したもの）を１本読んだことがあり、ものすごく参考になりました。
