実験計画法（２）分散分析

分散分析とは
一元配置分散分析
二元配置分散分析
参考文献

分散分析とは

分散分析は実験計画法の1つであり、実験で得られたデータより、処置（因子や水準）が結果に与える効果の検定を行うもの。基本的な考え方は平均値の差の検定と同じ。平均値の差の検定では、2群（処置あり、なし）のデータの平均値が等しいという帰無仮説についてt検定を行うが、分散分析では、因子や水準ごとのデータのばらつき（分散）を使いF検定を行う。帰無仮説や検定方法については後述する。

一元配置分散分析

一つの因子について複数の水準下での実験データを比較する分析を一元配置と呼ぶ。ここではある大きな農園におけるコーヒーの収穫量を最大化する問題を考える。コーヒー豆の品種Aを因子、3段階の水準（標高1000m、標高1200m、標高1500m）を設定する。

下表は、それぞれの標高で4回ずつ収穫を行った際に得られた収穫量データである。

f:id:good_na_life:20210829125952p:plain — 表1

構造式（モデル）

標高の水準を $i$ 、収穫の回数を $j$ （回目）とすると、コーヒーの収穫量は $y_{ij}$ となる。標高の水準 $i$ によって $y_{ij}$ の平均値に違いがあるかどうかを検定したいので、まず $y_{ij}$ の標高ごとの平均を $θ_i$ とする。一元配置モデルは

${ \begin{align} y_{ij} = θ_i + ε_{ij} , \quad i &= 1,\cdots ,I, \quad j = 1,\cdots, n_i, \\ ε_{ij} \sim \mathcal{N}(0,σ^2) \tag{1} \end{align} }$

表1では標高が3水準なので $I=3$ 、収穫回数がどの標高においても4回なので $n_i = n= 4$ 、実験データ数は $N=nI=12$ となる。
標高（水準）ごとの平均 $θ_i$ の平均を $μ$ （一般平均）とする。

${ \begin{align} μ = \frac{1}{I} \sum_{i=1}^{I} θ_i \tag{2} \end{align} }$

標高の違いが収穫量に与える効果を $α_i =θ_i - μ$ 、つまり $θ_i$ と一般平均 $μ$ との差とすると*1、(1)は

${ \begin{align} y_{ij} &= μ + α_i + ε_{ij} \tag{3} \end{align} }$

と表される。(3)のモデルを推定し、1つでも $α_i≠0$ となる $α_i$ が得られれば標高によって $y_{ij}$ の平均は異なることが示される。同様に、 $y_{ij}$ のデータ全体の平均を $\bar{y}$ とし、標高ごとの平均 $\bar{y_i}$ は以下のように表すことができる。

${ \begin{align} \bar{y} &= \frac{1}{N} \sum_{i=1}^{I} \sum_{j=1}^{n_i} y_{ij}\\ &= \frac{1}{I} \sum_{i=1}^{I} \bar{y_i}\tag{4} \end{align} }$

${ \begin{align} \bar{y_i} &= \frac{1}{n} \sum_{j=1}^{n_i} y_{ij} \tag{5} \end{align} }$

とおくと、 $μ,θ_i,α_i$ は次のように推定される。

${ \begin{align} μ \approx \bar{y}, \quad θ_i \approx \bar{y_i}, \quad α_i \approx \bar{y_i} - \bar{y} \tag{6} \end{align} }$

帰無仮説

帰無仮説は「標高によってコーヒー豆の収穫量の平均値に差がない」ことなので、(1)において

${ \begin{align} H_0 : θ_1 = θ_2 = \cdots = θ_I \tag{7} \end{align} }$

もしくは(3)において

${ \begin{align} H_0 : α_1 = α_2 = \cdots = α_I = 0 \tag{8} \end{align} }$

が成立するかどうかを調べればよい。

(8)を調べるには、パラメータ $α_i$ について

${ \begin{align} \sum_{i=1}^{I} α_i^2 = 0 \tag{9} \end{align} }$

