ANOVA君/平方和のタイプ/タイプⅢ平方和の計算法 - 井関龍太のページ

ホーム   編集 凍結 差分 バックアップ 添付 コピー 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS

ANOVA君/平方和のタイプ/タイプⅢ平方和の計算法

Last-modified: 2014年04月21日 (月) 20:33:29 (1666d)
Top > ANOVA君 > 平方和のタイプ > タイプⅢ平方和の計算法

タイプⅢ平方和の計算法

引き続き,タイプⅡと同じデータ例を使って,今度はタイプⅢの計算法について説明します。

     s  A  B    y
1   s1 a1 b1  9.8
2   s2 a1 b1 12.8
3   s3 a1 b1 11.2
4   s4 a1 b2 20.4
5   s5 a1 b2 14.7
6   s6 a1 b2 18.1
7   s7 a2 b1 26.7
8   s8 a2 b1 26.9
9   s9 a2 b2 19.5
10 s10 a2 b2 28.9
11 s11 a3 b1 27.3
12 s12 a3 b1 36.9
13 s13 a3 b2 25.9

平方和を計算するまでの手順は,基本的にこれまでと同じなのですが,1つ注意する点があります。
それは,計画行列の作り方です。
これまでは,特に気にせずデフォルトの方法を使ってきました。
しかし,Rにおける計画行列の作り方には,5つのオプションがあります。
このオプションは,以下の2つの方法で切り替えることができます。

デフォルトでは,contr.treatmentという方法が使われています。
これを例えば,contr.sumというオプションに変えるには,

> model.matrix(~ A + B + A:B, dat, contrasts = list(A = "contr.sum", B = "contr.sum"))

のように,各要因にどのようなcontrastsをしてほしいかを指定します。
交互作用については,要因のかけ算で作られるので指定する必要はありません。

もう1つの方法は,

> options(contrasts = c("contr.sum", "contr.sum"))

として,オプションの設定を変更することです。
このようにしておくと,以後,関数の使用時にいちいち指定しなくても,常にcontr.sumが使われます。
もとの設定を残したい場合は,

> options()$contrasts

と入力すると,その時点でのcontrastsの設定が表示されるので,これを適当な変数に代入しておいて,作業が終わった後で設定を戻すとよいと思われます。

それでは,5つのcontrastsオプションの違いを見てみましょう。

#contr.treatment
> model.matrix(~ A + B + A:B, dat)
   (Intercept) Aa2 Aa3 Bb2 Aa2:Bb2 Aa3:Bb2
1            1   0   0   0       0       0
2            1   0   0   0       0       0
3            1   0   0   0       0       0
4            1   0   0   1       0       0
5            1   0   0   1       0       0
6            1   0   0   1       0       0
7            1   1   0   0       0       0
8            1   1   0   0       0       0
9            1   1   0   1       1       0
10           1   1   0   1       1       0
11           1   0   1   0       0       0
12           1   0   1   0       0       0
13           1   0   1   1       0       1
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$A
[1] "contr.treatment"

attr(,"contrasts")$B
[1] "contr.treatment"
#contr.helmert
> model.matrix(~ A + B + A:B, dat)
   (Intercept) A1 A2 B1 A1:B1 A2:B1
1            1 -1 -1 -1     1     1
2            1 -1 -1 -1     1     1
3            1 -1 -1 -1     1     1
4            1 -1 -1  1    -1    -1
5            1 -1 -1  1    -1    -1
6            1 -1 -1  1    -1    -1
7            1  1 -1 -1    -1     1
8            1  1 -1 -1    -1     1
9            1  1 -1  1     1    -1
10           1  1 -1  1     1    -1
11           1  0  2 -1     0    -2
12           1  0  2 -1     0    -2
13           1  0  2  1     0     2
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$A
[1] "contr.helmert"

attr(,"contrasts")$B
[1] "contr.helmert"
#contr.sum
> model.matrix(~ A + B + A:B, dat)
   (Intercept) A1 A2 B1 A1:B1 A2:B1
1            1  1  0  1     1     0
2            1  1  0  1     1     0
3            1  1  0  1     1     0
4            1  1  0 -1    -1     0
5            1  1  0 -1    -1     0
6            1  1  0 -1    -1     0
7            1  0  1  1     0     1
8            1  0  1  1     0     1
9            1  0  1 -1     0    -1
10           1  0  1 -1     0    -1
11           1 -1 -1  1    -1    -1
12           1 -1 -1  1    -1    -1
13           1 -1 -1 -1     1     1
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$A
[1] "contr.sum"

attr(,"contrasts")$B
[1] "contr.sum"
#contr.poly
> model.matrix(~ A + B + A:B, dat)
   (Intercept)           A.L        A.Q        B.L       A.L:B.L    A.Q:B.L
