Rで手を動かす分散分析（1元配置）

はじめに
①Rのaov関数を用いて簡単に1元配置分散分析を行う方法
②Σの計算を自分でちゃんとやる方法
③行列を使った方法

はじめに

本記事では3つの分散分析を行う方法についてまとめます。

①aov関数を用いた簡単な方法
②Σの計算を自分でちゃんとやる方法
③行列を使った方法

実際に手を動かすことで、分散分析の理解につながれば幸いです。すべての方法で、図1の実験結果を用います。要因がAのみで4水準(1,2,3,4)、実験結果がyです。

①Rのaov関数を用いて簡単に1元配置分散分析を行う方法

#まず、データを用意して...
a <- factor(c(rep(1,5),rep(2,5),rep(3,5),rep(4,5)))
y <- c(575,542,530,539,570,
       565,593,590,579,610,
       600,651,610,637,629,
       725,700,715,685,710)

#aov関数で分散分析を実行
df <- data.frame(a,y)
summary(aov(y~a, df))

#実行結果
> summary(aov(y~a, df))
            Df Sum Sq Mean Sq F value   Pr(>F)    
a            3  66871   22290    66.8 2.88e-09 ***
Residuals   16   5339     334                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>

以上です。実行結果で、aに***がついていて有意であることが分かります。

②Σの計算を自分でちゃんとやる方法

次に、教科書に書いてある公式に従って、Σの計算をちゃんとやってみようと思います（図1）。

図1で、今回は水準数が4なのでa=4, 水準ごとの繰返しは5回なので、n=5になります。

図1の意味の解説は教科書に譲るとして、実際にRで計算してみます。

#Σの形で考える
ST <- sum((df$y-mean(df$y))^2) #STの計算
yA <- by(df$y, INDICES = df$a, FUN = mean)#aごとの平均を計算
n <- 5#繰り返し
SA <- n*sum((yA-mean(df$y))^2)#SAの計算
SE <- ST-SA#SEの計算

#実行結果（平方和はaov関数を用いた結果とあってますね）
> ST
[1] 72209.75
> SA
[1] 66870.55
> SE
[1] 5339.2
> 

#平均平方の計算
MSA <- SA/3
MSE <- SE/(20-1-3)

#実行結果（平均平方もあってます）
> MSA
[1] 22290.18
> MSE
[1] 333.7
> 
#F値の計算 (これも上と一致)
> MSA/MSE
[1] 66.79707 
> 

#p値: P(F>F0)
> pf(MSA/MSE, 3, 20-1-3, lower.tail = F)
[1] 2.882866e-09
> 
#p値<0.05より有意水準5%で主効果Aは有意

平方和、平均平方、Ｆ値、p値、すべてaov関数を使った計算と一致しました。

③行列を使った方法

最後に、行列を用いた方法を説明します。あまり教科書には載ってないですが、この方法で分散分析を考えると、もう少し分かった気になります。

まず図1を今回の分散分析のモデルを数式化すると、

yij=µ+ai+εij

という形になります。ir=μ+ai+εij

yir=μ+ai+εij

ir=μ+ai

+εij

それで、図2のように行列を定義すると、下のような表し方が出来ます。

そしてここで、射影行列という行列を持ち出します。計画行列Xの射影行列Pは次のようにして計算されます。

このPを用いて、射影行列の性質から

という関係が成り立ちます。Iは単位行列です。これにYの転置とYをかけて

という式が成立するんですね。これはどこかでみた形でST=SA+SRの平方和分解に相当します。第一項めは平均値になるので、左辺に移行して

によって平方和を計算します。

以下コードです。

#射影行列を返す関数
P <- function(M){
  M <- as.matrix(M)
  invMtM <- solve(t(M)%*%M)
  projection <- M%*%invMtM%*%t(M)
  return(projection)
}

#次のXを用意
> X
   a1 a2 a3 a4 r1 r2 r3 r4 r5   y
1   1  0  0  0  1  0  0  0  0 575
2   1  0  0  0  0  1  0  0  0 542
3   1  0  0  0  0  0  1  0  0 530
4   1  0  0  0  0  0  0  1  0 539
5   1  0  0  0  0  0  0  0  1 570
6   0  1  0  0  1  0  0  0  0 565
7   0  1  0  0  0  1  0  0  0 593
8   0  1  0  0  0  0  1  0  0 590
9   0  1  0  0  0  0  0  1  0 579
10  0  1  0  0  0  0  0  0  1 610
11  0  0  1  0  1  0  0  0  0 600
12  0  0  1  0  0  1  0  0  0 651
13  0  0  1  0  0  0  1  0  0 610
14  0  0  1  0  0  0  0  1  0 637
15  0  0  1  0  0  0  0  0  1 629
16  0  0  0  1  1  0  0  0  0 725
17  0  0  0  1  0  1  0  0  0 700
18  0  0  0  1  0  0  1  0  0 715
19  0  0  0  1  0  0  0  1  0 685
20  0  0  0  1  0  0  0  0  1 710

#平方和の計算
ST <- t(X$y)%*%X$y
Xa <- X[1:4]
X$y <- X$y-mean(X$y) 
SA<- t(X$y)%*%P(Xa)%*%X$y
SE <- ST-SA

#実行結果
> ST
         [,1]
[1,] 72209.75
> SA
         [,1]
[1,] 66870.55
> SE
       [,1]
[1,] 5339.2