ビッグデータ 公衆衛生学 医学 統計学

【Rで実践】医療データ分析入門:カプランマイヤー曲線の描き方

臨床研究やリアルワールドデータ分析において、「新薬Aは従来薬Bよりも患者さんの生存期間を延ばすのか?」「ある遺伝子変異があると、病気の再発までの時間は短くなるのか?」といった問いに答えることは非常に重要です。このような「イベント(死亡、再発など)が発生するまでの時間」を分析する手法を生存時間解析と呼びます。

生存時間解析は非常に強力なツールですが、「統計学は苦手…」「プログラミングなんてやったことがない」と感じる方も多いかもしれません。

しかし、ご安心ください。本記事では、無料で利用できる統計解析ソフトウェア「R」を使い、生存時間解析の最も基本的な手法であるカプランマイヤー(Kaplan-Meier)法ログランク(Log-rank)検定を、誰でも実践できるようにステップ・バイ・ステップで徹底解説します。

この記事を最後まで読めば、以下のことができるようになります。

  • 生存時間解析の基本的な考え方を理解する。
  • RとRStudioを準備し、解析に必要な環境を整える。
  • サンプルデータを使って、カプランマイヤー曲線を描画し、論文に載せられるレベルの綺麗なグラフを作成する。
  • ログランク検定を行い、グループ間の生存時間に統計的な差があるか判断する。

少し読み応えがありますが、医療データ分析の第一歩を踏み出していきましょう。

1. 生存時間解析とカプランマイヤー曲線とは?

まず、具体的な操作に入る前に、基本的な概念を少しだけおさらいしましょう。

生存時間データとは?

通常のデータと少し違う、生存時間データが持つ2つの重要な要素が「時間」と「打ち切り」です。

  • イベント発生までの時間 (Time-to-event):研究の観察開始時点(例:治療開始日、診断日)から、特定のイベント(例:死亡、病気の再発、心筋梗塞の発症)が起こるまでの時間のことです。これを「生存時間」と呼びます。
  • 打ち切り (Censoring):これが生存時間解析を少し複雑に、そして面白くする概念です。研究の観察期間中に、イベントが発生しなかった、あるいはイベントが確認できなかったケースを「打ち切り」と呼びます。

    例えば、以下のようなケースが該当します。

    • 研究期間が終了した時点で、まだ生存している。
    • 引っ越しや自己都合で、追跡(フォローアップ)が不可能になった。
    • 交通事故など、研究対象のイベントとは異なる原因で死亡した。

「打ち切り」になったからといって、そのデータを無視してはいけません。例えば、5年間観察してイベントが起こらなかった患者さんは、「少なくとも5年間は生存していた」という非常に重要な情報を持っています。生存時間解析では、この「打ち切り」データを適切に扱うことで、より正確な分析が可能になります。

カプランマイヤー法とは?

カプランマイヤー法は、この「打ち切り」データを考慮しながら、時間経過とともにイベントを発生せずに生存している個体の割合(生存確率)がどのように変化していくかを推定する手法です。

その結果をグラフにしたものがカプランマイヤー曲線です。

  • 縦軸は生存確率(1.0 = 100% から始まり、時間とともに低下)。
  • 横軸は時間(日数、月数、年数など)。

曲線はイベント(例:死亡)が発生するたびに階段状にガクンと下がります。「打ち切り」があった時点では、生存している人数は減りますが、生存確率は変わりません(曲線は下がりません)。この曲線を見ることで、「治療開始から3年後の生存率は約80%だな」といった情報を視覚的に読み取ることができます。

2. Rの準備をしよう

それでは、いよいよ実践の準備です。今回は、無料で高機能な統計ソフトウェア「R」と、それを快適に使うための統合開発環境「RStudio」を使用します。

なぜRを使うのか?

  • 無料: 何より費用がかかりません。
  • 高機能: 世界中の研究者が開発した豊富なパッケージ(拡張機能)により、最新の解析手法もすぐに利用できます。
  • 可視化に強い: ggplot2などのパッケージを使えば、論文掲載レベルの美しいグラフを自在に作成できます。
  • 再現性: コードを保存しておけば、誰でも、いつでも同じ解析を再現できます。

RとRStudioのインストール

まだインストールしていない方は、以下の公式サイトからダウンロードしてインストールしてください。手順は多くの解説サイトで紹介されているので、ここでは割愛します。必ずRを先にインストールしてから、RStudioをインストールしてください。

解析に必要なパッケージのインストール

Rの機能を拡張するために、「パッケージ」をインストールします。今回は、生存時間解析に必須の3つのパッケージを使います。

  • survival: 生存時間解析のコアとなる伝統的なパッケージ。Surv()survfit()といった基本関数を含みます。
  • survminer: カプランマイヤー曲線を美しく、情報豊かに描画するためのパッケージ。ggsurvplot()関数が非常に強力です。
  • tidyverse: データ整形や作図に便利なパッケージ群。今回は主にggplot2のテーマ機能を利用します。