が成立するかどうか、推定値を用いて表すと、

${ \begin{align} \sum_{i=1}^{I} (\bar{y_i} - \bar{y})^2 = 0 \tag{10} \end{align} }$

標高の水準 $i$ の平均値と全体の平均値の差の平方和がゼロであるかどうかを調べればよい。

検定統計量の算出

この検定を行うために、平方和の分解を考える。データ全体の偏差平方和をTSSとし、その内訳としてBSS（群間変動）、WSS（郡内変動）を定義する。ここでの変動は偏差平方和と同義。

${ \begin{align} TSS &= \sum_{i=1}^{I} \sum_{j=1}^{n_i} (y_{ij} - \bar{y})^2　\\ BSS &= \sum_{i=1}^{I} n_i(\bar{y_{i}} - \bar{y})^2 \\ WSS &= \sum_{i=1}^{I} \sum_{j=1}^{n_i} (y_{ij} - \bar{y_i})^2 \tag{11} \end{align} }$

帰無仮説の $\sum_{i=1}^{I}(\bar{y_{i}} - \bar{y})^2$ が0から離れるにつれてBSS（群間変動）の値は大きくなるため、検定にBSS（群間変動）を利用することができる。次のような検定統計量を用いる。

${ \begin{align} F = \frac{BSS/ (I-1)} {WSS/ (N-I)} \tag{12} \end{align} }$

(12)は群間変動と郡内変動の平均平方*2の比であり、不偏分散の比と同様にＦ分布を用いて検定することができる。

帰無仮説のもとで
${ \begin{align} F \sim F(I-1,N-I) \tag{13} \end{align} }$

$F$ の値が $F_{(I-1,N-I,α)}$ を超えるとき（十分に値が大きいとき）は帰無仮説は棄却される。これらの値を一覧表にまとめたものが分散分析表である。

	平方和	自由度	平均平方	F統計量
水準（標高）	$BSS$	$I-1$	$BSS/(I-1)$	$F = \frac{BSS/(I-1)}{WSS/(N-I)}$
誤差	$WSS$	$N-I$	$WSS/(N-I)$	-
合計	$TSS$	$nI-1$	-	-

長々と書いてしまったが、一元配置分散分析は、BSS（群間変動）とWSS（郡内変動）の比を求め、BSS（群間変動）が十分に大きければ水準を変えたことによる効果がある、と結論づける。構造式の誤差項が正規分布に従うと想定しており、正規分布に従う確率変数の二乗であるBSSとWSSはカイ二乗分布に従う。したがってBSSとWSSの比はF分布に従う。

水準ごとの平均値の点推定

分散分析表を用いて、(6)に書いたように、標高ごとの平均 $θ_i$ の近似値を推定することが可能である。点推定値は、

${ \begin{align} \hat{θ}_i = \bar{y}_i \tag{14} \end{align} }$

となり、 $\bar{y}_i$ の95%信頼区間は、群内変動（誤差）の平均平方（分散）、自由度 $N-1$ の $t$ 分布表を用いて、以下に表す。

${ \begin{align} \bar{y}_i \pm t_{α = 0.025, df={N-1}} \sqrt{\frac{1}{N}\frac{WSS}{N-I}} \tag{15} \end{align} }$

（補足）TSS = BSS + WSSの証明

