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

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

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

Last-modified: 2013年09月28日 (土) 22:16:13 (1871d)
Top > ANOVA君 > 平方和のタイプ > タイプⅠ平方和の計算法

各平方和の計算法

平方和のタイプの違いを具体的に説明するには,その計算法を述べることが近道であるように思われます。
しかし,純粋に数式のみを使って説明することは私には荷が重いので,Rを使って平方和を計算する手順に沿って話をしたいと思います。
ここでは,Rを使ったタイプⅠ,タイプⅡ,タイプⅢ平方和の計算法を簡単に述べることとします。
(ただし,以下で述べる手続きは,実際にANOVA君で用いているものとは若干異なります。基本的な部分は同じですが,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
14 s14 a3 b2 20.5

sは被験者のラベル,A,Bは各要因のラベル,yは従属変数を表します。
3×2の被験者間計画のデータで,非釣り合い型になっています。
(被験者間計画なのは,説明を簡単にするためです。)

このデータにdatという名前を付けてRに読み込みます(一行目はヘッダとして認識させます)。

> dat <- read.table("clipboard", header = T)

タイプⅠ平方和の計算法

タイプⅡ,タイプⅢの平方和を計算するには,一般線形モデルを用いて平方和を計算する方法を知る必要があります。
そのための基礎として,まず,最小二乗法によるタイプⅠ平方和の計算手順について述べます。

最初に,計画行列(design matrix)を作ります。
この計画行列は,重回帰分析におけるダミー変数に当たるものです。
各要因の各水準,要因間の交互作用を数値で表します。
例えば,上のデータについての計画行列は以下のように表すことができます。

  (Intercept) Aa1 Aa2 Aa3 Bb1 Bb2 Aa1:Bb1 Aa1:Bb2 Aa2:Bb1 Aa2:Bb2 Aa3:Bb1 Aa3:Bb2
1           1   1   0   0   1   0       1       0       0       0       0       0
2           1   1   0   0   1   0       1       0       0       0       0       0
3           1   1   0   0   1   0       1       0       0       0       0       0
4           1   1   0   0   0   1       0       1       0       0       0       0
5           1   1   0   0   0   1       0       1       0       0       0       0
6           1   1   0   0   0   1       0       1       0       0       0       0
7           1   0   1   0   1   0       0       0       1       0       0       0
8           1   0   1   0   1   0       0       0       1       0       0       0
9           1   0   1   0   0   1       0       0       0       1       0       0
10          1   0   1   0   0   1       0       0       0       1       0       0
11          1   0   0   1   1   0       0       0       0       0       1       0
12          1   0   0   1   1   0       0       0       0       0       1       0
13          1   0   0   1   0   1       0       0       0       0       0       1
14          1   0   0   1   0   1       0       0       0       0       0       1

各データ(y)について,それが該当する要因の該当する水準の列には1,該当しない水準の列には0を当てています。
交互作用の列は,対応する各水準の列の値同士をかけたものになっています。

ただし,重回帰分析と同じで,このままの行列で計算をしようとすると,解の出ない行が出てきます。
例えば,ある行のa1水準が0になるか1になるかは,a2水準の列とa3水準の列を見れば完全に予測できます。
このような状況では,重回帰分析でも解は出ません。

このような場合はこのような場合で対処法があるのですが,推定できない部分が出ないように最初からダミー変数の列を減らした行列を作る方法があります。
このタイプの行列はreduced matrixと呼ばれています。

Rでは,model.matrix関数を使うと簡単に計画行列を作ることができます。

> 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
14           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"

model.matrix関数には,model.matrix(モデル式,データフレーム)のように変数を投入します。
モデル式は,データフレームのヘッダとして用いた文字列を使って作ります。
「~」は“=”に当たるもので,model.matrixの場合には常に付ける必要があります。
「+」を使って変数をつなぎます。
「:」は変数間の交互作用を表します。

model.matrixで作られる計画行列は,最初からreduced matrixになっています(むしろ,この関数ではreducedさせない方が面倒です)。
reduced matrixでもreducedしていない計画行列と同じ解が出せる上に,行列のサイズが小さくなるので計算の負荷を減らすことができます。
そこで,以下では,基本的にreduced matrixの使用に基づいて話を進めます。

では,この計画行列をdesmatという名前を付けた変数に代入しておきます。

> desmat <- model.matrix(~ A + B + A:B, dat)

以下でも,同様に,「<-」を使って適当な名前の変数に中間結果を保存することがあります。
この計画行列にデータの従属変数部分(y)を取り付けて,拡大行列を作ります。

> exmat <- cbind(desmat, dat$y)
> exmat
   (Intercept) Aa2 Aa3 Bb2 Aa2:Bb2 Aa3:Bb2     
1            1   0   0   0       0       0  9.8
2            1   0   0   0       0       0 12.8
3            1   0   0   0       0       0 11.2
4            1   0   0   1       0       0 20.4
5            1   0   0   1       0       0 14.7
6            1   0   0   1       0       0 18.1
7            1   1   0   0       0       0 26.7
8            1   1   0   0       0       0 26.9
9            1   1   0   1       1       0 19.5
10           1   1   0   1       1       0 28.9
11           1   0   1   0       0       0 27.3
12           1   0   1   0       0       0 36.9
13           1   0   1   1       0       1 25.9
14           1   0   1   1       0       1 20.5

さらに,この拡大行列を転置した行列ともとの行列をかけ算します。
Rでは,行列の転置にt関数,行列のかけ算には「%*%」を用いるので,以下のようなコードになります。

> t(exmat) %*% exmat
            (Intercept) Aa2   Aa3 Bb2 Aa2:Bb2 Aa3:Bb2       
(Intercept)        14.0   4   4.0   7     2.0     2.0  299.6
Aa2                 4.0   4   0.0   2     2.0     0.0  102.0
Aa3                 4.0   0   4.0   2     0.0     2.0  110.6
Bb2                 7.0   2   2.0   7     2.0     2.0  148.0
Aa2:Bb2             2.0   2   0.0   2     2.0     0.0   48.4
Aa3:Bb2             2.0   0   2.0   2     0.0     2.0   46.4
                  299.6 102 110.6 148    48.4    46.4 7195.1 

ただし,転置行列とのかけ算は,crossprodという関数を使えば一発でできます。

> promat <- crossprod(exmat)#今度は変数に代入しています
> promat
            (Intercept) Aa2   Aa3 Bb2 Aa2:Bb2 Aa3:Bb2       
(Intercept)        14.0   4   4.0   7     2.0     2.0  299.6
Aa2                 4.0   4   0.0   2     2.0     0.0  102.0
Aa3                 4.0   0   4.0   2     0.0     2.0  110.6
Bb2                 7.0   2   2.0   7     2.0     2.0  148.0
Aa2:Bb2             2.0   2   0.0   2     2.0     0.0   48.4
Aa3:Bb2             2.0   0   2.0   2     0.0     2.0   46.4
                  299.6 102 110.6 148    48.4    46.4 7195.1

これで積和行列ができました。
この行列が平方和を算出するためのベースになります。

まず,積和行列の右下の数値に注目してください。
「7195.1」とあります。
この値は,分散の総和に当たります。
分散分析では,この分散の総和から,各要因が説明できる分散を分解していくことになります。
そのための方法の1つとして,最小二乗法を用いることができます。

では,切片によって説明できる分散だけを分解してみましょう。
積和行列(promat)を見てください。
切片(Intercept)は,一列目によって表現されています(積和行列は対角行列なので一行目と言っても同じです)。
そこで,この行列の左上から1×1の要素だけを取り出します。
[ ]を使い,数値を指定することで行列の一部だけを取り出すことができます。

> promat[1, 1]
[1] 14

この値からデータの分散を説明します。
自分で行列計算を行って解を出すこともできますが,Rでは,lm.fit関数を使うと簡単に最小二乗解が得られます。
lm.fit(計画行列, データベクトル)とすると,以下のようにずらりと解が得られます。

> lm.fit(promat[1, 1, drop = FALSE], promat[1, 7])
$coefficients
(Intercept) 
       21.4 

$residuals
[1] 0

$effects
(Intercept) 
      299.6 

$rank
[1] 1

$fitted.values
[1] 299.6

$assign
NULL

$qr
$qr
            (Intercept)
(Intercept)          14

$qraux
[1] 14

$pivot
[1] 1

$tol
[1] 1e-07

$rank
[1] 1

attr(,"class")
[1] "qr"

$df.residual
[1] 0

ここで,drop = FALSEと加えてあるのは,行列から要素を1個だけ取り出したときに,データの形式が行列でなくなることを防ぐためです。
lm.fit関数は,計画行列の部分に行列形式のデータを代入しないとうまく動作しません。
また,「promat[1, 7]」の部分は,積和行列のy変数に当たる部分(7列目)の1行目を指定しています。

出力の中で重要なのは,coefficients(係数)の部分です。
以下では,係数だけを一挙に取り出すため,次のように指定することにします。

> lm.fit(promat[1, 1, drop = FALSE], promat[1, 7])$coef
(Intercept) 
   21.4 

この値をもとのyの値をかけ,総和から引きます。

> rss1 <- promat[7, 7] - lm.fit(promat[1, 1, drop = FALSE], promat[1, 7])$coef * promat[1, 7]
> rss1
(Intercept) 
   783.66 

これが切片の残差平方和(resudual sum of square:RSS)になります。
これは,この要因(今回は切片)のみでは説明できない分散になります。

同様に,要因Aの残差平方和を計算します。
要因Aは積和行列の2,3列目で表現されます。
そこで,切片の1列目を加えて,1~3列目を指定します(1:3)。
今回は,最小二乗解を出すと,3つの値が得られます。

> lm.fit(promat[1:3, 1:3], promat[1:3, 7])$coef
(Intercept)   Aa2      Aa3 
   14.50    11.00    13.15 

この値は,それぞれ7列目の1~3行番目に相当する値とかけてから合計します。

> sum(lm.fit(promat[1:3, 1:3], promat[1:3, 7])$coef * promat[1:3, 7])
[1] 6920.59

この結果をやはり総和から引くことで要因Aのための残差平方和が得られます。

> rss2 <- promat[7, 7] - sum(lm.fit(promat[1:3, 1:3], promat[1:3, 7])$coef * promat[1:3, 7])
> rss2
[1] 274.51

(ちなみに,reducedしていない計画行列を用いた場合,A要因の3つの水準のうちどこか1つは解が出ないことになります。
解が出なかった部分には「NA」と表示されます。
この場合は,「NA」を0に置き換えて上と同様の計算をすれば,reduced matrixを使った場合と同じ解が得られます。)

さて,それでは,上の2つの残差平方和から要因Aの効果を計算してみます。
rss1は切片のみ,rss2は切片+要因Aの効果で説明できない分散を推定したものですから,この2つのRSSを引けば,要因A単独で説明される部分が得られるはずです。
この関係性はR(A|μ)のように表現できます。
この式は,μ(切片)を統制してAの効果を推定する,というような意味です。
実際に計算すると,以下の値が得られます。

> rss1 - rss2
(Intercept) 
   509.15 

これが要因AのタイプⅠ平方和になります。

同様に,要因Bと交互作用A×Bのための残差平方和も計算してみます。

> rss3 <- promat[7, 7] - sum(lm.fit(promat[1:4, 1:4], promat[1:4, 7])$coef * promat[1:4, 7])
> rss3
[1] 273.5843
> rss4 <- promat[7, 7] - sum(lm.fit(promat[1:6, 1:6], promat[1:6, 7])$coef * promat[1:6, 7])
> rss4
[1] 125.8133

要因Bのときは切片,要因Aの列を含み,交互作用A×Bのときはすべての効果の列を含むことに注意してください。
ここで,R(B|μ,A)はrss2 - rss3,R(A×B|μ,A,B)はrss3 - rss4によって得られることになります。

> rss2 - rss3#要因Bの平方和
[1] 0.9257143
> rss3 - rss4#要因A×Bの平方和
[1] 147.7710

さらに,誤差平方和を計算する必要があります。
誤差平方和については,積和行列のy変数に対応する部分を除いたすべての行列を用いて残差平方和を計算することで得られます。
実は,これはrss4に相当するので,rss4の値をそのまま使えばOKです。

これで,このデータについて必要なすべての平方和が得られました。
値を表にまとめました。

A509.15
B0.9257143
AxB147.7710
Error125.8133
Total783.66

これが授業などで習う,手計算による方法と同じ値の平方和,すなわち,タイプⅠ平方和になります。
タイプⅠ平方和は,この例のような釣り合い型データに対しては適切です。

しかし,非釣り合い型の場合には,問題が起こります。
仮に例として挙げたデータをs1~s13までしかないものと見なして,13名のデータでタイプⅠ平方和を計算したとします。
そうすると,変数の投入順序によって結果が変わるのです。
投入順によって,以下のように,要因A,Bの平方和の値が違ってきます。

要因A→Bの順B→Aの順
A576.44578.68
B2.780.54
AxB92.3392.33
Error111.23111.23
Total783.66783.66

この場合,どちらの順序による結果が適切であるかを決めることはできません。
釣り合い型の場合には,どちらの順序でも結果は同じになります。
このようなわけで,非釣り合い型計画でも投入順序によって結果が変わらない方法として,タイプⅡやタイプⅢの平方和が考案されたというわけです。

Rの関数では・・・

もともとRに入っている関数でタイプⅠ平方和による分散分析を行うには,まず,anova関数を用いる方法があります。
この関数は,lm関数の計算結果を引数として取ります。
例えば,以下のようにすれば,簡単に分散分析表を得ることができます。

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

ただし,この関数では,被験者内要因を含む分析を行うのが少々面倒になります。
anova関数では,F比を計算するときの誤差項を指定できないようなのです。
そのため,要因Aの平均平方をs×Aの平均平方で割ってほしいのに,両方ともs×A×Bに相当する平均平方で割られてしまうということが起こります。

そこで,より使い勝手のよい関数としてaov関数がお勧めです。
この関数は,lm関数を呼び出す必要がなく,誤差項を指定できます。
例えば,2要因の混合要因計画(AsBデザイン)であば,以下のように指定すれば,適切な結果を出力してくれます。

> aov(y ~ A + B + A:B + Error(s:A + s:A:B), dat)

問題は,どちらの関数も,タイプⅠ以外の結果を計算してくれないことです。
ANOVA君は,もともと,js-STARやSPSS形式のデータを整形してaov関数に投入する関数として構想されていたのですが,この点に不満があったことから迷走が始まりました。

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