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

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

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

Last-modified: 2013年09月28日 (土) 22:15:31 (1871d)
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

積和行列を作るところまでは,タイプⅠの計算と同じです。

> dat <- read.table("clipboard", T)
> desmat <- model.matrix(~ A + B + A:B, dat)
> exmat <- cbind(desmat, dat$y)
> promat <- crossprod(exmat)
> promat
            (Intercept) Aa2  Aa3   Bb2 Aa2:Bb2 Aa3:Bb2        
(Intercept)        13.0   4  3.0   6.0     2.0     1.0  279.10
Aa2                 4.0   4  0.0   2.0     2.0     0.0  102.00
Aa3                 3.0   0  3.0   1.0     0.0     1.0   90.10
Bb2                 6.0   2  1.0   6.0     2.0     1.0  127.50
Aa2:Bb2             2.0   2  0.0   2.0     2.0     0.0   48.40
Aa3:Bb2             1.0   0  1.0   1.0     0.0     1.0   25.90
                  279.1 102 90.1 127.5    48.4    25.9 6774.85

rss1についても,計算法は同じになります。

> 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)
要因B=R(B|μ,A)
要因A×B=R(A×B|μ,A,B)

まず,要因Aの平方和について計算してみます。

> rss2a <- promat[7, 7] - sum(lm.fit(promat[c(1,4), c(1,4)], promat[c(1,4), 7])$coef * promat[c(1,4), 7])
> rss2b <- promat[7, 7] - sum(lm.fit(promat[1:4, 1:4], promat[1:4, 7])$coef * promat[1:4, 7])
> rss2a - rss2b
[1] 578.6848

rss2aは,Aとその上位の効果(すなわち,A×B)を除いたモデル,すなわち,μ(1列目),B(4列目)のみを含むモデルの残差平方和です。
c(1, 4)とすることで,1列目と4列目をベクトル化して,2,3列目を含めないようにしています。
一方,rss2bは,このモデルにAの列のみ加えたものです。
これらの残差平方和の差分が要因AのタイプⅡ平方和です。

同様に,要因Bと要因A×Bの平方和も計算してみます。

#要因Bの平方和
> rss3a <- promat[7, 7] - sum(lm.fit(promat[1:3, 1:3], promat[1:3, 7])$coef * promat[1:3, 7])
> rss3b <- promat[7, 7] - sum(lm.fit(promat[1:4, 1:4], promat[1:4, 7])$coef * promat[1:4, 7])
> rss3a - rss3b
[1] 2.779298

#要因A×Bの平方和
> rss4a <- promat[7, 7] - sum(lm.fit(promat[1:4, 1:4], promat[1:4, 7])$coef * promat[1:4, 7])
> rss4b <- promat[7, 7] - sum(lm.fit(promat[1:6, 1:6], promat[1:6, 7])$coef * promat[1:6, 7])
> rss4a - rss4b
[1] 92.33404

誤差平方和は,rss4bと同じです。

> rss4b
[1] 111.2333

まとめると,タイプⅡ平方和の計算結果は以下のようになります。

要因A→Bの順B→Aの順
A578.68578.68
B2.782.78
AxB92.3392.33
Error111.23111.23
Total782.79782.79

右側には,要因Bを先に投入した場合の結果を示しています。
どちらの計算結果もまったく同じになります。

このように見ていくとわかるように,タイプⅡ平方和は,タイプⅠ平方和の計算手順から,各残差平方和の計算の際に投入する変数の種類と引き算する相手を変更したものであることがわかります。

注意したいのは,検討したい要因のすべての上位の項を除く必要があるということです。
例えば,3要因の分散分析(ABCsデザイン)では,各平方和の算出法は以下のようになります。

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

要因Aの場合には,Aを含むすべての効果,すなわち,A×B,A×C,A×B×Cのすべての効果を統制する必要があります。
ANOVA君のファイルには,この作業を行うための関数としてelematch関数を用意しています。

Rの関数では・・・

Rの関数でも,タイプⅡ平方和を計算することができます。
まずは,drop1関数を使う方法があります。
この関数を使う場合,1つ1つの効果についての平方和を個別に計算する必要があります(場合によっては,一回の適用で複数の適切な平方和が得られることもありますが,高次の要因計画では,そのような場合は少なくなります)。
例えば,要因AのタイプⅡ平方和は以下のように指定することで得られます。

> drop1(lm(y ~ A + B, dat), "A")

このとき投入するモデルは,上で述べたのと同じように,上位の効果を含まないモデルにする必要があります。
この方法では,効果の数だけモデルを変えつつdrop1関数を適用する必要があるので,高次の要因計画では非常に時間がかかります。
anovakun 2.0.0はこの方法を用いていたために計算時間が長くなっていました。

drop1関数よりも簡単にタイプⅡ平方和を計算する方法もあります。
そのためには,「car」パッケージをインストールする必要があります。
carパッケージを読み込んでいれば,「Anova」という関数が使えます。
ややこしいですが,これはanova関数とは別物です。
Rは大文字と小文字を別のものとして区別しています。
これは変数についても同様です。

Anova関数は,anova関数と同様に,lm関数の結果を引数として取ります。

> Anova(lm(y ~ A + B + A:B, dat))

これで簡単にタイプⅡ平方和を計算できます。

このAnova関数は,anova関数とは違い,うまく使うと誤差項を指定でき(被験者内要因を含むデータを適切に分析でき),Greenhouse-GeisserとHuynh-Feldtのεによって調整した結果も同時に得ることができるようです(ヘルプを参照してください)。

ただ,私には今ひとつこの関数の使い方がよくわからず,うまく被験者内要因を指定できません……。
また,計算結果を事後の分析に利用する方法もよくわかりません。
carパッケージをインストールしておく必要があるというのも,初心者には不親切かもしれません。
そのようなもろもろの理由から,平方和も自分で計算してしまえばいいのでは?と思ってしまったために,ANOVA君は今のような形になりました。

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