2015年9月7日追記
中澤先生のツイートで知りましたが、信頼区間をプロットする関数を含んでいるパッケージがすでに存在するようです。
http://www.inside-r.org/packages/cran/plotrix/docs/plotCI
こちらを使うことをお勧めします。
誰がこんな関数使うんだよシリーズ！
以下をコンソールにコピペしてください。
PlotCI <- function(TB,las=0, xlab="", ylab="" ,type="p",pch=1){
#基本のプロット
plot(TB[,1], ylim=c(min(TB),max(TB))
,type=type,axes=F,xlab=xlab,ylab=ylab,pch=pch)
axis(2,las=las)
axis(1,1:nrow(TB),rownames(TB),las=las)
#信頼区間のバーを描く
arrows(1:nrow(TB), TB[,2], 1:nrow(TB), TB[,3], length=.05, angle=90, code=3,col="grey")
}
使い方
関数本体
PlotCI (table , type, las, xlab, ylab, pch)
- table
- 1列目に点推定値（平均や回帰係数など）、2列目に信頼区間の上限、3列目に下限を入れた、n ×3行の表（ちなみに行列でもいける）
- type
- プロットの仕方。plot()関数に準じる。デフォルトは、type="p"
- las
- 見出しの向き。デフォルトは las = 0
- ylab/xlab
- Y軸・X軸のラベル。デフォルトでは何も描かない
- pch
- プロット点の形状。デフォルトは "○"
使用例
> #回帰係数っぽいもの乱数で10個出す
> b <- rnorm(10)
>
> #標準誤差っぽいものを乱数で10個出す
> se <- 0.5*rnorm(10)
>
> #95%信頼区間を求めるため、標準誤差に1.96をかけたものを回帰係数に足したり引いたり
> #ちなみに、1.96というのは自由度無限大と仮定してます（笑）
> 上限 <- b+se*1.96
> 下限 <- b-se*1.96
>
> # n列×3行の表にまとめる
> X <- cbind(b, 上限 , 下限)
>
> #行に適当に名前を付ける
> rownames(X) <- LETTERS[1:10]
>
> #表はこんな感じです
> X
b 上限 下限
A 0.4570674 -0.06099424 0.9751291
B 1.1403102 0.50063212 1.7799883
C -0.3636855 -1.13794614 0.4105752
D -0.7133917 -1.79882090 0.3720375
E 1.5478229 0.16446766 2.9311782
F -0.5235648 -0.29019635 -0.7569333
G 0.7012099 1.92979771 -0.5273778
H 1.0280072 1.74644329 0.3095711
I 0.8187337 -0.23824256 1.8757099
J 0.6582293 -0.79155813 2.1080168
>
>
> #プロットする
> PlotCI(X)
> abline(h=0,lty=2) #ついでに y=0 の線を引く
> 
さらに・・・・
別のペアを重ね書きする PointsCI
つまり、points()関数をベースにした信頼区間図示付け足し用です。
関数本体
以下をコピペ
PointsCI <- function(TB,hpos=0.1 , type="p",pch=2,lty=1){
#基本のプロット
points(x=1:nrow(TB)+hpos,y=TB[,1]
,type=type,pch=pch,lty=lty)
#信頼区間のバーを描く
arrows(1:nrow(TB)+hpos, TB[,2], 1:nrow(TB)+hpos, TB[,3], length=.05, angle=90, code=3,col="grey")
}
使用法
PointsCI(TB, hpos , pch, type)
- hpos
- 重ねるプロットを水平方向にどれだけずらすか。デフォルトは、hpos = 0.1
あとはPlotCI() と同様。
使用例
> # jitter()で、b と se にゆらぎを加え、微妙に違う b と se をゲット
> b2 <- jitter(b,100)
> se2<- jitter(se,100)
> 上限2 <- b2+se2*1.96
> 下限2 <- b2-se2*1.96
> X2 <- cbind(b2, 上限2 , 下限2)
> X2
b2 上限2 下限2
[1,] 0.54787711 -0.16677589 1.2625301
[2,] 1.52170702 1.01240594 2.0310081
[3,] -0.67233139 -1.88613187 0.5414691
[4,] -0.82281372 -2.22684888 0.5812214
[5,] 1.57784357 0.04915702 3.1065301
[6,] -0.38026040 -0.27633554 -0.4841853
[7,] 1.37210506 2.51883530 0.2253748
[8,] 0.38760203 0.61830278 0.1569013
[9,] 0.59942531 0.06430740 1.1345432
[10,] -0.06895293 -1.03655354 0.8986477
>
> #最初の図に書き加える
> PointsCI(X2, pch=2,hpos =.3, type="p")
>
