一元配置分散分析
分散分析：3標本以上の比較方法
以前の項で2標本のt検定を複数回行うと第1種の過誤を犯す確率が高まることを紹介した。それゆえに、もし3標本以上の比較で、第1種の過誤を犯す確率を5％以下に抑えたい場合、今までの検定方法のままでは問題があり、ひと工夫必要となる。そんなときに真っ先に習う方法が分散分析（Analysis of variance、ANOVA）と呼ばれる方法だ。分散分析は1回の検定で3標本以上の平均値の違いを検出することができる。分散分析を使うことができる条件は、すべての標本が等分散の正規分布に従い、独立であることである。t検定と同様である。この時の帰無仮説は、すべての標本の平均値が等しいというものだ。もし、この帰無仮説が棄却されるのであれば、どれかの標本は異なる平均値を持つと言える。ただし、あらかじめ述べておくと、どの標本が異なる平均値をもつかまではわからない。それを知るためには、別の解析が必要になる。が、まずは分散分析について理解を深めよう。
分散分析の検定統計量F
「～検定」というような名前ではないが、分散分析は統計的仮説検定の一種だ。なのでこれまでと同様に検定統計量を算出し、それを帰無分布と比較するという手続きは変わらない。分散分析の検定統計量はFと呼ばれ、このFが従うF分布と比較することで帰無仮説の判定を行う。Fは後述するように分散の比であり、分散分析に限らず分散の比を使う検定は統計量Fを求めることがある。例えば、等分散性を確認する時も統計量Fを求め、この時の検定の名前はF検定と呼ばれる。したがって、分散分析はF検定の仲間、ということもできる。
では、分散分析の手続きを紹介しよう。手順は多いが、一つ一つの手続きは単純なので、恐れずに手順を追ってほしい。今回は例として、X、Y、Zの3標本があるときを紹介する。
まず真っ先にやるべきなのは、各標本の平均、全体の平均、そして自由度と呼ばれる値を計算することだ。これさえ終われば、あとは手続きに従って、Fを計算すればよい。手順1～3は、平均平方和を求める、と書いているが、これは今までの言葉で言うなら、各標本平均や全体平均を基準として、様々な分散を求めていると言える。比較する標本を今回は要因factor、と呼んでいる。しばしばこれらの標本は何かの処理をしたグループであることが多い。例えば、標本Xは無処理だが、Yには薬剤A、Ｚには薬剤Bを処理しているみたいな感じだ。そこで、平均値に影響を与える要素、という意味で要因と呼ばれる。ここで注意だが、標本X、Y、Zと3つ存在するが、要因はこれらをひとまとめにして1つ、である。標本が複数あることを示したいときは、水準levelという言葉を用いる。今回であれば、標本が3つあるので1要因の中に3水準存在することになる。今回のような1要因の分散分析を一元配置分散分析One-way ANOVAと呼ぶ。今回は取り扱わないが、2要因以上の分散分析も存在する。例えば、1つ目の要因は薬剤処理（上記の例なら3水準）、2つ目の要因は生物学的な性別（男女の2水準）、のような感じである。2要因になれば二元配置分散分析Two-way ANOVAだ。標本の数もさらに増えて、3水準 × 2水準 = 6標本存在するような状態となる。さらに、残差residual、と呼ばれる新しい言葉も出てきたが、これは要因では説明できない標本の差分、という意味で、まさに要因で説明される部分の残りということだ。言葉だけではわかりにくいので、次の図を見て理解を深めてほしい。全体平均の平方和が、標本平均に対する平方和と残差平方和で構成されていることもわかるだろう。
手順4では、統計量Fを2つの分散の比率として求めている。このように分散の比率を使って平均値の差を検定するので分散分析と呼ばれている。さてFの意味するところであるが、このFは標本間の平均値の差が大きいほど大きな値を示す。これも直上の図に示した。直上の図の（A）と（B）では残差の平方和の大きさを変えず、標本XとYの平均だけを変更している。2標本の平均の差が大きいほど、残差平方和に対する標本平均を使った平方和の割合が大きくなっていることがわかるだろう。つまり、（A）では要因（標本Xであるか、Yであるか）は平均値に影響を与えているので、平均値の差が大きくなり、Fも大きくなる。逆に、（B）では要因は平均値にほとんど影響を与えておらず、ゆえにFが小さくなる。
言い方を変えれば、標本の数値のばらつきは、要因の影響を受けている、ともいえる。このような状態を、しばしば統計界隈で、要因がばらつきを説明している、と表現することがある。こういう表現をすることもあり、ここでの要因のことを説明変数（独立変数）なんて呼ぶ。一方で、要因に説明される側のこと（ここでは数値そのもの）を被説明変数（従属変数）と呼ぶ。ただし、説明と言っているが、あくまで形式上であり、真の因果関係を保証してはいないので注意が必要である。初めて出くわしたときは全く意味が分からないと思うので、ここで述べておく。
この統計量Fは、すべての標本が正規分布に従い、かつ平均が等しいときに従う分布、F分布の形が知られている。あとは、t検定の時と同じ手続きだ。F分布の端から5%の領域は、すべての標本平均が等しいときにはめったに生じないFの値だ。もし、分散分析で5%の領域に含まれるようなFが算出されたら、それは帰無仮説が正しいと考えることに無理があった、つまり標本の平均値に違いがある、と考えるわけだ。
分散分析の危険率
さて、ここからは、毎回恒例のRを使ったシミュレートだ。分散分析におけるF分布をシミュレーションで作ってみるとともに、分散分析の危険率を検討しよう。そもそもの話、複数回検定を行うと危険率が上がるので、別の方法を導入しようということで紹介したのが分散分析だ。分散分析の危険率が高かったら話にならない。では、平均値の等しい３つの標本に対して分散分析を行う状況を考える。すべての標本が平均50、標準偏差20、サンプルサイズ30にそろえてある。
------------------------------------------------------
nA <- 30
nB <- 30
nC <- 30
meA <- 50
meB <- 50
meC <- 50
s <- 20
p <- NULL
Fvalue <- NULL
for (i in 1:10000){
mA <- rnorm(nA, mean = meA, sd = s)
mB <- rnorm(nB, mean = meB, sd = s)
mC <- rnorm(nC, mean = meC, sd = s)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC))) #d列に標本の数値、g列に標本名を入れている。
fit1 <- aov(d$d ~ d$g) #分散分析を行うコマンド、被説明変数~説明変数と記述する。
fit2 <- summary(fit1)
mT <- mean(c(mA,mB,mC))
m_mA <- mean(mA)
m_mB <- mean(mB)
m_mC <- mean(mC)
dfT <- nA + nB + nC - 1 #全体の自由度
dfF <- 3 - 1 #要因の自由度
dfR <- dfT - dfF #残差の自由度
sT <- sum((c(mA,mB,mC) - mT)^2)/dfT #全体平均に対する平均平方和
sF <- sum(((m_mA - mT)^2)*nA, ((m_mB - mT)^2)*nB, ((m_mC - mT)^2)*nC)/dfF #標本平均に対する平均平方和
sR <- (sT*dfT - sF*dfF)/dfR #残差平方和
Fvalue <- c(Fvalue, sF/sR) #検定統計量F
p <- c(p, fit2[[1]]$`Pr(>F)`[1]) #分散分析のp値を取り出す
}
histn(Fvalue, xlab = "統計量F", ylab = "頻度") #図1の描画
overdraw(abline(v = qf(df1 = dfF, df2 = dfR, p = 0.95), col = "red"))
histn(p, xlab = "P value", ylab = "頻度") #図2の描画
sum(Fvalue > qf(df1 = dfF, df2 = dfR, p = 0.95)) #F分布は比をとる2つの分散の自由度を変数にとる（df1とdf2）
## [1] 486
sum(p < 0.05)
## [1] 486
------------------------------------------------------
図1 F分布。赤線（qf(df1 = dfF, df2 = dfR, p = 0.95)）は帰無仮説が正しいときの95%区間
図2 分散分析のp値の分布
図1はt検定を解説した時と同様に、帰無仮説が正しいときの検定統計量が従う分布である。赤線よりも右側の、「極端に大きな値」が出た場合は帰無仮説を正しいと考えることに無理があると判断する。図2は10000回、分散分析をした時のp値の分布である。F分布をもとに判断した場合（sum(Fvalue > qf(df1 = dfF, df2 = dfR, p = 0.95))）も、分散分析のp値で判断した場合（sum(p < 0.05)）も危険率はともに約5%であり、当初の危険率の問題をクリアできている。繰り返すが、分散分析は１回の検定で3標本以上の平均値の差があるかを判断できるので、複数回検定の問題が生じない。
分散分析の検出力
続いて分散分析の検出力を検討しよう。危険率ばかりが低く、検出力も低いのでは、これはこれで話にならない。さて、帰無仮説はすべての標本の平均値が等しい、というものだった。3標本（平均値をA、B、Cとする）の場合は、A = B = Cを考える。これに対する対立仮説は、どれかの標本の平均値が異なっている、となる。3標本の場合、平均値を大きいものからA、B、Cとしたとき、次の3パターンである：（1）A = B > C、（2）A > B = C、（3）A > B > C。ここで平均値の大きなものから名前を付けたのは、要因の順番には意味がないと考えているからである。例えば、標本Xには薬剤を0、Yには10、Zには100みたいな、対立仮説に要因の順序が反映される可能性がある場合は、別の手段が必要になる。今回はシミュレーションなので、要因の順序はどうでもよく、平均値の大小関係だけで分類した。対立仮説さえ分かれば、あとは帰無仮説の検証の時と同様にシミュレートできる。まずは、A = B > Cの場合だ。
------------------------------------------------------
nA <- 30
nB <- 30
nC <- 30
meA <- 70
meB <- 70
meC <- 50
s <- 20
p <- NULL
Fvalue <- NULL
for (i in 1:10000){
mA <- rnorm(nA, mean = meA, sd = s)
mB <- rnorm(nB, mean = meB, sd = s)
mC <- rnorm(nC, mean = meC, sd = s)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
fit1 <- aov(d$d ~ d$g)
fit2 <- summary(fit1)
mT <- mean(c(mA,mB,mC))
m_mA <- mean(mA)
m_mB <- mean(mB)
m_mC <- mean(mC)
dfT <- nA + nB + nC - 1
dfF <- 3 - 1
dfR <- dfT - dfF
sT <- sum((c(mA,mB,mC) - mT)^2)/dfT
sF <- sum(((m_mA - mT)^2)*nA, ((m_mB - mT)^2)*nB, ((m_mC - mT)^2)*nC)/dfF
sR <- (sT*dfT - sF*dfF)/dfR
Fvalue <- c(Fvalue, sF/sR)
p <- c(p, fit2[[1]]$`Pr(>F)`[1])
}
histn(Fvalue, xlab = "統計量F", ylab = "頻度") #図3の描画
overdraw(abline(v = qf(df1 = dfF, df2 = dfR, p = 0.95), col = "red"))
histn(p, xlab = "P value", ylab = "頻度") #図4の描画
sum(Fvalue > qf(df1 = dfF, df2 = dfR, p = 0.95))
## [1] 9822
sum(p < 0.05)
## [1] 9822
------------------------------------------------------
図3 F分布。A = B > Cの時。赤線（qf(df1 = dfF, df2 = dfR, p = 0.95)）は帰無仮説が正しいときの95%区間
図4 分散分析のp値の分布。A = B > Cの時。
次に、A > B = Cの場合だ。
------------------------------------------------------
nA <- 30
nB <- 30
nC <- 30
meA <- 70
meB <- 50
meC <- 50
s <- 20
p <- NULL
Fvalue <- NULL
for (i in 1:10000){
mA <- rnorm(nA, mean = meA, sd = s)
mB <- rnorm(nB, mean = meB, sd = s)
mC <- rnorm(nC, mean = meC, sd = s)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
fit1 <- aov(d$d ~ d$g)
fit2 <- summary(fit1)
mT <- mean(c(mA,mB,mC))
m_mA <- mean(mA)
m_mB <- mean(mB)
m_mC <- mean(mC)
dfT <- nA + nB + nC - 1
dfF <- 3 - 1
dfR <- dfT - dfF
sT <- sum((c(mA,mB,mC) - mT)^2)/dfT
sF <- sum(((m_mA - mT)^2)*nA, ((m_mB - mT)^2)*nB, ((m_mC - mT)^2)*nC)/dfF
sR <- (sT*dfT - sF*dfF)/dfR
Fvalue <- c(Fvalue, sF/sR)
p <- c(p, fit2[[1]]$`Pr(>F)`[1])
}
histn(Fvalue, xlab = "統計量F", ylab = "頻度") #図5の描画
overdraw(abline(v = qf(df1 = dfF, df2 = dfR, p = 0.95), col = "red"))
histn(p, xlab = "P value", ylab = "頻度") #図6の描画
sum(Fvalue > qf(df1 = dfF, df2 = dfR, p = 0.95))
## [1] 9800
sum(p < 0.05)
## [1] 9800
------------------------------------------------------
図5 F分布。A > B = Cの時。赤線（qf(df1 = dfF, df2 = dfR, p = 0.95)）は帰無仮説が正しいときの95%区間
図6 分散分析のp値の分布。A > B = Cの時。
最後に、A > B > Cの場合だ。
------------------------------------------------------
nA <- 30
nB <- 30
nC <- 30
meA <- 70
meB <- 60
meC <- 50
s <- 20
p <- NULL
Fvalue <- NULL
for (i in 1:10000){
mA <- rnorm(nA, mean = meA, sd = s)
mB <- rnorm(nB, mean = meB, sd = s)
mC <- rnorm(nC, mean = meC, sd = s)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
fit1 <- aov(d$d ~ d$g)
fit2 <- summary(fit1)
mT <- mean(c(mA,mB,mC))
m_mA <- mean(mA)
m_mB <- mean(mB)
m_mC <- mean(mC)
dfT <- nA + nB + nC - 1
dfF <- 3 - 1
dfR <- dfT - dfF
sT <- sum((c(mA,mB,mC) - mT)^2)/dfT
sF <- sum(((m_mA - mT)^2)*nA, ((m_mB - mT)^2)*nB, ((m_mC - mT)^2)*nC)/dfF
sR <- (sT*dfT - sF*dfF)/dfR
Fvalue <- c(Fvalue, sF/sR)
p <- c(p, fit2[[1]]$`Pr(>F)`[1])
}
histn(Fvalue, xlab = "統計量F", ylab = "頻度") #図7の描画
overdraw(abline(v = qf(df1 = dfF, df2 = dfR, p = 0.95), col = "red"))
histn(p, xlab = "P value", ylab = "頻度") #図8の描画
sum(Fvalue > qf(df1 = dfF, df2 = dfR, p = 0.95))
## [1] 9376
sum(p < 0.05)
## [1] 9376
------------------------------------------------------
図7 F分布。A > B > Cの時。赤線（qf(df1 = dfF, df2 = dfR, p = 0.95)）は帰無仮説が正しいときの95%区間
図8 分散分析のp値の分布。A > B > Cの時。
今回の解析では、標本が標準偏差20の正規分布から、サンプルサイズ30でサンプリングされるとき、標本の平均値の差が最大20程度あると、検出力が90%以上となり、申し分ない状態となる。もちろん、平均値の差が小さくなるほど検出力は落ちるので、あらゆる分析で90%以上の検出力が担保されるという意味ではないので注意しよう。このように分散分析を使えば、3標本以上の平均値の比較において、すべての標本の平均値が同じというわけではない、という判断を下すことができる。
Rで行う分散分析
上記のような背景のもと、分散分析を行ってみよう。以下のようなサンプルサイズ30の標本が3つあるとする。ともに同じ標準偏差の正規分布からデータを生成していると考えてもらってよい。有意水準は5%としよう。まずは、解析の前にデータの様子を確認だ。
------------------------------------------------------
mA <- c(79, 66, 47, 17, 78, 41, 83, 91, 9, 43, 64, 41, 58, 53, 62,
55, 30, 46, 64, 34, 26, 47, 55, 23, 51, 51, 26, 67, 40, 70)
mB <- c(37, 64, 70, 41, 43, 38, 46, 62, 21, 64, 31, 73, 99, 46, 29,
53, 19, 68, 70, 51, 103, 30, 61, 7, 46, 69, 40, 39, 53, 49)
mC <- c(60, 50, 80, 57, 59, 84, 78, 50, 51, 68, 68, 58, 69, 104, 61,
74, 69, 99, 46, 91, 60, 75, 82, 65, 47, 81, 91, 70, 71, 63)
nA <- length(mA)
nB <- length(mB)
nC <- length(mC)
d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))
boxplotn(d$d ~ d$g, ylab = "data", jitter.method = "swarm") #図9の描画
------------------------------------------------------
図9 標本A、B、Cの描画
これを見ると標本Cでは平均値が大きそうな様子がわかる。では、実際に違いそうかを解析してみよう。
------------------------------------------------------
fit1 <- aov(d$d ~ d$g) #分散分析
fit2 <- summary(fit1)
fit2
## Df Sum Sq Mean Sq F value Pr(>F)
## d$g 2 7007 3503 9.617 0.000168 ***
## Residuals 87 31692 364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------
まず、aov()関数で分散分析を行っている。その解析結果をsummary()関数を通すことで、解析結果の要約を確認することができる。fit2には要約を格納しているので、fit2と打って中身を確認しよう。d$gは要因のことで、Residualsは残差のことである。Dfは自由度のことで、要因と残差の自由度は以下の方法で計算できる。
------------------------------------------------------
dfT <- nA + nB + nC - 1
dfF <- 3 - 1
dfR <- dfT - dfF
dfT #全体の自由度
## [1] 89
dfF #要因の自由度
## [1] 2
dfR #残差の自由度
## [1] 87
------------------------------------------------------
次のSum SqとMean Sqは平方和（Sum of Square）と平均平方和（Mean of sum of square）のことで以下のように計算する。
------------------------------------------------------
mT <- mean(c(mA,mB,mC))
m_mA <- mean(mA)
m_mB <- mean(mB)
m_mC <- mean(mC)
ST <- sum((c(mA,mB,mC) - mT)^2) #全体の平方和
SF <- sum(((m_mA - mT)^2)*nA, ((m_mB - mT)^2)*nB, ((m_mC - mT)^2)*nC) #要因の平方和
sT <- sum((c(mA,mB,mC) - mT)^2)/dfT #全体の平均平方和
sF <- sum(((m_mA - mT)^2)*nA, ((m_mB - mT)^2)*nB, ((m_mC - mT)^2)*nC)/dfF #要因の平均平方和
SR <- (sT*dfT - sF*dfF) #残差の平方和
sR <- (sT*dfT - sF*dfF)/dfR #残差の平均平方和
ST
## [1] 38698.89
SF
## [1] 7006.689
SR
## [1] 31692.2
sT
## [1] 434.819
sF
## [1] 3503.344
sR
## [1] 364.2782
------------------------------------------------------
次にF valueは以下のように計算する。
------------------------------------------------------
Fvalue <- sF/sR
Fvalue
## [1] 9.617223
------------------------------------------------------
最後にPr(>F)はいわゆるp値のことで、F分布を利用して以下のように計算できる。
------------------------------------------------------
pf(df1 = dfF, df2 = dfR, q = Fvalue, lower.tail = F)
## [1] 0.0001684784
------------------------------------------------------
p値が有意水準を下回ったので、3標本のいずれかは平均値が違いそうだ、とわかる。何かしら要因の効果があったと、言及するだけであれば分散分析で事足りる。しかし、普通はこれらの標本の中で、どの標本ペアの比較に差があるといっていいのか、を知りたいものである。分散分析ではそこまで言及することはできない。次回はさらにその先の解析を紹介しよう。