統計ER

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

非劣性検定とは?どうやるの?

Toukei Consul Banner

KH Coder Consul Banner

いままでの方法と比べて、 格段にいいとか、 際立っていいとか、 じゃなくてもいい場合がある。

ダメじゃなければいい。 劣っていなければいい。

劣っていなければいいという検定方法が、 非劣性(ひれっせい)検定だ。

じゃ、どういうときが劣っていないって言うのか?

非劣性マージン:劣っているとは?劣っていないとは?

劣っていないとはどういう状態を言うのか?

劣っているとはどう状態なのか?

あなたは聞いて回ったことがあるだろうか?

臨床医は知っていると思うだろうか?

実は誰もわからないのだ。

非劣性検定で一番の問題はこの点だ。

「これ」以上劣っていたら、劣っているとする。 という閾値(しきいち)を決める。

「これ」が非劣性マージンと呼ばれる幅だ。

だが、非劣性マージンは何が正解か誰もわからない。

非劣性検定の非劣性マージンの決め方

参考になるスライドセットを見つけた。

www.slideshare.net

非劣性マージンを決めるときに必要なのは、 標準治療とプラセボの効果の差と、 保持率とか維持率と呼ばれる数値。

標準治療の効果と保持率=維持率の掛け算で、 非劣性マージンを決める。

保持率=維持率とは何か?

標準治療のプラセボに対する効果の 何割は効果があるか?ということ。

どんなに効果がなくても、 標準治療の効果の半分は効果がある。 もしくは効果があってほしい。

新しい治療はまだ効果がよくわからないから、 最悪でも標準治療の効果のどのくらいは効果がある ということはわからない。

だから、保持率=維持率は、あいまいだ。

主観。
というより、もはや願望。

どこまでいっても議論が尽きない理由がわかる。

非劣性試験の最大の弱点だ。

現実は、たくさん試験が行われていて、 何らかの理屈をつけて適切な試験とみなされている。

その結果、効果は変わらないが、 副作用が少ない治療法が生まれている。

どこまでも怪しさが残るが、 全否定するほどでもない。

これが 非劣性マージン、
非劣性検定、
非劣性試験だ。

母比率の非劣性検定はどうやる?

以下の条件で、下記のスクリプトで計算する。

使用する統計ソフトはRだ。

計算に使う関数を自作することができる。

今回は自作関数を使う。

自作関数はfunction()で作成する。

function()のあとに{}で関数の内容を指定する。

自作関数を使うときは、
関数名 <- function(){}
とする。

自作関数の全部を実行して、 そのあとに関数名()で使う。

臨床的に意味のある最小の差(DELTA)が 非劣性マージンだ。

どんなに最悪でもそこは下回らないという限界点。

  • pA.hat: 試験薬の有効率
  • pB.hat: 標準薬の有効率
  • pB.star: 帰無仮説(H0:pA=pB-DELTA)の下でのpB

  • DELTA: 臨床的に意味のある最小の差→10%(0.1)とする

  • rA: 試験薬群の有効症例数→64人と37人とする
  • nA: 試験薬群のサンプルサイズ→128とする
  • rB: 標準薬群の有効症例数→57人と39人とする
  • nB: 標準薬群のサンプルサイズ→127とする

STATISTICSが検定統計量で、 漸近的に正規近似できる 検定統計量だ。

non.inferiorが自作関数名だ。

non.inferior()で使う。

non.inferior <- function(nA, nB, rA, rB, DELTA, conf.level=0.9){
 METHOD <- "Non-Inferiority Test (Dunnett-Gent)"
 pB.star <- (rA+rB+nA*DELTA)/(nA+nB)
 se.delta.hat <- sqrt(((pB.star-DELTA)*(1-pB.star+DELTA))/nA+(pB.star*(1-pB.star))/nB)
 pA.hat <- rA/nA
 pB.hat <- rB/nB
 STATISTIC <- (pA.hat-(pB.hat-DELTA))/se.delta.hat
 names(STATISTIC) <- "Z"
 PVAL <- 1-pnorm(STATISTIC)
 se <- sqrt(pA.hat*(1-pA.hat)/nA+pB.hat*(1-pB.hat)/nB)
 delta <- pA.hat-pB.hat
 WIDTH <- qnorm(1-(1-conf.level)/2)
 CINT <- c(delta - WIDTH*se, delta + WIDTH*se)
 attr(CINT, "conf.level") <- conf.level
 RVAL <- list(statistic=STATISTIC, p.value=as.numeric(PVAL),
 estimate=c(pA=pA.hat, pB=pB.hat, pBstar=pB.star), conf.int=CINT, method=METHOD)
 class(RVAL) <- "htest"
 return(RVAL)
}

