rpy2 を使って R の wilcox.exact を使う（scipy の マンホイットニーU検定 の実装が 少サンプル数非対応 なので、代わりに）

追記（2020/2/25）

マンホイットニーのU検定よりBrunner-Munzel検定の方が良いそうです。

そのため、下のwilcox.exactをbrunnermunzel.test（片群のサンプルサイズが10程度以下ならbrunnermunzel.permutation.test）に置き換えて利用すると良いはずです。

参考

奥村晴彦 先生の解説：
Brunner-Munzel検定

他の検定（t検定、Welchのt検定、U検定）との比較：
マイナーだけど最強の統計的検定 Brunner-Munzel 検定 - ほくそ笑む

経緯

Jupyter notebook でマンホイットニーのU検定をしようと思いました。U検定を実装しているライブラリを探してみた所、scipyで実装されていました。具体的にはscipy.stats.mannwhitneyuという名前で実装されています。

このライブラリを利用しようと思い、引数等を確認するために公式Docを見た所、Notesに次の様な注記がありました。

Use only when the number of observation in each sample is > 20

訳すと、「2群それぞれのサンプルサイズが 20超 の場合のみ利用可能」となります。

なぜこのような条件が付されているのか調べた所、どうやらこのライブラリは「標準正規分布による近似をして検定する方法」のみを実装している事が分かりました。

この事はソースコードとWikipediaを見て分かりました。

そして、scipy の Github を調べると、どうやら片群のサンプルサイズが20以下の場合の実装はプルリク中のようでした。

さて、私は片群のサンプルサイズが20以下の場合の検定を行いたかったので、別のライブラリを探す必要が生まれました。

そこで、Python ライブラリをざっと調べて見たのですが、良さそうなものがパッとは見つかりませんでした。

どうしようか考えた所、「Rの実装を使えばいいのでは」と思いました。

そこで、RでU検定を実装したライブラリを探した所、wilcox.exactというライブラリが見つかりました。

このライブラリは、データ科学便覧で紹介されています。

データ科学便覧の解説によると、サンプルサイズによって自動で「標準正規分布による近似」をするか「p値を正確に計算」するかを決定してくれる事が分かりました。そしてこの機能により、scipyの様なサンプルサイズの制約はない事が分かりました。

この様に、良さそうなライブラリがRで見つかったので、Pythonカーネル上でRを使えるようにするライブラリrpy2を使って、wilcox.exactを使う事にしました。

方法

1. rpy2を install します

terminal

pip install rpy2

2. exactRankTests を install します

terminal

# R を起動します
$ R

R（teminal上）

# exactRankTests を install
> install.packages("exactRankTests", repos="http://cran.ism.ac.jp/")

# 終了
> q()

3. Jupyter notebook 上で次のコードを実行して動作確認をします

Input（Jupyter　notebook）

# %R が使えるようにする
%load_ext rpy2.ipython

# U検定が入っているライブラリをロード
%R library(exactRankTests)

# Python の配列を引数にとる U検定メソッド を作成
def u_test(x, y):
    
    # ndarray に変換
    # ndarrayなら 「%R -i」 した時に R で 「c()」 した時の形と同じになる
    x = np.array(x)
    y = np.array(y)

    # R側に x と y を読み込む
    %R -i x
    %R -i y
    
    # U検定
    %R output = wilcox.exact(x=x,y=y,paired=F)
    
    # Python 側に出力
    %R -o output
    
    # ListVector型で出力されるので、dictionary 型に変換
    result = dict(zip(output.names, map(list,list(output))))
    
    return result

# 試してみる
x = np.array([0.15165101, 0.47173771, 0.46735185, 0.67045667, 0.17489188])
y = x * 10
u_test(x, y)

Output（Jupyter　notebook）

{'statistic': [0.0],
 'pointprob': [0.003968253968253968],
 'p.value': [0.007936507936507936],
 'null.value': [0.0],
 'alternative': ['two.sided'],
 'method': ['Exact Wilcoxon rank sum test'],
 'data.name': ['x and y']}

p.valueがp値です。

参考

rpy2の使い方

• Jupyterで1つのノートブック内にRとPythonを共存させる - Qiita