【python】scipyでt検定する方法まとめ

概要

　いっっっつも使い方を忘れて調べているので、自分で備忘録を書いておくことにしました。

概要
t検定の概要
使う関数
引数と注意点
やってみる
まとめ
statsmodelsでやる方法

t検定の概要

　t検定と一口に言っても実際はいろいろありますが、今回やるのは2群の標本の平均に差異があるかどうかの検定です。帰無仮説は「両者の平均に差はない」、対立仮説は「両者の平均に差がある」です。

　詳しいことはwikipediaとかを見てください（手抜き）。

t検定 - Wikipedia

使う関数

　scipyのt検定を行う関数としては、

scipy.stats.ttest_ind
scipy.stats.ttest_rel

　の2つがあります。ttest_indは対応のないt検定、ttest_relは対応のあるt検定で使えます。

　使い所が多いのは対応のないt検定を行うttest_indの方なので、こちらだけ取り扱います。

引数と注意点

　いろいろあります。

scipy.stats.ttest_ind(a, b, axis=0, equal_var=True, nan_policy='propagate')

　a,bは普通にデータの入った1次元配列を渡して使うことが多いでしょう。axisという引数があることから想像が付く通り、多次元配列でも渡せるようです。使ったことはありません。

　equal_var=Trueだとスチューデントのt検定、equal=var=Falseだとウェルチのt検定です。これは等分散かどうかに関わらずウェルチのt検定で良いという話題があるので、Falseを指定してやると良いと思います。

等分散か否かに関わらずウェルチの t 検定を使う (べきである) - Qiita

　他の引数はあまり重要ではないので、説明を省略します。

　結果は(t統計量, p値)というtupleっぽいオブジェクト*1で返ります。p値が設定した有意水準（たとえば0.05）より小さいときに有意差があったと言えます（不慣れだと毎回「どっちだっけ？」と思うポイント）。

　また気になる点として、t検定は母集団が正規分布に従うというけっこうきつい仮定を置いています。しかし、実はあまり気にする必要はないという議論もあります。

実際には $X$ が正規分布でなくても， $n$ が大きければ中心極限定理により $\overline{X}$ は正規分布に近づくので，この検定は母集団が正規分布かどうかには鈍感です。データの分布が正規分布かどうかの検定をしてから t検定を行う必要はまったくありません。
t検定

　そういうことらしいです。

　なお、scipyのt検定に片側検定のオプションはありません。両側検定の結果から計算するか、他のライブラリ（statmodelsなど）でやることになります。両側検定の結果から計算する場合は、

t, p = stats.ttest_ind(male, female, equal_var=True)
pval3 = p
pval2 = p / 2.0
pval1 = 1.0 - pval2
if t < 0.0:
    w = pval2
    pval2 = pval1
    pval1 = w
pythonのpandasによる簡単な統計処理：第３回 F検定，t検定その他

　みたいな感じになるようです。

やってみる

　まず適当な確率分布のオブジェクトを生成する。

>>> from scipy import stats
>>> d1 = stats.norm(loc=0, scale=1)
>>> d2 = stats.norm(loc=1, scale=1)

　N(0, 1)とN(1, 1)です。

　標本はこんな感じで取れます。

>>> d1.rvs(10)
array([ 0.18711344,  0.3534579 , -0.52046706,  0.47855615, -0.51033227,
        0.70266909,  0.19253524,  0.28232341,  1.24373963, -0.70771188])

　参考：
scipyで確率分布のサンプルと確率密度関数を生成する - 静かなる名辞

　念のために正しいパラメータになっているか確かめます。

>>> np.mean(d1.rvs(1000))
0.031232377764520782
>>> np.var(d1.rvs(1000))
0.9921694887020086
>>> np.mean(d2.rvs(1000))
0.97389464006143
>>> np.var(d2.rvs(1000))
1.0268883977332324

　大丈夫そうなのでt検定します。有意水準0.05とします（先に決めるのがルールなので・・・）。

　最初は標本サイズ3でやってみます。

>>> a = d1.rvs(3)
>>> b = d2.rvs(3)
>>> a
array([-1.29621283,  0.42129238, -0.13701242])
>>> b
array([ 0.81419163,  1.21399486, -1.40737252])
>>> stats.ttest_ind(a, b, equal_var=False)
Ttest_indResult(statistic=-0.5672127490081906, pvalue=0.6064712381602329)

　pvalue=0.6064712381602329で、0.05より圧倒的におおきいので有意差なしということになります。サンプルサイズが少なすぎて有意差が出せないのです。

　10まで増やしてみます。

>>> c = d1.rvs(10)
>>> d = d2.rvs(10)
>>> stats.ttest_ind(c, d, equal_var=False)
Ttest_indResult(statistic=-2.8617115251996275, pvalue=0.011903486818782736)

　今度は出ました。ただし何回かやると有意になったりならなかったりするので、出方のばらつき次第で変わる可能性があります。

　サンプルサイズの見積もりは以下の方法があるそうです。

$n=16\frac{s^2}{d^2}$
幾つデータが必要か？―平均値の差の検定 | ブログ | 統計WEB
　※引用者注：
　 $s$ は標準偏差、 $d$ は期待される二群間の差
　上式は有意水準5%の設定の場合に80%の検出力になるサンプル数
　（詳細はリンク先で読んでください）

　今回は標準偏差、二群間の差ともに1という簡単な設定なので、16サンプルあれば良い計算です。

>>> sum(stats.ttest_ind(
...       d1.rvs(16), d2.rvs(16), equal_var=False)[1] < 0.05
...       for _ in range(100))
78

　よさそうです。

まとめ

　簡単なことなのですが、割とやり方を忘れやすいので書き記しました。これで今後は忘れないで済むでしょう（フラグ）。

statsmodelsでやる方法

　statsmodelsでもできるので、こちらの記事も参考にしてください。片側検定が簡単にできるなどのメリットがあります。

【python】statsmodelsでt検定する方法 - 静かなる名辞

*1:厳密にはTtest_indResultという型だが、tupleの派生クラスで実質的にtupleとして扱えるので気にしなくて良い。中身はともにfloat