34. 生存分析 (Survival Analysis)

Rを使った分析の準備

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

library(corrplot)
library(jtools)
library(margins)
library(patchwork)
library(prediction)
library(ROCR)
library(survival)
library(survminer)
library(stargazer)
library(tidyverse)

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 = 死亡せず or 打ち切り、1 = 死亡
x	Nonmaintained = 維持化学療法なし、Maitained = 維持化学療法あり

• time が少ない順に並べ替えてみる

• 患者12から患者3まで9人の患者の状況、time と status ごとに整理してみる

• 23 人の患者の中で死亡が判明しているのは 8 名・・・死亡率は \(\frac{8}{23}\) なのか？
  →　患者3が13日目に打ち切りされている
  →　患者3が打ち切り後に「死亡したのか」「死亡しなかったのか」不明
  →　この事実を確認しないまま死亡率を計算すると・・・死亡率は \(\frac{8}{23}\)
  
• しかし、患者3が打ち切り後に「死亡している」場合　→　実際の死亡率は\(\frac{9}{23}\)
  →　死亡率\(\frac{8}{23}\)という結果は、実際の死亡率 \(\frac{9}{23}\) を過小評価してしまうことになる
  
• 打ち切りされたケースを無視すると、死亡率の過小評価を防ぐため
  →　カプランマイヤー曲線の縦軸は「累積の生存率」で計算している

「累積の生存率」とは

• 下の図は、横軸が生存時間 (week)、縦軸が「累積の生存率」を表す
  
• その日ごとの生存時間を掛け合わせる事によって、その日の終わりまでに生存している累積の確率を計算できる
  
• 例えば、初日の生存確率・・・ \(\frac{23}{23}\)= 1
  
• 5日目の生存確率・・・\(\frac{21}{23}\) = 0.91
  
• 5日目までの累積の生存確率・・・ 1×0.91 = 0.91

• リスク人数 (n_risk)
  →　まだ死亡を経験せず、かつまだ観察が打ち切られていない患者数

初日（Time = 0）では誰も死亡していない

→　死亡確率は 0
→　生存確率 = \(\frac{23}{23}\) = 1

5日目には2名の患者（患者12と13）が死亡

→　生き残ったのは 23名中21名
→　生存確率 = \(\frac{21}{23}\) = 0.91

\[累積の生存時間（5日目） = 「生存時間（初日）」× 「生存時間（5日目）」\\ = 1×0.91 = 0.91\]

8日目には2名の患者（患者14と15）が死亡

→　生き残ったのは 21名中19名
→　生存確率 = \(\frac{19}{21}\) = 0.9047　（分母が21になっていることに注意）
→　生存確率 = 1 - 0.095 = 0.905
\[累積の生存時間（8日目） = 「生存時間（初日）」× 「生存時間（5日目）×「生存時間（8日目）\\ = 1 × 0.91 × 0.905 = 0.823\]

• 実際の「累積の生存率」をチェックしてみる

# KM推定（全体）
fit_all <- survfit(Surv(time, status) ~ 1, data = aml)

# サマリーを表示（timeごとの生存確率）
summary(fit_all)

Call: survfit(formula = Surv(time, status) ~ 1, data = aml)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    5     23       2   0.9130  0.0588       0.8049        1.000
    8     21       2   0.8261  0.0790       0.6848        0.996
    9     19       1   0.7826  0.0860       0.6310        0.971
   12     18       1   0.7391  0.0916       0.5798        0.942
   13     17       1   0.6957  0.0959       0.5309        0.912
   18     14       1   0.6460  0.1011       0.4753        0.878
   23     13       2   0.5466  0.1073       0.3721        0.803
   27     11       1   0.4969  0.1084       0.3240        0.762
   30      9       1   0.4417  0.1095       0.2717        0.718
   31      8       1   0.3865  0.1089       0.2225        0.671
   33      7       1   0.3313  0.1064       0.1765        0.622
   34      6       1   0.2761  0.1020       0.1338        0.569
   43      5       1   0.2208  0.0954       0.0947        0.515
   45      4       1   0.1656  0.0860       0.0598        0.458
   48      2       1   0.0828  0.0727       0.0148        0.462

