MENU

【無料プレゼント付き】学会発表・論文投稿に必要な統計を最短で学ぶことができる無料メルマガ

30項目以上の背景因子と非線形リスクをどう解析すべきか?Random Survival Forest(RSF)とrandomForestSRCパッケージによる解決法

「先生!多変量Cox比例ハザードモデルを組もうとしたら、比例ハザード性(PH仮定)が複数の変数で崩れてしまい、単一のハザード比(HR)が出せなくなりました!おまけに背景因子が30項目以上あって、ステップワイズ法を入れてもモデルが過学習で崩壊してしまいます……」

多変量Cox回帰を用いて予後予測や危険因子探索を行っていると、こうした「手詰まり」に直面することが少なくない。

臨床データやオミクスデータが高次元化・多項目化する現代において、従来のCox回帰が前提とする「比例ハザード性」や「単一の線形関係」にデータを無理やり当てはめることには限界がある。

しかし、「ハザード比(HR)が出ないから機械学習は医学論文で使えない」と思い込んで諦める必要は全くない。

右側打切り(Censoring)に対応した機械学習手法である Random Survival Forest(ランダムサバイバルフォレスト / RSF) を用いれば、PH仮定の破綻や複雑な非線形関係をクリアした上で、変数重要度(VIMP)Minimal Depth(最小深さ)Partial Dependence Plot(部分依存プロット) という3つの強力な武器によって、Cox回帰以上に柔軟かつ説得力を持って危険因子を証明できる。

本記事では、RSFの直感的な仕組みから、機械学習の初学者に向けた基本用語解説、標準パッケージ randomForestSRC による実践手順、そして論文でそのまま使える英文テンプレートまでを徹底的に解説する。

>>もう統計で悩むのは終わりにしませんか? 

↑1万人以上の医療従事者が購読中

目次

1. なぜ多変量Cox回帰は限界を迎えるのか?(PH仮定破綻と高次元データの壁)

医学研究において最も広く使われている多変量解析手法はCox比例ハザードモデルである。しかし、患者背景の複雑化に伴い、従来のCox回帰では対応しきれない「3つの致命的な限界」が顕在化している。

1. 比例ハザード性(PH仮定)の破綻

Cox回帰は「観察期間を通じて、2群間のハザード比(HR)が常に一定である」という極めて厳格な前提条件(比例ハザード性)に基づいている。しかし、追跡期間が長くなったり、治療薬の初期効果と長期効果が異なったりする場合、ハザード比は時間とともに変化する。

共変量が20〜30項目を超えると、どれか複数の変数でPH仮定が破綻し、妥当な単一のHRを算出できなくなる。

※補足(関連記事):

比例ハザード性の破綻時の対策やRMST解析については Cox回帰での比例ハザード性破綻とRMST を参照されたい。

2. 高次交互作用と非線形性の無視

臨床現場では「単独のマーカー上昇」ではなく、「70歳以上の高血圧患者でのみ、特定の炎症マーカーが高くなると死亡リスクが跳ね上がる」といった条件付きのリスク上昇(高次交互作用)が頻繁に起こる。

従来のCox回帰でこれを評価するには、研究者が手動で「交互作用項(Interaction terms)」をモデルに組み込む必要があるが、あらかじめ全ての組み合わせを想定することは事実上不可能である。

3. 多重共線性と過学習(Overfitting)の罠

共変量が背景項目や検査値など数十項目に及ぶ場合、変数間の相関(多重共線性)によりモデルが不安定になる。ステップワイズ法などで無理やり変数選択を行うと、サンプルサイズに対してモデルが複雑になりすぎ、未知のデータに対する予測能を失う「過学習(Overfitting)」を引き起こす。

【コラム1】臨床医のための機械学習用語集(基本編)

