# ------------------------------------------------------------------------------ # Program: twinMulAceCon.R # Author: Hermine Maes # Date: 03 04 2014 # # Multivariate Twin Analysis model to estimate causes of variation # Matrix style model - Raw data - Continuous data # # Modified: Eero Vuoksimaa 5/22/2015 # mxExpectationNormal used instead of mxFIMLObjective which was used in older versions of OpenMx # CE Model added # constraining of genetic and environmental correlations to zero added # -------|---------|---------|---------|---------|---------|---------| setwd("Z:Documents/1_OpenMx_Finland") # Load Library require(OpenMx) require(psych) source("myFunctions.R") source("GenEpiHelperFunctions.R") # -------------------------------------------------------------------- # PREPARE DATA # Load Data fsiq <- read.table ("FT12_iq_z.txt", header=T, na=-99) describe(fsiq, skew=F) # check the file dim(fsiq)[1] #number of cases dim(fsiq)[2] #number of variables names(fsiq)[1:2] #zygocity and typing, values are are constant within family names(fsiq)[3:10] #phenotypes that can have different values within family # Select Variables for Analysis Vars <- c('zresheight','zresweight') #Vars <- c('zresheight','zresfsiq') #Vars <- c('zresweight','zresfsiq') nv <- 2 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Data for Analysis mzmData <- subset(fsiq, typing==1, selVars) mzfData <- subset(fsiq, typing==2, selVars) dzmData <- subset(fsiq, typing==3, selVars) dzfData <- subset(fsiq, typing==4, selVars) osdzData <- subset(fsiq, typing==5, selVars) describe(mzmData, skew=F) describe(mzfData, skew=F) describe(dzmData, skew=F) describe(dzfData, skew=F) describe(osdzData, skew=F) dim(mzmData) dim(mzfData) dim(dzmData) dim(dzfData) dim(osdzData) # Generate Descriptive Statistics colMeans(mzmData,na.rm=TRUE) colMeans(mzfData,na.rm=TRUE) colMeans(dzmData,na.rm=TRUE) colMeans(dzfData,na.rm=TRUE) colMeans(osdzData,na.rm=TRUE) cov(mzmData) cov(mzfData,use="complete") cov(dzmData,use="complete") cov(dzfData,use="complete") cov(osdzData,use="complete") cor(mzmData,use="complete") cor(mzfData,use="complete") cor(dzmData,use="complete") cor(dzfData,use="complete") cor(osdzData,use="complete") # ------------------------------------------------------------------------------ # PREPARE SATURATED MODEL # Saturated Model # Set Starting Values svMe <- c(0,0) # start value for means svVa <- valDiag(ntv,0.3) # start values for variances lbVa <- valODiag(ntv,.0001,-10) # lower bounds for covariances # -------|---------|---------|---------|---------|---------|---------| ## Algebra for expected Mean Matrices in MZ & DZ twins meanMZM <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labFull("meMZM",1,ntv), name="expMeanMZM" ) meanMZF <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labFull("meMZF",1,ntv), name="expMeanMZF" ) meanDZM <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labFull("meDZM",1,ntv), name="expMeanDZM" ) meanDZF <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labFull("meDZF",1,ntv), name="expMeanDZF" ) meanOSDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labFull("meOSDZ",1,ntv), name="expMeanOSDZ" ) ## Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covMZM <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVa, lbound=lbVa, labels=labSymm("vaMZM",ntv), name="expCovMZM" ) covMZF <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVa, lbound=lbVa, labels=labSymm("vaMZF",ntv), name="expCovMZF" ) covDZM <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVa, lbound=lbVa, labels=labSymm("vaDZM",ntv), name="expCovDZM" ) covDZF <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVa, lbound=lbVa, labels=labSymm("vaDZF",ntv), name="expCovDZF" ) covOSDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVa, lbound=lbVa, labels=labSymm("vaOSDZ",ntv), name="expCovOSDZ" ) # Data objects for Multiple Groups dataMZM <- mxData( observed=mzmData, type="raw" ) dataMZF <- mxData( observed=mzfData, type="raw" ) dataDZM <- mxData( observed=dzmData, type="raw" ) dataDZF <- mxData( observed=dzfData, type="raw" ) dataOSDZ <- mxData( observed=osdzData, type="raw" ) objMZM <- mxExpectationNormal( covariance="expCovMZM", means="expMeanMZM", dimnames=selVars ) objMZF <- mxExpectationNormal( covariance="expCovMZF", means="expMeanMZF", dimnames=selVars ) objDZM <- mxExpectationNormal( covariance="expCovDZM", means="expMeanDZM", dimnames=selVars ) objDZF <- mxExpectationNormal( covariance="expCovDZF", means="expMeanDZF", dimnames=selVars ) objOSDZ <- mxExpectationNormal( covariance="expCovOSDZ", means="expMeanOSDZ", dimnames=selVars ) fitFunction <- mxFitFunctionML() # Combine Groups modelMZM <- mxModel( "MZM", meanMZM, covMZM, dataMZM, objMZM, fitFunction ) modelMZF <- mxModel( "MZF", meanMZF, covMZF, dataMZF, objMZF, fitFunction ) modelDZM <- mxModel( "DZM", meanDZM, covDZM, dataDZM, objDZM, fitFunction ) modelDZF <- mxModel( "DZF", meanDZF, covDZF, dataDZF, objDZF, fitFunction ) modelOSDZ <- mxModel( "OSDZ", meanOSDZ, covOSDZ, dataOSDZ, objOSDZ, fitFunction ) minus2ll <- mxAlgebra( MZM.objective+ MZF.objective+ DZM.objective+ DZF.objective+ OSDZ.objective, name="minus2sumloglikelihood" ) obj <- mxFitFunctionAlgebra( "minus2sumloglikelihood" ) ciCov <- mxCI( c('MZM.expCovMZM','MZF.expCovMZF','DZM.expCovDZM','DZF.expCovDZF','OSDZ.expCovOSDZ' )) ciMean <- mxCI( c('MZM.expMeanMZM','MZF.expMeanMZF','DZM.expMeanDZM','DZF.expMeanDZF','OSDZ.expMeanOSDZ' )) twinSatModel <- mxModel( "twinSat", modelMZM, modelMZF, modelDZM, modelDZF, modelOSDZ, minus2ll, obj, ciCov, ciMean ) # ------------------------------------------------------------------------------ # RUN SATURATED MODEL # Run Saturated Model #twinSatFit <- mxRun( twinSatModel, intervals=F ) #twinSatSumm <- summary( twinSatFit) # Generate Saturated Model Output #round(twinSatFit$MZM.expMeanMZM@values,2) #round(twinSatFit$MZF.expMeanMZF@values,2) #round(twinSatFit$DZM.expMeanDZM@values,2) #round(twinSatFit$DZF.expMeanDZF@values,2) #round(twinSatFit$OSDZ.expMeanOSDZ@values,2) #round(twinSatFit$MZM.expCovMZM@values,2) #round(twinSatFit$MZF.expCovMZF@values,2) #round(twinSatFit$DZM.expCovDZM@values,2) #round(twinSatFit$DZF.expCovDZF@values,2) #round(twinSatFit$OSDZ.expCovOSDZ@values,2) # ------------------------------------------------------------------------------ # RUN SATURATED SUBMODELS # Constrain expected Means and Variances to be equal across twin order #eqMeVaTwinModel <- mxModel( twinSatModel, name="eqM&Vtwins" ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labFull("meMZM",1,ntv), # free=TRUE, values=svMe, newlabels=labFull("meMZM",1,nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labFull("meMZF",1,ntv), # free=TRUE, values=svMe, newlabels=labFull("meMZF",1,nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labFull("meDZM",1,ntv), # free=TRUE, values=svMe, newlabels=labFull("meDZM",1,nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labFull("meDZF",1,ntv), # free=TRUE, values=svMe, newlabels=labFull("meDZF",1,nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labFull("meOSDZ",1,ntv), # free=TRUE, values=svMe, newlabels=labFull("meOSDZ",1,nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labDiag("vaMZM",ntv), # free=TRUE, values=diag(svVa), newlabels=labDiag("vaMZM",nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labDiag("vaMZF",ntv), # free=TRUE, values=diag(svVa), newlabels=labDiag("vaMZF",nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labDiag("vaDZM",ntv), # free=TRUE, values=diag(svVa), newlabels=labDiag("vaDZM",nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labDiag("vaDZF",ntv), # free=TRUE, values=diag(svVa), newlabels=labDiag("vaDZF",nv) ) #eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=labDiag("vaOSDZ",ntv), # free=TRUE, values=diag(svVa), newlabels=labDiag("vaOSDZ",nv) ) # #eqMeVaTwinFit <- mxRun( eqMeVaTwinModel, intervals=F ) #eqMeVaTwinSumm <- summary( eqMeVaTwinFit ) #mxCompare(twinSatFit, eqMeVaTwinFit) # Constrain expected Means and Variances to be equal across twin order and zygosity #eqMeVaZygModel <- mxModel( eqMeVaTwinModel, name="eqM&Vzyg" ) #eqMeVaZygModel <- omxSetParameters( eqMeVaZygModel, labels=c(labFull("meMZM",1,nv),labFull("meDZM",1,nv),labFull("meMZF",1,nv),labFull("meDZF",1,nv),labFull("meOSDZ",1,nv)), # free=TRUE, values=svMe, newlabels=labFull("me",1,nv) ) #eqMeVaZygModel <- omxSetParameters( eqMeVaZygModel, labels=c(labDiag("vaMZM",nv),labDiag("vaDZM",nv),labDiag("vaMZF",nv),labDiag("vaDZF",nv),labDiag("vaOSDZ",nv)), # free=TRUE, values=diag(svVa), newlabels=labDiag("va",nv) ) #eqMeVaZygFit <- mxRun( eqMeVaZygModel, intervals=F ) #eqMeVaZygSumm <- summary( eqMeVaZygFit ) #mxCompare(eqMeVaTwinFit, eqMeVaZygFit) # Print Comparative Fit Statistics #SatNested <- list(eqMeVaTwinFit, eqMeVaZygFit) #mxCompare(twinSatFit, SatNested) # ------------------------------------------------------------------------------ # PREPARE GENETIC MODEL # ------------------------------------------------------------------------------ # Cholesky Decomposition ACE Model # ------------------------------------------------------------------------------ # Set Starting Values svMe <- c(0,0) # start value for means svPa <- valDiag(nv,.3) # start values for parameters on diagonal lbPa <- valLUDiag(nv,.0001,-10,NA) # lower bounds for parameters on diagonal # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, labels=labLower("a",nv), lbound=lbPa, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, labels=labLower("c",nv), lbound=lbPa, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, labels=labLower("e",nv), lbound=lbPa, name="e" ) # Matrices generated to hold A, C, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) covP <- mxAlgebra( expression=A+C+E, name="V" ) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labFull("me",1,nv), name="expMean" ) covMZM <- mxAlgebra( expression= rbind( cbind(V, A+C), cbind(A+C, V)), name="expCovMZM" ) covMZF <- mxAlgebra( expression= rbind( cbind(V, A+C), cbind(A+C, V)), name="expCovMZF" ) covDZM <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovDZM" ) covDZF <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovDZF" ) covOSDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovOSDZ" ) # Data objects for Multiple Groups dataMZM <- mxData( observed=mzmData, type="raw" ) dataMZF <- mxData( observed=mzfData, type="raw" ) dataDZM <- mxData( observed=dzmData, type="raw" ) dataDZF <- mxData( observed=dzfData, type="raw" ) dataOSDZ <- mxData( observed=osdzData, type="raw" ) objMZM <- mxExpectationNormal( covariance="expCovMZM", means="expMean", dimnames=selVars ) objMZF <- mxExpectationNormal( covariance="expCovMZF", means="expMean", dimnames=selVars ) objDZM <- mxExpectationNormal( covariance="expCovDZM", means="expMean", dimnames=selVars ) objDZF <- mxExpectationNormal( covariance="expCovDZF", means="expMean", dimnames=selVars ) objOSDZ <- mxExpectationNormal( covariance="expCovOSDZ", means="expMean", dimnames=selVars ) fitFunction <- mxFitFunctionML() # Combine Groups pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, meanG ) modelMZM <- mxModel( pars, covMZM, dataMZM, objMZM, fitFunction, name="MZM" ) modelMZF <- mxModel( pars, covMZF, dataMZF, objMZF, fitFunction, name="MZF" ) modelDZM <- mxModel( pars, covDZM, dataDZM, objDZM, fitFunction, name="DZM" ) modelDZF <- mxModel( pars, covDZF, dataDZF, objDZF, fitFunction, name="DZF" ) modelOSDZ <- mxModel( pars, covOSDZ, dataOSDZ, objOSDZ, fitFunction, name="OSDZ" ) minus2ll <- mxAlgebra( expression=MZM.objective + MZF.objective + DZM.objective + DZF.objective + OSDZ.objective , name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) CholAceModel <- mxModel( "CholACE", pars, modelMZM, modelMZF, modelDZM, modelDZF, modelOSDZ, minus2ll, obj ) # ------------------------------------------------------------------------------ # RUN GENETIC MODEL # Run Cholesky Decomposition ACE model CholAceFit <- mxRun(CholAceModel) #mxCompare(twinSatFit,CholAceFit) #mxCompare(eqMeVaZygFit,CholAceFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices # ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) CholACEpathMatrices <- c("iSD %*% a","iSD %*% c","iSD %*% e") CholACEpathLabels <- c("stPathA","stPathC","stPathE") formatOutputMatrices(CholAceFit,CholACEpathMatrices,CholACEpathLabels,Vars,4) # ACE Covariance Matrices & Proportions of Variance Matrices ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE") formatOutputMatrices(CholAceFit,ACEcovMatrices,ACEcovLabels,Vars,4) # ACE Correlation Matrices ACEcorMatrices <- c("solve(sqrt(I*A)) %&% A","solve(sqrt(I*C)) %&% C","solve(sqrt(I*E)) %&% E") ACEcorLabels <- c("corA","corC","corE") formatOutputMatrices(CholAceFit,ACEcorMatrices,ACEcorLabels,Vars,4) # ------------------------------------------------------------------------------ # Cholesky Decomposition AE Model # ------------------------------------------------------------------------------ CholAeModel <- mxModel( CholAceFit, name="CholAE" ) CholAeModel <- omxSetParameters( CholAeModel, labels=labLower("c",nv), free=FALSE, values=0 ) CholAeFit <- mxRun(CholAeModel) mxCompare(CholAceFit, CholAeFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices # ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) CholAEpathMatrices <- c("iSD %*% a","iSD %*% c","iSD %*% e") CholAEpathLabels <- c("stPathA","stPathC","stPathE") formatOutputMatrices(CholAeFit,CholAEpathMatrices,CholAEpathLabels,Vars,4) # ACE Covariance Matrices & Proportions of Variance Matrices AEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") AEcovLabels <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE") formatOutputMatrices(CholAeFit,AEcovMatrices,AEcovLabels,Vars,4) # ACE Correlation Matrices AEcorMatrices <- c("solve(sqrt(I*A)) %&% A","solve(sqrt(I*E)) %&% E") AEcorLabels <- c("corA","corE") formatOutputMatrices(CholAeFit,AEcorMatrices,AEcorLabels,Vars,4) # ------------------------------------------------------------------------------ # Cholesky Decomposition CE Model # ------------------------------------------------------------------------------ CholAcModel <- mxModel( CholAceFit, name="CholAC" ) CholAcModel <- omxSetParameters( CholAcModel, labels=labLower("a",nv), free=FALSE, values=0 ) CholAcFit <- mxRun(CholAcModel) mxCompare(CholAceFit, CholAcFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices # ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) CholACpathMatrices <- c("iSD %*% a","iSD %*% c","iSD %*% e") CholACpathLabels <- c("stPathA","stPathC","stPathE") formatOutputMatrices(CholAcFit,CholACpathMatrices,CholACpathLabels,Vars,4) # ACE Covariance Matrices & Proportions of Variance Matrices ACcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACcovLabels <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE") formatOutputMatrices(CholAcFit,ACcovMatrices,ACcovLabels,Vars,4) # ACE Correlation Matrices ACcorMatrices <- c("solve(sqrt(I*C)) %&% C","solve(sqrt(I*E)) %&% E") ACcorLabels <- c("corC","corE") formatOutputMatrices(CholAcFit,ACcorMatrices,ACcorLabels,Vars,4) # ------------------------------------------------------------------------------ # set cor_A = 0 # ------------------------------------------------------------------------------ CholAeModel_noAcor <- mxModel( CholAeFit, name="CholAE_no_a_cor" ) CholAeModel_noAcor <- omxSetParameters( CholAeModel_noAcor, labels=labLower("a",nv), free=c(TRUE,FALSE,TRUE), values=c(.5,0,.5)) CholAe_noAcorFit <- mxRun(CholAeModel_noAcor) mxCompare(CholAeFit, CholAe_noAcorFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices # ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) CholAE_noAcorpathMatrices <- c("iSD %*% a","iSD %*% c","iSD %*% e") CholAE_noAcorpathLabels <- c("stPathA","stPathC","stPathE") formatOutputMatrices(CholAe_noAcorFit,CholAE_noAcorpathMatrices,CholAE_noAcorpathLabels,Vars,4) # ACE Covariance Matrices & Proportions of Variance Matrices AE_noAcorcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") AE_noAcorcovLabels <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE") formatOutputMatrices(CholAe_noAcorFit,AE_noAcorcovMatrices,AE_noAcorcovLabels,Vars,4) # ACE Correlation Matrices AE_noAcorcorMatrices <- c("solve(sqrt(I*A)) %&% A","solve(sqrt(I*E)) %&% E") AE_noAcorcorLabels <- c("corA","corE") formatOutputMatrices(CholAe_noAcorFit,AE_noAcorcorMatrices,AE_noAcorcorLabels,Vars,4) # ------------------------------------------------------------------------------ # set cor_E = 0 # ------------------------------------------------------------------------------ CholAeModel_noEcor <- mxModel( CholAeFit, name="CholAE_no_e_cor" ) CholAeModel_noEcor <- omxSetParameters( CholAeModel_noEcor, labels=labLower("e",nv), free=c(TRUE,FALSE,TRUE), values=c(.5,0,.5)) CholAe_noEcorFit <- mxRun(CholAeModel_noEcor) mxCompare(CholAeFit, CholAe_noEcorFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices # ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) CholAE_noEcorpathMatrices <- c("iSD %*% a","iSD %*% c","iSD %*% e") CholAE_noEcorpathLabels <- c("stPathA","stPathC","stPathE") formatOutputMatrices(CholAe_noEcorFit,CholAE_noEcorpathMatrices,CholAE_noEcorpathLabels,Vars,4) # ACE Covariance Matrices & Proportions of Variance Matrices AE_noEcorcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") AE_noEcorcovLabels <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE") formatOutputMatrices(CholAe_noEcorFit,AE_noEcorcovMatrices,AE_noEcorcovLabels,Vars,4) # ACE Correlation Matrices AE_noEcorcorMatrices <- c("solve(sqrt(I*A)) %&% A","solve(sqrt(I*E)) %&% E") AE_noEcorcorLabels <- c("corA","corE") formatOutputMatrices(CholAe_noEcorFit,AE_noEcorcorMatrices,AE_noEcorcorLabels,Vars,4) # ------------------------------------------------------------------------------ # set cor_A = 0 & cor_E = 0 # ------------------------------------------------------------------------------ CholAeModel_noAEcor <- mxModel( CholAeFit, name="CholAE_no_ae_cor" ) CholAeModel_noAEcor <- omxSetParameters( CholAeModel_noAEcor, labels=labLower("a",nv), free=c(TRUE,FALSE,TRUE), values=c(.5,0,.5)) CholAeModel_noAEcor <- omxSetParameters( CholAeModel_noAEcor, labels=labLower("e",nv), free=c(TRUE,FALSE,TRUE), values=c(.5,0,.5)) CholAe_noAEcorFit <- mxRun(CholAeModel_noAEcor) mxCompare(CholAeFit, CholAe_noAEcorFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices # ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) CholAE_noAEcorpathMatrices <- c("iSD %*% a","iSD %*% c","iSD %*% e") CholAE_noAEcorpathLabels <- c("stPathA","stPathC","stPathE") formatOutputMatrices(CholAe_noAEcorFit,CholAE_noAEcorpathMatrices,CholAE_noAEcorpathLabels,Vars,4) # ACE Covariance Matrices & Proportions of Variance Matrices AE_noAEcorcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") AE_noAEcorcovLabels <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE") formatOutputMatrices(CholAe_noAEcorFit,AE_noAEcorcovMatrices,AE_noAEcorcovLabels,Vars,4) # ACE Correlation Matrices AE_noAEcorcorMatrices <- c("solve(sqrt(I*A)) %&% A","solve(sqrt(I*E)) %&% E") AE_noAEcorcorLabels <- c("corA","corE") formatOutputMatrices(CholAe_noAEcorFit,AE_noAEcorcorMatrices,AE_noAEcorcorLabels,Vars,4)