#		Author: Frühling Rijsdijk
#		MODEL:	MV quantitative & qualitative sex-dif Sex-Limitation script (ACE Correlated Factors model)
#		Correlation Approach to ensure that the order of the variables does NOT affect the ability 
#		of the model to account for the DZOS data.
#		Ref: Neale et al., Multivariate genetic analysis of sex-lim and G x E interaction, Twin Research & Human Genetics, 2006 
#		Restrictions: Assumes means and variances can be equated across birth order within zygosity groups



# ------------------------------------------------------------------------------------------------------------------------------------------------------
# Variabels: T1SEX T2SEX T1AGE T2AGE T1M	T2M	T1F	T2F	T1T	T2T	ZYG ZYGSEX (1=MZ, 2=DZ) (0=male, 1=female), (1=mzm, 2=mzf, 3=dzm, 4=dzf, 5=dos)
# Three informants (Mother (M), father (F) and twins (T)), continuous measure.


setwd("C:/")   					#insert the path to your working directory. Remember Windows path format backslash must be changed to R format (single slash or double backslash )


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


# Read in data file 
# -----------------------------------------------------------------------

mv_data <- read.table("datafile.DAT ", header=TRUE, na.strings='-1.000')    #insert the name of your data file


nvar <- 3    					# number of variables
tnvar <- 2*nvar 					# number of variables*max family size
nlower <-tnvar*(tnvar+1)/2 			# number of free elements in a lower matrix tnvar*tnvar


MVdata <- mv_data[,c("T1M","T2M", "T1F", "T2F", "T1T", "T2T", "ZYG", "ZYGSEX", "T1SEX", "T2SEX", "T1AGE", "T2AGE")]

names(MVdata)

attach(MVdata)



Vars <- c('M','F','T')
selVars <- paste("T",c(rep(1,nvar),rep(2,nvar)),Vars,sep="")
print(selVars)
 
summary (MVdata)

summary (selVars)

mzmData <- as.data.frame (subset (MVdata, ZYGSEX ==1, selVars))
mzfData <- as.data.frame (subset (MVdata, ZYGSEX ==2, selVars))
dzmData <- as.data.frame (subset (MVdata, ZYGSEX ==3, selVars))
dzfData <- as.data.frame (subset (MVdata, ZYGSEX ==4, selVars))
dzoData <- as.data.frame (subset (MVdata, ZYGSEX ==5, selVars))

summary (mzmData)
summary (mzfData)
summary (dzmData)
summary (dzfData)
summary (dzoData)


# 1 Nonscalar Sex Limitation #####################################################################################
# Quantitative sex Diff & Qualitative sex Diff for 		A
# In this first model there are male and female paths, plus male and female Ra, Rc and Re between variables
# Male-Female correlations in DZO group between 		A factors FREE
# ----------------------------------------------------------------------------------------------------------------

ACENonScSLCoraModel <- mxModel("ACENonScSLCora",
    	mxModel("ACE",
    	# Matrices a, c, and e to store a, c, and e path coefficients
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("am11","am22", "am33"), name="am" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("cm11","cm22", "cm33"), name="cm" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("em11","em22", "em33"), name="em" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("af11","af22", "af33"), name="af" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("cf11","cf22", "cf33"), name="cf" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("ef11","ef22", "ef33"), name="ef" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("ram21", "ram31", "ram32"), lbound=-.9999, ubound=.9999, name="Ram" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("rcm21", "rcm31", "rcm32"), lbound=-.9999, ubound=.9999, name="Rcm" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("rem21", "rem31", "rem32"), lbound=-.9999, ubound=.9999, name="Rem" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("raf21", "raf31", "raf32"), lbound=-.9999, ubound=.9999, name="Raf" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("rcf21", "rcf31", "rcf32"), lbound=-.9999, ubound=.9999, name="Rcf" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("ref21", "ref31", "ref32"), lbound=-.9999, ubound=.9999, name="Ref" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=.5, lbound=-.5, ubound=.5, name="CorMFa" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=FALSE, values=1, name="CorMFc" ),

    	# Matrices A, C, and E compute variance components											
        	mxAlgebra( expression=am %*% (Ram) %*% t(am), name="Am" ),
        	mxAlgebra( expression=cm %*% (Rcm) %*% t(cm), name="Cm" ),
        	mxAlgebra( expression=em %*% (Rem) %*% t(em), name="Em" ),
        	mxAlgebra( expression=af %*% (Raf) %*% t(af), name="Af" ),
        	mxAlgebra( expression=cf %*% (Rcf) %*% t(cf), name="Cf" ),
        	mxAlgebra( expression=ef %*% (Ref) %*% t(ef), name="Ef" ),

    	# Algebra to compute total variances and standard deviations (diagonal only)
        	mxAlgebra( expression=Am+Cm+Em, name="Vm" ),
        	mxAlgebra( expression=Af+Cf+Ef, name="Vf" ),
        	mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
        	mxAlgebra( expression=solve(sqrt(I*Vm)), name="iSDm"),
        	mxAlgebra( expression=solve(sqrt(I*Vf)), name="iSDf"),


    # Matrix & Algebra for expected means vector
        	mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 18, name="Mm" ),
        	mxAlgebra(  cbind(Mm,Mm), name="expMeanm"),
        	mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 18, name="Mf" ),
        	mxAlgebra(  cbind(Mf,Mf), name="expMeanf"),
        	mxAlgebra(  cbind(Mm,Mf), name="expMeanOs"),

    # Algebra for expected variance/covariance matrix in MZ

		mxAlgebra(  rbind  ( cbind(Am+Cm+Em , 	Am+Cm),
                       			cbind(Am+Cm   , 		Am+Cm+Em)), name="expCovMZm" ),

		mxAlgebra(  rbind  ( cbind(Af+Cf+Ef , 	Af+Cf),
                                     cbind(Af+Cf   , 	Af+Cf+Ef)), name="expCovMZf" ),

    # Algebra for expected variance/covariance matrix in DZ

		mxAlgebra(  rbind  ( cbind(Am+Cm+Em     , 			0.5%x%Am+Cm),
                                     		cbind(0.5%x%Am+Cm, 		Am+Cm+Em)),  name="expCovDZm" ), 

		mxAlgebra(  rbind  ( cbind(Af+Cf+Ef     ,		0.5%x%(Af) + Cf ),
                                 		cbind(0.5%x%(Af) + Cf , 	Af+Cf+Ef)),  name="expCovDZf" ), 

# Algebra for expected variance/covariance matrix in DZOS

		mxAlgebra ( rbind ( cbind(Am+Cm+Em	,   am %*%(CorMFa%*%Raf)%*% t(af) + cm %*%(CorMFc%*%Rcf)%*% t(cf) ),
				 cbind(  af %*%(CorMFa%*%Ram)%*% t(am) + cf %*%(CorMFc%*%Rcm)%*% t(cm)	, 	Af+Cf+Ef) ), name="expCovOS"),



# Compute Variance Components
    		mxAlgebra( expression=rbind(cbind(Am/Vm,Cm/Vm,Em/Vm),cbind(Af/Vf,Cf/Vf,Ef/Vf)), name="VarComp")    


    		mxModel("MZm",
        		mxData( observed=mzmData, type="raw" ),
        		mxFIMLObjective( covariance="ACE.expCovMZm", means="ACE.expMeanm", dimnames=selVars )  ),
    		mxModel("DZm", 
        		mxData( observed=dzmData, type="raw" ),
        		mxFIMLObjective( covariance="ACE.expCovDZm", means="ACE.expMeanm", dimnames=selVars )  ),
    		mxModel("MZf",
        		mxData( observed=mzfData, type="raw" ),
        		mxFIMLObjective( covariance="ACE.expCovMZf", means="ACE.expMeanf", dimnames=selVars )  ),
    		mxModel("DZf", 
        		mxData( observed=dzfData, type="raw" ),
        		mxFIMLObjective( covariance="ACE.expCovDZf", means="ACE.expMeanf", dimnames=selVars ) ),
    		mxModel("OS", 
        		mxData( observed=dzoData, type="raw"),
 			mxFIMLObjective( covariance="ACE.expCovOS", means="ACE.expMeanOs", dimnames=selVars )  ),
    		mxAlgebra( expression=MZm.objective + DZm.objective + MZf.objective + DZf.objective + OS.objective, name="modelfit" ),
    		mxAlgebraObjective("modelfit")
