ログランク検定のサンプルサイズ計算を R で行う方法
ログランク検定のサンプルサイズ計算の前に ログランク検定はどうやるか？
ログランク検定は、時間経過とともにイベントが起きていくデータの群間比較をする方法。
患者さんの死亡をイベントとしたデータ、病気が再発することをイベントとしたデータ、病気が発症することをイベントとしたデータ、などが扱える。
病気ばかりではなく、時間経過とともに起きてくる事柄なら、なんでも扱える。
サンプルサイズ計算の前に、ログランク検定はどうやるか？を見てみる。
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の比例ハザードモデルのサンプルサイズ計算を使う。
ログランク検定も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
ログランク検定のサンプルサイズ計算をエクセルで
Coxの比例ハザードモデル用だが、上記の通りログランク検定のサンプルサイズ計算に使える。
よければどうぞ。
まとめ
R でログランク検定を行うには、survivalパッケージのsurvdiff()を使う。
そのときイベント発生を表すデータとイベント発生までの時間、群分けの変数が必要だ。
ログランク検定のサンプルサイズ計算は、Coxの比例ハザードモデルのサンプルサイズ計算と同一のスクリプトを使って行える。
計算に必要なのは両群の生存確率だ。
