ggplot2でヒストグラムに経験分布を重ねて描く

2ヶ月に1回くらい同じようなことをやろうとして都度調べている気がするのでまとめておく。

はじめに

そもそもヒストグラムに経験分布を重ねて描く必要はないし、おそらく描かない方が良い。ヒストグラムと経験分布のようにスケールの異なるプロットを1つにまとめてしまうと、軸と値の対応が分かりづらくなる。ggplot2でも異なるスケールの2つのプロットを1つにまとめて表示することは推奨されていないので、そういうことをやろうとすると非直感的な操作が必要となる。

グラフは分けて描いて、必要であれば並べれば良い。そうしたことを行う簡単な方法はいくつかあるが、patchworkパッケージ（The Composer of Plots • patchwork）を使用した例を示す。

library(ggplot2)
library(patchwork)

g <- ggplot(iris, aes(x = Sepal.Length, fill = Species, color = Species))
p1_hist <- g + geom_histogram(position = "identity", alpha = .5) +
  theme(
    axis.text.x = element_blank(), 
    axis.ticks.x = element_blank(), 
    axis.title.x = element_blank()
  ) +
  labs(title = "Histogram")
p2_ecdf <- g + stat_ecdf() + labs(title = "ecdf",y = "F(x)")

p1_hist / p2_ecdf

装飾等にこだわらなければ正味4〜5行程度の記述でプロットが作成できるし、複雑な手順もなく覚えやすい。また、プロット中の値について誤解される恐れもない。

したがって以下の方法は表現としてはあまり推奨できない。ただ、紙面の都合で重ねたい場合もあるだろうし、パレート図のように既に世の中に広まってしまっている複合グラフもある。そういったものを作成したい場合には役立つだろう。

ポイント

ヒストグラムに経験分布を重ねて描くための主なポイントは2つある。

経験分布の計算
第2軸の設定

1. 経験分布の計算

まず、経験分布を計算する必要がある。単に経験分布のみをプロットするのであれば、stat_ecdf()を使うという手がある。しかし、経験分布のy軸スケールは0〜1なので、大抵の場合ヒストグラムのレンジと合わない。

ggplot2ではプロットに使える座標はあくまで1種類なので、ヒストグラムに経験分布を重ねたければ、ヒストグラムのスケールに合うように経験分布の値を調整する必要がある。具体的には、stat_bin()のcomputed variablesを利用して累積和の計算と高さの調整を行うといった操作を行う。*1

2. 第2軸の設定

次に第2軸の設定が必要となる。第2軸はscale_y_continuous()等のsec.axis引数に詳細を指定することで表示できる。詳細の指定はsec_axis()関数などを用いて組み立てる。sec_axis()には第1軸の値を変換するための関数を指定できるので、これにより値を変換して第2軸のラベルを作成する。*2

今回は経験分布のスケールを作りたいので、0〜データの最大値を0〜1に変換した第2軸が表示できればよい。そのためには第1軸の値をデータの最大値で割れば良いので、「データの最大値」をいかに調達するかが課題となる。

方法1. countの最大値を保持しておく

以前Qiitaに同じような内容を書いた。

qiita.com

書き方が古かったりするので、データセットをirisに変えて再度書いた。

library(ggplot2)
bins <- 30

ggplot(iris, aes(x = Sepal.Length)) +
  geom_histogram(bins = bins) +
  stat_bin(
    aes(y = 
          after_stat(
            cumsum(count / sum(count) * (m_cnt <<- max(count)))
          )
    ), 
    geom = "line", bins = bins
  ) +
  scale_y_continuous(
    name = "count",
    sec.axis = sec_axis(
      trans = ~ . / m_cnt,
      name = "ecdf"
    )
  )

ここでは、after_stat()を使ってstat_bin()のcumputed variablesであるcountから累積和を計算し、さらに合計値と最大値を使ってスケールを調整している。同時にcountの最大値を覚えておいて、右側y軸の範囲を0〜1に調整するのに使用している。

この方法はグローバル環境にm_cntという変数が生じてしまうという欠点がある。

また、グループに分けるようなプロットではうまくいかない。例えば以下のようにSpeciesごとにヒストグラムを作成してみる。

ggplot(iris, aes(x = Sepal.Length, fill = Species, color = Species)) +
  geom_histogram(bins = bins, alpha = .5, position = "identity") +
  stat_bin(
    aes(y = 
          after_stat(
            cumsum(count / sum(count) * (m_cnt <<- max(count)))
          )
    ), 
    geom = "line", bins = bins
  ) +
  scale_y_continuous(
    name = "count",
    sec.axis = sec_axis(
      trans = ~ . / m_cnt,
      name = "ecdf"
    )
  )

すると、以下のような結果になってしまう。

コードと結果は省略するが、facet_gridなどでグラフを別々にしてみても同様である。これはafter_stat()内の計算がグループ別に分かれていないことによる。

方法2. 累積値をグループごとに計算する + データの最大値を第1軸の値から計算する

