Rを使った分析の準備

ここで使うRのパッケージは次のとおり。

library(tidyverse)
library(survival)
library(survminer)

1. 生存分析 (survival analysis) とは

生存時間解析とも呼ばれる
生物の死や機械システムの故障など、1つの事象（event）が発生するまでの予想される期間を分析する統計学の一分野
分野ごとに呼び名が異なる

工学	信頼性分析
経済学	継続時間分析
社会学	イベント履歴分析

1.1 生存分析でわかること

ある時間を過ぎて生存する人々の割合はどのくらいか
生き残った人々のうち、彼らはどのくらいの割合で死亡（または故障するか）
複数の死因または故障を考慮に入れることができるか
特定の状況または特性は、生存の確率をどのように増加（または減少）させるのか

1.2 目的に応じた生存分析の手法

医療を事例として

知りたいこと	適切な手法	概要
治療期間の全体傾向	カプラン＝マイヤー（KM）生存曲線	生存曲線を用いて、治療終了までの患者数の推移（生存率）を視覚化。
男女間での死亡率の差	ログランク検定	2つ以上の群で生存曲線が統計的に異なるかを検定
年齢などが死因に与える影響	Cox比例ハザードモデル	ハザード（死亡の危険度）に影響を与える要因を同時に分析可能。
Coxモデルの各説明変数の有意性	尤度比検定・ワルド検定	Cox回帰の係数の有意性を判定。尤度比検定はモデル全体、ワルド検定は各変数ごとの影響を見る。

1.3 生存分析で使う一般的な用語

事象（Event）	死亡、疾患の発生、疾患の再発、回復、またはその他の興味ある経験
時間（Time）	観察期間の開始（手術や治療の開始など）から、(i) 事象の発生、または (ii) 試験の終了、または (iii) 連絡が途絶えたり研究から離脱するまでの時間
打ち切り (Censoring)	個人の生存時間に関するいくらかの情報を持っている時、生存時間が正確にわからない場合に起こる
生存関数: S(t)	ある被験者が時間 t より長く生存する確率

1.4 データ

ここでは R に付属している survialパッケージに入っている「急性骨髄性白血病生存データセット」を使う

library(survival)

survival パッケージに含まれている aml.csvを使う
これは急性骨髄性白血病生存データセット: Acute Myelogenous Leukemia survival data
aml.csv データを読み込む

aml <- read_csv("data/aml.csv")

aml データに含まれている変数を確認

DT::datatable(aml)

time	生存時間または打ち切り時間 (weeks)
status	がんの再発：0 = 事象なし（打ち切り）、1 = 事象あり（再発）
x	治療群：Nonmaintained（維持化学療法なし）、Maitained（あり）

12番の患者（`time = 5`、`status = 1`) の意味

最も生存期間が短い患者（5週間）
この患者の治療は5週目まで続き (time = 5)
この患者はがんが再発した (status = 1)

11番の患者（`time = 161`、`status = 0`) の意味

最も生存期間が長い患者（161週間）
この患者の治療は 161週目まで続き (time = 161)
この患者の治療は打ち切られている (status = 0)
治療が打ち切られた = 11番の患者に事象がなかった（aml癌の再発がなかった）

2番の患者（`time = 13`、`status = 1`) の意味

この患者の治療は 13週目まで続き (time = 13)
この患者はがんが再発した (status = 1)

3番の患者（`time = 13`、`status = 0`) の意味

この患者の治療は 13週目まで続き (time = 13)
この患者の治療は打ち切られている (status = 0)
この患者は「13週目まで」がんの再発はなかった
この患者は「13週目まで」しか研究に参加していない
研究が終了する13週間前に遅れて登録したため、13週しか観察ででなかったかもしれない
もしくは、この患者は研究の初期に登録されたが、追跡調査を受けなかったか、研究を辞退したのかも

◎注意：次のコマンドを使うと、複数の変数 (time, status, x) の値を同時に指定して確認できる

reactable::reactable(aml,
  filterable = TRUE,  # 検索可能に設定
  defaultPageSize = 23) # 表示行数の限界を指定

ここで知りたいこと

維持療法を受けた患者が、維持療法を受けていない患者に比べて再発が遅くなるかどうか

