Triad sou.

ある仮定のもとでのKaplan–Meier推定量の正確な分布について

統計

はじめに

この記事はあまり実用的ではないです。
調べた内容を忘れてしまうともったいないような気がしたので、一応メモしておくことにしました。

生存関数のKaplan–Meier推定量$\hat{S}(t)$は非常によく使われていて、$\sqrt{n}\{\hat{S}(t)-S(t)\} \to^d -S(t)W(t)$が成り立つことはよく知られています。
ただし、$W(t)$は$\mathrm{E}[W(t)]=0, \mathrm{Var}[W(t)]=\int_0^t \{\Pr(U>s)S(s)\}^{-1} \mathrm{d}\Lambda(s)$なるブラウン運動、$U$は打切り時間の確率変数です。
この性質から、Kaplan–Meier推定量のpoint-wiseの信頼区間は正規近似を使って構成することができます。
ただし、そのまま正規近似すると$n$が小さかったり生存関数が0や1に近いところで近似精度が悪く、変換を行った統計量の方が正規近似の精度が良いことが知られています（Borgan & Liestøl, 1990など）。
シミュレーションの結果などから、より保守的な二重対数変換（$\log\{-\log S(t)\}$）と分散安定化変換っぽい逆正弦平方根変換（$\mathrm{arc sin}\sqrt{S(t)}$）を利用した近似信頼区間が推奨されています。
普通はこの変換推定量を利用した近似信頼区間が用いられますが、正確な分布に関する研究もいくつかあり、その一つを調べた結果をメモしています。

Chang (1996)

この論文では、Coxの比例ハザードモデルとは別の比例ハザードモデルを考え、その仮定のもとでKaplan–Meier推定量の正確な分布を導いたものでした。

生存時間モデルの基本の仮定

$T_i, U_i$（for $i=1, \ldots, n$）をそれぞれ生存時間と打切り時間の確率変数とし、$T_i$は分布$F$にしたがい$U_i$は分布$G$にしたがいこれらがi.i.d.であるという仮定と、$T_i$と$U_i$の間の独立性を仮定します。
そして$X_i=\min(T_i, U_i)$, $\delta_i=I(T_i \leq U_i)$が観測可能である（生存時間$T_i$は常に観測できるわけではない）というモデルを考えます。
また、確率$K(t)=\Pr(X \leq t)$を考え、この経験推定量を$K_n(t)=n^{-1}\sum_{i=1}^n I(X_i \leq t)$とします。
Kaplan–Meier推定量は生存関数$\bar{F}(t)=1-F(t)$の推定量であり、
\[
\bar{F}_n(t)=\prod_{i=1}^{nK_n(t)} c_{in}^{\delta_{(i)}}I(X_{(n)} \geq t),
\] と表わすことができます。
$X_{(i)}$は$X_{(1)} < X_{(2)} < \ldots < X_{(n)}$となるように$X_i$を並べ替えたもので、$\delta_{(i)}$は$\delta_i$を同じ順番に並べ替えたもので、$c_{in}=(n-i)/(n-i+1)$とします。

比例ハザードモデル

上の設定のもとで、比例ハザードモデル
\[
\bar{G}=1-G=(\bar{F})^{\beta}, \beta>0
\] を仮定します。
この仮定のもとでは、$\delta_1,\ldots,\delta_n$と$X_1,\ldots,X_n$の間の独立性が成り立ちます（Allen, 1963）。

正確な分布

比例ハザードモデルのもとでは$\gamma=\Pr(T \leq U)=1/(1+\beta)$が成り立ちます。
$T$と$U$の独立性から、$\bar{K}(t)=1-K(t)=\bar{F}(t)\bar{G}(t)$が成り立ちます。
また、$X\leq t$を満たす人数を$A_n(t)=n K_n(t)$とし、$t$以下の死亡の人数を$D_n(t)=\sum_{i=1}^{A_n(t)}\delta_{(i)}$とします。
このとき、
\[
\Pr(A_n(t)= q)=\binom{n}{q} \{K(t)\}^q\{\bar{K}(t)\}^{n-q},
\] であり（独立に確率$K(t)$で起こる事象が$n$人中$q$人で観測される確率より）、
\[
\Pr(D_n(t)= i \mid A_n(t)= q)=\binom{q}{i} \gamma^i (1-\gamma)^{q-i},
\] が成り立ちます（$\delta_1,\ldots,\delta_n$と$X_1,\ldots,X_n$の間の独立性より、独立に確率$\gamma$で起こる事象が$q$人中$i$人で観測される確率より）。
すると、Kaplan–Meier推定量が時点$t$で$x$以下となる確率$H_n^t(x)=\Pr(\bar{F}_n(t) \leq x)$（$0 < x \leq 1$）は
\[
\begin{aligned}
H_n^t(x)&=
\sum_{q=0}^n \Pr(A_n(t)= q) \Pr(\bar{F}_n(t) \leq x \mid A_n(t)= q)\\&=
\sum_{q=0}^n \Pr(A_n(t)= q) \left[\sum_{i=0}^q \Pr(D_n(t)= i \mid A_n(t)= q)\Pr(\bar{F}_n(t) \leq x \mid A_n(t)= q, D_n(t)=i) \right]
\end{aligned}
\] となります（周辺分布の定義より）。
ここで、$\Pr(\bar{F}_n(t) \leq x \mid A_n(t)= q, D_n(t)=i)$は
\[
\Pr\left(\prod_{j=1}^q c_{jn}^{\delta_{(j)}} \leq x \mid A_n(t)= q, D_n(t)=i\right)=
\Pr\left(\sum_{j=1}^q \delta_{(j)} \log c_{jn} \leq \log x \mid A_n(t)= q, D_n(t)=i\right),
\] から計算ができます。
具体的には、$q$と$i$が固定されているため、$\sum_{j=1}^q \delta_{(j)}=i$（$j=1,\ldots,q$）を満たすような$\delta_{(j)}$の並べ替え分布を作り、$\sum_{j=1}^q \delta_{(j)} \log c_{jn} \leq \log x$を満たす割合を計算すれば良いです。
自分は以下のコードで計算しました。

