R言語を使いWeb教材「ハンバーガーショップで学ぶ楽しい統計学」1の「分散分析（1要因）」について進めていきます（第6回）。

分散分析をしよう
それぞれのお店のポテトの評価データを入力する。
> ワクワク<-c(80,75,80,90,95,80,80,85,85,80,90,80,75,90,85,85,90,90,85,80)

> モグモグ<-c(75,70,80,85,90,75,85,80,80,75,80,75,70,85,80,75,80,80,90,80)

> パクパク<-c(80,80,80,90,95,85,95,90,85,90,95,85,98,95,85,85,90,90,85,85)

> 全体<-c(ワクワク,モグモグ,パクパク)
> お店文字<-c(rep("ワクワク",length(ワクワク)),rep("モグモグ",length(モグモグ)),rep("パクパク",length(パクパク)))
> お店<-factor(お店文字)
> お店

> 全体データ<-cbind(ワクワク,モグモグ,パクパク) #3つの群( ワクワク,モグモグ,パクパク )のデータを横方向につないで行列をつくる
> 全体データ
ワクワク モグモグ パクパク
[1,] 80 75 80
[2,] 75 70 80
[3,] 80 80 80
[4,] 90 85 90
[5,] 95 90 95
[6,] 80 75 85
[7,] 80 85 95
[8,] 85 80 90
[9,] 85 80 85
[10,] 80 75 90
[11,] 90 80 95
[12,] 80 75 85
[13,] 75 70 98
[14,] 90 85 95
[15,] 85 80 85
[16,] 85 75 85
[17,] 90 80 90
[18,] 90 80 90
[19,] 85 90 85
[20,] 80 80 85
まずは平均を求める。
> 群平均<-colMeans(全体データ) #colMeans()は列ごとの平均を求める関数
> 群平均
ワクワク モグモグ パクパク
84.00 79.50 88.15
> 全平均<-mean(全体)
> 全平均
[1] 83.88333
次に標準偏差を求めてみる。
> sqrt(mean((ワクワク-mean(ワクワク))^2))
[1] 5.385165
> sqrt(mean((モグモグ-mean(モグモグ))^2))
[1] 5.454356
> sqrt(mean((パクパク-mean(パクパク))^2))
[1] 5.341114
> sd(全体)*sqrt((length(全体)-1)/length(全体))
[1] 6.447459

帰無仮説、対立仮説を立てよう
帰無仮説H0: 3つのお店のポテトの評価（母集団）の平均のどの組み合わせにおいても差はない。
対立仮説H1: 3つのお店のポテトの評価（母集団）の平均の少なくともひとつの組み合わせに差がある。

平均からのずれ（分散）を分析するから分散分析
データは、全体の平均からのズレ＝群間のズレ＋群内のズレ（残差）を考える必要がある。

群間のずれと郡内のずれを比較する
もし、群内のズレに比べて、群間のズレが大きければ、「母集団の平均に差がない」という帰無仮説を棄却することになる。
逆に、群内のズレに比べて、群間のズレが小さければ、「母集団の平均に差がない」という帰無仮説を採択することになる。