non.inferior(nA=128, nB=127, rA=64+37, rB=57+39, DELTA=0.1)

結果は以下の通り。

p値は片側で出している。

非劣性は最低この点は下回らないという仮説で、 片側しか興味がないため、 興味がある側の片側p値を出す。

90%信頼区間を出力しているのも、 下限5%のラインを見るためだ。

非劣性マージンー10%(-0.1)よりも、 信頼区間の下限が大きく-0.05のため、 非劣性と言えるわけだ。

> non.inferior(nA=128, nB=127, rA=64+37, rB=57+39, DELTA=0.1)

        Non-Inferiority Test (Dunnett-Gent)

data:  
Z = 2.5561, p-value = 0.005293
90 percent confidence interval:
 -0.05314985  0.11946383
sample estimates:
       pA        pB    pBstar 
0.7890625 0.7559055 0.8227451 

ただし、有効性100%の場合は計算できない。

> non.inferior(nA=188, nB=91, rA=188, rB=91, DELTA=0.1)

        Non-Inferiority Test (Dunnett-Gent)

data:  
Z = NaN, p-value = NA
90 percent confidence interval:
 0 0
sample estimates:
      pA       pB   pBstar 
1.000000 1.000000 1.067384 

Warning message:
In sqrt(((pB.star - DELTA) * (1 - pB.star + DELTA))/nA + (pB.star *  :
  NaNs produced

母比率の非劣性の検定を最尤推定量で行うには?

正規近似は100%を超えてしまうときは使えない。

現実には、両群とも有効率100%がありうる。

有効率100%を非劣性検定するには、 最尤推定量に基づく検定を利用する必要がある。

pB.star: 帰無仮説(H0:pA=pB-DELTA)の下でのpB

の計算が前節の正規近似と異なり、ややこしい。

non.inferior.likelihood <- function(nA, nB, rA, rB, DELTA, conf.level=0.9){
 METHOD <- "Non-Inferiority Test (Likelihood Method)"

 pA.hat <- rA/nA
 pB.hat <- rB/nB

 a <- nA+nB
 b <- -1*(nB+nA+rB+rA+DELTA*(nA+2*nB))
 c <- nB*DELTA^2+DELTA*(2*rB+nA+nB)+rB+rA
 d <- -rB*DELTA*(1+DELTA)
 v <- b^3/(27*a^3)-(b*c)/(6*a^2)+d/(2*a)
 u <- sign(v)*sqrt(b^2/(9*a^2)-c/(3*a))
 w <- (pi+acos(v/u^3))/3
 pB.star <- 2*u*cos(w)-b/(3*a)

 se.delta.hat <- sqrt(((pB.star-DELTA)*(1-pB.star+DELTA))/nA+(pB.star*(1-pB.star))/nB)

 STATISTIC <- (pA.hat-(pB.hat-DELTA))/se.delta.hat
 names(STATISTIC) <- "Z"
 PVAL <- 1-pnorm(STATISTIC)

 se <- sqrt(pA.hat*(1-pA.hat)/nA+pB.hat*(1-pB.hat)/nB)
 delta <- pA.hat-pB.hat
 WIDTH <- qnorm(1-(1-conf.level)/2)
 CINT <- c(delta - WIDTH*se, delta + WIDTH*se)
 attr(CINT, "conf.level") <- conf.level

 RVAL <- list(statistic=STATISTIC, p.value=as.numeric(PVAL),
 estimate=c(pA=pA.hat, pB=pB.hat, pBstar=pB.star), conf.int=CINT, method=METHOD)
 class(RVAL) <- "htest"
 return(RVAL)
}

non.inferior.likelihood(nA=128, nB=127, rA=64+37, rB=57+39, DELTA=0.1)

結果は以下の通り。

p値がほんの少し大きくなり、 帰無仮説の下でのpB(pBstar)は ほんの少し小さく見積もられている。

> non.inferior.likelihood(nA=128, nB=127, rA=64+37, rB=57+39, DELTA=0.1)

        Non-Inferiority Test (Likelihood Method)

data:  
Z = 2.5181, p-value = 0.0059
90 percent confidence interval:
 -0.05314985  0.11946383
sample estimates:
       pA        pB    pBstar 
0.7890625 0.7559055 0.8129256 

最尤推定量に基づく方法なら有効性100%でも検定できる。

> non.inferior.likelihood(nA=188, nB=91, rA=188, rB=91, DELTA=0.1)

        Non-Inferiority Test (Likelihood Method)

data:  
Z = 4.5704, p-value = 2.434e-06
90 percent confidence interval:
 0 0
sample estimates:
    pA     pB pBstar 
     1      1      1 

最尤推定量に基づく方法で検定しておけば、 いつでも問題ない。

母比率の非劣性検定 エクセル計算機

2021年10月9日追記

母比率の非劣性検定(最尤推定量に基づく方法)を計算するエクセル計算機を作ったので、よければどうぞ。

割合の非劣性検定 計算機【エクセル計算機】 | HHA SHOP

使い方解説動画も、よければどうぞ。

youtu.be

平均値の非劣性検定はどうやる?

変化量の平均値のような場合、 非劣性検定はどうやるか?

帰無仮説H_0と対立仮説H_1は以下の通りだ。

 \displaystyle H_0: \mu_A \leqq \mu_B - \Delta

 \displaystyle H_1: \mu_A > \mu_B - \Delta

臨床上意味のある最小の差を引いて、 同じかそれよりも劣る、が帰無仮説H_0

対立仮説 H_1は臨床上意味のある最小の差を引くと、 少なくともそのレベルは上回る効果がある。

non.inferiority.mean <- function(nA, xbarA, sdA, nB, xbarB, sdB, Delta, conf.level=.9){
 data.name <- sprintf("\nN_A = %s, Mean_A = %s, SD_A = %s\nN_B = %s, Mean_B = %s, SD_B = %s\nDelta = %s", nA, xbarA, sdA, nB, xbarB, sdB, Delta)

 METHOD <- "Non-inferiority test (Mean)"
 df <- nA-1+nB-1
 names(df) <- "df"
 SE <- sqrt((1/nA+1/nB)*(((nA-1)*sdA^2+(nB-1)*sdB^2)/df))
 STATISTIC <- (xbarA - (xbarB - Delta))/SE
 names(STATISTIC) <- "T"
 PVAL <- 1 - pt(STATISTIC, df)

 delta <- xbarA - xbarB
 WIDTH <- qt(1-(1-conf.level)/2, df)
 CINT <- c(delta-WIDTH*SE, delta+WIDTH*SE)
 attr(CINT, "conf.level") <- conf.level

 RVAL <- structure(list(statistics=STATISTIC, parameter=df,
 p.value=as.numeric(PVAL), data.name=data.name,
 conf.int=CINT, method=METHOD))
 class(RVAL) <- "htest"
 return(RVAL)
}

non.inferiority.mean(nA=46, xbarA=34.5, sdA=32.02, nB=44, xbarB=29.7, sdB=28.42, Delta=7)

結果は以下の通り。

平均の差は4.8で、 これはどんなに悪くても-7よりは小さくはならないはず というのが証明したい仮説だ。

計算結果、片側p値は0.034で、有意水準片側5%とすると統計学的有意。

90%信頼区間(片側5%)の下限は-5.8で、 -7より大きいため、 最悪でも-7よりは悪くない。

よって、非劣性が証明された。

> non.inferiority.mean(nA=46, xbarA=34.5, sdA=32.02, nB=44, xbarB=29.7, sdB=28.42, Delta=7)

        Non-inferiority test (Mean)

data:  
N_A = 46, Mean_A = 34.5, SD_A = 32.02
N_B = 44, Mean_B = 29.7, SD_B = 28.42
Delta = 7
T = 1.8459, df = 88, p-value = 0.03413
90 percent confidence interval:
 -5.826446 15.426446

平均値の非劣性検定 エクセル計算機

2021年10月10日追記

平均値の非劣性検定を計算するエクセル計算機を作ったので、よければどうぞ。

平均値の非劣性検定 計算機【エクセル計算機】 | HHA SHOP

使い方解説動画も、よければどうぞ。

youtu.be

まとめ

非劣性の検定を母比率の場合と母平均の場合に分けて紹介した。

母比率の場合は最尤推定量の基づく方法がベターだ。

非劣性検定は、非劣性マージンの設定が問題となる。

非劣性マージンが適切に設定され、 サンプルサイズが適切に見積もられてはじめて 非劣性検定が意味を成す。

参考書籍

丹後俊郎著 無作為化比較試験 朝倉書店
5. 臨床的同等性・非劣性の評価