R による分散分析(一元配置)

分散分析(ANOVA)は、2 つ以上の実験群の母平均値が同じかどうかを検定する手法の一つである。R では aov 関数および anova 関数などが分散分析を行う関数である。

このページではシミュレーションデータを生成して、分散分析を行う例を説明する。シミュレーションデータは 4 つの実験群からなり、実験群には 10 個体のデータが含まれるものとする。4 つの実験群をそれぞれ A、B、C、D とする。A、B、C の 3 群は平均 100、分散 10 の正規分布からデータをサンプリングする。一方、D 群は平均 140、分散 10 の正規分布からデータをサンプリングする。

A <- rnorm(10, mean = 100, sd = sqrt(10))
B <- rnorm(10, mean = 100, sd = sqrt(10))
C <- rnorm(10, mean = 100, sd = sqrt(10))
D <- rnorm(10, mean = 140, sd = sqrt(10))

d <- data.frame(A = A, B = B, C = C, D = D)
head(d)
##           A         B         C        D
## 1  98.52897  99.32576 101.61557 138.4413
## 2 100.58735  95.54787 103.35295 138.9097
## 3  94.29980  99.15260  97.86004 135.7413
## 4 101.57001 102.17695  97.74807 140.2136
## 5  98.29363 101.43227 101.45534 139.3579
## 6 100.44516  98.56665 104.61613 138.7526

分散分析で利用しやすいように、データ形を変換する。変更後のデータは 2 列で、1 列目が実験群の名前、2 列目にはその実験群の各個体の観測データである。

library(reshape)

x <- melt(d, variable_name = "group")
head(x)
##    group     value
## 1 GroupA  98.52897
## 2 GroupA 100.58735
## 3 GroupA  94.29980
## 4 GroupA 101.57001
## 5 GroupA  98.29363
## 6 GroupA 100.44516

aov 関数

サンプルデータを用意できたところで、aov 関数を用いて分散分析を行う例を示す。

各実験群の母平均が同じかどうかということは、各個体の値はその実験群の影響をうけているかどうかと解釈することができる。つまり、各個体の値 value を実験群 group で説明できるかどうか、ということに言い換えることができる。そこで aov に代入する式として value ~ group とする。以下の例では group は数値ではなく、水準であることを明示するために factor(group) と変換してから利用している。

res <- aov(value ~ factor(group), data = x)

summary(res)
##               Df Sum Sq Mean Sq F value Pr(>F)
## factor(group)  3  12545    4182   440.6 <2e-16 ***
## Residuals     36    342       9
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

factor(group) の行を見ると、Pr(>F) の値が 2×10-16 と非常に小さくなっている。この分散分析の結果として、A、B、C、D の 4 つの群の母平均が同じではないことがわかる。ただし、分散分析ではどの群の母平均が大きかったか、または小さかったかを判定することはできない。

anova 関数

サンプルデータを用意できたところで、anova 関数を用いて分散分析を行う例を示す。anova 関数は一般線形モデルの解析結果から分散分析の結果を計算する関数である。そのため、anova 関数を利用する前に、一般線形モデルによる解析を行う必要がある。

各実験群の母平均が同じかどうかということは、各個体の値はその実験群の影響をうけているかどうかと解釈することができる。つまり、各個体の値 value を実験群 group で説明できるかどうか、ということに言い換えることができる。そこで aov に代入する式として value ~ group とする。以下の例では group は数値ではなく、水準であることを明示するために factor(group) と変換してから利用している。

l <- lm(value ~ factor(group), data = x)
res <- anova(l)

res
## Analysis of Variance Table
## 
## Response: value
##               Df  Sum Sq Mean Sq F value    Pr(>F)
## factor(group)  3 12544.7  4181.6   440.6 < 2.2e-16 ***
## Residuals     36   341.7     9.5
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anovalm の実行結果が aov の実行結果と少し異なっている。これは aov がタイプ I の平行和を、lm はタイプ II の平行和を利用して解析に用いたためである。仮に、以下のように lm 関数にダミー変数を与えると、実行結果が aov 関数と一致するようになる。

l <- lm(value ~ 0 + factor(group), data = x)
res <- anova(l)

res
## Analysis of Variance Table
## 
## Response: value
##               Df Sum Sq Mean Sq F value    Pr(>F)
## factor(group)  4 492146  123036   12964 < 2.2e-16 ***
## Residuals     36    342       9
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1