RStudioを起動し、コンソール(通常は左下のウィンドウ)に以下のコマンドをコピー&ペーストして、Enterキーを押してください。

R

install.packages(c("survival", "survminer", "tidyverse"))

黒い文字でインストール状況が表示され、最後に「... successfully unpacked and MD5 sums checked」のようなメッセージが出れば成功です。一度インストールすれば、次回からはこの作業は不要です。

3. Rで実践!データ解析の手順

準備が整いました。ここからは、具体的なコードと共に分析を進めていきましょう。

Step 1: サンプルデータの準備と読み込み

今回は、survivalパッケージに最初から含まれているlungというデータセットを使います。これは、進行肺がん患者の生存に関するデータで、生存時間解析の練習用として非常に有名です。

まず、インストールしたパッケージをRのセッションに読み込みます。これを「ライブラリの読み込み」と呼びます。

R

# ライブラリの読み込み
library(survival)
library(survminer)
library(tidyverse)

次に、lungデータを変数に読み込み、中身を確認してみましょう。

R

# lungデータを読み込み、内容の最初の6行を表示
data(lung)
head(lung)

コンソールに以下のような結果が表示されるはずです。

  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  66   1       0       90        90       NA      15
4   5   210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0

たくさんの変数がありますが、今回注目するのは以下の3つです。

  • time: 生存日数 (イベント発生または打ち切りまでの時間)
  • status: イベント発生状況 (1 = 打ち切り, 2 = 死亡)
  • sex: 性別 (1 = 男性, 2 = 女性)

このままでは少し分かりにくいので、データを少しだけ整形しましょう。

statusは「0 = 打ち切り, 1 = イベント(死亡)」とするのが一般的です。また、sexも数字のままではなく、「Male」「Female」という因子(カテゴリ変数)に変換しておくと、グラフが分かりやすくなります。

R

# データの整形
lung_cleaned <- lung %>%
  mutate(
    # status: 1(打ち切り)を0に、2(死亡)を1に変換
    status = ifelse(status == 2, 1, 0),
    # sex: 1を"Male", 2を"Female"に変換
    sex = factor(sex, levels = c(1, 2), labels = c("Male", "Female"))
  )

# 整形後のデータを確認
head(lung_cleaned)

これでデータの下準備は万端です。

Step 2: 生存時間オブジェクトの作成

次に、survivalパッケージに「このデータは生存時間データですよ」と教えるために、Surv()関数を使って特殊なオブジェクトを作成します。

引数には、Surv(時間, イベント) の順番で指定します。

R

# 生存時間オブジェクトを作成
Surv(lung_cleaned$time, lung_cleaned$status)

実行すると、[1] 306 455 1010+ 210 883 1022+ ... のように、数字の羅列が表示されます。数字の後に + が付いているものが「打ち切り」データを示しています。

Step 3: カプランマイヤー法によるモデルの構築

いよいよカプランマイヤー法で生存確率を計算します。これにはsurvfit()関数を使います。

~(チルダ)という記号を使いますが、これは「何によって分けるか」を指定する記号だと思ってください。

  • 全体の生存曲線: Surv(時間, イベント) ~ 1
  • グループ別(例:性別)の生存曲線: Surv(時間, イベント) ~ グループ変数

今回は、性別(Male vs Female)で生存曲線がどう違うかを見たいので、後者を使います。

R

# 性別で層別化したカプランマイヤー法を実行
fit <- survfit(Surv(time, status) ~ sex, data = lung_cleaned)

# 結果の要約を表示
print(fit)

print(fit)を実行すると、以下のような結果が表示されます。

Call: survfit(formula = Surv(time, status) ~ sex, data = lung_cleaned)

          n events median 0.95LCL 0.95UCL
sex=Male  138    112    270     212     310
sex=Female 90     53    426     348     550

この結果から、

  • 男性 (sex=Male) は138人中112人がイベント(死亡)を経験し、生存期間の中央値(50%の人が生存している期間)は270日。
  • 女性 (sex=Female) は90人中53人がイベントを経験し、生存期間の中央値は426日。

であることが読み取れます。女性の方が生存期間が長い傾向にありそうですね。これを視覚的に確認しましょう。

Step 4: カプランマイヤー曲線の描画

いよいよハイライトです。survminerパッケージのggsurvplot()関数を使って、美しく分かりやすいグラフを作成します。

まずは、最もシンプルな形で描いてみましょう。

R

# 最もシンプルなカプランマイヤー曲線
ggsurvplot(fit, data = lung_cleaned)