機械学習に馴染みのない読者向けに、本記事で登場する基本概念を臨床的なイメージで解説する。

  • 決定木(Decision Tree):「年齢 > 65 か? Yes/No」$\rightarrow$「CRP > 2.0 か? Yes/No」のように、YES/NO の質問(分岐)を繰り返して患者をリスクごとにグループ分けするフローチャート。
  • アンサンブル学習(Ensemble Learning):1人の医師(1本の決定木)の診断に頼るのではなく、1000人の医師(1000本の決定木)に少しずつ異なる患者データを見せて判定させ、その「多数決(平均)」で最終的なリスク予測を弾き出す手法。「森(Forest)」と呼ぶのは、無数の「木(Tree)」が集まっているためである。
  • ブートストラップ標本(Bootstrap Sample): 元の臨床データから重複を許してランダムに患者を抽出(復元抽出)し、少しずつ顔ぶれの違う「学習用データセット」を何パターンも作成する手法。

2. Random Survival Forest (RSF) の仕組み:打切りデータを決定木で扱う技法

通常のRandom Forestは「がん発生(1)/ なし(0)」のような二値分類や連続値の予測を行う手法であり、追跡期間(Time)や「打切り(Censoring)」が存在する生存時間解析にはそのまま適用できない。

この問題を解決したのが Random Survival Forest(RSF) である。

打切りを正しく扱う「Log-rank split rule(ログランク分割ルール)」

RSFが生存時間データを正しく扱える秘密は、決定木を伸ばしていく際の「分岐ルール」にある。

一般的な決定木がデータの不純度(Gini係数など)を下げるように分けるのに対し、RSFでは各ノード(分岐点)で Log-rank split rule(ログランク分割ルール) を採用する。

つまり、「患者群をどのような条件(例: 血圧 140 以上 vs 未満)で2分割したときに、2群間のログランク統計量が最も大きくなるか(生存曲線が最も大きく離れるか)」を検索し、分岐点を作成する。

打切りデータを正しく考慮したログランク検定に基づいて樹形図(Survival Tree)を構築するため、打切りを含む臨床データからでも正確な累積ハザードや予測生存率をアンサンブル(集約)することができる。

>>もう統計で悩むのは終わりにしませんか? 

↑1万人以上の医療従事者が購読中

3. ハザード比(HR)が出なくても査読者を黙らせる「3つの武器」

「RSFはハザード比(HR)が出力されないから、医学論文(査読)で使えないのではないか?」という質問をよく受ける。

しかし、これは完全な誤解である。査読者が本質的に知りたいのは「単一のHRという数値」そのものではなく、「どの因子が真の予後予測因子か」「その因子がリスクにどう影響しているか」「モデルの予測精度は高いか」 の3点である。

RSFでは、Cox回帰のHRに代わる「3つの強力な可視化ツール」を用いることで、査読者を十分に納得させるロジックを構築できる。

武器1:VIMP(Variable Importance / 変数重要度)

特定の変数(例: Age)の値をランダムにシャッフル(破棄)したときに、モデル全体の予測精度がどれくらい悪化するかを測定する指標。

もし、ある変数を破棄したときに予測精度が著しく低下(エラー率が上昇)した場合、その変数は「予後予測に極めて重要な因子である」と判定できる。

武器2:Minimal Depth(最小深さ)

決定木の中で、その変数が「どれくらい浅い(初期の)階層」で分割に使われたかを示す指標。

予後を決定づける強固なリスク因子ほど、決定木の上部(浅い階層)で頻繁に選択される。VIMPとMinimal Depthの双方で上位に入った因子は、極めて確実性の高いリスク因子として主張できる。

武器3:Partial Dependence Plot(PDP / 部分依存プロット)

他のすべての変数の影響を一定に保った(周辺化した)上で、特定のリスク因子(例: D-dimer値、CRP値、年齢)の変化に伴い、アンサンブル死亡リスクや予測生存率がどう推移するかを描画するプロット。

「検査値が〇〇を超えると急激に死亡リスクが跳ね上がる」といった非線形な閾値効果(Threshold effect)を視覚的に提示できるため、単一のHR値よりも臨床的に価値の高い情報を提供できる。

