### Test Analysis, Version 2.2
### edited by ABHW
### May 20, 2011
### Using Correlation Input, simplified matrix creation
###

#########
#########
#### Bivariate Analysis
#########
#########


require(OpenMx)
source("GenEpiHelperFunctions.R")

Vars <- c("traita", "traitb")                         ## Variable names here
Vars

nv <- 2
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")   
ntv <- nv*2

mzCov <- matrix(c(	1.0,	.26,	-.89,	.17,
					.26,	1.0,	.15,	.37,
					-.89,	.15,	1.0,	.29,
					.17,	.37,	.29,	1.0), 
      4,                                          ## Enter the total number of variables here - eg 2 vars * 2 twins
	dimnames=list(selVars,selVars))             ## Enter the dimnames here - use the selVars Command above

dzCov <- matrix(c(	1.0,	.26,	.27,	-.16,
					.26,	1.0,	.07,	.27,
					.27,	.07,	1.0,	.09,
					-.16,	.27,	.09,	1.0), 
      4,                                          ## Enter the total number of variables here - eg 2 vars * 2 twins
	dimnames=list(selVars,selVars))             ## Enter the dimnames here - use the selVars Command above




# Fit Multivariate ACE Model using Cholesky Decomposition with Correlations
# -----------------------------------------------------------------------

multACEModel <- mxModel("multACE",
    mxModel("ACE",
    # Matrices a, c, and e to store a, c, and e path coefficients
	mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("a11", "a21","a22"), name="a" ),
    mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("c11", "c21","c22"), name="c" ),
	mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("e11", "e21","e22"), name="e" ),
    
# Matrices A, C, and E compute variance components
        mxAlgebra( expression=a %*% t(a), name="A" ),
        mxAlgebra( expression=c %*% t(c), name="C" ),
        mxAlgebra( expression=e %*% t(e), name="E" ),
    
# Algebra to compute total variances and standard deviations (diagonal only)
        mxAlgebra( expression=A+C+E, name="V" ),
        mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
        mxAlgebra( expression=solve(sqrt(I*V)), name="isd"),
   
## Note that the rest of the mxModel statements do not change for bivariate/multivariate case
   
# Matrix & Algebra for expected means vector
#        mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= .5, name="M" ),
#        mxAlgebra( expression= cbind(M,M), name="expMean"),
    
# Algebra for expected variance/covariance matrix in MZ
        mxAlgebra( expression= rbind  ( cbind(A+C+E , A+C),
                                        cbind(A+C   , A+C+E)), name="expCovMZ" ),
    
# Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
        mxAlgebra( expression= rbind  ( cbind(A+C+E     , 0.5%x%A+C),
                                        cbind(0.5%x%A+C , A+C+E)),  name="expCovDZ" )         
),
  
mxModel("MZ",
        mxData( observed=mzCov, type="cor", numObs = 356),
        mxMLObjective( covariance="ACE.expCovMZ", dimnames=selVars )
    ),
    
mxModel("DZ", 
        mxData( observed=dzCov, type="cor", numObs = 240),
        mxMLObjective( covariance="ACE.expCovDZ", dimnames=selVars ) 
    ),
    mxAlgebra( expression=MZ.objective + DZ.objective, name="2sumll" ),
    mxAlgebraObjective("2sumll")
   )

multACEFit <- mxRun(multACEModel, intervals=F)
multACESumm <- summary(multACEFit)
multACESumm

# Fit AE model
# ---------------------------------------------------------------------
multAEModel <- mxModel(multACEFit, name="multAE",
    mxModel(multACEFit$ACE,
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=FALSE, values=0, name="c" ) # drop c at 0
    )
)
multAEFit <- mxRun(multAEModel)
multAESumm <- summary(multAEFit)
multAESumm
LL_AE <- mxEval(objective, multAEFit)
LL_AE


# Fit CE model
# --------------------------------------------------------------------
multCEModel <- mxModel(multACEFit, name="multCE",
    mxModel(multACEFit$ACE,
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=FALSE, values=0, name="a" ) # drop a at 0
    )
)
multCEFit <- mxRun(multCEModel)
multCESumm <- summary(multCEFit)
multCESumm
LL_CE <- mxEval(objective, multCEFit)
LL_CE


# Fit E model
# ---------------------------------------------------------------------
multEModel <- mxModel(multAEFit, name="multE",
    mxModel(multACEFit$ACE,
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=FALSE, values=0, name="a" ), # drop a at 0
        mxMatrix( type="Full", nrow=nv, ncol=nv, free=FALSE, values=0, name="c" )
    )
)
multEFit <- mxRun(multEModel)
multESumm <- summary(multEFit)
multESumm
LL_E <- mxEval(objective, multEFit)
LL_E


# Print Comparative Fit Statistics ACE models
# ---------------------------------------------------------------------
multACENested <- list(multAEFit, multCEFit, multEFit)
tableFitStatistics(multACEFit,multACENested)