箱ひげ図複数カテゴリ描写問題

また嵌ったので共有したいんや。

起こったこと

ggplot2で箱ひげ図＋αで色々書こうと思ったらちょっと工夫が必要だった。
でもそこまで説明するのに長くなりそうなので、
まずは箱ひげ図問題について書く。

ある連続変数の分布を、異なるカテゴリごとに可視化するのにオーソドックスな方法の一つが箱ひげ図である。
例えば異なる地域5つの地域の小学1年生の身長の分布、みたいなのを見たいとき。
でこの箱ひげ図ってのは実はいくつか問題があって、

よく経験されるのは
1. いくつかの表示の意味が誤解されている。
2. サンプル数についての情報を与えない
3. 分布が単峰性であるかのような誤解を与えやすい。
だね。
個人的にはこれが堅牢な可視化とは思っていない。

表示の意味について

require(tidyverse)
data(diamonds)
sd <- 114514

ggsave2 <- function(plot, wid=9, hei=9){
  plot_name <- deparse(substitute(plot))
  file_name <- paste(plot_name, ".png", sep = "",collapse = "")
  ggsave(filename = file_name,plot = plot,device = "png",width = wid, height = hei,dpi = 300,units = "cm")
}



dmd <- diamonds %>%
  .[sample(x = 1:nrow(.), size = 500, replace = F),]

ってしとくよ。
ggsave2はggplot2のプロットオブジェクトをそのマッマ保存する関数だ。
diamondはそのままだとちょっと多いので減らした。

で問題の箱ひげ図だけど、こうする。

pbox1 <- ggplot()+theme_light()+
  geom_boxplot(data = dmd,
               aes(x=cut, y=carat))
ggsave2(pbox1)

図のあれこれの意味がわかるだろうか。
真ん中の箱が第1四分位、第2四分位（中央値）、第3四分位なのはまあいいよね。
問題はその上下の縦線と、その先にある点だ。

すぐに答えられた君は偉い。
ワイは詰まった。
でマニュアルは以下だ。
https://cloud.r-project.org/web/packages/ggplot2/index.html

boxplotのマニュアルにはこんな感じ。

The lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles).
This differs slightly from the method used by the boxplot() function, and may be apparent with
small samples. See boxplot.stats() for for more information on how hinge positions are calculated
for boxplot().
The upper whisker extends from the hinge to the largest value no further than 1.5 * IQR from the
hinge (where IQR is the inter-quartile range, or distance between the first and third quartiles). The
lower whisker extends from the hinge to the smallest value at most 1.5 * IQR of the hinge. Data
beyond the end of the whiskers are called "outlying" points and are plotted individually.
In a notched box plot, the notches extend 1.58 * IQR / sqrt(n). This gives a roughly 95% confidence
interval for comparing medians. See McGill et al. (1978) for more details.

縦線はヒンジから1.5*IQR以内の最大数まで伸びてるってことらしいね。
これを超えたものがoutlyingとしてドットになっている。
outlyingがない場合には、縦線の終着点が最大値、最小値ということだろうね。

これはggplot2の方言なのかな？
ともかくこれは誤解しやすい。

サンプル数情報

さっきの図で、君はサンプル数をどのように捉えたかな？

dmd_smr <- dmd %>%
  group_by(cut) %>%
  summarise(n=n())
dmd_smr %>%
  DT::datatable()

ってするとこんな感じ。

意外!それは不均一！

ということで箱ひげ図はサンプル数についての誤解を与えやすい。

変な誤解を与えて仕事を増やしたくないので、ともかくサンプル数をオーディエンスに与えたい。
かんたんには上に棒グラフを積むことだろうか。

p_bar <- ggplot()+theme_light()+
  geom_bar(data = dmd_smr,
               aes(x=cut, y=n),
           stat = "identity")
ggsave2(p_bar)

ただ、いちいちFigをあわせてレイアウトを調整するのって結構面倒くさい。

分布の形態について

下程極端なデータセットはあまりないかもしれないけど、ともかくこんなのでどうか。

v1_0 <- rnorm(n=100, mean = 0, sd = 0.5)
v1_2 <- rnorm(n=100, mean = 2, sd=0.5)
v2_1 <- rnorm(n=200, mean = 1, sd = 1.5)

df_demo <- data.frame(val = c(v1_0, v1_2, v2_1),
                      label = c(rep("twopeak", 200),
                                rep("onepeak", 200)))
p_box2 <- ggplot()+theme_light()+
  geom_boxplot(data = df_demo,
               aes(x=label, y=val))
ggsave2(p_box2)

これだけ見るとピークは一つみたいに見えてしまう。

p_hist <- ggplot()+theme_light()+
  geom_histogram(data = df_demo,
               aes(x=val, fill=label),
               position = "identity", color="black", alpha = 0.7, bins=30)+
  coord_cartesian(expand = F)
ggsave2(p_hist, wid=13)

ってやると

twopeak側はピークが2つあって、中央値とかで2群を比較するのはそもそもどうなのかという話になる。

でもこれは群が3つ以上になるとヒストグラムを別につけるのも面倒になってくる。

問題まとめ

というわけで

1. 複数カテゴリグループにおける、連続変数の分布を、なるべく文句の出ないように可視化したい
2. 仕事を増やしたくない
3. 分布についての統計値を示したい
4. 各群のサンプル数をある程度示したい
5. 仕事を増やしたくない
6. 分布の形について情報を与えたい（さもなくば平然とt検定しろとか言われる）
7. 仕事を増やしたくない
8. 仕事を増やしたくない

という要請を考えたとき、まあ色々方法はあると思うんだけど、
サンプル数が200以下ぐらいなら個人的には以下のようにしている。

解決策

boxplot + jitterplot

pboxjit <- ggplot()+theme_light()+
  geom_boxplot(data = dmd,
               aes(x=cut, y=carat), outlier.shape = NA)+
  geom_jitter(data = dmd,
               aes(x=cut, y=carat),
              color="red", size=0.7)

ggsave2(pboxjit)

今回はnがちょっと多いのでごちゃぁってしているけど、
一般的な臨床研究のサンプル数なんかはそこまでイカないことも多くって、
とりあえずこれで示すことが多い。
そうすると、t検定しろとかよくわからない指令が降りることは少なくて、
この分布を踏まえて、どのような次の検証が可能かという議論に入れることが多い、気がしている。

現場からは以上です。
ではなくて、もうちょっとだけ続くんじゃよ。