生存曲線が群間で等しいかどうかを検定するときにはLog-rank検定が使われる。
これは「群間でハザード関数が等しい」という帰無仮説のもとで期待されるイベント数と、実際観察されたイベント数が乖離しているかどうかを検定する。
ここではsurvivalパッケージのveteranデータセットを改変して使う。
library(survival) library(tidyverse) d <- veteran %>% mutate(age_cat = if_else((age<60),"under60","over60"), age_cat = factor(age_cat), trt = factor(trt, labels=c("standard","test")), prior = factor(prior, labels=c("N","Y")))
> head(d,3) trt celltype time status karno diagtime age prior age_cat 1 standard squamous 72 1 60 7 69 N over60 2 standard squamous 411 1 70 5 64 Y over60 3 standard squamous 228 1 60 3 38 N under60
survidiff()を使ってLog-rank検定を行う
Surv(time, event)で生存時間データを作って従属変数とする。
> survdiff(Surv(time, status) ~ trt, data=d) Call: survdiff(formula = Surv(time, status) ~ trt, data = d) N Observed Expected (O-E)^2/E (O-E)^2/V trt=standard 69 64 64.5 0.00388 0.00823 trt=test 68 64 63.5 0.00394 0.00823 Chisq= 0 on 1 degrees of freedom, p= 0.9
デフォルトだとP値の桁数が少ないので、必要ならprint()で桁数を指定する（指定した桁数より2つほど少なくなる理由は不明）。
> print(survdiff(Surv(time, status) ~ trt, data=d), digit=10) Call: survdiff(formula = Surv(time, status) ~ trt, data = d) N Observed Expected (O-E)^2/E (O-E)^2/V trt=standard 69 64 64.50019666 0.003879006813 0.008227343202 trt=test 68 64 63.49980334 0.003940117750 0.008227343202 Chisq= 0 on 1 degrees of freedom, p= 0.92772723
Log-rank検定以外の方法を使う
Log-rank検定は比例ハザード性が成立しているときに最も検出力が高くなるが、ハザード比が多少時間で変化しても構わない。
どの期間の群間差に注目しているかに応じて重み付けをしたLog-rank検定（weighted Log-rank test）を行うこともできる。
観察期間中ハザードが一定にならない状況の例としては、
- 治療の効果が短期間で、時間の経過とともに有益性が薄まる
- 積極的治療では初期は害が目立っても後半で有益性がみられる
- 保存的治療ではある程度の時間が経過してようやく有益性がみられる
などがある。
代表的な手法と重みの付け方は次のとおり。
- Wilcoxon検定：イベントが発生したときのnumber at risk
に応じて重み付け
- Taron-Ware検定：
に応じて重み付け
- Peto検定：イベントが発生した時の生存確率のKaplan-Meier推定値
に応じて重み付け
- Fleming-Harrington検定（
-class）*1：
に応じて重み付け
- Harrington-Fleming検定（
-class）*2：
に応じて重み付け
型Fleming-Harrington検定では、
- 前半部分に重みを付けたければ
- 中盤部分に重みを付けたければ
- 後半部分に重みを付けたければ
- 全体に均等に重みを付けたければ
（Log-rank検定と同じ）
となり自由度が高い。
survidiff()ではrhoを指定することで 型Harrington-Fleming検定を行うことができる。
rho=0：通常のLog-rank検定rho=1：Peto-Peto検定（Gehan-Wilcoxon検定のPeto-Peto修正版）
> survdiff(Surv(time, status) ~ trt, data=d, rho=1) # Peto-peto Call: survdiff(formula = Surv(time, status) ~ trt, data = d, rho = 1) N Observed Expected (O-E)^2/E (O-E)^2/V trt=standard 69 32.2 35.4 0.279 0.871 trt=test 68 35.2 32.1 0.308 0.871 Chisq= 0.9 on 1 degrees of freedom, p= 0.4
P値の桁数を増やしたいときはprint()を使う。
層別Log-rank検定を行う
交絡調整の手段として層別化がある。
「対象をいくつかの層に分け、層内で群間比較し、結果を統合する」という方法で、Log-rank検定においても層別解析することが可能。
層別変数をstrata()で指定すればよい。
> survdiff(Surv(time, status)~ trt + strata(age_cat), data=d) Call: survdiff(formula = Surv(time, status) ~ trt + strata(age_cat), data = d) N Observed Expected (O-E)^2/E (O-E)^2/V trt=standard 69 64 64.6 0.00610 0.0131 trt=test 68 64 63.4 0.00622 0.0131 Chisq= 0 on 1 degrees of freedom, p= 0.9
層別化はweighted Log-rank testへも適用できる。
3つ以上の群間で生存曲線を比較する
単に3つ以上のカテゴリーを持つ変数を使うだけ。
strata, rhoについても同様に指定できる。
> survdiff(Surv(time, status) ~ celltype + strata(age_cat), rho=1, data=d) Call: survdiff(formula = Surv(time, status) ~ celltype + strata(age_cat), data = d, rho = 1) N Observed Expected (O-E)^2/E (O-E)^2/V celltype=squamous 35 13.61 20.1 2.10 4.66 celltype=smallcell 48 28.46 19.2 4.43 8.88 celltype=adeno 27 16.08 11.0 2.35 3.89 celltype=large 27 9.78 17.6 3.47 7.33 Chisq= 18.6 on 3 degrees of freedom, p= 3e-04
おわりに
- 迷ったときはLog-rank検定でよい
参考資料
- 生存時間解析について自己学習できる詳細な解説付きテキスト
*1:Fleming Thomas R, Harrington David P. A class of hypothesis tests for one and two sample censored survival data Communications in Statistics-Theory and Methods. 1981;10:763–794.
*2:Harrington, D. P. and Fleming, T. R. (1982). A class of rank test procedures for censored survival data. Biometrika69, 553-566.
