はじめに

業務上、生存時間分析を扱うこと~~になった~~ したので、少しづつ勉強しています。
今回は良い機会なので、生存時間分析の基本的な事項と、生存確率を推定する手法であるKaplan−Meier推定法およびその検定までをまとめてみました。

生存時間分析とは

ある時点から、死亡や故障、退会などのイベントが起こるまでの時間を推定する分析手法です。
その名の通り投薬の延命効果や、疾患発症後の余命の推定などに用いられるようですが、工学分野でも信頼性工学として故障までの時間の推定に用いられています。
また、マーケティング分野ではサブスクリプションビジネスにおけるチャーンの予測にも用いられているようです。

ここで気をつけるべきは、生存時間分析では、そのイベントが起こるかどうかではなく、イベントが起こるまでの時間、すなわち生存時間をモデリングしている点です。

生存時間

生存時間は、観測開始した時点から、注目しているイベントが発生した時点までの時間として定義されます。会員性サービスの場合はユーザが利用を開始した時点から解約した時点までの時間とすることができますし、医療へ応用する場合は、ある患者に投薬した時点から死亡するまでの時間を生存時間として定義できます。
これら時間は確定的なものではなく分布するので、確率変数として扱います。

打ち切りについて

生存時間分析では、打ち切りという概念が登場します。打ち切りとは、イベントが発生する前に何らかの理由でデータが取得できなくなってしまった場合のことを指します。例えば肺がんによる死亡に関するデータの場合は、違う病院に移った、事故など明らかに他の原因で死亡した、などがそれに該当します。

打ち切りの種類

打ち切りにはいくつか種類があるので、ここにまとめておきます。

右側打ち切り(right censoring)
観測期間が終了した時点でイベントがまだ発生していない場合
左側打ち切り(left censoring)
観測期間開始より前にすでにイベントが発生している様な場合
区間打ち切り
観測期間中のうち、ある期間内にイベントが発生していることはわかっているが、正確な時刻はわからない場合

マーケティング用途で生存時間分析を行う場合、右側打ち切り以外はほとんど考える必要は無いと思います。

生存関数

生存時間がどのように分布しているかは、生存関数を用いて表します。

生存時間
確率変数としての生存時間を記号 $T$ として定義
累積分布関数:
無作為抽出された観測対象の生存時間が時間t以下である確率

F(t) = {\rm Pr}(T \le t)

生存関数:
ある観測対象について時間 $t$ よりも長い生存時間が得られる確率

S(t) = {\rm Pr}(T > t)

任意の$t$について、累積分布関数と生存関数の和は1となります。つまり$F(t) = 1- S(t)$が成り立ちます。
ほとんどの場合、イベントがどれだけ早く起こるかよりも、生存時間がどれだけ長くなるか、つまり解約や死亡するタイミングがどれだけ遅いかに興味があることが多いため、以降では生存関数の推定に焦点を当てます。

Kaplan-Meier推定法

Kaplan-Meier推定法は、打ち切りを含む全ての観測データを利用して生存関数を推定出来る手法であり、ほとんどの統計ソフトウェアで算出できる標準的な推定方法のようです。この推定法では、各時刻における生存確率をそれぞれ計算し、それらを掛け合わせることで全体の生存関数を推定します。

Kaplan-Meier推定法の定義

時刻$t$における生存確率$\hat S(t)$は次のように定義される。

\hat S(t) = \prod_{t_i \le t} (1- \frac{d_i}{n_i})

ここで$n_{i}$は時刻$t_i$においてイベントが発生しうる人数を表します。$d_{i}$は$n_i$のうちイベントが実際に発生した人数を表します。$n_i$が表している観測対象の集合はリスク集合と呼ばれる場合もあります。

具体的な例として、サブスクリプションの音楽サービスについて考えてみます。
10月から1月までの4ヶ月間観測をしているとき、合計で5人のユーザがそれぞれのタイミングで利用を開始しました。その間にユーザBとユーザDは解約してしまいましたが、他のユーザは1月でも継続して利用しています。

Kaplan−Meier推定法を適用する際は、それぞれのユーザの利用開始月を0ヶ月目として、下図のように揃えて考えます。ユーザBとユーザDはそれぞれ利用開始から1ヶ月目と2ヶ月目に解約したことが分かっています。また、ユーザEは1月から４月までの4ヵ月利用利用していたことが分かっています。一方ユーザAは11月から、ユーザCは1月kから利用を開始しており、2月以降(それぞれ利用開始から3ヶ月目以降、1ヶ月目以降)の利用状況はわかりません。この場合これらのユーザは打ち切りとして扱います。

この例において2ヶ月目の生存関数の値$\hat S(2)$を求めてみると、次のように計算することができます。

