生存時間分析生存関数の標準誤差（Greenwoodの公式）

f:id:kimoppy126:20201006235307p:plain

前回の記事では様々な定義による生存関数を見ていきました。
今回はそんな生存関数について、重要な統計量である標準偏差を算出する方法について学んでいきます。
標準誤差を推定することで信頼区間を形成したり、種々の仮説検定に用いることができます。

前回の記事はこちらから↓ www.kimoton.com

生存関数の信頼区間は、特定の生存期間ごとに信頼区間が算出されることから、時点ごとの信頼区間（pointwise confidential interval）と呼ばれたりします。

前準備

前の記事でご紹介した通り、 $\hat{p_j} = (n_j-d_j)/n_j$ としたとき、生存関数のカプラン・マイヤー推定値は、下記の式で表せます。

$$ \hat{S}(t) = \prod_{j=1}^k{\hat{p_j}} $$

これを対数変換することにより、

$$ \log{\hat{S}(t)} = \sum^k_{j=1}{\log{\hat{p_j}}} $$ となります。よって、$\log{\hat{S(t)}}$の分散は、 $$ var\{\log{\hat{S}(t)}\} = \sum^k_{j=1}{var\{\log{\hat{p_j}}\}} $$ となります。また、$p_j$ を真の生存確率とすると、$n_j$, $p_j$ は二項分布に従うため、 $n_j-d_j$ の分散 $var(n_j-d_j)$ は、$n_j$, $p_j$ を用いて下記のように表せます。 $$ var(n_j-d_j) = n_jp_j(1-p_j) $$ $\hat{p_j} = (n_j-d_j)/n_j$ より、 $\hat{p_j}$ の分散 $var\{\hat{p_j}\}$ は、 $$ \begin{align} var\{\hat{p_j}\} & = var(n_j-d_j)/n_j^2\\ & = p_j(1-p_j) / n_j \end{align} $$

となります。この値は、 $\hat{p_j}(1-\hat{p_j}) / n_j$ によって推定できます。

デルタ法

デルタ法では、連続で滑らかな確率関数 $f(X)$ において、確率変数 $X$ の $\mu$ 周りの1次のTaylor展開を考えます。

$$ f(X) \simeq f(\mu) + (X - \mu) f'(\mu)\ $$

$$ \begin{align} var\{f(X)\} & \simeq var\{X - \mu\} \times {f'(\mu)}^2\\ & \simeq var\{X\} \times {f'(\mu)}^2 \end{align} $$

デルタ法による関数の分散の推定量は、

