Rで解析：ggplot2の文法のまま3Dプロット「ggcube」パッケージ

本パッケージは、3D対応のgeomとstat、coord_3dコマンドの座標系でggplot2の文法のまま3Dプロットが可能です。coord_3d()はpitch・roll・yawの3軸まわりの回転、透視投影と平行投影の切り替え、光源の向きや強さを指定する照明効果を1つの関数で制御することが可能です。さらに、animate_3dコマンドでは、プロットが回転するGIFや動画を書き出すことが可能です。

パッケージバージョンは0.2.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

東京大学のデータサイエンティスト育成講座 ~Pythonで手を動かして学ぶデ―タ分析~ | 塚本邦尊, 山田典一, 大澤文孝, 中山浩太郎, 松尾 豊[協力]

Amazonで塚本邦尊, 山田典一, 大澤文孝, 中山浩太郎, 松尾 豊[協力]の東京大学のデータサイエンティスト育成講座 ~Pythonで手を動かして学ぶデ―タ分析~。アマゾンならポイント還元本が多数。

www.amazon.co.jp

パッケージのインストール

下記コマンドを実行してください。library(“ggcube”)を実行すると自動で「ggplot2」パッケージが読み込まれます。

# パッケージのインストール
install.packages("ggcube")
# パッケージの読み込み（依存のggplot2も同時に読み込まれる）
library("ggcube")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

ggcubeの関数は、大きく4つの役割に分かれています。座標系を決めるcoord_3d()、データを立体に変換して描くgeom_bar_3d()やgeom_smooth_3d()などの3D用geom、z軸の目盛りや見た目を整えるscale_z_continuous()やテーマ要素、そして完成した図を回転させるanimate_3d()です。この順に組み合わせていくのが基本的な流れになります。

3D図形の中に注釈用の点・文字・線分を埋め込む：annotate_3dコマンド

実行すると、指定した3D座標に注釈が配置された図が得られます。ggplot2のannotate()が別レイヤーになるのに対し、annotate_3dで作った注釈は本体の図形と一緒に奥行き順でソートされるため、面の裏側に回り込んだ注釈が図形で正しく隠されます。

### データ作成：中央が盛り上がった地形を模した格子データ #####
terrain <- expand.grid(x = seq(-3, 3, length.out = 40),
                       y = seq(-3, 3, length.out = 40))

# 中心からの距離に応じて高さを減衰させる
terrain$z <- 30 * exp(-(terrain$x^2 + terrain$y^2) / 4)
########

# ベースの3Dプロットを作成
p_annot <- ggplot(terrain, aes(x, y, z)) +
  coord_3d(light = light(mode = "hsl", direction = c(-1, 0, 1)))

# 山頂に赤い点を1つ置く
p_annot + geom_surface_3d(
  annotate = annotate_3d("point", x = 0, y = 0, z = 30,
                         color = "red", size = 4))

# ベクトルを渡すと、1回の呼び出しで複数の点を配置できる
p_annot + geom_surface_3d(
  annotate = annotate_3d("point", x = seq(-2, 2, 1), y = 0,
                         z = 30 * exp(-seq(-2, 2, 1)^2 / 4),
                         color = "green", size = 3))

# 点・文字・線分を組み合わせて山頂を指し示す
p_annot + geom_surface_3d(
  annotate = list(
    annotate_3d("point", x = 0, y = 0, z = 30,
                color = "red", size = 4),
    annotate_3d("text", x = 0, y = 0, z = 30,
                label = "Rは進化する", color = "red",
                vjust = -1.2, fontface = "bold"),
    annotate_3d("segment", x = 0, y = 0, z = 30,
                xend = 0, yend = 0, zend = 0,
                color = "red", linewidth = 1)))

回転・透視投影・縦横比を指定して3D座標系を作る：coord_3dコマンド

ggcubeで作図する際の中心となる関数です。既存のggplotオブジェクトに足すだけで立体表示になり、視点の回転、遠近感の強さ、立方体の縦横比、どの面を描くか、軸ラベルをどの辺に置くを、このコマンドで指定できます。

# ベースの3Dプロットを作成
p_coord <- ggplot(terrain, aes(x, y, z)) +
  geom_surface_3d(color = NA) +
  scale_fill_viridis_c() +
  theme(legend.position = "none")

# 既定の視点で描画する
p_coord + coord_3d()

# 回転をすべて0にすると、x-y平面を真上から見た状態になる
p_coord + coord_3d(pitch = 0, roll = 0, yaw = 0)

# 3軸の回転を組み合わせて任意の角度から見る
p_coord + coord_3d(pitch = 25, roll = -50, yaw = 40)

