生存時間分析 ログランク検定

f:id:kimoppy126:20201011181656p:plain

2群以上のデータを比較する際に、最も簡単な方法として2つの生存関数を同じ軸の上にプロットすることが挙げられます。

from matplotlib import pyplot as plt
from lifelines.datasets import load_dd
from lifelines import KaplanMeierFitter

data = load_dd()
data.head()
T = data["duration"]
E = data["observed"]
ax = plt.subplot(111)

# 民主主義か否か
dem = (data["democracy"] == "Democracy")

# カプラン・マイヤー推定量を使用
kmf = KaplanMeierFitter()
# 民主主義データに関して生存関数を推定
kmf.fit(T[dem], event_observed=E[dem], label="Democratic Regimes")
kmf.plot(ax=ax)
# 非民主主義データに関して生存関数を推定
kmf.fit(T[~dem], event_observed=E[~dem], label="Non-democratic Regimes")
kmf.plot(ax=ax)

plt.title("Lifespans of different global regimes");

この際、2群の生存関数の推定値に差が見られた場合の解釈は2つ考えられます。

  1. 2群の間には本当に意味のある差があって、それゆえ両群の生存関数に差が見られた
  2. 2群の間には本当に意味のある差は存在せず、観測された差は単なる偶然変動である

これらの2つの解釈を区別するための手法が「仮説検定」となります。

仮説検定とは

観測されたデータと帰無仮説と呼ばれる特定の仮説がどの程度一致しているかを評価する手順です(例:「2群の間に差はない」、「生存時間と説明変数間に関係がない」等)。

仮説検定の流れ

  1. 帰無仮説及び対立仮説、有意水準を定義する
  2. 観測データから検定統計量を構築する
  3. 帰無仮説が正しい際に検定統計量がその値と同じかそれより極端な値が得られる確率(P-value)を求める
  4. P-valueの大きさが有意水準より大きいか小さいかによって帰無仮説を棄却できるか否かを判断する

ログランク検定

ログランク検定は、死亡または別の重要なイベントである可能性があるイベントについて、2つのサンプルからの生存分布を比較するためのノンパラメトリック仮説検定です。
2つのグループの個人の生存経験を比較するために、臨床試験でよく使用されます。

たとえば、薬物Aが薬物Bと比較して生存期間を延長するかどうかをテストするとします。
このテストでは、違いがあるかどうかがわかりますが、違いのサイズに関する詳細や信頼区間は提供されません。
2標本のt検定等の他の有意差検定では、時系列を跨いで差異を考慮することはできないため、生存時間分析ではログランク検定が用いられます。

ログランク検定では、下記の3つの強力な仮定を置いています。

  • 打ち切り(正確な生存期間がわからない場合に起こります)がイベントとは無関係に生じる
  • 研究対象の被験者について、観察開始時点での生存確率が等しい
  • イベントが特定の時間に生じる

イベントの発生時刻を $\tau_1 \lt\tau_2 \lt ... \lt \tau_k$ とします。

$Y_i (\tau_j)$ = $\tau_j$ におけるグループ$i$ のエクスポージャ数

$Y (\tau_j)$ = $Y_0 (\tau_j) + Y_1 (\tau_j)$

$d_{ij}$ = $\tau_j$ におけるグループ$i$ のアウトカム数

$d_{ij}$ = $d_{0j} + d_{1j}$

帰無仮説 $H_0$: 2群の生存時間の観測値に差がない($h_0(\cdot) = h_1(\cdot)$)

帰無仮説 $H_0$ の下で、死亡率 $P[d_{j} = d]$ は超幾何分布に従うため、期待値($E_j$)、分散($V_j$)は下記の式で表せます。

