先日線形回帰分析に関して、基本的でありながら実務家でも間違った理解をしていることが多い点について解説したんですが、t検定についても、間違った情報をもとにした質問をされることがあるので解説します。
とはいえ私も勉強中なので、間違いがあれば指摘いただけるとありがたいです。
ちなみにサンプルサイズ、効果量、検出力についてはいい情報があるので触れていません。本記事ではもっと基本的なトピックを扱います。サンプルサイズ周辺で興味がある方はTJO氏の記事や書籍:サンプルサイズの決め方がおすすめです。特に「サンプルサイズの決め方」はとても良い本だったのでおすすめです。
本記事のトピック
1.正規分布の仮定は必要？不必要？
前提
t検定の解説に「サンプルサイズが大きければ正規性を無視して良い」と書かれていたり、逆に「サンプルサイズが大きいとZ検定で良いので、小さいサンプルでこそt検定」という記述を見かけました。他にも、正規分布を前提として解説している書籍もあるのですが、実際のところどうなのでしょうか？式の成り立ちや実験を通じて確認してみましょう。
結論
先に結論です。
母集団の分布が正規分布から乖離している場合は、その乖離度合いや歪み1の強さによって検定にバイアスがかかるが、サンプルサイズが大きくなると頑健になる。従って、サンプルサイズが大きいときはZ検定でもt検定ほとんど同じ結果になるのでどちらでも良く、小サンプルの場合は正規分布の仮定が必要。
注意点は、乖離度合いや歪みが強ければ強いほど必要なサンプルサイズが増加するので、一概に何件あれば無視して良いとは言えないこと。
解説
実験 ランダムサンプリングでp値の分布を確認
帰無仮説のもとであればどのような仮説検定であってもp値は[0, 1]の一様分布になります。2
逆にもし一様分布でなければその検定にはバイアスがかかってしまうことになります。
この性質を使って、様々な分布で何度も検定を行い p値 がどのように分布するのか確認してみます。
グラフは上段にp値のヒストグラム、下段にp値の累積分布を表示しています。
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats as st
from statsmodels.distributions.empirical_distribution import ECDF
1標本検定
ここからはt検定を何度も行うので、以下の関数を使い回すことにします。
サンプルサイズを変えながら、1度の実験で10,000回t検定を行います。
def ttest_1sample(n, dist):
pvalues = []
for i in range(10000):
x = dist.rvs(size=n)
ttest = st.ttest_1samp(a=x, popmean=dist.mean())
pvalues.append(ttest.pvalue)
return pvalues
def plot(idx, pvalues, n):
plt.subplot(2, 4, idx + 1)
plt.title(f'n={n}')
plt.hist(pvalues, bins=40)
plt.subplot(2, 4, idx + 5)
ecdf = ECDF(pvalues)
plt.plot(ecdf.x, ecdf.y, label='累積経験分布')
plt.plot([0, 1], [0,1], color='gray', linestyle = "dashed", label='累積理論分布')
plt.legend()
def compare_ttest(dist):
plt.figure(figsize=(20, 6))
for i, n in enumerate([10, 30, 100, 200]):
pvalues = ttest_1sample(n, dist)
plot(i, pvalues, n)
正規分布
まずは皆さんご存知正規分布です。理論分布は左右対称の釣鐘型の分布です。
norm = st.norm(loc=0, scale=1)
x = np.linspace(-2, 2, 100)
plt.plot(x , norm.pdf(x));
図 理論分布
この理論分布から何度もランダムサンプリングして仮説検定を行います。
この理論分布は平均が0なので、帰無仮説を0とした1標本t検定を何度も行い、得られるp値の分布を確認します。
norm = st.norm(loc=0, scale=1)
compare_ttest(dist=norm)
図 p値の分布
サンプルサイズを変えながら実験を行いましたが、どのようなサンプルサイズであっても一様に分布しているようですね。
ガンマ分布
次にガンマ分布を使って少し歪みのある分布で実験します。
gamma = st.gamma(a=3)
x = np.linspace(0, 10, 100)
plt.plot(x , gamma.pdf(x));
図 理論分布
また同じように何度も検定を行い、p値の分布を確認してみます。先ほどと同じように理論分布の平均を帰無仮説に設定しています。
サンプルサイズが10で左側に山があり、帰無仮説を棄却しやすくなっています。それでも、サンプルサイズが30ぐらいから偏りがなくなっています。ただ、同じ条件で何度か試すと、サンプルサイズ30件でもバイアスがかかることがありました。みなさんもぜひ試してみてください。
指数分布
次はさらに歪みの強い指数分布です。
expon = st.expon()
x = np.linspace(0, 6, 100)
plt.plot(x , expon.pdf(x));
図 理論分布
また実験します。
サンプルサイズが大きくなるに従って一様分布に近づいているようですが、ガンマ分布よりも偏りが強いようです。
2標本検定
ではここからは2標本検定を実験するので新たに関数を定義します。
なお検定は分散が等しくなくても使えるWelchのt検定を用います。
def ttest_2sample(n, dist1, dist2):
pvalues = []
for i in range(10000):
x = dist1.rvs(size=n)
y = dist2.rvs(size=n)
ttest = st.ttest_ind(x, y, equal_var=False)
pvalues.append(ttest.pvalue)
return pvalues
def compare_ttest(dist1, dist2):
plt.figure(figsize=(20, 3))
for i, n in enumerate([10, 30, 100, 200]):
pvalues = ttest_2sample(int(n/2), dist1, dist2)
plt.subplot(1, 4, i + 1)
plt.title(f'n={n}')
plot(i, pvalues, n);
正規分布 vs 正規分布
norm = st.norm(loc=0, scale=1)
x = np.linspace(-2, 2, 100)
compare_ttest(dist1=norm, dist2=norm)
サンプルサイズによらず一様
ガンマ分布 vs 正規分布
gamma = st.gamma(a=3)
mu = gamma.mean()
norm = st.norm(loc=mu)
x = np.linspace(-2, 10, 100)
plt.plot(x , gamma.pdf(x), label='gamma')
plt.plot(x , norm.pdf(x), label='norm')
plt.legend()
compare_ttest(dist1=gamma, dist2=norm);
理論分布
p値の分布
このぐらいの違いだとほとんど影響がなさそうです。
指数分布 vs 正規分布
expon = st.expon()
mu = expon.mean()
norm = st.norm(loc=mu)
x = np.linspace(0, 10, 100)
plt.plot(x , expon.pdf(x), label='expon')
plt.plot(x , norm.pdf(x), label='norm')
plt.legend();
compare_ttest(dist1=expon, dist2=norm)
理論分布
p値の分布
案外安定していますね。パラメーターを変えるとかなり結果が変わるのですが、個人的には想像したよりも頑健な印象を受けました。
ガンマ分布 vs 指数分布
gamma = st.gamma(a=3)
lam = gamma.mean()
expon = st.expon(scale=lam)
x = np.linspace(0, 10, 100)
plt.plot(x , gamma.pdf(x), label='gamma')
plt.plot(x , expon.pdf(x), label='expon')
plt.legend();
compare_ttest(dist1=expon, dist2=expon)
理論分布
p値の分布
サンプルサイズが小さいときに棄却しにくい方にバイアスがかかっているようです。
実験まとめ
母集団が正規分布の場合はバイアスのない検定ができることがわかりました。
母集団の理論分布によってバイアスのかかり方は異なりますが、歪んだ分布であっても全体的にサンプルサイズが大きくなると安定する傾向が見られました。
定義を確認してみる
実験結果を参考にしながら、t分布の理解を深めるために数式も確認してみます。なおここからの解説は東大本ほぼそのままです。
t分布の定義は...$U$ と $V$ を互いに独立な確率変数とし、$U \sim N(0, 1)$、$V \sim \chi^2(m)$ とする。
このとき、以下の分布を自由度 $m$ の $t$ 分布と言う。
T = \frac{U}{\sqrt{V/m}}
検定に用いる統計量は、
$X_1, X_2,,,X_n$を標本データ、$\bar X$を標本平均として
T = \frac{\bar{X} - \mu}{\sqrt{s^2/n}} , \biggl(s^2 = \frac{1}{n-1}\sum_{i=1}^n(X_i - \bar{X})^2\biggr)
このように計算した値を用いています。$\mu$に帰無仮説となる値を代入すれば、$t$値からp-valueを計算していつも通りに検定を行えばOKですね。2つの式を突合して、
U = \bar{X} - \mu, \\
{\sqrt{V/m}} = {\sqrt{s^2/n}}
のように理解してしまいそうになりますが、グッと堪えてもう少し掘り下げます。
$U \sim N(0, 1)$ ですから $(\bar{X} - \mu)$ ではないんですね。ここで$\bar{X}$の分布を考えると、これは中心極限定理より
\bar{X} \sim N(\mu, \sqrt{\sigma^2/n})
なので、これを標準化すると
\frac{\bar{X} - \mu}{\sqrt{\sigma^2/n}} \sim N(0, 1)
となるので、
U = \frac{\bar{X} - \mu}{\sqrt{\sigma^2/n}}
が正しい定義です。残った$\sqrt{V/m}$については
V = \sqrt{\frac{(n-1)s^2}{\sigma^2}},\\
m = n-1
とすると, $V$ は自由度(n-1)の$\chi^2$分布となり、以下のように$\sqrt{V/m}$を定義できました。
\sqrt{V/m}=\sqrt{\left.\frac{(n-1)s^2}{\sigma^2}\right/(n-1)}
ここまで得られた式を組み合わせてやっと、当初の式に戻ってきます。
T = \left. \frac{\bar{X} - \mu}{\sqrt{\sigma^2/n}}\right/\sqrt{\left.\frac{(n-1)s^2}{\sigma^2}\right/(n-1)} = \frac{\bar{X} - \mu}{\sqrt{s^2/n}}
定義の解釈まとめ
さて、改めてこの式を振り返ると$V$ は$\chi^2$分布なので、ここに正規分布の仮定が入っています。3
よくあるt検定の解説で正規分布を仮定するのはここから来ています。しかし、実験で示した通り正規分布でなくてもサンプルサイズが大きくなると検定のバイアスは低減します。
歪んだ分布は、t分布が標準正規分布に漸近するよりも漸近する「速度」が遅いので、サンプルが少ない間はバイアスがかりますが、t値の$s^2$はサンプルサイズが大きくなると結局母分散$\sigma^2$に漸近するので、バイアスのない検定ができるようになる。というわけですね。
以上 後編へ続きます。
-
歪み:分布の左右非対称の度合い。skewnessとも。 ↩
-
証明が気になる方はこちら「新装改訂版 現代数理統計学」
こちらもにも「現代数理統計学の基礎」 ↩ -
$\chi^2$分布は独立に標準正規分布に従う確率変数の二乗和が従う分布 ↩