Rで生物系の研究者がよく使うような検定をしつつ論文っぽい図を描きたい。
(R3.6.2)を使用
Rでは「統計解析」と「図示」という側面があると思うが、今回はそれらを繋ぐことを目指す。
具体的には、統計解析によって有意に差があるとわかっても、普通のやり方では有意差マーク（"*"など）は図に表示されないので、
今回は、有意差があったときのアスタリスク（*）を「有意な度合いに従って」つけることを目標にやってみる。
##サンプルデータの取得
Rにもとから入ってるirisっていうデータを今回は使う。
#サンプルデータの取得
data <- iris
#入ったかの確認
head(data)
t検定
をやりたいのだが、irisデータセットは3群なので、ちょっとt検定から離れて2群だけ取り出す作業をする。
そのためにdplyrを読み込む（もともと2群データを持ってる人は飛ばしていい）。
行の抽出はdplyrパッケージのfilter関数とかでできます。
#パッケージのインストールと読み込み
if (!require(dplyr))install.packages('dplyr')
library(dplyr)
#setosaとversicolorだけ取り出したやつをdata2とする
data2 <- data %>% filter(Species %in% c("setosa", "versicolor"))
やっとt検定
t.test(Sepal.Length ~ Species, data = data2, var.equal=T)
出力結果↓（p < 0.05であることがわかる。）
Two Sample t-test
data: Sepal.Length by Species
t = -6.9334, df = 48, p-value = 9.304e-09
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.1042354 -0.6077646
sample estimates:
mean in group setosa mean in group versicolor
5.024 5.880
p値だけを取り出す。
せっかくなので、p値をcsvファイルにでも保存しておく（研究ではこういう作業が重要だと思う）。
#tにさっきの検定を格納
t <- t.test(Sepal.Length ~ Species, data = data2, var.equal=T)
#p値の取り出し
tp <- t$p.value
#csvファイルへの出力
write.csv(tp, file = "t.qiita.csv")
##p値によって適切なアスタリスク（*）を返す関数を作ってしまう
有意差があるやつだけ手動でつけてもいいのだが、p < 0.05（*）やp < 0.001（***）のようにアスタリスクの数を変えて有意な度合いを表現したい時に、めんどくさいし、Rにやらせたほうがミスが起きにくいのかなと思う。
てことで関数を作った。
#sigという関数を作り、pが以下の条件のときに適切な数だけ（*）を返すような関数を入れてみた
# p > 0.1 "", 0.1 ≧ p > 0.05 ".", 0.05 ≧ p > 0.01 "*", 0.01 ≧ p > 0.001 "**", 0.001 > p "***"
sig <- function(a) {
if (a > 0.1) {
return("")
} else {
if ((a <= 0.1)&&(a > 0.05)) {
return(".")
} else {
if ((a <= 0.05)&&(a > 0.01)) {
return("*")
} else {
if ((a <= 0.01)&&(a > 0.001)) {
return("**")
} else return("***")
}
}
}
}
#できるか確認してみる
sig(tp)
#結果↓（tp < 0.001なので合っている）
[1] "***"
追記：このsig関数について、もっとスタイリッシュに書く方法を、@doikoji 様が記事を書いてくださっているので、参考にしてください。
Rの多重ifブロックを良い感じに改良してみる
##論文っぽい箱ひげ図を描く
プロットするためにggplot2パッケージを（インストールし、）読み込む
if (!require(ggplot2))install.packages('ggplot2')
library(ggplot2)
基本のプロットgと、それを論文っぽくカスタマイズしたg2を作る
→ 参考にさせていただいたサイト
#基本
g <- ggplot(data2, aes(x=Species, y=Sepal.Length, fill=Species)) + #data2でx軸Species、y軸Sepal.Length、色Species
geom_boxplot(lwd=3) + #画像にすると線が細くなりがちなので太めに（デフォルトは1）
theme(panel.grid = element_blank()) #いったん枠線を無くす
g #確認
論文っぽくカスタマイズ
g2 <- g +
scale_fill_manual(values = c("gray40", "white")) + #論文って白黒のときあると思うので、色を指定
labs(y="Sepal length (cm)", x=NULL) + #y軸の名前いい感じに。x軸は凡例だけあればいいかなと思うので消す。
theme(axis.text=element_text(size=37, face = "bold", colour = "black"), #軸の数字の大きさ、字体、色
axis.title=element_text(size=40,face="bold", colour = "black"), #軸のタイトルの大きさ、字体、色
axis.line = element_line(size = 1.2, colour = "black"), #軸の大きさ、色
panel.grid.major = element_blank(), #以下、補助線や背景を消した
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position = "none"#, #凡例消した。つけたい時はここの先頭に#つけて下の4行の#を取って適宜調整する
#legend.position=c(0.80,0.80),
#legend.background = element_rect(fill = "white", colour = "black"), #凡例の背景や色
#legend.title = element_text(size = 37, face = "bold"), #凡例タイトルや中身のテキストの大きさや字体の調整
#legend.text = element_text(size = 37, face = "bold")
)
g2 #確認
ここで、グラフにタイトルを付け加えたい。
annotate関数で付けることができるが、座標を指定しなくてはならず、y軸用のデータの数値の大きさによって毎回高さを変えるのは面倒である。
そこで、y軸データの大きさに依存したyRoofを定義する。
→ 参考にさせていただいたサイト
#y軸データの最大値の1.2倍の値の小数第一位を四捨五入
yRoof <- round(max(data2$Sepal.Length)*1.2,1)
yRoofを使ってタイトルをつける。
g2 <- g2 +
annotate("text", x=1.5, y=yRoof, #タイトルをつけた。1.5は1個めと2個目の箱ひげの間を指す。
label="Sepal length",
vjust=0.5,
size=15,
fontface="bold")
何もしなければは有意差マーク（*）等は出てこない。
t検定で有意差がなかったときはこのあたりで保存すればいいと思う。
画像の保存は、ggplot系ならggsave関数でできる。plot系はまた別。
（この～系って呼び名は私が勝手につけたけど、理にかなってるとは思う）
ggsave(filename = "p.qiita.png", plot = g2, width = 9, height = 9, dpi = 300)
有意差がある場合、先ほど作ったsigとyRoofを使って、アスタリスクをちょうどいい位置に配置する。
やってみたらyRoofの0.9～0.93倍あたりがちょうどよかったのでこうしてみた。
g2 <- g2 +
geom_segment(x=1, y=yRoof*0.9, xend=2, yend=yRoof*0.9, size=1) + #geom_segmentで有意差マークの下にある棒を描く
annotate("text", x=1.5, y=yRoof*0.93, label=print(sig(tp)), size=15) #さっき作ったsig関数にp値であるtpを入れる
画像を保存
ggsave(filename = "p.qiita.png", plot = g2, width = 9, height = 9, dpi = 300)
##geom_signifについて（2020.08.02追記）
私が考えるようなことはもう誰かが考えている笑
ggsignifというパッケージがあって、インストールすることで今回やったようなことをできるgeom_signifが使えるようである。
やっていることは同じだが、操作性はこちらのほうがよさそうなので、今後はこっちも使ってみようと思う↓参考
https://note.com/eiko_dokusho/n/nbc9f39e8d31e
(2020.08.04)
やってみた。
ggsignifは環境によってはinstall.packagesでインストールできないらしいので、そのへんは、こちらのサイトを参照。
if (!require(ggsignif))install.packages('ggsignif')
library(ggsignif)
#上のやり方で有意差マークを入れる前（geom_segmentを使う前）のg2を使う
g3 <- g2 +
geom_signif(
y_position=yRoof*0.9, #高さの設定
xmin=1.0, #線の始点の設定
xmax=2.0, #終点の設定
annotation= sig(tp) #さっき作ったsigを使う,
tip_length=0.05, #棒の下向きの線の長さ
textsize=10 #*の大きさ)
このやり方では、geom_segmentとannotateというように2回書かなくていいので便利かもしれない。
私の作ったsig関数を入れることで、より汎用性の高いコードが仕上がると思う。
##まとめ
今回はよく使うt検定からの、論文っぽい図作成という王道的なことをやった。
Rを始めた頃は、*の場所（座標）を逐一指定していたが、有意な度合い（例えば多重比較でのa,b,ab,...など）によって、有意差マークを変えたり、高さを変えたりするのは面倒だと思うので、いい感じにまとめることができたのではないかと思う。
本当はyRoofを定義せずともマークをつけるところまでやってくれる関数も作ってみたいが、そこまでの熱意はなかった笑
今度はTwo-way ANOVAなどにこれを応用したものを投稿していきたい。