# 視点を近づけて遠近感を強調する
p_coord + coord_3d(dist = 1)

# 透視投影を切ると、平行な線が平行のまま描かれる
p_coord + coord_3d(persp = FALSE)

# データの比率を保ったまま描画する
p_coord + coord_3d(scales = "fixed")

# z軸だけを2倍の長さに引き伸ばす
p_coord + coord_3d(ratio = c(1, 1, 2))

# 立方体の6面すべてを描く（手前の面は既定で20％の不透明度）
p_coord + coord_3d(panels = "all")

# 面を名前で選んで描画する
p_coord + coord_3d(panels = c("xmin", "ymin", "zmin"))

# 軸ラベルを置く辺を手動で指定する
p_coord + coord_3d(xlabels = c("ymax", "zmax"),
                   zlabels = c("xmax", "ymin"))

# 軸ラベルの自動回転を止めて水平に固定する
p_coord + coord_3d(rotate_labels = FALSE)

# 少し引きの構図にして余白を広げる
p_coord + coord_3d(zoom = 0.85)

・コード例の上から4つ目までのプロットを紹介

立方体の面と軸ラベルの見た目をthemeで整える：cube_themingコマンド

cube_themingは「ggcube」パッケージが追加するテーマ要素のヘルプページです。?cube_theming()で実行します。そうすると、ヘルプページが表示されます。

下記にヘルプページに記載のオプションをggplot2の「theme」コマンドで使った例を紹介します。

# 6面すべてを描いたうえで、面とラベルの体裁を細かく指定する
ggplot(terrain, aes(x, y, z)) +
  geom_surface_3d(fill = "steelblue", color = "white", linewidth = 0.15) +
  coord_3d(panels = "all") +
  theme(
    # 奥側の面の塗りと枠線
    panel.background = element_rect(fill = "gray95", color = "gray40"),
    # 手前の面はデータが隠れないよう強めに透過させる
    panel.foreground = element_rect(alpha = 0.15),
    # 手前の面のグリッド線
    panel.grid.foreground = element_line(color = "gray70", linewidth = 0.25),
    # 目盛りラベルの色
    axis.text = element_text(color = "gray20"),
    # z軸の目盛りラベルだけ色を変える
    axis.text.z = element_text(color = "firebrick"),
    # marginで軸タイトルの余白を広げる
    axis.title = element_text(margin = margin(t = 20)),
    # x軸タイトルだけ色を変える
    axis.title.x = element_text(color = "steelblue4"))

＜おすすめのRに関する書籍です＞

The R Tips 第3版 データ解析環境Rの基本技・グラフィックス活用集

R は統計計算とグラフィックスのための言語・環境で、R の熟練度に合わせて読んでいただけるよう、2部構成にしています。初心者向けにはR入門編、統計や数値計算、プログラミング技法やデータハンドリング等、グラフィックスなど実行したい人はR-Tips 編で学ぶことができます。 著者: 舟尾暢男

www.amazon.co.jp

度数集計やビニングを自動で行い3D棒グラフを描く：geom_bar_3dコマンド

データを集計した3次元の棒グラフをプロットするコマンドです。x軸とy軸が離散値なら組み合わせごとの度数を数え、連続値なら2次元のビンに区切って集計し、片方だけが連続値なら群を保ったままビンを設定します。ggplot2のgeom_bar()とgeom_histogram()のコマンドを3Dでプロットします。集計済みの高さを使いたい場合はgeom_col_3d()のコマンドを使います。

# データ作成：店舗と曜日の組み合わせ
set.seed(1)
visits <- data.frame(
  store = sample(c("Sapporo", "Asahikawa", "Hakodate"), 600, replace = TRUE),
  weekday = sample(c("Mon", "Wed", "Fri", "Sun"), 600, replace = TRUE))

# 離散×離散なので、組み合わせごとの度数がそのまま柱の高さになる
ggplot(visits, aes(store, weekday)) +
  geom_bar_3d(fill = "steelblue") +
  coord_3d()

# データ作成：相関のある2つの計測値
set.seed(2)
meas <- data.frame(sensor_a = rnorm(1500))
# sensor_aに比例する成分とノイズを足してsensor_bを作る
meas$sensor_b <- 0.6 * meas$sensor_a + rnorm(1500, 0, 0.8)

# 連続×連続なので2次元ビンに区切られ、3Dヒストグラムになる
ggplot(meas, aes(sensor_a, sensor_b)) +
  geom_bar_3d(bins = c(18, 18), fill = "darkorange") +
  coord_3d()

