$\gdef \vec#1{\boldsymbol{#1}} \\ \gdef \rank {\mathrm{rank}} \\ \gdef \det {\mathrm{det}} \\ \gdef \Bern {\mathrm{Bern}} \\ \gdef \Bin {\mathrm{Bin}} \\ \gdef \Mn {\mathrm{Mn}} \\ \gdef \Cov {\mathrm{Cov}} \\ \gdef \Po {\mathrm{Po}} \\ \gdef \HG {\mathrm{HG}} \\ \gdef \Geo {\mathrm{Geo}}\\ \gdef \N {\mathrm{N}} \\ \gdef \LN {\mathrm{LN}} \\ \gdef \U {\mathrm{U}} \\ \gdef \t {\mathrm{t}} \\ \gdef \F {\mathrm{F}} \\ \gdef \Exp {\mathrm{Exp}} \\ \gdef \Ga {\mathrm{Ga}} \\ \gdef \Be {\mathrm{Be}} \\ \gdef \NB {\mathrm{NB}} \\ \gdef \indep {\mathop{\perp\!\!\!\!\perp}} \\ \gdef \tr {\mathrm{tr}}$
ここからはいよいよ生存時間解析ですか。
いくつか生存時間$T$についての関数がでてきたかと思いますが、なかなか整理することができていいません。。
ここでは生存時間$T$についての関数として5種類を扱います。
この5種類のうち1種類が決まれば、残り4種類は自動的に定ります。
5種類の関数の関係性を整理して覚えましょう。
1. 生存時間$T$についての関数
生存時間$T$についての関数として以下の5種類を扱います。
①生存関数$S(t)$
②分布関数$F(t)$
③確率密度関数$f(t)$
④ハザード関数$h(t)$
⑤累積ハザード関数$H(t)$
以下では各関数について、その定義と直感的意味合いを確認していきます。
1-1. 生存関数$S(t)$
（定義）
$$\begin{aligned} S(t) &= {\small (時間tまでイベントが発生していない確率)} \\ &= Pr \{ T \gt t \} \end{aligned}$$
（直感的意味合い）
イベントを『死亡』とすると、
$$\begin{aligned} S(t) &= {\small (時間tまで『死亡』が発生していない確率)} \\ &= {\small (時間tまで生存している確率)} \end{aligned}$$となります。
生存時間解析では生存関数$S(t)$についての検定を行うため、生存関数$S(t)$は生存時間解析の主役とも捉えられます。
1-2. 分布関数$F(t)$
（定義）
$$\begin{aligned} F(t) &= {\small (時間tまでにイベントが発生している確率)} \\ &= Pr \{ T \leqq t \} \end{aligned}$$
（直感的意味合い）
イベントを『死亡』とすると、
$$\begin{aligned} F(t) &= {\small (時間tまで『死亡』が発生している確率)} \\ &= {\small (時間tまでに死亡した確率)} \end{aligned}$$となります。
（補足）
定義から明らかですが、$S(t)$と$F(t)$には、
$$\begin{aligned} F(t) = 1-S(t) \end{aligned}$$の関係式が成立します。
1-3. 確率密度関数$f(t)$
（定義）
$$\begin{aligned} f(t) &= {\small (積分してF(t)となる関数)} \\ &= F'(t) \end{aligned}$$
（直感的意味合い）
$f(t)$についてはややイメージしがたいですが、これは$f(t)$が確率関数ではなく確率密度関数であるゆえです。イベントを『死亡』とすると、$[t,t+dt]$に死亡する確率は、$f(t)dt$、と表されます。
1-4. ハザード関数$h(t)$
（定義）
$$\begin{aligned} h(t) = \lim_{\varDelta t \to 0} \frac{ Pr \{ t \leqq T \leqq T + \varDelta t | T \geqq t \} }{\varDelta t} \end{aligned}$$
（直感的意味合い）
$h(t)$についても$f(t)$同様にややイメージしがたいですが、これは$h(t)$が『ある時点でのイベントの発生速度』であるゆえです。イベントを『死亡』とすると、時点$t$まで生存したという条件の下で$[t,t+dt]$に死亡する確率は、$h(t)dt$、と表されます。
1-5. 累積ハザード関数$H(t)$
（定義）
$$\begin{aligned} H(t) &= {\small (ハザード関数h(s)をs:0 \sim tまで累積したもの)} \\ &= \int_{0}^{t} h(s) ds \end{aligned}$$
（直感的意味合い）
『$H(t) \leftrightarrow h(t)$』の関係は『$F(t) \leftrightarrow f(t)$』の関係と同様です。
2. 生存時間$T$についての関数たちの関係性
（ここでは生存時間$T$が連続変数である場合を扱います）
生存時間$T$についての関数たちの関係性は以下のFig1の様にまとめられます。
重要なことは、いずれか1種類が決まれば残り4種類は自動的に定まる、ということです。
Fig1.
この関係性のうち、$S(t)$と$h(t)$の関係性を除いてはすでに『1. 生存時間$T$についての関数』で示しているので、以下$S(t)$と$h(t)$の関係性を示します。
（$S(t)$と$h(t)$の関係性）
$$\begin{aligned} h(t) &= \lim_{\varDelta t \to 0} \frac{ Pr \{ t \leqq T \leqq T + \varDelta t | T \geqq t \} }{\varDelta t} \\ &{\scriptsize(定義より)} \\[10px] &= \lim_{\varDelta t \to 0} \frac{1}{\varDelta t} \cdot \frac{ Pr \{ (t \leqq T \leqq T + \varDelta t) \cap (T \geqq t) \} } { Pr \{ T \geqq t \} } \\ &{\scriptsize(条件つき確率の変形より)} \\[10px] &= \lim_{\varDelta t \to 0} \frac{1}{\varDelta t} \cdot \frac{ Pr \{ t \leqq T \leqq T + \varDelta t \} } { Pr \{ T \geqq t \} } \\[10px] &= \lim_{\varDelta t \to 0} \frac{1}{\varDelta t} \cdot \frac{ S(t)-S(t + \varDelta t) }{ S(t) } \\[10px] &= -\frac{S'(t)}{S(t)} \\[10px] &= -\frac{d}{dt} \log S(t) \end{aligned}$$
（おまけ1. $S(t)$と$H(t)$の関係性）
$$\begin{aligned} H(t) &= \int_{0}^{t} h(s) ds \\[10px] &= [ -\log S(t)]_{0}^{t} \\[10px] &= -\log S(t) + \log S(0) \\[10px] &= -\log S(t) + \log 1 \\ &{\scriptsize(S(0) = 1であるため)} \\[10px] &= -\log S(t) \end{aligned}$$
（おまけ2. $f(t),h(t),S(t)$の関係式）
$$\begin{aligned} h(t) &= -\frac{S'(t)}{S(t)} \\[10px] &= \frac{F'(t)}{S(t)} \\ &{\scriptsize(F(t) = 1-S(t)よりF'(t)=-S'(t))} \\[10px] &= \frac{f(t)}{S(t)} \\[10px] \Rightarrow f(t) &= h(t) S(t) \end{aligned}$$
おまけ1,2はこんな関係性もあるのかと頭の片隅においておくだけで構いません。
Fig1さえ覚えていれば簡単に導出することができるからです。
例題1.
生存時間$T \sim \Exp(\lambda)$である時、即ち、$T$の確率密度関数$f(t)$が
$$\begin{aligned} f(t) = \lambda e^{-\lambda t} \end{aligned}$$である時、生存関数$S(t)$、分布関数$F(t)$、ハザード関数$h(t)$、累積ハザード関数$H(t)$を求めよ。
解答.
$$\begin{aligned} F(t) &= \int_{0}^{t} f(s) ds \\[10px] &= \int_{0}^{t} \lambda e^{-\lambda s} ds \\[10px] &= [-e^{-\lambda s}]_{0}^{t} \\[10px] &= 1-e^{-\lambda t} \\[30px] S(t) &= 1-F(t) \\[10px] &= e^{-\lambda t} \\[30px] h(t) &= -\frac{d}{dt} \log S(t) \\[10px] &= -\frac{d}{dt} (-\lambda t) \\[10px] &= \lambda \\[30px] H(t) &= \int_{0}^{t} h(s) ds \\[10px] &= \lambda t \end{aligned}$$
生存時間$T$が指数分布に従うならば、ハザード関数$h(t)=Const$となります。
逆にハザード関数$h(t)=Const$ならば、生存時間$T$は指数分布に従います。
これを示すには$h(t)=Const$とした時、生存時間$T$の確率密度関数$f(t)$が指数分布の形になってることを示せばよいです。
なおイベントを「死亡」とした場合、生存時間$T$が指数分布に従う（ハザードが一定）という設定は、健康状態が安定している人のよいモデルとなります。
例題2.
生存時間$T$がワイブル分布に従う時、即ち、$T$の確率密度関数$f(t)$が
$$\begin{aligned} f(t) = \lambda p (\lambda t)^{p-1} e^{-(\lambda t)^p} \end{aligned}$$である時、生存関数$S(t)$、分布関数$F(t)$、ハザード関数$h(t)$、累積ハザード関数$H(t)$を求めよ。
解答.
（いきなり$F(t) = \int_{0}^{t} f(s) ds$に$f(s)$を代入して積分するのは大変なので準備をします）
まず、
$$\begin{aligned} \frac{d}{dt} e^{-(\lambda t)^p} &= e^{-(\lambda t)^p} \cdot \frac{d}{dt} (-(\lambda t)^p) \\[10px] &= e^{-(\lambda t)^p} \cdot (-\lambda^p) \cdot p t^{p-1} \\[10px] &= – \lambda p (\lambda t)^{p-1} \cdot e^{-(\lambda t)^p} \end{aligned}$$であるから、
$$\begin{aligned} (f(t)=) \lambda p (\lambda t)^{p-1} \cdot e^{-(\lambda t)^p} = -\frac{d}{dt} e^{-(\lambda t)^p} \end{aligned}$$である。
よって、
$$\begin{aligned} F(t) &= \int_{0}^{t} f(s) ds \\[10px] &= \int_{0}^{t} -\frac{d}{ds} e^{-(\lambda s)^p} ds \\[10px] &= -[e^{-(\lambda s)^p}]_{0}^{t} \\[10px] &= 1-e^{-(\lambda t)^p} \\[30px] S(t) &= 1-F(t) \\[10px] &= e^{-(\lambda t)^p} \\[30px] h(t) &= -\frac{d}{dt} \log S(t) \\[10px] &= -\frac{d}{dt} (-(\lambda t)^p) \\[10px] &= \lambda^p p \cdot t^{p-1} \\[30px] H(t) &= \int_{0}^{t} h(s) ds \\[10px] &= \int_{0}^{t} \lambda^p p \cdot s^{p-1} ds \\[10px] &= \lambda^p p [\frac{s^p}{p}]_{0}^{t} \\[10px] &= (\lambda t)^p \end{aligned}$$
ハザード関数$h(t) = \lambda^p p \cdot t^{p-1}$の形から$h(t)$は、
⚫︎$p > 1$：単調増加
⚫︎$p = 1$：一定
⚫︎$p < 1$：単調減少
となります。
病気や治療介入の種類によってはこの様なモデルがよく当てはまります。
3. 打ち切り(censoring)
実際にデータを収集してみると、データの完全な観測が困難となることがありえます。
その様な場合を『打ち切りが生じている』と表現します。
例.
患者さんを3年間フォローアップする臨床研究において、患者さんが引っ越したためフォローアップから外れてしまった場合、それ以降の状態は不明となります。（右側打ち切り）
生存時間データにはしばしば打ち切りが生じるため、打ち切りを考慮した尤度関数が必要となります。 以下それを求めていきますが、やや発展的なので初見の方は飛ばしてもよいかと思います。
3-1. 生存時間データの尤度関数
以下の様に定式化します。
・$X_i$：生存時間（$\theta$：$X_i$の分布に関するパラメータ）
・$V_i$：打ち切り時間
・$T_i = \min(X_i, V_i)$ （右側打ち切りを想定）
・$C_i = \vec 1 \{ V_i \leqq X_i \}$ （打ち切りがなければ$0$、あれば$1$）
また仮定として、
・$\{ (T_i, C_i) \}_i$：ある同時密度関数に従い$i$について独立
・$X_i \indep V_i$
とします。
病院でデータを収集していたとして、状態が悪い方（死亡リスクの高い方）が継続的に受診することができずに打ち切りになってしまっていたとします。この時、打ち切りにならない方となる方で死亡リスクが異なってしまいます。上記の仮定は『打ち切りは死亡リスクによらずランダムに生じる』というものなので、こういった状況ではこの仮定は破られています。
なお$(T,C)$の同時密度関数を$f_{T, C}$と表し、$X, V$の確率密度関数、累積分布関数をそれぞれ$f_X(\bullet ; \theta), f_V(\bullet), F_X(\bullet ; \theta), F_V(\bullet)$と表します。
すると観測データは$(T_i, C_i)$となり、この同時密度関数をパラメータ$\theta$に関する関数とみると尤度関数となります。
先に結論を書くと、尤度$L(\theta)$は、
$$\begin{aligned} L(\theta) &= \prod_{i=1}^{n} f_{T, C} (t_i, c_i; \theta) \\ &\propto \prod_{i=1}^{n} f_X (t_i; \theta)^{1-C_i} S(t_i; \theta)^{C_i} \\ &{\scriptsize(\thetaについて比例)} \end{aligned}$$となります。
打ち切りが生じなければ確率密度関数$f_X (t_i; \theta)$を尤度の計算に用いており、これは自然ですね。一方で打ち切りを生じた場合には生存関数$S(t_i; \theta)$を尤度の計算に持ちいており、これは今までと違うところですね。
（証明）
$Pr\{ T \in A, C = 0 \}, Pr\{ T \in A, C = 1 \}$について、
$$\begin{aligned} Pr\{ T \in A, C = 0 \} &= Pr\{ X \in A, X \leqq V \} \\[10px] &= \int_{A} ( \int_{x}^{\infty} f_V(v) dv ) f_X (x ; \theta) dx \\[10px] &= \int_{A} ( 1-F_V(x) ) f_X (x ; \theta) dx \\[20px] Pr\{ T \in A, C = 1 \} &= Pr\{ V \in A, V \lt X \} \\[10px] &= \int_{A} ( \int_{v}^{\infty} f_X (x ; \theta) dx ) f_V (v) dv \\[10px] &= \int_{A} ( 1-F_X (v ; \theta) ) f_V (v) dv \\[10px] &= \int_{A} S (v ; \theta) f_V (v) dv \end{aligned}$$となる。
よって、まとめると同時密度関数$f_{T, C}$は、
$$\begin{aligned} f_{T, C} (t, c) &= [ (1-F_V (t)) f_X (t ; \theta) ]^{1-C} [ S(t ; \theta) f_V (t) ]^{C} \end{aligned}$$となるが、$\theta$：$X_i$の分布に関するパラメータ、であるため、
$$\begin{aligned} f_{T, C} (t, c; \theta) &= [ (1-F_V (t)) f_X (t; \theta) ]^{1-C} [ S(t; \theta) f_V (t) ]^{C} \\[10px] &\propto f_X (t; \theta)^{1-C} S(t; \theta)^{C} \\ &{\scriptsize(\thetaについて比例。ここではVについての分布パラメータに興味はない。)} \end{aligned}$$となっている。
以上より、$n$個の互いに独立な観測$(t_i, C_i)$が得られた時の尤度$L(\theta)$は、
$$\begin{aligned} L(\theta) &= \prod_{i=1}^{n} f_{T, C} (t_i, c_i; \theta) \\[10px] &\propto \prod_{i=1}^{n} f_X (t_i; \theta)^{1-C_i} S(t_i; \theta)^{C_i} \\ &{\scriptsize(\thetaについて比例)} \end{aligned}$$となる。
まとめ.
- 生存時間まわりの5種類の関数の関係性を確認した。いずれか1種類が定まれば残り4種類は自動的に決まる。
- 打ち切りを考慮した時の尤度は、生存時間の確率密度関数と生存時間関数を以て表現される。
