サンプルサイズブック
サンプルサイズブック
―コピペで使える！簡単便利！―
統計ソフト R で使えるサンプルサイズ計算のスクリプト集
2
3
まえがき
本書は、ブログ「一日いち R」に掲載した、サンプルサイズ計算のためのスクリプトを集め
たものである。スクリプト本体はブログに掲載したものと寸分たがわず同じものであるが、
記載の整理とスクリプトや結果の説明を加筆してある。また、ブログに掲載した順番より、
わかりやすい順番に並べ替えてある。
本書の利用方法はいたってシンプルである。下記のような太字のスクリプトを R のコンソ
ールにコピー＆ペーストする。
myPsize <- function(p,delta){
n <- 4*(1-p)/(p*delta^2)
c("N at least"=n)
}
次に、文中に出てくる結果セクションのプロンプト（>マーク）がついている先頭行を R の
エディタにコピー＆ペーストして、プロンプトは削除し、パラメータを希望に合わせて変
え、Run する。下記の場合 0.08 や 0.2 を変更する。
> myPsize(0.08,0.2)
参考文献中の式をスクリプトとして忠実に表現し、例題の解答とあっていることを確認し
ているが、誤りを見つけた場合はご連絡いただければ幸いである。ブログ「一日いち R」の
コメント欄に記載していただければと思う。
サンプルサイズを計算したいのだが、手ごろなソフトウェアや R のパッケージが見つから
ないと悩んでいる諸氏のお役にたてれば幸いである。
なお、本書に掲載していない検定・試験デザインのサンプルサイズ計算は、随時追加する
予定であるので期待していただきたい。
2014 年 6 月 30 日
一日いち R 管理人
4
目次
母比率の推定···························································5
平均値の差の検定 ·················································· 6
ANOVA―誤差分散がわかる場合―··········································7
ANOVA―誤差分散がわからない場合― ·····································9
Dunnett 検定·························································12
傾向性の検定··························································15
相関係数の検定························································17
ロジスティック回帰モデル···············································19
Cox 比例ハザードモデル·················································21
割合の非劣性検定（Dunnett-Gent の方法）
···································24
割合の非劣性検定（最尤推定量に基づく方法）
·································26
平均値の非劣性検定····················································29
Cox 比例ハザードモデルの非劣性検定······································31
クロスオーバー試験····················································35
Matched case-control study··············································37
5
母比率の推定
相対精度 delta で，母比率を推定する際のサンプルサイズ計算。One arm の試験に応用可
能。
myPsize <- function(p,delta){
n <- 4*(1-p)/(p*delta^2)
c("N at least"=n)
}
8%を相対精度 20%（これを 95%信頼区間とする）で推定するためには、1150 例必要。1150
例で 8%（92 例）だった場合、95%信頼区間は 6.5-9.7%。±1.6%(8%の 20%)程度であるこ
とを binom.test()で確認している。
> myPsize(0.08,0.2)
N at least
1150
> binom.test(1150*0.08, 1150)
Exact binomial test
data: 1150 * 0.08 and 1150
number of successes = 92, number of trials = 1150, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.06497587 0.09721197
sample estimates:
probability of success
0.08
参考文献
丹後俊郎著 古川俊之監修
新版 医学への統計学
9.1.2 母比率の推定 pp.184-185
朝倉書店
6
平均値の差の検定
２つの母平均の差の検定（いわゆる t 検定）を行って解析する場合のサンプルサイズ計算。
power.t.test(power=, delta=)
有意水準両側 5%（デフォルト）
、検出力 80%、delta（２つの母平均の差を母標準偏差で割
ったもの）が 0.7 とした場合のサンプルサイズを見積もった場合、一群 34 人必要と計算さ
れる。
> power.t.test(power=.8, delta=.7)
Two-sample t test power calculation
n = 33.02467
delta = 0.7
sd = 1
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group
delta については、参考文献の p.186 にある「慣例的性質」についての記述が役に立つ。 以
下 d は delta と同じ。
「d の推定は困難である場合が多いので、慣例的性質を使い目安をつ
けるのも一つの方法」とのこと。
1) 小さな差を検出したければ d=0.1～0.2
2) 中位な差を検出したければ d=0.4～0.5
3) 大きな差を検出したければ d=0.8～0.9
参考文献
丹後俊郎著 古川俊之監修
新版
医学への統計学
朝倉書店
9. 標本の大きさの求め方 9.2.1 二つの母平均の差の検定 pp.185-192
7
ANOVA―誤差分散がわかる場合―
連続変数を 3 群以上の群間で比較する分散分析のためのサンプルサイズ計算。検出力 80%・
90%、有意水準 5%・1%それぞれの組み合わせ（4 通り）に対応している。
k：群の数
d：群間で最大の差（母平均の差の最大値）
sigma0：誤差の母分散の平方根
anova.sample.size <- function(k, d, sigma0, alpha=0.05, power=.8){
DELTA <- d^2/(2*sigma0^2)
phiA <- k-1
if (power==.8 & alpha==.05) {lambda <- 4.860+3.584*sqrt(phiA)}
else if (power==.9 & alpha==.05) {lambda <- 7.049+4.244*sqrt(phiA)}
else if (power==.8 & alpha==.01) {lambda <- 7.736+4.551*sqrt(phiA)}
else if (power==.9 & alpha==.01) {lambda <- 10.439+5.213*sqrt(phiA)}
n <- lambda/DELTA
METHOD <- "ANOVA Sample Size Calculation"
structure(list(n=n, delta=d, sd=sigma0, sig.level=alpha, power=power,
method=METHOD),class="power.htest")
}
以下は、4 群比較で、最大の 2 群間の差を 2、誤差分散の平方根を 1、有意水準を 5％とし
て、検出力を 80％または 90％とすると、以下のように、それぞれ一群 6 例、8 例が必要と
計算される。
> anova.sample.size(k=4, d=2, sigma0=1)
ANOVA Sample Size Calculation
n = 5.533835
delta = 2
sd = 1
sig.level = 0.05
power = 0.8
8
> anova.sample.size(k=4, d=2, sigma0=1, power=.9)
ANOVA Sample Size Calculation
n = 7.199912
delta = 2
sd = 1
sig.level = 0.05
power = 0.9
参考文献
永田靖著 サンプルサイズの決め方
朝倉書店 pp.129-164
9.3 (誤差分散が既知の場合の)サンプルサイズ設計方法 pp.136-139
9
ANOVA―誤差分散がわからない場合―
誤差分散がわからない場合の分散分析のサンプルサイズ計算。わからないなりにも、わか
っているとして、まず既知の場合でサンプルサイズ見積もりをする。そのｎで未知の場合
の検出力を計算し、未知の場合のｎを求める。ｎを１ずつずらして検出力を計算して、微
調整をするような感じ。
k：群の数
n：サンプルサイズ
d：群間で最大の差（母平均の差の最大値）
sigma0：誤差母分散の推定値の平方根
anova.power.michi <- function(k, n, d, sigma0, alpha=.05){
DELTA <- d^2/(2*sigma0^2)
phiA <- k-1
phiE <- k*(n-1)
lambda <- n*DELTA
cA <- (phiA+2*lambda)/(phiA+lambda)
phiA.star <- (phiA+lambda)^2/(phiA+2*lambda)
w <- qf(p=alpha, df1=phiA, df2=phiE, lower.tail=FALSE)
u
<-
(sqrt(w/phiE)*sqrt(2*phiE-1)-sqrt(cA/phiA)*sqrt(2*phiA.star-1))/(sqrt(c
A/phiA+w/phiE))
power <- 1-pnorm(u)
METHOD <- "ANOVA Power Calculation Variance unknown"
structure(list(n=n, delta=d, sd=sigma0, sig.level=alpha, power=power,
method=METHOD), class="power.htest")
}
以下は、例として、3 群で、一群 5 例、最大の群間差を 2、誤差分散を 2 としたとき、検出
力が約 41％であると計算されている。
> anova.power.michi(k=3, n=5, d=2, sigma0=sqrt(2))
ANOVA Power Calculation Variance unknown
10
n = 5
delta = 2
sd = 1.414214
sig.level = 0.05
power = 0.4088796
誤差分散がわからない場合は、以下のように検討していく。まず誤算分散がわからないな
りにも仮に 1 だとして、4 群で、最大の群間差を 2 として、有意水準を 5％、検出力を 90％
とした場合のサンプルサイズを計算すると、以下の通り一群 8 例と計算される。
> anova.sample.size(k=4, d=2, sigma0=1, power=.9)
ANOVA Sample Size Calculation
n = 7.199912
delta = 2
sd = 1
sig.level = 0.05
power = 0.9
次に、一群 8 例の場合と 9 例の場合の検出力を、誤差分散がわからない時の計算方法で検
討する。
そうすると 9 例の時のほうが 90％を上回り、検出力 90％で検定したいのであれば、
一群 9 例のほうが適切であることがわかる。
> anova.power.michi(k=4, n=c(8,9), d=2, sigma0=1)
ANOVA Power Calculation Variance unknown
n = 8, 9
delta = 2
sd = 1
sig.level = 0.05
power = 0.8955750, 0.9338735
また、ほかの例として、5 群比較で、最大の群間差を 3、誤差分散を 3、有意水準 1％、検
出力 80％としたときのサンプルサイズ計算を行うと、一群 12 例と計算される。
11
> anova.sample.size(k=5, d=3, sigma0=sqrt(3), alpha=.01)
ANOVA Sample Size Calculation
n = 11.22533
delta = 3
sd = 1.732051
sig.level = 0.01
power = 0.8
同様に、誤差分散がわからないとして、一群 12 例と 13 例の時の検出力を比べてみると、
13 例の時が 80％を超えることがわかる。つまり検出力 80％で検定したいのであれば、一
群 13 例が適切であると言える。
> anova.power.michi(k=5, n=c(12,13), d=3, sigma0=sqrt(3), alpha=.01)
ANOVA Power Calculation Variance unknown
n = 12, 13
delta = 3
sd = 1.732051
sig.level = 0.01
power = 0.7754018, 0.8240690
参考文献
永田靖著 サンプルサイズの決め方
朝倉書店
10 １元配置分散分析－誤差分散が未知の場合 pp.129-164
10.2 検出力の計算方法 pp.150-155
10.3 サンプルサイズの設計方法 pp.156-158
12
Dunnett 検定
Dunnett（ダネット）検定は、対照群と複数の試験群の平均値を比較するときに、多重性を
考慮した検定方法。このセクションの参考文献である Liu の論文中の数値表を用いて計算
する。データフレーム dat.tab に数値表を格納して計算している。
dat.tab
<-
read.table("http://homepage3.nifty.com/ttyotani/ichirdat/110921_liu_199
7_table2_nu40_rev.txt")
sample.size.dunnett <function(
nu=40,sig.level=c(.1,.05,.01),power=c(.75,.80,.85,.90,.95),
treat=c(2:10),Delta,sd=1
)
{
dat.tab.abs <- subset(dat.tab, alpha==sig.level & Power==power & k==treat)
n <- (dat.tab.abs$sqrt.n.delta/Delta*sd)^2
n0 <- sqrt(treat)*n
NOTE <- "n is each number in *treatment* groups"
METHOD <- "Sample size determination for Dunnett test"
structure(
list("n0 (control)" = n0, "n (treatment)" = n,
"Treatment groups" = treat,
"Diff. to be rejected" = Delta,
SD = sd,
sig.level = sig.level,
power = power, alternative = "two.sided", note = NOTE,
method = METHOD),
class = "power.htest")
}
有意水準両側 5％、
検出力 80％、
試験群 3 群、
帰無仮説の棄却が期待される平均値の差が 1、
標準偏差が 1（デフォルト）とすると、対照群が 43 人、試験群各群 25 人（試験群合計 75
13
人）必要と計算される。
> sample.size.dunnett(sig.level=.05,power=.8,treat=3,Delta=1)
Sample size determination for Dunnett test
n0 (control) = 42.71438
n (treatment) = 24.66116
Treatment groups = 3
Diff. to be rejected = 1
SD = 1
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is each number in *treatment* groups
平均値の差を 0.5 と小さくすると、対照群が 171 人、試験各群 99 人必要と、必要例数が増
加する。
> sample.size.dunnett(sig.level=.05,power=.8,treat=3,Delta=0.5)
Sample size determination for Dunnett test
n0 (control) = 170.8575
n (treatment) = 98.64462
Treatment groups = 3
Diff. to be rejected = 0.5
SD = 1
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is each number in *treatment* groups
平均値の差を 1 に戻し、標準偏差を 2 倍の 2 にしても上記と同じ結果が得られる。
14
> sample.size.dunnett(sig.level=.05,power=.8,treat=3,Delta=1,sd=2)
Sample size determination for Dunnett test
n0 (control) = 170.8575
n (treatment) = 98.64462
Treatment groups = 3
Diff. to be rejected = 1
SD = 2
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is each number in *treatment* groups
参考文献
Liu W. On sample size determination of Dunnett's procedure for comparing several
treatments with a control. J Stat Plan Inference 1997;62:255-261.
永田 靖
ダネットの方法におけるサンプルサイズの決め方について 医薬安全性研究会
会報 2005;50(May):51-59.
15
傾向性の検定
3 群以上の比較における、平均値の傾向性検定のサンプルサイズ計算。母数とスコアの相関
係数を考える。合計が 0 になるようにスコアを与えるのがポイント。例えば、K=3(3 グル
ープ)の場合は c(-1,0,1)、K=4 の場合は c(-3,-1,1,3):等間隔（間隔は 2）など。ANOVA から
計算される誤差分散が必要であるが、実際の場面では、これを見つけるのが難しいように
思う。
score：線形対比（ベクトル）
mu：母平均（ベクトル）
sigma.sq：ANOVA から計算される誤差分散
samplesize.mean.trend.test <- function(mu, sigma.sq, score, sig.level=.05,
power=.8){
Za <- qnorm(sig.level/2, lower.tail=FALSE)
Zb <- qnorm(power)
n <- ((Za+Zb)^2 * sigma.sq * sum(score^2))/(sum(score*mu)^2)
list(n=n, mu=mu, sigma.sq=sigma.sq, score=score, sig.level=sig.level,
power=power)
}
3 群比較として、母平均が－2、－3.5、－3.5 であるとして、誤差分散が 70×0.362、スコ
アを－2、1、1、有意水準 5％、検出力 80％とすると、以下の通り、一群 48 人必要と計算
される。
>
samplesize.mean.trend.test(mu=c(-2,-3.5,-3.5),
score=c(-2,1,1))
$n
[1] 47.47002
$mu
[1] -2.0 -3.5 -3.5
$sigma.sq
[1] 9.072
sigma.sq=70*0.36^2,
16
$score
[1] -2 1 1
$sig.level
[1] 0.05
$power
[1] 0.8
参考文献
丹後俊郎著 無作為化比較試験
朝倉書店
3.7 傾向性検定-量反応関係の検出
17
相関係数の検定
相関係数を求める際のサンプルサイズ計算。検定の帰無仮説は、母相関係数が 0（ゼロ）で
あるため、統計学的有意でも母相関係数が 0 ではないと検証されるに過ぎない。つまり、
相関係数の検定が統計学的有意であることと、そのサンプルの相関係数が実質的に意味の
あるものかどうかは、必ずしも一致しない。相関係数が真の意味があるかどうかの判定は
以下の表を参考に。
表 相関係数の解釈
相関係数の絶対値
解釈
0.0～0.2
ほとんど相関関係がない
0.2～0.4
やや相関関係がある
0.4～0.7
かなり相関関係がある
0.7～1.0
強い相関関係がある
引用：ピアソンの積率相関係数
http://aoki2.si.gunma-u.ac.jp/lecture/Soukan/pearson.html
samplesize.cor.test <- function(r, sig.level=.05, power=.8){
Za <- qnorm(sig.level/2, lower.tail=FALSE)
Zb <- qnorm(power)
C <- 0.5*log((1+r)/(1-r))
N <- ((Za+Zb)/C)^2 + 3
c(N=N, r=r, alpha=sig.level, Power=power)
}
相関係数 0.2、0.3、0.4、0.5、0.6、0.7 を、有意水準 5％、検出力 80％で、統計学的有意
に検出するためには、それぞれ、194 例、85 例、47 例、30 例、20 例、14 例が必要と算出
される。
> samplesize.cor.test(0.2)
N
r
alpha
Power
193.9680
0.2000
0.0500
0.8000
> samplesize.cor.test(0.3)
N
r
alpha
Power
84.92781 0.30000 0.05000
0.80000
18
> samplesize.cor.test(0.4)
N
r
alpha
Power
46.73161 0.40000 0.05000
0.80000
> samplesize.cor.test(0.5)
N
r
29.0123
0.5000
alpha
Power
0.0500 0.8000
> samplesize.cor.test(0.6)
N
r
alpha
Power
19.33641 0.60000 0.05000
0.80000
> samplesize.cor.test(0.7)
N
r
alpha
Power
13.43442 0.70000 0.05000
0.80000
また、相関係数 0.4 を有意水準 1％、検出力 80％で検出する場合、有意水準 5％、検出力
90％で検出する場合は、それぞれ、69 例、62 例が必要と計算される。
> samplesize.cor.test(0.4, sig.level=.01)
N
r
alpha
Power
68.07171 0.40000 0.01000
0.80000
> samplesize.cor.test(0.4, power=.9)
N
r
alpha
Power
61.54422 0.40000 0.05000
0.90000
参考文献
Stephen B. Hulley & Steven R. Cummings, 木原正博監訳、現代疫学研究会訳
医学的研究のデザイン
研究の質を高める疫学的アプローチ
メディカル・サイエンス・インターナショナル
Appendix 13-C 相関係数を分析に用いる場合に必要なサンプルサイズ
19
ロジスティック回帰モデル
ロジスティック回帰モデルにおいて、連続変数を用いたサンプルサイズ計算。1SD 上昇し
た時の発生確率を想定して計算する。
p0：平均レベルにおける発生確率
p1：平均＋1SD レベルにおける発生確率
sample.size.logistic <- function(p0, p1,
sig.level=.05, power=.8,alternative=c("two.sided","one.sided")){
alternative <- match.arg(alternative)
tside <- switch(alternative, one.sided=1, two.sided=2)
odds0 <- p0/(1-p0)
odds1 <- p1/(1-p1)
theta <- odds1/odds0
lambda <- log(theta)
Za <- qnorm(sig.level/tside, lower.tail=FALSE)
Zb <- qnorm(power)
delta <- (1+(1+lambda^2)*(exp(5*lambda^2/4)))/(1+exp(-lambda^2/4))
n <- ((Za+Zb*exp(-lambda^2/4))^2*(1+2*p0*delta))/(p0*lambda^2)
NOTE <- "n is size of entire cohort"
METHOD <- "Sample size, logistic reg. with cont. var."
structure(list(n = n,
"Prob. in average level" = p0,
"Prob. in ave.+SD level" = p1,
sig.level = sig.level,
power = power,
alternative = alternative, note = NOTE,
method = METHOD), class = "power.htest")
}
20
平均レベルの発生確率を 8％、平均＋1SD レベルの発生確率を 8％の 1.5 倍＝12％、検出力
90％、片側検定とすると、コホート全体の人数は 613 人必要と計算される。
> sample.size.logistic(p0=0.08, p1=0.08*1.5, power=.9, alternative="one.sided")
Sample size, logistic reg. with cont. var.
n = 612.0427
Prob. in average level = 0.08
Prob. in ave.+SD level = 0.12
sig.level = 0.05
power = 0.9
alternative = one.sided
NOTE: n is size of entire cohort
検出力 80％、両側検定にすると、コホート全体で 569 人必要と計算される。
> sample.size.logistic(p0=0.08, p1=0.08*1.5)
Sample size, logistic reg. with cont. var.
n = 568.7537
Prob. in average level = 0.08
Prob. in ave.+SD level = 0.12
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is size of entire cohort
参考文献
Alan Agresti 著
渡邉裕之・菅波秀規・吉田光宏・角野修司・寒水孝司・松永信人 共訳
カテゴリカルデータ解析入門 サイエンティスト社
Hsieh FY. Sample size tables for logistic regression. Stat Med 1989;8(7):795-802.
21
Cox 比例ハザードモデル
Cox の比例ハザードモデルにおけるサンプルサイズ計算。Freedman の方法と Shoenfeld
の方法がある。
t：観察期間（年）
S1：試験群の生存率
S0：対照群の生存率
dF：Freedman の方法で算出した各群の死亡数
dS：Shoenfeld の方法で算出した各群の死亡数
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 (Shoenfeld)" = dS,
"Number(Shoenfeld)" = nS,
"Survival trt" = S1,
22
"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 として、有意水準 5％、
検出力 80％とすると、Freedman の方法では各群 142 人、Shoenfeld の方法では各群 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 (Shoenfeld) = 36.27976
Number(Shoenfeld) = 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
23
参考文献
浜田＆藤井 生存時間解析における症例数設計 第 22 回 日本 SAS ユーザー会総会およ
び研究発表会 論文集 2003 年 7 月 31 日～8 月 1 日
東京
24
割合の非劣性検定（Dunnett-Gent の方法）
割合がエンドポイントの場合の非劣性検定のサンプルサイズ計算(Dunnett-Gent の方法)。
以下の設定のほか、有意水準を片側 5%、検出率 80%をデフォルト設定としている。
pA：試験薬の有効率
pB：標準薬の有効率
pB.bar：帰無仮説 H0: pA=pB-DELTA の下で推定された pB.star の漸近値
delta：pA-pB
DELTA：臨床的に意味のある最小の差（「有効率が下回る差はこれよりも小さい」が証
明したい仮説）
non.inferior.sample.size <- function(pA, pB, DELTA, power=.8, sig.level=.05,
alternative="one.sided"){
alternative <- match.arg(alternative)
tside <- switch(alternative, one.sided = 1, two.sided = 2)
delta <- pA-pB
pB.bar <- pB+(delta+DELTA)/2
R <- sqrt((pB.bar-DELTA)*(1-pB.bar+DELTA)+pB.bar*(1-pB.bar))
S <- sqrt(pA*(1-pA)+pB*(1-pB))
Za <- qnorm(sig.level/tside, lower.tail=FALSE)
Zb <- qnorm(power)
n <- ((Za*R+Zb*S)/(delta+DELTA))^2
NOTE <- "n is number in *each* group"
METHOD <- "Non Inferiority Test Sample Size Calculation (Dunnett-Gent)"
structure(list(n = n, pA = pA, pB = pB, DELTA=DELTA, sig.level = sig.level,
power = power, alternative = alternative, note = NOTE,
method = METHOD), class = "power.htest")
}
例として、試験薬が有効である被験者の割合を 0.813（81.3％）
、標準薬が有効である被験
25
者の割合を 0.741（74.1％）と見積もり、臨床的に意味のある最小の差を 0.1（10％）とし
て、検出力 90％とすると、一群 100 例必要と計算される。
> non.inferior.sample.size(pA=0.813, pB=0.741, DELTA=.1, power=.9)
Non Inferiority Test Sample Size Calculation (Dunnett-Gent)
n = 99.17305
pA = 0.813
pB = 0.741
DELTA = 0.1
sig.level = 0.05
power = 0.9
alternative = one.sided
NOTE: n is number in *each* group
参考文献
丹後俊郎著 無作為化比較試験
朝倉書店
5. 臨床的同等性・非劣性の評価
5.3 独立な母比率の差の検定 pp.100-105
26
割合の非劣性検定（最尤推定量に基づく方法）
割合がエンドポイントの場合の非劣性検定のサンプルサイズ計算で、最尤推定量に基づく
方法。以下の設定のほか、有意水準片側 5%、検出率 80%をデフォルト設定としている。
pA：試験薬の有効率
pB：標準薬の有効率
pB.bar：帰無仮説 H0: pA=pB-DELTA の下で推定された pB.star の漸近値
delta：pA-pB
DELTA：臨床的に意味のある最小の差
non.inferior.sample.size.likelihood <- function(pA, pB, DELTA, power=.8,
sig.level=.05, alternative="one.sided"){
alternative <- match.arg(alternative)
tside <- switch(alternative, one.sided = 1, two.sided = 2)
delta <- pA-pB
a <- 2
b <- -2*pB-2-3*DELTA-delta
c <- DELTA^2+2*(1+pB)*DELTA+2*pB+delta
d <- -pB*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)
R <- sqrt((pB.star-DELTA)*(1-pB.star+DELTA)+pB.star*(1-pB.star))
S <- sqrt(pA*(1-pA)+pB*(1-pB))
Za <- qnorm(sig.level/tside, lower.tail=FALSE)
Zb <- qnorm(power)
n <- ((Za*R+Zb*S)/(delta+DELTA))^2
NOTE <- "n is number in *each* group"
27
METHOD <- "Non Inferiority Test Sample Size Calculation (Likelihood
Method)"
structure(list(n = n, pA = pA, pB = pB, DELTA=DELTA, sig.level = sig.level,
power = power, alternative = alternative, note = NOTE,
method = METHOD), class = "power.htest")
}
Dunnett-Gent の方法と同様に、試験薬が有効である被験者の割合を 0.813（81.3％）
、標準
薬が有効である被験者の割合を 0.741（74.1％）と見積もり、臨床的に意味のある最小の差
を 0.1（10％）として、検出力 90％とすると、一群 102 例必要と計算される。
> non.inferior.sample.size.likelihood(pA=0.813, pB=0.741, DELTA=.1, power=.9)
Non Inferiority Test Sample Size Calculation (Likelihood Method)
n = 101.4188
pA = 0.813
pB = 0.741
DELTA = 0.1
sig.level = 0.05
power = 0.9
alternative = one.sided
NOTE: n is number in *each* group
例として、2/72=2.8%と 2/132=1.5%を比べる場合、1.5％が 2.8％に比べて 0.02（2%）は小
さくないという仮説を検証するには一群約 5000 例必要となる。このような試験設定はあま
り現実的でないと言える。
> non.inferior.sample.size.likelihood(pA=2/132, pB=2/72, DELTA=.02)
Non Inferiority Test Sample Size Calculation (Likelihood Method)
n = 4982.974
pA = 0.01515152
pB = 0.02777778
28
DELTA = 0.02
sig.level = 0.05
power = 0.8
alternative = one.sided
NOTE: n is number in *each* group
参考文献
丹後俊郎著 無作為化比較試験
朝倉書店
5. 臨床的同等性・非劣性の評価
5.3 独立な母比率の差の検定 pp.100-105
29
平均値の非劣性検定
平均値を用いた非劣性検定のサンプルサイズ計算。
Delta：非劣性マージン
Delta1：対立仮説 H1 における真の差
sample.size.non.inferiority.mean <- function(sig.level=.05, power=.8,
Delta, Delta1, sd, alternative="one.sided"){
alternative <- match.arg(alternative)
tside <- switch(alternative, one.sided=1, two.sided=2)
d <- (Delta + Delta1)/sd
Za <- qnorm(sig.level/tside, lower.tail=FALSE)
Zb <- qnorm(power)
n <- 2*((Za+Zb)/d)^2
NOTE <- "n is number in *each* group"
METHOD <- "Non Inferiority Test Sample Size Calculation (Mean)"
structure(list(n = n, "Non-inferiority margin" = Delta,
"True difference"=Delta1, SD = sd, sig.level = sig.level,
power = power, alternative = alternative, note = NOTE,
method = METHOD), class = "power.htest")
}
非劣性マージンを 7、真の差を 34.5-29.7＝4.8、SD を 30 として、有意水準片側 5％、検出
力 80％では、一群 80 人必要と計算される。
> sample.size.non.inferiority.mean(Delta=7, Delta1=34.5-29.7, sd=30)
Non Inferiority Test Sample Size Calculation (Mean)
n = 79.92389
Non-inferiority margin = 7
True difference = 4.8
SD = 30
30
sig.level = 0.05
power = 0.8
alternative = one.sided
NOTE: n is number in *each* group
参考文献
丹後俊郎著 無作為化比較試験
朝倉書店
5. 臨床的同等性・非劣性の評価
5.5 独立な母平均の差の検定 pp.113-117
31
Cox 比例ハザードモデルの非劣性検定
Cox の比例ハザードモデルを用いた非劣性検定のサンプルサイズ計算。Freedman と
Shoenfeld の方法がある。
t：観察期間（年）
S1：試験群の生存率
S0：対照群の生存率
dF：Freedman の方法で算出した各群の死亡数
dS：Shoenfeld の方法で算出した各群の死亡数
Delta.S0：非劣性マージン
sample.size.cox.noninf <- function(t,S1,S0,Delta.S0,
alternative=c("one.sided","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
H0.prime <- -1*log(S0-Delta.S0)/t
HR.prime <- H1/H0.prime
Za <- qnorm(sig.level/tside, lower.tail=FALSE)
Zb <- qnorm(power)
dF <- (Za+Zb)^2*(HR.prime+1)^2/(2*(HR.prime-1)^2)
nF <- dF/(((1-S1)+(1-S0))/2)
dS <- 2*(Za+Zb)^2/((log(HR.prime))^2)
nS <- dS/(((1-S1)+(1-S0))/2)
NOTE <- "n is number in *each* group"
METHOD <- "Sample Size of Cox Model (noninferiority)"
structure(
list(
"Death (Freedman)" = dF,
"Number (Freedman)" = nF,
32
"Death (Shoenfeld)" = dS,
"Number(Shoenfeld)" = nS,
"Survival trt(S1)" = S1,
"Survival ctl(S0)" = S0,
"Noninf. margin"= Delta.S0,
"Hazard trt(H1)" = H1,
"Hazard ctl(H0)" = H0,
"H0 w/ margin" = H0.prime,
"Hazard ratio(HR)" = HR,
"HR w/ margin" = HR.prime,
"Follow up(y)" = t,
sig.level = sig.level,
power = power,
alternative = alternative,
note = NOTE,
method = METHOD),
class = "power.htest")
}
追跡期間を 5 年間、試験群の生存率を 0.8、対照群の生存率を 0.65 として、非劣性マージ
ンを 0.065（つまり対照群の生存率の 10％）、有意水準両側 5％、検出力 80％とすると、
Freedman の方法では各群 84 人、Shoenfeld の方法では各群 75 人必要と計算される。
> sample.size.cox.noninf(t=5, S1=0.8, S0=0.65, Delta.S0=0.065,
+ alternative="two.sided")
Sample Size of Cox Model (noninferiority)
Death (Freedman) = 23.09411
Number (Freedman) = 83.97858
Death (Shoenfeld) = 20.42905
Number(Shoenfeld) = 74.28746
Survival trt(S1) = 0.8
Survival ctl(S0) = 0.65
Noninf. margin = 0.065
Hazard trt(H1) = 0.04462871
33
Hazard ctl(H0) = 0.08615658
H0 w/ margin = 0.1072287
Hazard ratio(HR) = 0.5179954
HR w/ margin = 0.4162012
Follow up(y) = 5
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group
追跡期間 10 年間、試験群と対照群の生存率をともに 0.95 として、非劣性マージンを 0.05
（つまり対照群の生存率の約 5％）
、有意水準片側 5％、検出力 80％とすると、Freedman
の方法では各群 520 人、Shoenfeld の方法では各群 478 人必要と計算される。
> sample.size.cox.noninf(t=10, S1=0.95, S0=0.95, Delta.S0=0.05,
+ alternative="one.sided")
Sample Size of Cox Model (noninferiority)
Death (Freedman) = 25.95087
Number (Freedman) = 519.0175
Death (Shoenfeld) = 23.86386
Number(Shoenfeld) = 477.2773
Survival trt(S1) = 0.95
Survival ctl(S0) = 0.95
Noninf. margin = 0.05
Hazard trt(H1) = 0.005129329
Hazard ctl(H0) = 0.005129329
H0 w/ margin = 0.01053605
Hazard ratio(HR) = 1
HR w/ margin = 0.486836
Follow up(y) = 10
sig.level = 0.05
power = 0.8
alternative = one.sided
34
NOTE: n is number in *each* group
参考文献
浜田＆藤井 生存時間解析における症例数設計 第 22 回 日本 SAS ユーザー会総会およ
び研究発表会 論文集 2003 年 7 月 31 日～8 月 1 日
東京
35
クロスオーバー試験
連続変数をエンドポイントにしたクロスオーバー試験のサンプルサイズ計算。個人差指数
が必要である。個人間分散と個人内分散の比であり、通常 1 より大きいはずであるが、こ
の特定が難しいと思われる。先行研究でわかれば問題ないが、わからない場合は、1 で計算
すれば、控えめで無難かと思う。
delta：effect size
sigma：標準偏差
theta：個人差指数
sigma.w：個人内標準偏差
rep.level：繰り返し計算を行う条件。デフォルトでは計算結果の差の割合が 10-8 未満
になるまで繰り返す。
samplesize.cross
<-
function(alpha,
power,
delta,
rep.level=10^-8){
p.t.t <- power.t.test(power=power, delta=delta/sigma)
np <- p.t.t$n
nc <- np/(2*(1+theta**2))
nc.approx <- nc
sigma.w <- sqrt(sigma**2/(1+theta**2))
nc0 <- nc
nc1 <- nc+5
yy <- matrix(0,8,1)
i <- 0
diff <- abs(nc1-nc0)/nc0
while (diff > rep.level){
i <- i+1
nc <- nc0
ta <- qt(p=1-alpha/2, df=2*nc-2)
tb <- qt(p=power, df=2*nc-2)
nc <- ((ta+tb)/(delta/sigma.w))**2
nc1 <- nc
diff <- abs(nc1-nc0)/nc0
yy[i,1] <- nc1
nc.exact <- nc1
sigma,
theta,
36
nc0 <- nc1
}
print(c("Alpha"=alpha, "Power"=power, "Delta"=delta, "Sigma"=sigma))
print(c("n(parllel)"=np,
"n(cross-over)"=nc.approx,
"n(cross-over
exact)"=nc.exact))
print(c(yy))
}
有意水準 5％、検出力 80％、effect size を 10、SD を 25 として、個人差指数を 1.5 とする
と一群 17 例必要。ちなみに並行群間試験だと一群 100 例必要。また、個人差指数を 1.0 と
すると一群 26 例必要。幾分控えめになる。
> samplesize.cross(alpha=.05, power=.8, delta=10, sigma=25, theta=1.5)
Alpha Power Delta Sigma
0.05 0.80 10.00 25.00
n(parllel)
n(cross-over) n(cross-over exact)
99.08057
15.24316
16.12026
[1] 16.18707 16.11550 16.12060 16.12024 16.12026 16.12026 16.12026 0.00000
> samplesize.cross(alpha=.05, power=.8, delta=10, sigma=25, theta=1)
Alpha Power Delta Sigma
0.05 0.80 10.00 25.00
n(parllel)
n(cross-over) n(cross-over exact)
99.08057
24.77014
25.53465
[1] 25.56813 25.53323 25.53471 25.53464 25.53465 25.53465 0.00000 0.00000
参考文献
丹後俊郎著 無作為化比較試験
クロスオーバー試験 p.182
朝倉書店
37
Matched case-control study
Matched case-control study のサンプルサイズ計算。
p0：コントロールの曝露割合
rr：期待されるオッズ比
ctrl：ケース 1 に対するコントロールの人数
sample.size.matched.ccs <- function(p0, rr, sig.level=.05, power=.8,
ctrl=1){
za <- qnorm(1-sig.level/2)
zb <- qnorm(power)
p <- rr/(1+rr)
m <- (za/2 + zb*sqrt(p*(1-p)))**2/(p-1/2)**2
p1 <- p0*rr/(1+p0*(rr-1))
q0 <- 1-p0
q1 <- 1-p1
pe <- p0*q1+p1*q0
pair <- (ctrl+1)*(m/pe)/(2*ctrl)
names(pair) <- rr
print(pair)
}
コントロールの曝露割合を 0.3、期待されるオッズ比を 2、検出力 90％、ケースとコントロ
ールの比を 1:1（デフォルト）とすると、ケースとコントロールそれぞれ 187 人必要と計算
される。
> sample.size.matched.ccs(p0=.3, rr=2, power=.9)
2
186.4132
ケースとコントロールの比を 1:2 とすると、必要なケースの人数が 140 人まで減少する。
この場合コントロールは 140 人×2＝280 人必要になる。
> sample.size.matched.ccs(p0=.3, rr=2, power=.9, ctrl=2)
2
38
139.8099
参考文献
Schlesselman JJ, Case-Control Studies, Oxford University Press, Oxford NY, 1982
6.6 Sample size and power for pair-matched studies pp.160-162
© Copyright 2024