1.5 全体のカプラン=マイヤー生存曲線

生存関数 S(t) は、被験者が時間 t よりも長く生存する確率
S(t) は、理論的には滑らかな曲線であるが、通常はカプラン=マイヤー(KM)曲線を用いて推定する

# パッケージ読み込み
library(survival)
library(survminer)

# データ読み込み
aml <- read_csv("data/aml.csv")  # アップロードしたaml.csvを読み込む

# 生存オブジェクト作成
# 再発までの時間: time、イベント発生: status、治療: x（Maintained or Nonmaintained）
fit <- survfit(Surv(time, status) ~ x, data = aml)

# KM曲線の描画
km_1 <- ggsurvplot(fit, data = aml,
           pval = TRUE,              # p値（log-rank test）
           conf.int = TRUE,          # 信頼区間を表示
           risk.table = TRUE,        # リスクテーブルを表示
           xlab = "Time to Relapse",
           ylab = "Survival Probability",
           legend.title = "Therapy",
           legend.labs = c("Maintained", "Non-maintained"),
           palette = c("blue", "red"),
           ggtheme = theme_minimal())
km_1

結果解釈のポイント

KM曲線の線が上にあるほど、再発せずに生存している割合が高い
青線（Maintained）が赤線（Non-maintained）より上であれば、維持療法は再発を遅らせている可能性がある
ログランク検定で p-value < 0.05 なら、治療群間に統計的に有意な差がある

結果の解釈

x軸は 0（観察開始時）から最後に観察された時点までの時間
y軸は、生存している被験者の割合
時間がゼロの時点では、100%の被験者が事象なしで（つまり、一人もがんが再発せず）生存している
階段状の実線（青色と赤色）は事象発生の進行を示す
青色は「維持化学療法あり (Maintained)」の場合
赤色は「維持化学療法なし (Nonaintained)」の場合
垂直方向の落ち込みは事象が発生したことを示す
5週目に2人、事象が発生　⇒　最初の「落ち込み」↓
8週目に2人、事象が発生 ⇒　2回目の「落ち込み」↓
9週目に1人、事象が発生 ⇒　3回目の「落ち込み」↓
横線と交差する「＋」・・・この時点で患者が打ち切られたことを示す
13週目の「+」・・・維持化学療法ありの患者が打ち切り
16週目の「+」・・・維持化学療法なしの患者が打ち切り
時間が経過するにつれて、グラフが右下がりになる
⇒　つまり、生存確率が下がる
⇒　維持化学療法なし (Nomaintained)の方が生存確率が大きく下がる
⇒　維持化学療法あり(Maintained)の方がより生存確率が高い
を示している。amlデータ表では、5人の被験者がそれぞれ13、16、28、45、161週目で打ち切られた。KM生存曲線には、これらの打ち切られた観察に対応する5つの目盛り線がある。

1.6 ログランク検定

ログランク検定 (log-rank test) とは「生存率の差の検定」のこと
2 つ以上のグループの生存期間を比較する
ここでは維持療法群（Maintained）と非維持療法群（Non-maintained）での生存率の差についてのログランク検定を使用
ログランク検定の帰無仮説・・・両治療群の生存率が同じ
それぞれの各時点で生存している被験者の期待数を、各事象の時間に治療群内で危険（risk）を抱えている被験者の数に合わせて調整
各治療群で観察された「事象数」(Observed) が「期待数」(Expected) と有意に異なるかどうかを判定する
検定は、カイ二乗分布に基づいてなされる
ログランク検定統計量が大きければ、治療群間の生存期間に差があることの証拠
ログランク検定統計量は、自由度が 1 のカイ二乗分布に近似し、p値はカイ二乗分布を使用して計算される
上で示したカプラン=マイヤー生存曲線には p-value = 0.065 とログランクの検定結果が示されている
⇒　帰無仮説は棄却できない
⇒　治療群間に統計的に有意な差があるとはいえない
念のため、survdiff()関数を使って、log-rank検定結果をもう少し詳しく確認してみる

survdiff(Surv(time, status) ~ x, data = aml)

Call:
survdiff(formula = Surv(time, status) ~ x, data = aml)

                 N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained    11        7    10.69      1.27       3.4