# , mxCI (c('ACE.VarComp'))
)


#ACENonScSLCoraFit <- mxRun(ACENonScSLCoraModel, intervals=T)
ACENonScSLCoraFit <- mxRun(ACENonScSLCoraModel)


ACENonScSLCoraSumm <- summary(ACENonScSLCoraFit)
ACENonScSLCoraSumm

# Generate Quantitative sex Diff & Qualitative sex Diff for A ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(ACENonScSLCoraFit)
expectedMeansCovariances(ACENonScSLCoraFit)
tableFitStatistics(ACENonScSLCoraFit)
   ACEpathMatrices <- c("ACE.am","ACE.cm","ACE.em","ACE.iSDm","ACE.iSDm %*% ACE.am","ACE.iSDm %*% ACE.cm","ACE.iSDm %*% ACE.em", "ACE.af","ACE.cf","ACE.ef","ACE.iSDf","ACE.iSDf %*% ACE.af","ACE.iSDf %*% ACE.cf","ACE.iSDf %*% ACE.ef")
   ACEpathLabels <- c("pathEst_am","pathEst_cm","pathEst_em","sdm","stPathEst_am","stPathEst_cm","stPathEst_em", "pathEst_af","pathEst_cf","pathEst_ef","sdf","stPathEst_af","stPathEst_cf","stPathEst_ef")
   ACEcovMatrices <- c("ACE.Am", "ACE.Cm", "ACE.Em","ACE.Vm","ACE.Am/ACE.Vm", "ACE.Cm/ACE.Vm","ACE.Em/ACE.Vm", "ACE.Af", "ACE.Cf", "ACE.Ef","ACE.Vf","ACE.Af/ACE.Vf", "ACE.Cf/ACE.Vf","ACE.Ef/ACE.Vf")
   ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")

formatOutputMatrices(ACENonScSLCoraFit,ACEpathMatrices,ACEpathLabels,4)
formatOutputMatrices(ACENonScSLCoraFit,ACEcovMatrices,ACEcovLabels,4)

# Generate Table of Parameter Estimates using mxEval

pathEstimatesACE <- round(mxEval(cbind(ACE.am,ACE.cm,ACE.em,ACE.af,ACE.cf,ACE.ef), ACENonScSLCoraFit),4)
    rownames(pathEstimatesACE) <- c('V1pathEstimates', 'V2pathEstimates', 'V3pathEstimates')
    colnames(pathEstimatesACE) <- c('am1','am2','am3','cm1','cm2','cm3','em1','em2', 'em3','af1','af2', 'af3','cf1','cf2', 'cf3', 'ef1','ef2', 'ef3')				
pathEstimatesACE

varComponentsACE <- round(mxEval(cbind(ACE.Am/ACE.Vm,ACE.Cm/ACE.Vm,ACE.Em/ACE.Vm,ACE.Af/ACE.Vf,ACE.Cf/ACE.Vf,ACE.Ef/ACE.Vf), ACENonScSLCoraFit),2)
    rownames(varComponentsACE) <- c('V1varcomp','V2varcomp', 'V3varcomp')
    colnames(varComponentsACE) <- c('h2m1','h2m2', 'h2m3', 'c2m1','c2m2','c2m3', 'e2m1','e2m2','e2m3', 'h2f1','h2f2','h2f3', 'c2f1','c2f2', 'c2f3', 'e2f1', 'e2f2', 'e2f3')
varComponentsACE




# 2 Nonscalar Sex Limitation #####################################################################################
# Quantitative sex Diff & Qualitative sex Diff for 		C
# In the second model there are male and female paths, plus male and female Ra, Rc and Re between variables
# Male-Female correlations in DZO group between 		C factors 
# -------------------------------------------------------------------------------------------------------------------

ACENonScSLCorcModel <- mxModel (ACENonScSLCoraFit, name="ACENonScSLCorc",
    mxModel(ACENonScSLCoraFit$ACE,
    #  Specify the Male-Female C correlations:
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, lbound=-1, ubound=1, name="CorMFc" ),
    #  FIX the Male-Female A correlations:
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=FALSE, values=.5, name="CorMFa" )
            )
)
ACENonScSLCorcFit <- mxRun(ACENonScSLCorcModel, intervals=T)
ACENonScSLCorcSumm <- summary(ACENonScSLCorcFit)
ACENonScSLCorcSumm



# Generate Quantitative sex Diff & Qualitative sex Diff for C ACE  Output
# -----------------------------------------------------------------------

parameterSpecifications(ACENonScSLCorcFit)
expectedMeansCovariances(ACENonScSLCorcFit)
tableFitStatistics(ACENonScSLCorcFit)
   ACEpathMatrices <- c("ACE.am","ACE.cm","ACE.em","ACE.iSDm","ACE.iSDm %*% ACE.am","ACE.iSDm %*% ACE.cm","ACE.iSDm %*% ACE.em", "ACE.af","ACE.cf","ACE.ef","ACE.iSDf","ACE.iSDf %*% ACE.af","ACE.iSDf %*% ACE.cf","ACE.iSDf %*% ACE.ef")
   ACEpathLabels <- c("pathEst_am","pathEst_cm","pathEst_em","sdm","stPathEst_am","stPathEst_cm","stPathEst_em", "pathEst_af","pathEst_cf","pathEst_ef","sdf","stPathEst_af","stPathEst_cf","stPathEst_ef")
   ACEcovMatrices <- c("ACE.Am", "ACE.Cm", "ACE.Em","ACE.Vm","ACE.Am/ACE.Vm", "ACE.Cm/ACE.Vm","ACE.Em/ACE.Vm", "ACE.Af", "ACE.Cf", "ACE.Ef","ACE.Vf","ACE.Af/ACE.Vf", "ACE.Cf/ACE.Vf","ACE.Ef/ACE.Vf")
   ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")

formatOutputMatrices(ACENonScSLCoraFit,ACEpathMatrices,ACEpathLabels,4)
formatOutputMatrices(ACENonScSLCoraFit,ACEcovMatrices,ACEcovLabels,4)


# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(ACE.am,ACE.cm,ACE.em,ACE.af,ACE.cf,ACE.ef), ACENonScSLCorcFit),4)
    rownames(pathEstimatesACE) <- c('V1pathEstimates', 'V2pathEstimates', 'V3pathEstimates')
    colnames(pathEstimatesACE) <- c('am1','am2','am3','cm1','cm2','cm3','em1','em2', 'em3','af1','af2', 'af3','cf1','cf2', 'cf3', 'ef1','ef2', 'ef3')				
pathEstimatesACE

varComponentsACE <- round(mxEval(cbind(ACE.Am/ACE.Vm,ACE.Cm/ACE.Vm,ACE.Em/ACE.Vm,ACE.Af/ACE.Vf,ACE.Cf/ACE.Vf,ACE.Ef/ACE.Vf), ACENonScSLCorcFit),2)
    rownames(varComponentsACE) <- c('V1varcomp','V2varcomp', 'V3varcomp')
    colnames(varComponentsACE) <- c('h2m1','h2m2', 'h2m3', 'c2m1','c2m2','c2m3', 'e2m1','e2m2','e2m3', 'h2f1','h2f2','h2f3', 'c2f1','c2f2', 'c2f3', 'e2f1', 'e2f2', 'e2f3')
varComponentsACE




# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
tableFitStatistics(ACENonScSLCoraFit, ACENonScSLCorcFit)




# 3 Nonscalar Sex Limitation (Preliminary model) ##########################################################################################################
# Quantitative Sex Diff but NO Qualitative sex Diff
# In the third model there are male and female paths, plus male and female Ra, Rc and Re between variables
# but Male-Female correlations in DZO group between 	A & C factors are fixed to their theoretical values of .5 and 1, respectively 
# From the previous model, we simply fix the C cor to 1 as well
# -----------------------------------------------------------------------------------------------------------------------------------

