臨床研究やリアルワールドデータ分析において、「新薬Aは従来薬Bよりも患者さんの生存期間を延ばすのか?」「ある遺伝子変異があると、病気の再発までの時間は短くなるのか?」といった問いに答えることは非常に重要です。このような「イベント(死亡、再発など)が発生するまでの時間」を分析する手法を生存時間解析と呼びます。
生存時間解析は非常に強力なツールですが、「統計学は苦手…」「プログラミングなんてやったことがない」と感じる方も多いかもしれません。
しかし、ご安心ください。本記事では、無料で利用できる統計解析ソフトウェア「R」を使い、生存時間解析の最も基本的な手法であるカプランマイヤー(Kaplan-Meier)法とログランク(Log-rank)検定を、誰でも実践できるようにステップ・バイ・ステップで徹底解説します。
この記事を最後まで読めば、以下のことができるようになります。
- 生存時間解析の基本的な考え方を理解する。
- RとRStudioを準備し、解析に必要な環境を整える。
- サンプルデータを使って、カプランマイヤー曲線を描画し、論文に載せられるレベルの綺麗なグラフを作成する。
- ログランク検定を行い、グループ間の生存時間に統計的な差があるか判断する。
少し読み応えがありますが、医療データ分析の第一歩を踏み出していきましょう。
Table of Contents
1. 生存時間解析とカプランマイヤー曲線とは?
まず、具体的な操作に入る前に、基本的な概念を少しだけおさらいしましょう。
生存時間データとは?
通常のデータと少し違う、生存時間データが持つ2つの重要な要素が「時間」と「打ち切り」です。
- イベント発生までの時間 (Time-to-event):研究の観察開始時点(例:治療開始日、診断日)から、特定のイベント(例:死亡、病気の再発、心筋梗塞の発症)が起こるまでの時間のことです。これを「生存時間」と呼びます。
- 打ち切り (Censoring):これが生存時間解析を少し複雑に、そして面白くする概念です。研究の観察期間中に、イベントが発生しなかった、あるいはイベントが確認できなかったケースを「打ち切り」と呼びます。
例えば、以下のようなケースが該当します。
- 研究期間が終了した時点で、まだ生存している。
- 引っ越しや自己都合で、追跡(フォローアップ)が不可能になった。
- 交通事故など、研究対象のイベントとは異なる原因で死亡した。
「打ち切り」になったからといって、そのデータを無視してはいけません。例えば、5年間観察してイベントが起こらなかった患者さんは、「少なくとも5年間は生存していた」という非常に重要な情報を持っています。生存時間解析では、この「打ち切り」データを適切に扱うことで、より正確な分析が可能になります。
カプランマイヤー法とは?
カプランマイヤー法は、この「打ち切り」データを考慮しながら、時間経過とともにイベントを発生せずに生存している個体の割合(生存確率)がどのように変化していくかを推定する手法です。
その結果をグラフにしたものがカプランマイヤー曲線です。
- 縦軸は生存確率(1.0 = 100% から始まり、時間とともに低下)。
- 横軸は時間(日数、月数、年数など)。
曲線はイベント(例:死亡)が発生するたびに階段状にガクンと下がります。「打ち切り」があった時点では、生存している人数は減りますが、生存確率は変わりません(曲線は下がりません)。この曲線を見ることで、「治療開始から3年後の生存率は約80%だな」といった情報を視覚的に読み取ることができます。
2. Rの準備をしよう
それでは、いよいよ実践の準備です。今回は、無料で高機能な統計ソフトウェア「R」と、それを快適に使うための統合開発環境「RStudio」を使用します。
なぜRを使うのか?
- 無料: 何より費用がかかりません。
- 高機能: 世界中の研究者が開発した豊富なパッケージ(拡張機能)により、最新の解析手法もすぐに利用できます。
- 可視化に強い:
ggplot2
などのパッケージを使えば、論文掲載レベルの美しいグラフを自在に作成できます。 - 再現性: コードを保存しておけば、誰でも、いつでも同じ解析を再現できます。
RとRStudioのインストール
まだインストールしていない方は、以下の公式サイトからダウンロードしてインストールしてください。手順は多くの解説サイトで紹介されているので、ここでは割愛します。必ずRを先にインストールしてから、RStudioをインストールしてください。
- Rのダウンロード: The Comprehensive R Archive Network (CRAN)
- RStudioのダウンロード: Posit | RStudio a new name for RStudio (「RStudio Desktop」のFree版をダウンロード)
解析に必要なパッケージのインストール
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の
survival
とsurvminer
パッケージを使えば、これらの一連の分析を簡単かつ美しく実行できること。
これだけの知識とスキルがあれば、多くの臨床研究の論文で使われている生存時間解析を理解し、あなた自身のデータで実践する準備が整ったと言えるでしょう。
次のステップとしては…
- ご自身のデータで、今回学んだコードを試してみましょう。変数の名前などを書き換えるだけで応用できるはずです。
- 今回は「性別」という一つの変数だけで比較しましたが、「年齢」や「治療法」など、複数の要因が生存期間に影響を与えることも考えられます。そうした多変量の影響を考慮した解析手法がCox比例ハザードモデルです。
今後の記事では、このCox比例ハザードモデルについても解説していく予定です。ぜひ、引き続きリアルワールドエビデンスブログをチェックしてください!