CODEWITHSUNDEEP
rstatistics
trailblazer_1
trailblazer_1

Reputation: 133

how to calculate power value

I have simulated data in R and conducted a 2-sample t-test (2-sided hypothesis). I assumed an alpha value of 0.05, calculated p-values to determine if I accept or reject the null hypothesis. I need to now calculate power for each data set and plot. I need help on how to calculate power when alpha is known. More info on my data: I have 100 datasets of sample- rnorm(100,mean1,sd1) for each sample1 and sample 2.

Example data:

sample1 <- replicate(100,rnorm(100,0,1))
sample2 <- replicate(100,rnorm(100,5,3))

Upvotes: 0

Views: 561

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226027

I think you mean something like this:

extract a vector of p-values by running the test for each set of samples

alpha <- 0.05
nsim <- 100
set.seed(101)
sample1 <- replicate(nsim, rnorm(100,0,1))
sample2 <- replicate(nsim, rnorm(100,5,3))
pvalues <- rep(NA, nsim)
for (i in 1:nsim) {
   tt <- t.test(sample1[,i], sample2[,i]
   pvalues[i] <- tt$p.value
}

compute power (proportion of the time that p<alpha)

You can do this with sum(pvalues<alpha)/nsim, but

mean(pvalues<alpha)

is slightly more compact. (pvalues < alpha is TRUE or FALSE; applying a numeric function like sum() or mean() to it converts FALSE to 0 and TRUE to 1; taking the mean computes the proportion of values for which pvalues < alpha, which is exactly the observed probability of rejecting the null hypothesis → power.

As @neilfws points out, you can find the power for a t-test much more efficiently with power.t.test() (however, simulation is still useful for handling unusual cases that don't fit into the standard framework).

Upvotes: 3

Related Questions