$t$ 検定

はじめに

平均値の検定をする場合、高校では正規分布による検定しか扱いません。例えば、母集団が分散 1 の正規分布だということはわかっているけれども、平均値がわかりません。実際にデータを3個とったところ、値は 1、2、3 でした。平均値が 0 であるという帰無仮説は棄却できるでしょうか。

1、2、3 の平均値は 2 です。分散 1 を仮定すると、3個の平均値の分散は 1/3 です。母集団の平均が 0 だとすると、3個の平均値は平均 0、分散 1/3 の正規分布をするはずです。分散 1/3 ということは、標準偏差 $\sqrt{1/3} \approx 0.577$ です。1、2、3 の平均値 2 は、この標準偏差の $2 / \sqrt{1/3} \approx 3.46$ 倍も離れています。標準偏差の 1.96 倍以上離れる確率は 0.05 ですから、このようなことが起きる確率は 0.05 より小さく、帰無仮説は水準 0.05 で棄却できます。

この論理そのものも難しいですが、平均がわからないのに分散だけわかっているというのは何とも不思議です。現場では平均がわからなければ分散もわからないのが一般的です。

平均も分散もわからないけれども、中心極限定理から、サンプルの平均値はほぼ正規分布をすることがわかっているという一般の場合に、サンプルの平均値が 0 かどうか、あるいは2つのサンプルの平均値が等しいかどうかを検定するのが $t$ 検定です。

こういった統計計算をするための標準的なツールは R ですが、やっぱり Python でも同じことをしたいですよね。Python では SciPy を使うのが一般的でしたが、あまり便利ではありません。以下では Pingouin という新しめのパッケージを使っています。pip install pingouin でインストールできます。念のため SciPy による方法も書いておきました（pip install scipy）。なお、$t$ 分布の密度関数などを求めるには SciPy の scipy.stats.t を使います。

1標本の $t$ 検定

10人の患者にある睡眠薬を飲ませたところ、睡眠時間がそれぞれ次の時間だけ増えました（Arthur R. Cushny and A. Roy Peebles, The Journal of Physiology 32, 501-510 (1905)）：

x = [1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4]

この睡眠薬は効果があるといえるでしょうか。

帰無仮説は、効果はない（上の x の平均値は 0 である）です。10個の平均値は中心極限定理から正規分布で近似できそうですが、分散については何もわかりません。

こういうときに使うのが、1標本の $t$ 検定です（関数の名前は2標本と同じですが、2番目の引数に 0 を指定します）。

import pingouin as pg

pg.ttest(x, 0)

               T  dof alternative     p-val        CI95%   cohen-d    BF10     power
T-test  3.679916    9   two-sided  0.005076  [0.9, 3.76]  1.163692  10.715  0.904364

$p$ 値は 0.005 程度ですので、効果があると言えそうです。睡眠時間の増加の 95% 信頼区間は [0.90, 3.76] 時間です。

念のため、SciPy を使う方法：

from scipy import stats

stats.ttest_1samp(x, 0)

TtestResult(statistic=3.6799158947951884, pvalue=0.005076132649772408, df=9)

95% 信頼区間：

stats.ttest_1samp(x, 0).confidence_interval()

ConfidenceInterval(low=0.897677539412931, high=3.762322460587068)

参考までに、R では次のようになります：

> x = c(1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4)
> t.test(x)

	One Sample t-test

data:  x
t = 3.6799, df = 9, p-value = 0.005076
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.8976775 3.7623225
sample estimates:
mean of x 
     2.33

2標本の $t$ 検定（対応がない場合）

2組の患者たちに痛みのレベルを報告してもらったところ，次のような結果を得ました (T. Lumley, Biometrics 52, 354-361 (1996))。

x = [1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1]
y = [3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4]

これらの平均値に差があるといえるでしょうか。

このような問題では、両群の分散が等しいと仮定しない $t$ 検定（Welch の検定）を使うのが推奨です（R のデフォルトはそうなっています）。

pg.ttest(x, y, correction=True)

               T        dof alternative     p-val           CI95%   cohen-d   BF10     power
T-test -2.948615  15.915569   two-sided  0.009479  [-2.36, -0.38]  1.255397  6.556  0.846951

$p$ 値は 0.009 程度で、両群の平均値に差があると言えそうです。両群の平均値の差の 95% 信頼区間は [-2.36, -0.38] です。

pg.ttest() のデフォルトは、両群のサイズ（個数）が等しければ分散が等しいと仮定し、そうでなければ Welch の方法を用います。オプション correction=True で常に Welch の方法を用います。