カプラン・マイヤー曲線のポイント ・「打ち切り」によって抜け落ちた人を母数から抜いて生存確率を計算する
→　抜け落ちた人が「もし抜け落ちなかったらどれくらいイベントが起きるか」を計算できる
→　リスクを計算する上での過小評価を回避できる

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

• カプラン＝マイヤー生存曲線 (Kaplan–Meier survival curve)
  
• 生存関数 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

結果解釈のポイント

• 縦軸 (Survival Probability) は「累積の生存率」
  
• それぞれの日に生き延びる生存率をかけ算して「累積の生存率」を計算
  
• KM曲線の線が上にあるほど、死亡せずに生存している割合が高い
  
• 青線（Maintained）が赤線（Non-maintained）より上であれば、維持療法は死亡を遅らせている可能性がある
  
• ログランク検定で p-value < 0.05 なら、治療群間に統計的に有意な差がある

結果の解釈

• x軸は 0（観察開始時）から最後に観察された時点までの時間
  

• y軸は、生存している被験者の割合（累積の生存率）
  

• 時間がゼロの時点では、被験者全員が事象なしで（=1人も死亡せず）生存している
  

• 階段状の実線（青色と赤色）は事象発生の進行を示す
  

• 青色は「維持化学療法あり (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 は慎重に検討する必要あり

1.7 時間を指定しての生存確率の比較

• 長期追跡データではイベント数が減り、末期の推定は不安定になる
  →　解析対象を「臨床的に重要な時点（例えば：40週まで）」に制限することで、安定した比較が可能になる
  
• この場合の解釈は「40週までのフォローアップに限定すれば、群間の生存曲線に統計的差があるか？」

40週までの生存率

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

# --- KM推定（群別） ---
fit <- survfit(Surv(time, status) ~ x, data = aml)

# --- 40週時点の生存確率を取得 ---
s40 <- summary(fit, times = 40)
df40 <- data.frame(
  group = sub("^x=", "", s40$strata),  # "x=Maintained" → "Maintained"
  surv  = s40$surv
)

# ラベル（百分率表示）と色を準備（ggsurvplotのpaletteに合わせる）
df40$label <- paste0(df40$group, ": ", sprintf("%.1f%%", df40$surv * 100))
df40$col   <- ifelse(df40$group == "Maintained", "blue", "red")

# --- プロット作成 ---
p <- ggsurvplot(
  fit, data = aml,
  conf.int = TRUE,
  risk.table = TRUE,
  xlab = "Time to Relapse (weeks)",
  ylab = "Survival Probability",
  legend.title = "Therapy",
  legend.labs = c("Maintained", "Non-maintained"),
  palette = c("blue", "red")
)

# 40週の縦線・ポイント・ラベルを重ねる
p$plot <- p$plot +
  geom_vline(xintercept = 40, linetype = 2) +
  geom_point(
    data = df40,
    aes(x = 40, y = surv),
    color = df40$col, size = 3, inherit.aes = FALSE
  ) +
  geom_text(
    data = df40,
    aes(x = 40, y = surv, label = label),
    vjust = -0.8, color = df40$col, fontface = "bold",
    inherit.aes = FALSE
  )

p

• 上の図は二つの部分から構成されている

上部の曲線 (Survival Probability & Time to Relapse: weeks)

• 40週までに再発せずに生存している確率の推定値
  ・維持療法を受けた患者 Maintained: 0.36 (36.8%)
  ・維持療法を受けない患者 Notmintained: 0.194 (19.4%)

下部のリスクテーブル (Therapy & Number at risk)

• Therapy という因子の水準は2つ：Maintained, Non-maitained
  

• Number at risk: 各時点でまだイベントが起きておらず、死亡せず残っている人数
  

• t = 0 の時点では、療法を受けた患者 (11人)＋受けない患者（12人）= 23人
  

• t = 40 の時点では、療法を受けた患者 (3人)＋受けない患者（2人）= 5人

• survdiff()関数を使って、40週までの log-rank検定結果を確認してみる
  （注：survdiff() はデータフレームを第一引数として受け取る関数ではないので、formula と data = ... を明示的に指定する必要がある）

aml_40 <- aml |>      
  filter(time <= 40) 

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

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

                 N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained     8        6      8.3     0.639       1.6
x=Nonmaintained 10        9      6.7     0.793       1.6

 Chisq= 1.6  on 1 degrees of freedom, p= 0.2 

• 分析結果の表示：

N	サンプル数	→ 40週までのサンプル総数は 18
Observed	死亡数	→ 40週までの死亡総数は 15
(O-E)^2/E	期待値 (Expected)
(O-E)^2/V	カイ二乗検定統計量

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

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

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

• 期待される死亡数は 6.7人だが、実際には 9人も死亡
  
• 想定より多くの死亡（9人）が観察

さて、結果は・・・・

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

1.8 半分が生き残るまでの期間

生存時間中央値 (Median Survival Time)

• 生存確率が50%になる time の値のこと
  ・維持療法を受けた患者 Maintained: 31 weeks
  ・維持療法を受けない患者 Notmintained: 23 weeks
  →　維持療法を受けない患者の方が 7 週間早く生存者が半分になる

library(survival)
library(survminer)
library(ggplot2)

# --- Kaplan–Meier 推定 ---
fit <- survfit(Surv(time, status) ~ x, data = aml)

# --- 中央生存時間を取り出す ---
tb <- summary(fit)$table
df_med <- data.frame(
  group  = gsub("^x=", "", rownames(tb)),   # → "Maintained" / "Nonmaintained"
  median = as.numeric(tb[, "median"]),
  stringsAsFactors = FALSE
)

# NA（中央値が未到達）の群を除外
df_med2 <- subset(df_med, !is.na(median))

# 各群ごとに分割
df_m <- subset(df_med2, group == "Maintained")
df_n <- subset(df_med2, group == "Nonmaintained")   # ←注意：ハイフンなし

# --- KM曲線を描画 ---
p <- ggsurvplot(
  fit, data = aml,
  conf.int = TRUE, risk.table = TRUE,
  xlab = "Time to Relapse (weeks)", ylab = "Survival Probability",
  legend.title = "Med.Time(weeks)",
  legend.labs = c("Maintained", "Non-maintained"),
  palette = c("blue", "red"),
  xlim = c(0, 50)
)

# --- 横線・交点のドット＆数値ラベル ---
p$plot <- p$plot +
  geom_hline(yintercept = 0.5, linetype = 2) +
  # Maintained: ドット＆数値（青）
  geom_point(data = df_m, aes(x = median, y = 0.5),
             inherit.aes = FALSE, color = "blue", size = 3) +
  geom_text(data = df_m, aes(x = median, y = 0.55, label = sprintf("%.0f", median)),
            inherit.aes = FALSE, color = "blue", fontface = "bold", hjust = -0.3) +
  # Nonmaintained: ドット＆数値（赤）
  geom_point(data = df_n, aes(x = median, y = 0.5),
             inherit.aes = FALSE, color = "red", size = 3) +
  geom_text(data = df_n, aes(x = median, y = 0.55, label = sprintf("%.0f", median)),
            inherit.aes = FALSE, color = "red", fontface = "bold", hjust = -0.3)


p

• 「生存時間中央値が群間で有意に異なるかどうか」を調べたいとき
  →　中央値そのものに直接 t 検定のような方法はない
  →　中央値の差 (= 8) に加えて「生存曲線に差があるか」を log-rank test で統計的に検定する

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 

さて、結果は・・・・

半分が生き残るまでの期間差の検定 ・p-value = 0.07 は、有意水準 0.05 をわずかに上回っており、統計的に有意とはいえない
⇒「維持療法を受けた群と受けない群が半分生き残るまでの期間に差がある」と断定できる証拠は十分ではない

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

• 例えば「年齢や性別を考慮しても、維持療法を受けた患者のほうが、死亡が遅れる傾向はあるのか？」を知りたいとする
  
• しかし「KM法では 1 変数（グループ）しか扱えない」
  ⇒　多変量解析ができない

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

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

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

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

変数の変更 (sex) => female

• sex: 性別（1＝男性、2＝女性）
  → sex が 女性なら1、男性なら 0 変更する → female = 1（女性の場合）、female = 0（男性の場合）という変数を作る

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

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

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

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

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

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

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

# KM曲線の描画
km_1 <- ggsurvplot(fit, data = lung,
           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, death) ~ age + 
    female + 
    ph.ecog + 
    ph.karno + 
    pat.karno + 
    meal.cal + 
    wt.loss,
  data = lung)

