library(survival) # Kaplan-Meier 生存解析など
library(cmprsk) # 競合リスク (cumulative incidence) 解析
set.seed(123) # 再現性確保のための乱数シード
n <- 100 # 全体サンプル数
group <- sample(0:1, n, replace = TRUE) # group=0 と group=1 をランダムに割り振り
# 競合リスクを想定した「再発までの時間」と「他原因死までの時間」をランダム生成
# 例として、group=1 の方がやや再発・死亡リスクが高いように設定
# (rateパラメータが大きいほど、より早くイベントが起きやすい)
time_recur <- rexp(n, rate = 0.15 + 0.05 * group) # 再発の時間
time_other <- rexp(n, rate = 0.05 + 0.03 * group) # 他原因死の時間
# フォローアップを12ヶ月(1年)と仮定
censor_time <- 12
# 実際の観察時間 = min(再発時間, 他原因死時間, フォローアップ打ち切り)
obs_time <- pmin(time_recur, time_other, censor_time)
# イベントの種類を付与
# - 再発が先に起きた場合 : cause = 1
# - 他原因死が先に起きた場合 : cause = 2
# - いずれも起こらず12ヶ月で打ち切り : cause = 0 (検閲)
cause <- ifelse(obs_time == time_recur & obs_time < censor_time, 1,
ifelse(obs_time == time_other & obs_time < censor_time, 2, 0))
# データフレームにまとめる
df <- data.frame(
id = 1:n,
group = factor(group, labels = c("Group0","Group1")),
obs_time = obs_time,
cause = cause
)
# 他原因死(cause=2) や イベントなし(cause=0) は打ち切り扱いとする
df$event_km <- ifelse(df$cause == 1, 1, 0)
df # 確認用
#> id group obs_time cause event_km
#> 1 1 Group0 0.65779577 2 0
#> 2 2 Group0 3.06076085 2 0
#> 3 3 Group0 6.85246410 1 1
#> 4 4 Group1 1.42334230 1 1
#> 5 5 Group0 10.42034592 1 1
#> 6 6 Group1 0.07032614 2 0
#> 7 7 Group1 0.05459213 2 0
#> 8 8 Group1 0.49284656 1 1
#> 9 9 Group0 1.86923441 1 1
#> 10 10 Group0 1.97191306 1 1
#> 11 11 Group1 0.64449270 2 0
#> 12 12 Group1 4.62013612 1 1
#> 13 13 Group1 8.21202263 1 1
#> 14 14 Group0 10.79941845 1 1
#> 15 15 Group1 12.00000000 0 0
#> 16 16 Group0 10.14366416 1 1
#> 17 17 Group1 1.90007102 1 1
#> 18 18 Group0 1.59008645 1 1
#> 19 19 Group0 3.10991616 1 1
#> 20 20 Group0 0.28180968 1 1
#> 21 21 Group0 2.13117933 1 1
#> 22 22 Group1 0.85813206 2 0
#> 23 23 Group0 3.81708736 1 1
#> 24 24 Group0 1.44037002 1 1
#> 25 25 Group0 11.16263531 2 0
#> 26 26 Group0 9.17562608 2 0
#> 27 27 Group1 3.25329632 2 0
#> 28 28 Group1 7.19726313 1 1
#> 29 29 Group0 11.54102656 1 1
#> 30 30 Group1 6.22391650 1 1
#> 31 31 Group0 9.75533707 1 1
#> 32 32 Group1 7.57030682 2 0
#> 33 33 Group0 0.03066084 1 1
#> 34 34 Group1 5.54382734 1 1
#> 35 35 Group1 1.27407664 2 0
#> 36 36 Group0 7.94668671 1 1
#> 37 37 Group0 7.43285802 1 1
#> 38 38 Group0 0.44917261 1 1
#> 39 39 Group0 0.48346744 2 0
#> 40 40 Group1 5.94493508 2 0
#> 41 41 Group0 1.73297405 1 1
#> 42 42 Group1 9.28461144 1 1
#> 43 43 Group1 2.31609810 1 1
#> 44 44 Group0 1.57357164 1 1
#> 45 45 Group0 7.88066281 1 1
#> 46 46 Group0 0.39780914 1 1
#> 47 47 Group0 2.68825628 1 1
#> 48 48 Group1 1.41289674 2 0
#> 49 49 Group0 2.77720430 1 1
#> 50 50 Group0 5.02145603 1 1
#> 51 51 Group1 0.94343194 1 1
#> 52 52 Group0 5.84569247 1 1
#> 53 53 Group0 1.26691880 1 1
#> 54 54 Group0 6.52718804 1 1
#> 55 55 Group0 2.15583885 1 1
#> 56 56 Group1 6.60238777 1 1
#> 57 57 Group1 1.59230324 1 1
#> 58 58 Group0 10.70044864 1 1
#> 59 59 Group1 0.72867058 1 1
#> 60 60 Group0 2.89564104 2 0
#> 61 61 Group0 0.20039636 1 1
#> 62 62 Group1 3.22928082 2 0
#> 63 63 Group1 0.99896122 1 1
#> 64 64 Group0 1.50708732 2 0
#> 65 65 Group0 2.56875413 2 0
#> 66 66 Group1 4.21521789 1 1
#> 67 67 Group0 2.32501851 1 1
#> 68 68 Group0 4.10762035 2 0
#> 69 69 Group0 2.67804117 1 1
#> 70 70 Group0 7.33801156 1 1
#> 71 71 Group1 6.64461586 1 1
#> 72 72 Group0 4.27935099 1 1
#> 73 73 Group0 1.30555810 1 1
#> 74 74 Group0 3.04525904 1 1
#> 75 75 Group0 2.48500109 1 1
#> 76 76 Group1 4.51968878 2 0
#> 77 77 Group1 6.37013860 1 1
#> 78 78 Group0 7.20990100 1 1
#> 79 79 Group1 1.48479580 1 1
#> 80 80 Group1 0.42080367 1 1
#> 81 81 Group1 2.78753937 2 0
#> 82 82 Group1 9.84126127 1 1
#> 83 83 Group0 3.85054179 2 0
#> 84 84 Group1 8.04696072 1 1
#> 85 85 Group1 2.62447110 2 0
#> 86 86 Group1 0.46826157 1 1
#> 87 87 Group0 2.16423523 1 1
#> 88 88 Group0 7.81377446 2 0
#> 89 89 Group1 1.33055360 1 1
#> 90 90 Group0 1.49967285 2 0
#> 91 91 Group1 2.24554346 1 1
#> 92 92 Group1 6.44453240 1 1
#> 93 93 Group0 1.18859375 1 1
#> 94 94 Group1 2.09680840 2 0
#> 95 95 Group1 0.31070399 1 1
#> 96 96 Group0 3.80388895 1 1
#> 97 97 Group0 9.82292331 2 0
#> 98 98 Group1 6.55371952 1 1
#> 99 99 Group0 8.63248644 1 1
#> 100 100 Group0 3.11540473 1 1
# Survオブジェクトを作成 (time = obs_time, event=event_km)
km_fit <- survfit(Surv(obs_time, event_km) ~ group, data = df)
# プロット(Kaplan-Meier 生存曲線 → 再発しない確率の推移)
# 通常の慣習では「生存確率(再発しない確率)」で描かれるので、
# 『1 - 曲線』を見れば「累積再発率」のイメージになります。
plot(
km_fit,
conf.int = FALSE,
xlim = c(0, 12),
xlab = "Months",
ylab = "Probability of being recurrence-free",
col = c("red","blue"),
main = "Kaplan-Meier: no-recurrence"
)
legend("bottomleft", legend = c("Group0","Group1"), col=c("red","blue"), lty=1)
解釈 上記プロットは “再発が起こっていない確率” の推移を示します。 再発率(= 1 - 生存確率)を見たい場合は、1 - 生存確率 の形で見ることもできますが、 ここでは典型的な “生存曲線” として描画しました。 event == 2(他原因死)は 再発せずに打ち切られた として扱われています。 そのため、実際より再発率がやや高く推定される可能性があります。 2. 競合リスク解析による累積発症曲線 一方、累積発症曲線 (cumulative incidence) では、他原因死を「再発のリスクから外れる(再発できない)」と明示的に扱います。 そのため、「再発率」をより実態に近い形で推定できるようになります。
# cmprsk パッケージの cuminc() を用いる
# ここでは単純に '再発 (1)' と '他原因死(2)' を2つの競合リスクとして解析
# グループ別に解析したいため、グループをサブセットして cuminc() を実行
# fstatus: 1→再発, 2→他原因死, cencode=0→検閲
fit_ci <- cuminc(ftime = df$obs_time, fstatus = df$cause,
group = df$group, cencode = 0)
# cuminc() の結果オブジェクトは、各グループ・各"cause"の累積発症率を含む
# ここでは cause=1(再発) のみを図示したいので、plot() の引数で指定する
# plot.cuminc のデフォルトでは cause=1,2 両方の曲線を描くので、
# "再発(1)" だけを表示するには 'curves = "1"' を指定する
plot(
fit_ci,
xlim = c(0,12),
ylim = c(0,1),
lty = c(1,1,2,2),
col = c("red","blue", "red", "blue"), # Group0, Group1
main = "CI curve: recurrence (1) and death (2)",
xlab = "Months",
ylab = "Cumulative incidence of recurrence",
curves = "1" # 「再発(cause=1)」だけ描画
)
解釈 cuminc() の出力では、それぞれのイベント (ここでは 1=再発, 2=他原因死) の累積発症曲線を同時に計算しています。 グラフでは、青色の線が「再発」の累積発症曲線、赤色の線が「他原因死」の累積発症曲線。 Kaplan–Meier 法で「再発率」を見たときよりも、累積発症曲線の再発率はやや低めになる傾向が多いです。 なぜなら「他原因死」が起こった患者は、再発を経験するリスクから除外されるためです。 3. 二つの解析結果の違いを確認 Kaplan–Meier 法 他原因死は「再発なしの打ち切り」として扱い、再発率を過大推定する場合がある。 全生存 (OS) など「単一のイベント(死)」が注目される場面では有用。 累積発症曲線 (Competing Risks) 他原因死をしっかり競合イベントとして扱い、再発できない状況をモデルに反映。 再発率や特定死因の発生率など、競合リスクが多い状況で適切にイベント率を推定できる。
このR Markdownファイルでは、同じ仮想データを使って
Kaplan–Meier 法による再発曲線 競合リスク解析 (累積発症曲線) を比較しました。
K-M 法では、再発以外のイベントが「打ち切り」となるため、興味のあるイベントの確率を過大推定するリスクがあります。 累積発症曲線(競合リスク解析)では、他原因死を正しく除外しながら再発率を推定するので、より実態に近い推定が可能です。 競合リスクが無視できない状況(例:再発前に他原因死が頻繁に起こるがん研究など)では、累積発症曲線を用いた解析を検討するとよいでしょう。