# How to determine required difference for significance in R

I wrote a small script that generates a distribution with mean m1 and a second distribution with mean m2, where m2 is in the vicinity diff of m1.

m1   <- 2.50
diff <- 0.1
m2   <- seq(m1 - diff, m1 + diff, by=diff/15)
nsamples    <- 250

d <- data.frame("m1"=rep(m1, length(m2)), "m2"=m2, "pv"=rep(0, length(m1)))

for(i in 1:nrow(d)) {
d$pv[i] <- d$pv[i] + t.test(
rnorm(n=nsamples, mean=d$m1[i], sd=0.1)-rnorm(n=nsamples, mean=d$m2[i], sd=0.1)
)$p.value }  I'm interested now in how the p-value behaves depending on diff and nsamples, so I plotted it for the given parameters: Am I interpreting correctly (i.e. does my code do the right thing?) that from diff around 2.455 and 2.555 the p-value reaches significance? If its not too unrelated, is the way I'm setting up the d data.frame to accumulate the p-values done in a good way? Do I really have to prepare the m1 column with all the same value 2.5 in advance? • aside from your question, remember to correct for multiple comparisons. Aug 16, 2016 at 19:27 ## 1 Answer Your code can be simplified at many points. I'll go them trough step by step. m1 <- 2.50 diff <- 0.1 m2 <- seq(m1 - diff, m1 + diff, by=diff/15)  If you plug in m1 and diff you get m2 <- seq(2.4, 2.6, by = 0.1/15). I generaly prefer using length instead of by as long as you don't net this exact stepwidth (and your question doesn't gives the Impression that this is the case). So I would write m2 <- seq(2.4, 2.6, length = 31). nsamples <- 250 # is fine. d <- data.frame("m1"=rep(m1, length(m2)), "m2"=m2, "pv"=rep(0, length(m1))) # if you are only interesered in the plot shown, you don't need a data.frame.  About the loop: for(i in 1:nrow(d)) { d$pv[i] <- d$pv[i] + t.test( rnorm(n=nsamples, mean=d$m1[i], sd=0.1)-rnorm(n=nsamples, mean=d$m2[i], sd=0.1) )$p.value
}


The Expression d$pv[i] <- d$pv[i] + t.test(...) is futile as you initialise d$pv[i] with 0, so you calculate d$pv[i] <- 0 + t.test(...). Just write d$pv[i] <- t.test(...). Instead of building the difference from two independent normal variable draws, you can draw from the distribution of the differences: rnorm(n = nsamples, mean = d$m1[i] - d$m2[i], sd = sqrt(0.1^2 + 0.1^2)). If we recapitulate what d$m1 and d$m2 is we see that d$m1[i] - d$m2[i] is 2.5 - seq(2.4, 2.6, length = 31) which is nothing else than seq(-0.1, 0.1, length = 31). So what I would do is the following: nsamples <- 250 diff <- 0.1 new <- seq(-diff , diff, length = 31) p <- array(dim = length(new)) for(i in 1:length(new)) { p[i] <- t.test(rnorm(n=nsamples, mean = new[i], sd = sqrt(2)/10))$p.value
}
plot(x = new, y = p, type = "l")


Edit 2016-08-26 after request in comments: When one wants to avoid the loop, one could use sapply:

p <- sapply(new, function(x) t.test(rnorm(n=nsamples, mean = x, sd = sqrt(2)/10))$p.value)  • good points. Is there a way to avoid the for loop? I ask because many operations in R are vectorized, so I wonder if I can use this. Aug 25, 2016 at 19:05 • @TMOTTM, yes - see update (it is p <- sapply(new, function(x) t.test(rnorm(n=nsamples, mean = x, sd = sqrt(2)/10))$p.value)). Aug 26, 2016 at 8:27