# 空のビンも高さ0の柱として表示する
ggplot(meas, aes(sensor_a, sensor_b)) +
  geom_bar_3d(bins = c(12, 12), drop = FALSE) +
  coord_3d()

# データ作成：3系統の測定値
set.seed(3)
lots <- data.frame(
  lot = rep(c("L1", "L2", "L3"), each = 200),
  value = c(rnorm(200, 5, 1), rnorm(200, 7, 1.5), rnorm(200, 6, 0.7)))

# 離散×連続の組み合わせでは、連続側だけがビンに使用される
ggplot(lots, aes(x = lot, y = value, fill = lot)) +
  geom_bar_3d(bins = 20, width = c(0.6, 1)) +
  coord_3d(scales = "fixed", ratio = c(1, 1, 0.1))

# 度数の代わりに密度をz軸に割り当てる
ggplot(lots, aes(x = lot, y = value, z = after_stat(density))) +
  geom_bar_3d(bins = 20, width = c(0.6, 1)) +
  coord_3d()

実行例：詳細は上記コマンド例を確認してください。

回帰モデルの当てはめ結果を3D曲面として描く：geom_smooth_3dコマンド

散布した(x, y, z)データに平滑化モデルを当てはめ、その予測値を格子上で評価した3次元曲面がプロットできます。ggplot2のgeom_smooth()の3Dで、loess、lm、glm、gamに対応しています。信頼区間を追加の曲面として重ねたり、元のデータ点や残差線を一緒にプロット可能です。

### データの作成：温度と湿度から収量が決まる栽培試験#####
set.seed(123)
crop <- data.frame(
  temp = runif(120, 10, 30),
  humid = runif(120, 30, 90))
# 最適点から離れるほど収量が落ちる二次の関係にノイズを加える
crop$yield <- 50 - 0.08 * (crop$temp - 22)^2 -
  0.01 * (crop$humid - 60)^2 + rnorm(120, 0, 2)
########

# 曲面を重ねる土台となるプロットを作成
p_smooth <- ggplot(crop, aes(temp, humid, yield)) +
  coord_3d() +
  scale_fill_viridis_c()

# 既定のloessによる平滑化曲面を描く
p_smooth + geom_smooth_3d()

# 元のデータ点と、点から曲面へ下ろした残差線を重ねる
p_smooth + geom_smooth_3d(points = TRUE, residuals = TRUE,
                          point_color = "red", alpha = 0.8)

# 二次の項を含むモデル式を明示して線形モデルを当てはめる
p_smooth + geom_smooth_3d(method = "lm", n = 12,
                          formula = z ~ poly(x, 2) + poly(y, 2))

# loessの平滑化の強さをspanで調整し、照明効果を加える
p_smooth + geom_smooth_3d(method = "loess",
                          method.args = list(span = 0.5),
                          fill = "steelblue", color = "white", n = 20,
                          light = light(direction = c(0, -1, 1)))

# 95％信頼区間を上下の曲面として重ねて描く
p_smooth + geom_smooth_3d(aes(fill = after_stat(level)),
                          method = "lm", color = "black", se = TRUE,
                          level = 0.95, se_alpha = 0.6, n = 10) +
  scale_fill_manual(values = c("tomato", "gray40", "steelblue"))

# 標準誤差そのものを塗り分けに使い、予測の不確実さを可視化する
p_smooth + geom_smooth_3d(aes(fill = after_stat(se)), n = 15) +
  scale_fill_viridis_c(option = "magma")

# データが存在しない隅への外挿を避けるため、凸包の内側だけを描く
p_smooth + geom_smooth_3d(method = "lm", domain = "chull", n = 15)

# データ範囲の外側まで格子を広げて、明示的に外挿させる
p_smooth + geom_smooth_3d(method = "lm", n = 12,
                          xlim = c(5, 35), ylim = c(20, 100))

・95％信頼区間を上下の曲面として重ねて描く：他の実行例は省略

z軸の目盛りやラベル、範囲、変換を指定する：scale_z_continuousコマンド

scale_x_continuous()やscale_y_continuous()と同じ感覚で、目盛りの位置、ラベルの文字列、軸の範囲、対数や逆順といった変換を指定できます。

### データの作成：製品の価格・重量・連続使用時間#####
set.seed(7)
product <- data.frame(
  price = round(runif(60, 1000, 5000)),
  weight = round(runif(60, 100, 900)),
  hours = round(runif(60, 5, 40)))
########

# z軸の目盛り位置とラベル、下限を明示して描画する
ggplot(product, aes(price, weight, hours)) +
  geom_point_3d(color = "steelblue") +
  scale_z_continuous(
    name = "Battery life",
    breaks = c(10, 20, 30, 40),
    labels = c("10h", "20h", "30h", "40h"),
    limits = c(0, NA)) +
  coord_3d()