hnt <- function(t, x, n, gamma, kt) {
  ghat <- gamma
  prob <- kt^n
  for (q in 0:(n - 1)) {
    if (q == 0) {
      gx <- 1
    } else {
      gx <- 0
      bs <- e1071::bincombinations(q)
      grp <- rowSums(bs)
      bs[,1] <- bs[,1]*log((n - 1)/n)
      if (q > 1) {
        for (j in 2:q) {
          bs[,1] <- bs[,1] + bs[,j]*log((n - j)/(n - j + 1))
        }
      }
      hx <- aggregate(x = data.frame(prob = as.numeric(bs[,1] <= log(x))),
                      by = data.frame(grp), FUN = mean)
      for (i in 0:q) {
        gx <- gx + choose(q, i)*ghat^i*(1 - ghat)^(q - i)*hx[i + 1, 2]
      }
    }
    prob <- prob + gx*choose(n, q)*kt^q*(1 - kt)^(n - q)
  }
  prob
}

信頼区間

正確な分布に基づく信頼区間についてはこの論文には書かれていません。
$H_n^t(x)$を求めるためには$K(t)$と$\gamma$が必要ですが、実際にはデータから求めた推定値$K_n(t)=n^{-1}\sum_{i=1}^n I(X_i \leq t)$と$\hat{\gamma}=n^{-1}\sum_{i=1}^n \delta_i$しか得られないため、これらをプラグインしたものを使ってシミュレーションをしてみました。
また、時点$t$における両側$100(1-\alpha/2)$%信頼区間は$\alpha/2=H_n^t(x_l) < H_n^t(x) < H_n^t(x_u)=1-\alpha/2$から求められるので、これを満たす$x_l$と$x_u$を数値的に探さなくてはなりません。
以下は、下側信頼区間を計算する目的で$x_l$を二分探索する適当なコードです（計算が重くなりすぎるので、精度は適当です）。

sch <- function(prob, t, n, gamma, kt) {
  if (hnt(t, 0.5, n, gamma, kt) < prob) {
    if (hnt(t, 0.75, n, gamma, kt) < prob) {
      a <- 1
      b <- 0.75
    } else {
      a <- 0.75
      b <- 0.5
    }    
  } else {
    if (hnt(t, 0.25, n, gamma, kt) < prob) {
      a <- 0.5
      b <- 0.25
    } else {
      a <- 0.25
      b <- 0
    }
  }
  while (abs(a - b) >= 0.002) {
    if (hnt(t, round((a + b)/2, 3), n, gamma, kt) < prob) {
      b <- round((a + b)/2, 3)
    } else {
      a <- round((a + b)/2, 3)
    }
  }
  b
}

シミュレーションの結果は載せませんが、比例ハザードモデルが成り立つ条件下では、近似信頼区間よりもかなり保守的ではあるものの正確にType I errorが$\alpha$以下になっているような結果にはなりました。
計算方法も適当であり、大きい$q$についてe1071::bincombinations(q)を呼ぶとメモリが足りなくなるので、$n=25$ぐらいまでが計算の限界でした。
特に何の参考にもなりませんが、こういう方法もあるのだなということが分かって勉強になりました。

参考文献

  • Borgan Ø, Liestøl K. A note on confidence intervals and bands for the survival function based on transformations. Scandinavian Journal of Statistics 1990; 17: 35–41.
  • Chang MN. Exact distribution of the Kaplan–Meier estimator under the proportional hazards model. Statistics & Probability Letters 1996; 28: 153–157.
  • Allen WR. A note on conditional probability of failure when hazards are proportiona. Operations Research 1963; 11: 658–659.