統計ER

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

Cox比例ハザードモデルのサンプルサイズ計算はどうやる?

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

コックス Cox の比例ハザードモデルで解析する場合のサンプルサイズ計算の方法。

いくつか方法があるのでご紹介。

Coxの比例ハザードモデルのサンプルサイズ計算その1

S0がコントロール、S1が治療群、それぞれの生存率。

dが、両群合わせて合計の死亡者数。

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

NOTE <- "d is total number of death in *both* of groups"
METHOD <- "Sample Size Calculation of Cox Propotional Hazard Model"
structure(
list("Total number of death" = d,
"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")
}

治療群の生存率が0.6、コントロールの生存率が0.5のとき、両群合わせて死亡者数が343例必要と計算される。

> sample.size.cox(S1=0.6, S0=0.5)

     Sample Size Calculation of Cox Propotional Hazard Model 

 Total number of death = 342.267
Surv. rate of trt. arm = 0.6
Surv. rate of ctl. arm = 0.5
          Hazard ratio = 0.7369656
             sig.level = 0.05
                 power = 0.8
           alternative = two.sided

NOTE: d is total number of death in *both* of groups

最近の抗がん剤は寿命延長が期待できるものも出てきたので、例えば5年でS1=0.8、S0=0.3というように、中央値に届かないアームも存在する。

0.8 vs 0.3となると、死亡症例数は17例で済み、ハザード比は0.19という極めて低い値が想定される。

> sample.size.cox(S1=0.8, S0=0.3)

     Sample Size Calculation of Cox Propotional Hazard Model 

 Total number of death = 16.6165
Surv. rate of trt. arm = 0.8
Surv. rate of ctl. arm = 0.3
          Hazard ratio = 0.1853394
             sig.level = 0.05
                 power = 0.8
           alternative = two.sided

NOTE: d is total number of death in *both* of groups

参考書籍

医学統計学シリーズ3 Cox比例ハザードモデル 中村 剛 著 朝倉書店
3.8 必要sample sizeの計算法

Coxの比例ハザードモデルのサンプルサイズ計算その2

統計ソフトRを使って自前でスクリプトに起こした関数を使う方法

Freedmanの方法とSchoenfeldの方法。

こちらは観察期間を設定してハザードを計算してくれるスクリプトだが、ハザード比は観察期間が異なっても同じで、サンプルサイズも同じである。

比例ハザードモデルの「比例」の由来通り、観察期間にかかわらず比例関係が一定で、観察期間が変わっても、必要なサンプルサイズは変わらない。

  • t: 観察期間(年)
  • S1 と S0: 治療群とコントロール群の生存率
  • dF: Freedmanの方法による各群必要な死亡症例数
  • dS: Schoenfeldの方法による各群必要な死亡症例数
sample.size.cox <- function(t,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))
H1 <- -1*log(S1)/t
H0 <- -1*log(S0)/t
HR <- H1/H0
Za <- qnorm(sig.level/tside, lower.tail=FALSE)
Zb <- qnorm(power)
dF <- (Za+Zb)^2*(HR+1)^2/(2*(HR-1)^2)
nF <- dF/(((1-S1)+(1-S0))/2)
dS <- 2*(Za+Zb)^2/((log(HR))^2)
nS <- dS/(((1-S1)+(1-S0))/2)

NOTE <- "n is number in *each* group"
METHOD <- "Sample Size Calculation of Cox Propotional Hazard Model"
structure(
list(
"Death (Freedman)" = dF,
"Number (Freedman)" = nF,
"Death (Schoenfeld)" = dS,
"Number(Schoenfeld)" = nS,
"Survival trt" = S1,
"Survival ctl" = S0,
"Hazard trt" = H1,
"Hazard ctl" = H0,
"Hazard ratio" = HR,
"Follow up(y)" = t,
sig.level = sig.level,
power = power,
alternative = alternative,
note = NOTE,
method = METHOD),
class = "power.htest")
}

5年の観察期間で、治療群が生存率0.8、コントロール群が生存率0.65と想定すると、必要死亡症例数は、Freedmanの方法で各39例、Schoenfeldの方法で各37例と計算される。

全体の必要症例数はFreedmanで各142例、Schoenfeldで各132例。

> sample.size.cox(t=5, S1=0.8, S0=0.65)

     Sample Size Calculation of Cox Propotional Hazard Model 

  Death (Freedman) = 38.92388
 Number (Freedman) = 141.5414
