Zenn
🧑‍💻

ggplot2でboxplotに多重比較の結果を記入する

2020/12/09に公開
R
ggplot
anova
tech

調べて意外と出てこなかったのでまとめました．

つくりたいもの

こういうboxplotを作りたい．

box上部に書かれたアルファベットは，多重比較検定の結果を表している．

多重比較検定の結果…？

  • abcとつける例のアレ
    例えばある５つの群（Group1~5）があったとき，ANOVAを実行後に各群間の有意差をみると
    Group1とGroup2，Group2と4にそれぞれ有意差がなかった．
    その際は次のように書く．
    Group1 a
    Group2 ab
    Group3 c
    Group4 b
    Group5 d

同じアルファベットが振られた群同士には有意差がない事を示している．
詳しくは以下のサイトを参照してもらいたい．
多重比較の結果を表のみで記す、例のabc…の割り振り方 | 日々是独想 - 日々の徒然なることを独り想う。

abc振りは大変．．．

ルールさえわかってしまえば自分でabc振ってplotに適当に書き込めば良くない？

たしかに上述サイトの通りの作業を行えば手動でも振ることができるが，群数が増えると作業量も増え頭も混乱してくる．

あれ，，えっとAとBは有意，，でそれでCと．．．Eは有意，あ違うかえっと．．．

という作業を多重比較検定の結果を見ながら行わなければいけない．

達成したいこと

そこで今回は

  • abc振りを自動でやってくれて
  • ついでにplotもしてくれる
    コードを書いていく．

Rのmulticompパッケージをつかう

このパッケージを使うと簡単に上の対応を作ってくれる．
検定に使うのはTukeyHSDだ．今回は例としてirisのデータセットを使って実践する．

  • ANOVA実行結果のresをglhtに渡す．
  • さらにcld()に投げて返ってきた結果をmltvに格納しておく．
library(multcomp)
res =aov(Sepal.Length ~ Species, data=iris)
summary(res)
tuk = glht(res,linfct=mcp(Species="Tukey"))
mltv = cld(tuk,decreasing=F)
  • 結果
res =aov(Sepal.Length ~ Species, data=iris)
> summary(res)
             Df Sum Sq Mean Sq F value Pr(>F)    
Species       2  63.21  31.606   119.3 <2e-16 ***
Residuals   147  38.96   0.265                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> tuk = glht(res,linfct=mcp(Species="Tukey"))
> mltv = cld(tuk,decreasing=F)
> mltv
    setosa versicolor  virginica 
       "a"        "b"        "c"

今回は全てに有意差があったためバラバラのアルファベットが振られた つまらないな

一応Tukeyをかけて確認

TukeyHSD(res)

Tukey multiple comparisons of means
   95% family-wise confidence level

Fit: aov(formula = Sepal.Length ~ Species, data = iris)

$Species
                     diff       lwr       upr p adj
versicolor-setosa    0.930 0.6862273 1.1737727     0
virginica-setosa     1.582 1.3382273 1.8257727     0
virginica-versicolor 0.652 0.4082273 0.8957727     0

やっぱもっと面白いデータを使えばよかった．

stat_summary()を使ってboxplotに反映

ここで得られた結果をお手軽にboxplotに乗せるには，ggplotメンバーのstat_summary関数を使うのがいい．
stat_summaryはエラーバーやアノテーションをいい感じにつけてくれる．

テスト

stat_summaryの引数でgeom = 'text'を指定してやると任意の文字列を書き込める．

annos = c('A','B','C')
plot1 = 
  ggplot(data = iris, aes(x = Species, y = Sepal.Length,color = Species))+
  geom_boxplot()+
  stat_summary(geom = 'text', label =annos, fun.y = max, vjust = -1)+
  theme_classic()
print(plot1)


fun.y = maxを指定することで，boxの上部に表示される．便利．

stat_summaryの挙動は確認できたので，ここに比較検定の結果を渡してやる．

RstudioのView()を使って格納箇所を簡単に特定

先程のmltvはデータセットなので，対象のベクトルをいい感じに抽出してやらなければいけない．
ここでRstudioのViewを使うことで簡単にデータの構造を把握できる．

View(mltv)


マーカを引いた箇所に欲しいデータが確認できる．
カーソルを合わせポップアップをクリックするとコード内で使用できる形でアドレスが挿入される．非常に便利だ．

右端に薄っすらとボタンが出ています．

このような形でコピーされる．

mltv[["mcletters"]][["Letters"]]

中身はラベル付きベクトルなのでそのままstat_summaryにわたすことができる．

完成

abc = mltv[["mcletters"]][["Letters"]]
plot2 = 
  ggplot(data = iris, aes(x = Species, y = Sepal.Length,color = Species))+
  geom_boxplot()+
  stat_summary(geom = 'text', label =abc, fun.y = max, vjust = -1)+
  theme_classic()
print(plot2)

以上です．

Discussion