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


# What is the Power of A the Univariate Model as C increases?
modA1 <- acePow(add = .33, com = .1, Nmz = 1000, Ndz = 1000)
modA2 <- acePow(add = .33, com = .3, Nmz = 1000, Ndz = 1000)
modA3 <- acePow(add = .33, com = .5, Nmz = 1000, Ndz = 1000)

par(mfrow = c(1,2), mar = c(4, 5, 2, 1) + 0.1)
powerPlot(maxN = 1500, Wncp = c(modA1$WncpA, modA2$WncpA, modA3$WncpA))
legText <- c("A = .33 & C = .1", "A = .33 & C = .3", "A = .33 & C = .5")
legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Additive Genetic Variance", cex = 1.5)

# What is the Power of C the Univariate Model as A increases?
modC1 <- acePow(add = .1, com = .33, Nmz = 1000, Ndz = 1000)
modC2 <- acePow(add = .3, com = .33, Nmz = 1000, Ndz = 1000)
modC3 <- acePow(add = .5, com = .33, Nmz = 1000, Ndz = 1000)

#par( mar = c(4, 5, 2, 1) + 0.1)
powerPlot(maxN = 1500, Wncp = c(modC1$WncpC, modC2$WncpC, modC3$WncpC))
legText <- c("A = .1 & C = .33", "A = .3 & C = .33", "A = .5 & C = .33")
legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Common Environmental Variance", cex = 1.5)