Rのggplot2でグラフを描いたときに有意差の記号、*とかそういうのをグラフ上に書こうとしたときの方法について、私的まとめ。 いつもすぐに忘れてググってる。
グラフにannotation描く動機がなければ興味のない話になっている。自分はこうやってるよというのを晒すことでより良い知見にたどり着きたい。
あと必要なときにどこにやったのかわからないRのスクリプト探すより自分のブログをコピペしたほうが早いというのもある。
- ggplotで描いたグラフに有意差のマークを書く方法
- ggplotのgeom_textとgeom_segmentを使って有意差のマークをグラフ上に描く
- ggsignif::geom_signifを使って有意差のマークをグラフ上に描く
- ggpubr::stat_compare_meansを使って有意差のマークをグラフ上に描く
- 雑感
ggplotで描いたグラフに有意差のマークを書く方法
ggplotで書いたグラフに有意差のマークを書く方法は自分の把握している範囲で3つある。
- ggplotで描いたグラフを画像ファイルとして出力したあとに、画像処理用のソフトをつかって画像ファイルに書き込む。
- ggplot上でggplotの関数の
geom_text
とgeom_segment
をつかって書き込む。 ggsignif::geom_signif
やggpubr::stat_compare_means
などの別パッケージを使ってggplotで描いた(あるいはggpubr
で描いた)グラフに書き込む。
覚えることが少ないのはggplotで描いたグラフを画像ファイルとして出力したあとに、別のソフトをつかって有意差のマークを書き込む方法だが、解像度が微妙に合っていなくて、グラフや有意差のマークのどっちかがボケたり、グラフと線の位置をピッタリ合わせようとすると意外とストレスになる。本稿では触れない。
本稿では、残りの2つの方法について簡素にまとめていく。
個人的な所感としては、データの形がハマればggpubr::stat_compare_means
が楽だが、データがはまらない場合はggsignif::geom_signif
が使いやすいかと思う。データがそれぞれのパッケージでうまく処理できない形だと、位置指定をしなければいけないのでちょっと面倒。
ちなみに、ggpubr::stat_compare_means
は内部でggsignif::geom_signif
を使っている。
ggplot2
の用意しているgeom_text
やgeom_segment
を使って描く方法は、要するにグラフ上にせっせと領域指定して罫線ひいているだけなのでだいたい場合で万能。ただし、位置合わせが面倒だったり、コードが長くなる。
ggplotのgeom_text
とgeom_segment
を使って有意差のマークをグラフ上に描く
geom_text
はグラフ上の任意の箇所に文字(記号)を描く関数で、geom_segment
は始点と終点を指定して罫線を引くggplot2
で用意されている関数。
これらをつかってグラフ上にお絵かきしていく。
まずはRの組み込みデータのIRISデータを使って種別のがく片の長さの箱ひげ図を描く
# ggplot2インストールしていなければinstall.packages("ggplot2")でインストールする。 library(ggplot2) ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot()
これは見るからに差があるが、一元配置分散分析で確認すると
summary(aov(Sepal.Length~Species, data = iris))
## 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
TukeyHSD(aov(Sepal.Length~Species, data = iris))
## 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
とテューキーの多重比較でもすべての群間で有意な差があるということで、ざっとp\<0.001ということで***を各群間に書き込みたくなる。
はじめにsetosaとversicolor間で描いてみる。
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # setosaとverisicolorの2群の間(x = 1.5, y = 7.5)に有意差のマークをつける geom_text(x = 1.5, y = 7.7, label = "***") + # setosa(x = 1の座標)のグラフ上ちょっと上(y = 7.5から7.2に)から縦線をちょっと引いて geom_segment(x = 1, xend = 1, y = 7.5, yend = 7.2) + # setosa(x = 1)からversicolor(x = 2)にさっき引いた直線(y = 7.5)とつながるように線を引いて geom_segment(x = 1, xend = 2, y = 7.5, yend = 7.5) + # Versicolor(x = 2)でグラフの上にちょっと直線をsetosa上に描いたのと同じ位置(y = 7.5から7.2)に描く geom_segment(x = 2, xend = 2, y = 7.5, yend = 7.2)
縦線つけないで横線と米印だけなら2行で済むが、縦線もつけてコを横倒しした形にするとなると4行はコードを書く必要がある。
3群すべてで描くと同様の手順(ただしこれだけではだめなコードを示す)で、
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # setosaとverisicolor間 geom_text(x = 1.5, y = 7.7, label = "***") + geom_segment(x = 1, xend = 1, y = 7.5, yend = 7.2) + geom_segment(x = 1, xend = 2, y = 7.5, yend = 7.5) + geom_segment(x = 2, xend = 2, y = 7.5, yend = 7.2) + # vesicolorとvirginica間 geom_text(x = 2.5, y = 8.5, label = "***") + geom_segment(x = 2, xend = 2, y = 8, yend = 8.3) + geom_segment(x = 2, xend = 3, y = 8.3, yend = 8.3) + geom_segment(x = 3, xend = 3, y = 8, yend = 8.3) + # setosaとvirginica間 geom_text(x = 2, y = 9, label = "***") + geom_segment(x = 1, xend = 1, y = 8.4, yend = 8.7) + geom_segment(x = 1, xend = 3, y = 8.7, yend = 8.7) + geom_segment(x = 3, xend = 3, y = 8.4, yend = 8.7)
となり、ちゃんと描けない。右上のあたりに線が見切れているのでわかるかもしれないが、描画範囲内に有意差のマークが収まっていないから。
ということで、ylim
で表示範囲を調整する。
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # setosaとverisicolor間 geom_text(x = 1.5, y = 7.7, label = "***") + geom_segment(x = 1, xend = 1, y = 7.5, yend = 7.2) + geom_segment(x = 1, xend = 2, y = 7.5, yend = 7.5) + geom_segment(x = 2, xend = 2, y = 7.5, yend = 7.2) + # vesicolorとvirginica間 geom_text(x = 2.5, y = 8.5, label = "***") + geom_segment(x = 2, xend = 2, y = 8, yend = 8.3) + geom_segment(x = 2, xend = 3, y = 8.3, yend = 8.3) + geom_segment(x = 3, xend = 3, y = 8, yend = 8.3) + # setosaとvirginica間 geom_text(x = 2, y = 9, label = "***") + geom_segment(x = 1, xend = 1, y = 8.4, yend = 8.7) + geom_segment(x = 1, xend = 3, y = 8.7, yend = 8.7) + geom_segment(x = 3, xend = 3, y = 8.4, yend = 8.7) + # yの範囲を有意差のマークがすべて入るように設定 ylim(c(4.25, 9.1))
でなんとなくいい感じなる。あとは図のサイズなどもろもいじるといい感じになる。
この方法の強みとしては、グラフ上で線や米印をつけたい箇所の座標がわかれば簡単にできる。書きたいところに書きたい長さで書きたいものが簡単にかける。 あとはオプションで色や線の種類も簡単に変えられる。
弱みは入力するものがちょっと多いこと、後述するお絵かき座標用のデータを準備するなどで工夫もできなくない。
あとは、グラフ上の座標を線一本ずつ、米印の位置を入力したり調整したりしないといけない。後述する有意差のマークを付けるためのパッケージをつかった関数を使うと、うまくハマれば位置も簡単に変数名で指定できるし、p値を勝手に計算してくれる。 自分でお絵描きする方法だと計算しなければいけないというのが一手間だが、関数が計算してくれるのも個人的に一長一短があると思っている(この点は最後に書く)。
あとは、上記の書き方だともともと描画される範囲外にマークを書こうとするとylim
で範囲指定しなければいけないけど、少しあとで書くaes()
で描く方法だと範囲指定がいらない。
今回のIRISデータみたいに位置がわかりやすいと簡単だが、ちょっと複雑だとやりにくい。例えばRの組み込みデータセットのToothGrowthで、ビタミンCの投与量の箱ひげ図を投与方法で層化したグラフを描くと、
ggplot(ToothGrowth, aes(x = factor(dose), y = len, fill = supp)) + geom_boxplot()
となる。
着目する視点次第だけれど今回は微調整が面倒な例として上げたいので、同一投与量内での投与方法での効果の差に注目したときに、
with(ToothGrowth, t.test(len[dose == 0.5 & supp == "OJ"], len[dose == 0.5 & supp == "VC"]))
##
## Welch Two Sample t-test
##
## data: len[dose == 0.5 & supp == "OJ"] and len[dose == 0.5 & supp == "VC"]
## t = 3.1697, df = 14.969, p-value = 0.006359
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.719057 8.780943
## sample estimates:
## mean of x mean of y
## 13.23 7.98
with(ToothGrowth, t.test(len[dose == 1 & supp == "OJ"], len[dose == 1 & supp == "VC"]))
##
## Welch Two Sample t-test
##
## data: len[dose == 1 & supp == "OJ"] and len[dose == 1 & supp == "VC"]
## t = 4.0328, df = 15.358, p-value = 0.001038
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2.802148 9.057852
## sample estimates:
## mean of x mean of y
## 22.70 16.77
with(ToothGrowth, t.test(len[dose == 2 & supp == "OJ"], len[dose == 2 & supp == "VC"]))
##
## Welch Two Sample t-test
##
## data: len[dose == 2 & supp == "OJ"] and len[dose == 2 & supp == "VC"]
## t = -0.046136, df = 14.04, p-value = 0.9639
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.79807 3.63807
## sample estimates:
## mean of x mean of y
## 26.06 26.14
となるので、doseが0.5と1では有意な差があるが、2では有意な差がないということでIRIS同様にコを横倒しにしたのに米印を書きたいと思う。
このときに、dose0.5のOJ側のx座標は微妙な位置にあるということ。書きながら微調整していくしか無い。
例としては
ggplot(ToothGrowth, aes(x = factor(dose), y = len, fill = supp)) + geom_boxplot() + # dose0.5 geom_text(x = 1, y = 25, label = "**") + geom_segment(x = 0.81, xend = 0.81, y = 22, yend = 24) + geom_segment(x = 0.81, xend = 1.19, y = 24, yend = 24) + geom_segment(x = 1.19, xend = 1.19, y = 22, yend = 24) + # dose1 geom_text(x = 2, y = 31, label = "**") + geom_segment(x = 1.81, xend = 1.81, y = 28, yend = 30) + geom_segment(x = 1.81, xend = 2.19, y = 30, yend = 30) + geom_segment(x = 2.19, xend = 2.19, y = 28, yend = 30)
となる、x座標が微妙なところになっているので微調整しながら妥協点を見つけていく。一箇所見つければあとは対称な位置にあったりするので簡単だけど。
あとはこれの色やテーマや諸々を調整するといい感じの図になる。
ちなみに、ggplotで図を書くときにggplot()
のところにデータやエステティックを書く派と、グラフの形を指定したときに書く派がいると思うが、上記までの書き方だと描けない(次に解決策を書く)。
ggplot() + geom_boxplot(aes(x = Species, y = Sepal.Length), iris) + geom_text(x = 1.5, y = 7.7, label = "***") + geom_segment(x = 1, xend = 1, y = 7.5, yend = 7.2) + geom_segment(x = 1, xend = 2, y = 7.5, yend = 7.5) + geom_segment(x = 2, xend = 2, y = 7.5, yend = 7.2)
geom_text
とgeom_segment
のそれぞれでエステティックを書く必要がある。aes()
の中身はaes
書かなかった上のときと同じでOK。
ggplot() + geom_boxplot(aes(x = Species, y = Sepal.Length), iris) + geom_text(aes(x = 1.5, y = 7.7, label = "***")) + geom_segment(aes(x = 1, xend = 1, y = 7.5, yend = 7.2)) + geom_segment(aes(x = 1, xend = 2, y = 7.5, yend = 7.5)) + geom_segment(aes(x = 2, xend = 2, y = 7.5, yend = 7.2))
この方法だと、範囲指定しなくてもエステティックの要素でグラフの描画範囲をいい塩梅に調整してくれるので、
ggplot() + geom_boxplot(aes(x = Species, y = Sepal.Length), iris) + # setosaとverisicolor間 geom_text(aes(x = 1.5, y = 7.7, label = "***")) + geom_segment(aes(x = 1, xend = 1, y = 7.5, yend = 7.2)) + geom_segment(aes(x = 1, xend = 2, y = 7.5, yend = 7.5)) + geom_segment(aes(x = 2, xend = 2, y = 7.5, yend = 7.2)) + # vesicolorとvirginica間 geom_text(aes(x = 2.5, y = 8.5, label = "***")) + geom_segment(aes(x = 2, xend = 2, y = 8, yend = 8.3)) + geom_segment(aes(x = 2, xend = 3, y = 8.3, yend = 8.3)) + geom_segment(aes(x = 3, xend = 3, y = 8, yend = 8.3)) + # setosaとvirginica間 geom_text(aes(x = 2, y = 9, label = "***")) + geom_segment(aes(x = 1, xend = 1, y = 8.4, yend = 8.7)) + geom_segment(aes(x = 1, xend = 3, y = 8.7, yend = 8.7)) + geom_segment(aes(x = 3, xend = 3, y = 8.4, yend = 8.7))
となる。
ただし、毎回aes()
と書くのがちょっと手間。
ちなみにエステティックで描いていればグラフの描画範囲が調整されるので、ggplot()
内にデータやエステティックを書く派も、
ggplot(iris, aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # setosaとverisicolor間 geom_text(aes(x = 1.5, y = 7.7, label = "***")) + geom_segment(aes(x = 1, xend = 1, y = 7.5, yend = 7.2)) + geom_segment(aes(x = 1, xend = 2, y = 7.5, yend = 7.5)) + geom_segment(aes(x = 2, xend = 2, y = 7.5, yend = 7.2)) + # vesicolorとvirginica間 geom_text(aes(x = 2.5, y = 8.5, label = "***")) + geom_segment(aes(x = 2, xend = 2, y = 8, yend = 8.3)) + geom_segment(aes(x = 2, xend = 3, y = 8.3, yend = 8.3)) + geom_segment(aes(x = 3, xend = 3, y = 8, yend = 8.3)) + # setosaとvirginica間 geom_text(aes(x = 2, y = 9, label = "***")) + geom_segment(aes(x = 1, xend = 1, y = 8.4, yend = 8.7)) + geom_segment(aes(x = 1, xend = 3, y = 8.7, yend = 8.7)) + geom_segment(aes(x = 3, xend = 3, y = 8.4, yend = 8.7))
でylim()
を使わずとも描画範囲を調子してくれる。しかし、打鍵量はaes()を4回を群の数だけ多く書くのでylim()
で範囲指定したほうが書く分量は少ない。
エステティックで書く方法の旨味は範囲を自動調整してくれることよりも、座標をデータフレームで与えて描いてもらえることにあると思う。いままで、お絵描きするための座標をそれぞれでベタ打ちしていたが、座標用のデータを作ってそれを参照して作ることもできる。これならば、思った図にならなかっときに直す箇所がベタ打ちしているときよりも少なく、かつ直し忘れが少なくなるかもしれない。まずは、setosaとversicolorの2群間のみに印をつけてみる。
d.setosa <- data.frame(x = 1, xend = 2, y = 7.2, yend = 7.5, label = "***", label.x = 1.5, label.y = 7.7) ggplot(iris, aes(x = Species, y = Sepal.Length)) + geom_boxplot() + geom_text(aes(x = label.x, y = label.y, label = label), d.setosa) + geom_segment(aes(x = x, xend = x, y = y , yend = yend), d.setosa) + geom_segment(aes(x = x, xend = xend, y = yend, yend = yend), d.setosa) + geom_segment(aes(x = xend, xend = xend, y = y, yend = yend), d.setosa)
この方式なら、3群間にしてもgeom_text
やgeom_segment
の入力回数も減らせるし、コードもスッキリする。
d.signif <- data.frame(x = c(1, 2, 1), xend = c(2, 3, 3), y = c(7.2, 8, 8.4), yend = c(7.5, 8.3, 8.7), label = c("***", "***", "***"), label.x = c(1.5, 2.5, 2), label.y = c(7.7, 8.5, 9)) ggplot(iris, aes(x = Species, y = Sepal.Length)) + geom_boxplot() + geom_text(aes(x = label.x, y = label.y, label = label), d.signif) + geom_segment(aes(x = x, xend = x, y = y , yend = yend), d.signif) + geom_segment(aes(x = x, xend = xend, y = yend, yend = yend), d.signif) + geom_segment(aes(x = xend, xend = xend, y = y, yend = yend), d.signif)
ということスッキリかけるようになった。
このデータフレームを使う方法は次のように、fill
じゃなくてfacet_grid
で特定の要因で層化する(パネルを並べるイメージ?)ようなグラフのときに更に生きてくる。具体例として、
ggplot(ToothGrowth, aes(x = supp, y = len)) + geom_boxplot() + facet_grid(cols = vars(dose))
というありがちな感じの、グラフでdose0.5で**をつけようと思ったときに、次のようにするとうまく行かない。
ggplot(ToothGrowth, aes(x = supp, y = len)) + geom_boxplot() + facet_grid(cols = vars(dose)) + geom_text(aes(x = 1.5, y = 31, label = "**")) + geom_segment(aes(x = 1, xend = 1, y = 28, yend = 30)) + geom_segment(aes(x = 1, xend = 2, y = 30, yend = 30)) + geom_segment(aes(x = 2, xend = 2, y = 28, yend = 30))
なんでそうなってしまうかというと、パネルごとにx座標1と2が存在するので、それぞれのパネルごとに同じ位置に描画されてしまう。みんな同じだったらうまーって感じなのだが、今回はdose2で有意な差がなかったので困る。さらにそもそも位置も悪い。
そこでデータを与えてエステティックを描く手法が生きてくる。データにdoseを入れてパネルを選択するようにするとうまくいく。
d.signif <- data.frame(x = c(1, 1), xend = c(2, 2), y = c(28, 30), yend = c(30, 32), label = c("**", "**"), label.x = c(1.5, 1.5), label.y = c(31, 33), dose = c(0.5, 1.0)) ggplot(ToothGrowth, aes(x = supp, y = len)) + geom_boxplot() + facet_grid(cols = vars(dose)) + geom_text(aes(x = label.x, y = label.y, label = label), d.signif) + geom_segment(aes(x = x, xend = x, y = y , yend = yend), d.signif) + geom_segment(aes(x = x, xend = xend, y = yend, yend = yend), d.signif) + geom_segment(aes(x = xend, xend = xend, y = y, yend = yend), d.signif)
他のパッケージの用意している関数を使う場合ではggplot2
の仕様が変更になって使えなくなる可能性があるけれど(そうそうないと思うけど)、ggplot2
が自前で準備している関数なので仕様変更で使えなくなるという心配がほぼないのもいいところ。
ggsignif::geom_signif
を使って有意差のマークをグラフ上に描く
ggsignif
はggplotで描いたグラフに有意差のマークを書くためのgeom_signif
関数を提供するパッケージ。
github上に開発版があるが、CRANの安定版と変わらない。Githubを見ると今後機能追加はしていかない方針の様子。
install.packages("ggsignif")
# ggplotとggsignifのパッケージを読み込む library(ggplot2) library(ggsignif)
前の説と同じIRISデータでとりあえずversicolorと、virginicaの2群間に有意差の米印をつけてみる。
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # geom_signifで印をつける geom_signif(comparisons = list(c("versicolor", "virginica")), map_signif_level=TRUE)
位置は指定しなければグラフの上となる。線の位置もいい塩梅に、グラフの描画範囲もいい塩梅に調整してくれる。
geom_signif
はcomparisons
の引数で比較したい2群を名前でリストで渡して指定する。map_signif_level
の引数でTRUEとすると米印、FALSE(デフォルト)だとp値が直接数値で描画される。
デフォルトだとp値の計算にはwilcox.test
が使われる。testの引数で使いたい2群比較でlistでp.valueを持つ検定方法を指定するとよいので例えばt検定使いたかったら、
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # geom_signifで印をつける geom_signif(comparisons = list(c("versicolor", "virginica")), map_signif_level=TRUE, test = "t.test")
とするといい。
米印だと違いがわからないので、米印じゃないバージョンで2つやってみる。
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + geom_signif(comparisons = list(c("versicolor", "virginica")))
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # geom_signifで印をつける geom_signif(comparisons = list(c("versicolor", "virginica")), test = "t.test")
ちゃんと結果が違う。
これはlibrary()
で適当な2群間比較の検定を読み込んでいると、その他の検定もつかえる。
例えば、Rの組み込み関数のwilcox.testはタイがあると正しいp値を計算できないということが知られているので、代替案としてexactRankTests::wilcox.exactを使いたかったら
library(exactRankTests)
## Package 'exactRankTests' is no longer under development.
## Please consider using package 'coin' instead.
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # geom_signifで印をつける geom_signif(comparisons = list(c("versicolor", "virginica")), test = "wilcox.exact")
とするといい。今回の例だとタイが無いので普通にやった結果と変わらないけれど、特殊な検定を使いたかったらこの方式で計算するのがいい。library名::関数
では残念ながら読めないのでtest = "exactRankTests::wilcox.exact"
と指定してもうまく行かない。
2群間の比較しかできないので多重比較の結果で行けないので微妙なところがあるのは否めない。
annotationの引数で入れる文字(記号)を指定できる。
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + # geom_signifで印をつける geom_signif(comparisons = list(c("versicolor", "virginica")), annotations = "+++")
高さはデフォルトだとグラフの上部で決まった位置になっているので、3群を何も設定しないで米印をつけようとすると
ggplot(data = iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + geom_signif(comparisons = list( c("setosa", "versicolor"), c("versicolor", "virginica"), c("setosa", "virginica")), map_signif_level=TRUE)
とみんな同じ高さになってしまう。離れたところどうし(4群以上で1-2群、3-4群だけ印をつけるなど)なら問題ないのだけれども、こんな感じでくっついてしまう。
y_position
の引数でそれぞれ高さをかぶらないよう調整するといい。さらに調整した高さで勝手に描画範囲がいい塩梅に調整される。
ggplot(iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + # これはデフォルトで入る高さにしておく geom_signif(comparisons = list(c("versicolor", "virginica")), map_signif_level=TRUE) + # これはy7.5の高さにしておく geom_signif(comparisons = list(c("setosa", "versicolor")), map_signif_level = TRUE, y_position = 7.5) + # これはy8.5の高さにしておく geom_signif(comparisons = list(c("setosa", "virginica")), map_signif_level = TRUE, y_position = 8.5)
この書き方だとgeom_signifを複数回書かなければいけないので面倒だが、
ggplot(data = iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + geom_signif(comparisons = list( c("setosa", "versicolor"), c("versicolor", "virginica"), c("setosa", "virginica")), map_signif_level=TRUE, y_position = c(7.5, 8.1, 8.5))
とy_psitionにc()で高さのベクトルをcomparisonsの順番に合わせて与えても同じ図が出力される。これならgeom_signifを複数書かなくてもいいが、デフォルトの高さで入れられなくなるので、全て任意の位置を指定しなければいけない。
基本的にはコを左90°に回した形になっているが、tip_lengthでマイナス数値を指定するとコを右90°に回した形にできる。あとは、vjustで米印の位置を調整でき、マイナス数値で上にプラス数値で下につくようになる。この要領で、
ggplot(iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + geom_signif(comparisons = list(c("versicolor", "virginica")), map_signif_level=TRUE) + geom_signif(comparisons = list(c("setosa", "versicolor")), map_signif_level = TRUE, y_position = 7.5) + # これを下向きにする geom_signif(comparisons = list(c("setosa", "virginica")), map_signif_level = TRUE, # 位置をy4にする y_position = 4, # 足?的なやつの長さを-0.03にして上向きに線が出るようにする tip_length = -0.03, # 米印の位置を横線より下に出るようにする vjust = 3)
でできる。
ただし線の下に米印をつけるように設定すると、あまり遠くに設定すると描画範囲からはみ出てしまう。あくまで自動調整してくれる描画範囲はもともとデフォルトで描画されるだろう範囲のみとなっている。 必要があればylim()
で調整する必要がある。
しかし、fillなどでグラフを層化したときには、項目を選んで選択するということがうまくできなくなる。
ggplot(ToothGrowth, aes(x = factor(dose), y = len, fill = supp)) + geom_boxplot() + geom_signif(comparisons = list(c("OJ", "VC")))
## Warning: Computation failed in `stat_signif()`:
## TRUE/FALSE が必要なところが欠損値です
となり意図した図が描けない(マークが付いていない図しか出力されない)ので、せっせと座標を指定して、計算も自前でやってannonations
の引数で入れてやる。
ggplot(ToothGrowth, aes(x = factor(dose), y = len, fill = supp)) + geom_boxplot() + # これは2つに分けてもいい # annonationsで入れたい文字列を入れる、今回はわかりやすいように米印とプラスでかき分けている geom_signif(annotations = c("**", "++"), # yの位置を調整する。最初の方が上のannonationsの最初の選択したほう、以下同じ。 y_position = c(23, 30), # 有意差のコの字の始点のx座標を指定する。 xmin = c(0.81, 1.81), # 有意差のコの字の終点のx座標を指定する。 xmax = c(1.19, 2.19))
足の長さ的なのも調整できるが、デフォルトでもいい感じなるし、米印の位置も線の始点と終点のちょうど真ん中で勝手に調整してくれるので、ggplot2のgeom_textやgeom_segementでやるよりも入力項目が少なくて済むことが多いと思う。
ちなみにfacet_grid
で分けた場合はこうなる。
ggplot(ToothGrowth, aes(x = supp, y = len)) + geom_boxplot() + geom_signif(comparisons = list(c("OJ", "VC")), map_signif_level = TRUE, test = "t.test") + facet_grid(cols = vars(dose))
ちなみにその2として、doseの0.5と1とで比較したいと言う場合にはこうなる。
ggplot(ToothGrowth, aes(x = factor(dose), y = len, fill = supp)) + geom_boxplot() + geom_signif(comparisons = list(c("0.5", "1")), map_signif_level = TRUE, test = "t.test")
他にもコの字のグラフに伸びる足のような線をtip_length長くしたり、いろいろ調整することができるのでオプション駆使すればかなり融通は効く。
ggplot(data = iris, mapping = aes(x = Species, y = Sepal.Length)) + geom_boxplot() + geom_signif(comparisons = list(c("setosa", "versicolor")), map_signif_level=TRUE, # tip_lengthで長さを調整できる、最初に描いたほうが左側、次が右側になる。値を一つだけにすると両方とも同じ長さで調整できる tip_length = c(0.5, 0.2))
しかしそんなに微調整したいんだったら、geom_segment
やgeom_text
でお絵描きして調整したらいいんじゃないかとも思わなくもない。
エステティックも使えるけれど、どんどん操作が込み入っていくので今回はここまで。基本的にはgeom_textやgeom_segmentでできるようなことと似ているが、ちょっと違うところもあって微妙に使いづらい気がする、基本的にはgeom_signifで準備された引数で処理していくのがいいと個人的には思う。
ggpubr::stat_compare_means
を使って有意差のマークをグラフ上に描く
ggpubr
はggplot2ベースでいい感じのグラフを簡単に記載できるようにした関数群を含むパッケージとなっている。
githubに開発版があるが、CRANに安定版があるので
install.packages("ggpubr")
でインストールする。
# ライブラリで読み込む library(ggplot2) library(ggpubr)
ggpubr
でグラフを描く関数はggplot2
とはやや記法が異なるが、描いたグラフはggplot2
で描いたグラフと同じようにgeom_*
で装飾できる。逆に、ggpubr
で用意されている装飾用の関数のstat_compare_means
などはggplot2
で描いた図も装飾できる。
ggpubr::stat_compare_means
はggsignif::geom_signif
を使っており、ggsignif::geom_signif
の拡張的な感じ。依存関係でggpubr
をインストールするとggsignif
もインストールされ、ggpubr
をロードするとgeom_signif
を使えるようになるし。
ちなみに、ggpubr
は依存関係がggsignif
に比べて多いのでそういうのを気にするなら選択肢として微妙かもしれない。
とりあえず、ggplot2で描いたIRISデータで3群それぞれに米印を付ける。
ggplot(iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + stat_compare_means(comparisons = list( c("setosa", "versicolor"), c("versicolor", "virginica"), c("setosa", "virginica")), # labelで"p.signif"とすると米印、デフォルトだとp値のままで出力される。 label = "p.signif")
ggsignif
ではうまく行かなかったコードに近いものでもきれいに描ける。
デフォルトだと米印がggsignif
よりも一個多い。
comparisonsのリストは米印をつけたい下からの順番で指定しないと見た目変な描画になってしまう。
ggplot(iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + stat_compare_means(comparisons = list( c("setosa", "versicolor"), c("setosa", "virginica"), c("versicolor", "virginica")), label = "p.signif")
デフォルトで用いられるのはwilcox.test
なので変えたかったらmethod
の引数で他の2群比較の検定を指定するといい。ggsignif
同様、library
で読み込んでいるとその他のライブラリで提供される検定も使える。
例えばt検定を使いたかったら
ggplot(iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + stat_compare_means(comparisons = list( c("setosa", "versicolor"), c("versicolor", "virginica"), c("setosa", "virginica")), label = "p.signif", method = "t.test")
とする。
米印はsymnum.args
で調整する。ggsignif
だともし出力結果で制御したかったらif文を書かなければいけないと思うのでその点は楽(実用的かどうかは微妙)。
ggplot(iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + stat_compare_means(comparisons = list( c("setosa", "versicolor"), c("versicolor", "virginica"), c("setosa", "virginica")), label = "p.signif", symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")))
これで、ggsignif
と同じような出力結果に近づく。
そもそも、**じゃなく†使うよ(あるいは他の記号も含めて)、p値0.1から印つけたいよという派も簡単に設定できる。
ggplot(iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + stat_compare_means(comparisons = list( c("setosa", "versicolor"), c("versicolor", "virginica"), c("setosa", "virginica")), label = "p.signif", symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("§", "‡", "†", "*", "ns")))
ただし、§なんかは縦長なのでデフォルトだと位置が微妙なので調整が必要。
印の位置はvjust
で調整できる。マイナス数値で上に、プラス数値で下に移動する。
ggplot(iris, aes(x=Species, y=Sepal.Length)) + geom_boxplot() + stat_compare_means(comparisons = list( c("setosa", "versicolor"), c("versicolor", "virginica"), c("setosa", "virginica")), label = "p.signif", symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 1), symbols = c("§", "‡", "†", "*", "ns")), # ほんの僅かに上に上げる。 vjust = -0.2)
あとは基本的にはできることは、ggsignifとよく似ている。
面白いのは対応のある検定を指定できる点。組み込みデータのsleepデータで対応のあるデータをプロットして線でつないだ箱ひげ図を描いてみる。ggplot
で描くのがちょっと手がかかるのでggpubr
のggpaird
で描く。
ggpaired(data = sleep, x = "group", y = "extra", line.color = "gray", xlab = "group", ylab = "extra")
簡単にそれっぽいのが描ける。 ここで、
with(sleep, t.test(extra[group == 1], extra[group == 2]))
##
## Welch Two Sample t-test
##
## data: extra[group == 1] and extra[group == 2]
## t = -1.8608, df = 17.776, p-value = 0.07939
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.3654832 0.2054832
## sample estimates:
## mean of x mean of y
## 0.75 2.33
with(sleep, t.test(extra[group == 1], extra[group == 2], paired = TRUE))
##
## Paired t-test
##
## data: extra[group == 1] and extra[group == 2]
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.4598858 -0.7001142
## sample estimates:
## mean of the differences
## -1.58
となるので、対応のあるt検定ではちゃんと有意な差があるので、そのことを米印をつけたいが、
ggpaired(data = sleep, x = "group", y = "extra", line.color = "gray", xlab = "group", ylab = "extra") + stat_compare_means(comparisons = list(c("1", "2")), label = "p.signif", method = "t.test")
で有意差なしとなってしまう。そこでparied
をTRUEにすると対応のある検定をしてくれる。
ggpaired(data = sleep, x = "group", y = "extra", line.color = "gray", xlab = "group", ylab = "extra") + stat_compare_means(comparisons = list(c("1", "2")), label = "p.signif", method = "t.test", paired = TRUE)
雑感
p値を計算してくれる機能を持つ有意差の米印マークをつけてくれる関数geom_signif
やstat_compare_means
について、その関数が計算できる範囲の検定でうまくはまっているときはうれしいし(たとえば、純粋に2群比較だったり、sleepデータで対応のある検定できるときなど)気持ちがいいが、IRISデータの例ように本来多重比較すべきものも単に2群比較になってしまうというあまりうまくない展開になってしまってくる。
たまたま米印上は同じ結果になっていると図としての成果物としては問題ないのかもしれないけれど、個人的には危うさをはらんでいるなと思う。
最終成果物として、もうこれ以降絶対変更しないというものだったらいいけど、後でデータが変わったりしたときに多重比較してないがゆえなどの不適切な米印がつきかねないなと思っている。
基本的に米印を付ける段階なんて言うのは最後に検定結果を報告する段階にあると思うけど、膨大な量じゃなければgeom_signifやstat_compare_meansが用意している検定結果からマークをつけるのではなく、人力で米印を決めるか、判定間違いやつけ忘れを防ぐためにも自動化したいのであれば少し手間だけれど別途if文で米印のlabelsを決めてやった方が安全だと思う。
関係のない雑感として、ggplot()
だと引数の一番最初がdataで、mappingが次だけど、geom_*
だとmappingの方が先なので明示しない場合にggplot()
と同じように書くとaestheticおかしいよとエラーがでるのが地味にやりにくい。
ほかの関数的にも、dataがaestheticのあとなのは妥当なんだけど、普段ggplot()にいろいろ書いている派なので混乱する。 ggplot()
はggplot()
でパイプでデータを受けとる設計なのかもしれないのでこれはこれで妥当なんだろうなぁと思った。
あとは極めてどうでもいい雑感としてstat_compare_means
はデフォルトがwilcox.testなのでmeanじゃなくmedianじゃないか(それも正確にはどうなのという話もあるが)、そもそもmeanとかいらなかったのでは、あえてつけたくてどっちもいいよというならお茶を濁す程度にmとかしておけばと思った。 関数名はパッケージ作者の専権事項だし、こんな機能を使おうと思うユーザーが勘違いして平均値比較の検定(t.testがデフォルト)と思いこむこともそうなのでどうでもいいけど。