統計ER

R, EZR, SPSS, KH Coder を使ったデータ分析方法を紹介するブログ。ニッチな内容が多め

ログランク検定はどんな計算をしているのか?

生存時間解析に登場する群間比較検定のログランク検定。

いったいどんな計算をしているのか?

生存時間自体は使っていないことを知っているだろうか?

ログランク検定とは?

ログランク検定とは、以下のようなカプランマイヤー曲線で描かれるような生存時間データの群間比較する検定である。

f:id:toukeier:20220321082911p:plain

実際どのような計算をしているのだろうか?

ログランク検定はどのような計算をしているのか?

ログランク検定の帰無仮説は母集団におけるイベントの観測値と期待値に差がないというものである。

検定統計量は、イベントが起きた時点ごとの観測値と期待値の差分の合計を、時点ごとの分散の合計の平方根、つまり標準偏差的な数値で割ったものになる。

添え字等を省くと以下のような式になる。

 \displaystyle Z = \frac{\sum {(O - E)}}{\sqrt {\sum {V}}}

詳しくはWikipedia参照。

ログランク検定 - Wikipedia

いわゆるカイ2乗検定の検定統計量にとても似ている。

こちらも添え字等を省くと以下のような式になる。

 \displaystyle  \chi^{2}= \sum {\frac{(O - E)^{2}}{E}}

EZRでは、カイ2乗値として表示され、上記の検定統計量の2乗と一致する。

本稿の最後あたりで確認するのでお楽しみに。

カイ2乗検定の詳細はWikipedia参照。

カイ二乗検定 - Wikipedia

ログランク検定は生存時間自体を使っていない!?

上記のログランク検定の検定統計量の主体はイベントの観測値、つまり発生したイベント数である。

検定統計量の式で省略した添え字は、イベントが観測された時点の順番と、グループの番号である。

つまり、検定統計量に「時間」そのものは入っていない。

生存時間の平均値の差のようなものを計算しているのだと思っていたのなら、それは誤解だ。

生存時間解析ではあるが、生存時間そのものは使っていない。

ログランク検定はどのような計算をしているのか?

では、ログランク検定はどのような計算をしているのだろうか?

サンプルデータを用いて確認しよう

サンプルデータを用いて計算を再現してみようと思う。

サンプルデータは MASS パッケージの gehan データセットだ。

白血病の患者さんを6-メルカプトプリン(6-MP)で治療した群とコントロール群とで比較した試験のデータである。

上記のカプランマイヤー曲線は、このデータをプロットしたものである。

ログランク検定をしてみると以下のような結果が出力される。

EZRで、「統計解析」→「生存期間の解析」→「Logrank検定」を選択して以下のように指定する。

f:id:toukeier:20220321090517p:plain

ログランク検定の結果は以下のとおりである。

f:id:toukeier:20220321090630p:plain

> (res <- survdiff(Surv(time,cens==1)~treat, data=gehan, rho=0, na.action = 
+   na.omit))
Call:
survdiff(formula = Surv(time, cens == 1) ~ treat, data = gehan, 
    na.action = na.omit, rho = 0)

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

 Chisq= 16.8  on 1 degrees of freedom, p= 0.00004 

データの加工

gehanはこのようなデータである。

pairになっているが、今回はpairの要素は考慮しない。

timeがイベントもしくは打ち切りまでの観察時間(週)である。

censが1の時はイベント(死亡)発生、0は打ち切り(生存のまま観察終了)である。

treatが6-MP治療群かコントロール群かである。

f:id:toukeier:20220321093608p:plain
gehan dataset

先ほどカプランマイヤー曲線を描かせて、ログランク検定を行った結果を利用して、時点ごとのイベント数、打ち切り数をまとめたデータセット gehan.O12 を作成する。

gehan.O12 <- with(km, data.frame(treat=c(rep("6-MP",16),rep("CTRL",12)),time,n.risk,n.event,n.censor))

出来上がりがこちら。

timeがイベントもしくは打ち切り発生までの時間、n.riskがtime時点に試験に残っている各群の患者数、n.eventがtime時点でイベントが発生した人数、n.censorが打ち切りとなった人数である。

f:id:toukeier:20220321094647p:plain

この6-MP群とコントロール群群が縦に積み重なっているものを横に並べたいのだが、EZRでは実施方法がわからず、いったんCSV形式で出力し、エクセルで以下のように加工して、EZRに再度取り込んだ。

f:id:toukeier:20220321094941p:plain

ここで、timeは両群通してイベントか打ち切りが発生した時点までの時間、N1jはtime時点に試験に残っている6-MP群の患者数、O1jはtime時点でイベントが発生した人数、N2jはtime時点のコントロール群の患者数、O2jはtime時点の打ち切り例の人数である。

このデータをもとに以下のような加工をすると、期待値Eと分散Vが計算できる。添え字は各群(1が6-MP群、2がコントロール群)を示している。

#####計算式を入力して新たな変数を作成する#####
gehan_for_logrank$Nj <- with(gehan_for_logrank, N1j+N2j)
#新しい変数 Nj を作成しました。
#####計算式を入力して新たな変数を作成する#####
gehan_for_logrank$Oj <- with(gehan_for_logrank, O1j+O2j)
#新しい変数 Oj を作成しました。
#####計算式を入力して新たな変数を作成する#####
gehan_for_logrank$E1j <- with(gehan_for_logrank, N1j*Oj/Nj)
#新しい変数 E1j を作成しました。
#####計算式を入力して新たな変数を作成する#####
gehan_for_logrank$E2j <- with(gehan_for_logrank, N2j*Oj/Nj)
#新しい変数 E2j を作成しました。
#####計算式を入力して新たな変数を作成する#####
gehan_for_logrank$V1j <- with(gehan_for_logrank, 
  E1j*(Nj-Oj)/Nj*(Nj-N1j)/(Nj-1))
#新しい変数 V1j を作成しました。
#####計算式を入力して新たな変数を作成する#####
gehan_for_logrank$V2j <- with(gehan_for_logrank, 
  E2j*(Nj-Oj)/Nj*(Nj-N2j)/(Nj-1))
#新しい変数 V2j を作成しました。

f:id:toukeier:20220321095533p:plain

詳しい計算方法は、Wikipedia参照。

ログランク検定 - Wikipedia

検定統計量Zの計算とp値

ここまで準備が出来たら検定統計量Zを計算する。

> Z <- with(gehan_for_logrank, sum(O1j-E1j)/sqrt(sum(V1j, na.rm=T)))

> Z
[1] -4.097919

このZ値を用いて、標準正規分布を使ってp値を計算する。

pnorm()で、下側の確率が計算できるので、2倍して両側検定とする。

結果としてp値は0.00004となった。

> 2*pnorm(Z)
[1] 0.00004168809

これは先ほどのEZRのメニューからログランク検定を行ったときのp値と同じである。

f:id:toukeier:20220321090630p:plain

ちなみにZを2乗すると、カイ2乗値 Chisq=16.8 と一致する。

> Z^2
[1] 16.79294

これでログランク検定がどのような計算をしているかの全貌がわかった。

まとめ

生存時間解析で、群間を比較する方法として有名なログランク検定は、どのような検定で、どのような計算をしているかを確認してみた。

イベントまたは打ち切りの人数と発生した時点の情報は使用するものの、発生までの時間(生存時間)そのものを使ってはいない。

これで生存時間解析がどのような解析を行っているのか、理解が深まれば幸いである。

解説動画

youtu.be

参考サイト

ログランク検定 - Wikipedia

カイ二乗検定 - Wikipedia

おすすめ書籍

EZR公式マニュアル