$$ \widehat{var}{f(X)} \simeq \widehat{var}{X} \times {f'(\hat{\mu})}^2 $$

例として、関数 $\ln{(X)}$ を考えると、

$$ \ln{(X)} \simeq \ln{(\mu)} + (X-\mu)\frac{1}{\mu} $$

より、デルタ法による分散の推定量は、

$$ \widehat{var}{\ln{(X)}} \simeq \widehat{var}{X} \times \frac{1}{\hat{\mu}^2} $$

となります。

デルタ法を用いた標準誤差の導出（Greenwood's formula）

準備が整いました。それでは標準偏差を求めていきましょう。

デルタ法より、

$$ \begin{align} \widehat{var}\{\log{(\hat{p_j})}\} & \simeq \frac{1}{\hat{p_j}^2} \times \widehat{var}\{\hat{p_j}\} \\ & = \frac{1}{\hat{p_j}^2} \times \hat{p_j}(1-\hat{p_j}) / n_j \\ & = \frac{d_j}{n_j(n_j-d_j)} \end{align} $$

ここから、上述した $\log{\hat{S}(t)}$ の分散推定量は、

$$ \begin{align} \widehat{var}\{\log{\hat{S}(t)}\} & = \sum^k_{j=1}{\widehat{var}\{\log{(\hat{p_j})}\}} \\ & = \sum^k_{j=1}{\frac{d_j}{n_j(n_j-d_j)}} \end{align} $$

となります。

再度デルタ法により、

$$ \widehat{var}{\log{\hat{S}(t)}} \simeq = \frac{1}{[\hat{S(t)}]^2} \times \widehat{var}{\hat{S}(t)} $$

ゆえに、

$$ \begin{align} \widehat{var}\{\hat{S}(t)\} & = {[\hat{S(t)}]^2} \times \widehat{var}\{\log{\hat{S}(t)}\} \\ & = {[\hat{S(t)}]^2} \times \sum^k_{j=1}{\frac{d_j}{n_j(n_j-d_j)}} \end{align} $$

標準誤差の推定値は下記となります（Greenwood's formula）。

$$ \widehat{se}\{\hat{S}(t)\} \simeq {[\hat{S(t)}]} \times \sqrt{\sum^k_{j=1}{\frac{d_j}{n_j(n_j-d_j)}}} $$

補対数対数変換（二重対数変換）

上述した信頼区間は、対称な形状をしていることから、生存関数の推定値が0または1に近いとき、信頼区間が区間 $(0, 1)$ の外側に外れてしまうという問題が生じます。

このような問題に対する解決法は以下の2通りになります。

信頼限界が1以上の場合は1に、0以下の場合には0に置き換える。
$\hat{S}(t)$を変換し、変換後の値について信頼区間を計算、得られた信頼区間を逆変換することで $S(t)$ の信頼区間を得る

例として、補対数対数変換 $\log{{-\log{S(t)}}}$ による変換を用いた標準偏差の導出を行います。

デルタ法より、

$$ var\{\log(-X)\} \simeq var\{X\} \times \frac{1}{X^2} $$

ここで、 $X=\log{\hat{S}(t)}$ とおくと、

$$ \begin{align} var\{\log(-\log{\hat{S}(t)})\} & \simeq var\{\log{\hat{S}(t)}\} \times \frac{1}{\log{\hat{S}(t)}^2} \\ & =\sum^k_{j=1}{\frac{d_j}{n_j(n_j-d_j)}} \times \frac{1}{\log{\hat{S}(t)}^2} \end{align} $$

この値の平方根が$\log{{-\log{S(t)}}}$ の標準偏差となります。

この上限、下限を $c_1, c_2$ とすると、$exp(-exp(\hat{c_1})), exp(-exp(\hat{c_2}))$ により逆変換してあげることで、 $(-\infty, \infty)$ の間になるように $\hat{S}(t)$ を変換することができます。

計算は省略しますが、$z_{\alpha/2}$ を標準正規分布の上側 $\alpha/2$点とすると、$100(1-\alpha)\%$信頼区間は、

$$ \hat{S}(t)^{\exp{\{\pm z_{\alpha/2}se\{ \log(-\log{\hat{S}(t)})\}\}}} $$

となります（Greenwood’s Exponential formula）。

lifelinesを含む多くの統計ソフトウェアでは、この方法で信頼区間が算出されます。

lifelinesで見る信頼区間

前回同様lifelinesパッケージで提供されているテストデータ（世界中の政治指導者の生涯データ）を用いて生存時間関数の推定について見ていきます。

from lifelines.datasets import load_dd

data = load_dd()
data.head()
#       ctryname  cowcode2  politycode un_region_name  ...         regime start_year duration observed
# 0  Afghanistan       700       700.0  Southern Asia  ...       Monarchy       1946        7        1
# 1  Afghanistan       700       700.0  Southern Asia  ...  Civilian Dict       1953       10        1
# 2  Afghanistan       700       700.0  Southern Asia  ...       Monarchy       1963       10        1
# 3  Afghanistan       700       700.0  Southern Asia  ...  Civilian Dict       1973        5        0
# 4  Afghanistan       700       700.0  Southern Asia  ...  Civilian Dict       1978        1        0
# 
# [5 rows x 12 columns]

T = data["duration"] # 在職期間
E = data["observed"] # 退職が観測されたか否か

カプラン・マイヤー推定量をプロットすると、信頼区間は自動で描画されます。
有意水準 $\alpha$ の値は、KaplanMeierFitter の引数で指定することができます。

from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter(alpha=0.05)
kmf.fit(durations=T, event_observed=E)
kmf.plot()

上限・下限の値は confidence_interval_ 属性に格納されています。

kmf.confidence_interval_.head()
#      KM_estimate_lower_0.95  KM_estimate_upper_0.95
# 0.0                1.000000                1.000000
# 1.0                0.700522                0.741841
# 2.0                0.578805                0.624308
# 3.0                0.487205                0.534126
# 4.0                0.395233                0.442242