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

# What is the Power of A and C in the Univariate Model if the MZ and DZ groups are of different sizes?
modN1 <- acePow(add = .33, com = .33, Nmz = 5000, Ndz = 1000)
modN2 <- acePow(add = .33, com = .33, Nmz = 3000, Ndz = 3000)
modN3 <- acePow(add = .33, com = .33, Nmz = 1000, Ndz = 5000)

par(mfrow = c(1,2),  mar = c(4, 5, 2, 1) + 0.1)
powerPlot(maxN = 1500, Wncp = c(modN1$WncpA, modN2$WncpA, modN3$WncpA))
legText <- c("MZ:DZ = 5:1", "MZ:DZ = 1:1", "MZ:DZ = 1:5")
legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Additive Genetic Variance", cex = 1.5)

powerPlot(maxN = 1500, Wncp = c(modN1$WncpC, modN2$WncpC, modN3$WncpC))
mtext("Common Environmental Variance", cex = 1.5)