x=Nonmaintained 12       11     7.31      1.86       3.4

 Chisq= 3.4  on 1 degrees of freedom, p= 0.07

分析結果の表示：

N	サンプル数
Observed	再発数
(O-E)^2/E	期待値 (Expected)
(O-E)^2/V	カイ二乗検定統計量

x = Maintained 群（維持療法あり）

期待される再発数は10.69人だが、実際には7人だけ再発
→ 想定より少ない再発（7人）が観察
⇒　治療の効果がある可能性あり

x = Nonmaintained 群（維持療法なし）

期待される再発数は 7.31人だが、実際には 11人も再発
想定より多くの再発（11人）が観察

さて、結果は・・・・

ログランク検定結果・p-value = 0.07 は、有意水準 0.05 をわずかに上回っており、統計的に有意とはいえない
⇒　「維持療法を受けた患者のほうが、再発が遅れる」という傾向はある程度あるが、統計的に確実とは言えない

注意:

被験者 23 人というサンプルサイズは少なすぎる
⇒　治療群間の差を検出する力はほとんどない
⇒　カイ二乗検定は漸近近似法に基づいている
⇒　サンプルサイズが小さい場合は p-value は慎重に検討する必要あり
例えば、ここで「年齢や性別を考慮しても、維持療法を受けた患者のほうが、再発が遅れる傾向はあるのか？」を知りたいとする
しかし「KM法では 1 変数（グループ）しか扱えない」
⇒　多変量解析ができない

解決策　⇒　Cox比例ハザードモデルを使う

KM法とCox比例ハザードモデルの違い

状況	なぜKMでは不十分か？	Coxモデルを使う理由
🔸 複数の要因（共変量）を同時に扱いたい	KM法では 1 変数しか扱えない	Coxモデルは多変量解析が可能
🔸 量的変数を含めたい	KM法はカテゴリ変数しか扱えない	Coxモデルは連続変数も扱える
🔸 交絡要因を統制したい（性別や年齢など）	KM法は統制ができない	Coxモデルで交絡を調整できる
🔸 ハザード比（リスクの大きさ）を知りたい	KM法は「生存率」のみ	Coxモデルはハザード比を出力
🔸 予測モデルを作りたい	KM法は予測に使いにくい	Coxモデルは予測に応用可能

2. Cox比例ハザードモデルを使った分析

2.1 データの準備

ここでは R に付属している survialパッケージに入っている「肺がんデータセット」を使う

library(survival)

survival パッケージに含まれているデータを確認
次のいずれか1つのコマンドで確認する

data(package = "survival") 
help(package = "survival")

この中の lung（肺がんデータ: NCCTG Lung Cancer Data）を使う
lung データを読み込む

names(lung)

 [1] "inst"      "time"      "status"    "age"       "sex"       "ph.ecog"  
 [7] "ph.karno"  "pat.karno" "meal.cal"  "wt.loss"

inst	: 医療機関コード
time	: 生存期間（日数）
status	: 打ち切りステータス（1＝打ち切り、2＝死亡）
age	: 年齢（歳）
sex	: 性別（1＝男性、2＝女性）
ph.ecog	: 医師によるECOGパフォーマンススコア（0＝無症状、1＝症状ありだが完全に歩行可能、2＝1日のうち半分未満をベッドで過ごす、3＝1日のうち半分以上をベッドで過ごすが完全に寝たきりではない、4＝完全に寝たきり）
ph.karno	: 医師によるカルノフスキー・パフォーマンススコア（0＝最悪、100＝最良）
pat.karno	: 患者自身によるカルノフスキー・パフォーマンススコア
meal.cal	:食事から摂取したカロリー量
wt.loss	: 過去6か月間の体重減少（ポンド）

変数の変更 (status)

status: 打ち切りステータス（1＝打ち切り、2＝死亡）
→ 生存分析では event 発生 = 1 （つまり死亡）と指定
→ status が 2（死亡の場合）が 1、それ以外は 0 と変更する

変数の変更 (sex)

sex: 性別（1＝男性、2＝女性）
→ sex が女性なら1、男性なら 0 変更する

lung <- lung |> 
  mutate(status = if_else(status == 2, 1, 0)) |> 
  mutate(female = if_else(sex == 2, 1, 0))