要は累積値の計算がグループごとに行われないのが問題なので、グループごとに計算してやればよい。また、データの最大値に対応する値は、実は第1軸を変換する最に用いる値から計算できる。

まずコードと結果を示す。この方法だと事前の集計等は不要で、グローバル環境に変数も出てこず、ggplot()で完結する。

ggplot(iris, aes(x = Sepal.Length, fill = Species, color = Species)) +
  geom_histogram(bins = bins, position = "identity", alpha = .5) +
  stat_bin(
    aes(
      y = after_stat(
        ave(
          count, 
          group,
          FUN = function(.){cumsum(.) / sum(.) * max(count)} 
          )
      )
    ),
    geom = "line",
    position = "identity",
    bins = bins
  ) +
  scale_y_continuous(
    name = "count",
    sec.axis = sec_axis(
      name = "ecdf",
      trans = ~ ./ (max(.) / 1.05),
      labels = scales::percent
    )
  )

ポイントはstat_bin()で使用しているave()関数と、sec_axis()のtransの変換式に指定している./max(.) * 1.05という式にある。

`ave()`によるグループ別の累積値の計算

ave()はRが標準で読み込むstatsパッケージに含まれている関数で、グループごとに何らかの関数を適用できる。戻り値は第一引数の数値ベクトルと同じ長さのベクトルとなる（SQLで言うところのWindow関数と同様の挙動をする）ので、ggplot2などとは組み合わせやすい。

R: Group Averages Over Level Combinations of Factors

ave()をstat_bin()内部でafter_stat()と組み合わせて使用すると、グループごとに累積和を計算するといったようなことができる。該当箇所を再掲する。

after_stat(
        ave(
          count, 
          group,
          FUN = function(.){cumsum(.) / sum(.) * max(count)} 
          )
      )

groupは明示的にaes()の中で指定していないが、after_stat()内で利用できる。中身はfillとcolorの指定に対応している。もしfillやcolorを指定しておらずグループが存在していなくともgroup変数は利用可能なので、この書き方をしておけば問題が起こりにくい。

FUN =には各グループに適用する関数を指定する必要がある。なお、引数名は省略できない（グループ変数の指定が任意個可能なため）。ここでは無名関数をその場で定義して与えているが、別途定義したものを与えても良い。

`. / (max(.) / 1.05)`による第2軸の値の計算

最初の方法では、countの最大値を覚えておいて、その値で第1軸の値を割ることで「0〜countの最大値」から「0〜1」への変換を行った。ここでは「countの最大値」の代わりに「max(.) / 1.05」を用いている。

まず、transに指定する式が実際には何をしているかを説明する。transに指定する式(formula)は第1軸の値を変換するものであるが、実際には第1軸のレンジを1000等分したベクトルが式に基づいて定義された関数に渡され、第2軸のラベルに用いる値が計算されている。そして、「第1軸のレンジ」は連続値の場合、デフォルトでデータの最大、最小から5%だけ拡張された範囲となっている。なお、これはexpand引数で制御でき、ヘルプにも記載がある（Position scales for continuous data (x & y) — scale_continuous • ggplot2）。

つまり、max(.)は描画範囲の上端の値に対応し、それはデフォルトだとcountの最大値の1.05倍である。つまりmax(.) / 1.05がcountの最大値に対応する。この値で割れば、「0〜countの最大値」を「0〜1」に変換することができる。

他

ave()を使う方法は色々と応用ができる。

例えばstat_density()と組み合わせれば、密度関数の積分値を図示できる。

ggplot(iris, aes(x = Sepal.Length, fill = Species, color = Species)) +
  geom_density(alpha = .5) +
  stat_density(
    aes(
      y = after_stat(
        ave(
          density,
          group,
          FUN = function(.) {cumsum(.) / sum(.) * max(density)}
        )
      )
    ),
    position = "identity",
    geom = "line"
  ) +
  scale_y_continuous(
    name = "density",
    sec.axis = sec_axis(
      trans = ~ ./max(.) * 1.05,
      name = "area under curve"
    )
  )

また、ecdfは不要で単に累積カウントを比較したいという場合にも役立つ。

ggplot(diamonds, aes(x = price, color = cut)) +
  stat_count(
    aes(y = after_stat(ave(count, group, FUN = cumsum))),
    position = "identity",
    geom = "line"
  ) + theme_minimal() +
  labs(y = "Cumulative count")

参考

ave()を使う方法は以下で見つけた。

r - ggplot faceted cumulative histogram - Stack Overflow

*1:stat_ecdf()だと（当然ではあるが）computed variablesからヒストグラムの高さを計算できないので、stat_bin()を使わなければならない。なお、利用可能なcomputed variablesは関数のヘルプを見ると確認できる。

*2:第2軸が設定できるといっても、これは複数の異なるスケールのグラフを重ねて描くのに適した機能ではない。第2軸は第1軸の単純な複製（そのためにdup_axis()という関数も用意されている）や、別の表現（摂氏と華氏、異なるタイムゾーン下での時刻など）を示すために用いるのがおそらく本来の用法で、そのような操作であれば容易にできる。