\begin{eqnarray}
 \hat S(2) & = &   \prod_{t_i \le 2} (1- \frac{d_i}{n_i})\\
   & = & (1 - \frac{d_0}{n_0}) (1-\frac{d_1}{n_1})(1-\frac{d_2}{n_2})\\
& = & (1- \frac{0}{5})(1-\frac{1}{4})(1-\frac{1}{3})\\
& = & 1 \times \frac{3}{4} \times \frac{2}{3}\\
& = & 0.5
\end{eqnarray}

生存関数の比較

投薬した群とそうでない群、もしくは年代や地域別の生存関数の間に差があるかどうかを確かめるためには検定を行う必要があります。Kaplan-Meier推定法で時刻ごとの生存率を掛け合わせて生存関数を求めたように、検定統計量を求める場合も、各時刻の統計量の和を計算します。
それぞれの検定統計量は時刻ごとの、次のような分割表に基づいて計算されます。ここで添え字のiは時刻が$t_i$であることを表しています。

イベント\群	1	0	合計
発生	$d_{1i}$	$d_{0i}$	$d_i$
未発生	$n_{1i} - d_{1i}$	$n_{0i} - d_{0i}$	$n_i - d_i$
リスク集合	$n_{1i}$	$n_{0i}$	$n_i$

分割表の差を検定する場合は$\chi^2$検定を用いますが、生存関数の差を検定するときも$\chi^2$検定と同じ様な流れで計算します。

計算手順
1. 期待度数を求める
2. データと期待度数の差を求める
3. 統計量Qを求める
4. 統計量Qが自由度1の$\chi^2$分布に従うことを利用してp値を求める
5. p値を解釈する

期待度数を求める
期待度数とは、各群でイベントの発生率に差がないとしたとき、きっとこうなるだろうというという数値を指します。この期待度数とデータの間に大きな乖離があったとしたら、各群の間の生存確率に差がありそうだとみなせそうです。
ここではそれぞれの期待度数を次のようにおいてみます。

イベント\群 1 0 合計

発生 $\hat e_{1i}$ $\hat e_{0i}$ $d_i$

未発生 $\hat e'_{1i}$ $\hat e'_{0i}$ $n_i - d_i$

リスク集合 $n_{1i}$ $n_{0i}$ $n_i$

ただし、すべての期待度数を求める必要はなく、ここで言う$\hat e_{1i}$もしくは$\hat e_{0i}$のどちらかを求めればよいです。

例として$\hat e_{1i}$を求めてみると次のようになります。
```
\begin{eqnarray}
\hat e_{1i} &=& n_{1i} \frac{n_i}{d_i}
\end{eqnarray}
```
データと期待度数の差およびイベント発生数の分散を求め、時刻について総和を取る
次に,$d_{1i}$と$\hat e_{1i}$の差$O_i$と、$d_{1i}$の分散$V_1$をもとめ、時刻について総和を取ります。$O_i$は次のように表すことができます。
```
O_1 = \sum_{i=1}^{m}(d_{1i} - \hat e_{1i})
```
全体のイベント発生数およびリスク集合が定数であると仮定すると$d_{1i}$は超幾何分布に従います。このことを利用すると、その分散 $\hat v_1$は
```
\hat v_1 = \frac{n_{1i}n_{0i} d_i(n_i - d_i) }{n_i^2(n_i - 1)}
```
で求める事ができます。これを用いると$V_1$は
```
V_1 = \sum_{i=1}^{m} \hat v_1
```
となります。
統計量$Q$をもとめる
求めた$O_1$と$V_1$を用いて統計量Qを計算します。
```
Q = \frac{O_1^2}{V_1}
```
統計量$Q$が自由度1の$\chi^2$分布に従うことを利用してp値を求める
p値を求めます。
手順3で求めた$Q$は自由度1の$\chi^2$分布に従うようですので、$p$値は次のように求められます。
```
p = \rm{Pr} \left(\chi^2(1) \ge Q \right) 
```
p値を解釈する

伝統的にはp値が0.05もしくは0.01を下回れば有意差があるとみなすことが多いようです。
ただし、p値が大きかったからと言って差がないと断定することはできませんし、p値が小さければ万事OKということでも無いことには注意が必要です。

ここで説明した検定方法をログランク検定(log-rank test)と呼びます。

終わりに

統計ムズカシイ

参考文献

http://www.math.s.chiba-u.ac.jp/~wang/suvival.pdf
https://www.ism.ac.jp/editsec/toukei/pdf/54-2-445.pdf
https://qiita.com/saltcooky/items/409329485be499a5b270
http://racco.mikeneko.jp/Kougi/10a/IS/IS07pr.pdf
https://note.com/maxwell/n/nc78c55afe944#DVcZz
https://logics-of-blue.com/chi-squared-test/