2024/3/18: 同じ機能でコードを書き直しました。
0. 目次
PythonからRのいくつかの統計関数を使えるようにしました。自分でルーティンに使う用。GitHubで公開しています。
前回は環境構築の方法について書きました。
Pythonからよく使うRの統計関数を簡単に呼び出せるようにしてみた (多重比較・Rpy2)
今回は関数の使い方について書いていきます。
- 準備
- TukeyHSD test
- Dunnett's test
- Steel-Dwass test
- Fisher's Exact test
- 計算結果をテキストファイルで保存
多重比較そのものについてはあまり詳しく解説しないので、このあたりのページを参考にしてください。
統計検定を理解せずに使っている人のために III
基本解説→ポストホックテストとしての多重比較検定
1. 準備
rpy2_wrapperディレクトリに移動し、JupyterかPythonを起動します。
Rクラスのインスタンス生成
from rpy2_wrapper import R
#Create R instance
R_inst = R()
csv形式のデータファイルをインポート
読み込みはcsvファイルから行います。
#Import csv file and transform it to both Pandas and R dataframes
R_inst.import_csv("sample_data.csv")
出力結果
column names:
['category', 'data.1', 'data.2', 'data.3', 'data.4', 'data.5']
Import was succeeded
ここで自動的にPandasデータフレームとRデータフレームが生成されます。
Pandasデータフレームは.p_dfに代入されています。
#Pandas dataframe
print(type(R_inst.p_df))
R_inst.p_df
出力結果
<class 'pandas.core.frame.DataFrame'>
| category | data 1 | data 2 | data3 | data 4 | data 5 | |
|---|---|---|---|---|---|---|
| 0 | A | 0 | 3 | 5 | 8 | 12 |
| 1 | A | 1 | 4 | 6 | 10 | 14 |
| 2 | A | 2 | 6 | 7 | 11 | 15 |
| 3 | A | 0 | 2 | 4 | 8 | 13 |
| 4 | B | 2 | 3 | 4 | 4 | 5 |
| 5 | B | 0 | 1 | 2 | 2 | 3 |
| 6 | B | 1 | 2 | 2 | 3 | 4 |
| 7 | B | 1 | 1 | 2 | 3 | 5 |
| 8 | B | 2 | 2 | 3 | 4 | 5 |
| 9 | C | 2 | 3 | 4 | 5 | 6 |
| 10 | C | 2 | 4 | 5 | 6 | 7 |
| 11 | C | 1 | 2 | 4 | 5 | 6 |
| 12 | C | 0 | 3 | 3 | 4 | 5 |
一方、Rデータフレームは.r_dfに代入されています。
#R dataframe
print(type(R_inst.r_df))
出力結果
<class 'rpy2.robjects.vectors.DataFrame'>
表はMarkdownでうまく表示できなかったので省略。
2. TukeyHSD test
TukeyHSD testを実施します。多重比較で通常、One-way ANOVAの後にやるPost-hoc testの１つです。ざっくり言うと、全てのカテゴリーに対して1対1の比較を行うための方法です。内容については以下のページが詳しいです。
チューキー・クレーマー検定は，チューキー検定 (Tukey test) とかチューキーの範囲検定 (Tukey range test) とかチューキーのHSD検定 (Tukey honestly significant difference test) 等とも呼ばれる多重比較の検定法のひとつであり，多群のデータ中における各2群間の平均値の差についての検定を行う方法である．
rpy2_wrapperでは、ANOVAとTukeyHSDが同時に実行・出力されるように実装しているので、事前にANOVAを実施する必要は必ずしもありません。
実行する際、ラベルとデータのカラム名をそれぞれ引数で指定します。ここではラベルは"category"、データは"data.5"に入っています。カラム名はimport_csvを実行時に出力されるのでそちらを参照してください(PandasからRに変換するときに文字化けすることがあるので注意)。
#Perform TukeyHSD
R_inst.tukeyHSD(data_column="data.5",label_column="category")
実行結果は以下のように、Rで実行したときと同じ形式でJupyter notebookに出力されます。
-------------Result of ANOVA-------------
Analysis of Variance Table
Response: data.5
Df Sum Sq Mean Sq F value Pr(>F)
label_data 3 242.366 80.789 85.096 7.973e-11 ***
Residuals 18 17.089 0.949
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-------------Result of TukeyHSD test-------------
diff lwr upr p adj
B-A -9.1000000 -10.9473264 -7.2526736 2.573011e-10
C-A -7.5000000 -9.4472530 -5.5527470 1.350361e-08
D-A -8.6111111 -10.2659580 -6.9562643 1.042004e-10
C-B 1.6000000 -0.2473264 3.4473264 1.033965e-01
D-B 0.4888889 -1.0471250 2.0249028 8.051780e-01
D-C -1.1111111 -2.7659580 0.5437357 2.637763e-01
p adjがP値。
3. Dunnett's test
TukeyHSDは1対1の検定を全てのカテゴリーに対して行いましたが、Dunnett's testではカテゴリーのうちの1つを対照区として、それと他のカテゴリーとの1対1の比較を行います。事前にANOVAはやらなくても良いらしいですが、このコードを使うと実行時に一緒にANOVAの結果も出ちゃいます。
R にてダネット検定 (Dunnett test) を行う．ダネット検定は，ひとつの対照群 (コントロール群) と2つ以上の処理群からなる多群のデータ中における対照群と処理群の各2群間の平均値の差についての検定を行う方法である．
対照区は指定できないので、カテゴリー名を工夫する必要がありそうです。
同ページ内からの引用を以下に載せておきます。
水準名 (ここでの C，X，Y，Z) をアルファベット順 (水準名が数字の場合，降順) で並べたときの最も若いものを対照群にするので，そこに注意しなければ正確な対比較が行えない．
早速実行してみます。
#Dunnett's test
rinst.dunnett_test(data_column=data_column,label_column=label_column)
-------------Result of ANOVA-------------
Analysis of Variance Table
Response: data.5
Df Sum Sq Mean Sq F value Pr(>F)
label_data 3 242.366 80.789 85.096 7.973e-11 ***
Residuals 18 17.089 0.949
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-------------Result of Dunnett's test-------------
Estimate Std. Error t value Pr(>|t|)
B - A -9.100000 0.653622 -13.922412 0.000000e+00
C - A -7.500000 0.688978 -10.885681 2.231247e-10
D - A -8.611111 0.585519 -14.706797 9.992007e-16
Pr(>|t|)がP値。
地味に、この関数の実装が一番大変でした。
4. Steel-Dwass test
最後にSteel-Dwass test。TukeyHSD testのノンパラメトリック版だそうです。正直勉強不足で理解しきれていないですが、とりあえず実装をば。
方法として、3つ設定可能ですが、デフォルトの設定で実装しています。順列が少ないときはExact、多い時はMonte Carloが自動的に選択されるそうです。詳しくは以下を参照のこと。
pSDCFlig: Dwass, Steel, Critchlow, Fligner
Either "Exact", "Monte Carlo", or "Asymptotic", indicating the desired distribution. When method=NA, "Exact" will be used if the number of permutations is 10,000 or less. Otherwise, "Monte Carlo" will be used.
以下のページに説明があるのですが、モンテカルロについてはあまり詳しく紹介されていないのでよくわかりません。時間があるときに突っ込んで勉強したいところ。
多重比較 Steel-Dwass 正規近似と正確検定： U 検定が基準
データが多いと固まるので、少なめのデータを読み込み直して実行します。
import_data = "sample_data3.csv"
rinst3 = R()
rinst3.import_csv(import_data)
rinst3.p_df.head()
出力結果
column names:
['category', 'data.1', 'data.2', 'data.3', 'data.4', 'data.5']
Import was succeeded
| category | data 1 | data 2 | data 3 | data 4 | data 5 | |
|---|---|---|---|---|---|---|
| 0 | A | 0 | 3 | 5 | 8 | 12 |
| 1 | A | 1 | 4 | 6 | 10 | 14 |
| 2 | A | 2 | 6 | 7 | 11 | 15 |
| 3 | A | 0 | 2 | 4 | 8 | 13 |
| 4 | B | 2 | 3 | 4 | 4 | 5 |
実行してみます。
data_column="data.5"
label_column="category"
#Steel-Dwass test with Monte Carlo
rinst3.steeldwass(data_column=data_column,label_column=label_column)
rinst3.close()
-------------Result of Steel-Dwass test-------------
Ties are present, so p-values are based on conditional null distribution.
Group sizes: 4 5 4
Using the Monte Carlo (with 10000 Iterations) method:
For treatments A - B, the Dwass, Steel, Critchlow-Fligner W Statistic is -3.5233.
The smallest experimentwise error rate leading to rejection is 0.012 .
For treatments A - C, the Dwass, Steel, Critchlow-Fligner W Statistic is -3.2856.
The smallest experimentwise error rate leading to rejection is 0.0444 .
For treatments B - C, the Dwass, Steel, Critchlow-Fligner W Statistic is 3.0895.
The smallest experimentwise error rate leading to rejection is 0.0676 .
ちょっと時間がかかります。
5. Fisher's exact test
最後に、Fisher's exact testです。クロス表で独立性の検定に使う手法です。
直接確率検定は正確確率検定と訳されることも多い．...<中略>... 本検定法は分割表 (クロス集計表) における2つ以上のカテゴリーが独立であるかどうかを調べる手法である． ...<中略>... 同様の検定を行うものにカイ二乗検定 (適合度検定・独立性の検定) があるが，分割表のセルの期待値に10未満のものがある場合はカイ二乗検定は避け，フィッシャーの直接確率検定を用いるほうが良いとされる．
こちらはクロス表を読み込んでそのまま実行するようにしたので、Import_csvするときにカテゴリーのカラムをインデックスに指定してください。
import_data = "sample_data2.csv"
rinst2 = R()
rinst2.import_csv(import_data,index_col=0)
rinst2.p_df
column names:
['A', 'B', 'C', 'D']
Import was succeeded
| A | B | C | D | |
|---|---|---|---|---|
| positive | 2 | 8 | 9 | 4 |
| negative | 12 | 5 | 2 | 10 |
では、実行してみます。引数はなしでOK。
#Fisher's Exact Test
rinst2.fisher_test()
出力結果
-------------Result of Fisher's Exact test-------------
Fisher's Exact Test for Count Data
data:
p-value = 0.002194
alternative hypothesis: two.sided
以上です。Rと同じように出力できています。
6. 計算結果をテキストファイルで保存
テキストファイルとして保存できるようにしました。
各コードの引数にsave_file="{name}.txt"を突っ込んでもらえれば、指定のファイルの最後に追記されます。
次回は実装したコードを見ながら、Rpy2自体の使い方を軽く紹介する予定です。
一般的なRpy2の使い方について、記事にしました。(2018/7/29)
Rpy2の使い方の備忘録