トレンド検定のサンプルサイズ計算はどうやる?

にほんブログ村 科学ブログ 数学へ

トレンド検定のサンプルサイズ計算。 トレンドが見られればOKという試験ならこちら。

トレンドとはカテゴリが大きくなると、 平均値や割合が大きくなるとか小さくなるとか、 傾向や相関があることを指す。

統計ソフトRでやってみた。

平均値の場合

母数とスコアの相関係数を考える。

ポイントは、合計が0になるようにスコアを与えること。

c: 線形対比(ベクトル)←これがポイント。

K=3(3グループ)の場合はc(-1,0,1) K=4の場合はc(-3,-1,1,3):等間隔(間隔は2)

mu: 母平均(ベクトル)

sigma.sq: ANOVAから計算される誤差分散(下記分散分析表の黄色ハイライト部分)

Df Sum Sq Mean Sq F value Pr(>F)
tension 2 2034 1017.1 7.206 0.00175 **
Residuals 51 7199 141.1

スクリプトは以下の通り。

samplesize.mean.trend.test <- function(mu, sigma.sq, score, sig.level=.05, power=.8){
METHOD <- "Sample size determination for trend test of the means"
NOTE <- "n is each number in groups"
Za <- qnorm(sig.level/2, lower.tail=FALSE)
Zb <- qnorm(power)
n <- ((Za+Zb)^2 * sigma.sq * sum(score^2))/(sum(score*mu)^2)
structure(list(n=n, mu=mu, sigma.sq=sigma.sq, score=score, sig.level=sig.level, power=power, method=METHOD, note=NOTE), class = "power.htest")
}

 

参考書籍の例題は、以下のように計算される。

例題は分散分析表でなくSE(標準誤差)の値がわかるタイプ。

症状スコアの変化量は、-2, -3.5, -3.5と想定されるとする。

相関を見る線形対比は-2, 1, 1とする。

合計でゼロになるようにすることと、 症状スコアの変化量を表している数字にすること。

三群のうち、-2の群だけ違うわけなので、 その場合は-2, 1, 1というスコアが妥当。 (変化量と割り当てるスコアがともに-2なのはたまたまなので。念のため)

誤差分散は最大のSE=0.36の二乗に三群の症例数の真ん中あたり70例をかける。

SEはこんなふうに計算される。

$$ SE = \sqrt{\frac{分散}{n}} $$

よって二乗して症例数をかけると分散に戻る。

$$ \left(\sqrt{\frac{分散}{n}}\right)2n = \frac{分散}{n}n = 分散 $$

結果として一群48例必要と計算される。

> samplesize.mean.trend.test(mu=c(-2,-3.5,-3.5), sigma.sq=70*0.36^2, score=c(-2,1,1))

     Sample size determination for trend test of the means 

              n = 47.47002
             mu = -2.0, -3.5, -3.5
       sigma.sq = 9.072
          score = -2, 1, 1
      sig.level = 0.05
          power = 0.8

NOTE: n is each number in groups

 

もう一つの例として、 warpbreaksデータセットで計算してみると、 各群11例と計算される。

warpbreaksのデータでトレンド検定してみた結果はこちら(記事最下段)。

実際のデータは各群18例なので、 トレンド検定のためには少し多めだが悪くはない。

> samplesize.mean.trend.test(mu=c(36.4, 26.4, 21.7), sigma.sq=141.1, score=c(1,0,-1))

     Sample size determination for trend test of the means 

              n = 10.25015
             mu = 36.4, 26.4, 21.7
       sigma.sq = 141.1
          score = 1, 0, -1
      sig.level = 0.05
          power = 0.8

NOTE: n is each number in groups

 

割合の場合

割合のベクトルはprop、 線形対比ベクトルはscore、 でそれぞれインプットする。

