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

# What is the Power to detect Additive Genetic Variance in a Common Pathway Twin Model?

# Simulate data for the Common Pathway Model 
ComPathDat <- comPathSim(aL=.2, cL=.2, eL=.6, load = c(.9,.8,.7,.6,.5), 
                         specA = c(.2, .2, .2, .2,.2), specC = c(.2, .2, .2, .2,.2), specE = c(.2, .2, .2, .2,.2), 
						 Nmz = 2000, Ndz = 2000)

# Put the data into MZ and DZ objects
mzData <- ComPathDat$mzData	
dzData <- ComPathDat$dzData	
	
# Set up objects for the analysis
selVars <- colnames(mzData)            
nf  <- 1       # number of factors
ntv <-  dim(mzData)[2]
nv <- ntv/2
nl        <- 1
vars <- paste("V", 1:nv, sep = "")

# Create Mean Matrices
meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=0, labels=paste("mean",vars,sep="_"), name="meanG" )

# Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s)
pathAl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels="aPath", lbound=.00001, name="al" )
pathCl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels="cPath", lbound=.00001, name="cl" )
pathEl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels="ePath", lbound=.00001, name="el" )

# Matrix and Algebra for constraint on variance of latent phenotype
covarLP   <- mxAlgebra( expression=al %*% t(al) + cl %*% t(cl) + el %*% t(el), name="CovarLP" )
varLP     <- mxAlgebra( expression= diag2vec(CovarLP), name="VarLP" )
unit      <- mxMatrix( type="Unit", nrow=nl, ncol=1, name="Unit")
varLP1    <- mxConstraint( expression=VarLP == Unit, name="varLP1")

# Matrix f for factor loadings on latent phenotype
pathFl    <- mxMatrix( type="Full", nrow=nv, ncol=nl, free=TRUE, values=.2, labels=paste("load",1:nv,sep = ""), name="fl" )

# Matrices as, cs, and es to store a, c, and e path coefficients for specific factors
pathAs    <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.4, labels=paste("as",1:nv, sep = ""), lbound=.00001, name="as" )
pathCs    <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.4, labels=paste("cs",1:nv, sep = ""), lbound=.00001, name="cs" )
pathEs    <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=paste("es",1:nv, sep = ""), lbound=.00001, name="es" )

# Matrices A, C, and E compute variance components
covA      <- mxAlgebra( expression=fl %&% (al %*% t(al)) + as %*% t(as), name="A" )
covC      <- mxAlgebra( expression=fl %&% (cl %*% t(cl)) + cs %*% t(cs), name="C" )
covE      <- mxAlgebra( expression=fl %&% (el %*% t(el)) + es %*% t(es), name="E" )

# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP      <- mxAlgebra( expression= A+C+E, name="V" )
covMZ     <- mxAlgebra( expression= A+C, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" )
expCovMZ  <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ  <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )

# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML     <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars      <- list(meanG, pathAl, pathCl, pathEl, covarLP, varLP, unit, pathFl, pathAs, pathCs, pathEs, covA, covC, covE, covP)
modelMZ   <- mxModel( name="MZ", pars, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( name="DZ", pars, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )

# Build & Run the Full Model 
modelCP   <- mxModel( "ComPath", pars, varLP1, modelMZ, modelDZ, multi )
fullCP     <- mxRun(modelCP)
summary(fullCP)

# Build & Run the Reduced Model 
redCP <- omxSetParameters(fullCP, labels = "aPath", value = 0, free = FALSE, name = "reduced")
redCP     <- mxRun(redCP)
summary(redCP)

# Calculate the Weighted NCP
neg2LL <- redCP$output$fit - fullCP$output$fit
totN   <- dim(mzData)[1]+dim(dzData)[1]
Wncp <- neg2LL / totN
Wncp

powerPlot(maxN = 2000, Wncp = c( 0.003413346, 0.00928087, 0.03469791)) 
legText <- c("A = .2, C = .2, & E = .6", "A = .3, C = .2, & E = .5", "A = .4, C = .2, & E = .4")
legend("bottomright", legend = legText, col = 1:3, ncol = 1, cex = 1, lwd = 3, box.lwd = 0)
mtext("Power to detect Additive Genetic \n Variance in a Latent Factor", cex = 1.5, line = 1)