[2023-06-25] SciPy 1.11.0 が出たのでそれに合わせて修正しました。
平均値の検定をする場合、高校では正規分布による検定しか扱いません。例えば、母集団が分散 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 による方法も書いておきました。なお、$t$ 分布の密度関数などを求めるには SciPy の scipy.stats.t を使います。
Pingouin は、英語の penguin(ペンギン)に似ていますが、19世紀に絶滅したオオウミガラスという海鳥を意味するフランス語です(英語 great auk、学名 Pinguinus impennis)。フランス語の発音は[pɛ̃-ɡwɛ̃](パングアンに近い)ですが、英語圏の人はピングインと読んでいるようです。日本でもピングインのほうが読みやすいと思います。
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] です。
例えば信頼区間だけ欲しい場合には pg.ttest(x, 0)["CI95%"][0]
とします。
念のため、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. 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 の方法を用います。
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
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