はじめに
数理統計学の本を読んでいたらWelchの$t$検定が出てきたのですが、そういえばWelchの$t$検定の導出って知らないなぁ、と思いました。本には導出までは書かれていなかったためWikipediaを見てみるとWelchさんの書いた原著論文が参考文献として記載されていたので、その論文を参考に導出をまとめておきます。
母平均の検定
問題設定
2つの正規母集団からのサンプルサイズ$m,n$の独立な標本をそれぞれ
X_1, \ldots,X_m \sim N(\mu_X, \sigma_X^2)
Y_1, \ldots,Y_n \sim N(\mu_Y, \sigma_Y^2)
とします。また、各標本の平均と分散を以下のように定義します。
\bar{X} = \frac{1}{m}\sum_{i=1}^{m}X_i\hspace{3pt}, \hspace{15pt} \bar{Y} = \frac{1}{n}\sum_{i=1}^{n}Y_i
s_X^2 = \frac{1}{m-1}\sum_{i=1}^{m}(X_i-\bar{X})^2 \hspace{3pt}, \hspace{15pt} s_Y^2 = \frac{1}{n-1}\sum_{i=1}^{n}(Y_i-\bar{Y})^2
これらのデータを用いて2つの正規母集団の母平均が異なるかどうかを調べることにします。ここで考える検定問題は
H_0: \mu_X=\mu_Y\hspace{8pt}vs\hspace{8pt}H_1: \mu_X\neq\mu_Y
であり、これを有意水準$\alpha$で検定します。
母分散が既知の場合
\bar{X} - \bar{Y}\sim N(\mu_X-\mu_Y, \frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n})
が成り立つので、帰無仮説$H_0$が正しいという仮定の下で、
Z=\frac{\bar{X} - \bar{Y}}{\sqrt{\frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n}}} \sim N(0,1)
となります。母分散$\sigma_X^2$および$\sigma_Y^2$が既知であれば$Z$の値を実際に計算することができます。そこで、標準正規分布の上側$100α$%点を$z_\alpha$とすると、$|Z|>z_{\alpha/2}$の時に帰無仮説$H_0$を棄却する検定が有意水準$\alpha$の検定になります。
母分散が未知だが等しい場合
母分散が未知の場合は上記の$Z$の値を計算することはできないので、母分散の代わりに母分散の推定値を用いて検定を行うことになります。そこでまずは、
\frac{(m-1)s_X^2}{\sigma_X^2} \sim \chi^2(m-1)\hspace{3pt}, \hspace{15pt}\frac{(n-1)s_Y^2}{\sigma_Y^2} \sim \chi^2(n-1)
のようにそれぞれ自由度$m-1$及び$n-1$の$\chi^2$分布に従う確率変数を作ります。$\chi^2$分布の再生性より、これらの和は
\frac{(m-1)s_X^2}{\sigma_X^2} +\frac{(n-1)s_Y^2}{\sigma_Y^2} \sim \chi^2(m+n-2)
のように自由度$m+n-2$の$\chi^2$分布に従います。一方で、帰無仮説$H_0$が正しいという仮定の下で$Z$は標準正規分布に従うので、$t$分布の定義1より
T=\frac{\frac{\bar{X}-\bar{Y}}{\sqrt{\frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n}}}}{\sqrt{\frac{\frac{(m-1)s_X^2}{\sigma_X^2}+\frac{(n-1)s_Y^2}{\sigma_Y^2}} {m+n-2}}}
は自由度$m+n-2$の$t$分布に従います。ごちゃごちゃした式ですが、母分散が等しい、つまり$\sigma_X^2=\sigma_Y^2$を仮定できる場合は上記の式をきれいに整理することができて
T=\frac{\bar{X}-\bar{Y}}{\sqrt{\frac{(m-1)s_X^2+(n-1)s_Y^2}{m+n-2}}\sqrt{\frac{1}{m}+\frac{1}{n}}} \sim t(m+n-2)
となります。分母の一つ目のルートの中身はよく「プールした分散」とか呼ばれているやつです。このように整理すると$T$には未知母数が含まれていないので実際に値を計算することができます。そこで、自由度$k$の$t$分布の上側$100α$%点を$t_\alpha(k)$とすると、$|T|>t_{\alpha/2}(m+n-2)$の時に帰無仮説$H_0$を棄却する検定が有意水準$\alpha$の検定になります。
母分散について仮定をおけない場合
現実的には母分散について何も仮定を置けない、という状況が最も多いのではないかと思います。しかし、母分散について何も仮定を置けない場合は、上記の$T$の値を計算することはできません。なので2群の母分散が等しいか否かに関わらず母平均の検定を行う方法が欲しくなります。
とりあえず等分散の場合と同様に、母分散は未知なので母分散の代わりに母分散の推定値を用いて検定を行う、という方針に従って、$Z$の母分散$\sigma_X^2,\sigma_Y^2$を母分散の推定値である$s_X^2,s_Y^2$に置き換えた統計量$v$を考えてみます。
v=\frac{\bar{X} - \bar{Y}}{\sqrt{\frac{s_X^2}{m}+\frac{s_Y^2}{n}}}
$Z$は帰無仮説$H_0$が正しいという仮定の下で標準正規分布に従いましたが、$v$はもはや標準正規分布には従いません。検定を作るには$v$の分布を知る必要がありますが、一体どんな分布なのでしょうか？実は統計量$v$は以下の自由度$f$の$t$分布で近似できることが知られています。
f=\frac{\Bigl(\frac{s_X^2}{m}+\frac{s_Y^2}{n}\Bigr)^2} {\frac{\bigl(\frac{s_X^2}{m}\bigr)^2}{m-1}+\frac{\bigl(\frac{s_Y^2}{n}\bigr)^2}{n-1}}
したがって、$|v|>t_{\alpha/2}(f)$の時に帰無仮説$H_0$を棄却する検定が近似的に有意水準$\alpha$の検定になります。これがWelchの$t$検定と呼ばれている方法です。これは一体どこから来たのでしょうか？
Welchのt検定の導出
それではWelchの$t$検定を導出していきましょう。以下の議論では帰無仮説$H_0$が正しいというということを仮定しています。まずは、
Z=\frac{\bar{X} - \bar{Y}}{\sqrt{\frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n}}} \sim N(0,1)
が成り立つので、$\bar{X} - \bar{Y}$は標準正規分布に従う確率変数$Z$を用いて、以下のように書くことができます。
\bar{X}-\bar{Y}=Z\sqrt{\frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n}}
また、
\chi_X^2=\frac{(m-1)s_X^2}{\sigma_X^2} \sim \chi^2(m-1)\hspace{3pt}, \hspace{15pt}\chi_Y^2=\frac{(n-1)s_Y^2}{\sigma_Y^2} \sim \chi^2(n-1)
が成り立つので、$s_X^2$および$s_Y^2$は$\chi^2$分布に従う確率変数を用いてそれぞれ以下のように書くことができます。
s_X^2 = \frac{\chi_X^2\sigma_X^2}{m-1} \hspace{3pt}, \hspace{15pt} s_Y^2=\frac{\chi_Y^2\sigma_Y^2}{n-1}
これらをWelchの$t$検定で用いる統計量$v$に代入して式を整理すると
\begin{align}
v&=\frac{\bar{X} - \bar{Y}}{\sqrt{\frac{s_X^2}{m}+\frac{s_Y^2}{n}}} \\
&= \frac{Z\sqrt{\frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n}}}{\sqrt{\frac{\chi_X^2\sigma_X^2}{m(m-1)}+\frac{\chi_Y^2\sigma_Y^2}{n(n-1)}}} \\
&=\frac{Z}{\sqrt{a\chi_X^2+b\chi_Y^2}}
\end{align}
と書き直すことができます。ここで、$a$と$b$はそれぞれ以下のような値です。
a=\frac{\frac{\sigma_X^2}{m(m-1)}}{\frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n}} \hspace{3pt}, \hspace{15pt}b=\frac{\frac{\sigma_Y^2}{n(n-1)}}{\frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n}}
次に、書き直した統計量$v$の分布を知るために$w=a\chi_X^2+b\chi_Y^2$とおいて、$w$がどんな分布に従うかを調べることにします。$a$と$b$のどちらかが$0$、もしくは$a=b$であれば$w$が$\chi^2$分布を定数倍した分布（つまりガンマ分布）に従うことがすぐにわかります。しかし、それ以外の場合は$w$の分布は簡単な形にはならない（らしい）ので、別の分布で近似しようと思います。じゃあどんな分布で近似するかですが、$w$は$\chi^2$分布の定数倍を複数個足したものなので、1つの$\chi^2$分布の定数倍でいい感じに近似できる気がしてきました。そこで、$w$の分布を形状母数$f/2$、尺度母数$2g$のガンマ分布$Ga(f/2, 2g)$で近似することにしましょう2。ここで問題になるのは、ガンマ分布のパラメータに含まれている$f$と$g$をどうやって決めるのか、という点です。ここで、$w$の正確な分布はわからなくても期待値と分散は以下のように計算できることを利用します。
E[w]=E[a\chi_X^2+b\chi_Y^2]=aE[\chi_X^2]+bE[\chi_Y^2]=a(m-1)+b(n-1)
V[w]=V[a\chi_X^2+b\chi_Y^2]=a^2V[\chi_X^2]+b^2V[\chi_Y^2]=2\{a^2(m-1)+b^2(n-1)\}
そこで、近似に用いる$Ga(f/2, 2g)$の期待値と分散が$w$の期待値と分散に一致するように$f$と$g$を決めることにします。$Ga(f/2, 2g)$の期待値と分散はそれぞれ$gf$、$2g^2f$なので、連立方程式
\begin{align}
&gf=a(m-1)+b(n-1) \\
&2g^2f=2\{a^2(m-1)+b^2(n-1)\}
\end{align}
が成り立つように$f$と$g$を決めます。上記の連立方程式を解くと
\begin{align}
&f=\frac{\{a(m-1)+b(n-1)\}^2}{a^2(m-1)+b^2(n-1)} \\
&g=\frac{a^2(m-1)+b^2(n-1)}{a(m-1)+b(n-1)}
\end{align}
となります。この$f$と$g$を用いて$w$の分布を$Ga(f/2, 2g)$で近似します。
ここで、$w$が$Ga(f/2, 2g)$に近似的に従うとすれば$w/g$という確率変数は自由度$f$の$\chi^2$分布に近似的に従います。$Z$は標準正規分布に従うので、$t$分布の定義より、
\begin{align}
\frac{Z}{\sqrt{\frac{\frac{w}{g}}{f}}}&=\frac{Z}{\sqrt{\frac{a\chi_X^2+b\chi_Y^2}{gf}}}=\frac{Z}{\sqrt{a\chi_X^2+b\chi_Y^2}} \sim t(f)
\end{align}
のように、Welchの$t$検定で用いる統計量$v$が自由度$f$の$t$分布に近似的に従うことが分かりました。あとは、自由度$f$を具体的に求めることができれば検定を構成することができます。分数がたくさん出てきて面倒ですが、上記の$f$の計算式に$a$と$b$の具体的な形を代入して整理すると、
\begin{align}
f=\frac{\{a(m-1)+b(n-1)\}^2}{a^2(m-1)+b^2(n-1)}
=\frac{\Bigl(\frac{\sigma_X^2}{m}+\frac{\sigma_Y^2}{n}\Bigr)^2} {\frac{\bigl(\frac{\sigma_X^2}{m}\bigr)^2}{m-1}+\frac{\bigl(\frac{\sigma_Y^2}{n}\bigr)^2}{n-1}}
\end{align}
と書くことができます。しかし、この式には母分散が含まれているため真の自由度は求めることができません。そこで母分散を母分散の推定量で置き換えた
\hat{f}=\frac{\Bigl(\frac{s_X^2}{m}+\frac{s_Y^2}{n}\Bigr)^2} {\frac{\bigl(\frac{s_X^2}{m}\bigr)^2}{m-1}+\frac{\bigl(\frac{s_Y^2}{n}\bigr)^2}{n-1}}
を$f$の推定量として用いることで自由度を算出します。
ここまでの議論より、帰無仮説$H_0$が正しいという仮定の下で
v=\frac{\bar{X} - \bar{Y}}{\sqrt{\frac{s_X^2}{m}+\frac{s_Y^2}{n}}} \sim t(\hat{f})
が近似的に成立するので、$|v|>t_{\alpha/2}(\hat{f})$の時に帰無仮説$H_0$を棄却する検定が近似的に有意水準$\alpha$の検定になります。これでWelchの$t$検定を導出することができました！
近似・推定量の性質の確認
分布の近似精度の確認
Welchの$t$検定の導出では$w=a\chi_X^2+b\chi_Y^2$の分布をガンマ分布$Ga(f/2, 2g)$で平均と分散が一致するように近似しました。この近似がどれくらいの精度なのか気になったので、$w$から乱数を発生させて作成したヒストグラムに$Ga(f/2, 2g)$の確率密度関数を重ねて雰囲気を見てみることにします。
$m=n$かつ$\sigma_X^2=\sigma_Y^2$であれば$w$の分布と$Ga(f/2, 2g)$は一致するので、気になるのは$m\neq n$、$\sigma_X^2\neq\sigma_Y^2$のケースです。まずは、$m\neq n$の場合を確認してみます。
上のグラフを見ると$m=5, n=10$の場合は$w$の分布と$Ga(f/2, 2g)$は少しずれがあるみたいです。また$m$と$n$が等しくなくてもサンプルサイズが増えると近似精度は良くなっています。次に$\sigma_X^2\neq\sigma_Y^2$の場合を見てみます。
分散が等しくないない場合はちょっと分布がずれているようにも見えますが、大きな乖離はなさそうです。この近似は不等分散については頑健なんですかね。$m\neq n$でサンプルサイズが小さいと近似精度が落ちるようですが、ある程度サンプルサイズが大きければ良く近似できていると考えてよさそうです。
自由度の推定量について
Welchの$t$検定では検定統計量$v$が近似的に従う$t$分布の真の自由度$f$を知ることはできないため、代わりに$f$の推定量である$\hat{f}$を用いて検定を行います。$\hat{f}$がどんな感じの推定量なのか気になりました。
まず、$s_X^2$及び$s_Y^2$はそれぞれ$\sigma_X^2,\sigma_Y^2$に確率収束するので、$\hat{f}$は$f$に確率収束します。したがって、$\hat{f}$は$f$の一致推定量になっていることが分かります。
次に不偏性については期待値を計算するのは面倒そうだったので次のようにシミュレーションで確認してみます。サンプルサイズを固定して、$\sigma_Y^2/\sigma_X^2$の値を変えながら正規乱数を発生させて、$\sigma_Y^2/\sigma_X^2$の値ごとにたくさんの$\hat{f}$を算出します。その$\hat{f}$の平均をとることで$\hat{f}$の期待値を調べ、真の自由度$f$と比較します。自由度そのものを見ても検定にどんな影響があるかよく分からないので、推定した上側2.5%点$t_{0.025}(\hat{f})$と真の上側2.5%点$t_{0.025}(f)$も一緒にプロットすることにします。
グラフを見ると$\hat{f}$の期待値は真の自由度$f$と異なる場合があるため、$\hat{f}$は一様に不偏な推定量ではなさそうです。サンプルサイズが小さい時ほど$\hat{f}$の期待値と真の自由度$f$の乖離は大きく、サンプルサイズが大きくなると乖離が小さくなっていることが確認できます。とはいえ、上側2.5%点を見るとサンプルサイズが小さくても値に大きな差はないですね。
サンプルサイズの比を変えても同じような感じですね。
終わりに
今回のように2群の母分散について仮定を置かずに母平均の検定を行うにはどうすればよいのか、という問題はBehrens–Fisher問題と呼ばれています。Welchの$t$検定はBehrens–Fisher問題の解決策ということですね。ただし、Welchの$t$検定が唯一の方法というわけではなく別の方法も提案されています。Wikipediaに書いてありました。単に母平均の検定といってもだいぶ奥が深いなぁ。
参考文献
- Welch, B. L. “The Significance of the Difference Between Two Means When the Population Variances Are Unequal.” Biometrika, vol. 29, no. 3/4, 1938, pp. 350–62. JSTOR, https://doi.org/10.2307/2332010. Accessed 23 Jul. 2022.
- Welch, B. L. “The Generalization of `Student’s’ Problem When Several Different Population Variances Are Involved.” Biometrika, vol. 34, no. 1/2, 1947, pp. 28–35. JSTOR, https://doi.org/10.2307/2332510. Accessed 23 Jul. 2022.
シミュレーションで用いてプログラム: