大学院の授業でRを用いたサンプルサイズ(必要症例数)の算出方法のシミュレーションを学んだので、拡張して遊んでみた。
サンプルサイズの求め方はいくつかあるが、一例として
α(有意水準)、β(検出力)、SD(標準偏差)、興味のある効果量の差、を決めることで算出することができる。
通常はα=5%,β=80-90%として、SDは先行研究や手元の自施設データから決めていく事が多い。興味のある効果量の差は臨床的に重要と思われる最小の効果量を定めることで決定していく。
実際にはシミュレーションベースで必要症例数を決定していくことも多いらしい。
ということで症例数とSDを色々変えてみたときのβの変化の様子をシミュレーションで表現してみた。
＜シチュエーション＞
等分散の正規分布のデータを2つ用意。
2群の差を10で固定し、有意水準α=5%としたときの2群の平均値差が0か否か、つまり
H0: 平均値差=0
H1: 平均値差≠0
という仮説検定を考える。この場合はt検定が適用になる。
差を10で固定した2つの正規分布に従うデータをランダムに生成するので、当然2群間に差はあるのだが、同じ試行を何百回と繰り返していくと、確率的なゆらぎによって有意差を検出できない状況も生じてくる。
仮に100回検定をした場合に、80回は差を検出できたとき、検出力は80%といえる。
サンプルサイズとデータが従う正規分布の標準偏差（データのばらつき）を変えたときに、検出力βがどう変わるのか、をシミュレーションしていく。
下準備。
正規分布に従うデータを2個生成し、差が10となるようにする。
上記2群のデータから1000回t検定を繰り返すシミュレーションを実施し、有意水準5%で有意差を検出できた回数/全試行回数からβ(power)を算出するような自作関数を作成。
sim<-
function(x,y=10){ #x:症例数,y:SD, default=10
set.seed(1)
count<-0
n<-1000
for(i in 1:n){
a1<-rnorm(x,mean=120,sd=y)
a2<-rnorm(x,mean=130,sd=y)
p<-t.test(a1,a2,var.equal=T)$p.value
p01<-ifelse(p<0.05,1,0)
count<-count+p01
}
power<-100*(count/n)
return(power)
}
次に、症例数、SD(データのばらつき)を変えたときに検出力がどのように変化するのか可視化。今回は症例数を2~50までの範囲で設定。
#症例数と検出力、SDの関連性をグラフ化
res_sd<-function(y){ #y:SD
result<-
matrix(nrow=49,ncol=3,
dimnames=list(NULL,c("n","power(%)","SD"))
)
for(x in 2:50){
result[x-1,1]<-x
result[x-1,2]<-sim(x,y)
result[x-1,3]<-y
}
return(result)
}
#SDごとに結果を結合
results<-
res_sd(1) %>%
rbind(
res_sd(5),res_sd(10),res_sd(20),res_sd(50),res_sd(100)
)
#結果の図示
results %>%
as.data.frame() %>%
mutate(SD=as.factor(SD)) %>%
ggplot(aes(x=n,y=`power(%)`,color=SD))+
geom_point()+
geom_line()
横軸にサンプルサイズ(n)、縦軸に検出力(power)としてSDごとに色分けしてグラフを図示。
- サンプルサイズが増えるにつれ検出力も増大
- SDが小さいほうが検出力が上がりやすい
ということがわかる。
概ね予想通りの見た目であったが、想像よりもSDの影響が大きいんだなという印象。
SDが2群の差よりも大きい(SD: 20~)と、検出力を保つためには一気にサンプルサイズがたくさん必要になるっぽい。
簡単なシミュレーションだったが、R codingの良い練習になった。
