分散分析（analysis of variance; ANOVA）は誤差の変動に着目して、その変動を誤差変動と各要因の変動に分解して、各要因の効果を調べる仮説検定の一つである。分散分析は、データの構造によって一次元配置分散分析（one-way ANOVA）や多次元配置分散分析などに分類することができるが、どちらも変動を、誤差変動と各要因の変動に分解することによって仮説検定を行う方法である。このページでは、一元配置分散分析を説明する。一元配置分散分析（以下、分散分析）は、複数の実験群の母平均が同じであるかどうかを検定する仮説検定の一つである。
母平均の 2 群間比較を例に分散分析を説明すると次のようになる。A 群と B 群を考えて、それぞれの平均を \( \overline{X_{A}} \) と \( \overline{X_{B}} \) とする。また、全データの平均を総平均 \( \frac{\overline{X_{A}} + \overline{X_{B}}}{2} \) とする。このとき、A 群に含まれるあるデータ \( X_{A1} \) が総平均からズレ（総変動）を考えたとき、その総変動を \( X_{A1} \) が A 群平均からズレ（群内変動）と A 群平均が総平均からのズレ（群間変動）に分解することができる。
群内変動と群間変動の相対的な大きさに着目すると、群内変動が群間変動よりも著しく大きければ、総変動はほとんど群内変動に影響される。つまり、全データをみたとき、そのばらつきは、ほとんど群内変動によるものとみなすことができる。したがって、データ全体が 1 つの実験群からなるとみることができるようになる。逆に、群間変動が郡内変動よりも著しく大きければ、データ全体のばらつきは、ほとんど群間変動によるものとみなすことができる。つまり、データ全体のばらつきは、ほとんど A 群平均と総平均の差に由来する。よって、A 群平均とそれ以外の実験群の平均に有意差が見られるといえる。
3 群以上の母平均の比較の場合も同様に考えることができる。すなわち、データの総変動を群内変動と群間変動に分解し、両者の相対的な大小を比べる。
分散分析と F 検定
比較対象の実験群が r 群あり、実験群 i で観測されたデータ数を ni とする。実験群 i の j 番目の観測データを xij と表す。このとき、全データに対して、総平均からの変動を計算し、その和を求める。すなわち、総平均に対する誤差平方和を計算する。
\[ S_{T} = \sum_{i=1}^{r}\sum_{j=1}^{n_{i}}\left( x_{ij} - \overline{x_{\cdot \cdot}} \right)^{2} \]この総変動を群内変動と群間変動に分解していく。
\[ \begin{eqnarray} S_{T} &=& \sum_{i=1}^{r}\sum_{j=1}^{n_{i}}\left( x_{ij} - \overline{x_{i\cdot}} + \overline{x_{i\cdot}} - \overline{x_{\cdot \cdot}} \right)^{2} \\ &=& \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} \left( x_{ij} - \overline{x_{i\cdot}} \right)^{2} + \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} 2 \left( x_{ij} - \overline{x_{i\cdot}} \right) \left( \overline{x_{i\cdot}} - \overline{x_{\cdot \cdot}} \right) + \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} \left( \overline{x_{i\cdot}} - \overline{x_{\cdot \cdot}} \right)^{2} \\ &=& \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} \left( x_{ij} - \overline{x_{i\cdot}} \right)^{2} + 2 \sum_{i=1}^{r} \left( \overline{x_{i\cdot}} - \overline{x_{\cdot \cdot}} \right) \sum_{j=1}^{n_{i}} \left( x_{ij} - \overline{x_{i\cdot}} \right) + \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} \left( \overline{x_{i\cdot}} - \overline{x_{\cdot \cdot}} \right)^{2} \\ &=& \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} \left( x_{ij} - \overline{x_{i\cdot}} \right)^{2} + 2 \sum_{i=1}^{r} \left( \overline{x_{i\cdot}} - \overline{x_{\cdot \cdot}} \right) \cdot 0 + \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} \left( \overline{x_{i\cdot}} - \overline{x_{\cdot \cdot}} \right)^{2} \\ &=& \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} \left( x_{ij} - \overline{x_{i\cdot}} \right)^{2} + \sum_{i=1}^{r}\sum_{j=1}^{n_{i}} \left( \overline{x_{i\cdot}} - \overline{x_{\cdot \cdot}} \right)^{2} \\ &=& S_{w} + S_{b} \end{eqnarray} \]このように、総変動が群内変動 Sw と群間変動 Sb に分解することができた。次に、群内変動 Sw と群間変動 Sb の相対的な大きさを比較する。
Sw と Sb の相対的な大きさを調べたいが、Sw と Sb を計算する際に使用したデータ数が異なっている。データ数が多くなれば、誤差の二乗の和も大きくなるので、このままでは両者を比較できない。そこで、両者をそれぞれのデータ数で割れば、1 データあたりの Sw と 1 データあたりの Sb が求まる。このように補正した Sw/(n - r) と Sb/(r - 1) は直接比較できるようになる。また、このように補正を行うと、両者は分散そのものとみなすことができる。
正規分布の仮定下において、誤差の二乗の和として計算された分散がカイ二乗分布に従う。また、カイ二乗分布に従う 2 つの確率変数の比は、F 分布に従うことが導かれる。このことから、Sw/(n - r) と Sb/(r - 1) に着目すると、Sw/(n - r) は自由度 n - r のカイ二乗分布に従い、Sb/(r - 1) は自由度 r - 1 のカイ二乗分布に従う、そして、両者の比は F 分布に従う。
以上により、Sw/(n - r) と Sb/(r - 1) が同じ値であるかどうかを検定するには F 検定 を用いることができる。つまり、分散分析は、総変動を群内変動と群間変動に分解した上で、群内変動と群間変動の比の分布を利用して、群間に差があるかどうかを検定する方法である。
\[ \frac{\frac{S_{b}}{r-1}}{\frac{S_{w}}{n-r}} \sim \mathcal{F}(r-1, n-r) \]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このデータ構造を R で解析しやすいように、次にように変形する。変更後のデータは 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.44516aov 関数
各実験群の母平均が同じかどうかということは、各個体の値はその実験群の影響をうけているかどうかと解釈することができる。つまり、各個体の値 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 ‘ ’ 1factor(group) の行を見ると、Pr(>F) の値が 2×10-16 と非常に小さくなっている。この分散分析の結果として、A、B、C、D の 4 つの群の母平均が同じではないことがわかる。ただし、分散分析ではどの群の母平均が大きかったか、または小さかったかを判定することはできない。
anova 関数
anova 関数を使うとき、線形モデルを作成する関数と同時に使用する必要がある。手順として、観測値を目的変数として、実験群名を説明変数として線形モデルを構築し、次にその線形モデルを anova 関数に与える。
lm_model <- lm(value ~ factor(group), data = x)
res <- anova(lm_model)
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 ‘ ’ 1aov 関数と anova 関数の違い
aov の実行結果と lm + anova の実行結果がやや異なっている。これは aov がタイプ I の平行和を、lm はタイプ II の平行和を利用して解析しているために起因する。仮に、以下のように lm 関数にダミー変数を与えると、実行結果が aov 関数と一致するようになる。
lm_model <- lm(value ~ 0 + factor(group), data = x)
res <- anova(lm_model)
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