対応のあるt検定について
- 2群間の平均値の差を比較する検定(下図参照)
- 2群間の平均値が独立とはいえない、つまりデータに対応がある(従属である)場合に用いる。
(例:ある患者たちの1ヶ月の体重の推移をみる場合など)
コード
Rのデフォルトデータセット「sleep」を用いる。
(薬剤1と薬剤2の2種類の睡眠薬を10人の被験者にそれぞれ投与した際に、増加した睡眠時間を記録したもの)
2種類の薬剤を同一被験者に投与しているので、その値には対応があるといえる。
データは正規性がある(正規分布している)と仮定する。
有意水準を0.05、帰無仮説を「2群のデータに差がないこと」として検定する。
import numpy as np
from scipy import stats
import seaborn as sns
#用いるデータ(それぞれ正規分布に従うと仮定)
A = np.array([0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0])
B = np.array([1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4])
print('平均値(A):',round(A.mean(),2))
print('平均値(B):',round(B.mean(),2))
"""output
平均値(A): 0.75
平均値(B): 2.33
"""
#試しにバイオリンプロットを作成
sns.violinplot(x=['A' for i in range(len(A))]+['B' for i in range(len(B))], y=np.concatenate([A, B]))
result = stats.ttest_rel(A, B)
print('p-value: ', round(result.pvalue, 4))
# >>> p-value: 0.0028
p < 0.05なので帰無仮説は棄却され、2群のデータ間には差があるといえる。
検定の方法詳細
手順
- 仮説の設定
帰無仮説(H0):「2群間のデータに差がない」と仮定する。
対立仮説(H1):「2群間のデータに差がある」と仮定する。 - 有意水準 $\alpha$ を設定する。
- 統計量$t$を算出する。
- 自由度 $\nu$ を算出する。
- 自由度 $\nu$ の $t$ 分布を用いて、p値を算出し判定する。
p > $\alpha$ :帰無仮説は棄却できない(データに差があると言えない)
p < $\alpha$ :帰無仮説を棄却する(データに差があると予想できる)
$$ t = \frac{\bar{d}-\mu}{\sqrt{\frac{s^2}{n}}} $$
$$ \nu = n-1 $$
$\bar{d}$:差の平均
$\mu $:差の母平均(差が0かどうかを検定する場合は0とおく)
$s^2$:差の不偏分散
$n$:サンプルサイズ
詳細コード
以下の2パターンで検証した
- t値が棄却域に入るか否かで判定
- t値からp値を求めて判定
1. t値が棄却域に入るか否かで判定
import numpy as np
from scipy import stats
import seaborn as sns
#用いるデータ(それぞれ正規分布に従うと仮定)
A = np.array([0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0])
B = np.array([1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4])
#サンプルサイズ
n = len(A)
#それぞれの差を計算する
delta = np.array([a-b for a, b in zip(A,B)])
#差が0か否かを検定するため μ=0とおく
mu = 0
#差の平均
d = delta.mean()
#差の普遍分散s^2
s = delta.var(ddof=True)
#統計量t と 自由度ν を算出する
t = (d - mu)/np.sqrt(s/n)
nu = n -1
print('t:',round(t,3))
print('nu:',nu)
"""output
t: -4.062
nu: 9
"""
# 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.833
上限点: 1.833
"""
求めた統計量t(-4.062)が棄却域に入っているので帰無仮説が棄却され、2群のデータに差があると予想できる。
2. t値からp値を求めて判定
import numpy as np
from scipy import stats
import seaborn as sns
#用いるデータ(それぞれ正規分布に従うと仮定)
A = np.array([0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0])
B = np.array([1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4])
#サンプルサイズ
n = len(A)
#それぞれの差を計算する
delta = np.array([a-b for a, b in zip(A,B)])
#差が0か否かを検定するため μ=0とおく
mu = 0
#差の平均
d = delta.mean()
#差の普遍分散s^2
s = delta.var(ddof=True)
#統計量t と 自由度ν を算出する
t = (d - mu)/np.sqrt(s/n)
nu = n -1
print('t:',round(t,3))
print('nu:',nu)
"""
t: -4.062
nu: 9
"""
# 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.5)
plt.xlim(-6,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))
"""
p-value(lower side): 0.0014
p-value(upper side): 0.9986
"""
t値付近を拡大してみてみると...
p値のうち小さい方(負側からt値までの面積 | 0.0014 | 黄緑)が、有意水準0.05よりも小さいので帰無仮説は棄却され、2群間のデータには差があると予想される。
参考
リンク
リンク