1            1 -7.071068e-01  0.4082483 -0.7071068  5.000000e-01 -0.2886751
2            1 -7.071068e-01  0.4082483 -0.7071068  5.000000e-01 -0.2886751
3            1 -7.071068e-01  0.4082483 -0.7071068  5.000000e-01 -0.2886751
4            1 -7.071068e-01  0.4082483  0.7071068 -5.000000e-01  0.2886751
5            1 -7.071068e-01  0.4082483  0.7071068 -5.000000e-01  0.2886751
6            1 -7.071068e-01  0.4082483  0.7071068 -5.000000e-01  0.2886751
7            1 -9.073264e-17 -0.8164966 -0.7071068  6.415766e-17  0.5773503
8            1 -9.073264e-17 -0.8164966 -0.7071068  6.415766e-17  0.5773503
9            1 -9.073264e-17 -0.8164966  0.7071068 -6.415766e-17 -0.5773503
10           1 -9.073264e-17 -0.8164966  0.7071068 -6.415766e-17 -0.5773503
11           1  7.071068e-01  0.4082483 -0.7071068 -5.000000e-01 -0.2886751
12           1  7.071068e-01  0.4082483 -0.7071068 -5.000000e-01 -0.2886751
13           1  7.071068e-01  0.4082483  0.7071068  5.000000e-01  0.2886751
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$A
[1] "contr.poly"

attr(,"contrasts")$B
[1] "contr.poly"
#contr.SAS
> model.matrix(~ A + B + A:B, dat)
   (Intercept) Aa1 Aa2 Bb1 Aa1:Bb1 Aa2:Bb1
1            1   1   0   1       1       0
2            1   1   0   1       1       0
3            1   1   0   1       1       0
4            1   1   0   0       0       0
5            1   1   0   0       0       0
6            1   1   0   0       0       0
7            1   0   1   1       0       1
8            1   0   1   1       0       1
9            1   0   1   0       0       0
10           1   0   1   0       0       0
11           1   0   0   1       0       0
12           1   0   0   1       0       0
13           1   0   0   0       0       0
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$A
[1] "contr.SAS"

attr(,"contrasts")$B
[1] "contr.SAS"

reducedされているので少しわかりにくいですが,contr.helmert,contr.sum,contr.polyは,各要因内のすべての水準の列を足すと0になるように計画行列を構成しています。

タイプⅢの平方和は,これらの,列を足して0になる計画行列を使った場合にのみ,適切に計算できます。
他の計画行列を使った場合には,数値は計算できますが,おそらく,あまり意味のない値が出てきます。
つまり,タイプⅢ平方和では,使用する計画行列によって計算結果が変わるのです。
タイプⅠやタイプⅡでは,基本的に,このようなことはありません。
これが「タイプⅢ平方和は特定の仮定に依存する」ということの意味だと思います。
以下では,contr.sumによる計画行列を使って,タイプⅢ平方和を計算してみます。

contrastsオプションの指定を除くと,積和行列,rss1まではこれまでと同じです。

> dat <- read.table("clipboard", header = T)
> options(contrasts = c("contr.sum", "contr.sum"))
> desmat <- model.matrix(~ A + B + A:B, dat)
> exmat <- cbind(desmat, dat$y)
> promat <- crossprod(exmat)
> promat
            (Intercept)   A1   A2   B1 A1:B1 A2:B1        
(Intercept)        13.0  3.0  1.0  1.0  -1.0  -1.0  279.10
A1                  3.0  9.0  3.0 -1.0   1.0   1.0   -3.10
A2                  1.0  3.0  7.0 -1.0   1.0   1.0   11.90
B1                  1.0 -1.0 -1.0 13.0   3.0   1.0   24.10
A1:B1              -1.0  1.0  1.0  3.0   9.0   3.0  -57.70
A2:B1              -1.0  1.0  1.0  1.0   3.0   7.0  -33.10
                  279.1 -3.1 11.9 24.1 -57.7 -33.1 6774.85 
> rss1 <- promat[7, 7] - lm.fit(promat[1, 1, drop = FALSE], promat[1, 7])$coef * promat[1, 7]
> rss1
(Intercept) 
782.7877 

タイプⅢ平方和では,検討したい効果以外のすべての効果を含んだモデルとそれに検討したい効果も含めたモデル(つまり,すべての効果を含めたモデル)の差分を取ります。

そこで,各平方和は以下のように求められます。

要因A=R(A|μ,B,A×B)
要因B=R(B|μ,A,A×B)
要因A×B=R(A×B|μ,A,B)

具体的には,次のようになります。

> rss2 <- promat[7, 7] - sum(lm.fit(promat[c(1,4:6), c(1,4:6)], promat[c(1,4:6), 7])$coef * promat[c(1,4:6), 7])
> rss0 <- promat[7, 7] - sum(lm.fit(promat[1:6, 1:6], promat[1:6, 7])$coef * promat[1:6, 7])
> rss2 - rss0
[1] 505.1579

これが要因AのタイプⅢ平方和になります。
ここで,rss2は,全効果から要因Aを除いた場合の残差平方和,rss0は全効果を投入した場合の残差平方和です。

同様に,要因Bと要因A×Bについては以下のようになります。

> rss3 <- promat[7, 7] - sum(lm.fit(promat[c(1:3,5:6), c(1:3,5:6)], promat[c(1:3,5:6), 7])$coef * promat[c(1:3,5:6), 7])
> rss4 <- promat[7, 7] - sum(lm.fit(promat[1:4, 1:4], promat[1:4, 7])$coef * promat[1:4, 7])
> rss3 - rss0#要因Bの平方和
[1] 1.719298
> rss4 - rss0#要因A×Bの平方和
[1] 92.33404

