🎻

Rの組み込み関数でバイオリンプロットを書く

2022/12/24に公開

はじめに

R 言語においてバイオリンプロット（violin plot）を描くとき， ggplot2::geom_violin() を使うのが簡単です。紹介している記事もたくさん見つかります。全然それで構わないのですが，最近，個人的に「できる限り組み込み関数でなんとかする」ということに挑戦していて，かつバイオリンプロットを描きたいと思うことがありました。案外簡単にできたので，備忘録も兼ねて記事にすることにしました。

本記事は R Advent Calendar 2022 16 日目の記事になります（空いていたので後から埋めました）。

完成図

本記事のコードを使用すると，以下のようなバイオリンプロットを作成できるようになります。

順を追って作図

バイオリンプロットはデータのカーネル密度グラフを 90 度回転させ^[1]左右対称に描いたものです。R ではdensity()関数でカーネル密度推定ができるので，まずは，その返り値を素直に plot してみます。

カーネル密度のプロット

a <- rnorm(1000)
plot(density(a))

横軸が（カーネル密度の）データ点，縦軸がカーネル密度となるようなプロットが生成されました。当然ですが，バイオリンプロットの向きとは異なります。

カーネル密度の図を 90 度右に回転させる

図の回転^[2]に向けて，density()の返り値の構造を確認してみます。

str(density(a))
# List of 7
#  $ x        : num [1:512] -3.54 -3.53 -3.51 -3.5 -3.49 ...
#  $ y        : num [1:512] 4.09e-05 4.95e-05 5.96e-05 7.14e-05 8.53e-05 ...
#  $ bw       : num 0.223
#  $ n        : int 1000
#  $ call     : language density.default(x = a)
#  $ data.name: chr "a"
#  $ has.na   : logi FALSE
#  - attr(*, "class")= chr "density"

どうやら，7 つの要素からなるリストになっていて，$xの数値ベクトルが上記のプロットの x 軸，$yの数値ベクトルが上記の y 軸に対応しているようです。実際，これらをplot関数に具体的に指定すると，同じグラフが描画されます。

dens <- density(a)
plot(dens$x, dens$y, type = "l")

ということは，この x と y を入れ替えてplotに渡せば横軸がカーネル密度，縦軸がデータ点となる図が作成されます。

plot(dens$y, dens$x, type = "l")

左右対称なグラフにする

x 軸と y 軸を入れ替えたグラフは作成できたので，カーネル密度の値を負に（=左右）反転したグラフを一緒に描画すれば，バイオリンプロットが出来上がります。

plot(dens$y, dens$x, type = "l", xlim = c(-0.4, 0.4))
lines(-dens$y, dens$x)

plot()を実行する時点で，描画範囲が決まってしまうので，xlimを指定しておくのを忘れないようにしましょう。

2 つのグループのバイオリンプロットを描く

バイオリンプロットは，複数のグループのデータの分布をまとめて描きたいときに利用されると思います。グループが一つしかないのであれば，通常のカーネル密度グラフでいいですからね。

データ用の変数をもう一つ作成して，横に並べてみます。とりあえず，横方向に 1 だけずらした位置に描画してみます。2 つ目のバイオリンプロットを並べるときは，始めから lines() を使います（plot()すると新しい図が作成されてしまう^[3]）。

b <- rnorm(1000, 1, 0.5)
dens_b <- density(b)
# 以下の最初の2行は xlim 以外同じ
plot(dens$y, dens$x, type = "l", xlim = c(-0.4, 1.4))
lines(-dens$y, dens$x)
lines(1 - dens_b$y, dens_b$x)
lines(1 + dens_b$y, dens_b$x)

新しく追加したバイオリンプロットが元のものに重なっていますし，右側ははみ出してしまっています。バイオリンプロットにおいて横方向に表現されているカーネル密度は，分布感さえ伝わればよく，実際の値である必要はない（はず）です。そこで，カーネル密度の値$y に 1 未満の値をかけて幅を調整し見た目を整えます。

wd <- 0.5 # 補正値
plot(dens$y * wd, dens$x, type = "l", xlim = c(-0.4, 1.4))
lines(-dens$y * wd, dens$x)
lines(1 - dens_b$y * wd, dens_b$x)
lines(1 + dens_b$y * wd, dens_b$x)

良さげなバイオリンプロットになりました。

複数のグループを描画できるようにしつつ関数にまとめる

以上の処理を関数にまとめます。任意のグループ数のバイオリンプロットをかけるようにしています。

