Rでエラーバー付きのプロットを行う

例えば↓こんなサンプルデータがあったとしましょう。属するグループと観測値の組みたいな感じですね。

> x1
   group value
1      A  2.66
2      A  2.24
3      A  2.41
4      A  1.12
5      A  2.48
6      A  2.33
7      A  2.90
8      A  4.35
9      A  3.82
10     A  3.45
11     B  6.65
12     B  4.61
13     B  5.26
14     B  5.24
15     B  7.26
16     B  5.97
17     B  2.90
18     B  2.37
19     B  4.94
20     B  4.08
21     C  5.63
22     C  5.75
23     C  4.87
24     C  8.14
25     C  1.01
26     C  1.58
27     C  3.76
28     C  8.12
29     C  5.95
30     C  4.48

で、グループごとの平均やばらつきにどんな傾向があるか見るために、エラーバー付きでプロットしたいような状況だと考えてください。

こういうときは、gplotsパッケージのplotmeans関数が便利です。

install.packages("gplots")
library(gplots)
plotmeans(value~group, data=x1, connect=F, ylim=c(0,10))

plotmeansによるエラーバー付きプロット

ところで、エラーバーの長さは何を表しているのでしょうか？

> ?plotmeans

でマニュアルの引数のところを見ると、

p     confidence level for error bars. Defaults to 0.95.

と書いてあります。p= で指定することによって、エラーバーの信頼水準を指定することができ、そのデフォルト値は0.95であると書いてあるので、つまりエラーバーは信頼区間を表しているということですね。

確認してみましょう。

グループAだけ抽出して、t.test関数を使ってt検定をし、その95パーセント信頼区間を求めてみます。

> tt <- t.test(x1[x1$group=="A",]$value)
> tt

        One Sample t-test

data:  x1[x1$group == "A", ]$value
t = 9.6159, df = 9, p-value = 4.952e-06
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 2.122944  3.429056
sample estimates:
mean of x
    2.776

2.122944 ～ 3.429056 という結果が出ました。これを先ほどのplotmeansの結果に追記してみます。

abline(h=tt$conf.int[1], lty=2)
abline(h=tt$conf.int[2], lty=2)

plotmeansのエラーバーは信頼区間であることを確認

うん、見事に一致しました。

ということで、plotmeansにおけるエラーバーはt検定における信頼区間を現していて、デフォルトは0.95だけど、p= で指定できるということが分かってもらえたと思います。

で、以降は補足です。

データの形が違うってことがありますよね。例えば↓こんな感じ。

> x2
      A    B    C
1  2.66 6.65 5.63
2  2.24 4.61 5.75
3  2.41 5.26 4.87
4  1.12 5.24 8.14
5  2.48 7.26 1.01
6  2.33 5.97 1.58
7  2.90 2.90 3.76
8  4.35 2.37 8.12
9  3.82 4.94 5.95
10 3.45 4.08 4.48

こんな場合は、plotmeansの入力に合わせて、形を組み替えちゃえばいいと思います。

> # plotmeansの入力データの形式に直す
> x3 <- data.frame(value=unlist(x2), group=rep(names(x2), each=nrow(x2)))
> x3
    value group
A1   2.66     A
A2   2.24     A
A3   2.41     A
A4   1.12     A
A5   2.48     A
A6   2.33     A
A7   2.90     A
A8   4.35     A
A9   3.82     A
A10  3.45     A
B1   6.65     B
B2   4.61     B
B3   5.26     B
B4   5.24     B
B5   7.26     B
B6   5.97     B
B7   2.90     B
B8   2.37     B
B9   4.94     B
B10  4.08     B
C1   5.63     C
C2   5.75     C
C3   4.87     C
C4   8.14     C
C5   1.01     C
C6   1.58     C
C7   3.76     C
C8   8.12     C
C9   5.95     C
C10  4.48     C

行名が気に入らない場合は、↓こんな感じで消してしまいましょう。

> row.names(x3) <- NULL
> x3
   value group
1   2.66     A
2   2.24     A
3   2.41     A
4   1.12     A
5   2.48     A
6   2.33     A
7   2.90     A
8   4.35     A
9   3.82     A
10  3.45     A
11  6.65     B
12  4.61     B
13  5.26     B
14  5.24     B
15  7.26     B
16  5.97     B
17  2.90     B
18  2.37     B
19  4.94     B
20  4.08     B
21  5.63     C
22  5.75     C
23  4.87     C
24  8.14     C
25  1.01     C
26  1.58     C
27  3.76     C
28  8.12     C
29  5.95     C
30  4.48     C

以下は、さらに蛇足。

Rでエラーバーを表示させる方法といえば、私は↓この方法しか知らなかったんですね。

誤差範囲 | Rで棒グラフと折れ線グラフにエラーバーを付ける | バイオスタティスティクス


矢印を書くarrows関数というのがあるんですが、矢尻の角度を矢の本体の棒に対して90度にすることで、エラーバーっぽく見せるというやり方です。矢の長さを決めるのに必要な値は自分で事前に計算しておくという感じで、かなり手作り感覚あふれる描画方法です。

で、Rにはエラーバーを表示する関数はないのかなあと思っていたところで、↓この書籍で、gplotsパッケージとplotmeans関数の存在を知りました。

Rによる統計的検定と推定


ちなみに、この書籍の良いところは、それぞれの手法に対して、Rによる検定コマンドの紹介の前に、ちゃんと可視化の例が載っているところですね。

検定の前にちゃんとデータの可視化を行っている
「Rによる統計的検討と推定」より

上記は、2つの平均値（対応あり）の検定を紹介しているページですが、t.testによる検定の前に、matplotによる可視化でデータの傾向を確認しています。可視化をすっ飛ばして、いきなり検定するとかありえないですもんね。

話を戻すと、前述のリンク先（バイオスタティスティクス）のarrows関数による方法は、それでもやっぱり有用だったりします。

というのも、plotmeansによるエラーバーは信頼区間を現していましたが、エラーバーの長さって他の値で決めることもあるじゃないですか。1σとか。plotmeansではそれはできないのですが、arrowsを使う方法ならば、自分で値を計算するから、任意の指標でのエラーバー付けができますからね。これはこれでいろいろ使えて便利だと思います。