# ----------------------------------------------------------------------- # 5/16/2011 Changing values for thresholds on remission and changing matrix and # algebra in ACE model to accomodate separate male and female values # # TESTING WIDER STARTING VALUES TO SEE IF CAN GET RID OF CODE RED MESSAGE - 5/16/2011 OK # # TESTING DIFFERENCES BETWEEN GROUPS (5-groups: MZF, MZM, DZF, DZM, DZOS) # 5/6/2011 THIS WORKS WITH MOST RECENT R (2.13.0) and OPENMX (1.0.6-1581)BUILDS # # 4/28/2010 trying to do this in R - Cholesky w AUD and remission # Program: OrdinalMultivariateTwinAnalysis_MatrixRaw.R # Author: Sarah Medland (Adapted from Meike Bartels script which was based on script H.M.) # Date: 01-03-2010 # # Multivariate Twin Saturated model to estimate means and (co)variances across multiple groups # Multivariate Cholesky ACE model to estimate genetic and environmental sources of variance # Matrix style model input - Raw data input # # NOTE- SCRIPT NOT GENERALIZABLE. MUST REWRITE MATRIX AND ALGEBRA STATEMENTS FOR # ORDINAL VARIABLES WITH DIFFERENT NUMBERS OF LEVELS # ----------------------------------------------------------------------- setwd("S:/My Documents/AUSTRA/remission_relapse/remis_2011") require(OpenMx) source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R") # Prepare Data # Read data from csv file audrem.csv # Variables zyg89 female preg thrsh_a rem_a thrsh_b rem_b # ----------------------------------------------------------------------- audrem <-read.table ('audrem.csv', sep=',', na='', header=TRUE) #create subset with variables for analysis nthresh1 <- 4 nthresh2 <- 1 Naudrem <-data.frame(audrem$zyg89, audrem$thrsh_a, audrem$thrsh_b, audrem$rem_a, audrem$rem_b) names(Naudrem) <- c("zyg89", "thrsh1", "thrsh2", "rem1", "rem2") Naudrem[,c(2,3)] <- mxFactor(Naudrem[,c(2,3)], c(0 : nthresh1)) Naudrem[,c(4,5)] <- mxFactor(Naudrem[,c(4,5)], c(0 : nthresh2)) # tell R how many levels in the ordinal variables mxFactor(Naudrem$thrsh[,1:2], levels=c(0:4), ordered=TRUE) mxFactor(Naudrem$rem[,1:2], levels=c(0:1), ordered=TRUE) head(Naudrem) summary(Naudrem) str(Naudrem) table(Naudrem$zyg89) table(Naudrem$thrsh1) table(Naudrem$thrsh2) table(Naudrem$rem1) table(Naudrem$rem2) nv <- 2 ntv <- nv*2 maxthresh <- 4 selVars <- c('thrsh1', 'rem1', 'thrsh2', 'rem2') # name thresholds for dimname statement (to label rows and columns) thresh <-c("Th1", "Th2", "Th3", "Th4") mzfData <- subset(Naudrem, zyg89==1, selVars) dzfData <- subset(Naudrem, zyg89==3, selVars) mzmData <- subset(Naudrem, zyg89==2, selVars) dzmData <- subset(Naudrem, zyg89==4, selVars) dzosData <- subset(Naudrem, zyg89==5, selVars) summary(mzfData) summary(dzfData) summary(mzmData) summary(dzmData) summary(dzosData) table(mzfData$thrsh1) table(dzfData$thrsh1) table(mzmData$thrsh1) table(dzmData$thrsh1) table(dzosData$thrsh1) # Fit Multivariate Saturated Model # ----------------------------------------------------------------------- multTwinSatModel <- mxModel("multTwinSat", mxModel("MZF", mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="expCovMZF" ), mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, values=0, name="expMeanMZF" ), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=TRUE, values=c(-1.0, .01, .10, .20), lbound=c(-2.0, .00001, .01, .10), name="V1T"), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=c(TRUE,FALSE,FALSE,FALSE), values=c(-.10,NA,NA,NA), name="V2T"), # NOTE Inc means increments between theresholds mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Inc" ), mxAlgebra( expression= cbind((Inc %*% V1T),V2T,(Inc %*% V1T),V2T), dimnames=list(thresh,selVars), name="expThMZF"), mxData( observed=mzfData, type="raw" ), mxFIMLObjective( covariance="expCovMZF", means="expMeanMZF", dimnames=selVars, thresholds="expThMZF")), mxModel("DZF", mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="expCovDZF" ), mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, values=0, name="expMeanDZF" ), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=TRUE, values=c(-1.0, .01, .10, .20), lbound=c(-2.0, .00001, .01, .10), name="V1T"), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=c(TRUE,FALSE,FALSE,FALSE), values=c(-.10,NA,NA,NA), name="V2T"), mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Inc" ), mxAlgebra( expression= cbind((Inc %*% V1T),V2T,(Inc %*% V1T),V2T), dimnames=list(thresh,selVars), name="expThDZF"), mxData( observed=dzfData, type="raw" ), mxFIMLObjective( covariance="expCovDZF", means="expMeanDZF", dimnames=selVars, thresholds="expThDZF")), mxModel("MZM", mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="expCovMZM" ), mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, values=0, name="expMeanMZM" ), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=TRUE, values=c(-2.0, .01, .10, .20), lbound=c(-3.0, .00001, .01, .10), name="V1T"), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=c(TRUE,FALSE,FALSE,FALSE), values=c(.3,NA,NA,NA), name="V2T"), # NOTE Inc means increments between theresholds mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Inc" ), mxAlgebra( expression= cbind((Inc %*% V1T),V2T,(Inc %*% V1T),V2T), dimnames=list(thresh,selVars), name="expThMZM"), mxData( observed=mzmData, type="raw" ), mxFIMLObjective( covariance="expCovMZM", means="expMeanMZM", dimnames=selVars, thresholds="expThMZM")), mxModel("DZM", mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="expCovDZM" ), mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, values=0, name="expMeanDZM" ), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=TRUE, values=c(-2.0, .01, .10, .20), lbound=c(-3.0, .00001, .01, .10), name="V1T"), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=c(TRUE,FALSE,FALSE,FALSE), values=c(.3,NA,NA,NA), name="V2T"), mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Inc" ), mxAlgebra( expression= cbind((Inc %*% V1T),V2T,(Inc %*% V1T),V2T), dimnames=list(thresh,selVars), name="expThDZM"), mxData( observed=dzmData, type="raw" ), mxFIMLObjective( covariance="expCovDZM", means="expMeanDZM", dimnames=selVars, thresholds="expThDZM")), mxModel("DZOS", mxMatrix( type="Stand", nrow=ntv, ncol=ntv, free=TRUE, values=.5, name="expCovDZOS" ), mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, values=0, name="expMeanDZOS" ), #FEMALE VALUES mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=TRUE, values=c(-1.0, .01, .10, .20), lbound=c(-2.0, .00001, .01, .10), name="V1TF"), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=c(TRUE,FALSE,FALSE,FALSE), values=c(-.10,NA,NA,NA), name="V2TF"), #MALE VALUES mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=TRUE, values=c(-2.0, .01, .10, .20), lbound=c(-3.0, .00001, .01, .10), name="V1TM"), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=c(TRUE,FALSE,FALSE,FALSE), values=c(.30,NA,NA,NA), name="V2TM"), # MATRIX TO HOLD INCREMENTS AND ALGEBRA (FEMALES FIRST SINCE THEY'RE IN FIRST COLUMN IN DZOS PAIRS) mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Inc" ), mxAlgebra( expression= cbind((Inc %*% V1TF),V2TF,(Inc %*% V1TM),V2TM), dimnames=list(thresh,selVars), name="expThDZOS"), mxData( observed=dzosData, type="raw" ), mxFIMLObjective( covariance="expCovDZOS", means="expMeanDZOS", dimnames=selVars, thresholds="expThDZOS")), mxAlgebra( MZF.objective + DZF.objective + MZM.objective + DZM.objective + DZOS.objective, name="minus2sumll" ), mxAlgebraObjective("minus2sumll") ) multTwinSatFit <- mxRun(multTwinSatModel) multTwinSatSumm <- summary(multTwinSatFit) multTwinSatSumm # Generate Saturated Output # ----------------------------------------------------------------------- parameterSpecifications(multTwinSatFit) expectedMeansCovariances(multTwinSatFit) LL_Sat <- mxEval(objective, multTwinSatFit) tableFitStatistics(multTwinSatFit) # Fit Multivariate ACE Model using Cholesky Decomposition # ----------------------------------------------------------------------- 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,lbound=c(.00001, NA, .00001), label=c("af11", "af21", "af22"), name="af" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,lbound=c(.00001, NA, .00001), label=c("cf11", "cf21", "cf22"), name="cf" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,lbound=c(.00001, NA, .00001), label=c("ef11", "ef21", "ef22"), name="ef" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,lbound=c(.00001, NA, .00001), label=c("am11", "am21", "am22"), name="am" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,lbound=c(.00001, NA, .00001), label=c("cm11", "cm21", "cm22"), name="cm" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,lbound=c(.00001, NA, .00001), label=c("em11", "em21", "em22"), name="em" ), # Matrices A, C, and E compute variance components mxAlgebra( expression=af %*% t(af), name="Af" ), mxAlgebra( expression=cf %*% t(cf), name="Cf" ), mxAlgebra( expression=ef %*% t(ef), name="Ef" ), mxAlgebra( expression=am %*% t(am), name="Am" ), mxAlgebra( expression=cm %*% t(cm), name="Cm" ), mxAlgebra( expression=em %*% t(em), name="Em" ), # Algebra to compute total variances and standard deviations (diagonal only) mxAlgebra( expression=Af+Cf+Ef, name="Vf" ), mxAlgebra( expression=Am+Cm+Em, name="Vm" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*Vf)), name="isdf"), mxAlgebra( expression=solve(sqrt(I*Vm)), name="isdm"), # Constraint on variance of ordinal variables mxConstraint( Vf == I, name="ConstrainedVf"), mxConstraint( Vm == I, name="ConstrainedVm"), ## 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=ntv, free=FALSE, values= 0, name="expMean"), # Matrix for expected thresholds vectors in FEMALES mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=TRUE, values=c(-1.0, .01, .10, .20), lbound=c(-2.0, .00001, .01, .10), name="V1TF"), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=c(TRUE,FALSE,FALSE,FALSE), values=c(-.10,NA,NA,NA), name="V2TF"), # Matrix for expected thresholds vectors in MALES mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=TRUE, values=c(-2.0, .01, .10, .20), lbound=c(-3.0, .00001, .01, .10), name="V1TM"), mxMatrix(type="Full", nrow=maxthresh, ncol=1, free=c(TRUE,FALSE,FALSE,FALSE), values=c(.30,NA,NA,NA), name="V2TM"), # Matrix to hold incrtements mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Inc" ), # Algebra for thresholds vector in FEMALE PAIRS mxAlgebra( expression= cbind((Inc %*% V1TF),V2TF,(Inc %*% V1TF),V2TF), dimnames=list(thresh,selVars), name="expThreshF"), # Algebra for thresholds vector in MALE PAIRS mxAlgebra( expression= cbind((Inc %*% V1TM),V2TM,(Inc %*% V1TM),V2TM), dimnames=list(thresh,selVars), name="expThreshM"), # Algebra for expected thresholds vector in DZOS PAIRS (FEMALES ARE IN 1ST COLUMN) mxAlgebra( expression= cbind((Inc %*% V1TF),V2TF,(Inc %*% V1TM),V2TM), dimnames=list(thresh,selVars), name="expThreshFM"), # Algebra for expected variance/covariance matrix in MZs F and M mxAlgebra( expression= rbind ( cbind(Af+Cf+Ef , Af+Cf), cbind(Af+Cf , Af+Cf+Ef)), name="expCovMZF" ), mxAlgebra( expression= rbind ( cbind(Am+Cm+Em , Am+Cm), cbind(Am+Cm , Am+Cm+Em)), name="expCovMZM" ), # Algebra for expected variance/covariance matrix in DZs F and M and DZOS, note use of 0.5, converted to 1*1 matrix mxAlgebra( expression= rbind ( cbind(Af+Cf+Ef , 0.5%x%Af+Cf), cbind(0.5%x%Af+Cf , Af+Cf+Ef)), name="expCovDZF" ), mxAlgebra( expression= rbind ( cbind(Am+Cm+Em , 0.5%x%Am+Cm), cbind(0.5%x%Am+Cm , Am+Cm+Em)), name="expCovDZM" ), mxAlgebra( expression= rbind ( cbind(Af+Cf+Ef , 0.5%x%(af%*%t(am))+cf%*%t(cm)), cbind(0.5%x%(am%*%t(af))+cm%*%t(cf),Am+Cm+Em)), name="expCovDZOS" ) ), mxModel("MZF", mxData( observed=mzfData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovMZF", means="ACE.expMean", dimnames=names(mzfData), thresholds="ACE.expThreshF") ), mxModel("DZF", mxData( observed=dzfData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovDZF", means="ACE.expMean", dimnames=names(dzfData), thresholds="ACE.expThreshF") ), mxModel("MZM", mxData( observed=mzmData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovMZM", means="ACE.expMean", dimnames=names(mzmData), thresholds="ACE.expThreshM") ), mxModel("DZM", mxData( observed=dzmData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovDZM", means="ACE.expMean", dimnames=names(dzmData), thresholds="ACE.expThreshM") ), mxModel("DZOS", mxData( observed=dzosData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovDZOS", means="ACE.expMean", dimnames=names(dzosData), thresholds="ACE.expThreshFM") ), mxAlgebra( expression=MZF.objective + DZF.objective + MZM.objective + DZM.objective + DZOS.objective, name="minus2sumll" ), mxAlgebraObjective("minus2sumll") ) multACEFit <- mxRun(multACEModel) multACESumm <- summary(multACEFit) multACESumm # Generate Multivariate ACE Output # ----------------------------------------------------------------------- parameterSpecifications(multACEFit) expectedMeansCovariances(multACEFit) # Path Estimates Matxaf <- mxEval(ACE.af, multACEFit) Matxcf <- mxEval(ACE.cf, multACEFit) Matxef <- mxEval(ACE.ef, multACEFit) Matxaf Matxcf Matxef Matxam <- mxEval(ACE.am, multACEFit) Matxcm <- mxEval(ACE.cm, multACEFit) Matxem <- mxEval(ACE.em, multACEFit) Matxam Matxcm Matxem tableFitStatistics(multACEFit) # EQUATE M & F PARAMETERS TO TEST GENDER HETEROGENEITY # HAVEN'T FIGURED OUT HOW TO EQUATE THRESHOLDS YET # PROBLEM WITH HAVING MULTIPLE START VALUES AND FREE AND FIXED VALUES IN SAME MATRIX # WHEN TRYING TO EQUATE multNoSxACEModel <-mxModel(multACEFit, name="multNoSxACE", mxModel(multACEFit$ACE, # Change labels to equate M & F (no longer "af11" but "a11" etc for M & F) # EQUATE PATH ESTIMATES mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("a11", "a21", "a22"), name="af" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("c11", "c21", "c22"), name="cf" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("e11", "e21", "e22"), name="ef" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("a11", "a21", "a22"), name="am" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("c11", "c21", "c22"), name="cm" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("e11", "e21", "e22"), name="em" )) ) multNoSxACEFit <-mxRun(multNoSxACEModel) multNoSxACESumm <- summary(multNoSxACEFit) parameterSpecifications(multNoSxACEFit) expectedMeansCovariances(multNoSxACEFit) tableFitStatistics(multACEFit, multNoSxACEFit) # ESTIMATE GENDERS SEPARATELY, DROP PARAMETERS TO FIND MOST PARSIMONIOUS MODEL # FEMALES: DROPPING PARAMETERS FROM MultACE # DROP a 2 1 a 2 2 dropFAModel <-mxRename(multACEFit, "dropFA") dropFAModel$ACE.af@free[2,1] <-FALSE dropFAModel$ACE.af@values[2,1] <-0 dropFAModel$ACE.af@free[2,2] <-FALSE dropFAModel$ACE.af@values[2,2] <-0 dropFAFit <-mxRun(dropFAModel) dropFASumm <- summary(dropFAFit) dropFASumm tableFitStatistics(multACEFit, dropFAFit) # DROP c 1 1 and c 2 2 dropFCModel <-mxRename(multACEFit, "dropFC") dropFCModel$ACE.cf@free[1,1]<-FALSE dropFCModel$ACE.cf@values[1,1]<-0 # DROP c 1 1 #dropFCModel$ACE.cf@free[2,1]<-FALSE #dropFCModel$ACE.cf@values[2,1]<-0 # DROP c 1 2 dropFCModel$ACE.cf@free[2,2]<-FALSE dropFCModel$ACE.cf@values[2,2]<-0 # DROP c 2 2 dropFCFit <-mxRun(dropFCModel) dropFCSumm <- summary(dropFCFit) dropFCSumm tableFitStatistics(multACEFit, dropFCFit) Matxcf <- mxEval(ACE.cf, dropFCFit) Matxcf # DROP E 2 1 dropFEModel <-mxRename(multACEFit, "dropFE") dropFEModel$ACE.ef@free[2,1]<-FALSE dropFCModel$ACE.ef@values[2,1]<-0 # DROP E 2 1 dropFEFit <-mxRun(dropFEModel) dropFESumm <- summary(dropFEFit) dropFESumm tableFitStatistics(multACEFit, dropFEFit) # FEMALE FINAL MODEL DROPPING ALL NS PARAMETERS FemModel <- mxRename(multACEFit, "Fem") FemModel$ACE.af@free[2,1] <-FALSE FemModel$ACE.af@values[2,1] <-0 # DROP a 2 1 FemModel$ACE.af@free[2,2] <-FALSE FemModel$ACE.af@values[2,2] <-0 # DROP a 2 2 FemModel$ACE.cf@free[1,1]<-FALSE FemModel$ACE.cf@values[1,1]<-0 # DROP c 1 1 FemModel$ACE.cf@free[2,2]<-FALSE FemModel$ACE.cf@values[2,2]<-0 # DROP c 2 2 FemModel$ACE.ef@free[2,1]<-FALSE FemModel$ACE.ef@values[2,1]<-0 # DROP e 2 1 FemFit <-mxRun(FemModel) FemSumm <- summary(FemFit) FemSumm tableFitStatistics(multACEFit, FemFit) # MALES: DROPPING PARAMETERS FROM MultACE # DROP A 2 1 dropMAModel <-mxRename(multACEFit, "dropMA") dropMAModel$ACE.am@free[2,1] <-FALSE dropMAModel$ACE.am@values[2,1] <-0 dropMAFit <-mxRun(dropMAModel) dropMASumm <- summary(dropMAFit) dropMASumm tableFitStatistics(multACEFit, dropMAFit) # DROP C 2 1 C 2 2 dropMCModel <-mxRename(multACEFit, "dropMC") #dropMCModel$ACE.cm@free[1,1] <-FALSE #dropMCModel$ACE.cm@values[1,1] <-0 dropMCModel$ACE.cm@free[2,1] <-FALSE dropMCModel$ACE.cm@values[2,1] <-0 dropMCModel$ACE.cm@free[2,2] <-FALSE dropMCModel$ACE.cm@values[2,2] <-0 dropMCFit <-mxRun(dropMCModel) dropMCSumm <- summary(dropMCFit) dropMCSumm tableFitStatistics(multACEFit, dropMCFit) # DROP E 2 1 dropMEModel <-mxRename(multACEFit, "dropME") dropMEModel$ACE.em@free[2,1] <-FALSE dropMEModel$ACE.em@values[2,1] <-0 dropMEFit <-mxRun(dropMEModel) dropMESumm <- summary(dropMEFit) dropMESumm tableFitStatistics(multACEFit, dropMEFit) # FINAL MALE MODEL MaleModel <-mxRename(multACEFit, "Male") MaleModel$ACE.am@free[2,1] <-FALSE MaleModel$ACE.am@values[2,1] <-0 MaleModel$ACE.cm@free[2,1] <-FALSE MaleModel$ACE.cm@values[2,1] <-0 MaleModel$ACE.cm@free[2,2] <-FALSE MaleModel$ACE.cm@values[2,2] <-0 MaleModel$ACE.em@free[2,1] <-FALSE MaleModel$ACE.em@values[2,1] <-0 MaleFit <-mxRun(MaleModel) MaleSumm <- summary(MaleFit) MaleSumm tableFitStatistics(multACEFit, MaleFit) # FINAL DROPPING ALL NS PARAMETERS< MALE AND FEMALE FinalModel <-mxRename(multACEFit, "Final") FinalModel$ACE.af@free[2,1] <-FALSE FinalModel$ACE.af@values[2,1] <-0 # DROP a 2 1 FinalModel$ACE.af@free[2,2] <-FALSE FinalModel$ACE.af@values[2,2] <-0 # DROP a 2 2 FinalModel$ACE.cf@free[1,1]<-FALSE FinalModel$ACE.cf@values[1,1]<-0 # DROP c 1 1 FinalModel$ACE.cf@free[2,2]<-FALSE FinalModel$ACE.cf@values[2,2]<-0 # DROP c 2 2 FinalModel$ACE.ef@free[2,1]<-FALSE FinalModel$ACE.ef@values[2,1]<-0 # DROP e 2 1 FinalModel$ACE.am@free[2,1] <-FALSE FinalModel$ACE.am@values[2,1] <-0 FinalModel$ACE.cm@free[2,1] <-FALSE FinalModel$ACE.cm@values[2,1] <-0 FinalModel$ACE.cm@free[2,2] <-FALSE FinalModel$ACE.cm@values[2,2] <-0 FinalModel$ACE.em@free[2,1] <-FALSE FinalModel$ACE.em@values[2,1] <-0 FinalFit <-mxRun(FinalModel) FinalSumm <- summary(FinalFit) FinalSumm tableFitStatistics(multACEFit, FinalFit)