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

# What is the Power to detect Rg in the Bivariate Model?
rg1 <- bivPow(A1 = .3, A2 = .3, rg = .10, C1 = .33, C2 = .33, rc = .1, Nmz = 1000, Ndz = 1000)
rg2 <- bivPow(A1 = .3, A2 = .3, rg = .30, C1 = .33, C2 = .33, rc = .3, Nmz = 1000, Ndz = 1000)
rg3 <- bivPow(A1 = .3, A2 = .3, rg = .50, C1 = .33, C2 = .33, rc = .5, Nmz = 1000, Ndz = 1000)

rg4 <- bivPow(A1 = .4, A2 = .4, rg = .10, C1 = .33, C2 = .33, rc = .1, Nmz = 1000, Ndz = 1000)
rg5 <- bivPow(A1 = .4, A2 = .4, rg = .30, C1 = .33, C2 = .33, rc = .3, Nmz = 1000, Ndz = 1000)
rg6 <- bivPow(A1 = .4, A2 = .4, rg = .50, C1 = .33, C2 = .33, rc = .5, Nmz = 1000, Ndz = 1000)

rg7 <- bivPow(A1 = .5, A2 = .5, rg = .10, C1 = .33, C2 = .33, rc = .1, Nmz = 1000, Ndz = 1000)
rg8 <- bivPow(A1 = .5, A2 = .5, rg = .30, C1 = .33, C2 = .33, rc = .3, Nmz = 1000, Ndz = 1000)
rg9 <- bivPow(A1 = .5, A2 = .5, rg = .50, C1 = .33, C2 = .33, rc = .5, Nmz = 1000, Ndz = 1000)

par(mfrow = c(3,2),  mar = c(4, 5, 2, 1) + 0.1)
powerPlot(maxN = 750, Wncp = c(rg1$WncpA1, rg2$WncpA1, rg3$WncpA1))
legText <- c("A = .3 & rg = .1", "A = .3 & rg = .3", "A = .3 & rg = .5")
legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Additive Genetic Variance", cex = 1)

powerPlot(maxN = 2000, Wncp = c(rg1$WncpRg, rg2$WncpRg, rg3$WncpRg))
mtext("Genetic Correlation", cex = 1)

powerPlot(maxN = 750, Wncp = c(rg4$WncpA1, rg5$WncpA1, rg6$WncpA1))
legText <- c("A = .4 & rg = .1", "A = .4 & rg = .3", "A = .4 & rg = .5")
legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Additive Genetic Variance", cex = 1)

powerPlot(maxN = 2000, Wncp = c(rg4$WncpRg, rg5$WncpRg, rg6$WncpRg))
mtext("Genetic Correlation", cex = 1)

powerPlot(maxN = 750, Wncp = c(rg7$WncpA1, rg8$WncpA1, rg9$WncpA1))
legText <- c("A = .5 & rg = .1", "A = .5 & rg = .3", "A = .5 & rg = .5")
legend("right", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Additive Genetic Variance", cex = 1)

powerPlot(maxN = 2000, Wncp = c(rg7$WncpRg, rg8$WncpRg, rg9$WncpRg))
mtext("Genetic Correlation", cex = 1)