最近では、さまざまな統計解析がパッケージや生成AIを使うことで手軽に実行できるようになりました。そのため、「とにかく結果が出れば良い」という考えが強くなり、解析手法の本質を十分に理解せずに使われることも多いのではないでしょうか。
私自身は、理解すること自体を楽しみたいので、理論的な背景についても学ぶようにしています。しかし、本や論文を完全に読みこなすのは容易ではありません。そんなとき、理解を諦めるのではなく、少しでも近づくための方法を試します。その一つが「数値実験」です。仮定したモデルに基づいて数値データを生成し、それを分析することで、最初は難解に見えた数式の意味が徐々に明確になってくることがあります。
そこで今回は、生存時間解析の理解を深めるために、生存時間解析用の数値データを実際に生成してみたいと思います。
以下では、3 つの代表的なケースに分けて生存解析用の数値データをRで生成します。各ケースとも、実際のイベント（故障・死亡など）の発生時刻と、打ち切り時刻（観測終了時刻）を独立に生成しています。
イベントのフラグは、「1」がイベント (故障や死亡)発生，「0」が打ち切りを表しています。
共通設定
各モデルで使用するサンプル数は n = 1000 件、打ち切り時刻は 0～20 の一様分布から生成しています。
# 共通の設定 n <- 1000 # サンプル数（各モデルごとのデータ数） c_max <- 20 # 打ち切り時刻の上限（打ち切り時刻は一様分布 U(0, c_max) で生成）
ケース1: 指数分布モデル
モデルの考え方
指数分布は、一定のハザード率 lambda を持つ場合の故障（または死亡）モデルです。
生存関数は
となり、平均生存時間は 1/lambda です。
Rスクリプトについて
rexp(n, rate = lambda)でイベント発生時刻（実際の故障時刻）を生成します。- 打ち切り時刻は
runif(n, 0, c_max)により一様分布から生成。 - 観測時刻は
pmin()でイベント時刻と打ち切り時刻の小さい方を選び、イベントが観測されたかどうかはas.numeric(event_time_exp <= censor_time_exp)で判定しています。
# ### ケース1: 指数分布モデル ### # 共通の設定 n <- 1000 # サンプル数（各モデルごとのデータ数） c_max <- 20 # 打ち切り時刻の上限（打ち切り時刻は一様分布 U(0, c_max) で生成） # 指数分布のパラメタ lambda <- 0.1 # ハザード率（例：平均生存時間 10） # イベント発生時刻（実際の故障時刻）の生成 event_time_exp <- rexp(n, rate = lambda) # 打ち切り時刻を一様分布から生成 censor_time_exp <- runif(n, 0, c_max) # 観測時刻はイベント時刻と打ち切り時刻のうち小さい方 time_obs_exp <- pmin(event_time_exp, censor_time_exp) # イベントが観測されたかのフラグ（1: イベント発生, 0: 打ち切り） status_exp <- as.numeric(event_time_exp <= censor_time_exp) # データフレームにまとめる df_exp <- data.frame(time = time_obs_exp, status = status_exp)
Kaplan-Meier生存曲線の描画
survival パッケージを用いで Kaplan-Meier 曲線を描くサンプルを以下に示します。
# パッケージ "survival" のインストールと読み込みが必要です。 if(!require(survival)) install.packages("survival") library(survival) # Kaplan-Meier 曲線の描画例（例: 指数分布モデル） km_fit <- survfit(Surv(time, status) ~ 1, data = df_exp) plot(km_fit, xlab = "Time", ylab = "Survival Probability", col=4, main = "Kaplan-Meier Curve with Theoretical Curve (Exponential Model)") # 理論的な生存曲線（指数分布: S(t) = exp(-lambda * t)）を重ね描き curve(exp(-lambda * x), add = TRUE, col = 2, lwd = 2, lty = 2) # 凡例の追加 legend("topright", legend = c("Kaplan-Meier", "Theoretical (Exponential)"), col = c(4,2), lty = c(1, 2), lwd = c(1, 2))
ケース2: ワイブル分布モデル
モデルの考え方
ワイブル分布は、最弱リンクモデル（Weakest-Link Model） という物理的な故障メカニズムの考え方から導くことができます。このモデルは、システムが複数の独立した部品から成り、最も弱い部品が故障したときにシステム全体が壊れる という状況を表しています。
例えば、鎖（チェーン） を考えたとき、鎖全体の強度は最も弱い鎖の強度によって決まります。これは、ある負荷がかかったとき、鎖の中で最も脆い部分が壊れた瞬間に全体が切れるためです。
ワイブル分布は、形状パラメータと尺度パラメータで特徴付けられる分布で、 故障データの解析に広く用いられます。 生存関数は
です。
Rスクリプトについて
rweibull(n, shape, scale)を使用してイベント発生時刻を生成します。- その他はケース1と同様の手法で、打ち切り時刻や観測時刻、フラグを作成しています。
# ### ケース2: ワイブル分布モデル ### # 共通の設定 n <- 1000 # サンプル数（各モデルごとのデータ数） c_max <- 20 # 打ち切り時刻の上限（打ち切り時刻は一様分布 U(0, c_max) で生成） # ワイブル分布のパラメタ shape <- 2 # 形状パラメータ scale <- 5 # 尺度パラメータ # イベント発生時刻をワイブル分布から生成 event_time_weibull <- rweibull(n, shape = shape, scale = scale) # 打ち切り時刻を一様分布から生成 censor_time_weibull <- runif(n, 0, c_max) # 観測時刻と打ち切りフラグの作成 time_obs_weibull <- pmin(event_time_weibull, censor_time_weibull) status_weibull <- as.numeric(event_time_weibull <= censor_time_weibull) # データフレームにまとめる df_weibull <- data.frame(time = time_obs_weibull, status = status_weibull)
Kaplan-Meier生存曲線の描画
結果を確かめるためのRスクリプトの例は以下です。
# Kaplan-Meier 曲線の描画 library(survival) km_fit <- survfit(Surv(time, status) ~ 1, data = df_weibull) plot(km_fit, xlab = "Time", ylab = "Survival Probability", col=4, main = "Kaplan-Meier Curve with Theoretical Curve (Weibull Model)") # 理論的な生存曲線（Weibull分布の生存関数 S(t)=exp(-(t/scale)^shape)）を重ね描き curve(exp(- (x/scale)^shape), add = TRUE, col = 2, lwd = 2, lty = 2) # 凡例の追加 legend("topright", legend = c("Kaplan-Meier", "Theoretical (Weibull)"), col = c(4, 2), lty = c(1,2), lwd = c(1,2))
ケース3: 生存関数が与えられた場合
モデルの考え方
ここでは、例として生存関数を
と仮定します。
累積分布関数 (cumulative distribution function: CDF) は
となります。
逆関数法を用いて、
一様乱数 から
によりイベント時刻を生成します。逆関数法とは、累積分布関数 (確率分布関数)の逆関数を使って、一様乱数 を目的の確率密度関数に従うように変換する方法です。
Rスクリプトについて
データを生成するRスクリプトはこんな感じです。
# ### ケース3: 与えられた生存関数 S(t)=1/(1+t) の場合 ### # 共通の設定 n <- 1000 # サンプル数（各モデルごとのデータ数） c_max <- 20 # 打ち切り時刻の上限（打ち切り時刻は一様分布 U(0, c_max) で生成） # ランダムな一様乱数 u を生成 u <- runif(n) # 逆関数法によりイベント時刻を計算: t = u/(1-u) event_time_custom <- u / (1 - u) # 打ち切り時刻を一様分布から生成 censor_time_custom <- runif(n, 0, c_max) # 観測時刻と打ち切りフラグの作成 time_obs_custom <- pmin(event_time_custom, censor_time_custom) status_custom <- as.numeric(event_time_custom <= censor_time_custom) # データフレームにまとめる df_custom <- data.frame(time = time_obs_custom, status = status_custom)
Kaplan-Meier生存曲線の描画
結果を確かめるためのRスクリプトの例は以下です。
# Kaplan-Meier 曲線の描画 library(survival) km_fit_custom <- survfit(Surv(time, status) ~ 1, data = df_custom) plot(km_fit_custom, xlab = "Time", ylab = "Survival Probability", col=4, main = "Kaplan-Meier Curve with Theoretical Curve (Custom Model)") # 理論的な生存曲線 (S(t)=1/(1+t)) の重ね描き curve(1/(1+x), add = TRUE, col = 2, lwd = 2, lty = 2) # 凡例の追加 legend("topright", legend = c("Kaplan-Meier", "Theoretical (1/(1+t))"), col = c(4, 2), lty = c(1, 2), lwd = c(1, 2))
最後に一言
最近、町田康さんが YouTube に出演されている動画をよく見ます。町田康さんは、若いころ「町田町蔵」として音楽活動をされていました。憧れというほどではありませんが、私が若いころに最も好きだったアーティストです。
昔、新宿ロフトを拠点に音楽活動をしていたころ、クラブチッタ川崎でのライブの打ち上げで町田さんと直接お話しする機会がありました。当時の私は20歳くらいで、ミュージシャンを目指しながら、工事現場で働きつつバンド活動を続けていました。大学にはほとんど通わず、留年や休学を繰り返していました。
しかし、やがて自分には音楽の才能も文学の才能もないことを痛感し、「ただの勘違い野郎だった」と思い知らされました。その挫折をきっかけに、物理や数学の本を読み始め、大学にも通うようになりました。すると、自分にも理解できることがあると気づき、学ぶことの面白さを感じるようになりました。そして、なぜか研究の機会をいただき、気づけば20年以上、飽きずに研究を続けています。今になって、研究の世界が自分に合っていたのかもしれないと感じています。
そんな私ですが、最近町田康さんの話を聞いていると、意外な共通点があることに気づきました。それは、物事の背後にある構造を深く考えるという点です。
町田康さんの話には、物事の本質を見抜き、それをわかりやすい言葉で説明する力があります。
私の専門分野は異なりますが、物事の数理的な構造が見える気がしています。もちろん、町田康さんと自分を比べるのはおこがましいし、もしかすると、これは人生で 2 度目の大きな勘違いかもしれません。
しかし、私が人生で学んだことは、「勘違いかもしれない」と恐れずに、進み続けることが大切だということです。