ずれの平方和を計算する
全体のずれ、群間のずれ、群内のずれを計算する。
ずれとは、各データについて平均からの差を２乗して足したもの、つまり「平方和」となる。
したがって 全体の平方和、群間の平方和、群内の平方和を求める。
まずは全体の平方和を求める。
> 全体平均行列<-matrix(rep(全平均,60),nrow=20,ncol=3)
> 全体平均行列
[,1] [,2] [,3]
[1,] 83.88333 83.88333 83.88333
[2,] 83.88333 83.88333 83.88333
[3,] 83.88333 83.88333 83.88333
[4,] 83.88333 83.88333 83.88333
[5,] 83.88333 83.88333 83.88333
[6,] 83.88333 83.88333 83.88333
[7,] 83.88333 83.88333 83.88333
[8,] 83.88333 83.88333 83.88333
[9,] 83.88333 83.88333 83.88333
[10,] 83.88333 83.88333 83.88333
[11,] 83.88333 83.88333 83.88333
[12,] 83.88333 83.88333 83.88333
[13,] 83.88333 83.88333 83.88333
[14,] 83.88333 83.88333 83.88333
[15,] 83.88333 83.88333 83.88333
[16,] 83.88333 83.88333 83.88333
[17,] 83.88333 83.88333 83.88333
[18,] 83.88333 83.88333 83.88333
[19,] 83.88333 83.88333 83.88333
[20,] 83.88333 83.88333 83.88333
> 全体偏差<-全体-全体平均行列
> 全体偏差
[,1] [,2] [,3]
[1,] -3.883333 -8.883333 -3.883333
[2,] -8.883333 -13.883333 -3.883333
[3,] -3.883333 -3.883333 -3.883333
[4,] 6.116667 1.116667 6.116667
[5,] 11.116667 6.116667 11.116667
[6,] -3.883333 -8.883333 1.116667
[7,] -3.883333 1.116667 11.116667
[8,] 1.116667 -3.883333 6.116667
[9,] 1.116667 -3.883333 1.116667
[10,] -3.883333 -8.883333 6.116667
[11,] 6.116667 -3.883333 11.116667
[12,] -3.883333 -8.883333 1.116667
[13,] -8.883333 -13.883333 14.116667
[14,] 6.116667 1.116667 11.116667
[15,] 1.116667 -3.883333 1.116667
[16,] 1.116667 -8.883333 1.116667
[17,] 6.116667 -3.883333 6.116667
[18,] 6.116667 -3.883333 6.116667
[19,] 1.116667 6.116667 1.116667
[20,] -3.883333 -3.883333 1.116667
> 全体平方和<-sum(全体^2)
> 全体平方和
[1] 2494.183
次に群内平方和を求める。
> 群平均行列<-matrix(rep(群平均,20),nrow=20,ncol=3,byrow=TRUE)
> 群平均行列
[,1] [,2] [,3]
[1,] 84 79.5 88.15
[2,] 84 79.5 88.15
[3,] 84 79.5 88.15
[4,] 84 79.5 88.15
[5,] 84 79.5 88.15
[6,] 84 79.5 88.15
[7,] 84 79.5 88.15
[8,] 84 79.5 88.15
[9,] 84 79.5 88.15
[10,] 84 79.5 88.15
[11,] 84 79.5 88.15
[12,] 84 79.5 88.15
[13,] 84 79.5 88.15
[14,] 84 79.5 88.15
[15,] 84 79.5 88.15
[16,] 84 79.5 88.15
[17,] 84 79.5 88.15
[18,] 84 79.5 88.15
[19,] 84 79.5 88.15
[20,] 84 79.5 88.15
> 郡内<-全体-群平均行列
> 郡内
[,1] [,2] [,3]
[1,] -4 -4.5 -8.15
[2,] -9 -9.5 -8.15
[3,] -4 0.5 -8.15
[4,] 6 5.5 1.85
[5,] 11 10.5 6.85
[6,] -4 -4.5 -3.15
[7,] -4 5.5 6.85
[8,] 1 0.5 1.85
[9,] 1 0.5 -3.15
[10,] -4 -4.5 1.85
[11,] 6 0.5 6.85
[12,] -4 -4.5 -3.15
[13,] -9 -9.5 9.85
[14,] 6 5.5 6.85
[15,] 1 0.5 -3.15
[16,] 1 -4.5 -3.15
[17,] 6 0.5 1.85
[18,] 6 0.5 1.85
[19,] 1 10.5 -3.15
[20,] -4 0.5 -3.15
> 郡内平方和<-sum(郡内^2)
> 郡内平方和
[1] 1745.55
群内平方和は、不偏分散を使い各群に分けて計算しても求まる。
> ワクワク平方和<-var(ワクワク)*((length(ワクワク)-1))
> モグモグ平方和<-var(モグモグ)*((length(モグモグ)-1))
> パクパク平方和<-var(パクパク)*((length(パクパク)-1))
> 群内平方和<-ワクワク平方和+モグモグ平方和+パクパク平方和
> 群内平方和
[1] 1745.55
最後に群間平方和をもとめる。
> 群間<-群平均行列-全体平均行列
> 群間
[,1] [,2] [,3]
[1,] 0.1166667 -4.383333 4.266667
[2,] 0.1166667 -4.383333 4.266667
[3,] 0.1166667 -4.383333 4.266667
[4,] 0.1166667 -4.383333 4.266667
[5,] 0.1166667 -4.383333 4.266667
[6,] 0.1166667 -4.383333 4.266667
[7,] 0.1166667 -4.383333 4.266667
[8,] 0.1166667 -4.383333 4.266667
[9,] 0.1166667 -4.383333 4.266667
[10,] 0.1166667 -4.383333 4.266667
[11,] 0.1166667 -4.383333 4.266667
[12,] 0.1166667 -4.383333 4.266667
[13,] 0.1166667 -4.383333 4.266667
[14,] 0.1166667 -4.383333 4.266667
[15,] 0.1166667 -4.383333 4.266667
[16,] 0.1166667 -4.383333 4.266667
[17,] 0.1166667 -4.383333 4.266667
[18,] 0.1166667 -4.383333 4.266667
[19,] 0.1166667 -4.383333 4.266667
[20,] 0.1166667 -4.383333 4.266667
> 群間平方和<-sum(群間^2)
> 群間平方和
[1] 748.6333
群間の平方和は、各群ごとに値が同じになるので群内平均と全体平均の差の２乗に各群のデータの個数をかけたものを合計しても求まる。
> ワクワク平方和<-((mean(ワクワク)-全平均)^2)*length(ワクワク)
> モグモグ平方和<-((mean(モグモグ)-全平均)^2)*length(モグモグ)
> パクパク平方和<-((mean(パクパク)-全平均)^2)*length(パクパク)
> 群間平方和<-ワクワク平方和+モグモグ平方和+パクパク平方和
> 群間平方和
[1] 748.6333

