ライブラリの読み込み

library(survival)   # Kaplan-Meier 生存解析など
library(cmprsk)     # 競合リスク (cumulative incidence) 解析

1. データの作成(仮想)

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

2. Kaplan-Meier 法: “再発(cause=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) では、他原因死を「再発のリスクから外れる(再発できない)」と明示的に扱います。 そのため、「再発率」をより実態に近い形で推定できるようになります。

3. 競合リスク解析: 累積発症曲線で「再発率」を推定

# 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 法では、再発以外のイベントが「打ち切り」となるため、興味のあるイベントの確率を過大推定するリスクがあります。 累積発症曲線(競合リスク解析)では、他原因死を正しく除外しながら再発率を推定するので、より実態に近い推定が可能です。 競合リスクが無視できない状況(例:再発前に他原因死が頻繁に起こるがん研究など)では、累積発症曲線を用いた解析を検討するとよいでしょう。