# Power analysis for binomial test via simulation

Suppose a gene is mutated in 30% of the samples.

I want plot the number of samples required to see 30% of them mutated at various level of power.

I want to do this through a simulation so I generated different sample sizes and for each sample size I generated a 1000 trials of simulated number of mutated samples based on p.

p = 0.3
# from sample size 1:100, I generated 1000 binomial random variable
for (n in 1:100) {
k = rbinom(1000, n, p)
}


How do I construct the power analysis after I have generated some simulated data?

• Did you want to run some test? Eg, to show that the observed proportion is <.5? Or did you want to get a 95% CI of a given width & that would include .3 95% of the time? Mar 14, 2015 at 23:15
• I want to see how many samples would I need to see 30% samples mutated at 80% power. Mar 16, 2015 at 21:12
• What does that mean? Do you want to run a binomial test against a null of .5? Mar 16, 2015 at 21:15
• i think my null p0=.01 which a sample will be mutated in a gene in a background model, and my p1=.3 Mar 16, 2015 at 21:37
• i think my problem is quite similar to this question except my P0 is .01, so the minimum number of success to reject null should be some number greater than 1 i think. Mar 16, 2015 at 21:58

## 1 Answer

You may be able to figure out everything you need to know from my answer here: Simulation of logistic regression power analysis - designed experiments, which is quite comprehensive.

The basic procedure for simulating a power analysis is to:

1. Simulate data according to your preferred scenario (the alternative hypothesis).
2. Run the test you intend to use.
3. Do this many times, and see what percentage of the time your results are significant.
4. If you want to solve for the required $N$ to achieve a given power, do the above with different $N$s and find the one that yields the power you want.

In your case, you need to apply the binomial test with the null proportion set to $.01$.

Here's an example in R (note that your code had errors):

set.seed(8063)
p = 0.3
# from sample size 1:100, I generated 1000 binomial random variable
k = matrix(NA, nrow=1000, ncol=100)
for(n in 1:100){
k[,n] = rbinom(1000, n, p)
}

p.mat  = matrix(NA, nrow=1000, ncol=100)
for(i in 1:1000){
for(n in 1:100){
p.mat[i,n] = binom.test(x=k[i,n], n=n, p=0.01)$p.value } } power = c() power = apply(p.mat, 2, FUN=function(x){ mean(x<.05) }) min(which(power>.8)) # [1] 5 power[4] # [1] 0.763 power[5] # [1] 0.82  ### Brute force search: tl = NULL tl[paste("n", 1:10, sep="")] = list(NULL) pl = tl for(n in 1:10){ tl[[n]] = 0:n for(h in 0:n){ pl[[n]][h+1] = binom.test(tl[[n]][h+1], n, 0.01, "g")$p.value
}
}
msig = lapply(pl, function(x){ min(which(x<.05))-1 })
pow = c()
for(n in 1:10){
pow[n] = 1-sum(dbinom(0:(msig[[n]]-1), n, .3))
}
names(pow) = paste("n", 1:10, sep=" = ")
min(which(pow>.8))
# [1] 5
pow[4:5]
#   n = 4   n = 5
# 0.75990 0.83193

• Wouldnt it be binom.test(x=k[i,n], n=n, p=0.01, alternative="greater") in my case? Mar 17, 2015 at 20:06
• @kylelyk02, you could do that if you wanted. Mar 17, 2015 at 20:22
• i noticed that the power flucuate when sample is low, would it be safer to say power > .8 at power[9] instead of power[5]? Mar 17, 2015 at 20:51
• @kylelyk02, I don't follow that comment, but it occurs to me that your situation is so simple it can easily be found analytically by brute force search. I updated by answer. Note that the analytical results agree very closely w/ the original simulation. Mar 18, 2015 at 0:27
• here's what i meant by power fluctuate. As you can see there are occasional dips in the power as sample size increases. Mar 18, 2015 at 23:17