Goodな生活

経済学修士→環境コンサル→データサイエンス

実験計画法(2)分散分析

分散分析とは

分散分析は実験計画法の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の比であり、不偏分散の比と同様にF分布を用いて検定することができる。

帰無仮説のもとで
 {
\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分布に従う。


f:id:good_na_life:20210908224228p:plain

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

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

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

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

 {
\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:この式は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