2
$\begingroup$

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:

enter image description here

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?

$\endgroup$
1

1 Answer 1

0
$\begingroup$

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)
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – TMOTTM
    Aug 25, 2016 at 19:05
  • $\begingroup$ @TMOTTM, yes - see update (it is p <- sapply(new, function(x) t.test(rnorm(n=nsamples, mean = x, sd = sqrt(2)/10))$p.value)). $\endgroup$
    – Qaswed
    Aug 26, 2016 at 8:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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