5
$\begingroup$

Not sure if I'm misunderstanding the pwr::pwr.t.test() function and result.

set.seed(1299448)
g1 <- rnorm(35, 12, 3)
g2 <- rnorm(35, 14, 3)
t.test(g1, g2)
pwr.t.test(d = ((mean(as.numeric(g1), na.rm = T) -
                 mean(as.numeric(g2), na.rm = T)) /
                 sqrt(((sd(g1)^2)+ (sd(g2)^2))/2)),
           power = .80,
           sig.level = .05,
           alternative = "two.sided")

Results

Welch Two Sample t-test

data:  g1 and g2
t = -2.1738, df = 63.822, p-value = 0.03343


Two-sample t test power calculation 

          n = 59.10848
          d = 0.519647
  sig.level = 0.05
      power = 0.8
alternative = two.sided

NOTE: n is number in *each* group

So the effect is significant with a sample size of 70, but with such an effect size, a sample size of 120 would be required for 80% power.

Can someone ELI5? Sorry if this is more of a cross-validated question. I started originally thinking this was a programming/misunderstanding of the pwr package, but now I'm starting to think I just don't understand the sample size estimation/power well enough.

$\endgroup$
1
  • 2
    $\begingroup$ This seems much more about misunderstanding about statistical concepts than programming. I've voted to migrate to Cross Validated. $\endgroup$
    – MrFlick
    Dec 1, 2015 at 15:10

1 Answer 1

1
$\begingroup$

The issue is that the effect size for that particular realization of the simulation is different than the effect size of the model that generated it, and just because the effect is strong enough to produce a p-value less than 0.05 one time, doesn't mean that it's strong enough to create a p-value less than 0.05 80% of the time:

set.seed(1299448)
g1 <- rnorm(35, 12, 3)
g2 <- rnorm(35, 14, 3)
(eff <- ((mean(as.numeric(g1), na.rm = T) -
           mean(as.numeric(g2), na.rm = T)) /
          sqrt(((sd(g1)^2)+ (sd(g2)^2))/2)))
# [1] -0.519647

n.samp <- 10000
out70 <- rep(NA,n.samp)
out120 <- rep(NA,n.samp)

for(i in 1:n.samp) {
  g1 <- rnorm(35, 0)
  g2 <- rnorm(35, eff)
  out70[i] <- t.test(g1, g2)$p.value<0.05
  g1 <- rnorm(60, 0)
  g2 <- rnorm(60, eff)
  out120[i] <- t.test(g1, g2)$p.value<0.05
}

mean(out70)
# [1] 0.5682
mean(out120)
# [1] 0.8109
$\endgroup$

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.