はじめに
素直なやり方ではできなかったのでむりくりなんとかする黒魔術シリーズです。もっといい方法があれば教えてください。
この記事では、R言語で三角図の等高線図とヒートマップを作成します。
【前の内容】
【他の内容】
【目次】
ggplot2で三角グラフの等高線を作図したい
前回は、三角図や三角ダイアグラム(ternary diagram)などと呼ばれるグラフをggplot2パッケージを利用して作成した。今回は、三角図上に等高線図(contour map)とヒートマップ(heatmap)を描画したい。しかし、ggplot2で等高線を作図するには格子状の点(直交する点)を使う必要があり、三角座標に変換すると格子点にならない。そこで、三角座標を含めた2次元座標上の格子点を作成して、元の3次元座標に戻して目的の計算をすることで、三角図上に等高線を描画する。
三角図の座標の作図については前回の記事を参照のこと。
次のパッケージを利用する。
# 利用パッケージ library(tidyverse) library(gganimate) library(MCMCpack)
基本的にパッケージ名::関数名()の記法を使うので、パッケージを読み込む必要はない。ただし、作図コードがごちゃごちゃしないようにパッケージ名を省略しているので、ggplot2を読み込む必要がある。
また、magrittrパッケージのパイプ演算子%>%ではなく、ネイティブパイプ演算子|>を使っている。%>%に置き換えても処理できるが、その場合はmagrittrを読み込む必要がある。
三角座標からの再変換式の確認
まずは、三角座標上の点$\mathbf{y}$から3次元座標上の点$\mathbf{x}$に変換する計算式を確認する。
総和が1の3次元の変数(点)
に対して、次の式で2次元の変数(点)$\mathbf{y}$に変換する。
$\mathbf{y}$は、正三角形の座標上の点になるのであった。
ここでは、$\mathbf{y}$から$\mathbf{x}$に再変換する式を導出する。
途中式(クリックで展開)
式(1)を$x_3$について整理する。
この式を式(2)に代入して、$x_2$について整理する。
$y_1, y_2$から$x_2$を求められる式が得られた。
続いて、式(3)を式(1')に代入して、$x_3$について整理する。
$y_2$から$x_3$を求められる式が得られた。
さらに、$\mathbf{x}$の条件式を$x_1$について整理する。
$x_2, x_3$から$x_1$を求められる式が得られた。
式(3),(4),(5)を用いて、$\mathbf{x}$が得られる。
格子点$\mathbf{y}$を用意して、$\mathbf{x}$に変換した上で目的の計算を行うことで等高線を作図できる。
等高線図
三角座標上の等高線図を作成する。例として、ディリクレ分布の確率密度の等高線図を作成する。
$y_1$と$y_2$の値を作成する。
# 三角座標の値を作成 y_1_vals <- seq(from = 0, to = 1, length.out = 301) y_2_vals <- seq(from = 0, to = 0.5*sqrt(3), length.out = 300) head(y_1_vals); head(y_2_vals)
## [1] 0.000000000 0.003333333 0.006666667 0.010000000 0.013333333 0.016666667 ## [1] 0.000000000 0.002896406 0.005792812 0.008689218 0.011585624 0.014482030
$0 \leq y_1 \leq 1$の値をy_1_vals、$0 \leq y_2 \leq \frac{\sqrt{3}}{2}$の値をy_2_valsとして作成する。グラフが粗い場合や処理が重い場合は、y_1_valsの間隔(by引数)や要素数(length.out引数)を調整する。
$\mathbf{y}$の値を作成する。
# 格子点を作成 y_mat <- tidyr::expand_grid( y_1 = y_1_vals, y_2 = y_2_vals ) |> # 格子点を作成 as.matrix() # マトリクスに変換 head(y_mat)
## y_1 y_2 ## [1,] 0 0.000000000 ## [2,] 0 0.002896406 ## [3,] 0 0.005792812 ## [4,] 0 0.008689218 ## [5,] 0 0.011585624 ## [6,] 0 0.014482030
y_1_valsとy_2_valsの要素の全ての組み合わせ(格子状の点)をexpand_grid()で作成して、マトリクスに変換する。y_matはグラフ描画用の値で、各行が$\mathbf{y} = (y_1, y_2)$に対応する。
$\mathbf{x}$の値を作成する。
# 3次元変数に変換 x_mat <- tibble::tibble( x_2 = y_mat[, 1] - y_mat[, 2] / sqrt(3), x_3 = 2 * y_mat[, 2] / sqrt(3) ) |> # 元の座標に変換 dplyr::mutate( x_2 = dplyr::if_else(x_2 >= 0 & x_2 <= 1, true = x_2, false = as.numeric(NA)), x_3 = dplyr::if_else(x_3 >= 0 & x_3 <= 1 & !is.na(x_2), true = x_3, false = as.numeric(NA)), x_1 = 1 - x_2 - x_3, x_1 = dplyr::if_else(x_1 >= 0 & x_1 <= 1, true = x_1, false = as.numeric(NA)) ) |> # 範囲外の値をNAに置換 dplyr::select(x_1, x_2, x_3) |> # 順番を変更 as.matrix() # マトリクスに変換 head(x_mat)
## x_1 x_2 x_3 ## [1,] 1 0 0 ## [2,] NA NA NA ## [3,] NA NA NA ## [4,] NA NA NA ## [5,] NA NA NA ## [6,] NA NA NA
$\mathbf{y}$の値y_matから$\mathbf{x}$の値x_matに変換する。ただし、y_matは三角座標外の値を含み、その値を式(3-5)で計算すると総和が1でない点になる。範囲外の値(点)を取り除くのではなく、欠損値NAに置換する。
x_matの各行が$\mathbf{x} = (x_1, x_2, x_3)$に対応するように、マトリクスに変換する。
ディリクレ分布の確率密度を計算する。
# ディリクレ分布のパラメータを指定 alpha_k <- c(4, 2, 3) # ディリクレ分布の確率密度を計算 dens_df <- tibble::tibble( y_1 = y_mat[, 1], # x軸の値 y_2 = y_mat[, 2], # y軸の値 density = MCMCpack::ddirichlet(x = x_mat, alpha = alpha_k), # 確率密度 ) |> dplyr::mutate( fill_flg = !is.na(rowSums(x_mat)), density = dplyr::if_else(fill_flg, true = density, false = as.numeric(NA)) ) # 範囲外の値をNAに置換 dens_df
## # A tibble: 90,300 × 4 ## y_1 y_2 density fill_flg ## <dbl> <dbl> <dbl> <lgl> ## 1 0 0 0 TRUE ## 2 0 0.00290 NA FALSE ## 3 0 0.00579 NA FALSE ## 4 0 0.00869 NA FALSE ## 5 0 0.0116 NA FALSE ## 6 0 0.0145 NA FALSE ## 7 0 0.0174 NA FALSE ## 8 0 0.0203 NA FALSE ## 9 0 0.0232 NA FALSE ## 10 0 0.0261 NA FALSE ## # … with 90,290 more rows
グラフ用の値y_matと、計算用の値x_matにより求めた確率密度をデータフレームに格納する。
ディリクレ分布の確率密度は、MCMCpackパッケージのddirichlet()で計算できる。確率変数の引数xにx_mat、パラメータの引数alphaに設定したパラメータalpha_kを指定する。
パラメータの値を数式で表示するための文字列を作成する。
# パラメータラベル用の文字列を作成 param_text <- paste0( "list(", "alpha==(list(", paste0(alpha_k, collapse = ", "), "))", ", x==(list(x[1], x[2], x[3]))", ")" ) param_text
## [1] "list(alpha==(list(4, 2, 3)), x==(list(x[1], x[2], x[3])))"
ギリシャ文字などの記号を使った数式を表示する場合は、expression()の記法を使う。等号は"=="、複数の(数式上の)変数を並べるには"list(変数1, 変数2)"とする。(プログラム上の)変数の値を使う場合は、parse()のtext引数に指定する。
三角図上に確率密度の等高線を描画する。
# 等高線図を作成 ggplot() + geom_segment(data = ternary_grid_df, mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end), color = "gray50", linetype = "dashed") + # 三角図のグリッド線 geom_segment(data = ternary_axis_df, mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end), color = "gray50") + # 三角図の枠線 geom_text(data = ternary_ticklabel_df, mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v, angle = angle)) + # 三角図の軸目盛ラベル geom_text(data = ternary_axislabel_df, mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v), parse = TRUE, size = 6) + # 三角図の軸ラベル geom_contour_filled(data = dens_df, mapping = aes(x = y_1, y = y_2, z = density, fill = ..level..), alpha = 0.8) + # 確率密度の等高線 scale_x_continuous(breaks = c(0, 0.5, 1), labels = NULL) + # x軸 scale_y_continuous(breaks = c(0, 0.25*sqrt(3), 0.5*sqrt(3)), labels = NULL) + # y軸 coord_fixed(ratio = 1, clip = "off") + # アスペクト比 theme(axis.ticks = element_blank(), panel.grid.minor = element_blank()) + # 図の体裁 labs(title = "Ternary Contour Plot", subtitle = parse(text = param_text), fill = "density", x = "", y = "")
geom_contour()またはgeom_contour_filled()で等高線図を描画する。三角図外の点は欠損値NAなので描画されない(つまりデータの半分を捨てる非効率設計である)。
ヒートマップ
ついでに、ヒートマップを作成する。
# ヒートマップを作成 ggplot() + geom_segment(data = ternary_grid_df, mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end), color = "gray50", linetype = "dashed") + # 三角図のグリッド線 geom_segment(data = ternary_axis_df, mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end), color = "gray50") + # 三角図の枠線 geom_text(data = ternary_ticklabel_df, mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v, angle = angle)) + # 三角図の軸目盛ラベル geom_text(data = ternary_axislabel_df, mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v), parse = TRUE, size = 6) + # 三角図の軸ラベル geom_tile(data = dens_df, mapping = aes(x = y_1, y = y_2, fill = density, alpha = fill_flg)) + # 確率密度のヒートマップ scale_alpha_manual(breaks = c(TRUE, FALSE), values = c(0.8, 0), guide = "none") + # 透過 scale_fill_viridis_c() + # グラデーション #scale_fill_gradientn(colors = c("blue", "green", "yellow", "orange")) + # グラデーション scale_x_continuous(breaks = c(0, 0.5, 1), labels = NULL) + # x軸 scale_y_continuous(breaks = c(0, 0.25*sqrt(3), 0.5*sqrt(3)), labels = NULL) + # y軸 coord_fixed(ratio = 1, clip = "off") + # アスペクト比 theme(axis.ticks = element_blank(), panel.grid.minor = element_blank()) + # 図の体裁 labs(title = "Ternary Heatmap", subtitle = parse(text = param_text), fill = "density", x = "", y = "")
geom_tile()でヒートマップを描画する。alpha引数にfill_flg列を指定し、scale_alpha_manual()のbreaks引数にTRUE, FALSE(fill_flg列の値)、values引数に透過率を指定する。FALSEの場合の透過率を0にすると、三角図の範囲外の要素(density列が欠損値NAの要素)は描画されない。
グラデーションの色はscale_fill_viridis_c()やscale_fill_gradientn()などで設定できる。
同様に、散布図を使ってヒートマップを描画する。
# 散布図によるヒートマップを作成 ggplot() + geom_point(data = dens_df, mapping = aes(x = y_1, y = y_2, color = density, alpha = fill_flg)) + # 確率密度の散布図 scale_alpha_manual(breaks = c(TRUE, FALSE), values = c(0.5, 0), guide = "none") + # 透過 scale_color_viridis_c() + # グラデーション #scale_color_gradientn(colors = c("blue", "green", "yellow", "orange")) + # グラデーション geom_segment(data = ternary_grid_df, mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end), color = "gray50", alpha = 0.5, linetype = "dashed") + # 三角図のグリッド線 geom_segment(data = ternary_axis_df, mapping = aes(x = y_1_start, y = y_2_start, xend = y_1_end, yend = y_2_end), color = "gray50") + # 三角図の枠線 geom_text(data = ternary_ticklabel_df, mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v, angle = angle)) + # 三角図の軸目盛ラベル geom_text(data = ternary_axislabel_df, mapping = aes(x = y_1, y = y_2, label = label, hjust = h, vjust = v), parse = TRUE, size = 6) + # 三角図の軸ラベル scale_x_continuous(breaks = c(0, 0.5, 1), labels = NULL) + # x軸 scale_y_continuous(breaks = c(0, 0.25*sqrt(3), 0.5*sqrt(3)), labels = NULL) + # y軸 coord_fixed(ratio = 1, clip = "off") + # アスペクト比 theme(axis.ticks = element_blank(), panel.grid.minor = element_blank()) + # 図の体裁 labs(title = "Ternary Heatmap", subtitle = parse(text = param_text), fill = "density", x = "", y = "")
color引数を使って同様に処理する。
(グリッド線がうまく透けなかったので上から表示した。点の数を減らせば透けると思う。)
以上で、私が欲しいグラフが得られた。等高線図のアニメーションについては「ディリクレ分布の作図」を(このあと書くので)参照のこと。
おわりに
なんか数式がうじゃうじゃ出てきた。そんなつもりなかった。せめてもの気持ちでクリックしないと途中式が見えないようにしておきましたので許してください。
