Rで解析：フォレストプロットを作成「ggforestplotR」パッケージの紹介

複数の変数の効果量と信頼区間を一覧できるフォレストプロットは便利ですが、作図用のデータを整え、体裁の整った図に仕上げるには手間がかかります。本パッケージは、データを作図用に整えるコマンドが収録されているだけでなく、「ggplot2」パッケージのコマンドを適応できるので図の体裁を整えるのも簡単です。

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

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

データサイエンティストのための最新知識と実践 Rではじめよう! [モダン]なデータ分析

モダンな書き方でデータ分析を進めるための一冊。今どきのRの流儀を押さえたい人に。

www.amazon.co.jp

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

下記コマンドを実行してください。

# パッケージのインストール
install.packages("ggforestplotR")
# パッケージの読み込み
library("ggforestplotR")
# コマンド例に必要なパッケージ
# install.packages("broom")
library("broom")
# install.packages("tidyverse")
library("tidyverse")

コマンド例

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

フォレストプロットの作成手順は、「tidy_forest_model」と「bind_forest_models」を使用してモデルオブジェクトを整理し、次に「add_forest_table」や「add_split_table」を使って要約表や分割表を追加し、最後に「ggforestplot」です。

複数モデルの係数を1つのフォレストプロット用データに結合：bind_forest_modelsコマンド

コマンド例ではggplot2パッケージのコマンドで色・ラベル・体裁を整えています。

# 基本モデルを作成
fit1 <- lm(mpg ~ wt + hp, data = mtcars)
# 調整済みモデルを作成
fit2 <- lm(mpg ~ wt + qsec, data = mtcars)

# 2つのモデルの係数を1つのフォレストプロット用データに結合する（group列にモデル名が入る）
bound <- bind_forest_models(list(Base = fit1, Adjusted = fit2))

# フォレストプロットを作成（点・基準線・行の塗り分けはこの関数の引数で指定する）
p <- ggforestplot(
  bound,
  striped_rows = TRUE,      # 行を交互に塗り分けて見やすくする
  stripe_fill  = "grey96",  # 塗り分けの色
  point_size   = 2.8,       # 点（推定値）のサイズ
  ref_line     = 0,         # 効果なしを示す基準線（加法モデルなので0）
  ref_color    = "grey70"   # 基準線の色
)

# ggplot2パッケージのコマンドで色・ラベル・体裁を整える
p <- p +
  # モデル名（group）の色を手動で指定し、凡例タイトルも付ける
  scale_colour_manual(
    name   = "モデル",
    values = c(Base = "#2C7FB8", Adjusted = "#D95F0E")
  ) +
  # タイトルや軸ラベルを付ける
  labs(
    title    = "回帰係数のフォレストプロット",
    subtitle = "mtcars を用いた2モデルの比較",
    x        = "推定値（95%信頼区間）",
    caption  = "データ: mtcars"
  ) +
  # 全体の体裁を微調整する（theme_bwを上書きせず追記する）
  theme(
    plot.title       = element_text(face = "bold"),
    legend.position  = "top",
    panel.grid.minor = element_blank()
  )

# 図を表示
p

単一のモデルをフォレストプロット用データに整形：tidy_forest_modelコマンド

モデルの結果をわかりやすく整理し、フォレストプロットで可視化したい場合に便利です。

# broomパッケージが利用可能な場合のみ実行
if (requireNamespace("broom", quietly = TRUE)) {
  # 線形回帰モデルを作成
  fit <- lm(mpg ~ wt + hp + qsec, data = mtcars)
  # フォレストプロット用データに整形する
  tidy_forest_model(fit)

  # 乱数のシードを固定
  set.seed(123)
  # ロジスティック回帰用のサンプルデータを作成
  logit_data <- data.frame(
    age = rnorm(250, mean = 62, sd = 8),
    bmi = rnorm(250, mean = 28, sd = 4),
    treatment = factor(rbinom(250, 1, 0.45), labels = c("Control", "Treatment"))
  )
  # 線形予測子を計算
  linpred <- -9 + 0.09 * logit_data$age + 0.11 * logit_data$bmi +
    0.9 * (logit_data$treatment == "Treatment")
  # 線形予測子からイベントの発生有無を生成
  logit_data$event <- rbinom(250, 1, plogis(linpred))
  # ロジスティック回帰モデルをfitする
  logit_fit <- glm(event ~ age + bmi + treatment, data = logit_data, family = binomial())

  # オッズ比として表示するため指数変換して整形する
  tidy_forest_model(logit_fit, exponentiate = TRUE)
}

実行結果：

> tidy_forest_model(fit)
  term    estimate   conf.low   conf.high label group grouping separate_groups    n events      p.value