samplesize.prop.trend.test <- function(prop, score, sig.level=.05, power=.8){
METHOD <- "Sample size determination for trend test of the proportions"
NOTE <- "n is each number in groups"
p <- mean(prop)
R <- sqrt(p*(1-p)*sum(score^2))
S <- sqrt(sum(score^2*prop*(1-prop)))
Za <- qnorm(sig.level/2, lower.tail=FALSE)
Zb <- qnorm(power)
n <- ((Za*R+Zb*S)/sum(score*prop))^2
structure(list(n=n, proportion=prop, score=score, sig.level=sig.level, power=power, method=METHOD, note=NOTE), class = "power.htest")
}

 

例えば有効率が55%、75%、75%の三群と予想できるとする。

線形対比スコアは平均値の時と同じ-2, 1, 1とする。 55%の群だけ有効性が低いと見えるからだ。

このときにトレンド検定で有意になるには各群65例必要と計算される。

> samplesize.prop.trend.test(prop=c(0.55,0.75,0.75), score=c(-2,1,1))

     Sample size determination for trend test of the proportions 

              n = 64.66423
     proportion = 0.55, 0.75, 0.75
          score = -2, 1, 1
      sig.level = 0.05
          power = 0.8

NOTE: n is each number in groups

もう一つ、carDataライブラリのTitanicSurvivalというデータで検証してみよう。

carDataライブラリは最初からインストールされているので、 呼び出すだけで使えるようになる。

タイタニック号が沈没したとき生存した乗客の性別、年齢、客室クラスのデータ。

客室クラスで生存割合にトレンドがあるかどうかを見てみる。

割合のトレンド検定はCochran-Armitage検定で、 prop.trend.test()でできる。

トレンド検定はどうやる?参照

library(carData)

with(TitanicSurvival, table(passengerClass, survived))

tab <- with(TitanicSurvival, table(passengerClass, survived))

prop.table(tab,1)

prop.table(tab,1)[,2]*100

barplot(prop.table(tab,1)[,2]*100, xlab="Passenger Class", ylab="Survived (%)")

prop.trend.test(tab[,2], rowSums(tab), score=c(1,0,-1))

 

生存割合は、一等が62%、二等が43%、三等が26%だった。

グラフにすると見事にトレンドがみられる。

結果は統計学的有意でトレンドあり。

> library(carData)
> 
> with(TitanicSurvival, table(passengerClass, survived))
              survived
passengerClass  no yes
           1st 123 200
           2nd 158 119
           3rd 528 181
> 
> tab <- with(TitanicSurvival, table(passengerClass, survived))
> 
> prop.table(tab,1)
              survived
passengerClass        no       yes
           1st 0.3808050 0.6191950
           2nd 0.5703971 0.4296029
           3rd 0.7447109 0.2552891
> 
> prop.table(tab,1)[,2]*100
     1st      2nd      3rd 
61.91950 42.96029 25.52891 
> 
> barplot(prop.table(tab,1)[,2]*100, xlab="Passenger Class", ylab="Survived (%)")
> 
> prop.trend.test(tab[,2], rowSums(tab), score=c(1,0,-1))

        Chi-squared Test for Trend in Proportions

data:  tab[, 2] out of rowSums(tab) ,
 using scores: 1 0 -1
X-squared = 127.81, df = 1, p-value < 2.2e-16

 

もし仮にサンプルサイズ計算をすると、 一群29例と計算される。

タイタニックは試験ではなく悲劇の事故なので、 サンプルサイズ計算の対象ではないが、 もしもの話なので、あしからず。

> samplesize.prop.trend.test(prop=c(0.62,0.43,0.26), score=c(1,0,-1))

     Sample size determination for trend test of the proportions 

              n = 28.603
     proportion = 0.62, 0.43, 0.26
          score = 1, 0, -1
      sig.level = 0.05
          power = 0.8

NOTE: n is each number in groups

 

まとめ

トレンド検定のサンプルサイズ計算を、 平均値の場合と割合の場合に分けて紹介した。

三群以上で、群間比較は不要で、 トレンドだけ見られればいいという状況は少ないと思うが、 もし必要な場合はこの記事のスクリプトを使えば計算できる。

参考書籍

出典:丹後俊郎著 無作為化比較試験 朝倉書店 3.7 傾向性検定-量反応関係の検出

新版