${ \begin{align} TSS &= \sum_{i=1}^{I} \sum_{j=1}^{n} (y_{ij} - \bar{y})^2 \\ &= \sum_{i=1}^{I} \sum_{j=1}^{n} (y_{ij} - \bar{y_i} + \bar{y_i} - \bar{y})^2 \\ &= \sum_{i=1}^{I} \sum_{j=1}^{n} \left\{(y_{ij} - \bar{y_i})^2 + (\bar{y_i} - \bar{y})^2 +2(y_{ij} - \bar{y_i})(\bar{y_i} - \bar{y}) \right\} \\ &= \sum_{i=1}^{I} \sum_{j=1}^{n} (y_{ij} - \bar{y_i})^2 + \sum_{i=1}^{I} \sum_{j=1}^{n} (\bar{y_i} - \bar{y})^2 + 2 \sum_{i=1}^{I} (\bar{y_i} - \bar{y}) \sum_{j=1}^{n} (y_{ij} - \bar{y_i}) \\ &= \sum_{i=1}^{I} \sum_{j=1}^{n} (y_{ij} - \bar{y_i})^2 + \sum_{i=1}^{I} \sum_{j=1}^{n} (\bar{y_i} - \bar{y})^2 + 0 \\ &= BSS + TSS \tag{16} \end{align} }$

表1では簡単のために、各標高での収穫回数を同一（ $n_i=n=4$ ）としたが、必ずしも同一である必要はない。標高によって収穫回数*3が異なる場合にも（例えば表2では $n_1=3,n_2=4,n_3=2$ ）同様の検定を行うことができる。

f:id:good_na_life:20210829131324p:plain — 表2

二元配置分散分析

次に、品種と標高の2つの因子の組み合わせによる二元配置の実験を行う。1つの因子（A因子）はコーヒーの品種（アラビカ種かティピカ種）、もう1つの因子（B因子）は標高（1,000m、1,200m、1,500m）。先ほどと同じ農園で、標高の異なる3つの場所に、2種類の品種を栽培する。表3は標高3パターンと2品種の計6パターンの環境下において、それぞれ3回にわたってコーヒー豆を収穫した結果である。

f:id:good_na_life:20210829132319p:plain — 表3

構造式（モデル）

コーヒーの品種を $i$ 、標高を $j$ 、収穫の回数を $k$ （回目）とすると、コーヒーの収穫量 $y_{ijk}$ の二元配置モデルは

${ \begin{align} y_{ijk} = μ + α_i + β_j + (αβ)_{ij} + ε_{ijk}, \quad &i = 1,\cdots ,I, \quad j = 1,\cdots, J, \quad k = 1, \cdots , n_{ij} \\ & ε_{ijk} \sim \mathcal{N}(0,σ^2) \tag{17} \end{align} }$

$α_i$ はA因子の主効果、 $β_j$ はB因子の主効果といい、 $(αβ)_{ij}$ をA因子とB因子の交互作用効果と呼ぶ。
表3では $I=2, J=3, n_{ij} = n = 3$ であり、
全体のデータ数は $N = \sum_{i=1}^{I} \sum_{j=1}^{J} n_{ij} = nIJ=18$ となる。

$α_i,β_j,(αβ)_{ij}$ は

${ \begin{align} \sum_{i=1}^{I}α_i = \sum_{i=1}^{J} β_j = \sum_{i=1}^{I} (αβ)_{ij} = \sum_{j=1}^{J} (αβ)_{ij} = 0 \tag{18} \end{align} }$

を満たすパラメータである。検定にあたっては「要因A（品種）に影響力があるなら $α_i$ にゼロではないものがある」と考え、

${ \begin{align} \sum_{i=1}^{I}α_i^2 >0 \tag{19} \end{align} }$

同様に、「要因B（標高）に影響力があるなら $β_j$ にゼロではないものがある」と考えると、

${ \begin{align} \sum_{j=1}^{J}β_j^2 >0 \tag{20} \end{align} }$

「要因A（品種）、B（標高）に相乗効果（もしくは相殺効果）があるなら $(αβ)_{ij}$ にゼロではないものがある」と考えると

${ \begin{align} \sum_{j=1}^{J}(αβ)_{ij}^2 >0 \tag{21} \end{align} }$

ここで $\bar{y},\bar{y_i},\bar{y_j},\bar{y_{ij}}$ を以下のように定義すると、 $μ,α_i,β_j,(αβ)_{ij}$ の推定値を表すことができる*4。

