「統計学関連なんでもあり」の過去ログ--- 045

No.16387　R gplotsライブラリのplotmeansのエラーバー　　【jas】　2012/02/09(Thu) 05:07

お世話になります。Rを用いてグループ毎に平均にエラーバーをつけた図を作りたいと思います。下記の表をDATAに読み込み，gplotsライブラリのplotmeansを用いて，下のようにして作図しました。

library(gplots)
plotmeans(DATA$Value~DATA$Group, n.label=F, connect=F, cex=2.5, pch=19,
 barwidth=3, ci.label=TRUE, xlab="Group", ylab="Value")

この時のエラーバーの大きさは，p=で指定でき，confidence level for error bars. Defaults to 0.95とあります。そのため，この値は1.96 x SDだと思うのですが，たとえばA群は平均が4，標準偏差が1.581であり，エラーバーの上端は図のように5.963にはなりません。また，エラーバーを標準偏差にするためにp=0.682と設定しても，計算が合いません。
この考え方のどこが間違えているのでしょうか。よろしくお願いします。

Group	Value
A	2
A	4
A	3
A	5
A	6
B	7
B	6
B	8
B	9
B	12
B	7
C	15
C	16
C	18
C	19
C	21
C	22

No.16389　Re: R gplotsライブラリのplotmeansのエラーバー　　【青木繁伸】　2012/02/09(Thu) 08:31

plotmeans のプログラムを参照すれば疑問は解決出来るでしょう。該当部分は以下の通り。

        myvar <- function(x) var(x[!is.na(x)]) # 欠損値を除いて不偏分散を計算する関数の定義
        vars <- sapply(split(mf[[response]], mf[[-response]]), myvar) # 各群の不偏分散
        vars <- ifelse(vars < (minsd^2), (minsd^2), vars) # minsd はデフォルトで 0 なので，デフォルト時にはもとの vars のまま
        if (use.t) # t 分布を使うとき（デフォルト）
            ci.width <- qt((1 + p)/2, ns - 1) * sqrt(vars/ns) # 普通の定義
        else ci.width <- qnorm((1 + p)/2) * sqrt(vars/ns) # 正規分布を使うとき

実行結果は以下の通り

ns :
A B C 
5 6 6 
qt((1 + p)/2, ns - 1) :
       A        B        C 
2.776445 2.570582 2.570582 
sqrt(vars/ns) :
        A         B         C 
0.7071068 0.8724168 1.1180340 
vars :
       A        B        C 
2.500000 4.566667 7.500000 
minsd = 0 
ci.width :
       A        B        C 
1.963243 2.242619 2.873998 
> tapply(DATA$Value, DATA$Group, mean) # 各群の平均
        A         B         C 
 4.000000  8.166667 18.500000 
> tapply(DATA$Value, DATA$Group, sd) # 各群の標準偏差
       A        B        C 
1.581139 2.136976 2.738613

No.16395　Re: R gplotsライブラリのplotmeansのエラーバー　　【jas】　2012/02/09(Thu) 15:16

元のプログラムを参照するように指摘するだけではなく，コメントをつけていただきありがとうございました。よくわかりました。
エラーバーを標準偏差にするのもci.widthを修正しうまくいきました。
ありがとうございました。