simulated power analysis for difference in proportions

power sim logit.R

#  ------------------------------------------------------------------
#  |PROGRAM NAME: power sim logit.R
#  |DATE: 4/29/18
#  |CREATED BY: MATT BOGARD 
#  |PROJECT FILE: /Users/amandabogard/Google Drive/R Training
#  |----------------------------------------------------------------
#  | PURPOSE: compare modified simulated power analysis for difference in proportions from egap.org
#  | to results from the 'pwr' function in R
#  |----------------------------------------------------------------

# reference: http://egap.org/content/power-analysis-simulations-r 

# rm(list=ls()) # get rid of any existing data 

options("scipen" =100, "digits" = 4) # override R's tendency to use scientific notation

# cat("\f") # clear console

#-------------------------------------------------
# simulated power and sample size calculations
#-------------------------------------------------

possible.ns <- seq(from=100, to=500, by=10)     # The sample sizes we'll be considering
powers <- rep(NA, length(possible.ns))           # Empty object to collect simulation estimates
alpha <- 0.05                                    # Standard significance level
sims <- 500                                      # Number of simulations to conduct for each N

#### Outer loop to vary the number of subjects ####
for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  
  significant.experiments <- rep(NA, sims)         # Empty object to count significant experiments
  
  #### Inner loop to conduct experiments "sims" times over for each N ####
  for (i in 1:sims){
    Y0 <-  rbinom(n = N, 1,.5)                          # control potential outcome proportion
    Y1 <- rbinom(n = N, 1,.25)                          # treatment potential outcome proportion
    Z.sim <- rbinom(n=N, size=1, prob=.5)           # random treatment assignment
    Y.sim <- Y1*Z.sim + Y0*(1-Z.sim)                  # reveal outcomes according to assignment
    fit.sim <- glm( Y.sim ~ Z.sim,family="binomial")  # Do analysis (logistic regression)
    p.value <- summary(fit.sim)$coefficients[2,4]  # Extract p-values
    significant.experiments[i] <- (p.value <= alpha) # Determine significance according to p <= 0.05
  }
  
  powers[j] <- mean(significant.experiments)       # store average success rate (power) for each N
}
plot(possible.ns, powers, ylim=c(0,1))

tmp <- cbind(possible.ns, powers)

print(tmp) # n = 120 power ~ .80 (this matches what we get using Gpower)

# cleanup
rm("alpha", "fit.sim","i","j","N","p.value","possible.ns","powers","significant.experiments","sims",
   "tmp","Y.sim","Y0","Y1","Z.sim")

#--------------------------------------------
# compare this to the power function in R
#--------------------------------------------

library(pwr)

# reference: https://www.statmethods.net/stats/power.html

# syntax
# pwr.2p.test(h = , n = , sig.level =, power = )
# where h is the effect size and n is the common sample size in each group.

# h = 2*arcsin(sqrt(p1)) - 2* arcsin(sqrt(p2))

2*asin(sqrt(.5)) - 2*asin(sqrt(.25)) # calculate h = .5236

pwr.2p.test(h = .5236, n = , sig.level = .05, power = .80) # ~ 58*2 = 116 similar to results from simulation