rpy2_wrapperを使ってPythonコードで統計検定（TukeyHSD・Dunnett・Steel-Dwassなど)

2024/3/18: 同じ機能でコードを書き直しました。

https://qiita.com/MTNakata/items/657d71abaab5ad29515f

---

# 0. 目次

PythonからRのいくつかの統計関数を使えるようにしました。自分でルーティンに使う用。[GitHubで公開しています](https://github.com/MTNakata/rpy2_wrapper)。

前回は環境構築の方法について書きました。
[Pythonからよく使うRの統計関数を簡単に呼び出せるようにしてみた (多重比較・Rpy2)](https://qiita.com/MTNakata/items/129d334cea397a6ec0c3)

今回は関数の使い方について書いていきます。


1. 準備
2. TukeyHSD test
3. Dunnett's test
4. Steel-Dwass test
5. Fisher's Exact test
6. 計算結果をテキストファイルで保存

多重比較そのものについてはあまり詳しく解説しないので、このあたりのページを参考にしてください。

[統計検定を理解せずに使っている人のために III](https://www.jstage.jst.go.jp/article/kagakutoseibutsu/51/7/51_483/_pdf)
[基本解説→ポストホックテストとしての多重比較検定](http://www.shiga-med.ac.jp/~koyama/stat/com-ph.html)

# 1. 準備

`rpy2_wrapper`ディレクトリに移動し、JupyterかPythonを起動します。

### Rクラスのインスタンス生成

```python
from rpy2_wrapper import R
#Create R instance
R_inst = R()
```
### csv形式のデータファイルをインポート

読み込みはcsvファイルから行います。

```python
#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`に代入されています。

```python
#Pandas dataframe
print(type(R_inst.p_df))
R_inst.p_df
```

出力結果

```
    <class 'pandas.core.frame.DataFrame'>
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>category</th>
      <th>data 1</th>
      <th>data 2</th>
      <th>data3</th>
      <th>data 4</th>
      <th>data 5</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>A</td>
      <td>0</td>
      <td>3</td>
      <td>5</td>
      <td>8</td>
      <td>12</td>
    </tr>
    <tr>
      <th>1</th>
      <td>A</td>
      <td>1</td>
      <td>4</td>
      <td>6</td>
      <td>10</td>
      <td>14</td>
    </tr>
    <tr>
      <th>2</th>
      <td>A</td>
      <td>2</td>
      <td>6</td>
      <td>7</td>
      <td>11</td>
      <td>15</td>
    </tr>
    <tr>
      <th>3</th>
      <td>A</td>
      <td>0</td>
      <td>2</td>
      <td>4</td>
      <td>8</td>
      <td>13</td>
    </tr>
    <tr>
      <th>4</th>
      <td>B</td>
      <td>2</td>
      <td>3</td>
      <td>4</td>
      <td>4</td>
      <td>5</td>
    </tr>
    <tr>
      <th>5</th>
      <td>B</td>
      <td>0</td>
      <td>1</td>
      <td>2</td>
      <td>2</td>
      <td>3</td>
    </tr>
    <tr>
      <th>6</th>
      <td>B</td>
      <td>1</td>
      <td>2</td>
      <td>2</td>
      <td>3</td>
      <td>4</td>
    </tr>
    <tr>
      <th>7</th>
      <td>B</td>
      <td>1</td>
      <td>1</td>
      <td>2</td>
      <td>3</td>
      <td>5</td>
    </tr>
    <tr>
      <th>8</th>
      <td>B</td>
      <td>2</td>
      <td>2</td>
      <td>3</td>
      <td>4</td>
      <td>5</td>
    </tr>
    <tr>
      <th>9</th>
      <td>C</td>
      <td>2</td>
      <td>3</td>
      <td>4</td>
      <td>5</td>
      <td>6</td>
    </tr>
    <tr>
      <th>10</th>
      <td>C</td>
      <td>2</td>
      <td>4</td>
      <td>5</td>
      <td>6</td>
      <td>7</td>
    </tr>
    <tr>
      <th>11</th>
      <td>C</td>
      <td>1</td>
      <td>2</td>
      <td>4</td>
      <td>5</td>
      <td>6</td>
    </tr>
    <tr>
      <th>12</th>
      <td>C</td>
      <td>0</td>
      <td>3</td>
      <td>3</td>
      <td>4</td>
      <td>5</td>
    </tr>
  </tbody>
</table>
</div>


一方、Rデータフレームは`.r_df`に代入されています。

```python
#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の比較を行うための方法です。内容については以下のページが詳しいです。

[Rによるチューキー・クレーマー検定](https://data-science.gr.jp/implementation/ist_r_tukey_kramer_test.html)

>チューキー・クレーマー検定は，チューキー検定 (Tukey test) とかチューキーの範囲検定 (Tukey range test) とかチューキーのHSD検定 (Tukey honestly significant difference test) 等とも呼ばれる多重比較の検定法のひとつであり，多群のデータ中における各2群間の平均値の差についての検定を行う方法である．

rpy2_wrapperでは、ANOVAとTukeyHSDが同時に実行・出力されるように実装しているので、事前にANOVAを実施する必要は必ずしもありません。

実行する際、ラベルとデータのカラム名をそれぞれ引数で指定します。ここではラベルは"category"、データは"data.5"に入っています。カラム名は`import_csv`を実行時に出力されるのでそちらを参照してください(PandasからRに変換するときに文字化けすることがあるので注意)。

```python
#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によるダネット検定](https://data-science.gr.jp/implementation/ist_r_dunnett_test.html)

>R にてダネット検定 (Dunnett test) を行う．ダネット検定は，ひとつの対照群 (コントロール群) と2つ以上の処理群からなる多群のデータ中における対照群と処理群の各2群間の平均値の差についての検定を行う方法である．

対照区は指定できないので、カテゴリー名を工夫する必要がありそうです。

同ページ内からの引用を以下に載せておきます。

>水準名 (ここでの C，X，Y，Z) をアルファベット順 (水準名が数字の場合，降順) で並べたときの最も若いものを対照群にするので，そこに注意しなければ正確な対比較が行えない．

早速実行してみます。

```python
#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](https://rdrr.io/cran/NSM3/man/pSDCFlig.html)

>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 検定が基準](https://sites.google.com/site/fishermultiplecomparison/steel-dwass)

データが多いと固まるので、少なめのデータを読み込み直して実行します。

```python
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
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>category</th>
      <th>data 1</th>
      <th>data 2</th>
      <th>data 3</th>
      <th>data 4</th>
      <th>data 5</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>A</td>
      <td>0</td>
      <td>3</td>
      <td>5</td>
      <td>8</td>
      <td>12</td>
    </tr>
    <tr>
      <th>1</th>
      <td>A</td>
      <td>1</td>
      <td>4</td>
      <td>6</td>
      <td>10</td>
      <td>14</td>
    </tr>
    <tr>
      <th>2</th>
      <td>A</td>
      <td>2</td>
      <td>6</td>
      <td>7</td>
      <td>11</td>
      <td>15</td>
    </tr>
    <tr>
      <th>3</th>
      <td>A</td>
      <td>0</td>
      <td>2</td>
      <td>4</td>
      <td>8</td>
      <td>13</td>
    </tr>
    <tr>
      <th>4</th>
      <td>B</td>
      <td>2</td>
      <td>3</td>
      <td>4</td>
      <td>4</td>
      <td>5</td>
    </tr>
  </tbody>
</table>
</div>


実行してみます。

```python
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です。クロス表で独立性の検定に使う手法です。

[Rによるフィッシャーの直接確率検定](https://data-science.gr.jp/implementation/ist_r_fisher_exact_probability_test.html)

>直接確率検定は正確確率検定と訳されることも多い．...<中略>... 本検定法は分割表 (クロス集計表) における2つ以上のカテゴリーが独立であるかどうかを調べる手法である． ...<中略>... 同様の検定を行うものにカイ二乗検定 (適合度検定・独立性の検定) があるが，分割表のセルの期待値に10未満のものがある場合はカイ二乗検定は避け，フィッシャーの直接確率検定を用いるほうが良いとされる．

こちらはクロス表を読み込んでそのまま実行するようにしたので、`Import_csv`するときにカテゴリーのカラムをインデックスに指定してください。

```python
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
```



<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>A</th>
      <th>B</th>
      <th>C</th>
      <th>D</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>positive</th>
      <td>2</td>
      <td>8</td>
      <td>9</td>
      <td>4</td>
    </tr>
    <tr>
      <th>negative</th>
      <td>12</td>
      <td>5</td>
      <td>2</td>
      <td>10</td>
    </tr>
  </tbody>
</table>
</div>



では、実行してみます。引数はなしでOK。

```python
#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の使い方の備忘録](https://qiita.com/MTNakata/items/04a045f4beb0fa9e3a4d)