【コラム2】臨床医のための機械学習用語集(評価・解釈編)

  • Out-of-Bag (OOB) データ: ブートストラップ標本を作成する際、「たまたま選ばれなかった約37%の患者データ」 のこと。【なぜ約37%になるのか?】 サンプルサイズ$N$人のデータから「重複を許して$N$回ランダムに選ぶ(復元抽出)」とき、ある患者が1回の抽出で『選ばれない確率』は$(1 – 1/N)$である。これを$N$回繰り返したときに『一度も選ばれない確率』は$(1 – 1/N)^N$となる。 サンプルサイズ$N$が十分大きい臨床データでは、数学の性質上、この値はネイピア数の逆数$1/e \approx 0.368$(約36.8%$\approx$37%)へ収束する。RSFでは、この約37%の未学習データ(OOB)をテスト用データとして用いて精度検証を行うため、別途検証用データを用意しなくても過学習(Overfitting)のない客観的な予測精度を算出できる。
  • C-index (Harrell’s Concordance Index):「より高リスクと予測された患者が、実際に早くイベントを起こしたか」という一致率を示す指標(0.5〜1.0)。0.5でコイン投げと同等、0.7以上で良好な予測能、0.8以上で極めて高精度と判定される。

4. 【実践Rコード】randomForestSRC による RSFモデル構築と可視化(Step by Step)

ここからは、CRANで最も安定してメンテされている標準パッケージ randomForestSRC を用いて、RSFモデル構築から精度評価、要因抽出までを Step by Step で解説する。

Step 1: 疑似臨床データの作成と確認

患者200人分、追跡時間(time)、イベント(status)、背景因子10項目の疑似臨床データを生成する。

library(randomForestSRC)
library(tidyverse)

# 乱数シードの設定
set.seed(123)
n <- 200

# 疑似臨床データの作成
df_rsf <- tibble(
  id = 1:n,
  time = round(runif(n, min = 1, max = 60), 1),
  status = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.4, 0.6)),
  age = round(rnorm(n, mean = 65, sd = 10)),
  gender = factor(sample(c("Male", "Female"), n, replace = TRUE)),
  bmi = round(rnorm(n, mean = 24, sd = 4), 1),
  sbp = round(rnorm(n, mean = 135, sd = 15)),
  crp = round(rgamma(n, shape = 2, scale = 1.5), 2),
  eGFR = round(rnorm(n, mean = 70, sd = 15), 1),
  bnp = round(rlnorm(n, meanlog = 4.5, sdlog = 0.8), 1),
  statine = factor(sample(c("No", "Yes"), n, replace = TRUE)),
  smoking = factor(sample(c("Never", "Former", "Current"), n, replace = TRUE))
)

# 先頭行の確認
head(df_rsf)

Step 2: RSFモデルの構築と OOB C-index による精度評価

rfsrc() 関数を用いて ntree = 1000(1000本の決定木)のフォレストを構築し、Out-of-Bag (OOB) データに基づく C-index(予測精度)を出力する。

# RSFモデルの構築
rsf_fit <- rfsrc(
  Surv(time, status) ~ .,
  data = df_rsf %>% select(-id),
  ntree = 1000,
  splitrule = "logrank",
  importance = TRUE,
  seed = 123
)

# モデル概要の出力
print(rsf_fit)

# OOB C-index の確認 (1 - OOB Error Rate)
oob_c_index <- 1 - rsf_fit$err.rate[rsf_fit$ntree]
cat(sprintf("Out-of-Bag (OOB) C-index: %.3f\n", oob_c_index))

Step 3: 変数重要度(VIMP)と Minimal Depth の抽出・比較

vimp() および max.subtree() を用いて、重要変数を抽出・ランキング化する。

# VIMP (Variable Importance) の取得
vimp_res <- vimp(rsf_fit)$importance

# Minimal Depth (最小深さ) の抽出
max_sub_res <- max.subtree(rsf_fit)
min_depth_res <- max_sub_res$order[,1]

# VIMPとMinimal Depthの結果を統合
df_importance <- tibble(
  variable = names(vimp_res),
  vimp = as.numeric(vimp_res)
) %>%
  left_join(
    tibble(
      variable = names(min_depth_res),
      min_depth = as.numeric(min_depth_res)
    ),
    by = "variable"
  ) %>%
  arrange(desc(vimp))

print(df_importance)

【結果の読み方】VIMPの値が大きいほど予測能への貢献度が高く、Minimal Depthの値が小さい(0に近い)ほど決定木の根元近くで早期に分類に使われた重要な変数となる。双方の指標で上位に入った変数こそが、真に確実性の高い予後予測因子である。

