# Simulating p-values as a function of sample size

We are trying to prove a very subtle effect occurring in cells after a certain treatment. Let's assume that the measurements are normally distributed. Let's also assume the untreated cells have $\mu = 1$ and $\sigma = 0.1$ and the treated cells have $\mu = 1.1$ and $\sigma = 0.22$. The question is:

How large must the sample size be in order for the observed effect to be statistically significant ($\alpha = 0.05$)?

I know that very subtile effects require a larger sample size than more apparent effects, but how many? I'm still learning statistics, so please be patient with me. I tried to perform a little simulation in R. Assuming that you randomly pick $n$ samples from a normal distribution, I tried to calculate the mean p-value as a function of $n$.

Is this a correct way to do find the right sample size? Or am I completely off the track with this approach?

Code:

library(ggplot2)

ctrl.mean <- 1
ctrl.sd <- 0.1
treated.mean <- 1.1
treated.sd <- 0.22

# Function that repeats t-test a number of times (rpt) with given sample size, means and sds.
# Returns a list of p-values from the test

tsim <- function(rpt, n, mean1, sd1, mean2, sd2) {
x <- 0
ppool <- NULL
while (x <= rpt) {
ppool <- c(ppool, t.test(rnorm(n,mean1,sd1), y = rnorm(n,mean2,sd2))$p.value) x <- x + 1 } return(ppool) } # Iterate through sample sizes and perform the function # Returns data frame with list of mean p-values at a given sample size i <- 2 num <- 50 res <- NULL while (i <= num) { sim <- tsim(1000, i, ctrl.mean, ctrl.sd, treated.mean, treated.sd) res <- rbind(res, cbind(i, mean(sim), sd(sim))) i <- i + 1 } # Plot the result res <- as.data.frame(res) ggplot(res, aes(x=i, y=-log10(V2))) + geom_line() + geom_ribbon(aes(ymin=-log10(V2)-log10(V3), ymax=-log10(V2)+log10(V3)), alpha = 0.2) + annotate("segment", x = 6, xend = num, y = -log10(0.05), yend = -log10(0.05), colour = "red", linetype = "dashed") + annotate("text", x = 0, y=-log10(0.05), label= "p = 0.05", hjust=0, size=3) + annotate("segment", x = 6, xend = num, y = -log10(0.01), yend = -log10(0.01), colour = "red", linetype = "dashed") + annotate("text", x = 0, y=-log10(0.01), label= "p = 0.01", hjust=0, size=3) + annotate("segment", x = 6, xend = num, y = -log10(0.001), yend = -log10(0.001), colour = "red", linetype = "dashed") + annotate("text", x = 0, y=-log10(0.001), label= "p = 0.001", hjust=0, size=3) + xlab("Number of replicates") + ylab("-log10(p-value)") + theme_bw()  ## 1 Answer You have almost performed what is usually called a power analysis. I say almost, because what you usually measure in a power calculation is not the mean p-value, but rather the probability that, given the sample size and the hypothesised mean difference, you would get a p-value lower than say 0.05. You can make small changes to your calculations in order to get this probability, however. The following script is a modification of your script that calculates the power for sample sizes from 2 to 50: ctrl.mean <- 1 ctrl.sd <- 0.1 treated.mean <- 1.1 treated.sd <- 0.22 n_range <- 2:50 max_samples <- 50 power <- NULL p.theshold <- 0.05 rpt <- 1000 for(n in n_range) { pvals <- replicate(rpt, { t.test(rnorm(n,ctrl.mean, ctrl.sd), y = rnorm(n, treated.mean, treated.sd))$p.value
})
power <- rbind(power, mean(pvals < p.theshold) )
}

plot(n_range, power, type="l", ylim=c(0, 1))


The way I would read this graph goes like: "Given my assumptions of the two groups, the probability that I would find a significant effect at n = 30 is roughly 50%". Often an 80% chance of finding an actual effect is considered a high level of power. By the way, power analysis is generally considered a good thing. :)

• Thank you for your answer! I was surprised to find out that I did a kind-of power analysis without knowing. :-) Statistics is more intuitive than I thought. – Eekhoorn Oct 17 '13 at 8:43
• I would say it is often intuitive, but unfortunately hidden by a thick coating of arithmetics and mathematical notation... :) – Rasmus Bååth Oct 17 '13 at 8:46
• @RasmusBååth, hey, I´m looking for improving a code inspired on your reply. Could you please check it at codereview.stackexchange.com/questions/215958/… – Luis Mar 21 at 21:16