Death (Schoenfeld) = 36.27976
Number(Schoenfeld) = 131.9264
      Survival trt = 0.8
      Survival ctl = 0.65
        Hazard trt = 0.04462871
        Hazard ctl = 0.08615658
      Hazard ratio = 0.5179954
      Follow up(y) = 5
         sig.level = 0.05
             power = 0.8
       alternative = two.sided

NOTE: n is number in *each* group

治療群0.6、コントロール群0.5で計算すると、観察期間が1年でも5年でも同様で、

  • Freedmanでは死亡症例は各群172例、全体症例は各群381例必要、
  • Schoenfeldでは死亡症例は各群169例、全体症例は各群375例必要

と計算される。

>  sample.size.cox(t=1, S1=0.6, S0=0.5)

     Sample Size Calculation of Cox Propotional Hazard Model 

  Death (Freedman) = 171.1335
 Number (Freedman) = 380.2966
Death (Schoenfeld) = 168.5111
Number(Schoenfeld) = 374.4692
      Survival trt = 0.6
      Survival ctl = 0.5
        Hazard trt = 0.5108256
        Hazard ctl = 0.6931472
      Hazard ratio = 0.7369656
      Follow up(y) = 1
         sig.level = 0.05
             power = 0.8
       alternative = two.sided

NOTE: n is number in *each* group


> sample.size.cox(t=5, S1=0.6, S0=0.5)

     Sample Size Calculation of Cox Propotional Hazard Model 

  Death (Freedman) = 171.1335
 Number (Freedman) = 380.2966
Death (Schoenfeld) = 168.5111
Number(Schoenfeld) = 374.4692
      Survival trt = 0.6
      Survival ctl = 0.5
        Hazard trt = 0.1021651
        Hazard ctl = 0.1386294
      Hazard ratio = 0.7369656
      Follow up(y) = 5
         sig.level = 0.05
             power = 0.8
       alternative = two.sided

NOTE: n is number in *each* group

5年観察で、治療群が0.8、コントロール群が0.3の場合、
FreedmanとSchoenfeldの方法で、
それぞれ各群死亡症例は9例、6例、
全体症例は各群、19例、13例、
が必要と計算される。

> sample.size.cox(t=5, S1=0.8, S0=0.3)

     Sample Size Calculation of Cox Propotional Hazard Model 

  Death (Freedman) = 8.308251
 Number (Freedman) = 18.46278
Death (Schoenfeld) = 5.525171
Number(Schoenfeld) = 12.27816
      Survival trt = 0.8
      Survival ctl = 0.3
        Hazard trt = 0.04462871
        Hazard ctl = 0.2407946
      Hazard ratio = 0.1853394
      Follow up(y) = 5
         sig.level = 0.05
             power = 0.8
       alternative = two.sided

NOTE: n is number in *each* group

EZRを使った方法(2020年6月29日追記)

5年の観察期間で、治療群が生存率0.8、コントロール群が生存率0.65の条件で、EZRを使って計算してみると、Freedman式と同じ141例と計算された。

f:id:toukeier:20200629001805p:plain
EZRの設定(登録期間ゼロ、観察期間5年、5年生存率が0.8と0.65のとき)

観察研究のイメージで、登録は同時に行われることを想定して、登録期間はゼロとした。

> SampleHazard(0, 5, 5, 0.8, 0.65, 0.05, 0.80, 2, 1)
                               仮定
P1                              0.8
P2                             0.65
P1、P2の観察期間                  5
登録期間                          0
全研究期間                        5
αエラー                       0.05
                           両側検定
検出力                          0.8
N2とN1のサンプルサイズの比        1
                                   
必要サンプルサイズ         計算結果
N1                              141
N2                              141

治療群0.6、コントロール群0.5、観察期間5年で、EZRで同様に計算してみた。

f:id:toukeier:20200629002258p:plain
EZRの設定(観察期間5年、5年生存率を0.6と0.5とした)

上記のFreedman式の結果と同様に一群380例必要と計算された。

> SampleHazard(0, 5, 5, 0.6, 0.5, 0.05, 0.80, 2, 1)
                               仮定
P1                              0.6
P2                              0.5
P1、P2の観察期間                  5
登録期間                          0
全研究期間                        5
αエラー                       0.05
                           両側検定
検出力                          0.8
N2とN1のサンプルサイズの比        1
                                   
必要サンプルサイズ         計算結果
N1                              380
N2                              380