ACENonScSLModel <- mxModel (ACENonScSLCorcFit, name="ACENonScSL",
    mxModel(ACENonScSLCorcFit$ACE,
    #  Specify the Male-Female C correlations:
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=FALSE, values=1, name="CorMFc" ),
    #  FIX the Male-Female A correlations:
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=FALSE, values=.5, name="CorMFa" )
            )
)
ACENonScSLFit <- mxRun(ACENonScSLModel)
ACENonScSLSumm <- summary(ACENonScSLFit)
ACENonScSLSumm


# Generate non-scalar sex Diff ACE  Output
# -----------------------------------------------------------------------

parameterSpecifications(ACENonScSLFit)
expectedMeansCovariances(ACENonScSLFit)
tableFitStatistics(ACENonScSLFit)
   ACEpathMatrices <- c("ACE.am","ACE.cm","ACE.em","ACE.iSDm","ACE.iSDm %*% ACE.am","ACE.iSDm %*% ACE.cm","ACE.iSDm %*% ACE.em", "ACE.af","ACE.cf","ACE.ef","ACE.iSDf","ACE.iSDf %*% ACE.af","ACE.iSDf %*% ACE.cf","ACE.iSDf %*% ACE.ef")
   ACEpathLabels <- c("pathEst_am","pathEst_cm","pathEst_em","sdm","stPathEst_am","stPathEst_cm","stPathEst_em", "pathEst_af","pathEst_cf","pathEst_ef","sdf","stPathEst_af","stPathEst_cf","stPathEst_ef")
   ACEcovMatrices <- c("ACE.Am", "ACE.Cm", "ACE.Em","ACE.Vm","ACE.Am/ACE.Vm", "ACE.Cm/ACE.Vm","ACE.Em/ACE.Vm", "ACE.Af", "ACE.Cf", "ACE.Ef","ACE.Vf","ACE.Af/ACE.Vf", "ACE.Cf/ACE.Vf","ACE.Ef/ACE.Vf")
   ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")

formatOutputMatrices(ACENonScSLFit,ACEpathMatrices,ACEpathLabels,4)
formatOutputMatrices(ACENonScSLFit,ACEcovMatrices,ACEcovLabels,4)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(ACE.am,ACE.cm,ACE.em,ACE.af,ACE.cf,ACE.ef), ACENonScSLFit),4)
    rownames(pathEstimatesACE) <- c('V1pathEstimates', 'V2pathEstimates', 'V3pathEstimates')
    colnames(pathEstimatesACE) <- c('am1','am2','am3','cm1','cm2','cm3','em1','em2', 'em3','af1','af2', 'af3','cf1','cf2', 'cf3', 'ef1','ef2', 'ef3')			
pathEstimatesACE

varComponentsACE <- round(mxEval(cbind(ACE.Am/ACE.Vm,ACE.Cm/ACE.Vm,ACE.Em/ACE.Vm,ACE.Af/ACE.Vf,ACE.Cf/ACE.Vf,ACE.Ef/ACE.Vf), ACENonScSLFit),2)
    rownames(varComponentsACE) <- c('V1varcomp','V2varcomp', 'V3varcomp')
    colnames(varComponentsACE) <- c('h2m1','h2m2', 'h2m3', 'c2m1','c2m2','c2m3', 'e2m1','e2m2','e2m3', 'h2f1','h2f2','h2f3', 'c2f1','c2f2', 'c2f3', 'e2f1', 'e2f2', 'e2f3')
varComponentsACE



# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
tableFitStatistics(ACENonScSLCorcFit, ACENonScSLFit)




# 4 Scalar Sex Limitation ########################################################################################################
# Quantitative sex Diff but NO Qualitative sex Diff
# In the fourth model there are male and female paths, but one set of Ra, Rc and Re between variables (same for males and females)
# Male-Female correlations in DZO group between 	A & C factors are fixed to their theoretical values of .5 and 1, respectively 
# From the previous model, we simply use the same 'labeling' to equate male and female correlations
# -------------------------------------------------------------------------------------------------------------------------------------

ACEScSLModel <- mxModel (ACENonScSLFit, name="ACEScSL",
    mxModel(ACENonScSLFit$ACE,
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("ram21", "ram31", "ram32"), lbound=-.9999, ubound=.9999, name="Ram" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("rcm21", "rcm31", "rcm32"), lbound=-.9999, ubound=.9999, name="Rcm" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("rem21", "rem31", "rem32"), lbound=-.9999, ubound=.9999, name="Rem" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("ram21", "ram31", "ram32"), lbound=-.9999, ubound=.9999, name="Raf" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("rcm21", "rcm31", "rcm32"), lbound=-.9999, ubound=.9999, name="Rcf" ),
        	mxMatrix( type="Stand", nrow=nvar, ncol=nvar, free=TRUE, values=.5, label=c("rem21", "rem31", "rem32"), lbound=-.9999, ubound=.9999, name="Ref" )
            )
)
ACEScSLFit <- mxRun(ACEScSLModel)
ACEScSLSumm <- summary(ACEScSLFit)
ACEScSLSumm


# Generate Quantitative Scalar SL ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(ACEScSLFit)
expectedMeansCovariances(ACEScSLFit)
tableFitStatistics(ACEScSLFit)
   ACEpathMatrices <- c("ACE.am","ACE.cm","ACE.em","ACE.iSDm","ACE.iSDm %*% ACE.am","ACE.iSDm %*% ACE.cm","ACE.iSDm %*% ACE.em", "ACE.af","ACE.cf","ACE.ef","ACE.iSDf","ACE.iSDf %*% ACE.af","ACE.iSDf %*% ACE.cf","ACE.iSDf %*% ACE.ef")
   ACEpathLabels <- c("pathEst_am","pathEst_cm","pathEst_em","sdm","stPathEst_am","stPathEst_cm","stPathEst_em", "pathEst_af","pathEst_cf","pathEst_ef","sdf","stPathEst_af","stPathEst_cf","stPathEst_ef")
   ACEcovMatrices <- c("ACE.Am", "ACE.Cm", "ACE.Em","ACE.Vm","ACE.Am/ACE.Vm", "ACE.Cm/ACE.Vm","ACE.Em/ACE.Vm", "ACE.Af", "ACE.Cf", "ACE.Ef","ACE.Vf","ACE.Af/ACE.Vf", "ACE.Cf/ACE.Vf","ACE.Ef/ACE.Vf")
   ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")

formatOutputMatrices(ACEScSLFit,ACEpathMatrices,ACEpathLabels,4)
formatOutputMatrices(ACEScSLFit,ACEcovMatrices,ACEcovLabels,4)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(ACE.am,ACE.cm,ACE.em,ACE.af,ACE.cf,ACE.ef), ACEScSLFit),4)
    rownames(pathEstimatesACE) <- c('V1pathEstimates', 'V2pathEstimates', 'V3pathEstimates')
    colnames(pathEstimatesACE) <- c('am1','am2','am3','cm1','cm2','cm3','em1','em2', 'em3','af1','af2', 'af3','cf1','cf2', 'cf3', 'ef1','ef2', 'ef3')				
pathEstimatesACE

varComponentsACE <- round(mxEval(cbind(ACE.Am/ACE.Vm,ACE.Cm/ACE.Vm,ACE.Em/ACE.Vm,ACE.Af/ACE.Vf,ACE.Cf/ACE.Vf,ACE.Ef/ACE.Vf), ACEScSLFit),2)
    rownames(varComponentsACE) <- c('V1varcomp','V2varcomp', 'V3varcomp')
    colnames(varComponentsACE) <- c('h2m1','h2m2', 'h2m3', 'c2m1','c2m2','c2m3', 'e2m1','e2m2','e2m3', 'h2f1','h2f2','h2f3', 'c2f1','c2f2', 'c2f3', 'e2f1', 'e2f2', 'e2f3')
varComponentsACE






# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
tableFitStatistics(ACENonScSLCorcFit, ACEScSLFit)								


# 5 Homogeneity ########################################################################################################################
# NO Quantitative dif AND NO Qualitative sex Diff
# In the fifth model there is one set of parameters, same for males and females
# Male-Female correlations in DZO group between 	A & C factors are fixed to their theoretical values of .5 and 1, respectively 
# From the previous model, we simply use the same 'labeling' to equate male and female diagonal paths
# ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

