rstatixと仲良くなりたい

さいしょに

rstatixと仲良くなって様々な検定をささっと出来るようになりたいなあ、と思い少しかじった内容をまとめておきます。

rstatix概要

rstatixとはRのライブラリのひとつで、tidyverseのデザインと親和性のあるpipe-friendlyな感じに色々と検定する関数を提供します。
導入は単純にcranからインストールすれば良いですし、以下のようにpacmanを使っても良いです。

#パッケージ管理用のpacmanをインストール
if (!require(pacman)) {install.packages("pacman")} 

#必要なライブラリがインストールされていなければインストールしたうえで呼び出し
pacman::p_load(tidyverse, rstatix)

t検定のやりかた

t検定のやりかたを通してrstatixの特徴を見ていきます。

base関数の場合から説明

最初にbase関数でのt検定を見ていきます。base関数ではt.test()関数を使います。base関数でもパイプ関数を使ってデータを渡すことができますが、したのように.で指定する必要があります。なお、パイプ関数もmagrittrのほうの%>%では可能ですが、nativeな|>のほうでは指定できません。

#base
ToothGrowth %>%
  t.test(len ~ supp, data = .)

で、下ような結果が返ってきます。詳細が書かれているのはありがたいですが、取り扱いが不便です。

	Welch Two Sample t-test

data:  len by supp
t = 1.9153, df = 55.309, p-value = 0.06063
alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
95 percent confidence interval:
 -0.1710156  7.5710156
sample estimates:
mean in group OJ mean in group VC 
        20.66333         16.96333 

rstatixの場合

rstatixではt_test()関数を使います。"."のかわりに"_"になっていますね。dplyrの操作と同じ感じで使えるようになっているので、dataを指定する必要はありません。

#rstatix
ToothGrowth %>%
  t_test(len ~ supp)

下ような結果が返ってきます。tibbleで返ってきますので、とても取り扱いがしやすいです。

# A tibble: 1 × 8
  .y.   group1 group2    n1    n2 statistic    df      p
* <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>  <dbl>
1 len   OJ     VC        30    30      1.92  55.3 0.0606

また、以下のように3水準のものに対しても可能です。base関数のほうでは、2水準ではないのでだめですとエラーが出ますが、t_test()のほうは対応してくれます。

ToothGrowth %>%
  t_test(len ~ dose)

結果は次ように出力されます。デフォルトでは総当たりでt検定を行いp.adjにholmの方法で調整したp値が示されています。ありがたいですね。指定することでbonferroniなど別の方法による調整に変更可能ですが、そのあたりの違いについてはリンク先を参照してください。
https://www.med.osaka-u.ac.jp/pub/kid/clinicaljournalclub1.html

# A tibble: 3 × 10
  .y.   group1 group2    n1    n2 statistic    df        p    p.adj p.adj.signif
* <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>    <dbl> <chr>       
1 len   0.5    1         20    20     -6.48  38.0 1.27e- 7 2.54e- 7 ****        
2 len   0.5    2         20    20    -11.8   36.9 4.40e-14 1.32e-13 ****        
3 len   1      2         20    20     -4.90  37.1 1.91e- 5 1.91e- 5 ****    

まとめ

rstatixを使ってt検定を行う方法についてみてきました。rstatixはpipe-friendlyでささっとコードが書けますし、結果もtibbleで渡されるので非常に取り扱いが楽です。たとえばgroup_by()と組み合わせて検定を行い、p値を小さい順に並べて順番に見ていく、という機械的な取り扱いも簡単にできます。この方法は褒められるようなやり方ではないと思いますが、製造業などで問題がありそうな工程を抽出したいときにひとまず行うには都合がよく、まあまあな正答率を得られたりします（そのときはanovaを主に使いますが、概ね関数の使い方は今回書いたのと同じです）。
以上でrstatixの概要をt検定を通して書きましたが、別の関数でハマったところもあるので、今後はその辺も書いていきたいと思います。

参考