ウェルチのt検定（概要とpython実装）

t検定について

• 2群間の平均値の差を比較する検定
• 2群間の平均値が独立であり(データに対応がない)、かつ2群間に等分散性が仮定できない場合に用いる。（下図参照）
• 「多重性の問題」を避けるために、最近では等分散性の検定は行わず、等分散性の有無を問わない「ウェルチのt検定」を行った方が良いという考えが一般的になりつつあるらしい。

コード

違いに対応のないデータ、かつ不等分散であり、それぞれ正規分布に従うと仮定。
帰無仮説を「2群の平均値に差がない」と設定し、有意水準0.05で検定を行う。

import numpy as np
from scipy import stats
import seaborn as sns

#用いるデータ（違いに対応のないデータ、かつ不等分散であり、それぞれ正規分布に従うと仮定）
A = np.array([13.8, 10.2, 4.6, 10.0, 4.2, 16.1, 14.4, 4.9, 7.7, 11.4])
B = np.array([3.3, 2.6, 4.0, 4.7, 1.9, 2.9, 4.7, 5.3, 4.3, 3.0, 2.0])

#試しにバイオリンプロットを作成
sns.violinplot(x=['A' for i in range(len(A))]+['B' for i in range(len(B))],
y=np.concatenate([A, B]))
plt.savefig('ウェルチのt検定_バイオリンプロット.png')
result = stats.ttest_ind(A, B, equal_var=False)
print('p-value: ', round(result.pvalue, 4))
# >>> p-value: 0.0012

p < 0.05 なので帰無仮説は棄却され、2群間にの平均値には差があると言える。

検定の方法詳細

手順

1. 仮説の設定
   帰無仮説(H0)：「2群間の平均値に差がない」と仮定する。
   対立仮説(H1)：「2群間の平均値に差がある」と仮定する。
2. 有意水準 $\alpha$ を設定する。
3. 統計量 $t$ を算出する。
4. 自由度 $\nu$ を算出する。（小数点以下を切り捨てて整数にすることが多いとのこと）
5. 自由度 $\nu$ の $t$ 分布を用いて、p値を算出し判定する。
   p > α：帰無仮説は棄却できない（平均値が異なると言えない）
   p < α：帰無仮説を棄却する（平均値が異なると予想できる）

$$ t = \frac{\bar{X_1}-\bar{X_2}}{\sqrt{\frac{s_1^2}{N_1}+\frac{s_1^2}{N_1}}}\\　　$$

$$ \nu \approx \frac{(\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2})^2}{\frac{s_1^4}{n_1^2(n_1-1)}+\frac{s_2^4}{n_2^2(n_2-1)}}$$

$ \bar{X_i}$ ：標本平均
$s_i^2$ ：不偏分散
$n_i$ ：サンプルサイズ

詳細コード

以下の2パターンで検証した

1. t値が棄却域に入るか否かで判定
2. t値からp値を求めて判定

1. t値が棄却域に入るか否かで判定

import matplotlib.pyplot as plt

A = np.array([13.8, 10.2, 4.6, 10.0, 4.2, 16.1, 14.4, 4.9, 7.7, 11.4])
B = np.array([3.3, 2.6, 4.0, 4.7, 1.9, 2.9, 4.7, 5.3, 4.3, 3.0, 2.0])

#標本平均
X1 = A.mean()
X2 = B.mean()

#普遍分散
s1 = A.var(ddof=True)
s2 = B.var(ddof=True)

#サンプル数
n1 = len(A)
n2 = len(B)

#tとνを算出する
t = (X1-X2)/np.sqrt(s1/n1+s2/n2)
nu = int((s1/n1+s2/n2)**2 / (s1**2/(n1**2*(n1-1)) + s2**2/(n2**2*(n2-1))))
print('t:',t)
print('nu:',nu)
"""output
t: 4.426442804187721
nu: 10
"""

# x軸の等差数列を生成
X = np.linspace(start=-6, stop=6, num=100)

# t分布
t_sf = stats.t.pdf(x=X, df=nu)

