統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

ログランク検定のサンプルサイズ計算はどうやる?

kaiseki daiko banner

ログランク検定 log-rank testは、 時間経過とともにイベントが起きていくデータを検定する方法。

患者さんの死亡をイベントとしたデータ、
病気が再発することをイベントとしたデータ、
病気が発症することをイベントとしたデータ、
などが扱える。

病気ばかりではなく、 時間経過とともに起きてくる事柄なら、 なんでも扱える。

統計ソフトRでログランク検定のサンプルサイズ計算スクリプトを書いてみた。

ログランク検定はどうやるか?

サンプルサイズ計算の前に、 ログランク検定はどうやるか?を見てみる。

survivalパッケージのsurvdiff()を使う。

survivalパッケージはインストール済みなので、 library()で呼び出すだけでOK。

ovarianという卵巣がんの治療後生存時間の解析データ。

futimeが生存時間、fustatが死亡発生のデータ。rxが治療群のデータ。

library(survival)
survdiff(Surv(futime, fustat) ~ rx,data=ovarian)

結果はp=0.3で統計学的有意に異ならなかった。

> survdiff(Surv(futime, fustat) ~ rx,data=ovarian)
Call:
survdiff(formula = Surv(futime, fustat) ~ rx, data = ovarian)

      N Observed Expected (O-E)^2/E (O-E)^2/V
rx=1 13        7     5.23     0.596      1.06
rx=2 13        5     6.77     0.461      1.06

 Chisq= 1.1  on 1 degrees of freedom, p= 0.3 

もう一つの例は、MASSパッケージの中の、 gehanというデータ。

白血病の患者さんの寛解から再燃までのデータ。

timeが再燃までの時間、censが再燃したかどうか、 treatが6-メルカプトプリン(6-MP)で治療したか、対照だったか。

library(MASS)
survdiff(Surv(time, cens) ~ treat, data = gehan)

結果は、統計学的有意に6-MP治療群のほうが再燃が少なかった。

> library(MASS)
> survdiff(Surv(time, cens) ~ treat, data = gehan)
Call:
survdiff(formula = Surv(time, cens) ~ treat, data = gehan)

               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= 4e-05 

ログランク検定のサンプルサイズ計算はどうやるか?

ログランク検定のサンプルサイズ計算は、 実はCoxの比例ハザードモデルのサンプルサイズ計算を使うのだ。

toukeier.hatenablog.com

ログランク検定もCoxも生存時間解析と言いながら、 イベントが起きる順番しか使っていない。

あとはイベントが起きる比を使う。

サンプルサイズ計算はどちらにも流用できるということだ。

sample.size.cox <- function(S1,S0,alternative="two.sided",power=.8,
sig.level=.05){
alternative <- match.arg(alternative)
tside <- switch(alternative, one.sided=1, two.sided=2)

beta <- log(log(S1)/log(S0))
Za <- qnorm(sig.level/tside, lower.tail=FALSE)
Zb <- qnorm(power)
d <- ((Za+Zb)*(1+exp(beta))/(1-exp(beta)))^2
n <- d/(((1-S1)+(1-S0))/2)

NOTE <- "n is total number of subjects in *both* of groups"
METHOD <- "Sample Size Calculation of Cox Propotional Hazard Model"
structure(
list("Total number of death" = d,
"Total number of subjects" = n,
"Surv. rate of trt. arm" = S1,
"Surv. rate of ctl. arm" = S0,
"Hazard ratio" = exp(beta),
sig.level = sig.level,
power = power,
alternative = alternative,
note = NOTE,
method = METHOD),
class = "power.htest")
}

先ほどのovarianの結果をもとにサンプルサイズ計算をすると、 両群合わせて151例の死亡、326例の患者さんが必要だった。

各群13例の試験では統計学的有意にならなくて当然だ。

> sample.size.cox(S1=8/13, S0=6/13)

     Sample Size Calculation of Cox Propotional Hazard Model 

   Total number of death = 150.2536
Total number of subjects = 325.5495
  Surv. rate of trt. arm = 0.6153846
  Surv. rate of ctl. arm = 0.4615385
            Hazard ratio = 0.6279283
               sig.level = 0.05
                   power = 0.8
             alternative = two.sided

NOTE: n is total number of subjects in *both* of groups

gehanデータは、全部で8例の再燃症例があればよく、 必要なサンプルサイズは二群全体で11例と計算された。

> sample.size.cox(S1=(21-9)/21, S0=(21-21)/21)

     Sample Size Calculation of Cox Propotional Hazard Model 

   Total number of death = 7.84888
Total number of subjects = 10.98843
  Surv. rate of trt. arm = 0.5714286
  Surv. rate of ctl. arm = 0
            Hazard ratio = 0
               sig.level = 0.05
                   power = 0.8
             alternative = two.sided

NOTE: n is total number of subjects in *both* of groups

まとめ

統計ソフトRでログランク検定を行うには、 survivalパッケージのsurvdiff()を使う。

そのときイベント発生を表すデータと イベント発生までの時間、 群分けの変数が必要だ。

ログランク検定のサンプルサイズ計算は、 Coxの比例ハザードモデルのサンプルサイズ計算と同一の スクリプトを使って行える。

計算に必要なのは両群の生存確率だ。

エクセルでサンプルサイズ計算(2020年4月10日追記)

よければ以下からどうぞ。Coxの比例ハザードモデル用だが、上記の通りログランク検定のサンプルサイズ計算に使える。

https://happyhappygk.base.ec/items/26107287

参考書籍