F分布
これまで学んだことを整理します。
- 独立な2つの確率変数から合成された合成確率変数は、どんな分布に従うか？
（和、差、積、商、偏差平方和などの分布）
（例：正規分布とカイ二乗分布の平方根の商はt分布、偏差平方和はカイ二乗分布）
- 帰無仮説のもとで、合成確率変数の分布はどのような分布に従うか？
（統計量が有意差の範囲に入れば、帰無仮説は棄却されます）
- 帰無仮説の分布のもとで、対立仮説の分布の期待値はどこに位置するか？
（対立仮説の分布のもとで帰無仮説の棄却境界値はどこに位置するかが検出力を決定）
1 等分散の検定
独立な2つの確率変数\(x_A,~x_B\)
が出力した\(n_A\)個のサンプルAと \(n_B\)個のサンプルBをもとに、\(x_A,~x_B\)
の分散が等しいかどうか検定します。
どのような合成確率変数を用い、どのような分布になるかに注目します。
合成確率変数としては、独立な2つの確率変数の不偏分散の比を考えます。
等分散という帰無仮説のもとでは次のことが成立します。
\(N(\mu_1,\sigma^2)\) の\(n_A\)個の偏差平方和を\(S_A\)とすると、 \(\dfrac{S_A}{\sigma^2}\sim
\chi^2(n_A-1)\)
\(N(\mu_2,\sigma^2)\) の\(n_B\)個の偏差平方和を\(S_B\)とすると、 \(\dfrac{S_B}{\sigma^2}\sim
\chi^2(n_B-1)\)
共通の母分散は商の計算で約分されるので、
\(F=\dfrac{S_A/(n_A-1)}{S_B/(n_B-1)}\sim
F(n_A-1,n_B-1)\) 期待値：\(E(F) =
\dfrac{(n_B-1)}{(n_B-1)-2}\)
期待値が\(n_B\)だけによって決まることは驚くべきことですが、乱数によるシミュレーションで確かめてみましょう。
X = c()
for (i in 1:10000) {
A <- rnorm(4,1,1)
B <- rnorm(10,1,1)
p <-var(A)/var(B)
X <-c(X,p)
}
mean(X)## [1] 1.299261E <- (10-1)/(10-3);E## [1] 1.285714hist(X,prob=T)
curve(df(x, 3, 9),add=T)
公式による期待値とシミュレーションによる平均は、ほぼ一致しています。
分子の自由度を変えてシミュレーションしてみましょう。
X = c()
for (i in 1:10000) {
A <- rnorm(6,1,1)
B <- rnorm(10,1,1)
p <-var(A)/var(B)
X <-c(X,p)
}
mean(X)## [1] 1.293673E <- (10-1)/(10-3);E## [1] 1.285714hist(X,prob=T)
curve(df(x, 5, 9),add=T)
確かに期待値は、分母の自由度のみによって決まるようです。
シミュレーションによる分散比の分布と、var.test の結果が一致するか見てみましょう。
X = c()
for (i in 1:1000) {
A <- rnorm(20,1,1)
B <- rnorm(30,5,1)
p <-var(A)/var(B)
X <-c(X,p)
}
mean(X)## [1] 1.04734E <- (30-1)/(30-3);E## [1] 1.074074hist(X,prob=T)
curve(df(x, 19, 29),add=T)
r の関数 var.testで等分散検定をすると、
A <- rnorm(20,1,1)
B <- rnorm(30,5,1)
p <-var(A)/var(B);p## [1] 1.203964res <- var.test(A,B)
res$estimate## ratio of variances
## 1.203964res$p.value## [1] 0.6371428等分散仮説は棄却されず、妥当な検定結果と言えます。
2 分散分析
独立な2つの確率変数\(x_A,~x_B\)
が出力した\(n_A\)個のサンプルAと \(n_B\)個のサンプルBをもとに、\(x_A,~x_B\)
の平均が等しいかどうか検定します。
この判断は、\(x_A,~x_B\)
の分散が等しいという前提のもとで、初めて可能になります。
通常は、帰無仮説「分散が等しい」が棄却されないことを確認してからこの段階に進みます。
では、どのような合成確率変数を用いればよいでしょうか。
東京大学
大森宏 様 によれば、
帰無仮説「群間平均に差がない」のもとで、
\(E(S_e)=(n-a)\sigma^2\)
より、\(S_e\) の自由度は \((n-a)\) \(\therefore \dfrac{S_e}{\sigma^2}\sim
\chi^2(n-a)\)
\(E(S_A)=(a-1)\sigma^2\)より、\(S_A\) の自由度は \((a-1)\) \(\therefore \dfrac{S_A}{\sigma^2}\sim
\chi^2(a-1)\)
大森先生によれば、
・\(S_A/(a-1)\)
は、処理平均平方（treatment mean square）と呼ばれ、処理平均から求めた
誤差分散 σ2 の推定値
\(S_e/(n-a)\)
は、誤差平均平方と呼ばれ、（error mean square）誤差分散 σ2 の推定値
合成確率変数としては、
\(F=\dfrac{S_A/(a-1)}{S_e/(n-a)}\sim
F(a-1,n-a)\) これが分散分析のF値です。
このことをシミュレーションで確かめてみましょう。
SE = c()
SR = c()
X=c()
for (i in 1:1000) {
x1<-rnorm(30,5,1);x2<-rnorm(30,5,1)
se <- var(x1)*29+var(x2)*29 #群内変動
SE <- c(SE,se)
sr <- 30*(mean(x1)-mean(c(x1,x2)))^2+30*(mean(x2)-mean(c(x1,x2)))^2 #群間変動
SR <- c(SR,sr)
f <- (sr)/(se/58)#F値
X <- c(X,f)
}
mean(X)## [1] 1.028515hist(X,prob=T,breaks = "Scott")
curve(df(x, 1, 58),add=T)
mean(SE/58)## [1] 1.006167mean(SR/1)## [1] 0.9937356cor(SE,SR)## [1] -0.02239687\(E(S_e)=(n-a)\sigma^2=(60-2)\cdot
1\) \(E(S_A)=(a-1)\sigma^2=(2-1)\cdot 1\)
\(cor(SE,SR)\fallingdotseq 0\) より
SE,SR は無相関なので独立と予想され、
\(F=\dfrac{S_A/(a-1)}{S_e/(n-a)}\sim F(a-1,n-a)\)
２群が等分散であり、２群の平均も等しいとき、F値は自由度(a-1,n-a) のF分布に従います。ヒストグラムと理論上のグラフは一致しています。
r の関数 aov を適用すると、
A <- rnorm(30,5,1)
B <- rnorm(30,5,1)
x<-data.frame(group=factor(c(rep("A",30),rep("B",30))),value=c(A,B))
boxplot(value ~ factor(group), data = x)
res <- aov(value ~ factor(group), data = x)
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## factor(group) 1 0.34 0.3384 0.367 0.547
## Residuals 58 53.49 0.9222サンプル数を変えてみます。
群間平均に差がないとき、
SE = c()
SR = c()
F1 = c()
for (i in 1:1000) {
x1<-rnorm(20,5,1);x2<-rnorm(10,5,1)
se <- var(x1)*19+var(x2)*9 #群内変動
SE <- c(SE,se)
sr <- 20*(mean(x1)-mean(c(x1,x2)))^2+10*(mean(x2)-mean(c(x1,x2)))^2 #群間変動
SR <- c(SR,sr)
f <- (sr)/(se/28)#F値
F1 <- c(F1,f)
}
hist(F1,prob=T,breaks = "Scott")
curve(df(x, 1, 28),add=T)
mean(SE/28)## [1] 1.004944mean(SR/1)## [1] 0.9921942Fの出力値と理論上のグラフは、ほぼ一致しています。
間平均に差があるとき、
\(S_A/(a-1)\) は誤差分散 σ2
の推定値より大きくなります。
大森先生によれば、\(E(S_A)=\sum_i n_i
\alpha_i^2+(a-1)\sigma^2\) より、
\(E\left(\dfrac{S_A}{a-1}\right)=\dfrac{\sum_i n_i
\alpha_i^2}{a-1}+\sigma^2\)
SE = c()
SR = c()
F2 = c()
for (i in 1:1000) {
x1<-rnorm(20,6,1);x2<-rnorm(20,5,1)
se <- var(x1)*19+var(x2)*19 #群内変動
SE <- c(SE,se)
sr <- 20*(mean(x1)-mean(c(x1,x2)))^2+20*(mean(x2)-mean(c(x1,x2)))^2 #群間変動
SR <- c(SR,sr)
f <- (sr)/(se/38)#F値
F2 <- c(F2,f)
}
hist(F2,prob=T,breaks = "Scott")
curve(df(x, 1, 38),add=T)
mean(SE/38)## [1] 0.9986498mean(SR/1)## [1] 11.19774Fの出力値は、理論上の値より大きい値をとります。
参考までに、大森先生によれば
SR/1 の期待値は、\(20(0.5^2+0.5^2))+1=11\) となります。
hist(F2,prob=T,col="pink",breaks = "Scott")
hist(F1,prob=T,col="lightblue",breaks = "Scott",add=T)
水色は帰無仮説の場合のF値の分布でF(19,19)に従います。
ピンクは２群の平均差を１としたとき分布で、右に裾が大きいため値が大きくなってF分布にならず、帰無仮説は棄却され易いことがわかります。
つまり検定の検出力が高いことが予想されます。
組み込み関数 aov を用いた検定では、
A <- rnorm(20,6,1)
B <- rnorm(20,5,1)
x<-data.frame(group=factor(c(rep("A",20),rep("B",20))),value=c(A,B))
boxplot(value ~ factor(group), data = x)
res <- aov(value ~ factor(group), data = x)
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## factor(group) 1 4.04 4.042 4.656 0.0373 *
## Residuals 38 32.99 0.868
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#t.test(C,D,paired=FALSE,alternative ="two.sided",var.equal=T)帰無仮説「群間平均に差がない」は棄却されます。