これだけでも基本的なグラフは描けますが、ggsurvplot()の真価は豊富なカスタマイズオプションにあります。論文やプレゼンテーションでそのまま使えるような、情報量の多いグラフを作成してみましょう。

R

# 論文掲載レベルのカスタマイズを加えたカプランマイヤー曲線
ggsurvplot(
  fit,
  data = lung_cleaned,
  pval = TRUE,                 # ログランク検定のp値を表示
  conf.int = TRUE,             # 生存曲線の95%信頼区間を表示
  risk.table = TRUE,           # 時間経過毎のリスク集団の数(Number at Risk)を表示
  risk.table.col = "strata",   # リスクテーブルを層別(グループ別)に色分け
  linetype = "strata",         # 層別に線の種類を変更
  surv.median.line = "hv",     # 生存期間中央値の線を追加 (horizontal/vertical)
  ggtheme = theme_bw(),        # グラフのテーマを白黒ベースに変更
  palette = c("#E7B800", "#2E9FDF"), # 曲線の色を指定
  legend.title = "Sex",        # 凡例のタイトルを指定
  legend.labs = c("Male", "Female") # 凡例のラベルを指定
)

このコードを実行すると、非常に洗練されたグラフが出力されるはずです。

  • 実線が各グループの生存曲線を示し、薄い色の帯が95%信頼区間です。
  • グラフ右上にp値が表示されています。これが次のステップで解説するログランク検定の結果です。
  • グラフの下にはRisk Tableが表示されており、各時点で何人の患者がまだ追跡中か(リスクに晒されているか)が一目でわかります。これは論文で非常に好まれる形式です。
  • 点線で示された生存期間中央値も確認できます。

4. ログランク検定で群間比較

グラフを見て、「女性の方が男性より生存率が高そうだ」と感じましたが、それが統計的に意味のある「偶然とは言えない差」なのかを客観的に評価する必要があります。そのために使うのがログランク検定です。

ログランク検定とは?

ログランク検定は、「複数のグループ間で、生存時間の分布に差はない」という帰無仮説を検定する手法です。

簡単に言うと、描画した複数のカプランマイヤー曲線が、全体として統計的に有意に離れているかどうかを評価します。

p値が一般的に用いられる有意水準(例:0.05)よりも小さい場合、「帰無仮説は棄却され、グループ間の生存時間には統計的に有意な差がある」と結論づけることができます。

Rでの実行方法と結果の解釈

実は、先ほどのggsurvplot()pval = TRUEと指定したため、検定は既に実行され、p値がグラフに表示されています。今回の例では、p = 0.0013 のように表示されているはずです。

このp値は0.05よりも十分に小さいため、「男性と女性の生存期間には、統計的に有意な差がある」と結論できます。

検定の詳細な結果を見たい場合は、survdiff()関数を使います。

R

# ログランク検定の実行
surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung_cleaned)

# 検定結果の表示
print(surv_diff)

結果は以下のようになります。

Call:
survdiff(formula = Surv(time, status) ~ sex, data = lung_cleaned)

            n Observed Expected (O-E)^2/E (O-E)^2/V
sex=Male  138      112     91.6      4.55      10.3
sex=Female 90       53     73.4      5.68      10.3

 Chisq= 10.3  on 1 degrees of freedom, p= 0.001

Observedは実際に観測されたイベント数、Expectedはもし男女で差がないとした場合の期待イベント数です。この観測値と期待値のズレが、カイ二乗(Chisq)統計量として計算され、最終的にp= 0.001というp値が算出されています。

まとめ

お疲れ様でした!本記事では、生存時間解析の基本概念から、Rを用いた具体的な実践方法までを解説しました。

  • 生存時間解析は、「時間」と「打ち切り」を考慮したデータ分析手法であること。
  • カプランマイヤー法で生存曲線を推定し、視覚的に評価できること。
  • ログランク検定で、グループ間の生存時間の差を統計的に評価できること。
  • Rのsurvivalsurvminerパッケージを使えば、これらの一連の分析を簡単かつ美しく実行できること。

これだけの知識とスキルがあれば、多くの臨床研究の論文で使われている生存時間解析を理解し、あなた自身のデータで実践する準備が整ったと言えるでしょう。

次のステップとしては…

  • ご自身のデータで、今回学んだコードを試してみましょう。変数の名前などを書き換えるだけで応用できるはずです。
  • 今回は「性別」という一つの変数だけで比較しましたが、「年齢」や「治療法」など、複数の要因が生存期間に影響を与えることも考えられます。そうした多変量の影響を考慮した解析手法がCox比例ハザードモデルです。

今後の記事では、このCox比例ハザードモデルについても解説していく予定です。ぜひ、引き続きリアルワールドエビデンスブログをチェックしてください!

-ビッグデータ, 公衆衛生学, 医学, 統計学

© 2025 RWE