さらに教科書によっては、母集団が正規分布であると仮定すべき理由がなければ（あるいはまず正規分布かどうかの検定をして正規分布でなければ）$t$ 検定ではなくノンパラメトリックな方法（Wilcoxon(-Mann-Whitney) の検定など）を使えと書いてあると思います。しかし、ノンパラメトリックな方法も結局は順位和が漸近的に正規分布に近づくことを使ってしまっています。ノンパラメトリックな方法を使わなくても、標本平均は漸近的に正規分布に近づく（中心極限定理）ので、結局は同じ穴のムジナということもできます。さらに、現代においては $p$ 値よりもサンプル間の測定値の差の信頼区間が重視されますので、ノンパラメトリックな方法で $p$ 値だけを出す意味はあまりなさそうです。

SciPy では scipy.stats.ttest_ind を使って次のようにします（関数名の _ind は「独立」independent を意味します）：

from scipy import stats

stats.ttest_ind(x, y, equal_var=False)

TtestResult(statistic=-2.948615301572607, pvalue=0.009478901687159976, df=15.915568793934192)

stats.ttest_ind(x, y, equal_var=False).confidence_interval()

ConfidenceInterval(low=-2.3556088373583686, high=-0.38465090290137105)

statsmodels.stats.weightstats.ttest_ind を使うこともできます：

from statsmodels.stats.weightstats import ttest_ind

ttest_ind(x, y, usevar='unequal')

(-2.9486153015726067, 0.009478901687159983, 15.915568793934193)

参考までに、R では次のようになります：

> x = c(1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1)
> y = c(3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4)
> t.test(x, y)

	Welch Two Sample t-test

data:  x and y
t = -2.9486, df = 15.916, p-value = 0.009479
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.3556088 -0.3846509
sample estimates:
mean of x mean of y 
 1.357143  2.727273

Python で R と似た出力をする関数：

import numpy as np
from scipy import stats

def ttest(x, y):
    n1 = len(x)
    n2 = len(y)
    m1 = np.mean(x)
    m2 = np.mean(y)
    v1 = np.var(x, ddof=1)
    v2 = np.var(y, ddof=1)
    se = np.sqrt(v1 / n1 + v2 / n2)
    d = m1 - m2
    t = d / se
    df = (v1 / n1 + v2 / n2) ** 2 / ((v1 / n1) ** 2 / (n1 - 1) +
                                     (v2 / n2) ** 2 / (n2 - 1))
    p = 2 * stats.t.cdf(-abs(t), df)
    print(f"t = {t:.5g}, df = {df:.5g}, p-value = {p:.4g}")
    e = stats.t.ppf(0.975, df) * se
    print("95 percent confidence interval:")
    print(f" {d - e:.8g} {d + e:.8g}")

ttest(x, y)

t = -2.9486, df = 15.916, p-value = 0.009479
95 percent confidence interval:
 -2.3556088 -0.3846509

2標本の $t$ 検定（対応がある場合）

5人の成績が 1、2、3、4、5 だったのが、勉強したのでそれぞれ 3、5、4、7、6 になりました。勉強によって成績が伸びたと言えるでしょうか？

pg.ttest([1,2,3,4,5], [3,5,4,7,6], paired=True)

               T  dof alternative     p-val           CI95%   cohen-d   BF10     power
T-test -4.472136    4   two-sided  0.011056  [-3.24, -0.76]  1.264911  6.555  0.571609

$p$ 値は 0.01 ほどなので、成績が伸びたと言えそうです。

実は、これは一人一人の成績の伸びを1標本の $t$ 検定したのと本質的に同じです：

pg.ttest([2, 3, 1, 3, 1], 0)

               T  dof alternative     p-val         CI95%  cohen-d   BF10     power
T-test  4.472136    4   two-sided  0.011056  [0.76, 3.24]      2.0  6.555  0.908885

SciPy では stats.ttest_rel を使います（_rel は「関連する」related の意味です）：

stats.ttest_rel([1,2,3,4,5], [3,5,4,7,6])

TtestResult(statistic=-4.47213595499958, pvalue=0.011056493393450068, df=4)

stats.ttest_rel([1,2,3,4,5], [3,5,4,7,6]).confidence_interval()

ConfidenceInterval(low=-3.2416639982037667, high=-0.7583360017962335)

R では次のようになります：

> t.test(c(1,2,3,4,5), c(3,5,4,7,6), paired=TRUE)

	Paired t-test

data:  c(1, 2, 3, 4, 5) and c(3, 5, 4, 7, 6)
t = -4.4721, df = 4, p-value = 0.01106
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -3.241664 -0.758336
sample estimates:
mean difference 
             -2