#install.packages('clusterPower')
library(clusterPower)

Set the parameters

m0 <- 40
#m_sd <- 20 # leaving m as constant for now, although should vary 
p_t <- .88*.71 + (1-.88)*.41
p_c <- .88*.66 + (1-.88)*.23 # this is the overall employment level for the STAR report 
p_c <- p_c * seq(.5,1,.1)

Set a number of possible ICC values

icc <- seq(0.1,0.4,0.05)

Build the scenario dataset

d <- cbind(rep(p_c,times=length(icc)), rep(icc,each=length(p_c)))
d <- cbind(d, rep(p_t, nrow(d)), rep(m0, nrow(d)))

Sample size calcs

d <- cbind(d, apply(d,1,function(x) crtpwr.2prop(n=x[4], p1=x[1], p2=x[3], icc=x[2])))

Plot the results

Explore the effect of lower values of m

m1 <- c(40,30,20)
d <- cbind(d[rep(seq_len(nrow(d)), length(m1)), ], rep(m1, each=nrow(d)))

# Re-calculate the power 
d <- cbind(d, apply(d,1,function(x) crtpwr.2prop(n=x[6], p1=x[1], p2=x[3], icc=x[2])))

# Restict to 50 clusters per arm or fewer 
# d <- d[d[,7]<=50, ]