平方和がでそろった
ここで、群内の平方和＋群間の平方和が全体の平方和であることがわかる。
> 群内平方和+群間平方和
[1] 2494.183
> 全体平方和
[1] 2494.183

分散分析表をつくる
> 全体自由度<-length(全体)-1 #全体数-1
> 全体自由度
[1] 59
> 群内自由度<-(nrow(全体データ)-1)*(ncol(全体データ)) #(各群データ数-1)をすべての群について合計したもの
> 群内自由度
[1] 57
> 群間自由度<-ncol(全体データ)-1 #群の数-1
> 群間自由度
[1] 2

> 群間平均平方<-群間平方和/群間自由度 #検定統計量Fの分子
> 群間平均平方
[1] 374.3167

> 群内平均平方<-郡内平方和/群内自由度 #検定統計量Fの分母
> 群内平均平方
[1] 30.62368

> F<-群間平均平方/群内平均平方
> F
[1] 12.22311
Rの関数を使い一発で分散分析表をつくることもできる
まずは最も一般的 aov()関数を実行
> aov(全体~お店)
Call:
aov(formula = 全体 ~ お店)

Terms:
お店 Residuals
Sum of Squares 748.6333 1745.5500
Deg. of Freedom 2 57

Residual standard error: 5.533867
Estimated effects may be unbalanced

> summary(aov(全体~お店)) #より詳しく
Df Sum Sq Mean Sq F value Pr(>F)
お店 2 748.6 374.3 12.22 3.82e-05 ***
Residuals 57 1745.6 30.6
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
次に今回のような一元配置分散分析のみ実行可能なoneway.test()関数を実行
> oneway.test(全体~お店,var.equal=TRUE)

One-way analysis of means

data: 全体 and お店
F = 12.223, num df = 2, denom df = 57, p-value = 3.825e-05
最後に、数のモデルの比較など高度な分析に対応 するanova()関数を実行
> anova(lm(全体~お店))
Analysis of Variance Table

Response: 全体
Df Sum Sq Mean Sq F value Pr(>F)
お店 2 748.63 374.32 12.223 3.825e-05 ***
Residuals 57 1745.55 30.62
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

F分布をみる
まずは、5%棄却域を計算する。
> qf(0.05,群間自由度,群内自由度,lower.tail=FALSE)
[1] 3.158843
これを図示すると以下のようになる。
> curve(df(x,群間自由度,群内自由度),　0,13)
> abline(v=qf(0.05,群間自由度,群内自由度,lower.tail=FALSE))

検定統計量Fは12.22なので、5%有意水準で棄却域に入る。したがって帰無仮説は棄却される。
つまり「3つのお店のポテトの評価の平均には差がある」ということである。
これは、少なくともひとつの組み合わせの間で差があるが、 どれとどれの間に差があるかはわからない。

多重比較
どの組み合わせで差があるかということをさらに調べるためには、多重比較という手法を用いなければならない。
今回は多重比較の項ではないので詳しい説明は省くが、ここではTukeyという方法で多重比較を行う。
> #ワクワクとモグモグの多重比較
> q<-abs(mean(ワクワク)-mean(モグモグ))/sqrt(群内平均平方/nrow(全体データ))
> q
[1] 3.636627

> #モグモグとパクパクの多重比較
> q<-abs(mean(モグモグ)-mean(パクパク))/sqrt(群内平均平方/nrow(全体データ))
> q
[1] 6.990406

> #パクパクとワクワクの多重比較
> q<-abs(mean(パクパク)-mean(ワクワク))/sqrt(群内平均平方/nrow(全体データ))
> q
[1] 3.353778

> #多重比較の棄却域
> qtukey(0.95,3,群内自由度)
[1] 3.403189
「ワクワクとモグモグ」および「モグモグとパクパク」は5%の棄却域を上回ったので、帰無仮説は棄却され、これらの2点のお店のポテトの評価の平均には差があるといえる。一方、「ワクワク-パクパク」の検定統計量のみが棄却域を下回ったので帰無仮説が採択される。つまり、 これら2点のお店のポテトの評価の平均には差があるとはいえない。
Rの関数 TukeyHSD ()で一発検定。
> TukeyHSD(aov(全体~お店))
Tukey multiple comparisons of means
95% family-wise confidence level

Fit: aov(formula = 全体 ~ お店)

$`お店`
diff lwr upr p adj
モグモグ-パクパク -8.65 -12.8611413 -4.43885869 0.0000210
ワクワク-パクパク -4.15 -8.3611413 0.06114131 0.0542391
ワクワク-モグモグ 4.50 0.2888587 8.71114131 0.0336261
p値（p adj）からも「モグモグ-パクパク」および「ワクワク-モグモグ」は有意水準の5%より小さいことから、有意差ありと判断できる。
  1. Web教材「ハンバーガーショップで学ぶ楽しい統計学」第6章 分散分析（1要因）http://kogolab.chillout.jp/elearn/hamburger/chap6/sec0.html