箱ひげ図にp値と棒をいっぱいかく

箱ひげ図にp値とどの群を比較したか棒をつけている図を論文でよく見かけるが、Rにはそれを楽に描いてくれるような機能がない。
というのも、多重比較に関わるから難しいし、棒をつけたらつけたでゴチャゴチャするからなのだろう…
ただ、この図が描きたいという要望が多かったのでやってみる。
データはここからパクった。

effect <- c(80, 73, 80, 82, 74, 73, 68, 81, 85, 85, 93, 88)    # データ
group  <- c( 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3)    # 群
#素直に描く
boxplot(effect ~ group)

#本来なら、正規分布性をbartlett.testで確認するべきなんだろうけど
k1 <- kruskal.test(effect ~ group)

#ペアワイズなt検定の繰り返し
p1 <- pairwise.t.test(effect, group, p.adj = "bonf")
pv1 <- p1$p.value[lower.tri(p1$p.value, diag=TRUE)]

b1 <- boxplot(effect ~ group) #boxplotの統計量の各座標を取得しておく
boxplot(effect ~ group, ylim=c(65, 100))
segments(x0=min(group), y0=min(b1$stats[1, ]) - 1, x1=max(group)) #全体のp値
text(2, min(b1$stats[1, ]) - 1, paste("p =", round(k1$p.value, 3)), adj=c(0.5, 1.5))

segp <- combn(3, 2)
segy <- c(max(b1$stats[, 1:2]), 96, max(b1$stats[, 2:3])) + 1.5 #ここは手動で設定しておかないと、重なる場合がある
for(i in seq(ncol(segp))){
	segments(x0=segp[1, i], y0=segy[i], x1=segp[2, i]) #群間比較のp値
	segments(x0=segp[1, i], y0=segy[i], y1=segy[i]-1)
	segments(x0=segp[2, i], y0=segy[i], y1=segy[i]-1)
	text(mean(segp[, i]), segy[i] + 1, paste("p =", round(pv1[i], 3)), adj=c(0.5, 0))
}