記述統計を表示

DT::datatable(lung)

2.2 カプラン＝マイヤー生存曲線（男女別）

# パッケージ読み込み
library(survival)
library(survminer)

# データ読み込み
aml <- read_csv("data/lung.csv")  # アップロードしたlung.csvを読み込む

# 生存オブジェクト作成
# 再発までの時間: time、イベント発生: status、治療: x（Maintained or Nonmaintained）
fit <- survfit(Surv(time, status) ~ sex, data = lung)

# KM曲線の描画
km_1 <- ggsurvplot(fit, data = aml,
           pval = TRUE,              # p値（log-rank test）
           conf.int = TRUE,          # 信頼区間を表示
           risk.table = TRUE,        # リスクテーブルを表示
           xlab = "Time to Relapse",
           ylab = "Survival Probability",
           legend.title = "Therapy",
           legend.labs = c("Male Patient", "Female Patient"),
           palette = c("blue", "red"),
           ggtheme = theme_minimal())
km_1

グラフからわかること
・女性の曲線が上にある
⇒　男性グループと比較すると、女性グループは、長く生存している傾向あり
・ログランク検定（p値=0.0013）
⇒　日数と生存率に関する男女差は統計的に有意

2.4 Cox比例ハザードモデル（男女別）

ここで知りたいこと

性別が肺がんの生存期間に対してどのような影響を与えているか？
その際、患者の年齢、医師による患者へのパフォーマンススコア、接種カロリー量、過去6ヶ月の体重減少を統制したい
=> カプラン=マイヤー生存曲線ではなく、Cox比例ハザードモデルを使う

cox_model <- coxph(Surv(time, status) ~ age + 
    female + 
    ph.ecog + 
    ph.karno + 
    pat.karno + 
    meal.cal + 
    wt.loss,
  data = lung)

summary(cox_model)

Call:
coxph(formula = Surv(time, status) ~ age + female + ph.ecog + 
    ph.karno + pat.karno + meal.cal + wt.loss, data = lung)

  n= 168, number of events= 121 
   (60 observations deleted due to missingness)

                coef  exp(coef)   se(coef)      z Pr(>|z|)   
age        1.065e-02  1.011e+00  1.161e-02  0.917  0.35906   
female    -5.509e-01  5.765e-01  2.008e-01 -2.743  0.00609 **
ph.ecog    7.342e-01  2.084e+00  2.233e-01  3.288  0.00101 **
ph.karno   2.246e-02  1.023e+00  1.124e-02  1.998  0.04574 * 
pat.karno -1.242e-02  9.877e-01  8.054e-03 -1.542  0.12316   
meal.cal   3.329e-05  1.000e+00  2.595e-04  0.128  0.89791   
wt.loss   -1.433e-02  9.858e-01  7.771e-03 -1.844  0.06518 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

          exp(coef) exp(-coef) lower .95 upper .95
age          1.0107     0.9894    0.9880    1.0340
female       0.5765     1.7347    0.3889    0.8545
ph.ecog      2.0838     0.4799    1.3452    3.2277
ph.karno     1.0227     0.9778    1.0004    1.0455
pat.karno    0.9877     1.0125    0.9722    1.0034
meal.cal     1.0000     1.0000    0.9995    1.0005
wt.loss      0.9858     1.0144    0.9709    1.0009

Concordance= 0.651  (se = 0.029 )
Likelihood ratio test= 28.33  on 7 df,   p=2e-04
Wald test            = 27.58  on 7 df,   p=3e-04
Score (logrank) test = 28.41  on 7 df,   p=2e-04

上部の結果・・・logハザード比　⇒　そのまま解釈できない
下部の結果・・・ハザード比　⇒　解釈可能

female の結果（ハザード比）

指標	値	解釈
female の coef	-5.509e-01	logハザード比（女性であることの効果）
female の exp(coef)	0.5765	ハザード比（解釈可能）
female の p-value	0.00609	有意水準5%未満なので、統計的に有意な差あり

Cox比例ハザードモデルの結果からわかること (female)

ハザード比の結果

女性であること（female = 1）が死亡リスクに与える影響

