箱ひげ図のノッチの値を取得する

 各群におけるデータのばらつきを視覚的に分析したい場合、箱ひげ図が便利です。また、箱ひげ図にノッチ(切れ目)を入れると、各群の中央値の95%信頼区間が表示されます。そして、2つの群のノッチがオーバーラップしていなければ、それらの群の間に有意差があるとされます。*1 しかしながら、箱ひげ図を目で見るだけでは、ノッチの重なりが微妙な場合もあります。そのようなときは、各群の信頼区間を数値で確認するとよいでしょう。
 まず、お馴染みのirisデータセットを使って、箱ひげ図を描いてみます。以下の例では、Species (= iris[, 5]) によるSepal.Length (= iris[, 1]) の違いに注目しています。

# 箱ひげ図の描画
boxplot(iris[, 1] ~ iris[, 5], col = "pink", notch = TRUE)

 上記のスクリプトを実行すると、以下のような図が出力されます。

f:id:langstat:20150601224856p:plain

 そして、上記のスクリプトを変数に代入することで、箱ひげ図の描画に使われている様々な値を取得することができます。*2

# 変数に代入
x <- boxplot(iris[, 1] ~ iris[, 5], col = "pink", notch = TRUE)
# 変数の中身を表示
x

 上記のスクリプトを実行すると、以下のような図が結果されます。それを見ると、各群の信頼区間が$confに格納されていることが分かります。

$stats
     [,1] [,2] [,3]
[1,]  4.3  4.9  5.6
[2,]  4.8  5.6  6.2
[3,]  5.0  5.9  6.5
[4,]  5.2  6.3  6.9
[5,]  5.8  7.0  7.9

$n
[1] 50 50 50

$conf
         [,1]     [,2]     [,3]
[1,] 4.910622 5.743588 6.343588
[2,] 5.089378 6.056412 6.656412

$out
[1] 4.9

$group
[1] 3

$names
[1] "setosa"     "versicolor" "virginica" 

 因みに、信頼区間の行列 ($conf) に列名と行名をつけるには、以下のような処理を行います。

# 列名を指定
colnames(x$conf) <- x$names
# 行名を指定
rownames(x$conf) <- c("lwr", "upr")
# 結果を確認
x$conf

 上記のスクリプトを実行すると、以下のような表が結果されます。

      setosa versicolor virginica
lwr 4.910622   5.743588  6.343588
upr 5.089378   6.056412  6.656412

*1:(2015年6月3日追記)しかし、ノッチがオーバーラップしていても、有意差がある場合はあります。例えば、irisデータセットのversicolorとvirginicaにおけるSepal.Width (= iris[, 2]) の場合、このようにノッチがオーバーラップしていても、有意差はあります。信頼区間と検定に関しては、Belia, S., Fidler, F., Williams, J., and Cumming, G. (2005). Researchers misunderstand confidence intervals and standard error bars. Psychological Methods, 10(4), 389-396などを参照(PDFはこちら)。

*2:@kazutanさんのツイートで知りましたが、boxplot.statという関数もあります。ただ、全ての群の統計値をまとめて出してはくれないようです。