#パーセント点関数（ppf）でα=0.05の棄却域を計算する
upper_005 = stats.t.ppf(q=0.95, df=nu)
lower_005 = stats.t.ppf(q=0.05, df=nu)

# # 可視化
plt.plot(X, t_sf)
X_lower = np.linspace(start=-6, stop=lower_005, num=100)
X_upper = np.linspace(start=upper_005, stop=6, num=100)
t_lower = stats.t.pdf(x=X_lower, df=nu)
t_upper = stats.t.pdf(x=X_upper, df=nu)
base = np.zeros(100)
plt.fill_between(X_lower, base, t_lower, facecolor='lime', alpha=0.5)
plt.fill_between(X_upper, base, t_upper, facecolor='lime', alpha=0.5)

#求めたt値をプロット
plt.plot(t, 0.01, 's',label='t value')
plt.title('t distribution({})'.format(nu), fontsize=15)
plt.ylim(0,0.5)
plt.legend(fontsize=15)
plt.show()
print('下限点:',round(lower_005,3))
print('上限点:',round(upper_005,3))
"""output
下限点: -1.812
上限点: 1.812
"""

求めた統計量t（4.426）が棄却域に入っているので帰無仮説が棄却され、2群の平均値に差があると予想できる。

2. t値からp値を求めて判定

import matplotlib.pyplot as plt

#自由度と先ほど求めたf値
A_df = 9
B_df = 10
f = 0.771

# x軸の等差数列を生成
A = np.array([13.8, 10.2, 4.6, 10.0, 4.2, 16.1, 14.4, 4.9, 7.7, 11.4])
B = np.array([3.3, 2.6, 4.0, 4.7, 1.9, 2.9, 4.7, 5.3, 4.3, 3.0, 2.0])

#標本平均
X1 = A.mean()
X2 = B.mean()

#普遍分散
s1 = A.var(ddof=True)
s2 = B.var(ddof=True)

#サンプル数
n1 = len(A)
n2 = len(B)

#tとνを算出する
t = (X1-X2)/np.sqrt(s1/n1+s2/n2)
nu = int((s1/n1+s2/n2)**2 / (s1**2/(n1**2*(n1-1)) + s2**2/(n2**2*(n2-1))))
print('t:',t)
print('nu:',nu)
"""output
t: 4.426442804187721
nu: 10
"""

# x軸の等差数列を生成
X = np.linspace(start=-6, stop=6, num=100)

# t分布
t_sf = stats.t.pdf(x=X, df=nu)

#累積確率関数(cdf)と生存関数(sf)でp値を求める。
one_sided_pval1 = stats.t.cdf(t, nu)
one_sided_pval2 = stats.t.sf(t, nu)

# 可視化
plt.plot(X, t_sf)
X_lower = np.linspace(start=-6, stop=t, num=100)
X_upper = np.linspace(start=t, stop=6, num=100)
t_lower = stats.t.pdf(x=X_lower, df=nu)
t_upper = stats.t.pdf(x=X_upper, df=nu)
base = np.zeros(100)

#棄却域を塗りつぶす
plt.fill_between(X_lower, base, t_lower, facecolor='lime', alpha=0.5)
plt.fill_between(X_upper, base, t_upper, facecolor='yellow', alpha=0.5)

#求めたt値をプロット
plt.plot(t, 0.00009, 's',label='t value')
plt.title('t distribution({})'.format(nu), fontsize=15)
plt.ylim(0,0.005)
plt.xlim(3,6)
plt.legend(fontsize=15)
plt.show()
print('p-value(lower side):',round(one_sided_pval1,4))
print('p-value(upper side):',round(one_sided_pval2,4))
"""output
p-value(lower side): 0.9994
p-value(upper side): 0.0006
"""

t値付近を拡大してみてみると...

p値のうち小さい方（正側からt値までの面積 |0.0006|黄色）が、有意水準0.05よりも小さいので帰無仮説は棄却され、2群間の平均値には差があると予想される。

参考

• 2群間での統計検定手法まとめ | qiita
• Welchiのt検定 | 統計WEB
• ウェルチのt検定 | Wikipedia