【統計】生存時間解析　ログランク検定、一般化Wilcoxon検定（Gehan's Wilcoxon検定）等

ログランク検定、重み付きログランク検定たちについてのメモ。

【目次】

計算式等
ノート
プログラムコード
参考

前回のKaplan-meier 法に続く内容。
cochineal19.hatenablog.com

計算式等

先に計算方法をまとめておく。

ログランク検定は生存時間解析において、2群間の生存時間の差を検定する手法。
時点 t ごとにイベントの観測度数と期待度数のズレを見ており、下のクロス表にまとめることができる。
（※クロス表の黄色セル部分を見ている。）

■ 時点 t ごとのイベント有無のクロス表

		イベント		合計
		有	無	合計
群	A	$d_{A}$	$n_{A} - d_{A}$	$n_{A}$
群	B	$d_{B}$	$n_{B} - d_{B}$	$n_{B}$
全体		$d_{all}$	$n_{all} - d_{all}$	$n_{all}$

■ ログランク検定（Log-Rank test）

$\quad \chi ^{2} = \dfrac{ \sum w_{i} \left( O_{i} - E_{i} \right)^{2}}{\sum w_{i}^{2}\ Var[O_{i} - E_{i} ]}$

$\qquad O_{i} = d_{A}$

$\qquad E_{i} = n_{A} \dfrac{d_{all}}{n_{all}}$

$\qquad Var[O_{i} - E_{i} ] = \dfrac{n_{A} \cdot n_{B} \cdot d_{all}\left( n_{all} - d_{all} \right)}{n_{all}^{2} \left(n_{all}-1\right)}$

この式はマンテル・ヘンツェルの統計量（CMH統計量）であり、ログランク検定はマンテル・ヘンツェルの方法を用いていることが分かる（Cox-Mantel検定、Mantel Log Rank検定とも呼ばれる）。
cochineal19.hatenablog.com

■ ログランク検定の重み
生存時間解析はイベント・打ち切りの発生により観察期間の後期になるほどN数が少なくなる特徴がある。
このため各時点に重みを与える方法がいくつか考案されている。

手法	重み $w_{i}$	補足
Log-Rank (Mantel)	1	一定の重み。
Gehan's Wilcoxon	$n_{all}$	初期に強い重み。
Tarone-Ware	$\sqrt{n_{all}}$	初期に強い重み。
Peto-Peto	$\tilde{S} _{\left( t \right)} = \prod _{t_{i}＜t} \left(1-\dfrac{d_{all}}{n_{all}+1} \right)$	初期に少し重み。
Harrington-Fleming (a, b)	$\tilde{S} \left( t_{i-1} \right)^{a} [ 1-\tilde{S} \left( t_{i-1} \right) ]^{b}$	初期に少し重み。

f:id:cochineal19:20210328114133p:plain:w500

重みを見ての通り、観察期間の前半に差が出やすいものは重み付きのログランク検定の方が有意差が出やすく、後半に差が出やすいものは重みなしのログランク検定の方が有意差が出やすい。

簡単な歴史
* Cox (1953) : 二群の差をF検定で示す方法を提案（指数分布に従う場合）。
* Mantel and Haenszel (1959), Mantel (1966) : Coxの手法を拡張して、ログランク検定（重みなし）を提案。
* Gehan (1965) : Gehan-Wilcoxon検定を提案。
* Peto and Peto (1972) : Peto-Peto検定を提案。
* Tarone and Ware (1977) : Tarone-Ware検定を提案。
* Fleming and Harrington (1981), Harrington and Fleming (1982) : Harrington-Fleming検定を提案。

ノート

今回はSAS中心で記事を書く。（Rの関数は調べれていない。。。）

サンプルデータ

前回のKaplan-meier法の記事と同じく、RのMASSパッケージにあるgehanデータ（Remission Times of Leukaemia Patients）を用いる。

まずは手計算してみる

ログランク検定をExcelで計算してみると下図の通りになる。
（累積生存率 S(t) の計算は前回のKaplan-meier法の記事を参照。）

Log-rank (Mantel)
f:id:cochineal19:20210324000945p:plain

Gehan's Wilcoxon
f:id:cochineal19:20210324001220p:plain