# z軸を逆順にして、値が小さいほど上に来るようにする
ggplot(product, aes(price, weight, hours)) +
  geom_point_3d(color = "firebrick") +
  scale_z_continuous(transform = "reverse") +
  coord_3d()

# 目盛りの個数だけを指定して、位置は自動で決めさせる
ggplot(product, aes(price, weight, hours)) +
  geom_point_3d() +
  scale_z_continuous(n.breaks = 8) +
  coord_3d()

・z軸を逆順にして、値が小さいほど上に来るようにする：他の実行例は省略

＜おすすめのRに関する書籍です＞

Rクックブック 第2版 | J.D.Long, Paul Teetor, 大橋真也

AmazonでJ.D.Long, Paul Teetor, 大橋真也のRクックブック 第2版。アマゾンならポイント還元本が多数。J.D.Long, Paul Teetor, 大橋真也作品ほか、お急ぎ便対象商品は当日お届けも可能。またRクックブック 第2版もアマゾン配送商品なら通常配送無料。

www.amazon.co.jp

2次元のプロットを簡単に3D表現にする：aesコマンド

ggcubeのaes()は、aes(x, y, z)と3つ並べて書くだけでz軸を割り当てることができ、2次元のプロットを簡単に3D表現にすることが可能です。

実行結果のプロットは省略します。便利なコマンドです。

### データの作成：製品の価格・重量・連続使用時間#####
set.seed(7)
product <- data.frame(
  price = round(runif(60, 1000, 5000)),
  weight = round(runif(60, 100, 900)),
  hours = round(runif(60, 5, 40)))
########

# 2つだけ引数を渡した場合は、通常のggplot2と同じ2次元プロットになる
ggplot(product, aes(price, weight)) +
  geom_point()

# 3つ目の引数が自動的にz軸へ割り当てられる
ggplot(product, aes(price, weight, hours)) +
  geom_point() +
  coord_3d()

# z = を明示した書き方でも結果は同じ
ggplot(product, aes(x = price, y = weight, z = hours)) +
  geom_point() +
  coord_3d()

回転させたコマ送りをGIFや動画として書き出す：animate_3dコマンド

coord_3d()を使ったプロットの回転をフレーム間で補間しながら描画し、アニメーションとして出力します。pitch、roll、yawに長さ2以上の数値ベクトルを渡すと、その値をキーフレームとして線形補間されます。既定のGIF出力にはgifskiパッケージ、動画出力にはavパッケージが別途必要です。

なお、RStudioを使用するとViewerペインにプロットが表示されます。そのまま、Show in new windowボタンをクリックすると、ブラウザでプロットが表示されます。ブラウザからデスクトップにドラッグアンドドロップするとアニメーションgifの保存が可能です。保存はanim_save_3dコマンドでも可能です。

# GIF出力にはgifskiパッケージが必要
install.packages("gifski")

# アニメーションさせる地形プロットを作成
# anchor = "camera"を指定すると、回転しても光源が固定される
p_anim <- ggplot(terrain, aes(x, y, z)) +
  geom_surface_3d(fill = "steelblue", color = "steelblue",
                  linewidth = 0.2) +
  coord_3d(light = light(mode = "hsl", anchor = "camera",
                         direction = c(1, 1, 0))) +
  theme_void()

# z軸まわりに1周させる（実行時間を抑えるためフレーム数は少なめ）
animate_3d(p_anim, yaw = c(-30, 330), nframes = 24, fps = 12)

# yawで2周させながら、途中でrollを変化させる
animate_3d(p_anim, yaw = c(0, 720), roll = c(-90, -30, -90),
           nframes = 36, fps = 15)

# 開始と終了で静止させ、逆再生を付けてループさせる
animate_3d(p_anim, yaw = c(-30, 90), nframes = 24,
           start_pause = 5, end_pause = 5, rewind = TRUE)

# 作成したアニメーションをGIFとして保存
anim <- animate_3d(p_anim, yaw = c(0, 360), nframes = 24)
anim_save_3d(anim, file.path(tempdir(), "terrain_rotation.gif"))

＜おすすめのRに関する書籍です＞

初心者でもすぐにできるフリー統計ソフトEZR(Easy R)で誰でも簡単統計解析(改訂第2版) | 神田善伸

Amazonで神田善伸の初心者でもすぐにできるフリー統計ソフトEZR(Easy R)で誰でも簡単統計解析(改訂第2版)。アマゾンならポイント還元本が多数。

www.amazon.co.jp

この記事が誰かの役に立ちますように。