箱ひげ図（boxplot）は、表にまとめると視認性が悪い複数のサブサンプルの四分位と最大最小の値を、簡潔に見やすく表示することで、分布の特徴の比較が簡単にできるプロットです。今では中高でも教えているぐらい一般化しています。
1 箱ひげ図（boxplot）の構造
プロットする前にRの箱ひげ図がどういうものかを確認しておきましょう。箱ひげ図はサブサンプルごとに以下の図を作成します。
第1四分位（Q1; 25%パーセンタイル）、第2四分位（Q2; 中央値）、第3四分位（Q3; 25%パーセンタイル）は教科書通りですが、最大と最小が異常値（extreme dataもしくはoutlier）を除いて計算されます。
異常値は\(\{x \mid x \lt \text{Q1} -
\text{range}\cdot\text{IQR}\} \cup \{x \mid x \geq \text{Q3} +
\text{range}\cdot\text{IQR}\}\)で定義されます。rangeはデフォルトでは1.5で、0を指定すると\(\infty\)として解釈されます。\(\text{IQR}\)は\(\text{Q3}-\text{Q1}\)です。
2 箱ひげ図をプロットしてみる
習うより慣れろということで、プロットしてみましょう。
2.1 データセット
乱数から個体数4, サンプルサイズ160のパネル・データセットを作りましょう。
set.seed(1212)
df01 <- with(list(n = 160, n1 = 15, nc = 4), {
x <- runif(n, min=0, max=5)
z <- c(rep(1, n1), rep(2, n/2 - n1), rep(3, n/4), rep(4, n/4))
zv <- runif(nc, min=-5, max=5)
y <- 1 + 2*x + zv[z] + rnorm(n)
data.frame(y, x, z=factor(sprintf("group-%s", z)))
})2.2 シンプルにプロット
単純に使っても、概ね期待通りのモノが出力されます。
r_boxplot <- boxplot(y ~ z, data = df01)
上のプロットはモデル式でデータを指定しましたが、
with(df01, boxplot(
y[z=="group-1"],
y[z=="group-2"],
y[z=="group-3"],
y[z=="group-4"]))と、ベクトルを並べていくこともできます。ただし、namesを指定しないと、グループ名が入りません。
3 戻り値を確認
boxplotはinvisibleで、
- 正常値の最小・四分位・正常値の最大（
stats） - 観測数（
n） - 95%信頼区間の近似
notch\(= \pm 1.58 \text{IQR} / \sqrt{n}\)（conf） - 異常値（
out） - 異常値が属するグループ番号（
group） - グループ名（
names）
を要素として含むリストを戻します。
print(r_boxplot)$stats
[,1] [,2] [,3] [,4]
[1,] -1.67134921 -3.3336851 -5.0511667 -0.3469611
[2,] -0.02562462 -0.0885454 -0.5182996 1.6359863
[3,] 2.25986714 3.0815960 1.4288402 3.8609982
[4,] 4.47099636 4.9971122 3.5222448 6.5509501
[5,] 6.98958634 8.4444394 6.7823874 10.2391417
$n
[1] 15 65 40 40
$conf
[,1] [,2] [,3] [,4]
[1,] 0.4254515 2.084935 0.4194317 2.633142
[2,] 4.0942828 4.078257 2.4382487 5.088854
$out
numeric(0)
$group
numeric(0)
$names
[1] "group-1" "group-2" "group-3" "group-4"plot = FALSEをつけると、描画しないで戻り値を戻します。
4 箱ひげ図に載せる情報の追加/調整
表示する情報を増やしたり、調整したりできます。
4.1
notchを表示する
サブサンプルが同一かの目安に使うnotchを表示できます。なお、統計的仮説検定（e.g. t検定）と対応しているわけではなく、尤度原理にも即していないので、参照するときは注意深く行ないましょう。
r_boxplot <- boxplot(y ~ z, data = df01, notch = TRUE)
4.2 サンプルサイズに比例した太さにする
サブサンプルごとに観測数が違う場合、varwidth = TRUEをオプションにつけると、太さで表現することができます。
boxplot(y ~ z, data = df01, varwidth = TRUE)
4.3 異常値になる範囲を変更する
rangeで異常値とする外れ値を調整できます。0を指定すると、すべて正常値として取り扱います。
boxplot(y ~ z, data = df01, range=0.5)
4.3.1 異常値になる範囲（outlier）を注に入れる
異常値が無い場合は不要ですが、異常値は任意に指定できるrangeに依存するのでプロットする場合は、異常値の範囲を明示する必要があります。
例えば、以下のように書いておきましょう。
mtext(substitute(paste("Note)",
phantom(0),
"outlier"==paste(group("{", paste(x, "|", x>=Q3+range%.%IQR), "}"),
"∪",
group("{", paste(x, "|", x<Q1-range%.%IQR), "}"))), list(range = 1.5)),
1, 3, adj=0)5 箱ひげ図の見栄えを調整する
デフォルトでは素朴ですが、見栄えの自由度は意外に広いです。
5.1 縦ならび
業界慣習の縦横にあわせる事もできます。
boxplot(y ~ z, data = df01, horizontal = TRUE)
5.2 箱の太さと、線の長さを変える
箱の部分はboxwexを指定することで自動設定（0.8）よりも太くしたり細くしたりできます。
boxplot(y ~ z, data = df01, boxwex = 0.5)
ヒゲの先の線の長さは、staplewexで指定できます。boxwexに対する相対サイズです。
boxplot(y ~ z, data = df01, staplewex = 1.0)
5.3 事細かに表示を調整する
pars引数にリストで事細かに表示方法を指定できます。デフォルトのままでよいパラメーターは、リストの要素を消してしまえばよいです。グループごとに表示設定を変えたい場合は、グループの数と同じ長さのペクトル（もしくはグループの数の公約数と同じ長さのベクトル）をリストの値とすればよいです。
params <- list(
# 四分位の表示の設定
boxwex = 0.5,
boxlty = 1,
boxlwd = 1,
boxcol = "darkred",
boxfill = c("red", "pink", "purple", "yellow"),
# 中央値の表示の設定
medlty = NULL,
medlwd = NULL,
medpch = NULL,
medcex = NULL,
medcol = NULL,
medbg = NULL,
# ヒゲの部分（第1四分位未満/第2四分位以上から異常値にならない最小/最大値まで）の線の表示設定
whisklty = 2,
whisklwd = 1,
whiskcol = "blue",
# 最大/最小値の表示の設定
staplewex = c(0.5, 1, 0.8, 0.5),
staplelty = 1,
staplelwd = 2,
staplecol = "darkgreen",
# 異常値の表示の設定
outwex = NULL,
outlty = NULL,
outlwd = NULL,
outpch = 21,
outcex = 2,
outcol = "black",
outbg = "white"
)
boxplot(y ~ z, data = df01, pars = params)
5.4 グリッド線を描く
boxplotはpanel.firstが使えないので、グリッド線が箱ひげ図の上に描かれないように、低水準描画関数でプロット領域をつくって軸とグリッド線を描いた後に、add=TRUEをつけてboxplotを低水準描画関数として呼び出します。
plot.new()
xlim <- c(0.5, 0.5 + length(levels(df01$z)))
round_lim <- function(y, unit){
if(length(y)!=2) stop("Length of y must be two.")
y <- sort(y)
larger <- c(FALSE, TRUE)
sign(y)*ifelse(larger,
ifelse(y < 0, floor(abs(y) / unit), ceiling(abs(y) / unit)),
ifelse(y < 0, ceiling(abs(y) / unit), floor(abs(y) / unit))
)*unit
}
y_interval <- 3
ylim <- round_lim(range(df01$y), y_interval)
plot.window(xlim=xlim, ylim=ylim)
par(mar = c(4, 3, 0, 0))
axis(1, at = ceiling(xlim[1]):floor(xlim[2]),
labels = levels(df01$z), tick = FALSE)
axis(2, at=seq(ylim[1], ylim[2], y_interval), las=1, pos=xlim[1])
# grid関数だと目盛りにグリッド線をあわせられない為ablineで描く
abline(h = seq(ylim[1], ylim[2], y_interval), lty = 3, col="gray")
title(xlab = expression(z), ylab = expression(y))
boxplot(y ~ z, data = df01, add = TRUE, axes=FALSE)
6 まとめ
あまり箱ひげ図に親しみはなかったのですが、サブサンプルの基本統計量の確認作業に重宝しそうです。サブサンプル間の差異を検定したり、評価したりする方法ではない事を意識しながら、便利に使っていきましょう。boxplotで差がはっきり見えたら、統計的検定をしても、ベイズファクターを見ても、だいたい差異があると言えるような結果になると思いますが。