CODEWITHSUNDEEP
rloops
CadisEtRama
CadisEtRama

Reputation: 1111

loop for multiple power calculations with R

I need help writing a loop in R that will allow me to do a series of power calculations for different parameters. Alpha will always = 0.05 and N is the sample size that will take on various values so that: 

N= 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 1000.

H2 will also take on different values for each N analyzed, these are: 0.005, 0.01, 0. 02

So as an example, for N=1000, I would like to calculate the power with H2=0.005, H2=0.01 and H2= 0.02 at alpha=0.05. Then I want to do this analysis for all the values of N I mentioned above.

This is the code I am using below for just one run: 

N = 1000
alpha = 0.05
H2 = 0.005

threshold = qchisq(alpha, df = 1, lower.tail = FALSE)
power = pchisq(threshold, df = 1, lower.tail = FALSE, ncp = N * H2)

Can someone please help me turn this into a loop so I can get these results all at once in a structured table please? Thanks.

Upvotes: 0

Views: 1100

Answers (2)

Metrics
Metrics

Reputation: 15458

N <- c(1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000)
N<-as.list(rep(N,each=3)) # 3 is the length of original H2
H2 <- c(0.005, 0.01, 0.02)
H2<-as.list(rep(H2,10)) # 10 is the length of original N

myresult<-Map(function(x,y) cbind(x,y,power=pchisq(qchisq(0.05, df = 1, lower.tail = FALSE), df = 1, lower.tail = FALSE, ncp = x * y)),N,H2)
myout<-do.call(rbind,myresult)
colnames(myout)[1:2]<-c("N","H2")
> head(myout,10) 
         N    H2  power
 [1,] 1000 0.005 0.6088
 [2,] 1000 0.010 0.8854
 [3,] 1000 0.020 0.9940
 [4,] 2000 0.005 0.8854
 [5,] 2000 0.010 0.9940
 [6,] 2000 0.020 1.0000
 [7,] 3000 0.005 0.9721
 [8,] 3000 0.010 0.9998
 [9,] 3000 0.020 1.0000
[10,] 4000 0.005 0.9940

Upvotes: 0

Roland
Roland

Reputation: 132969

N <- c(1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000)
H2 <- c(0.005, 0.01, 0.02)
alpha <- 0.05

threshold <- qchisq(alpha, df = 1, lower.tail = FALSE)

paraComb <- expand.grid(N, H2)
ncp <- with(paraComb, Var1 * Var2)

setNames(cbind(paraComb,
               sapply(ncp, 
                      function(ncp) pchisq(threshold, 
                                           df = 1, lower.tail = FALSE, 
                                           ncp = ncp))
              ),
         c("N", "H2", "power"))

#        N    H2     power
# 1   1000 0.005 0.6087795
# 2   2000 0.005 0.8853791
# 3   3000 0.005 0.9721272
# 4   4000 0.005 0.9940005
# 5   5000 0.005 0.9988173
# <snip>

Upvotes: 1

Related Questions