ACEHomModel <- mxModel (ACEScSLFit, name="ACEHom",
    mxModel(ACEScSLFit$ACE,
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=2, label=c("am11","am22", "am33"), name="am" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("cm11","cm22", "cm33"), name="cm" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("em11","em22", "em33"), name="em" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=2, label=c("am11","am22", "am33"), name="af" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("cm11","cm22", "cm33"), name="cf" ),
        	mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("em11","em22", "em33"), name="ef" )		 
            )
)

ACEHomFit <- mxRun(ACEHomModel)
(ACEHomSumm <- summary(ACEHomFit))


# Generate Homogeneity ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(ACEHomFit)
expectedMeansCovariances(ACEHomFit)
tableFitStatistics(ACEHomFit)


   ACEpathMatrices <- c("ACE.am","ACE.cm","ACE.em","ACE.iSDm","ACE.iSDm %*% ACE.am","ACE.iSDm %*% ACE.cm","ACE.iSDm %*% ACE.em", "ACE.af","ACE.cf","ACE.ef","ACE.iSDf","ACE.iSDf %*% ACE.af","ACE.iSDf %*% ACE.cf","ACE.iSDf %*% ACE.ef")
   ACEpathLabels <- c("pathEst_am","pathEst_cm","pathEst_em","sdm","stPathEst_am","stPathEst_cm","stPathEst_em", "pathEst_af","pathEst_cf","pathEst_ef","sdf","stPathEst_af","stPathEst_cf","stPathEst_ef")
   ACEcovMatrices <- c("ACE.Am", "ACE.Cm", "ACE.Em","ACE.Vm","ACE.Am/ACE.Vm", "ACE.Cm/ACE.Vm","ACE.Em/ACE.Vm", "ACE.Af", "ACE.Cf", "ACE.Ef","ACE.Vf","ACE.Af/ACE.Vf", "ACE.Cf/ACE.Vf","ACE.Ef/ACE.Vf")
   ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")

formatOutputMatrices(ACEHomFit,ACEpathMatrices,ACEpathLabels,4)
formatOutputMatrices(ACEHomFit,ACEcovMatrices,ACEcovLabels,4)


# Generate Table of Parameter Estimates using mxEval
pathEstimatesACE <- round(mxEval(cbind(ACE.am,ACE.cm,ACE.em,ACE.af,ACE.cf,ACE.ef), ACEHomFit),4)
    rownames(pathEstimatesACE) <- c('V1pathEstimates', 'V2pathEstimates', 'V3pathEstimates')
    colnames(pathEstimatesACE) <- c('am1','am2','am3','cm1','cm2','cm3','em1','em2', 'em3','af1','af2', 'af3','cf1','cf2', 'cf3', 'ef1','ef2', 'ef3')				
pathEstimatesACE

varComponentsACE <- round(mxEval(cbind(ACE.Am/ACE.Vm,ACE.Cm/ACE.Vm,ACE.Em/ACE.Vm,ACE.Af/ACE.Vf,ACE.Cf/ACE.Vf,ACE.Ef/ACE.Vf), ACEHomFit),2)
    rownames(varComponentsACE) <- c('V1varcomp','V2varcomp', 'V3varcomp')
    colnames(varComponentsACE) <- c('h2m1','h2m2', 'h2m3', 'c2m1','c2m2','c2m3', 'e2m1','e2m2','e2m3', 'h2f1','h2f2','h2f3', 'c2f1','c2f2', 'c2f3', 'e2f1', 'e2f2', 'e2f3')
varComponentsACE





# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
tableFitStatistics(ACEScSLFit, ACEHomFit)								


#############################################################################################################
# IP and CP MODEL

# 6 Independent Pathways Heterogeneity model ##########################################################################################################	
# The A, C and E cor between variables are 1 in this model, for both males and females
# Quantitative sex Diff: We simply specify different paths (A, C and E factor loading) for males and females
# Plus different A, C and E Specific paths for males and females
# Qualitative sex dif: Does not need testing again 
# Male-Female correlations in DZO group between 	A & C factors are fixed to their theoretical values of .5 and 1, respectively 
# ----------------------------------------------------------------------------------------------------------------------------------

nf <- 1 
ACEmvhetIPModel <- mxModel("ACEmvhetIP",
	mxModel("ACEip",
    # Matrices a, c, and e to store a, c, and e path coefficients
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("amf1","amf2", "amf3"), name="amc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("cmf1","cmf2", "cmf3"), name="cmc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("emf1","emf2", "emf3"), name="emc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("aff1","aff2", "aff3"), name="afc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("cff1","cff2", "cff3"), name="cfc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("eff1","eff2", "eff3"), name="efc" ),

        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("ams1","ams2", "ams3"), name="ams" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("cms1","cms2", "cms3"), name="cms" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("ems1","ems2", "ems3"), name="ems" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("afs1","afs2", "afs3"), name="afs" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("cfs1","cfs2", "cfs3"), name="cfs" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("efs1","efs2", "efs3"), name="efs" ),

    # Matrices A, C, and E compute variance components
        mxAlgebra( amc %*% t(amc) + ams %*% t(ams), name="Am" ),
        mxAlgebra( cmc %*% t(cmc) + cms %*% t(cms), name="Cm" ),
        mxAlgebra( emc %*% t(emc) + ems %*% t(ems), name="Em" ),
        mxAlgebra( afc %*% t(afc) + afs %*% t(afs), name="Af" ),
        mxAlgebra( cfc %*% t(cfc) + cfs %*% t(cfs), name="Cf" ),
        mxAlgebra( efc %*% t(efc) + efs %*% t(efs), name="Ef" ),


    # Algebra to compute total variances and standard deviations (diagonal only)
	mxAlgebra( Am+Cm+Em, name="Vm" ),
	mxAlgebra( Af+Cf+Ef, name="Vf" ),
	mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
	mxAlgebra( solve(sqrt(I*Vm)), name="iSDm"),
	mxAlgebra( solve(sqrt(I*Vf)), name="iSDf"),

    # Matrix & Algebra for expected means vector
        mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 18, name="Mm" ),
        mxAlgebra(  cbind(Mm,Mm), name="expMeanm"),
        mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 18, name="Mf" ),
        mxAlgebra(  cbind(Mf,Mf), name="expMeanf"),
        mxAlgebra (expression=cbind(Mm,Mf), name="expMeanOs"),

# Algebra for expected variance/covariance matrix in MZ
	  mxAlgebra(  rbind  ( cbind(Am+Cm+Em , Am+Cm),
                                        cbind(Am+Cm   , Am+Cm+Em)), name="expCovMZm" ),

	  mxAlgebra(  rbind  ( cbind(Af+Cf+Ef , Af+Cf),
                                        cbind(Af+Cf   , Af+Cf+Ef)), name="expCovMZf" ),


# Algebra for expected variance/covariance matrix in DZ

	  mxAlgebra(  rbind  ( cbind(Am+Cm+Em     , 0.5%x%Am+Cm),
                                        cbind(0.5%x%Am + Cm, Am+Cm+Em)),  name="expCovDZm" ), 
	  mxAlgebra(  rbind  ( cbind(Af+Cf+Ef    , 0.5%x%Af + Cf),
                                        cbind(0.5%x%Af + Cf , Af+Cf+Ef)),  name="expCovDZf" ),
 
# Algebra for expected variance/covariance matrix in OS

	  mxAlgebra ( rbind ( cbind(Am+Cm+Em	,	0.5%x%(amc%*%t(afc))+ 0.5%x%(ams%*%t(afs)) +  cmc%*%t(cfc)+ cms%*%t(cfs)  ),
				 cbind(  0.5%x%(afc%*%t(amc))+ 0.5%x%(afs%*%t(ams)) +  cfc%*%t(cmc)+ cfs%*%t(cms)	, 	Af+Cf+Ef) ), name="expCovOS"),