summary(cox_model)

Call:
coxph(formula = Surv(time, death) ~ 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 → 統計的に有意ではない

3. ロジスティック回帰分析との違い

ロジスティック回帰分析（性別と肺がんの死亡）

library(tidyverse)

• データの読み込み

stargazer::stargazer(lung, 
  type = "html")

Statistic	N	Mean	St. Dev.	Min	Max
inst	227	11.088	8.303	1	33
time	228	305.232	210.646	5	1,022
status	228	1.724	0.448	1	2
age	228	62.447	9.073	39	82
sex	228	1.395	0.490	1	2
ph.ecog	227	0.952	0.718	0	3
ph.karno	227	81.938	12.328	50	100
pat.karno	225	79.956	14.623	30	100
meal.cal	181	928.779	402.175	96	2,600
wt.loss	214	9.832	13.140	-24	68
death	228	0.724	0.448	0	1
female	228	0.395	0.490	0	1

• death と age / death と female の散布図を描く（オッズの対数をとった値で）

• 文字化けを避けるため、Macユーザは次の2行を実行する

theme_set(theme_gray(base_size = 10, 
                     base_family = "HiraginoSans-W3"))

p_age <- ggplot(lung, aes(x = age, y = death)) + 
  geom_jitter(size = 1,
              alpha = 1/3,
              width = 0,
              height = 0.05) +
  geom_smooth(method = "glm", 
    color = "red",
    method.args = list(family = binomial(link = "logit"))) +
  labs(x = "年齢",
       y = "死亡確率")

p_female <- ggplot(lung, aes(x = female, y = death)) + 
  geom_jitter(size = 1,
              alpha = 1/3,
              width = 0.05,
              height = 0.05) +
  geom_smooth(method = "glm", 
    color = "red",
    method.args = list(family = binomial(link = "logit"))) +
  labs(x = "性別",
       y = "死亡確率")

p_age + p_female 

model_lung <- glm(death ~ age + female + ph.ecog + ph.karno + pat.karno + meal.cal + wt.loss, 
            data = lung, 
            family = binomial(link = "logit")) # 係数を「オッズの対数」に指定

• tidy() を使って、推定結果を確認する

broom::tidy(model_lung, conf.int = TRUE)

# A tibble: 8 × 7
  term         estimate std.error statistic p.value  conf.low conf.high
  <chr>           <dbl>     <dbl>     <dbl>   <dbl>     <dbl>     <dbl>
1 (Intercept) -2.28      3.30        -0.692  0.489  -8.80       4.24   
2 age          0.0210    0.0219       0.958  0.338  -0.0221     0.0643 
3 female      -0.940     0.392       -2.40   0.0164 -1.72      -0.177  
4 ph.ecog      1.05      0.480        2.19   0.0287  0.120      2.02   
5 ph.karno     0.0293    0.0271       1.08   0.279  -0.0245     0.0827 
6 pat.karno   -0.0148    0.0161      -0.914  0.361  -0.0471     0.0165 
7 meal.cal     0.000250  0.000502     0.498  0.618  -0.000700   0.00129
8 wt.loss     -0.00534   0.0147      -0.364  0.716  -0.0336     0.0247 

→　表示される回帰式の係数 estimate は「オッズの対数」 = log(odds)
= log-oddsとも呼ばれる

• 例えば、female の係数の推定値は-0.940
  
• これは log-odds = log(odds) = logit
  
• もし、これが重回帰分析の結果であれば、「選挙費用が 1 単位（100 万円）増えるごとに、死亡確率が 0.940 ポイント減る」という意味
  
• しかし、ロジスティック回帰分析の係数は、重回帰分析の係数と同じように解釈できない

log-odds は解釈しにくい → オッズ比に変換する

• model_lung の結果をオッズ比で表示させてみる

tidy(model_lung, 
  conf.int = TRUE, 
  exponentiate = TRUE)

# A tibble: 8 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)    0.102  3.30        -0.692  0.489  0.000151    69.4  
2 age            1.02   0.0219       0.958  0.338  0.978        1.07 
3 female         0.390  0.392       -2.40   0.0164 0.179        0.837
4 ph.ecog        2.86   0.480        2.19   0.0287 1.13         7.51 
5 ph.karno       1.03   0.0271       1.08   0.279  0.976        1.09 
6 pat.karno      0.985  0.0161      -0.914  0.361  0.954        1.02 
7 meal.cal       1.00   0.000502     0.498  0.618  0.999        1.00 
8 wt.loss        0.995  0.0147      -0.364  0.716  0.967        1.02 

Forest Plot による Odds Ratio の表示

# 結果を整形（オッズ比に変換）
results <- tidy(model_lung, conf.int = TRUE, exponentiate = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    term = recode(term,
                  "female" = "女性ダミー",
                  "age" = "年齢",
                  "ph.ecog" = "ECOGパフォーマンススコア",
                  "ph.karno" = "カルノフスキースコア（医師）",
                  "pat.karno" = "カルノフスキースコア（患者）",
                  "meal.cal" = "食事摂取したカロリー量",
                  "wt.loss" = "体重減少（ポンド）"),
    OR_label = sprintf("%.2f", estimate)   # ORを文字列化
  )

# forest plot + 数値ラベル
ggplot(results, aes(x = estimate, y = term)) +
  geom_point(size = 3, color = "blue") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "red") +
  geom_text(aes(label = OR_label), # ORを表示
    hjust = 0.5, 
    vjust = -2, 
    size = 4) +  
  labs(
    title = "Odds Ratios from Logistic Regression",
    x = "Odds Ratio (95% CI)",
    y = ""
  ) +
  theme_minimal(base_size = 14) + # ラベルが切れないように余白
  xlim(0, max(results$conf.high) * 1.2)   +
   theme_bw(base_family = "HiraKakuProN-W3")

