当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.$19$の「回帰分析その他」に関して演習問題を中心に解説を行います。
本章のまとめ
生存時間解析・ハザード関数
生存時間解析(survival time analysis)やハザード関数(hazard function)を理解するにあたっては、累積分布関数と確率密度関数を元に定義を確認することが重要である。
$0$以上の離散確率変数$T \geq 0$に対して累積分布関数を$F(t) = P(T \leq t)$、確率関数を$f(t) = P(T=t)$のように定義する。このとき、$F(t)$と$f(t)$には下記のような関係式が成立する。
$$
\large
\begin{align}
F(t) = P(T \leq t) = \sum_{x=0}^{t} P(T=x) = \sum_{x=0}^{t} f(x)
\end{align}
$$
同様に連続確率変数$T \geq 0$に関して累積分布関数$F(t)$と確率密度関数$f(t)$を考えるとそれぞれ下記のように表すことができる。
$$
\large
\begin{align}
f(t) &= \lim_{\epsilon \to 0} \frac{P(t \leq T \leq t+\epsilon)}{\epsilon} \\
F(t) &= P(T \leq t) = \int_{0}^{t} f(x) dx
\end{align}
$$
上記に対し、生存関数$S(t)$は累積分布関数で計測する事象が観測されなかった確率であるので、累積分布関数で表した事象に関する余事象を考えればよく、これは$S(t) = 1 – F(t)$で定義することができる。
また、$S(t)$と$f(t)$を用いてハザード関数$h(t)$は下記のように定義される。
$$
\large
\begin{align}
h(t) = \frac{f(t)}{S(t)} = \frac{f(t)}{1-F(t)}, \quad t \geq 0
\end{align}
$$
指数分布のハザード関数
ハザード関数の具体例を考える際に式がシンプルであることから指数分布がよく用いられる。パラメータ$\lambda$の指数分布$Ex(\lambda)$の累積分布関数$F(t)$と確率密度関数$f(t)$はそれぞれ下記のように表される。
$$
\large
\begin{align}
f(t) &= \lambda e^{-\lambda t}, \quad t \geq 0 \\
F(t) &= \int_{0}^{t} f(x) dx = \lambda \int_{0}^{t} e^{-\lambda x} dx \\
&= \left[ -e^{-\lambda x} \right]_{0}^{t} = 1 – e^{-\lambda t}, \quad t \geq 0
\end{align}
$$
上記に対して生存関数$S(t)$とハザード関数$h(t)$は下記のように導出できる。
$$
\large
\begin{align}
S(t) &= 1 – F(t) = 1 – (1 – e^{-\lambda t}) \\
&= e^{-\lambda t}, \quad t \geq 0 \\
h(t) &= \frac{f(t)}{S(t)} \\
&= \frac{\lambda e^{-\lambda t}}{e^{-\lambda t}} \\
&= \lambda
\end{align}
$$
上記より指数分布におけるハザード関数は$h(t)=\lambda$であり、$\lambda$は変数$t$に関して考えたとき変化しない定数であることに着目しておくと良い。
ハザード関数が条件付き確率と同様の式で表されると考えることで、ハザード関数は「これまで故障しなかったものがそのタイミングで故障する確率」のように解釈することもできる。これを「瞬間故障率」のように表すことがある。
・参考
赤本 章末課題$6.6$： 指数分布の瞬間故障率
ハザード関数と微分方程式
以下はワークブック$2$章の「確率分布と母関数」や「$1$級テキストの$2$章」の内容も含むが、$19$章の内容との関連が大きいので、当記事で取りまとめを行う。
$$
\large
\begin{align}
h(t) = \frac{f(t)}{S(t)} = \frac{f(t)}{1-F(t)}, \quad t \geq 0
\end{align}
$$
ハザード関数$h(x)$は上記のように定義されるが、$\log{S(x)}$の微分を考えることで下記のような関係式で表すこともできる。
$$
\large
\begin{align}
\frac{d}{dx} (\log{S(x)}) &= \frac{S(x)’}{S(x)} \\
&= \frac{(1-F(x))’}{S(x)} \\
&= \frac{-f(x)}{S(x)} = -h(x)
\end{align}
$$
上記を元に下記の微分方程式を$F(0)=0$に基づいて解くことでハザード関数$h(x)$から累積分布関数$F(x)$や確率密度関数$f(x)$の導出を行うことができる。
$$
\large
\begin{align}
\frac{d}{dx} (\log{(1-F(x))}) = -h(x)
\end{align}
$$
「指数関数」や「ワイブル分布」に関する具体例は下記で取り扱った。
・微分方程式を用いた指数分布・ワイブル分布の確率密度関数の導出
演習問題解説
問19.1
$[1]$
左打ち切りに関する$L=0$であることから、次のように尤度関数$L(\mathbf{\beta},\sigma)$を考えることができる。
$$
\large
\begin{align}
L(\mathbf{\beta},\sigma) &= \prod_{i \in \{i|y_i>L\}} \frac{1}{\sigma} \varphi \left( \frac{y_i-\mathbf{x}_i^{T} \mathbf{\beta}}{\sigma} \right) \prod_{i \in \{i|y_i \leq L\}} \Phi \left( \frac{L-\mathbf{x}_i^{T} \mathbf{\beta}}{\sigma} \right) \\
&= \prod_{i \in \{i|y_i>0\}} \frac{1}{\sigma} \varphi \left( \frac{y_i-\beta_0-\beta_1x_{i1}-\beta_2x_{i2}}{\sigma} \right) \\
& \times \prod_{i \in \{i|y_i \leq 0\}} \Phi \left( \frac{L-\beta_0-\beta_1x_{i1}-\beta_2x_{i2}}{\sigma} \right)
\end{align}
$$
$[2]$
AICはそれぞれ下記のように計算できる。
import numpy as np
logL = np.array([-262.3469, -261.9037, -261.9794, -261.8608])
num_param = np.array([4., 5., 5., 6.])
AIC = -2*logL + 2*num_param
print("AIC: {}".format(AIC))
・実行結果
> print("AIC: {}".format(AIC))
AIC: [ 532.6938 533.8074 533.9588 535.7216]
AICは小さい方が良いので、上記より、$1$つ目が良いことが確認できる。
問19.2
確率密度関数を$f(t)$、累積分布関数を$F(t)$とするとき、ハザード関数$h(t)$は下記のように定義される。
$$
\large
\begin{align}
h(t) = \frac{f(t)}{1-F(t)}
\end{align}
$$
ここで$f(t) = \lambda e^{-\lambda t}$であり、$F(t)$は下記のように導出できる。
$$
\large
\begin{align}
F(t) &= \int_{0}^{t} \lambda e^{-\lambda x} dx \\
&= \left[ -e^{-\lambda x} \right]_{0}^{t} \\
&= 1 – e^{-\lambda t}
\end{align}
$$
よって、ハザード関数$f(t)$は下記のように得られる。
$$
\large
\begin{align}
h(t) &= \frac{f(t)}{1-F(t)} \\
&= \frac{\lambda e^{-\lambda t}}{1-(1 – e^{-\lambda t})} \\
&= \frac{\lambda e^{-\lambda t}}{e^{-\lambda t}} = \lambda
\end{align}
$$
問19.3
$f(t|x)=(-\log{S(t|x)})’$より下記が成立する。
$$
\large
\begin{align}
f(t|x) &= (-\log{S(t|x)})’ \\
-\log{S(t|x)} &= \int_{0}^{t} f(u|x) du \\
-\log{S(t|x)} &= \int_{0}^{t} h_0(u) \exp(x^{T} \beta) du \\
&= \exp(x^{T} \beta) \int_{0}^{t} h_0(u) du \\
&= \exp(x^{T} \beta) H_0(t)
\end{align}
$$
上記の両辺の対数を取ることで下記が導出できる。
$$
\large
\begin{align}
-\log{S(t|x)} &= \exp(x^{T} \beta) H_0(t) \\
\log{(-\log{S(t|x)})} &= x^{T} \beta + \log{H_0(t)}
\end{align}
$$
また、上記で導出した式は式$(19.5)$に一致することも確認しておくとよい。
問19.4
$[1]$
式$(19.5)$の$\mathbf{x}$をダミー変数$x_1 \in \{0,1\}$を用いて、$\mathbf{x} = x_1$のように表すことを考える。このとき、式$(19.5)$に基づいて下記のような連立方程式が得られる。
$$
\large
\begin{align}
\log{(-\log(S(t|x_1=1)))} &= \beta + \log{H_0(t)} \\
\log{(-\log(S(t|x_1=0)))} &= \log{H_0(t)}
\end{align}
$$
上記で表した二つの式より$\log{H_0(t)}$を消去すると下記が得られる。
$$
\large
\begin{align}
\log{(-\log(S(t|x_1=1)))} = \beta + \log{(-\log(S(t|x_1=0)))} \quad (1)
\end{align}
$$
上記より比例ハザード性が成立する場合は、$x_1 \in \{0,1\}$で表した$2$群に対して生存関数の$2$重対数$\log{(-\log{S})}$が平行となると考えることができる。
ここで対象のグラフは生存関数のノンパラメトリック推定量に基づくカプランマイヤー推定量の曲線が平行していることから、グラフで表された検証結果に対して比例ハザード性を仮定することが妥当であると考えられる。
$[2]$
$[1]$で導出を行なった$(1)$は下記のように変形を行うことができる。
$$
\large
\begin{align}
\beta = \log{(-\log(S(t|x_1=1)))} – \log{(-\log(S(t|x_1=0)))}
\end{align}
$$
上記より二つのカプランマイヤー推定量の曲線の差は治療効果の大きさ$\beta$に一致すると考えることができる。
参考
・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1
・統計検定1級 統計応用 理工学の生存関数の出題例
https://www.hello-statisticians.com/toukei-kentei-1/stat_app/stat_certifi_1_app_sci_19_1.html