# Compute Variance Components
    	  mxAlgebra( expression=rbind(   cbind(Am/Vm,Cm/Vm,Em/Vm), cbind(Af/Vf,Cf/Vf,Ef/Vf)   ), name="VarComp") 

),		

    mxModel("MZm",
        mxData( observed=mzmData, type="raw" ),
        mxFIMLObjective( covariance="ACEip.expCovMZm", means="ACEip.expMeanm", dimnames=selVars )
    ),
    mxModel("DZm", 
        mxData( observed=dzmData, type="raw" ),
        mxFIMLObjective( covariance="ACEip.expCovDZm", means="ACEip.expMeanm", dimnames=selVars )
    ),
    mxModel("MZf",
        mxData( observed=mzfData, type="raw" ),
        mxFIMLObjective( covariance="ACEip.expCovMZf", means="ACEip.expMeanf", dimnames=selVars )
    ),
    mxModel("DZf", 
        mxData( observed=dzfData, type="raw" ),
        mxFIMLObjective( covariance="ACEip.expCovDZf", means="ACEip.expMeanf", dimnames=selVars )
    ),
    mxModel("OS", 
        mxData( observed=dzoData, type="raw" ),
        mxFIMLObjective( covariance="ACEip.expCovOS", means="ACEip.expMeanOs", dimnames=selVars ) 
    ),

    mxAlgebra(expression= MZm.objective + DZm.objective + MZf.objective + DZf.objective + OS.objective, name="modelfit" ),
    mxAlgebraObjective("modelfit")

)

   ACEmvhetIPFit <- mxRun(ACEmvhetIPModel)
   ACEmvhetIPSumm <- summary(ACEmvhetIPFit)
   ACEmvhetIPSumm


# Generate IP Heterogeneity ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(ACEmvhetIPFit)
expectedMeansCovariances(ACEmvhetIPFit)
tableFitStatistics(ACEmvhetIPFit)

   ACEpathMatricesC <- c("ACEip.amc","ACEip.cmc","ACEip.emc","ACEip.iSDm","ACEip.iSDm %*% ACEip.amc","ACEip.iSDm %*% ACEip.cmc","ACEip.iSDm %*% ACEip.emc", "ACEip.afc","ACEip.cfc","ACEip.efc","ACEip.iSDf","ACEip.iSDf %*% ACEip.afc","ACEip.iSDf %*% ACEip.cfc","ACEip.iSDf %*% ACEip.efc")
   ACEpathLabelsC <- c("pathEst_amc","pathEst_cmc","pathEst_emc","sdm","stPathEst_amc","stPathEst_cmc","stPathEst_emc", "pathEst_afc","pathEst_cfc","pathEst_efc","sdf","stPathEst_afc","stPathEst_cfc","stPathEst_efc")
   ACEpathMatricesS <- c("ACEip.ams","ACEip.cms","ACEip.ems","ACEip.iSDm","ACEip.iSDm %*% ACEip.ams","ACEip.iSDm %*% ACEip.cms","ACEip.iSDm %*% ACEip.ems", "ACEip.afs","ACEip.cfs","ACEip.efs","ACEip.iSDf","ACEip.iSDf %*% ACEip.afs","ACEip.iSDf %*% ACEip.cfs","ACEip.iSDf %*% ACEip.efs")
   ACEpathLabelsS <- c("pathEst_ams","pathEst_cms","pathEst_ems","sdm","stPathEst_ams","stPathEst_cms","stPathEst_ems", "pathEst_afs","pathEst_cfs","pathEst_efs","sdf","stPathEst_afs","stPathEst_cfs","stPathEst_efs")
   ACEcovMatrices <- c("ACEip.Am", "ACEip.Cm", "ACEip.Em","ACEip.Vm","ACEip.Am/ACEip.Vm", "ACEip.Cm/ACEip.Vm","ACEip.Em/ACEip.Vm", "ACEip.Af", "ACEip.Cf", "ACEip.Ef","ACEip.Vf","ACEip.Af/ACEip.Vf", "ACEip.Cf/ACEip.Vf","ACEip.Ef/ACEip.Vf")
   ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")


formatOutputMatrices(ACEmvhetIPFit,ACEpathMatricesC,ACEpathLabelsC,4)
formatOutputMatrices(ACEmvhetIPFit,ACEpathMatricesS,ACEpathLabelsS,4)
formatOutputMatrices(ACEmvhetIPFit,ACEcovMatrices,ACEcovLabels,4)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACEip <- round(mxEval(cbind(ACEip.amc,ACEip.cmc,ACEip.emc, ACEip.afc,ACEip.cfc,ACEip.efc),ACEmvhetIPFit),4)
    rownames(pathEstimatesACEip) <- c('V1cpath','V2cpath','V3cpath')
    colnames(pathEstimatesACEip) <- c('amc','cmc','emc', 'afc','cfc','efc')
pathEstimatesACEip
pathEstimatesACEipspecific <- round(mxEval(cbind(ACEip.ams,ACEip.cms,ACEip.ems, ACEip.afs,ACEip.cfs,ACEip.efs), ACEmvhetIPFit),4)
    rownames(pathEstimatesACEipspecific) <- c('V1spath','V2spath','V3spath')
    colnames(pathEstimatesACEipspecific) <- c('amsV1','amsV2','amsV3','cmsV1','cmsV2','cmsV3','emsV1','emsV2','emsV3', 'afsV1','afsV2','afsV3','cfsV1','cfsV2','cfsV3','efsV1','efsV2','efsV3')
pathEstimatesACEipspecific

varComponentsACEip <- round(mxEval(cbind(ACEip.Am/ACEip.Vm,ACEip.Cm/ACEip.Vm,ACEip.Em/ACEip.Vm, ACEip.Af/ACEip.Vf,ACEip.Cf/ACEip.Vf,ACEip.Ef/ACEip.Vf), ACEmvhetIPFit),4)
    rownames(varComponentsACEip) <- c('V1varcomp','V2varcomp','V3varcomp')
    colnames(varComponentsACEip) <- c('%Am1^2','%Am2^2', '%Am3^2','%Cm1^2','%Cm2^2','%Cm3^2','%Em1^2','%Em2^2','%Em3^2','%Af1^2','%Af2^2','%Af3^2','%Cf1^2','%Cf2^2','%Cf3^2','%Ef1^2', '%Ef2^2', '%Ef3^2')
varComponentsACEip




# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
tableFitStatistics(ACENonScSLFit, ACEmvhetIPFit)								


# 7 Independent Pathways Homogeneity model ##########################################################################################################	
# The A, C and E cor between variables are 1 in this model, for both males and females
# NO Quantitative sex Diff: We simply specify one set of parameters for males and females (by using same labels)
# Male-Female correlations in DZO group between 	A & C factors are fixed to their theoretical values of .5 and 1, respectively 
# ----------------------------------------------------------------------------------------------------------------------------------

ACEmvhomIPModel <- mxModel(ACEmvhetIPFit, name="ACEmvhomIP",
	mxModel(ACEmvhetIPFit$ACEip,
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("amf1","amf2", "amf3"), name="amc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("cmf1","cmf2", "cmf3"), name="cmc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("emf1","emf2", "emf3"), name="emc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("amf1","amf2", "amf3"), name="afc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("cmf1","cmf2", "cmf3"), name="cfc" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1, label=c("emf1","emf2", "emf3"), name="efc" ),

        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("ams1","ams2", "ams3"), name="ams" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("cms1","cms2", "cms3"), name="cms" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("ems1","ems2", "ems3"), name="ems" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("ams1","ams2", "ams3"), name="afs" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("cms1","cms2", "cms3"), name="cfs" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1, label=c("ems1","ems2", "ems3"), name="efs" )
	)
)

ACEmvhomIPFit <- mxRun(ACEmvhomIPModel)
ACEmvhomIPSumm <- summary(ACEmvhomIPFit)
ACEmvhomIPSumm