Step 4: Partial Dependence Plot (PDP) による非線形リスク推移の描画

上位に抽出された重要な連続変数について、数値の変化に伴う死亡リスクの非線形な変化を描画する。

# PDP(部分依存プロット)の描画
plot.variable(
  rsf_fit,
  xvar.names = c("eGFR", "sbp", "bnp"),
  partial = TRUE,
  plots.per.page = 3
)

5. 論文でそのまま使える英文テンプレート(Methods & Results)

論文の Methods および Results セクションにそのまま記述できる標準的な英文表現テンプレートを提示する。

Methods セクション英文テンプレート

“To account for potential non-proportional hazards, high-dimensional covariates, and non-linear interactions without making parametric assumptions, a Random Survival Forest (RSF) model was constructed using the randomForestSRC package in R. The forest was built with 1,000 survival trees using the log-rank splitting rule. Model performance was evaluated using the Out-of-Bag (OOB) Harrell’s concordance index (C-index) to prevent overfitting. Variable importance (VIMP) and minimal depth analysis were utilized to identify and rank key prognostic factors. The non-linear relationship between continuous clinical variables and mortality risk was visualized using Partial Dependence Plots (PDP).”

Results セクション英文テンプレート

“The Random Survival Forest model demonstrated robust predictive accuracy, achieving an Out-of-Bag C-index of [0.545]. Among the evaluated clinical factors, [eGFR, SBP, and BNP,] were identified as the most dominant prognostic factors according to both Variable Importance (VIMP) and minimal depth metrics. Partial Dependence Plots revealed a non-linear threshold effect, demonstrating a sharp increase in predicted mortality risk when [eGFR levels dropped below 60 mL/min/1.73 $\mathrm{m}^2$].”

6. まとめ&一括コピペ用Rスクリプト

本記事のポイントおさらい

  1. PH仮定破綻や多項目データにはRSFが超強力: Cox回帰の限界(比例ハザード性の破綻、非線形・高次交互作用、過学習)を自動で回避できる。
  2. 打切りデータもLog-rank splitで完璧に処理: 通常の決定木と異なり、ログランク分割ルールにより生存時間データを正確に集約する。
  3. HRが出なくても「3つの武器」で査読に対応: VIMP(重要度)、Minimal Depth(深さ)、PDP(部分依存プロット)を提示することで、Cox回帰以上に柔軟かつ実務的なリスク評価が可能。

一括実行用Rスクリプト

以下のコードをコピー&ペーストすることで、疑似データ生成からRSFモデル構築、OOB C-index算出、VIMP/Minimal Depth抽出、PDP描画までを一気通貫で走らせることができる。

# ==============================================================================
# Random Survival Forest (RSF) 一括実行スクリプト
# パッケージ: randomForestSRC, tidyverse
# ==============================================================================

if (!requireNamespace("randomForestSRC", quietly = TRUE)) install.packages("randomForestSRC")
if (!requireNamespace("tidyverse", quietly = TRUE)) install.packages("tidyverse")

library(randomForestSRC)
library(tidyverse)

# 1. 疑似臨床データの作成
set.seed(123)
n <- 200

df_rsf <- tibble(
  id = 1:n,
  time = round(runif(n, min = 1, max = 60), 1),
  status = sample(c(0, 1), size = n, replace = TRUE, prob = c(0.4, 0.6)),
  age = round(rnorm(n, mean = 65, sd = 10)),
  gender = factor(sample(c("Male", "Female"), n, replace = TRUE)),
  bmi = round(rnorm(n, mean = 24, sd = 4), 1),
  sbp = round(rnorm(n, mean = 135, sd = 15)),
  crp = round(rgamma(n, shape = 2, scale = 1.5), 2),
  eGFR = round(rnorm(n, mean = 70, sd = 15), 1),
  bnp = round(rlnorm(n, meanlog = 4.5, sdlog = 0.8), 1),
  statine = factor(sample(c("No", "Yes"), n, replace = TRUE)),
  smoking = factor(sample(c("Never", "Former", "Current"), n, replace = TRUE))
)

