require(OpenMx)
require(MASS)
mxOption(NULL, 'Default optimizer', 'NPSOL')
source("powerFun.R")



# What is the Power of A the Univariate Model for Binomial Data (relative to a continuous variable)
con  <- acePow( add  = .33, com  = .33, Nmz  = 1000, Ndz  = 1000)
bin5 <- acePowOrd( add  = .33, com  = .33,  percent = c(.5, .5), Nmz  = 1000, Ndz  = 1000)
bin4 <- acePowOrd( add  = .33, com  = .33,  percent = c(.4, .6), Nmz  = 1000, Ndz  = 1000)
bin3 <- acePowOrd( add  = .33, com  = .33,  percent = c(.3, .7), Nmz  = 1000, Ndz  = 1000)
bin2 <- acePowOrd( add  = .33, com  = .33,  percent = c(.2, .8), Nmz  = 1000, Ndz  = 1000)
bin1 <- acePowOrd( add  = .33, com  = .33,  percent = c(.9, .1), Nmz  = 1000, Ndz  = 1000)
bin05 <- acePowOrd( add  = .33, com  = .33,  percent = c(.05, .95), Nmz  = 1000, Ndz  = 1000)
 
powerPlot(maxN = 5000, Wncp = c(con$WncpA, bin5$WncpA, bin4$WncpA, bin3$WncpA, bin2$WncpA, bin1$WncpA, bin05$WncpA), wu = F)
legText <- c("Continuous", "Binary (Prevalence = .5)","Binary (Prevalence = .4)","Binary (Prevalence = .3)",
             "Binary (Prevalence = .2)", "Binary (Prevalence = .1)", "Binary (Prevalence = .05)")
legend("bottomright", legend = legText, col = 1:7, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)