Simulation for sample size calculation

3 views
Skip to first unread message

JonStat

unread,
Nov 13, 2009, 5:09:49 AM11/13/09
to MedStats
Dear all,

I often get asked for simple sample size calculations, e.g. for a
difference between two independent means or two independent
proportions, however in the pursuit of greater accuracy I'd like to be
able to run simulations to get my figures instead of relying on the
standard formula. The trouble is I'm not sure of a (gentle) source of
guidance on how to go about this, so your help would be greatfully
appreciated in pointing me in the right direction (as you very kindly
have done before)

Kind regards,

Jon

Steve Simon, P.Mean Consulting

unread,
Nov 13, 2009, 12:51:19 PM11/13/09
to meds...@googlegroups.com
JonStat wrote:

Simulations, of course, do not provide a more exact answer than a
formula, when a formula is available. The value of simulations is that
you can stress your assumptions a bit by adding some non-normality,
heterogeneity of variances, and so forth. Under those situations, the
formulas for power may be approximate or may not even exist.

Here's some code in R that shows how to get power when there is a half a
standard deviation shift and the sample size is 64 in each group.

alpha <- 0.05
alt <- "two.sided"
mu1 <- 0
mu2 <- 0.5
sigma1 <- 1
sigma2 <- 1
n1 <- 64
n2 <- 64
nreps <- 1000
pv <- rep(NA,nreps)
for (i in 1:nreps) {
x1 <- rnorm(n1,mu1,sigma1)
x2 <- rnorm(n2,mu2,sigma2)
pv[i] <- t.test(x1,x2,paired=FALSE,alternative=alt)$p.value
}
power <- sum(pv<alpha)/nreps
lower.ci.power <- power-1.96*sqrt(power*(1-power)/nreps)
upper.ci.power <- power+1.96*sqrt(power*(1-power)/nreps)

It's trivial, I'm sure, to do this in SAS also, but I stopped using SAS
ten years ago.
--
Steve Simon, Standard Disclaimer
The Monthly Mean is celebrating its first anniversary.
Find out more about the newsletter that dares
to call itself "average" at www.pmean.com/news

Frank Harrell

unread,
Nov 13, 2009, 1:50:12 PM11/13/09
to MedStats
For two proportions, see the bpower.sim function in the R Hmisc
package. For more advanced simulation (survival 2-group comparison
with delayed treatment effect, dropout, etc.) see the Hmisc spower
function.

Frank
Reply all
Reply to author
Forward
0 new messages