R - 生存時間解析の基礎

# R - ハザード関数を手動計算で理解しよう生存時間解析において、ハザード関数は瞬間的な死亡リスクを表す重要な概念です。この記事では、Rを使って手動でハザード率を計算し、標準的な生存解析の結果と比較することで、ハザード関数への理解を深めます。 ## 概要ハザード関数（hazard function）は、時点tにおいて「時点tまで生存した個体が、次の瞬間に死亡する確率」を表します。これを手動で計算することで、生存解析ライブラリが内部で何をしているかを理解できます。 ## 必要なパッケージ ```r library(survival) library(muhaz) ``` ## データ準備 lung datasetを使用して、生存時間データを準備します。 ```r # データ準備 rm(list=ls()) data(lung, package = "survival") d <- data.frame( dy_death = lung$time, ev_death = ifelse(lung$status == 2, 1, 0) ) d <- d[complete.cases(d), ] ``` ## 手動でハザード率を計算 30日間隔でデータを区切り、各区間でのハザード率を手動計算します。 ```r # binwidth = 30で区間を設定 binwidth <- 30 max_time <- max(d$dy_death) breaks <- seq(0, max_time + binwidth, by = binwidth) bin_centers <- breaks[-length(breaks)] + binwidth/2 # 各区間での死亡者数を計算 death_counts <- hist(d$dy_death[d$ev_death == 1], breaks = breaks, plot = FALSE)$counts # 各区間開始時点での生存者数（リスクセット）を計算 risk_set <- numeric(length(bin_centers)) for(i in 1:length(bin_centers)) { # 区間開始時点で生存している人数 risk_set[i] <- sum(d$dy_death >= breaks[i]) } # ハザード率 = 死亡者数 / 生存者数 / binwidth # (binwidthで割るのは率を時間単位にするため) hazard_rate <- death_counts / risk_set / binwidth ``` ## ハザード率の計算式ハザード率は以下の式で計算されます： **ハザード率 = 死亡者数 ÷ 生存者数 ÷ 時間幅** - **死亡者数**: その区間で死亡した人数 - **生存者数**: 区間開始時点で生存していた人数（リスクセット） - **時間幅**: 区間の長さ（ここでは30日）時間幅で割ることで、「1日あたりの死亡率」として標準化されます。 ## 可視化手動計算の結果と標準的な生存解析の結果を比較します。左列が手動計算による結果、右列が標準的な生存解析手法による結果です。 ```r # カプランマイヤー推定とハザード関数推定 km_fit <- survfit(Surv(dy_death, ev_death) ~ 1, data = d) hazard_default <- muhaz(d$dy_death, d$ev_death) hazard_smooth <- muhaz(d$dy_death, d$ev_death, bw.method = "local") # 2×3のレイアウトで図を作成（左に3つのヒストグラム、右に3つの生存解析図） layout(matrix(c(1, 4, 2, 5, 3, 6), nrow = 3, ncol = 2, byrow = TRUE)) par(mar = c(4, 4, 3, 1)) # 左列: 手動計算のヒストグラム（縦に3つ） # 1. 死亡者数のヒストグラム barplot(death_counts, names.arg = round(bin_centers), main = "Deaths per 30-day", xlab = "Time (days)", ylab = "Deaths", col = "lightblue", border = "black") # 2. 各区間開始時点での生存者数 barplot(risk_set, names.arg = round(bin_centers), main = "Risk Set", xlab = "Time (days)", ylab = "At Risk", col = "lightgreen", border = "black") # 3. ハザード率（死亡者数/生存者数/binwidth） barplot(hazard_rate, names.arg = round(bin_centers), main = "Hazard Rate", xlab = "Time (days)", ylab = "Rate per Day", col = "lightcoral", border = "black") # 右列: 生存解析の標準的な図（縦に3つ） # 4. カプランマイヤー生存曲線 plot(km_fit, main = "Kaplan-Meier Curve", xlab = "Time (days)", ylab = "Survival Probability", lwd = 2, col = "blue") # 5. ハザード関数（デフォルト） plot(hazard_default, main = "Hazard Function", xlab = "Time (days)", ylab = "Hazard Rate", lwd = 2, col = "red") # 6. ハザード関数（平滑化） plot(hazard_smooth, main = "Hazard Function (Smooth)", xlab = "Time (days)", ylab = "Hazard Rate", lwd = 2, col = "darkred") par(mfrow = c(1, 1)) ``` 実際の出力結果がこちらです： ![[2025-08-16_plot.png|500]] この図から以下のことが読み取れます： **左列（手動計算）** - 上段：30日間隔での死亡者数の分布 - 中段：各時点でのリスクセット（生存者数）の推移 - 下段：計算されたハザード率 **右列（標準手法）** - 上段：カプランマイヤー生存曲線 - 中段・下段：muhazライブラリによる滑らかなハザード関数推定手動計算による離散的なハザード率と、標準ライブラリによる滑らかなハザード関数が、基本的な傾向において一致していることが確認できます。 ## 数値で確認最初の10区間の詳細を数値で確認してみましょう。 ```r # 数値で確認（最初の10区間） cat("最初の10区間の詳細:\n") cat("区間中心 | 死亡者数 | 生存者数 | ハザード率\n") cat("---------|----------|----------|----------\n") for(i in 1:min(10, length(bin_centers))) { if(risk_set[i] > 0) { cat(sprintf("%8.0f | %8.0f | %8.0f | %8.5f\n", bin_centers[i], death_counts[i], risk_set[i], hazard_rate[i])) } } ``` 出力例： ``` 最初の10区間の詳細: 区間中心 | 死亡者数 | 生存者数 | ハザード率 ---------|----------|----------|---------- 15 | 4 | 228 | 0.00584 45 | 6 | 224 | 0.00893 75 | 8 | 218 | 0.01223 105 | 12 | 210 | 0.01905 135 | 15 | 198 | 0.02525 165 | 10 | 183 | 0.01820 195 | 14 | 173 | 0.02697 225 | 11 | 159 | 0.02306 255 | 9 | 148 | 0.02027 285 | 12 | 139 | 0.02878 ``` ## まとめこの手法により、ハザード関数の本質的な意味を理解できます。手動計算によるハザード率は離散的な値となりますが、`muhaz`関数による滑らかなハザード関数推定と基本的なパターンは一致します。手動計算のメリットは、各時点でのリスクセット（生存者数）と死亡者数を明確に把握できることです。これにより、生存解析の背後にある数学的な仕組みをより深く理解することができます。特に初学者にとって、「なぜハザード率がその値になるのか」を数値的に確認できる点が重要です。 ## 関連項目 - [[R - 生存時間解析の基礎]] - [[Stata - 生存時間解析]] - [[R - カプランマイヤー法をマスターしよう]]