検定

対応のあるt検定(概要とpython実装)

対応のある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群のデータ間には差があるといえる。

バイオリンプロット

検定の方法詳細

手順

  1. 仮説の設定
    帰無仮説(H0):「2群間のデータに差がない」と仮定する。
    対立仮説(H1):「2群間のデータに差がある」と仮定する。
  2. 有意水準 $\alpha$を設定する。
  3. 統計量$t$を算出する。
  4. 自由度 $\nu$を算出する。
  5. 自由度 $\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パターンで検証した

  1. t値が棄却域に入るか否かで判定
  2. 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群間のデータには差があると予想される。

参考

 

 

ottantacinqueをフォローする
実践ケモインフォマティクス

コメント