Goodな生活

Goodな生活

データサイエンスと弦楽器を探究する

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

分散分析とは

分散分析は、少ないデータで効果的な検証を組み立てる実験計画手法の1つ。医薬生物学や農事試験等で応用される。基本的な考え方は平均値の差の検定と同じである。平均値の差の検定では、2群のデータが平均値が異なるかどうかを検定するため、t検定を行う。分散分析では、3群以上のデータや1つの群に対して2つの要素を含むデータに対し、データのばらつき(分散)を元にF検定を行う。

一元配置

構造式(モデル)

農園ごとにコーヒー豆の収穫量に違いがあるのかどうかを検定する。下表は3つの農園でそれぞれ4回ずつコーヒー豆を収穫した時の収穫量を示している。

f:id:good_na_life:20210629061311p:plain
表1

農園の種類がi、各農園から収穫した回数をjとし、コーヒー豆の収穫量をy_{ij}とする。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=μ+α_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統計量
群間変動 I-1 BSS BSS/(I-1) F = \frac{BSS/(I-1)}{WSS/(N-I)}
郡内変動 N-I WSS WSS/(N-I) -
合計 nI-1 TSS - -

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{14}
\end{align}
}

なお表1では簡単のために、各農園での収穫回数を同一(n_i=n=4)としたが、必ずしも同一である必要はない。農園ごとに収穫回数が異なる場合にも(例えばn_1=3,n_2=4,n_3=2)同様の検定を行うことができる。

f:id:good_na_life:20210629064104p:plain
表2

二元配置分散分析

構造式(モデル)

次に二元配置分散分析のモデルを考える。ある農園には標高の異なる3つの畑があり、2種類の品種のコーヒー豆を栽培している。表3は、標高3パターンと2品種の計6パターンの環境下において、それぞれ3回にわたってコーヒー豆を収穫した結果である。品種と標高の2種類の因子が含まれているため、二元配置と呼ぶ。

f:id:good_na_life:20210629064505p: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{15}
\end{align}
}

α_i,β_jはそれぞれA因子主効果、B因子主効果といい、(αβ)_{ij}を交互作用効果と呼ぶ。
表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{16}
\end{align}
}

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

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

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

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

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

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

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

 {
\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{20}
\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{21}
\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{22}
\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{23}
\end{align}
}

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

次に検定統計量を算出する。
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{24}
\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{25}
\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{26}
\end{align}
}

残差の偏差平方和*4

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

全ての偏差平方和(総変動)*5

 {
\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{28}
\end{align}
}

分散分析表は以下のようになる。

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

帰無仮説と検定統計量

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

 {
\begin{align}
H_α : α_1 = α_2 = \cdots = α_I = 0 \\
H_β:β_1 = β_2 = \cdots = β_J = 0 \\
H_γ:(αβ)_{11} = \cdots = (αβ)_{IJ} = 0 \tag{29}
\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{30}
\end{align}
}

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

繰り返し

コーヒーの収穫回数n_iを繰り返し数と呼ぶ。繰り返しがないとき、つまりn_i=1のときには交互作用の影響は検定できない。\bar{y_{ij}} = \bar{y_i}となるためである。よってS_ε=0となりF統計量を作成できない。

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

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

*5:\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