この記事のKaplan-Meier推定量が何かわからない方はぜひこちらをご覧ください。
Pythonで学ぶ生存時間解析2 -Kaplan-Meier推定量
https://qiita.com/Goriwaku/items/767cdc2640e29fddf7cc
今回の記事では番外編ということで、Pythonは特に使用せずに、生存関数の信頼区間の構築や分散推定量について見ていきます。高度な内容になっている可能性があるため、わからなければコメントで聞いていただいてもかまいませんし、生存関数についてまだ自信がない方は上の記事をご覧になってから読まれることをお勧めします。
デルタ法による分散導出
今回は、Kaplan-Meier推定量の分散を導出することで信頼区間、および検定のための示唆を与えます。そのために、まずはデルタ法と呼ばれる分散の導出法について見ていきます。
一般に、単純な線形結合でない関数の分散の推定量を求めるとき、通常の手法では求めることが難しいです。そこで登場するのがこのデルタ法です。確率変数$X$に対して関数$f(X)$を考え、その分散についてデルタ法での導出を見ていきます。$f(X)$を$X$の期待値$\mu$の周りで1次の項までテイラー展開し、2次以降の項を無視すると以下のようになります。
f(X) = f(\mu) + (X-\mu)f'(\mu)
ここで、この式の両辺について分散を考えると、$f(\mu)$は定数、$Var(X-\mu)=\sigma^2$の下で、
Var(f(X)) = \sigma^2\{f'(\mu)\}^2
以上の議論より、関数の分散の推定量は
\hat{Var}(f(X))=\hat{\sigma}^2\{f'(\hat{\mu})\}^2
と表すことができます。これはデルタ法による関数の分散の導出です。
Kaplan-Meier推定量の分散
ここでは、先ほどのデルタ法を用いて、Kaplan-Meier推定量の分散の推定量を求めていきます。生存関数$\hat{S}(t)=\prod_{t_i\leq t}\frac{n_i - d_i}{n_i}$に対して、i.i.d.であることとそれぞれの時点におけるリスク集合の観測値が確率一定のベルヌーイ標本であると仮定します。この時、$\hat{p_i}=(n_i-d_i)/n_i$がこの確率の推定値であり、その分散推定量は$\hat{p_i}(1-\hat{p_i})/n_i$となります。したがって、デルタ法を用いて、生存関数の対数、および$\hat{p_i}$の対数の分散は、
\ln(\hat{S}(t))=\sum_{t_i\leq t}\ln\hat{p_i} \\
\hat{Var}[\ln(\hat{p_i})] = \frac{1}{\hat{p_i}^2}\frac{\hat{p_i}(1-\hat{p_i})}{n_i}=\frac{d_i}{n_i(n_i-d_i)}
ここで、各観測値が独立であるため、対数生存関数の分散推定量は、
\hat{Var}[\ln(\hat{S}(t))] = \sum_{t_i\leq t}\frac{d_i}{n_i(n_i-d_i)}
これに対して、$f(X)=\exp(X)$と考えてもう一度デルタ法を用いることで、Kaplan-Meier推定量の分散推定量を求めることができます。手法としては同じなので詳細の計算は省略しますが、これを解くことによって以下の生存関数の分散に関するGreenwoodの公式を得ます。
Var(\hat{S}(t))=(\hat{S}(t))^2\sum_{t_i \leq t}\frac{d_i}{n_i(n_i-d_i)}
以上により、生存関数の分散推定量を求めることができました。ここで詳しく言及はしませんが、計数過程の理論を使用すると、Kaplan-Meier推定量とその関数が漸近的に正規分布に従うことを示すことができます。しかし、生存時間解析においてはこの漸近正規性の前提条件となる十分なサンプルサイズが保証されないことも少なくありません。そこで、次の節で二重対数生存関数に基づく信頼区間の推定量について見ていきます。
二重対数生存関数に基づく信頼区間
二重対数生存関数$\ln[-\ln[\hat{S}(t)]]$は、値域が$-\infty$から$\infty$まであることから正規分布での近似を考えやすいというメリットがあります。前述したデルタ法を利用し、二重対数生存関数の分散推定量を求めると以下のようになります。
\hat{Var}\{\ln[-\ln(\hat{S}(t)]\}=\frac{1}{\ln(\hat{S}(t))}\sum_{t_i \leq t}\frac{d_i}{n_i(n_i-d_i)}
そして、正規近似を考えたときの二重対数生存関数の$100(1-\alpha)\%$信頼区間は、
\ln[-\ln(\hat{S}(t))] \pm z_{1-\frac{\alpha}{2}}\hat{SE}[\ln[-\ln(\hat{S}(t))]]
これに対してcを信頼区間の上限または下限とするとき、$\exp(-\exp(c))$を考えると、生存関数の信頼区間を構築することができます。この信頼区間は、50%程の右側打ち切りを含む25個ほどの小標本のケースでも十分に近似することができるということが研究によって示されているようです(Borgan and Leitol(1990), 私は読んだわけではありませんが興味がある方はぜひ読んでみてください。)。
ちなみに、ここでは触れませんが、二重対数変換以外にも対数変換、ロジット変換、逆制限平方変換などで信頼区間を構築している統計ソフトウェアもあるようです。もし皆さんが生存関数の信頼区間を統計ソフトに任せて構築するならば、ぜひどの変換による信頼区間なのか気にしてみてください。
Hall-Weller信頼幅
前節では、ある時点における生存関数の信頼区間について議論してきましたが、この節では生存関数全体に対する同時信頼幅について考えます。この信頼幅はKaplan-Meier推定量が定義可能な時間、すなわち最大非打ち切り時間$t_m$以下の値に限って使用できるということに注意してください。区間$[0, t_m]$における$100(1-\alpha)\%$信頼幅の端点は、
\ln[-\ln(\hat{S}(t))] \pm H_{\hat{a}, \alpha}\frac{1+n\hat{\sigma}^2}{\sqrt{n}|\ln\hat{S}(t)|}
となります。ただし、$H_{\hat{a}, \alpha}$は$\hat{a}=n\hat{\sigma}^2t_m/(1+n\hat{\sigma}^2t_m)$、$\hat{\sigma}^2$はKaplan-Meier推定量の分散推定量としたときのHall-Wellner分布表から得られる臨界値です。これも、先ほどの例と同様に$\exp(-\exp(c))$を考えると生存関数の同時信頼幅に変換することが可能です。Hall-Wellner信頼幅は各時点における信頼区間を含んでいるため、それを包括する形になっています。これらの概念を利用すると、生存関数の四分位点の信頼区間を考えることも可能になります。生存関数は階段関数のため、四分位点は必ずイベントが起こった時点になることに注意してください。また、一般に打ち切りを伴う生存時間の分布は右に裾が重いことが多く、単に記述統計として値を知りたいだけならば平均値よりも中央値のほうが適切である、ということも覚えておくとよいでしょう。
今回は、生存関数の推定量の導出を行いました。内容が結構重くなってしまいましたが、デルタ法は生存時間解析に留まらずほかの統計の分野でも使える知識だと思うので、これを機に考え方を知っておくのが良いと思います。
参考文献
生存時間解析入門 Hosmer DW, Lemeshow S, May S