ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rで生存時間解析】Kaplan-Meier (カプラン・マイヤー)法による生存曲線の推定

 Kaplan-Meier (カプラン・マイヤー)法は、生存曲線 (survival curve)、あるいは、生存関数と呼ばれる、生存確率の推移を推定するための方法です。生存曲線は、時刻ゼロで全員が生きていた状態から出発し、時間の経過とともに、残念ながら、亡くなってしまう方が発生する状況をグラフで表したものです。臨床研究では、予後予測の精度を示したり、治療・介入の効果を示したりする目的で、生存曲線が登場します。

 生存曲線を推定するためには、患者様がいつ亡くなったかという情報 (日付)を知る必要があります。しかし、患者様が、他の病院にうつったり、引っ越したりして、ご家族などとも連絡がつかなくなり、生存の有無が確認できなくなることがあります。そのように、追跡ができなくなったデータを、打ち切りデータ (censored data)と呼びます。Kaplan-Meier (カプラン・マイヤー)法は、打ち切りデータが存在する場合にも有効なので、医学統計で頻繁に使われます。

生存関数

 生存関数 \displaystyle{
S(t)
} は、時間  t までに対象が 生存している確率 を表します:

\displaystyle{
S(t) = P (T > t)

}

ここで、 T生存時間(故障時間、イベント発生時間) を表す確率変数です。

生存関数のグラフ (生存曲線)。グラフは次第に減少し、減少が大きいほど死亡率が高く、良くない状態であることを意味する。

カプラン・マイヤー推定量の導出

 カプラン・マイヤー法では、時系列データにおける生存確率を 段階的に更新しながら推定 します。

記号の定義

  • \displaystyle{
t _ 1 \lt t _ 2 \lt \dots \lt t _ k
} : イベント(故障・死亡)が発生した時点
  • \displaystyle{
d _ i
} : 時点 \displaystyle{
t _ i
} でイベント (死亡など)が発生した人数(個体数)
  • \displaystyle{
n _ i
} : 時点 \displaystyle{
t _ i
} 直前で 生存していた 人数
  • \displaystyle{
S(t)
} : 生存関数(生存確率)

変数の定義と生存割合

生存確率の逐次更新

 ある時点 \displaystyle{
 t _ i
} での生存確率は、直前の生存確率 に、イベントが起こらなかった確率(=その時点で生存していた人のうち生存し続けた割合)を掛けることで求められます:

\displaystyle{
S(t_i) = S(t_{i-1}) \times \left( \frac{n_i - d_i}{n_i} \right)
}

これは、それまでに生存していた割合

\displaystyle{
S(t _ {i-1})
}

に、その時点 \displaystyle{
t _ i
} での 生存割合(何人中、何人生存していたか)

\displaystyle{
\frac{n _ i - d _ i}{n _ i}
}

を掛けているだけです。

この関係を初期条件 \displaystyle{
S(0) = 1
} から繰り返し適用すると、カプラン・マイヤー推定量は以下の形になります:

\displaystyle{
\hat{S}(t) = \prod_{t_i \leq t} \left( \frac{n_i - d_i}{n_i} \right)
}

ここで、積は \displaystyle{
t
} 以下のすべてのイベント発生時点 \displaystyle{
t _ i
} について計算されます。

 推測統計の式の書き方に慣れていないと、「\displaystyle{
{S}(t)
}\displaystyle{
\hat{S}(t)
} を、何でわざわざ区別するの?」と疑問に感じるかもしれません。しかし、統計の考え方では、神様だけが知っている真の答え (正解) \displaystyle{
{S}(t)
} と私たちが限られた観測データを使ってその答えを推測したもの (推定量) \displaystyle{
\hat{S}(t)
} を区別して表現します。

 仮に、以下のようなデータがあるとします。〇が「発生したこと」を表します。

個体 生存時間 (月) イベント発生 打ち切り
A 2
B 3
C 5
D 5
E 6

このとき、各時点での \displaystyle{
n _ i
}(生存者数)と \displaystyle{
d _ i
}(イベント数)は次のようになります:

時間 \displaystyle{t _ i} 生存者 \displaystyle{n _ i} イベント数 \displaystyle{d _ i}
2 5 1
3 4 1
5 3 2

カプラン・マイヤー推定量を計算すると、以下のようになります。

\displaystyle{
S(2) = S(0) \times \frac{5-1}{5} = 1 \times \frac{4}{5} = 0.8
}
\displaystyle{
S(3) = S(2) \times \frac{4-1}{4} = 0.8 \times 0.75 = 0.6

}
\displaystyle{
S(5) = S(3) \times \frac{3-2}{3} = 0.6 \times \frac{1}{3} = 0.2
}

したがって、推定された生存関数は次のようになります。

時間 \displaystyle{t } 生存確率 \displaystyle{ S(t) }
\displaystyle{t \lt 2} 1.0
2 0.8
3 0.6
5 0.2
\displaystyle{t \gt 5} 0.2

カプラン・マイヤー曲線

 カプラン・マイヤーの生存曲線は、時間を横軸、推定された生存確率 \displaystyle{
S(t)
} を縦軸に取る 階段状のグラフ になります。生存イベントが発生するたびに確率が下がり、打ち切りデータ(右打ち切り)が発生した場合はそのまま継続されます。

