Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

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

追記(2020/2/25)

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

そのため、下のwilcox.exactbrunnermunzel.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の使い方

Masahiro_T
困ったときの再起動
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away