violinplot <- \(..., wd = 0.1, adj = 1) {
  dat <- list(...)
  n_dat <- length(dat)

  # 各グループのデータに対して density() を実行する
  l_dens <- lapply(dat, density, adjust = adj)
  # すべてのカーネル密度のデータ点の最大と最小を算出する
  # これが plot の y軸の範囲になる
  yl <- range(sapply(l_dens, \(dens) dens$x))

  # 土台となる描画範囲を作成する
  # xaxt = "n" は x軸のテキストを消している
  plot(1, type = "n", xlab = "", ylab = "", xaxt = "n",
       xlim = c(0.5, n_dat + 0.5), ylim = yl)
  # x軸をずらしながら，各グループのバイオリンプロットを描画
  for (i in seq_len(n_dat)) {
    dens <- l_dens[[i]]
    lines(i - dens$y * wd, dens$x)
    lines(i + dens$y * wd, dens$x)
  }
}

violinplot() 関数は以下のように使います。...には数値ベクトルの変数をカンマ区切りで好きなだけ指定することができます。軸のタイトルと，もはや無意味な x 軸の目盛はviolinplot()内で描画しないようにしているので，関数の外にそれぞれを設定する処理を追加しておきます。

violinplot(a, b, wd = 0.5, adj = 2)
axis(1, at = c(1, 2), labels = c("A", "B"))
title(ylab = "value")

実行結果はこのようになります。

いい感じですね。

色を塗りたい場合

個人的には，上の質素な感じでも実用に耐えると思ったのですが，バイオリンプロットの中を塗りたいという需要もありそうなので，そのための関数も作成してみました。

ほとんどの処理や発想はこれまでと同様ですが，塗りつぶすために，バイオリンプロットを図形（polygon）として描画しています。polygon() 関数の仕様に合わせるために，for 文の中がすこし変わっていてトリッキーな感じになっています。この説明は省略します。ややこしくなりそうですし，本記事のメインの内容からは逸れるからです。

violinplot_filled <- \(..., wd = 0.1, adj = 2, cols = NULL) {
  dat <- list(...)
  n_dat <- length(dat)
  # 指定された色がグループ数より少ないとき，足りない分だけRの色空間（？）からランダムに追加
  if (length(cols) < n_dat) {
    cols <- c(cols, sample(colors(), n_dat - length(cols)))
  }

  l_dens <- lapply(dat, density, adjust = adj)
  yl <- range(sapply(l_dens, \(dens) dens$x))

  plot(1, type = "n", xlab = "", ylab = "", xaxt = "n",
       xlim = c(0.5, n_dat + 0.5), ylim = yl)
  for (i in seq_len(n_dat)) {
    dens <- l_dens[[i]]
    polyx <- c(i - dens$y * wd, rev(i + dens$y * wd))
    polyy <- c(dens$x, rev(dens$x))
    polygon(polyx, polyy, col = cols[i])
  }
}

この関数を以下のように実行した結果が，本記事の最初に載せた図になります。

violinplot_filled(a, b, wd = 0.5, cols = palette.colors(2, "Classic Tableau"))
axis(1, at = c(1, 2), labels = c("A", "B"))
title(ylab = "value")

おわりに

はじめに述べたように，最近，筆者は組み込み関数縛りで R のコードを書いていて，バイオリンプロット以外の作図もいわゆる base plot でしているのですが，割とシンプルに作図できると感じています。最後に紹介した関数のせいでシンプルな印象がなくなっていそうですが，それまでの描画だけに限ればシンプルだったのではないでしょうか。R での作図といえばggplotがスタンダードな感じになっています。この現状はそれだけggplotが使いやすいからこそだと思うのですが^[4]，base plot ももう少し注目されて新しい記事が出たりしてもいいのかなと思いました^[5]。

最後までお読みいただきありがとうございました。いいね・サポートをいただけると大変励みになりますので，ぜひそちらもよろしくお願いします。

脚注

「90 度回転」は適切ではありません。直感的にわかりやすい気がしたので，この表現を用いました。後にも出てくるように，「x 軸と y 軸を入れ替える」が正しいです。 ↩︎
繰り返しになりますが，「回転」というのは適切ではありません。図を回転させると異なった図になります。 ↩︎
もちろん，plot(..., add = TRUE) で既存の図に重ねて描画することはできますが，軸を整えたりするのが面倒です。 ↩︎
かくいう，私自身もこれまでずっと ggplot で作図してきました。 ↩︎
base plot は ggplot よりも前から存在しているから紹介記事や書籍がすでに飽和している？？？ ↩︎