$$ \begin{align} E_j &= \frac{Y_1(\tau_j)}{Y(\tau_j)} d_j \\ V_j &= E_j \frac{Y(\tau_j) - d_j}{Y(\tau_j)}\frac{Y(\tau_j) - Y_1(\tau_j)}{Y(\tau_j) - 1} \\ &= \frac{Y_0(\tau_j) Y_1(\tau_j) d_j (Y(\tau_j) - d_j)}{Y(\tau_j)^2 (Y(\tau_j) - 1)} \\ \end{align} $$ ログランク検定では、観測値の期待値からの逸脱度の指標として、検定量 $Z$ が標準正規分布に、検定量 $Z^2$ は自由度1の$\chi^2$分布に従うことを検定します。 $$ Z = \frac{\sum^{K}_{j=1} (O_{i, j} - E_{i, j})}{\sqrt{\sum^{K}_{j=1} V_{i, j}}} $$

$O_j$: 各時点における観測された死亡者数

$E_j$: 各時点における死亡者数の期待値

この検定統計量は2群の順位(ランク)から導出されており、結果として得られる順位検定統計量がネルソン・アーレン推定量の対数(ログ)となっていることが検定名の由来となっています。

lifelinesでは、logrank_testとして実装されています。

from lifelines.statistics import logrank_test
# ログランク検定を実行
results = logrank_test(T[dem], T[~dem], E[dem], E[~dem], alpha=.99)
results.print_summary()

## <lifelines.StatisticalResult: logrank_test> ## t_0 = -1 ## null_distribution = chi squared ## degrees_of_freedom = 1 ## alpha = 0.99 ## test_name = logrank_test ## ## --- ## test_statistic p -log2(p) ## 260.47 <0.005 192.23

検定統計量($\chi^2$)、P-value、対数P-valueが出力されます。

加重ログランク検定

加重ログランク検定は、グループ間の有意差を比較したいが、特定のイベントをより重要(「加重」)としたい場合に使用されます。

通常のlog-rank検定では、すべての時刻における重要度が均一であると仮定していますが、これはハザード(死のリスク)が長期にわたって比較的一定であるか、または比例している場合に有効です。 加重ログランクテストでは、異なる時間のハザードの重み付けを行うことができます(例:生存の確率が時間の初めにはるかに大きく、終わりに次第に減少する場合など)

なお、重要度を表す重みは生存曲線からではなく、研究トピックに関する事前情報に基づいて選択する必要があります。
つまりデータセットを確認する前、分析設計の際に設定する必要があります。

加重ログランク検定の統計量は下記の通りとなります。

通常のログランク検定と同様に、観測値の期待値からの逸脱度の指標として、検定量 $Z$ が標準正規分布に、検定量 $Z^2$ は自由度1の$\chi^2$分布に従うことを検定します。

$$ Z = \frac{\sum^{K}_{j=1}w_j (O_{i, j} - E_{i, j})}{\sqrt{\sum^{K}_{j=1}w_j^2 V_{i, j}}} $$

$O_j$: 各時点における観測された死亡者数

$E_j$: 各時点における死亡者数の期待値

$w_j$: 各時点における重み

階層化ログランク検定

年齢、性別、体重などの変数を制御したい場合に使用します。たとえば、2つのグループと、性別(男性、女性)を制御して、4つの可能なグループ、つまり階層を与えることができます。

観測したいイベントの発症率にある説明変数の値が影響を与えている、もしくはある説明変数の値で分類した際にサンプル数が均一でない場合には層別ログランク検定が有用となります。

層$l$、グループ$j$の生存関数を$S_j^{(l)}$としたとき帰無仮説は下記になります。

帰無仮説 $H_0$:
全ての層について、2群の生存時間の観測値に差がない($S_0^{(l)}(\cdot) = S_1^{(l)}(\cdot) (l=1, ..., L)$)

層別ログランク検定では、層ごとにログランク検定を実行し、全ての統合的な検定量について有意性を判定します。

  1. データを $L$ グループに分割する(例:性別のみによって層別化する場合、$L=2$)
  2. 通常のログランク検定により、各層について $O,E,V$ を算出する
  3. 帰無仮説 $H_0$ の下で、下記の統計検定量$Z$が標準正規分布に従うことを検定する

検定統計量 $Z$ は下記となります。

$$ Z=\frac{\sum^L_{l=1}{(O^{(l)}-E^{(l)}})}{\sqrt{\sum^L_{l=1}{V^{(l)}}}} $$

参考