先日人に説明する事があって、案外知らない人が多いかも知らないと思い、筆を執ります。2サンプルt検定をN回も行うとtype 1 error の確率が有意水準\( \alpha \)より遥かに大きな値
$$\begin{eqnarray}
1- \left( 1-\alpha \right) ^N
\end{eqnarray}$$
になるのは良く知られた話です。1大量の平均値を比べる為に、分散分析という手法があります。これは変動の分解を元にして行いますが、回帰分析の特別な場合として捉える事が出来ます。
その解説をしようと思います。
分散分析
一番簡単な分散分析は、一元配置法(one way ANOVA)と呼ばれています。実務で得られる感じの規模小さめなデータにも使えるので、有名なのではないでしょうか。
以下のような(実験)データが得られているとします。
データがまとまっている表をそのまま行列と見なすと以下のようになります。
$$\begin{eqnarray}
X=
\begin{pmatrix}
x_{11} &\cdots & x_{1m} \\
x_{21} & \cdots & x_{2m} \\
\vdots & \vdots & \vdots \\
x_{n_1 1} & \cdots & x_{n_m m}
\end{pmatrix}
\end{eqnarray}$$
\(n_i \)で\(i \)列目のデータ数を表しています。必ずしも\(n_1 = \cdots =n_m \)となる必要はありません。全データ数は\(N\)としておきます。
全データの平均を\(\bar{x} \)とすると、全データの変動は以下のように書けます。
$$\begin{eqnarray}
S_T = \sum_{j=1} ^{m} \sum_{i=1}^{n_j} (x_{ij}-\bar{x} )^2
\end{eqnarray}$$
この変動が以下のように分解できるのがミソです。
$$\begin{eqnarray}
S_T &=& S_R +S_E \\
S_R &=& \sum_{j=1} ^{m} \sum_{i=1}^{n_j} ( \bar{x_{j}} -\bar{x} )^2 \\
S_E &=& \sum_{j=1} ^{m} \sum_{i=1}^{n_j} ( x_{ij} – \bar{x_{j}} )^2
\end{eqnarray}$$
ただし、\(\bar{x_j} \)はj列目のデータの平均値を表しています。
同じ列のデータの実験条件が等しいとすると、\(S_R \)は実験条件の変化によって説明できる変動になります。
$$\begin{eqnarray}
R^2= \frac{S_R}{S_T } = 1- \frac{S_E }{S_T}
\end{eqnarray}$$
が大きい時、実験の条件によって、データの値が変化していると思えそうです。
しかし、変動はデータが大きければ無限に大きくなるので、変動だけで判断するのは良くありません。自由度で割って、(不偏)分散で考えた方が良いです。2
$$\begin{eqnarray}
f_T &=& N-1 ,f_R = m-1 ,f_E = N-m \\
V_T &=& S_T /f_T ,V_R = S_R / f_R , V_E =S_E / f_E \\
R’^2&=& \frac{V_R}{V_T } = 1- \frac{V_E }{V_T}
\end{eqnarray}$$
これで目安が得られる訳ですが、統計的に有意かどうかを考えましょう。
列方向のデータは、分散が等しい正規分布から独立に得られているとします。
$$\begin{eqnarray}
x_{ij} \sim \mathcal{N}( \mu _j , \sigma ^{2} )
\end{eqnarray}$$
こうしておくと、
$$\begin{eqnarray}
F= \frac{V_R}{V_E}
\end{eqnarray}$$
は自由度\( (f_R , f_E ) \)のF分布に従います。
このF値と、有意水準を決めた統計的仮説検定から、実験条件の効果はあるのか確かめる事が出来ます。
これが一元配置法でした。
重回帰分析
回帰分析で分散分析を説明出来る事を示すために、簡単に重回帰分析の復習をします。
細かい計算などは、解説記事があるので、読んでみてください。
重回帰分析は、データ\(X\)と予測したい値\(y\)がある時、係数\(\beta \)を使って
$$\begin{eqnarray}
\hat{y}=X\beta
\end{eqnarray}$$
としてyを予測するモデルです。\( \beta \)は最小二乗法で求めます。つまり、
$$\begin{eqnarray}
L= \| y- \hat{y} \| ^2
\end{eqnarray}$$
を最小にする\(\beta \)を求めます。計算を省いて答えを書くと、
$$\begin{eqnarray}
\hat{\beta} = \left( X^{T} X \right) ^{-1} X^T y
\end{eqnarray}$$
となります。
予測したい値\(y\)の変動について、以下の等式が成り立つ事に注意しましょう。
$$\begin{eqnarray}
S_y &=& S_R + S_E \\
S_y&=& \sum ( y_i – \bar{y} )^2 \\
S_R&=& \sum (\hat{y_i} -\bar{y} )^2 \\
S_E &=& \sum (y_i – \hat{y_i} )^2
\end{eqnarray}$$
重回帰分析の精度は、予測したい値\(y \)の持つバラツキを\( \hat{y} \)がどれくらい説明出来るかで推定します。
$$\begin{eqnarray}
R^2 = \frac{S_R}{S_y} = 1- \frac{S_E}{S_y}
\end{eqnarray}$$
この値は、データの列が増えれば勝手に大きくなっていくので、自由度で補正して分散を使った量が大事です。
$$\begin{eqnarray}
R’^2= \frac{V_R}{V_T } = 1- \frac{V_E }{V_T}
\end{eqnarray}$$
\( \beta = ( \beta _0 , \cdots , \beta _m ) \)の内、意味のある係数はどれでしょうか。
この問いは、モデルを
$$\begin{eqnarray}
y&=& \hat{ \beta _0 } +\sum \hat{\beta _i} x_i + \epsilon\\
\epsilon &\sim & \mathcal{N} (0,1)
\end{eqnarray}$$
とすれば良いです。こうしておくことで、
$$\begin{eqnarray}
F=\frac{V_R}{V_E}
\end{eqnarray}$$
はF分布に従います。この量を使って以下の仮説検定を行います。
帰無仮説: \( \beta _i =0 \)
対立仮設: \( \beta _i \neq 0 \)
\( \beta _i =0 \)として、\(F \)の値が小さい時は3、\(\beta _i \)はあっても意味がないという事になります。
逆に、\( \beta _i =0 \)とした時にFが大きくなれば\( \beta _i \) に意味があるという事になります。
回帰分析で分散分析を再現する
これまでの準備を使って、回帰分析で分散分析を再現出来る事を見ましょう。
簡単の為に、1因子、2水準の場合の表を考えます。
表まとめられたデータを以下のよう\(X, y \)に分けられます。
$$\begin{eqnarray}
X= \begin{pmatrix}
1 & 0 \\
\vdots & \vdots \\
1 & 0\\
0 &1 \\
\vdots & \vdots \\
0 &1
\end{pmatrix}
,y= \begin{pmatrix}
x_{11} \\
\vdots \\
x_{1n_1 }\\
x_{21} \\
\vdots \\
x_{2 n2}
\end{pmatrix}
\end{eqnarray}$$
これを使って回帰分析をしましょう。\( \beta \)の推定に必要な行列を計算すると以下のようになります。\(X^T X \)はいつも対称行列になることに注意しましょう。
$$\begin{eqnarray}
X^T X &=& \begin{pmatrix}
N & n_1 &n_2 \\
n_1 &n_1 & 0 \\
n_2 & 0 & n_2
\end{pmatrix} \\
\left( X^T X \right) ^{-1} &=&\frac{1}{N-n_1 -n_2} \begin{pmatrix}
1 & -1 & -1 \\
-1 & \frac{N-n_2 }{n_1 } & -1 \\
-1 & -1 & \frac{N-n_1}{n_2}
\end{pmatrix}
\end{eqnarray}$$
これを使って、\(\beta \)と\(y\)の推定値は以下のようになります。
$$\begin{eqnarray}
\hat{\beta} &=& \begin{pmatrix}
0 \\
\bar{x_1} \\
\bar{x_2}
\end{pmatrix}
,\hat{y}=\begin{pmatrix}
\bar{x_1} \\
\vdots \\
\bar{x_1}\\
\bar{x_2}\\
\vdots \\
\bar{x_2}
\end{pmatrix}
\end{eqnarray}$$
この推定値を使って、変動を計算すると、
$$\begin{eqnarray}
S_y &=& S_R + S_E \\
S_y&=& \sum_{j=1}^{2} \sum_{i}^{n_{j}} ( x_{ij} – \bar{x_{ij} } )^2 \\
S_R&=&\sum_{j=1}^{2} \sum_{i}^{n_{j}}( \bar{x_j}-\bar{x_{ij} } )^2 \\
S_E &=&\sum_{j=1}^{2} \sum_{i}^{n_{j}} (x_{ij} – \bar{x_j} )^2
\end{eqnarray}$$
となり、分散分析を再現します。
回帰分析をすれば良いという事は、交互作用の有意性を計算したければ、モデルに
$$\begin{eqnarray}
x_i x_j
\end{eqnarray}$$
をつけ足せば良いだけという事です。
交互作用については、表ベースで考えると”交互作用の変動”を考えないといけませんが、回帰分析なら、交互作用を見たい列同士を掛けて、データの列方向に付け加えるだけなので、とても楽です。
まとめ
- 分散分析の説明をした
- 回帰分析の説明をした
- 回帰分析で分散分析を再現した