female の係数：0.39 （p-value = 0.0164）

• 「女性のオッズ」・・・\(odds_{female = 1}\)　　
• 「男性のオッズ」・・・\(odds_{female = 0}\)

\[OddRatios = \frac{odds_{female=1}}{odds_{female = 0}} = 0.39\]

→　\(odds_{female = 1}\) のオッズが \(odds_{female = 0}\) のオッズの 0.39倍になる
→　女性だと、死亡確率が下がる

統計的有意性

• female の p.value は 0.0164
  →　有意水準を 0.05 とすれば、「性別は死亡確率に影響がない」という帰無仮説を棄却できる
  → female は死亡確率に負の影響があり、その効果は統計的に有意と判断する

log-odds は解釈しにくい → 確率に変換する

3.1 性別ごとの予想死亡確率: predict()

→　「代表的な患者像での性差」を示したいとき
- age, ph.ecog, ph.karno, pat.karno, meal.cal, wt.loss を平均に固定

library(tidyverse)

# 平均 previous を固定して、female を 0 / 1 に
newdata <- data.frame(
  female = c(0, 1),
  age = mean(lung$age, na.rm = TRUE),
  ph.ecog = mean(lung$ph.ecog, na.rm = TRUE),
  ph.karno = mean(lung$ph.karno, na.rm = TRUE),
  pat.karno = mean(lung$pat.karno, na.rm = TRUE),
  meal.cal = mean(lung$meal.cal, na.rm = TRUE),
  wt.loss = mean(lung$wt.loss, na.rm = TRUE)
)