# Generate Homogeneity ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(ACEmvhomIPFit)
expectedMeansCovariances(ACEmvhomIPFit)
tableFitStatistics(ACEmvhomIPFit)

   ACEpathMatricesC <- c("ACEip.amc","ACEip.cmc","ACEip.emc","ACEip.iSDm","ACEip.iSDm %*% ACEip.amc","ACEip.iSDm %*% ACEip.cmc","ACEip.iSDm %*% ACEip.emc", "ACEip.afc","ACEip.cfc","ACEip.efc","ACEip.iSDf","ACEip.iSDf %*% ACEip.afc","ACEip.iSDf %*% ACEip.cfc","ACEip.iSDf %*% ACEip.efc")
   ACEpathLabelsC <- c("pathEst_amc","pathEst_cmc","pathEst_emc","sdm","stPathEst_amc","stPathEst_cmc","stPathEst_emc", "pathEst_afc","pathEst_cfc","pathEst_efc","sdf","stPathEst_afc","stPathEst_cfc","stPathEst_efc")
   ACEpathMatricesS <- c("ACEip.ams","ACEip.cms","ACEip.ems","ACEip.iSDm","ACEip.iSDm %*% ACEip.ams","ACEip.iSDm %*% ACEip.cms","ACEip.iSDm %*% ACEip.ems", "ACEip.afs","ACEip.cfs","ACEip.efs","ACEip.iSDf","ACEip.iSDf %*% ACEip.afs","ACEip.iSDf %*% ACEip.cfs","ACEip.iSDf %*% ACEip.efs")
   ACEpathLabelsS <- c("pathEst_ams","pathEst_cms","pathEst_ems","sdm","stPathEst_ams","stPathEst_cms","stPathEst_ems", "pathEst_afs","pathEst_cfs","pathEst_efs","sdf","stPathEst_afs","stPathEst_cfs","stPathEst_efs")
   ACEcovMatrices <- c("ACEip.Am", "ACEip.Cm", "ACEip.Em","ACEip.Vm","ACEip.Am/ACEip.Vm", "ACEip.Cm/ACEip.Vm","ACEip.Em/ACEip.Vm", "ACEip.Af", "ACEip.Cf", "ACEip.Ef","ACEip.Vf","ACEip.Af/ACEip.Vf", "ACEip.Cf/ACEip.Vf","ACEip.Ef/ACEip.Vf")
   ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")


formatOutputMatrices(ACEmvhomIPFit,ACEpathMatricesC,ACEpathLabelsC,4)
formatOutputMatrices(ACEmvhomIPFit,ACEpathMatricesS,ACEpathLabelsS,4)
formatOutputMatrices(ACEmvhomIPFit,ACEcovMatrices,ACEcovLabels,4)


# Generate Table of Parameter Estimates using mxEval
pathEstimatesACEip <- round(mxEval(cbind(ACEip.amc,ACEip.cmc,ACEip.emc, ACEip.afc,ACEip.cfc,ACEip.efc),ACEmvhomIPFit),4)
    rownames(pathEstimatesACEip) <- c('V1cpath','V2cpath','V3cpath')
    colnames(pathEstimatesACEip) <- c('amc','cmc','emc', 'afc','cfc','efc')
pathEstimatesACEip
pathEstimatesACEipspecific <- round(mxEval(cbind(ACEip.ams,ACEip.cms,ACEip.ems, ACEip.afs,ACEip.cfs,ACEip.efs), ACEmvhomIPFit),4)
    rownames(pathEstimatesACEipspecific) <- c('V1spath','V2spath','V3spath')
    colnames(pathEstimatesACEipspecific) <- c('amsV1','amsV2','amsV3','cmsV1','cmsV2','cmsV3','emsV1','emsV2','emsV3', 'afsV1','afsV2','afsV3','cfsV1','cfsV2','cfsV3','efsV1','efsV2','efsV3')
pathEstimatesACEipspecific

varComponentsACEip <- round(mxEval(cbind(ACEip.Am/ACEip.Vm,ACEip.Cm/ACEip.Vm,ACEip.Em/ACEip.Vm, ACEip.Af/ACEip.Vf,ACEip.Cf/ACEip.Vf,ACEip.Ef/ACEip.Vf), ACEmvhomIPFit),4)
    rownames(varComponentsACEip) <- c('V1varcomp','V2varcomp','V3varcomp')
    colnames(varComponentsACEip) <- c('%Am1^2','%Am2^2', '%Am3^2','%Cm1^2','%Cm2^2','%Cm3^2','%Em1^2','%Em2^2','%Em3^2','%Af1^2','%Af2^2','%Af3^2','%Cf1^2','%Cf2^2','%Cf3^2','%Ef1^2', '%Ef2^2', '%Ef3^2')
varComponentsACEip






# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
tableFitStatistics(ACEmvhetIPFit, ACEmvhomIPFit)	


# 8 Common Pathways model Heterogeneity model #############################################################################################################	
# There is only one set of A, C and E factors for all variables loading on one LF, for both males and females
# Quantitative sex Diff: We simply specify different paths (A, C and E factor loading) for males and females
# Plus different A, C and E Specific paths for males and females
# Qualitative sex dif: Does not need testing again 
# Male-Female correlations in DZO group between 	A & C factors are fixed to their theoretical values of .5 and 1, respectively 
# -------------------------------------------------------------------------------------------------------------------------------------


nf<-1
ACEmvhetCPModel <- mxModel("ACEmvhetCP",
    mxModel("ACEcp",
    # Matrices a, c, and e to store a, c, and e path coefficients
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1.6, label=c("amc1"), name="amc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1.6, label=c("cmc1"), name="cmc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1.6, label=c("emc1"), name="emc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1.6, label=c("afc1"), name="afc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1.6, label=c("cfc1"), name="cfc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1.6, label=c("efc1"), name="efc" ),

    # Algebra to compute variance for latent common factor and fix it to 1
        mxAlgebra(  amc %*% t(amc) + cmc %*% t(cmc) + emc %*% t(emc), name="CovarLPm" ),
        mxAlgebra(  diag2vec(CovarLPm), name="VarLPm" ),
        mxMatrix( type="Unit", nrow=nf, ncol=1, name="Unit"),
        mxConstraint (VarLPm == Unit),
	  mxAlgebra(  afc %*% t(afc) + cfc %*% t(cfc) + efc %*% t(efc), name="CovarLPf" ),
        mxAlgebra(  diag2vec(CovarLPf), name="VarLPf" ),
        mxMatrix( type="Unit", nrow=nf, ncol=1, name="Unit"),
        mxConstraint (VarLPf == Unit),

        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1.5, label=c("ldm1", "ldm2", "ldm3"), name="fm" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=1.5, label=c("ldf1", "ldf2", "ldf3"), name="ff" ),

        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1.4, label=c("ams1","ams2", "ams3"), name="ams" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1.4, label=c("cms1","cms2", "cms3"), name="cms" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1.5, label=c("ems1","ems2", "ems3"), name="ems" ),
  	  mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1.4, label=c("afs1","afs2", "afs3"), name="afs" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1.4, label=c("cfs1","cfs2", "cfs3"), name="cfs" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=1.5, label=c("efs1","efs2", "efs3"), name="efs" ),

    # Matrices A, C, and E compute variance components				
        mxAlgebra( fm %&% (amc %*% t(amc)) + ams %*% t(ams), name="Am" ),
        mxAlgebra( fm %&% (cmc %*% t(cmc)) + cms %*% t(cms), name="Cm" ),
        mxAlgebra( fm %&% (emc %*% t(emc)) + ems %*% t(ems), name="Em" ),
        mxAlgebra( ff %&% (afc %*% t(afc)) + afs %*% t(afs), name="Af" ),
        mxAlgebra( ff %&% (cfc %*% t(cfc)) + cfs %*% t(cfs), name="Cf" ),
        mxAlgebra( ff %&% (efc %*% t(efc)) + efs %*% t(efs), name="Ef" ),

    # Algebra to compute total variances and standard deviations (diagonal only)
        mxAlgebra( Am+Cm+Em, name="Vm" ),
	  mxAlgebra( Af+Cf+Ef, name="Vf" ),
        mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
        mxAlgebra( solve(sqrt(I*Vm)), name="iSDm"),
	  mxAlgebra( solve(sqrt(I*Vf)), name="iSDf"),



    # Algebra for expected means vectors
        mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= .18, label=c("Mmean1", "Mmean2", "Mmean3"),name="Mm" ),
        mxAlgebra(  cbind(Mm,Mm), name="expMeanm"),
        mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= .18, label=c("Fmean1", "Fmean2", "Fmean3"), name="Mf" ),
        mxAlgebra(  cbind(Mf,Mf), name="expMeanf"),
        mxAlgebra (expression=cbind(Mm,Mf), name="expMeanOs"),

    # Algebra for expected variance/covariance matrix in MZ
	  mxAlgebra(  rbind  ( cbind(Am+Cm+Em , Am+Cm),
                                        cbind(Am+Cm   , Am+Cm+Em)), name="expCovMZm" ),

	  mxAlgebra(  rbind  ( cbind(Af+Cf+Ef , Af+Cf),
                                        cbind(Af+Cf   , Af+Cf+Ef)), name="expCovMZf" ),


    # Algebra for expected variance/covariance matrix in DZ

	  mxAlgebra(  rbind  ( cbind(Am+Cm+Em     , 0.5%x%Am + Cm),
                                        cbind(0.5%x%Am + Cm, Am+Cm+Em)),  name="expCovDZm" ), 

	  mxAlgebra(  rbind  ( cbind(Af+Cf+Ef    , 0.5%x%Af + Cf),
                                        cbind(0.5%x%Af + Cf , Af+Cf+Ef)),  name="expCovDZf" ),
 