上の例の生存曲線

Rスクリプトの例

# サンプルデータの作成
# 10人の生存データ(time: 生存時間,status: イベント発生の有無)
time <- c(1, 2, 5, 5, 6, 7, 10, 12, 15, 18)  # 生存時間 (故障時間)
status <- c(1, 1, 1, 1, 0, 1, 0, 1, 0, 1)    # 1: イベント発生, 0: 打ち切り

# イベント発生時点を取得
event_times <- sort(unique(time[status == 1]))

# カプラン・マイヤー推定量の計算
n <- length(time)  # 総サンプル数
surv_prob <- 1     # 初期生存確率は1
S <- c(surv_prob)  # 生存確率のリスト
t_values <- c(0)   # 時間のリスト

for (t in event_times) {
    n_i <- sum(time >= t)  # その時点での生存者数
    d_i <- sum(time == t & status == 1)  # その時点でイベントが起こった人数
    surv_prob <- surv_prob * (n_i - d_i) / n_i  # 生存確率を更新
    S <- c(S, surv_prob)
    t_values <- c(t_values, t)
}

# 生存曲線のプロット
plot(t_values, S, type="s", col=2, lwd=3, ylim=c(0,1), las=1, 
     xlab="Time", ylab="Survival Probability", main="Kaplan-Meier Survival Curve")

信頼区間の推定

 カプラン・マイヤー推定量には、信頼区間を計算する方法もあります。よく使われるのが グリーンウッドの公式(Greenwood's formula) による標準誤差(SE)の推定です:

\displaystyle{
\text{Var} \left(\hat{S}(t) \right) = \hat{S}^2(t) \sum_{t_i \leq t} \frac{d_i}{n_i (n_i - d_i)}
}

これを用いて、\displaystyle{
95\%
} 信頼区間を次のように求めます:

\displaystyle{
\hat{S}(t) \pm 1.96 \times {\rm SE} \left(\hat{S}(t) \right) = \hat{S}(t) \pm 1.96 \times \sqrt{\text{Var} \left(\hat{S}(t) \right) }
}

Rスクリプトの例

# サンプルデータの作成
# 10人の生存データ(time: 生存時間,status: イベント発生の有無)
time <- c(1, 2, 5, 5, 6, 7, 10, 12, 15, 18)  # 生存時間 (故障時間)
status <- c(1, 1, 1, 1, 0, 1, 0, 1, 0, 1)    # 1: イベント発生, 0: 打ち切り

# イベントが発生時点を取得
event_times <- sort(unique(time[status == 1]))

# カプラン・マイヤー推定量の計算とグリーンウッドの公式
n <- length(time)  # 総サンプル数
surv_prob <- 1     # 初期生存確率
S <- c(surv_prob)  # 生存確率のリスト
t_values <- c(0)   # 時間のリスト
varS <- 0          # 分散の初期値
SE <- c(0)         # 標準誤差のリスト
lower_CI <- c(surv_prob) # 95% 信頼区間(下限)
upper_CI <- c(surv_prob) # 95% 信頼区間(上限)

# 分散計算用の累積和
sum_term <- 0

for (t in event_times) {
    n_i <- sum(time >= t)  # その時点での生存者数
    d_i <- sum(time == t & status == 1)  # その時点でのイベント数
    
    # カプラン・マイヤー生存確率の更新
    surv_prob <- surv_prob * (n_i - d_i) / n_i
    S <- c(S, surv_prob)
    t_values <- c(t_values, t)
    
    # グリーンウッドの公式で分散を計算
    sum_term <- sum_term + (d_i / (n_i * (n_i - d_i)))
    varS <- S[length(S)]^2 * sum_term
    SE_current <- sqrt(varS)
    SE <- c(SE, SE_current)
    
    # 信頼区間の計算(95% 信頼区間)
    lower_CI <- c(lower_CI, max(0, S[length(S)] - 1.96 * SE_current))  # 下限
    upper_CI <- c(upper_CI, min(1, S[length(S)] + 1.96 * SE_current))  # 上限
}

# 生存曲線と信頼区間のプロット
plot(t_values, S, type="s", col=2, lwd=3, ylim=c(0,1), las=1, 
     xlab="Time", ylab="Survival Probability", main="Kaplan-Meier Survival Curve")

# 信頼区間を描画(破線)
lines(t_values, lower_CI, type="s", col="pink", lwd=2, lty=2)  # 下限
lines(t_values, upper_CI, type="s", col="pink", lwd=2, lty=2)  # 上限

# 凡例の追加
legend("topright", legend=c("Survival Probability", "95% Confidence Interval"), 
       col=c(2, "pink"), lwd=c(3,2), lty=c(1,2))

まとめ

  • カプラン・マイヤー法は、生存関数 \displaystyle{
S(t)
}段階的な確率の積の形 で推定する方法。

  • 生存確率は、各時点での 生存割合 を掛け合わせることで求められる。

  • グラフは 階段状 になり、打ち切りデータ(censored data)も考慮できる。

  • 信頼区間の推定には グリーンウッドの公式 を利用する。

医学研究では、この方法を使って 生存曲線を描き、予後予測の精度、治療効果の比較 などに活用します。