全効果を投入したときの残差平方和は既に計算したので,これを再び用いています。
ANOVA君のアルゴリズムにおいて,タイプⅢ平方和の方がタイプⅡ平方和よりも計算量が少なくすむのは,このためです。
(タイプⅡのときは,引き算する残差平方和は必ずしも一貫しないので,そのつど適切なものを計算する必要があります。)

誤差平方和は,rss0と一致します。

> rss0
[1] 111.2333

まとめると,タイプⅢ平方和は以下のようになります。

要因A→Bの順B→Aの順
A505.16505.16
B1.721.72
AxB92.3392.33
Error111.23111.23
Total782.79782.79

投入順序を変えても結果は同じになっています。

以上のように,タイプⅢ平方和では,自身以外のすべての効果の影響を取り除いた上で,その効果によって説明される分散を計算しています。
一方,タイプⅡ平方和では,自身を含む上位の効果は予めモデルから除いてあります。
このため,一般には,(下位の効果については)タイプⅡ平方和の方が検出力が高くなります。

このことの当否を問うのはむずかしい問題です。
交互作用は,実際に存在する効果の組み合わせによって生じています。
そこで,交互作用があるのに,そこに含まれる要因の効果がまったくないと考えるのはおかしい,といった考え方ができます(Langsrud, 2003)。
その要因は,交互作用という形で(遠まわしに)影響を与えているではないか,ということです。
すると,主効果を検討するときにはその効果を含む交互作用はモデルに含めず,交互作用を検討するときには下位の効果の影響は考慮するという意味で,タイプⅡが適切であるという議論ができます。
タイプⅢでは,主効果の平方和を計算するときに,その主効果は含めないが上位の交互作用は含むモデルを使っています。

一方,主効果は主効果で交互作用とは別物であると考えるなら,主効果(や下位の効果)を検討するときに,他のすべての効果からの影響を統制することも適切であるかもしれません。
そうなると,タイプⅡよりもタイプⅢの方が妥当であるといえるかもしれません。

また,タイプⅠとの仮定の違い,セル内のデータ数の影響の有無の問題もあり,タイプⅡとタイプⅢのどちらが適切かを決めることは容易ではなくなっています。

Rの関数では・・・

タイプⅡの場合と同じく,carパッケージのAnova関数を使うことでタイプⅢ平方和を計算できます。
Anova関数では,タイプⅢを計算するためのオプションがあります。
これは,以下のように指定します。

> Anova(lm(y ~ A + B + A:B, dat), type = "III")

ただし,注意すべき点として,lm関数の部分を実行する際にcontrastsの設定を変更しておく必要があります。
上に述べたのと同じように,contr.sum,contr.helmert,contr.polyを使うと,結果はすべて一致し,適切なタイプⅢ平方和を得ることができます。
contr.treatment,contr.SASを使った場合には,それぞればらばらな(おそらくは)意味不明な値が出てきます(Rのデフォルトは,contr.treatmentになっていると思います)。
ANOVA君では,関数の中で一時的にcontrastsの設定をcontr.sumに変更してから計画行列を作り,作業が終わった後で設定を元に戻すようにしています。

参考文献

以上の計算手順は,主に,Milliken & Johnson(1992)に基づいています。
タイプⅡ,タイプⅢ平方和の計算法について,具体的な手順のレベルで説明している文献がなかなか見つからず苦労していたところ,この文献に出会うことができました。
積和行列を作るところまでは,高橋他(1989)にもくわしく説明されています。

Milliken & Johnson(1992)では,タイプⅢ平方和の計算法についてはそれほど具体的には述べていません。
この点については,Overall & Spiegel(1969)の記述が参考になりました。

高橋他(1989)を見ればわかるように,SASの計算法は,ここで紹介した方法とは若干異なります。
特に,タイプⅢ平方和については,ANOVA君の計算法は,SASのそれというよりは,carパッケージのAnova関数のそれに近いものになっています。
ただ,その計算結果は,複数の文献に載っているタイプⅡ,タイプⅢの計算例と一致しますし,SPSSによる出力とも一致しています(万一,丸め誤差以上の不一致が見られた場合には,お知らせください)。
SPSSは,SASともAnova関数とも別の計算法を使っている可能性があります。


Langsrud, Ø. (2003). ANOVA for unbalanced data:Use Type II instead of Type III sums of squares. Statistics and Computing, 13, 163-167.
Milliken, G. A., & Johnson, D. E. (1992). Analysis of messy data volume 1:Designed experiments. Chapman and Hall:New York.
Overall, J. E., & Spiegel, D. K. (1969). Concerning least squares analysis of experimental data. Psychological Bulletin, 72, 311-322.
高橋行雄・大橋靖雄・芳賀敏郎 (1989). SASによる実験データの解析. 東京大学出版会

TrackBack(0) | 外部リンク元 | このエントリーをはてなブックマークに追加