# 予測値と標準誤差
pred <- predict(model_lung, 
  newdata = newdata, 
  type = "response", 
  se.fit = TRUE)

# プロット用データ
plot_df <- newdata %>%
  mutate(
    fit = pred$fit,
    se = pred$se.fit,
    lower = fit - 1.96 * se,
    upper = fit + 1.96 * se,
    sex = factor(female, labels = c("male", "female"))
  )

# グラフ描画：p値を図中にテキストとして追加
plot_prediction <- ggplot(plot_df, aes(x = sex, y = fit)) +
  geom_point(size = 5, color = "black") +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
  geom_text(aes(label = paste0(round(fit * 100, 1), "%")), 
            hjust = -0.5, size = 5) +
  geom_text(
    data = data.frame(x = 1.5, y = max(plot_df$upper) + 0.05),
    aes(x = x, y = y, label = "p = 0.0164"),
    inherit.aes = FALSE,
    size = 5
  ) +
  labs(
    title = "性別による予測死亡確率",
    x = "患者の性別",
    y = "予測死亡確率"
  ) +
   theme_bw(base_family = "HiraKakuProN-W3")

plot_prediction

結果（予測死亡確率の比較）

・図は、性別によって予測される死亡確率の違いを示している
・6つの統制変数をそれぞれの平均値に固定している

