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


# What is the Power to detect Qualitative (different genes) and Quantitative (same genes but to different extents) Sex Limitation?

# Qualitative Sex Limitation
qual9 <- sexLimPower(am = .5, cm = .1, af = .5, cf = .1, rg = .9  , Nmzm = 5000, Nmzf = 5000, Ndzm = 5000, Ndzf = 5000, Ndzo = 10000)
qual7 <- sexLimPower(am = .5, cm = .1, af = .5, cf = .1, rg = .7  , Nmzm = 5000, Nmzf = 5000, Ndzm = 5000, Ndzf = 5000, Ndzo = 10000)
qual5 <- sexLimPower(am = .5, cm = .1, af = .5, cf = .1, rg = .5  , Nmzm = 5000, Nmzf = 5000, Ndzm = 5000, Ndzf = 5000, Ndzo = 10000)
qual3 <- sexLimPower(am = .5, cm = .1, af = .5, cf = .1, rg = .3  , Nmzm = 5000, Nmzf = 5000, Ndzm = 5000, Ndzf = 5000, Ndzo = 10000)
qual1 <- sexLimPower(am = .5, cm = .1, af = .5, cf = .1, rg = .1  , Nmzm = 5000, Nmzf = 5000, Ndzm = 5000, Ndzf = 5000, Ndzo = 10000)

# Quantitative Sex Limitation
quant1 <- sexLimPower(am = .2, cm = .2, af = .5, cf = .2, rg = 1, Nmzm = 5000, Nmzf = 5000, Ndzm = 5000, Ndzf = 5000, Ndzo = 5000)
quant2 <- sexLimPower(am = .3, cm = .2, af = .5, cf = .2, rg = 1, Nmzm = 5000, Nmzf = 5000, Ndzm = 5000, Ndzf = 5000, Ndzo = 5000)
quant3 <- sexLimPower(am = .4, cm = .2, af = .5, cf = .2, rg = 1, Nmzm = 5000, Nmzf = 5000, Ndzm = 5000, Ndzf = 5000, Ndzo = 5000)


par(mfrow = c(1,2),  mar = c(4, 5, 2, 1) + 0.1)

powerPlot(maxN = 15000, Wncp = c(qual9$WncpRg, qual7$WncpRg, qual5$WncpRg, qual3$WncpRg, qual1$WncpRg))
legText <- c("rg = .9", "rg = .7", "rg = .5", "rg = .3", "rg = .1")
legend("bottomright", legend = legText, col = 1:5, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Qualitative Sex Limitation", cex = 1)
# 41738 twin pairs required for 80% power for rg=.9

powerPlot(maxN = 5000, Wncp = c(quant1$WncpEqualA, quant2$WncpEqualA, quant3$WncpEqualA ))
legText <- c("Am = .2 & Af = .5", "Am = .3 & Af = .5", "Am = .4 & Af = .5")
legend("bottomright", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Quantitative Sex Limitation", cex = 1)