1   wt -4.35879720 -5.9006341 -2.81696034    wt  <NA>     <NA>            <NA> <NA>   <NA> 3.217222e-06
2   hp -0.01782227 -0.0485098  0.01286526    hp  <NA>     <NA>            <NA> <NA>   <NA> 2.441762e-01
3 qsec  0.51083369 -0.3888708  1.41053822  qsec  <NA>     <NA>            <NA> <NA>   <NA> 2.546284e-01
   std.error statistic
1 0.75270039 -5.790879
2 0.01498117 -1.189645
3 0.43922153  1.163043

> tidy_forest_model(logit_fit, exponentiate = TRUE)
                term estimate conf.low conf.high label group grouping separate_groups    n events
1                age 1.094012 1.056100  1.136280 age  <NA>     <NA>            <NA> <NA>   <NA>
2                bmi 1.109515 1.035159  1.192999 bmi  <NA>     <NA>            <NA> <NA>   <NA>
3 treatmentTreatment 1.905342 1.112625  3.299050 treatmentTreatment  <NA>     <NA>            <NA> <NA>   <NA>
       p.value  std.error statistic
1 1.371793e-06 0.01860639  4.829060
2 3.961390e-03 0.03606915  2.881220
3 1.981210e-02 0.27669215  2.329887

フォレストプロットに要約表を追加：add_forest_tableコマンド

# サンプルの係数データフレームを作成
coefs <- data.frame(
  term = c("Age", "BMI", "Treatment"),
  estimate = c(0.3, -0.2, 0.4),
  conf.low = c(0.1, -0.4, 0.2),
  conf.high = c(0.5, 0.0, 0.6),
  sample_size = c(120, 115, 98),
  p_value = c(0.012, 0.031, 0.004))

# フォレストプロットを作成
p <- ggforestplot(coefs, n = "sample_size", p.value = "p_value")

# プロットの左側に要約表を追加
add_forest_table(
  p,
  position = "left",
  columns = c("term", "n", "estimate", "p"),
  estimate_label = "Beta")

# 「+」演算子を使って、プロットの右側に要約表を追加する場合
ggforestplot(coefs, n = "sample_size", p.value = "p_value") +
  add_forest_table(
    position = "right",
    columns = c("term", "n", "estimate", "p"),
    estimate_label = "Beta")

・プロットの左側に要約表を追加

・「+」演算子を使って、プロットの右側に要約表を追加する場合

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

RとRStudioによる教育テストデータの分析

テストの得点データを扱うための一冊。教育評価の現場で必要になる分析をRとRStudioで実行する。

www.amazon.co.jp

フォレストプロットの前後に表を分割して配置：add_split_tableコマンド

# サンプルの係数データフレームを作成
coefs <- data.frame(
  term = c("Age", "BMI", "Treatment"),
  estimate = c(0.3, -0.2, 0.4),
  conf.low = c(0.1, -0.4, 0.2),
  conf.high = c(0.5, 0.0, 0.6),
  sample_size = c(120, 115, 98),
  p_value = c(0.012, 0.031, 0.004)
)

# フォレストプロットを作成
p <- ggforestplot(coefs, n = "sample_size", p.value = "p_value")

# プロットの左右に列名を指定して表を分割配置
add_split_table(
  p, left_columns = c("term", "n"),
  right_columns = c("estimate", "p"),
  estimate_label = "HR")

# 「+」演算子を使い、列番号を指定して表を分割配置する場合
# プロットの左右に列名を指定して表を分割配置と同じ
# ggforestplot(coefs, n = "sample_size", p.value = "p_value") +
#  add_split_table(
#    left_columns = c(1, 5),
#    right_columns = c(2, 6),
#    estimate_label = "HR")

・プロットの左右に列名を指定して表を分割配置

係数データからフォレストプロットを描画：ggforestplotコマンド

# サンプルの係数データフレームを作成
coefs <- data.frame(
  term = c("Age", "BMI", "Treatment"),
  estimate = c(0.10, -0.08, 0.34),
  conf.low = c(0.02, -0.16, 0.12),
  conf.high = c(0.18, 0.00, 0.56))

# 初期設定でフォレストプロット
ggforestplot(coefs)

# 交互の行ストライプと点スタイルを変えてフォレストプロット
ggforestplot(coefs, striped_rows = TRUE, point_shape = 17)

係数データフレームをフォレストプロット用の標準形式に変換：as_forest_dataコマンド

# サンプルデータを作成
raw <- data.frame(
  variable = c("Age", "BMI", "Treatment"),
  beta = c(0.10, -0.08, 0.34),
  lower = c(0.02, -0.16, 0.12),
  upper = c(0.18, 0.00, 0.56))

# 列名を、フォレストプロット対応の列名に変換する
as_forest_data(
  data = raw,
  term = "variable",
  estimate = "beta",
  conf.low = "lower",
  conf.high = "upper")

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

Rによる教育・言語・心理系のためのデータサイエンス入門 | 柳川浩三

教育・言語・心理の研究データを想定したデータサイエンス入門。分野特有のデータの扱い方に触れられる。

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