5
$\begingroup$

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$.

enter image description here

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()
$\endgroup$
9
$\begingroup$

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))

power analysis

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. :)

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Eekhoorn Oct 17 '13 at 8:43
  • $\begingroup$ I would say it is often intuitive, but unfortunately hidden by a thick coating of arithmetics and mathematical notation... :) $\endgroup$ – Rasmus Bååth Oct 17 '13 at 8:46
  • $\begingroup$ @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/… $\endgroup$ – Luis Mar 21 at 21:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.