${ \begin{align} \bar{y} &= \frac{1}{nIJ} \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n} y_{ijk} \\ &= \frac{1}{nIJ} \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n} (μ + α_i + β_j + (αβ)_{ij} + ε_{ijk}) \\ &= μ + \frac{1}{nIJ} \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n} (α_i + β_j + (αβ)_{ij} + ε_{ijk}) \\ &= μ \tag{22} \end{align} }$

したがって $μ \approx \bar{y}$ 。

${ \begin{align} \bar{y_i} &= \frac{1}{nJ}\sum_{j=1}^{J} \sum_{k=1}^{n} y_{ijk} \\ &= \frac{1}{nJ} \sum_{j=1}^{J} \sum_{k=1}^{n} (μ + α_i + β_j + (αβ)_{ij} + ε_{ijk}) \\ &= μ + α_i + \frac{1}{nJ} \sum_{j=1}^{J} \sum_{k=1}^{n} (β_j + (αβ)_{ij} + ε_{ijk}) \\ &= μ + α_i \tag{23} \end{align} }$

したがって $α_i \approx \bar{y_i} - \bar{y}$ 。

${ \begin{align} \bar{y_j} &= \frac{1}{nI}\sum_{i=1}^{I} \sum_{k=1}^{n} y_{ijk} \\ &= \frac{1}{nI} \sum_{i=1}^{I} \sum_{k=1}^{n} (μ + α_i + β_j + (αβ)_{ij} + ε_{ijk}) \\ &= μ + β_j+ \frac{1}{nI} \sum_{i=1}^{I} \sum_{k=1}^{n} (α_i + (αβ)_{ij} + ε_{ijk}) \\ &= μ + β_j \tag{24} \end{align} }$

したがって $β_j \approx \bar{y_j} - \bar{y}$ 。

${ \begin{align} \bar{y_{ij}} &= \frac{1}{n} \sum_{k=1}^{n} y_{ijk} \\ &= \frac{1}{n} \sum_{k=1}^{n} (μ + α_i + β_j + (αβ)_{ij} + ε_{ijk}) \\ &= μ + α_i +β_j + (αβ)_{ij} \tag{25} \end{align} }$

したがって $(αβ)_{ij} \approx \bar{y_{ij}} - \bar{y_i} - \bar{y_i} + \bar{y}$ 。

分散分析表

次に分散分析表を作成し、F検定のための検定統計量を算出する。
A因子の偏差平方和は

${ \begin{align} S_{α} &= \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n} (\bar{y_{i}} - \bar{y})^2 \\ &= nJ \sum_{i=1}^{I} (\bar{y_{i}} - \bar{y})^2 \tag{26} \end{align} }$

B因子の偏差平方和は

${ \begin{align} S_{β} &= \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n} (\bar{y_{j}} - \bar{y})^2 \\ &= nI \sum_{i=1}^{I} (\bar{y_{j}} - \bar{y})^2 \tag{27} \end{align} }$

交互作用の偏差平方和は

${ \begin{align} S_{αβ} &= \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n} (\bar{y_{ij}} - \bar{y_i} -\bar{y_j} + \bar{y})^2 \\ &= n \sum_{i=1}^{I} \sum_{j=1}^{J} (\bar{y_{ij}} - \bar{y_i} -\bar{y_j} + \bar{y})^2 \tag{28} \end{align} }$

残差の偏差平方和*5は

${ \begin{align} S_{ε} &= \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n} (\bar{y_{ijk}} - \bar{y_{ij}} )^2 \tag{29} \end{align} }$

全ての偏差平方和（総変動）*6は

${ \begin{align} S_{T} &= \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{n} (y_{ijk} - \bar{y})^2\\ &= S_α + S_β + S_{αβ} + S_ε \tag{30} \end{align} }$

したがって分散分析表は以下のようになる。

	平方和	自由度	平均平方	F統計量