# 2. RSFモデルの構築
cat("\n--- Random Survival Forest モデル構築中 ---\n")
rsf_fit <- rfsrc(
  Surv(time, status) ~ .,
  data = df_rsf %>% select(-id),
  ntree = 1000,
  splitrule = "logrank",
  importance = TRUE,
  seed = 123
)

# 3. 予測精度(OOB C-index)の評価
oob_c_index <- 1 - rsf_fit$err.rate[rsf_fit$ntree]
cat(sprintf("\nOut-of-Bag (OOB) C-index: %.3f\n", oob_c_index))

# 4. 変数重要度 (VIMP) の算出
vimp_res <- vimp(rsf_fit)$importance

# Minimal Depth (最小深さ) の抽出
max_sub_res <- max.subtree(rsf_fit)
min_depth_res <- max_sub_res$order[,1]

# VIMPとMinimal Depthの結果を統合
df_importance <- tibble(
  variable = names(vimp_res),
  vimp = as.numeric(vimp_res)
) %>%
  left_join(
    tibble(
      variable = names(min_depth_res),
      min_depth = as.numeric(min_depth_res)
    ),
    by = "variable"
  ) %>%
  arrange(desc(vimp))

print(df_importance)

cat("\n--- 変数重要度 (VIMP) と Minimal Depth ランキング ---\n")
print(df_importance)

# 5. PDP(部分依存プロット)の描画
plot.variable(
  rsf_fit,
  xvar.names = df_importance$variable[1:3], # 上位3変数をプロット
  partial = TRUE,
  plots.per.page = 3
)

おすすめ書籍

誰も教えてくれなかった 医療統計の使い分け〜迷いやすい解析手法の選び方が,Rで実感しながらわかる!

よかったらシェアしてね!
  • URLをコピーしました!
  • URLをコピーしました!

リサーチクエスチョン探し?データ分析?論文投稿?、、、で、もう悩まない!

第1章臨床研究ではなぜ統計が必要なのか?計画することの重要性
  • 推定ってどんなことをしているの?
  • 臨床研究を計画するってどういうこと?
  • どうにかして標本平均を母平均に近づけられないか?
第2章:研究目的をどれだけ明確にできるのかが重要
  • データさえあれば解析でどうにかなる、という考え方は間違い
  • 何を明らかにしたいのか? という研究目的が重要
  • 研究目的は4種類に分けられる
  • 統計専門家に相談する上でも研究目的とPICOを明確化しておく
第3章:p値で結果が左右される時代は終わりました
  • アメリカ統計協会(ASA)のp値に関する声明で指摘されていること
  • そうは言っても、本当に有意差がなくてもいいの…?
  • なぜ統計専門家はp値を重要視していないのか
  • 有意差がない時に「有意な傾向があった」といってもいい?
  • 統計を放置してしまうと非常にまずい
第4章:多くの人が統計を苦手にする理由
  • 残念ながら、セミナー受講だけで統計は使えません。
  • インプットだけで統計が使えない理由
  • どうやったら統計の判断力が鍛えられるか?
  • 統計は手段なので正解がないため、最適解を判断する力が必要
第5章:統計を使えるようになるために今日から何をすれば良いか?
  • 論文を読んで統計が使えるようになるための5ステップ
第6章:統計を学ぶために重要な環境
  • 統計の3つの力をバランスよく構築する環境

以下のボタンをクリックして、画面に出てくる指示に従って、必要事項を記入してください。

この記事を書いた人

医師、医学博士(専門は疫学および統計学)。10年以上の研究歴を有し、筆頭著者として10本の英文原著論文を執筆した実績を持つ。現在、外資系製薬企業にて15年以上にわたり、臨床開発やデータ解析、承認審査、薬価交渉戦略といった実務の最前線に従事。

「統計に悩むあなたを救いたい」をコンセプトに、ブログ「統計ER」やYouTube、SNSを通じて、医療従事者や研究者向けに「実務直結の統計ノウハウ」を配信。

学術的な厳密さと、製薬業界での高度な実務経験を融合させた独自の視点から、ブラックボックス化しない統計解析の普及を目指す。オンラインショップ「TKER SHOP」での計算ツール提供や、クラウドワークス等を通じた統計コンサルティング・解析代行(主にRを使用)も幅広く展開中。

コメント

コメントする

目次