geom_pointrangeで推定値と95%信頼区間の図示

右の図みたいに

いくつかの種に対してGLM等で得られた推定値と95%信頼区間を比較できる図を作りたいと思ったので作った。

・Rのggplot2にて描画

・とりあえず描画できればいいだろと思い、スクリプトは雑

・geom_pointrangeってのを使うとできる

・95%信頼区間の計算が別途必要

下みたいなデータセットがあればいける。

name	factor	estimate	SE	p_value	a	signif	order	max	min	est_round
種A	mowing	0.209	0.4691	0.6559		no	1	1.003603	-0.5856	0.21
種A	grazing	0.7708	0.4342	0.0758		no	1	1.506287	0.035314	0.77
種A	fire	0.8621	0.481	0.073		no	1	1.676861	0.04734	0.86
種B	mowing	1.8593	0.8313	0.0253	*	yes	2	3.26743	0.45117	1.86 *
種B	grazing	1.1811	0.5461	0.0305	*	yes	2	2.106133	0.256067	1.18 *
種B	fire	0.3572	0.4326	0.409		no	2	1.089976	-0.37558	0.36
種C	mowing	1.9794	0.7334	0.0069	**	yes	3	3.221698	0.737102	1.98 **
種C	grazing	-0.0081	0.4374	0.985		no	3	0.732807	-0.74901	-0.01
種C	fire	0.0986	0.4196	0.8141		no	3	0.809356	-0.61216	0.1

なぜか枠線が消えてるけど気にしない

なければ下のをRで読み込めばok

#適当なデータフレームを作る準備（あまり重要ではない）

nam<-rep(c("種A", "種B", "種C"),each=3)

fac <- rep(c("mowing", "grazing", "fire"),3)

est <- c(0.2090, 0.7708, 0.8621, 1.8593, 1.1811, 0.3572, 1.9794, -0.0081, 0.0986)

SE <- c(0.4691, 0.4342, 0.4810, 0.8313, 0.5461, 0.4326, 0.7334, 0.4374, 0.4196)

pvl <-c(0.6559, 0.0758, 0.0730, 0.0253, 0.0305, 0.4090, 0.0069, 0.9850, 0.8141)

aaa <-c("","","","*","*","","**","","")

sig <-c("no","no","no","yes","yes","no","yes","no","no")

ord <-rep(c(1,2,3),each=3)

#データフレーム作成

df <-data.frame(name=nam,factor=fac,estimate=est, SE=SE, p_value=pvl, a=aaa, signif=sig, order=ord)

#95%信頼区間の範囲をデータフレームに追加

# 自由度dfと区間pを決めるとqt値ってのが返ってくる。→詳細(外部リンク)

　qt<-qt(df=32, p=0.95)

# 信頼区間の最大値の列を追加

　df$max<-df$estimate+df$SE*qt

# 信頼区間の最小値の列を追加

　df$min<-df$estimate+df$SE*-qt

#推定値とアスタリスクを1つの文字としてまとめる

df$est_round <- paste(round(df$estimate, digits=2), df$a) # round()で推定値を丸める

#データフレームの確認

df　#たぶん上の表みたいになってると思う

#読み込みます

library(ggplot2)

#描画します

f<-ggplot()+　

#推定値0のところに灰色の破線を描画.

geom_hline(yintercept = 0, linetype=2, col="grey50")+

#後で90度回転させるからx軸に種名、y軸に推定値とか信頼区間の最大値・最小値を指定、x軸は何も指定しないと五十音順なのでreorder()で意味ある順序に.

geom_pointrange(data = df, aes(x = reorder(x = name, X = -order), y = estimate, ymin = min, ymax = max, col = signif))+

#ラベルの追加.vjustで高さ合わせ.

geom_text(data = df,aes(x = name, y = estimate, label = est_round), size = 4, vjust = -1) +

#90度回転

coord_flip()+

#今回は複数の要因があるので要因ごとに分けて描画. ※1要因ならいらない.

facet_wrap(~ factor)+

#x軸ラベル

xlab(' ')+

#y軸ラベル

ylab('Estimate')+

#グラフをシンプルにしたいなら.

theme_classic( )+

#凡例の場所"top","bottom","right","left"とか、今回はなし。

theme(legend.position = "none",

#facet_wrap()で出てくるフレームの枠線色と背景色の設定.

strip.background = element_rect(color = "white", fill = "snow3"),

#factorの文字の大きさの設定.

strip.text.x = element_text(size = 11),

#種名の文字の設定.

axis.text.y = element_text(colour = "black", size = 11),

#"Estimate"って文字の設定.

axis.text.x = element_text(colour = "black", size = 11))+

#col=で設定した色の指定（geom_pointrangeの中で設定してる）.

scale_color_manual(values = c("snow3", "black"))

#出力

f

下図のようなものが出力される。

以上