A因子	$S_α$	$I-1$	$\frac{S_α}{(I-1)}$	$F_α = \frac{S_α}{(I-1)} / \hat{σ}^2$
B因子	$S_β$	$J-I$	$\frac{S_β}{(N-I)}$	$F_β = \frac{S_β}{(N-I)} / \hat{σ}^2$
交互作用	$S_{αβ}$	$(I-1)(J-1)$	$\frac{S_{αβ}}{(I-1)(J-1)}$	$F_{αβ} =\frac{S_{αβ}}{(I-1)(J-1)} /\hat{σ}^2$
残差	$S_{ε}$	$(n-1)IJ$	$\hat{σ}^2 = \frac{S_{ε}}{(n-1)IJ}$	-
合計	$S_T$	$nIJ-1$	-	-	-

帰無仮説と検定統計量

帰無仮説はA因子（品種）については「品種の違いによって収穫量の平均値に差がない」、B因子（標高）については「標高の違いによって収穫量の平均値に差がない」、交互作用については「品種と標高における相互作用（相殺作用）による収穫量への影響はない」となる。
帰無仮説のもとで

${ \begin{align} H_α : α_1 = α_2 = \cdots = α_I = 0 \\ H_β：β_1 = β_2 = \cdots = β_J = 0 \\ H_γ：(αβ)_{11} = \cdots = (αβ)_{IJ} = 0 \tag{31} \end{align} }$

F検定統計量は
${ \begin{align} F_α &\sim F(I-1,(n-I)IJ) \\ F_β &\sim F(J-1,(n-1)IJ) \\ F_{αβ} &\sim F((I-1)(J-1),(n-1)IJ) \tag{32} \end{align} }$

検定から分かることは、帰無仮説を棄却できるかどうか、すなわち品種や標高によって収穫量に差があるのかどうかだけである。ここではどの品種や標高の組み合わせによって差が生じるのかまでは分からない。帰無仮説が棄却されたときに、もう少し細かい枠組みでの検定を設定し、その組み合わせによって差が生じるのかを調べることもできる。こうした検定を多重比較検定と呼ぶ。

参考文献

日本統計学会公式認定統計検定準1級対応統計学実践ワークブック

学術図書出版社

Amazon

現代数理統計学の基礎 (共立講座数学の魅力)

作者:達也, 久保川
共立出版

Amazon

29-1. 分散分析とは | 統計学の時間 | 統計WEB

*1:この式は $I$ 個のパラメータを $I+1$ 個のパラメータで表現いるため一意ではない。したがって $θ_i,α_i$ が推定できない。推定する（識別可能とする）ために $\sum_{i=1}^{I}α_i=0$ または $\sum_{i=1}^{I}n_{i} α_{i}=0$ と制約をおく。前者の制約下では $μ$ は $I$ 個の $θ_i$ の平均、後者の制約下では $θ_i$ の加重平均として推定される。

*2:平均平方は偏差平方和を自由度で割った値であるため標本分散と同義。

*3:コーヒーの収穫回数 $n_i$ を繰り返し数と呼ぶ。二元配置分析において繰り返しがないとき、つまり $n_i=1$ のときには交互作用の影響は検定できない。 $\bar{y_{ij}} = \bar{y_i}$ となり、 $S_ε=0$ つまりF統計量を作成できない。

*4:F検定の結果、交互作用 $(αβ)_{ij}$ が認められない場合、A因子とB因子それぞれの収穫量の平均値の高い同士を最適な水準として選択する。この場合の点推定値は $\bar{y_i*} + \bar{y_j*} - \bar{y}$

*5:実際の計算ではA因子、B因子、交互作用、総変動を計算した後、総変動から残差以外の平方和の差分を求めることで導出できる

*6: $\sum_{i} \sum_{j} \sum_{k} y^2_{ijk} - nIJ(\bar{y})^2$ としても計算できる。これは分散を求めるときに頻繁に使うテクニック。 $V[x]=\frac{1}{n}\sum_{i}(x_i-\bar{x}) = \frac{1}{n} (\sum_{i}x^2 - n\bar{x}^2)=E[x^2] - (E[x])^2$