# Algebra for expected variance/covariance matrix in OS
	  mxAlgebra ( expression=rbind ( cbind(Am+Cm+Em	, .5%x%(fm%*%amc%*%afc%*%t(ff)) + .5%x%(ams%*%t(afs))  +  (fm%*%cmc%*%cfc%*%t(ff)) + (cms%*%t(cfs)) ),
						    cbind(   .5%x%(ff%*%amc%*%afc%*%t(fm)) + .5%x%(afs%*%t(ams))  +  (ff%*%cmc%*%cfc%*%t(fm)) + (cfs%*%t(cms))    , Af+Cf+Ef) ), name="expCovOS"),
 

	# Compute Variance Components
    	 mxAlgebra( expression=rbind(cbind(Am/Vm,Cm/Vm,Em/Vm),
    	 								 cbind(Af/Vf,Cf/Vf,Ef/Vf)), name="VarComp")
  ),
    mxModel("MZm",
        mxData( observed=mzmData, type="raw" ),
        mxFIMLObjective( covariance="ACEcp.expCovMZm", means="ACEcp.expMeanm", dimnames=selVars )
    ),
    mxModel("MZf",
        mxData( observed=mzfData, type="raw" ),
        mxFIMLObjective( covariance="ACEcp.expCovMZf", means="ACEcp.expMeanf", dimnames=selVars )

    ),
    mxModel("DZm", 
        mxData( observed=dzmData, type="raw" ),
        mxFIMLObjective( covariance="ACEcp.expCovDZm", means="ACEcp.expMeanm", dimnames=selVars )
    ),
    mxModel("DZf", 
        mxData( observed=dzfData, type="raw" ),
        mxFIMLObjective( covariance="ACEcp.expCovDZf", means="ACEcp.expMeanf", dimnames=selVars )
    ),
    mxModel("OS", 
        mxData( observed=dzoData, type="raw" ),
        mxFIMLObjective( covariance="ACEcp.expCovOS", means="ACEcp.expMeanOs", dimnames=selVars )

    ),
    mxAlgebra( expression=MZm.objective + DZm.objective + MZf.objective + DZf.objective + OS.objective, 
	name="modelfit" ),
    mxAlgebraObjective("modelfit")
# , mxCI (c("ACEcp.VarComp"))
 , mxCI (c("ACEcp.VarComp"))
)
ACEmvhetCPFit <- mxRun(ACEmvhetCPModel)	
#ACEmvhetCPFit <- mxRun(ACEmvhetCPModel, intervals=T)
ACEmvhetCPSumm <- summary(ACEmvhetCPFit)
ACEmvhetCPSumm

# Generate Heterogeneity ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(ACEmvhetCPFit)
expectedMeansCovariances(ACEmvhetCPFit)
tableFitStatistics(ACEmvhetCPFit)


	ACEpathMatricesLP <- c("ACEcp.amc","ACEcp.cmc","ACEcp.emc", "ACEcp.afc","ACEcp.cfc","ACEcp.efc")
	ACEpathLabelsLP <- c("stParEst_amc","stPathEst_cmc","stPathEst_emc", "stParEst_afc","stPathEst_cfc","stPathEst_efc")
	ACEpathMatricesS <- c("ACEcp.ams","ACEcp.cms","ACEcp.ems","ACEcp.iSDm","ACEcp.iSDm %*% ACEcp.ams","ACEcp.iSDm %*% ACEcp.cms","ACEcp.iSDm %*% ACEcp.ems")
	ACEpathLabelsS <- c("pathEst_ams","pathEst_ems","pathEst_ems","sdm","stPathEst_ams","stPathEst_cms","stPathEst_ems")
	ACEcovMatrices <- c("ACEcp.Am", "ACEcp.Cm", "ACEcp.Em","ACEcp.Vm","ACEcp.Am/ACEcp.Vm", "ACEcp.Cm/ACEcp.Vm","ACEcp.Em/ACEcp.Vm", "ACEcp.Af", "ACEcp.Cf", "ACEcp.Ef","ACEcp.Vf","ACEcp.Af/ACEcp.Vf", "ACEcp.Cf/ACEcp.Vf","ACEcp.Ef/ACEcp.Vf")
	ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")

formatOutputMatrices(ACEmvhetCPFit,ACEpathMatricesLP,ACEpathLabelsLP,4)
formatOutputMatrices(ACEmvhetCPFit,ACEpathMatricesS,ACEpathLabelsS,4)
formatOutputMatrices(ACEmvhetCPFit,ACEcovMatrices,ACEcovLabels,4)


# Generate Table of Parameter Estimates using mxEval
pathEstimatesACEcommon <- round(mxEval(cbind(ACEcp.amc,ACEcp.cmc,ACEcp.emc, ACEcp.afc,ACEcp.cfc,ACEcp.efc),ACEmvhetCPFit),4)
    rownames(pathEstimatesACEcommon)  <- c('pathestimates')
    colnames(pathEstimatesACEcommon) <- c('amc','cmc','emc', 'afc','cfc','efc')
pathEstimatesACEcommon
pathEstimatesACEspecific <- round(mxEval(cbind(ACEcp.ams,ACEcp.cms,ACEcp.ems, ACEcp.afs,ACEcp.cfs,ACEcp.efs), ACEmvhetCPFit),4)
    rownames(pathEstimatesACEspecific)  <- c('V1Spath','V2cpath','V3Spath')
    colnames(pathEstimatesACEspecific) <- c('amsV1','amsV2','amsV3','cmsV1','cmsV2','cmsV3','emsV1','emsV2','emsV3', 'afsV1','afsV2','afsV3','cfsV1','cfsV2','cfsV3','efsV1','efsV2','efsV3')
pathEstimatesACEspecific

varComponentsACE <- round(mxEval(cbind(ACEcp.Am/ACEcp.Vm,ACEcp.Cm/ACEcp.Vm,ACEcp.Em/ACEcp.Vm, ACEcp.Af/ACEcp.Vf,ACEcp.Cf/ACEcp.Vf,ACEcp.Ef/ACEcp.Vf), ACEmvhetCPFit),4)
    rownames(varComponentsACE) <- c('V1varcomp','V2varcomp','V3varcomp')
    colnames(varComponentsACE) <- c('%Am1^2','%Am2^2', '%Am3^2','%Cm1^2','%Cm2^2','%Cm3^2','%Em1^2','%Em2^2','%Em3^2','%Af1^2','%Af2^2','%Af3^2','%Cf1^2','%Cf2^2','%Cf3^2','%Ef1^2', '%Ef2^2', '%Ef3^2')
varComponentsACE





# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
tableFitStatistics(ACENonScSLFit, ACEmvhetCPFit)	


# 9 Common Pathways model Homogeneity model #############################################################################################################	
# There is only one set of A, C and E factors for all variables loading on one LF, for both males and females
# NO Quantitative sex Diff: We simply specify one set of parameters for males and females (using same labels)
# Male-Female correlations in DZO group between 	A & C factors are fixed to their theoretical values of .5 and 1, respectively 
# -------------------------------------------------------------------------------------------------------------------------------------