・exp(coef) = 0.5765 → 女性の死亡リスクは男性の約57.6%
・女性は男性と比較して、死亡リスクが約42.4%低い
・ハザード比が 1 未満 → リスクが低い
・p-value = 0.00609 → 統計的に有意

age の結果（ハザード比）

指標	値	解釈
age の coef(coef)	1.065e-02	logハザード比（年齢の効果）
age の exp(coef)	1.0107	ハザード比（解釈可能）
age の p-value	0.35906	有意水準5%以上なので、統計的に有意ではない

Cox比例ハザードモデルの結果からわかること (age)

ハザード比の結果

年齢が死亡リスクに与える影響

・exp(coef) = 1.0107 →　年齢が1歳高くなるごとに、死亡リスクは約 1.07% 増加すると推定
・しかし、この差は統計的に有意ではない（p-value = 0.359）　　
・従って、他の変数を統制したこのモデルでは「年齢が高いほど死亡リスクが高い」とは断言できない
・ハザード比が 1 未満 → リスクが低い
・p-value = 0.35906 → 統計的に有意ではない

参考文献

Kieran Healy, DATA VISUALIZATION, Princeton, 2019
宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」
浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年
浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年
Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.
Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017

34. 生存分析 (Survival Analysis)

Masahiko Asano

2025-08-08

1. 生存分析 (survival analysis) とは

1.1 生存分析でわかること

1.2 目的に応じた生存分析の手法

医療を事例として

1.3 生存分析で使う一般的な用語

1.4 データ

12番の患者（`time = 5`、`status = 1`) の意味

11番の患者（`time = 161`、`status = 0`) の意味

2番の患者（`time = 13`、`status = 1`) の意味

3番の患者（`time = 13`、`status = 0`) の意味

ここで知りたいこと

1.5 全体のカプラン=マイヤー生存曲線

結果解釈のポイント

結果の解釈

1.6 ログランク検定

x = Maintained 群（維持療法あり）

x = Nonmaintained 群（維持療法なし）

さて、結果は・・・・

解決策　⇒　Cox比例ハザードモデルを使う

KM法とCox比例ハザードモデルの違い

2. Cox比例ハザードモデルを使った分析

2.1 データの準備

変数の変更 (status)

変数の変更 (sex)

2.2 カプラン＝マイヤー生存曲線（男女別）

2.4 Cox比例ハザードモデル（男女別）

ここで知りたいこと

female の結果（ハザード比）

ハザード比の結果

女性であること（female = 1）が死亡リスクに与える影響

age の結果（ハザード比）

ハザード比の結果

年齢が死亡リスクに与える影響

34. 生存分析 (Survival Analysis)

Masahiko Asano

2025-08-08

1. 生存分析 (survival analysis) とは

1.1 生存分析でわかること

1.2 目的に応じた生存分析の手法

医療を事例として

1.3 生存分析で使う一般的な用語

1.4 データ

12番の患者（time = 5、status = 1) の意味

11番の患者（time = 161、status = 0) の意味

2番の患者（time = 13、status = 1) の意味

3番の患者（time = 13、status = 0) の意味

ここで知りたいこと

1.5 全体のカプラン=マイヤー生存曲線

結果解釈のポイント

結果の解釈

1.6 ログランク検定

x = Maintained 群（維持療法あり）

x = Nonmaintained 群（維持療法なし）

さて、結果は・・・・

解決策 ⇒ Cox比例ハザードモデルを使う

KM法とCox比例ハザードモデルの違い

2. Cox比例ハザードモデルを使った分析

2.1 データの準備

変数の変更 (status)

変数の変更 (sex)

2.2 カプラン＝マイヤー生存曲線（男女別）

2.4 Cox比例ハザードモデル（男女別）

ここで知りたいこと

female の結果（ハザード比）

ハザード比の結果

女性であること（female = 1）が死亡リスクに与える影響

age の結果（ハザード比）

ハザード比の結果

年齢が死亡リスクに与える影響

12番の患者（`time = 5`、`status = 1`) の意味

11番の患者（`time = 161`、`status = 0`) の意味

2番の患者（`time = 13`、`status = 1`) の意味

3番の患者（`time = 13`、`status = 0`) の意味

解決策　⇒　Cox比例ハザードモデルを使う