Python
scipy
初心者
医療統計
BiSH

Pythonで行う基礎医学研究のための統計解析法

More than 1 year has passed since last update.

皆さんが普段一番よく使う統計学的検定法は何でしょうか?
私の場合は「t検定と多重比較検定」の2つが圧倒的な割合を占めています。

臨床研究を行うには他にも様々な手法や知識が必要ですが、私のように10匹前後のマウスを使って2-3群で比較を行うだけならば、この2つを押さえておけば当座をしのげます。
この記事ではPythonおよび数値解析用ライブラリであるSciPy(あるいはNumPy)を用いて、当座をしのぐための解析方法を簡単に紹介したいと思います。

Welchのt検定(独立した2群の平均値の比較)

"マウスに異なる薬剤を投与した時の血圧の変化を調べる実験をした"と仮定しましょう。
A薬を投与した8匹の血圧が「100, 96, 80, 92, 101, 99, 89, 93(mmHg)」であり、
B薬を投与した9匹の血圧が「83, 77, 86, 71, 70, 98, 80, 75, 74(mmHg)」だったとします。
A薬群の方が血圧が高そうに見えますが、それが偶然かどうかを検定してみましょう。

>>> from scipy import stats
>>> A = [100, 96, 80, 92, 101, 99, 89, 93]
>>> B = [83, 77, 86, 71, 70, 98, 80, 75, 74]
>>> stats.ttest_ind(A, B, equal_var = False) #等分散を仮定せずに検定

Ttest_indResult(statistic=3.770683261231587, pvalue=0.0018828551658725231)

このように、SciPyを使って簡単に結果(t値とp値)が得られました。
p値(=0.00188) が有意水準である0.05 (あるいは0.01)を下回りますので、(A薬群とB薬群の平均に差が無いという)帰無仮説が棄却されます。どうやらB薬を投与した方が血圧が有意に低くなるようです。

第3引数を省略すると等分散を前提とした検定(Studentのt検定)ができますが、等分散性に関わらずWelchのt検定をデフォルトで使用すべきとの記事がありました。
なお、t値とp値は下記のように返り値からunpack代入することもできます。

>>> t, p = stats.ttest_ind(A, B, equal_var = False)
>>> print("t_value is", round(t, 3), "and p_value is", round(p, 3))
t_value is 3.771 and p_value is 0.002

ちなみにfrom文を使わずに呼び出す場合、import scipy だけでは scipy.stats は呼び出せません。import scipy.stats と書かなければいけないので、初心者の方は注意しましょう。

一元配置分散分析法(One Way ANOVA)

先ほどのA薬群・B薬群のほかに、C薬群も加えてみましょう。
7匹のマウスの血圧が「90, 91, 92, 89, 88, 90, 91(mmHg)」だったとして、今度はこれら3群の血圧に有意な差が存在するかどうかを検定します。

>>> from scipy import stats
>>> A = [100, 96, 80, 92, 101, 99, 89, 93]
>>> B = [83, 77, 86, 71, 70, 98, 80, 75, 74]
>>> C = [90, 91, 92, 89, 88, 90, 91]
>>> stats.f_oneway(A, B, C)

F_onewayResult(statistic=10.437458422647644, pvalue=0.00071250192100598777)

やり方はt検定とほぼ同じですね。p値が0.0007ということで帰無仮説が棄却され、これら3群間に有意な差が存在することが分かります。比較する組み合わせは"A vs B"、"A vs C"、"B vs C"の3通りがありますから、これらのどの組み合わせに有意差が存在するのかを調べるため、引き続いて多重比較検定を行いましょう。

...その前に

厳密な話をすると、本来は分散分析の前に「各群の分散が等しいかどうか」をBartlett検定で確認しなくてはいけません。じつは本作例は母分散が等しくないので、ANOVAではなくKruskal-Wallis検定を行うのが正解です。そしてKruskal-Wallis検定で水準間に差が認められた場合にはSteel-Dwass法などのノンパラメトリック多重比較検定を行うべきなのですが、今回はPythonを活用することが主眼ですので何卒ご容赦下さい。

Tukey-Kramerの検定(多重比較検定)

Tukeyの検定はパラメトリック多重比較検定の一種で、分散分析で水準間に有意差が見られた場合に、各水準が等分散で正規分布に従っていると仮定して検定します。更に、各群のサンプル数が異なっていても比較できるように改良されたのがTukey-Kramer法だそうです。

PythonでTukey検定を行うには、StatsModelsというライブラリからstats.multicompモジュールを読み込む必要があります。

>>> from statsmodels.stats.multicomp import pairwise_tukeyhsd
>>> from statsmodels.stats.multicomp import MultiComparison

StatsModelsの公式ドキュメントによれば、関数はstatsmodels.stats.multicomp.pairwise_tukeyhsd(endog, groups, alpha=0.05) の形で呼び出すように書かれています。endogには「比較したいデータそのもの」、groupsは「それらのデータが属する群をインデックスで指定」して、いずれもnumpy.ndarrayの形で渡すようです。

少し不格好ですが、Pythonの通常のリスト型からNumPyのn次元配列に置き換えて関数に渡すコードを実直に書いてみました。

#!/usr/bin/env python
#-*-coding:utf-8-*-

import numpy as np
from statsmodels.stats.multicomp import pairwise_tukeyhsd

A = [100, 96, 80, 92, 101, 99, 89, 93]
B = [83, 77, 86, 71, 70, 98, 80, 75, 74]
C = [90, 91, 92, 89, 88, 90, 91]

def TukeyKramer(*args): #可変長引数で比較するリストを受け取る
    idx = 0
    data_arr = np.array([]) #全てのデータを一次元配列で列記する
    idx_arr = np.array([]) #データに対応するインデックスを一次元配列で列記する
    for lst in args:
        for nums in lst:
            if isinstance(nums, str) == False: #文字列(欠損値など)はskip
                data_arr = np.append(data_arr, nums)
                idx_arr = np.append(idx_arr, idx)
        idx += 1
    print(pairwise_tukeyhsd(data_arr, idx_arr)) #第3引数のdefaultはp=0.05

if __name__ == "__main__":
    TukeyKramer(A, B, C)

上記をシェルから実行すると、出力は以下のようになります。

Multiple Comparison of Means - Tukey HSD,FWER=0.05
==============================================
group1 group2 meandiff  lower    upper  reject
----------------------------------------------
 0.0    1.0   -14.4167 -22.7227 -6.1106  True
 0.0    2.0   -3.6071  -12.454   5.2397 False
 1.0    2.0   10.8095   2.1951   19.424  True
----------------------------------------------

1行目の"0.0(A薬群) vs 1.0(B薬群)" と 3行目の "1.0(B薬群) vs 2.0(C薬群)"で「rejectがTrue(帰無仮説を棄却)」になっていますので、ここにp=0.05未満の有意差があったことが判明しました。
なんとかPythonで多重比較検定を行うことができましたね。
記事が書けてホッとしたのでBiSHの動画でも見ようと思います。

終わりに

SPSSやJMPのような高価なソフトを使わなくても、私たちにはPythonがあります。

...とは言うものの、ネット上の解説記事は非常に少なく、英語オンリーの公式ドキュメントを読み込むのは骨が折れます。私はPythonも統計学も素人ではありますが、似たような境遇の方々にとって、この記事が少しでも支えになれば幸甚です。
引き続き、NumPyやPandasについても勉強を進めていきたいと思います。