ACEmvhomCPModel <- mxModel(ACEmvhetCPFit, name="ACEmvhomCP",
	mxModel(ACEmvhetCPFit$ACEcp,
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1, label=c("amc1"), name="amc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1, label=c("cmc1"), name="cmc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1, label=c("emc1"), name="emc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1, label=c("amc1"), name="afc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1, label=c("cmc1"), name="cfc" ),
        mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=1, label=c("emc1"), name="efc" ),

	# Algebra for expected means vectors
        #mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= .18, label=c("Mmean1", "Mmean2", "Mmean3"),name="Mm" ),
        #mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= .18, label=c("Mmean1", "Mmean2", "Mmean3"), name="Mf" ),
        

        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=.9, label=c("ldm1", "ldm2", "ldm3"), name="fm" ),
        mxMatrix( type="Full", nrow=nvar, ncol=nf, free=TRUE, values=.9, label=c("ldm1", "ldm2", "ldm3"), name="ff" ),

        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=.8, label=c("ams1","ams2", "ams3"), name="ams" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=.8, label=c("cms1","cms2", "cms3"), name="cms" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=.9, label=c("ems1","ems2", "ems3"), name="ems" ),
  	  mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=.8, label=c("ams1","ams2", "ams3"), name="afs" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=.8, label=c("cms1","cms2", "cms3"), name="cfs" ),
        mxMatrix( type="Diag", nrow=nvar, ncol=nvar, free=TRUE, values=.9, label=c("ems1","ems2", "ems3"), name="efs" )
	  )
	
)

ACEmvhomCPFit <- mxRun(ACEmvhomCPModel)
ACEmvhomCPSumm <- summary(ACEmvhomCPFit)
ACEmvhomCPSumm

# Generate Homogeneity CP ACE Output
# -----------------------------------------------------------------------
parameterSpecifications(ACEmvhomCPFit)
expectedMeansCovariances(ACEmvhomCPFit)
tableFitStatistics(ACEmvhomCPFit)
	ACEpathMatricesLP <- c("ACEcp.amc","ACEcp.cmc","ACEcp.emc", "ACEcp.afc","ACEcp.cfc","ACEcp.efc")
	ACEpathLabelsLP <- c("stParEst_amc","stPathEst_cmc","stPathEst_emc", "stParEst_afc","stPathEst_cfc","stPathEst_efc")
	ACEpathMatricesS <- c("ACEcp.ams","ACEcp.cms","ACEcp.ems","ACEcp.iSDm","ACEcp.iSDm %*% ACEcp.ams","ACEcp.iSDm %*% ACEcp.cms","ACEcp.iSDm %*% ACEcp.ems", "ACEcp.afs","ACEcp.cfs","ACEcp.efs","ACEcp.iSDf","ACEcp.iSDf %*% ACEcp.afs","ACEcp.iSDf %*% ACEcp.cfs","ACEcp.iSDf %*% ACEcp.efs")
	ACEpathLabelsS <- c("pathEst_ams","pathEst_cms","pathEst_ems","sdm","stPathEst_ams","stPathEst_cms","stPathEst_ems", "pathEst_afs","pathEst_cfs","pathEst_efs","sdf","stPathEst_afs","stPathEst_cfs","stPathEst_efs")
	ACEcovMatrices <- c("ACEcp.Am", "ACEcp.Cm", "ACEcp.Em","ACEcp.Vm","ACEcp.Am/ACEcp.Vm", "ACEcp.Cm/ACEcp.Vm","ACEcp.Em/ACEcp.Vm", "ACEcp.Af", "ACEcp.Cf", "ACEcp.Ef","ACEcp.Vf","ACEcp.Af/ACEcp.Vf", "ACEcp.Cf/ACEcp.Vf","ACEcp.Ef/ACEcp.Vf")
	ACEcovLabels <- c("covComp_Am","covComp_Cm","covComp_Em","Varm","stCovComp_Am","stCovComp_Cm","stCovComp_Em", "covComp_Af","covComp_Cf","covComp_Ef","Varf","stCovComp_Af","stCovComp_Cf","stCovComp_Ef")

formatOutputMatrices(ACEmvhomCPFit,ACEpathMatricesLP,ACEpathLabelsLP,4)
formatOutputMatrices(ACEmvhomCPFit,ACEpathMatricesS,ACEpathLabelsS,4)
formatOutputMatrices(ACEmvhomCPFit,ACEcovMatrices,ACEcovLabels,4)

# Generate Table of Parameter Estimates using mxEval
pathEstimatesACEcommon <- round(mxEval(cbind(ACEcp.amc,ACEcp.cmc,ACEcp.emc, ACEcp.afc,ACEcp.cfc,ACEcp.efc),ACEmvhomCPFit),4)
    rownames(pathEstimatesACEcommon)  <- c('pathestimates')
    colnames(pathEstimatesACEcommon) <- c('amc','cmc','emc', 'afc','cfc','efc')
pathEstimatesACEcommon
pathEstimatesACEspecific <- round(mxEval(cbind(ACEcp.ams,ACEcp.cms,ACEcp.ems, ACEcp.afs,ACEcp.cfs,ACEcp.efs), ACEmvhomCPFit),4)
    rownames(pathEstimatesACEspecific)  <- c('V1spathest','V2spathest','V3spathest')
    colnames(pathEstimatesACEspecific) <- c('amsV1','amsV2','amsV3','cmsV1','cmsV2','cmsV3','emsV1','emsV2','emsV3', 'afsV1','afsV2','afsV3','cfsV1','cfsV2','cfsV3','efsV1','efsV2','efsV3')
pathEstimatesACEspecific

varComponentsACE <- round(mxEval(cbind(ACEcp.Am/ACEcp.Vm,ACEcp.Cm/ACEcp.Vm,ACEcp.Em/ACEcp.Vm, ACEcp.Af/ACEcp.Vf,ACEcp.Cf/ACEcp.Vf,ACEcp.Ef/ACEcp.Vf), ACEmvhomCPFit),4)
    rownames(varComponentsACE) <- c('V1varcomp','V2varcomp','V3varcomp')
    colnames(varComponentsACE) <- c('%Am1^2','%Am2^2', '%Am3^2','%Cm1^2','%Cm2^2','%Cm3^2','%Em1^2','%Em2^2','%Em3^2','%Af1^2','%Af2^2','%Af3^2','%Cf1^2','%Cf2^2','%Cf3^2','%Ef1^2', '%Ef2^2', '%Ef3^2')
varComponentsACE






# Print Comparative Fit Statistics
# -----------------------------------------------------------------------
tableFitStatistics(ACEmvhetCPFit, ACEmvhomCPFit)	


# Print Comparative Fit Statistics

compare <- c(ACENonScSLCorcFit, ACENonScSLFit,  ACEScSLFit, ACEHomFit, ACEmvhetIPFit, ACEmvhomIPFit, ACEmvhetCPFit, ACEmvhomCPFit)
tableFitStatistics(ACENonScSLCoraFit, compare)	
compare <- c(ACEHomFit,  ACEmvhomIPFit,  ACEmvhomCPFit)
tableFitStatistics(ACEHomFit, compare)	



# Generates an output table of all submodels compared to previous (Nested models)
# --------------------------------------------------------------------------------------------------------------------------------

Nested.fit <- 	rbind(	mxCompare(ACENonScSLCoraFit, ACENonScSLFit),
					mxCompare(ACENonScSLCorcFit, ACENonScSLFit)[2,],
	   				mxCompare(ACENonScSLFit, ACEHomFit)[2,],
	  				mxCompare(ACENonScSLFit, ACEmvhetIPFit)[2,],
	  				mxCompare(ACEmvhetIPFit, ACEmvhomIPFit)[2,],
	  				mxCompare(ACENonScSLFit, ACEmvhetCPFit)[2,],
	  				mxCompare(ACEmvhetCPFit, ACEmvhomCPFit)[2,])#, #[,c(-4,-6)]
Nested.fit

write.csv(Nested.fit,"Model.Fitting.csv",row.names=FALSE)