1年にしても同じ結果である。ただし、1年観察期間として、1年生存率と読み替えて計算することになる。

f:id:toukeier:20200629002602p:plain
EZRの設定(観察期間1年、1年生存率が0.6と0.5と考えた場合)

> SampleHazard(0, 1, 1, 0.6, 0.5, 0.05, 0.80, 2, 1)
                               仮定
P1                              0.6
P2                              0.5
P1、P2の観察期間                  1
登録期間                          0
全研究期間                        1
αエラー                       0.05
                           両側検定
検出力                          0.8
N2とN1のサンプルサイズの比        1
                                   
必要サンプルサイズ         計算結果
N1                              380
N2                              380

5年観察で、治療群が0.8、コントロール群が0.3の場合、EZRで計算すると、以下の通り一群18例と計算された。

f:id:toukeier:20200629003123p:plain
EZRの設定(5年生存率が0.8と0.3)

> SampleHazard(0, 5, 5, 0.8, 0.3, 0.05, 0.80, 2, 1)
                               仮定
P1                              0.8
P2                              0.3
P1、P2の観察期間                  5
登録期間                          0
全研究期間                        5
αエラー                       0.05
                           両側検定
検出力                          0.8
N2とN1のサンプルサイズの比        1
                                   
必要サンプルサイズ         計算結果
N1                               18
N2                               18

統計ソフトRのパッケージ powerSurvEpi を使った場合(2020年9月29日追記)

まず統計ソフトRに powerSurvEpi パッケージをインストールする。

install.packages("powerSurvEpi")

使うときは library() で呼び出す。

library(powerSurvEpi)

ssizeCT.default() という関数を使う。サンプルサイズ計算に必要な数値は以下の通り。

  • power : 検出力
  • k : 実験群(新薬群、実薬群など)と対照群(従来薬群、プラセボ群など)の比
  • pE : 実験群のイベント発生割合
  • pC : 対照群のイベント発生割合
  • RR : ハザード比
  • alpha : 有意水準(指定しなければ0.05)

治療群が生存率0.8、コントロール群が生存率0.65のとき、pE=1-0.8=0.2、pC=1-0.65=0.35である。またハザード比は、log(0.8)/log(0.65)で計算できる。この時一群142例と計算される。Freedman式による計算だ。

> ssizeCT.default(power=0.8, k=1, pE=1-0.8, pC=1-0.65, RR=log(0.8)/log(0.65))
 nE  nC 
142 142 

治療群0.6、コントロール群0.5の場合は、pE=1-0.6=0.4, pC=1-0.5=0.5, RR=log(0.6)/log(0.5)である。この時一群381例必要と計算される。

> ssizeCT.default(power=0.8, k=1, pE=0.4, pC=0.5, RR=log(0.6)/log(0.5))
 nE  nC 
381 381 

治療群が0.8、コントロール群が0.3の場合は、pE=1-0.8=0.2, pC=1-0.3=0.7, RR=log(0.8)/log(0.3)となり、一群19例必要と計算される。

> ssizeCT.default(power=0.8, k=1, pE=0.2, pC=0.7, RR=log(0.8)/log(0.3))
nE nC 
19 19 

Freedman式かSchoenfeld式か

Freedmanに比べて、Schoenfeldのほうが、小さいサンプルサイズでよいと計算される。

小さいサンプルサイズで済むならそれに越したことはない。倫理的に最小人数が適切である。

しかし、Schoenfeldの式は、必要なイベント数を過小評価していると指摘されている。

結論として、比較的多くの症例を集めることが可能であれば、Freedman式の計算結果をお勧めする。

参考文献

浜田&藤井 生存時間解析における症例数設計
第22回 日本SASユーザー会総会および研究発表会 論文集
2003年7月31日~8月1日

エクセルでサンプルサイズ計算

よければ以下からどうぞ。

コックス比例ハザードモデルのサンプルサイズ計算【エクセルでサンプルサイズ】 | HHA SHOP

f:id:toukeier:20200618205654p:plain
https://happyhappygk.base.ec/items/26107287

コックス比例ハザードモデルのサンプルサイズ計算【エクセルでサンプルサイズ】【動画】(2020年6月19日追記)

エクセルファイルの使い方動画を公開した。よかったらどうぞ。

youtu.be

EZRでCoxのサンプルサイズ計算【無料統計ソフトEZRで簡単統計】【動画】(2020年10月1日追記)

EZRでは、登録期間を考慮できる。すぐれもの!

youtu.be