独立多群の差の検定
分散分析は独立多群に差があるかどうか検定します。統計量はF値であり、このためF分布に対する理解が基本となります。
http://hyoccori-toukei.la.coocan.jp/id-2/id-2/id7.html#2_分散分析
1 １元配置分散分析
このページは、市原清志先生の「バイオサイエンスの統計学」の問題を、統計ソフト”R”で解説するものです。 詳細については、市原先生のテキスト（ p152 ）を参照ください。
1.1 例5-1（p152）
データを読み込みます。
a1<-c(13,11,6)
a2<-c(11,10,7,7,5)
a3<-c(8,7,5,5,4,3,3)縦型のデータフレームに直します。
A1<-data.frame(A=rep("x1",3),B=a1)
A2<-data.frame(A=rep("x2",5),B=a2)
A3<-data.frame(A=rep("x3",7),B=a3)
data<-rbind(A1,A2,A3)
colnames(data) <- c("groupe","score")
data## groupe score
## 1 x1 13
## 2 x1 11
## 3 x1 6
## 4 x2 11
## 5 x2 10
## 6 x2 7
## 7 x2 7
## 8 x2 5
## 9 x3 8
## 10 x3 7
## 11 x3 5
## 12 x3 5
## 13 x3 4
## 14 x3 3
## 15 x3 3分散分析を行い、箱ヒゲ図でも表示します。
modelANOVA <- lm(data$score ~ data$groupe, data=data)
anova(modelANOVA)## Analysis of Variance Table
##
## Response: data$score
## Df Sum Sq Mean Sq F value Pr(>F)
## data$groupe 2 60 30 5 0.02634 *
## Residuals 12 72 6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1boxplot(score ~ groupe, data=data)
群間の変動が群内変動よりも大きく、群間に差があると判断できます。
1.2 演習5-1（p166）
a<-data.frame(A=rep("a",4),B=c(10,9,9,7));a## A B
## 1 a 10
## 2 a 9
## 3 a 9
## 4 a 7b<-data.frame(A=rep("b",5),B=c(8,7,6,6,5));b## A B
## 1 b 8
## 2 b 7
## 3 b 6
## 4 b 6
## 5 b 5c<-data.frame(A=rep("c",4),B=c(11,10,9,8));c## A B
## 1 c 11
## 2 c 10
## 3 c 9
## 4 c 8d<-data.frame(A=rep("d",3),B=c(6,4,4));d## A B
## 1 d 6
## 2 d 4
## 3 d 4e<-data.frame(A=rep("e",4),B=c(8,7,7,6));e## A B
## 1 e 8
## 2 e 7
## 3 e 7
## 4 e 6data1<-rbind(a,b,c,d,e);data1## A B
## 1 a 10
## 2 a 9
## 3 a 9
## 4 a 7
## 5 b 8
## 6 b 7
## 7 b 6
## 8 b 6
## 9 b 5
## 10 c 11
## 11 c 10
## 12 c 9
## 13 c 8
## 14 d 6
## 15 d 4
## 16 d 4
## 17 e 8
## 18 e 7
## 19 e 7
## 20 e 6boxplot(B ~ A, data=data1)
modelANOVA <- lm(B ~ A, data=data1)
anova(modelANOVA)## Analysis of Variance Table
##
## Response: B
## Df Sum Sq Mean Sq F value Pr(>F)
## A 4 52.933 13.2333 10.119 0.0003556 ***
## Residuals 15 19.617 1.3078
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12 シミュレーション
分散分析を理解するには、帰無仮説のもとでデータがどのようになるかを見ておく必要がある。
2.1 ３群の差がないとき
郡内平均がどの群も７、標準偏差がどの群も3.5とする。
この条件で各群のデータを乱数で生成する。
A1<-data.frame(A=rep("x1",3),B=rnorm(3,7,3.5))
A2<-data.frame(A=rep("x2",5),B=rnorm(5,7,3.5))
A3<-data.frame(A=rep("x3",7),B=rnorm(7,7,3.5))
data<-rbind(A1,A2,A3)
colnames(data) <- c("groupe","score")
data## groupe score
## 1 x1 4.6385325
## 2 x1 9.6048238
## 3 x1 10.1387121
## 4 x2 8.4156380
## 5 x2 8.7584708
## 6 x2 5.6372892
## 7 x2 3.0050673
## 8 x2 10.4601363
## 9 x3 4.4156972
## 10 x3 -0.3624942
## 11 x3 13.5033322
## 12 x3 3.1532504
## 13 x3 12.3818076
## 14 x3 4.9798915
## 15 x3 10.4311675modelANOVA <- lm(data$score ~ data$groupe, data=data)
anova(modelANOVA)## Analysis of Variance Table
##
## Response: data$score
## Df Sum Sq Mean Sq F value Pr(>F)
## data$groupe 2 3.02 1.5098 0.084 0.92
## Residuals 12 215.72 17.9765boxplot(score ~ groupe, data=data)
F値は有意にならなりません。
そもそもF値とはなんでしょうか？
\(F=\dfrac{S_A/(a-1)}{S_e/(n-a)}\sim
F(a-1,n-a)\)
\(S_A\) は群間変動、\(S_e\) は群内変動。\((a-1),~~(n-a)\)
はそれぞれの自由度です。
自由度は実質的なデータの個数です。
例えば、10個のデータがあり自由度10としましょう。
データ数が少ないからといって、同じデータを10回使って100個のデータとして分析しても自由度は100にはなりません。自由度は10のままです。
そのデータには10個の価値しかないからです。偏差平方和を自由度で割れば分散が得られます。東京大学
大森宏先生の説明では、
分子の \(S_A/(a−1)\) は、処理平均平方（treatment mean square）といい、処理平均から求めた 誤差分散 \(σ^2\) の推定値で群間分散と呼ばれます。
\(S_e/(n−a)\) は、誤差平均平方（error mean square）という誤差分散 \(σ^2\) の推定値で群内分散と呼ばれています。
F値は母分散の推定値の比なので値は１に近いのです。
帰無仮説 「３群に差がない」 のもとでは、群間分散も群内分散も母分散にほぼ一致するのでF値は１に近い値をとります。
上の操作を1000回行って確かめてみます。
SE = c()
SR = c()
X=c()
for (i in 1:1000) {
x1<-rnorm(3,7,3.5);x2<-rnorm(5,7,3.5);x3<-rnorm(7,7,3.5)
se <- var(x1)*2+var(x2)*4+var(x3)*6 #群内変動
SE <- c(SE,se)
sr <- 3*(mean(x1)-mean(c(x1,x2,x3)))^2+5*(mean(x2)-mean(c(x1,x2,x3)))^2+7*(mean(x3)-mean(c(x1,x2,x3)))^2 #群間変動
SR <- c(SR,sr)
f <- (sr/2)/(se/12)#F値
X <- c(X,f)
}
mean(X)## [1] 1.160657hist(X,prob=T,breaks = "Scott")
curve(df(x, 2, 12),add=T)
F値は１に近くなっている。
群間分散と群内分散の1000回の平均をみると
mean(SE/12)## [1] 12.24235#hist(SE/12)
mean(SR/2)## [1] 11.55941#hist(SR/2)
cor(SE,SR)## [1] -0.03649154群間分散、群内分散とも \(3.5^2\)
に近い値となっている。
相関係数がほぼ０であるのは、F値の式の分母と分子が独立であることを意味する。
2.2 ３群の差があるとき
郡内平均を 10, 8, 5 とし、標準偏差はどの群も3.5とする。
この条件で各群のデータを乱数で生成する。
SE = c()
SR = c()
X=c()
for (i in 1:1000) {
x1<-rnorm(3,10,3.5);x2<-rnorm(5,8,3.5);x3<-rnorm(7,5,3.5)
se <- var(x1)*2+var(x2)*4+var(x3)*6 #群内変動
SE <- c(SE,se)
sr <- 3*(mean(x1)-mean(c(x1,x2,x3)))^2+5*(mean(x2)-mean(c(x1,x2,x3)))^2+7*(mean(x3)-mean(c(x1,x2,x3)))^2 #群間変動
SR <- c(SR,sr)
f <- (sr/2)/(se/12)#F値
X <- c(X,f)
}
mean(X)## [1] 4.137319F値はF>4 となる。上のヒストグラムで見ると、F>4 は帰無仮説のもとでは確率的に起こらないと判断され、帰無仮説は棄却される。
このように、Rを用いて帰無仮説のもとでのF値の動きをシミュレーションすることっで、F値によって異常が検知できることが見えてくる。