【Forest Plot】Rでフォレストプロットの書き方：臨床研究のサブグループ解析

# 練習用データの作成
set.seed(123) # 再現性のためのシード
n <- 200 # 患者数

# 各変数を生成
PatientID <- 1:n
Age <- round(runif(n, 50, 80)) # 年齢: 50歳から80歳
Sex <- sample(c("Male", "Female"), n, replace = TRUE)
Hypertension <- sample(c(0, 1), n, replace = TRUE, prob = c(0.6, 0.4))
Diabetes <- sample(c(0, 1), n, replace = TRUE, prob = c(0.7, 0.3))
Group <- sample(c("NewDrug", "Placebo"), n, replace = TRUE)
Time <- round(runif(n, 1, 52)) # フォローアップ期間（週）
Event <- rbinom(n, 1, ifelse(Group == "NewDrug", 0.15, 0.25)) # イベント発生率の差を反映

# データフレームを作成
datasheet1 <- data.frame(PatientID, Age, Sex, Hypertension, Diabetes, Group, Time, Event)

#AgeGroup列を追加
datasheet1$AgeGroup <- ifelse(datasheet1$Age < 65, "<65", ">=65")

library(tidyverse)
library(survival)

perform_analysis <- function(data, group_var = NULL) {
  if (!is.null(group_var)) {
    # サブグループ解析
    results <- by(data, data[[group_var]], function(subgroup_data) {
      analyze_cox_model(subgroup_data, subgroup_data[[group_var]][1])
    })
    do.call(rbind, results)
  } else {
    # 全患者解析
    analyze_cox_model(data, "All patients")
  }
}

analyze_cox_model <- function(data, subgroup_name) {
  
  # サバイバルオブジェクトを作成(ここのTime, Eventは自分のデータの列名に応じて変更してください)
  surv_obj <- Surv(time = data$Time, event = data$Event)
  
  # Cox比例ハザードモデルを適用(ここのGroupは自分のデータの列名に応じて変更してください)
  cox_model <- coxph(surv_obj ~ Group, data = data)
  
  # HR, 95%CI, p値を取得
  hr <- exp(coef(cox_model))
  ci <- exp(confint(cox_model))
  pval <- summary(cox_model)$coefficients[5]
  
  # 各群のN数を計算(ここのGroup, "NewDrug", "Placebo"は自分のデータの列名や内容に応じて変更してください)
  counts <- table(data$Group)
  new_drug_count <- counts["NewDrug"]
  placebo_count <- counts["Placebo"]
  
  # 結果をデータフレームとして返す
  data.frame(
    Subgroup = subgroup_name,
    NewDrug = ifelse(is.na(new_drug_count), 0, new_drug_count),
    Placebo = ifelse(is.na(placebo_count), 0, placebo_count),
    HR = hr,
    CI_Lower = ci[1],
    CI_Upper = ci[2],
    p_value = pval
  )
}

# list()の中身は自分の調べたいサブグループの列名に応じて変更してください
subgroups <- list(
  AgeGroup = "AgeGroup",
  Sex = "Sex",
  Hypertension = "Hypertension",
  Diabetes = "Diabetes"
)

# datasheet1はデータシート名に応じて変更してください
subgroup_results <- lapply(names(subgroups), function(var) {
  cbind(Variable = var, perform_analysis(datasheet1, subgroups[[var]]))
})

all_patients_result <- cbind(Variable = "All patients", perform_analysis(datasheet1))

all_results <- do.call(rbind, c(list(all_patients_result), subgroup_results))

# 空白行を作成(適宜列名に応じて以下の内容は変更してください)
empty_rows <- data.frame(
  Variable = NA, Subgroup = NA, NewDrug = NA, Placebo = NA,
  HR = NA, CI_Lower = NA, CI_Upper = NA, p_value = NA,
  stringsAsFactors = FALSE
)

# 空白行を4つ複製(今回4因子でサブグループ解析したので4つ複製します)
empty_rows <- empty_rows[rep(1, 4), ]

# Subgroup列にリストの値を挿入(このsubgroupは【サブグループ解析の実行】で作成したリストに対応しています)
empty_rows$Subgroup <- unlist(subgroups)

# all_results に新しい行を追加
all_results <- rbind(all_results, empty_rows)

# 指定した行番号順(行をマニュアルで指定して並び替えます)
row_order <- c(1, 10, 2, 3, 11, 4, 5, 12, 6, 7, 13, 8, 9)

# 並び替え
all_results <- all_results[row_order, ]

# 最初の列と最後の列を同時に削除
all_results <- all_results[, -c(1, ncol(all_results))]

# 新しい列 se(ハザード比の標準誤差) を作成
all_results$se <- (log(all_results$CI_Upper) - log(all_results$HR)) / 1.96

# 結果を確認
view(all_results)

# 空白行を20個作成。もしフォレストプロットの横線が足りなければ30,40などに増やす
all_results$` ` <- paste(rep(" ", 20), collapse = " ")

# 描画用のHR(95% CI)を作成する
all_results$`HR (95% CI)` <- ifelse(is.na(all_results$se), "",
                           sprintf("%.2f (%.2f to %.2f)",
                                   all_results$HR, all_results$CI_Lower, all_results$CI_Upper))

# Subgroup列の名前を変更する(0, 1をNo, Yesに変更する)
all_results$Subgroup <- ifelse(all_results$Subgroup == 0, "No", 
                                ifelse(all_results$Subgroup == 1, "Yes", all_results$Subgroup))

# 見やすい用にSubgroupのインデント調整(Placeboに数値があればインデント"   "を入れるようにする)
all_results$Subgroup <- ifelse(is.na(all_results$Placebo), 
                      all_results$Subgroup,
                      paste0("   ", all_results$Subgroup))
                      
# NAが表示されないように消す（NewDrug, Placeboの列にあるNAを空白にする）
all_results[c("NewDrug", "Placebo")] <- lapply(all_results[c("NewDrug", "Placebo")], function(x) {
  x[is.na(x)] <- ""
  x
})

# 結果を確認
view(all_results)

p <- forest(all_results[,c(1:3, 8:9)],
            est = all_results$HR,
            lower = all_results$CI_Lower, 
            upper = all_results$CI_Upper,
            sizes = all_results$se,
            ci_column = 4,
            ref_line = 1,
            arrow_lab = c("NewDrug Better", "Placebo Better"),
            xlim = c(-1, 4),
            ticks_at = c(0, 1, 2, 3),
            footnote = "フォレストプロットが完成したよ！")

# Print plot
plot(p)