・男性の死亡確率は 80.7%
・女性の死亡確率は 61.9%
・p-value = 0.024 なので、5%有意水準で「女性候補は男性候補に比べて死亡確率が低い」が有意
・予測死亡確率の差は約 18.8%ポイント (80.7 − 61.9 = 18.8)
・同じ条件の人と比べると、女性は男性に比べて死亡する確率が有意に低い

• 独立変数がカテゴリカル変数の場合、分析の主要な関心が「性差が死亡確率に与える効果」であれば、結果の提示はこれで十分
  
• しかし「条件によって女性の効果がどう変わるか」が主要関心であれば、female = 0 と female = 1 それぞれの条件付き限界効果を示す必要がある

3.2 性別ごとの予想死亡確率: margins()

Average Marginal Effect (AME) を計算して「性別の効果（差分）」を示す

→　「全体平均での性差」を示したいとき　　

• カテゴリ変数（female）では「男性 → 女性に変わったとき、死亡確率がどれだけ変化するか」を各観測値で計算
  →　その平均をとる →　これが平均限界効果 (AME: Average Marginal Effect)
  →　「平均的に、女性は男性より死亡確率が・・%ポイント高い（低い）」と示せる
  →　係数の意味がさらに明快になる

# 必要パッケージ
library(margins)
library(tidyverse)

# --- p値の整形関数 ---
format_p <- function(x) {
  ifelse(x < 0.001,
         "p < 0.001",
         formatC(x, format = "e", digits = 3))
}
# --- 限界効果の計算 ---
mfx <- margins(model_lung, variables = "female")

summary(mfx)

 factor     AME     SE       z      p   lower   upper
 female -0.1660 0.0650 -2.5543 0.0106 -0.2933 -0.0386

• この AME (-0.1660) が示すこと:
  ⇒ 「男性 → 女性に変わったとき、死亡確率が平均的に 16.6%ポイント (-0.1660) 低下する」
• p-value ≈ 0.0106　→　統計的に有意

結果（限界効果の比較）

「性別の効果（差分）」がわかる

=「性別が変わることで死亡確率に与える平均的な影響」がわかる
・AME = -0.1660 → 女性だと平均的に死亡確率が約 16.6 %ポイント下がる
・男女別に「その効果が統計的に有意」(p-value = 0.0106)

結果のまとめ

・3.1の結果: predict()

・「平均的な患者 1 人」を仮定して男女の予測確率を比較
→　「女性は男性に比べて死亡確率が 18.8%ポイント低い」

・3.2の結果: margins()

・サンプル各観測の条件を尊重して、女性の時の確率変化の平均値を示す（AME、95%CI付き）
→　「女性だと平均的に死亡確率が約 16.6 %ポイント低い」

ロジスティック回帰分析は非線形なので、数値が一致しないのが普通

3.3 妥当性の検討

Cox比例ハザードモデル

目的	「生存時間」に基づいて死亡リスク（ハザード比）を推定
結果	女性は男性に比べて死亡リスクが約 42% 低い、ph.ecog や ph.karnoが有意
強み	死亡確率そのものではなく、ある時点まで生き残る確率（生存曲線）がわかる
妥当性	生存データ（time + censoring）があるなら、こちらが基本的に最も妥当

ロジスティック回帰の予測値プロット

目的	死亡を0/1の二値として扱い、説明変数から「死亡確率」を推定
結果	女性の死亡確率・・・61.9%、男性の死亡確率・・・80.7%
弱み	生存期間中に観察終了した人 (status = 1) を「死亡していない (death = 0)」と扱ってしまう
妥当性	生存データを無視しているため、Cox 比例ハザードよりは劣る

• 次のデータに status = 1、death = 0 と入力
  →　生存期間中に観察終了した人（status = 1）で「死亡していない (death = 0)」人が 63 人いることがわかる

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

どちらが妥当か？

• 生存時間データ（time, censoring）があるなら
  → Cox比例ハザードモデルがより妥当
  
• ロジスティック回帰は「ある時点までに死亡したか否か」しか見ていない
  →　時間情報を捨てている
  
• ただし、論文や発表で「分かりやすく男女差を示す図」を作るには、ロジスティック回帰による予測確率グラフも直感的で有用

参考文献

• Box-Steffensmeier, Janet M., and Bradford S. Jones (2004). Event History Modeling: A Guide for Social Scientists. Cambridge University Press.
• 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