Peto-Peto
f:id:cochineal19:20210324001336p:plain

公式どおりなので特記はないが、1行目のデータ（時点t=1）について式を追ってみる。

		イベント		合計
		有	無	合計
群	A	$d_{A}=0$	$n_{A} - d_{A}=21$	$n_{A}=21$
群	B	$d_{B}=2$	$n_{B} - d_{B}=19$	$n_{B}=21$
全体		$d_{all}=2$	$n_{all} - d_{all}=40$	$n_{all}=42$

時点t=1は、A群：イベント＝０例、B群：イベント＝２例である。
よって、

$\quad O_{1} = d_{A}=0$

$\quad E_{1} = n_{A} \dfrac{d_{all}}{n_{all}}= 21 \dfrac{2}{42} = 1$

$\quad Var[O_{1} - E_{1} ] = \dfrac{n_{A} \cdot n_{B} \cdot d_{all}\left( n_{all} - d_{all} \right)}{n_{all}^{2} \left(n_{all}-1\right)}=\dfrac{21 \times 21 \times 2\left( 42 - 2 \right)}{42^{2} \left(42 - 1\right)}=0.4878$

ここに重みを与えれば各手法となる：

Log-rank (Mantel)

$\quad w_{1}=1$

$\quad w_{1}\left(O_{1} - E_{1}\right) = 1 \times \left(0 - 1\right) = -1$

$\quad w_{1}^{2} Var[O_{1} - E_{1} ] = 1^{2} \times 0.4878 = 0.4878$

Gehan's Wilcoxon

$\quad w_{1}=n_{all,1}=42$

$\quad w_{1}\left(O_{1} - E_{1}\right) = 42 \times \left(0 - 1\right) = -42$

$\quad w_{1}^{2} Var[O_{1} - E_{1} ] = 42^{2} \times 0.4878 = 860.4878$

Peto-Peto

$\quad w_{1} = \tilde{S} _{\left( 1 \right)} = 1 - \dfrac{d_{all}}{n_{all}+1} = 1 - \dfrac{2}{42+1} = 0.953$

$\quad w_{1}\left(O_{1} - E_{1}\right) = 0.953 \times \left(0 - 1\right) = -0.953$

$\quad w_{1}^{2} Var[O_{1} - E_{1} ] = 0.953^{2} \times 0.4878 = 0.4435$

以上の処理を時点分行ってカイ二乗統計量を算出する。

SASでの実行

LIFETESTプロシジャのSTRATAステートメントのTESTオプションで指定する。

proc lifetest data=ADTTE plots=s(cl atrisk maxlen=30 outside) alpha=0.5 conftype=loglog;
    time AVAL * CNSR(0);
    strata TRT01A / test=(LOGRANK WILCOXON TARONE PETO FH(0.5,0.5));
    ods output SurvivalPlot = SP1;
    ods output Quartiles    = QT1;
run;

f:id:cochineal19:20210324002059p:plain:w300

Rでの実行

Log-rank (Mantel) についてsurvivalパッケージのsurvdiff関数で実行できる。
（※他の手法は調べれていない）

> library(survival)
> survdiff(Surv(time=AVAL, event=EVENT) ~ TRT01A, data=ADTTE)
Call:
survdiff(formula = Surv(time = AVAL, event = EVENT) ~ TRT01A, 
    data = ADTTE)

                N Observed Expected (O-E)^2/E (O-E)^2/V
TRT01A=6-MP    21        9     19.3      5.46      16.8
TRT01A=control 21       21     10.7      9.77      16.8

 Chisq= 16.8  on 1 degrees of freedom, p= 4e-05

プログラムコード

library(survival)
survdiff(Surv(time=AVAL, event=EVENT) ~ TRT01A, data=ADTTE)

SAS

proc lifetest data=ADTTE plots=s(cl atrisk maxlen=30 outside) alpha=0.5 conftype=loglog;
    time AVAL * CNSR(0);
    strata TRT01A / test=(LOGRANK WILCOXON TARONE PETO FH(0.5,0.5));
    ods output SurvivalPlot = SP1;
    ods output Quartiles    = QT1;
run;

Python