rm(list=ls()) #ls() require(OpenMx) require(psych) library(foreign) require(readstata13) library("OpenMx",lib="~/R/lib3.5.1") library("psych",lib="~/R/lib3.5.1") library("polycor",lib="~/R/lib3.5.1") library("readstata13",lib="~/R/lib3.5.1") library("mvtnorm", lib="~/R/lib3.5.1") # change optimizer ######## #mxOption(NULL, "Default optimizer", "NPSOL") mxOption( NULL, 'Default optimizer' ,'SLSQP' ) #mxOption(NULL, "Default optimizer", "CSOLNP") mxOption( NULL, "mvnRelEps", value= mxOption(NULL, key="mvnRelEps")/5) getOption('mxOptions')$mvnRelEps setwd("E:/奚玉?-前瞻?/新建文件?/usedata") twinData <- read.dta13("登记-18岁成年双胞胎-11391pair-随访-高血压同性别瘦正常超?.dta") twinData$age1 = scale(twinData$age1) twinData$age2 = scale(twinData$age2) twinData$Bmi1 = scale(twinData$bmi1) twinData$Bmi2 = scale(twinData$bmi2) ########################################################### ### select continuous variables varsc <-c('Bmi') nvc <- 1 # number of covariates ntvc <- nvc*2 # number of variables per pair # number of max thresholds = the number of category - 1 conVars <-paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="") #c('tea2g1', 'smk_cur1', 'tea2g2', 'smk_cur2' )# ### select binary variables varso <- c('hypertension') nth <- 1 nvo <- 1 # number of covariates ntvo <- nvo*2 # number of variables per pair ordVars <-paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") #c('tea2g1', 'smk_cur1', 'tea2g2', 'smk_cur2' )# ### select variables vars <-c('Bmi','hypertension') covars <- c("sex", "age") nv <- nvc+nvo # number of variables per twin ncv <- 2 # number of covariates ntv <- nv*2 # number of variables per pair selVars <-paste(vars,c(rep(1,nv),rep(2,nv)),sep="") #c('tea2g1', 'smk_cur1', 'tea2g2', 'smk_cur2' )# covVars <- paste(covars,c(rep(1,ncv),rep(2,ncv)),sep="") # Select Data for Analysis mzData <- subset(twinData, zy_sex1==1|zy_sex1==2, c(selVars,covVars)) dzData <- subset(twinData, zy_sex1==3|zy_sex1==4, c(selVars,covVars)) sum(is.na(mzData)) sum(is.na(dzData)) mzData <- na.omit(mzData) dzData <- na.omit(dzData) dim(mzData) #8454,3564 dim(dzData) # Generate Descriptive Statistics mzData <-cbind(mzData[,conVars],mzData[,covVars],mxFactor( x=mzData[,ordVars], levels=c(0:nth)) ) dzData <-cbind(dzData[,conVars],dzData[,covVars],mxFactor( x=dzData[,ordVars], levels=c(0:nth)) ) table(mzData$rehypertension1,mzData$rehypertension2) table(dzData$rehypertension1,dzData$rehypertension2) # Set Starting Values # Set Starting Values frMV <- c(TRUE,FALSE) #frCvD <- valDiagLU(frMV,T,T,ntv) # free status for diagonal, lower & upper elements of covariance matrix frCvD <- diag(frMV,ntv,ntv) # values for diagonal of covariance matrix frCvD[lower.tri(frCvD)] <- T# values for below diagonal elements frCvD[upper.tri(frCvD)] <- T# values for above diagonal elements frCvD frCv <- matrix(as.logical(frCvD),4) svMe <- c(0,0) # start value for means svVa <- c(.5,1) # start value for variances #valDiag <- function(valD,dim) { # valF <- diag(valD,dim,dim) # values for diagonal of covariance matrix # valF #} svVaC <- diag(svVa,ntv,ntv) lbVa <- .0001 # lower bound for variances lbVaC <- diag(lbVa,ntv,ntv) svLTh <- 0 # start value for first threshold svITh <- 1 # start value for increments svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bound for thresholds svb <- 0.5 # sv for beta,get from tryhardordinal svTh <- 1.5 # Create Labels labMeMZ <- paste("meanMZ",selVars,sep="") labMeDZ <- paste("meanDZ",selVars,sep="") labMeZ <- paste("meanZ",vars,sep="") labCvMZ <- paste("covMZ",rev(ntv+1-sequence(1:ntv)),rep(1:ntv,ntv:1),sep="") labCvDZ <- paste("covDZ",rev(ntv+1-sequence(1:ntv)),rep(1:ntv,ntv:1),sep="") labCvZ <- paste("covZ",rev(ntv+1-sequence(1:ntv)),rep(1:ntv,ntv:1),sep="") labVaMZ <- paste("covMZ",1:ntv,1:ntv,sep="") labVaDZ <- paste("covDZ",1:ntv,1:ntv,sep="") labVaZ <- paste("covZ",1:ntv,1:ntv,sep="") labThMZ <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="") labThDZ <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="") labThZ <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="") labdata <- paste("data.",covVars,sep = "")#cbind(paste("data.",covars[1],c(1,2),sep = ""),paste("data.",covars[2],c(1,2),sep = "")) #c("data.age1","data.age2","data.region_type1","data.region_type2") labeta <-cbind(paste("B",covars,"_",vars[1],sep = ""),paste("B",covars[2],"_",vars[2],sep = "")) #c('Bage_tea2g','Bregion_type_tea2g','Bage_smkcur', 'Bregion_type_smkcur') betaLabs_age <- paste("beta","Age",vars, sep="_") betaLabs_sex <- paste("beta","Sex",vars, sep="_") ########################################## # tea2g1 smkcur1 | tea2g2 smkcur2 (the sequence of setting labels ??? ) # tea2g1 VP1T1 | # smkcur1 CVT1* VP2T1 | #----------------------------------------- # tea2g2 WP1* CTT2T1 | VP1T2 # smkcur2 CTT1T2* WP2* | CVT2 VP2T2 ########################################## # VP1T1(variance of pheno1 of twin1) VP2T1(variance of pheno2 of twin1) VP1T2(variance of pheno1 of twin2) VP2T2(variance of pheno2 of twin2) # CVT1(covar pheno1 pheno2 of twin1) CVT2(covar pheno1 pheno2 of twin2) # WP1(covar pheno1 of twin1 twin2) WP2(covar pheno2 of twin1 twin2) # CTT1T2(pheno1 of twin1 and pheno2 of twin2) CTT2T1(pheno1 of twin2 and pheno2 of twin1) # Define definition variables to hold the Covariates # Matrix & Algebra for expected means (SND), Thresholds, effect of Age on Th and correlations modelMZ <- mxModel( "MZ", mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_age, name="bAge" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.age1","data.age2"), name="Age"), mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_sex, name="bSex" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex"), mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeMZ, name="meanG" ), mxAlgebra( expression= meanG + (Age %x% bAge) + (Sex %x% bSex), name="meanMZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labThMZ, name="thinMZ" ), mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ), mxAlgebra( expression= inc %*% thinMZ, name="threMZ" ), mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVaC, lbound=lbVaC, labels=labCvMZ, name="covMZ" ), mxData( observed=mzData, type="raw" ), mxExpectationNormal( covariance="covMZ", means="meanMZ", dimnames=selVars, thresholds="threMZ", threshnames=ordVars ), mxFitFunctionML() ) modelDZ <- mxModel( "DZ", mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_age, name="bAge" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.age1","data.age2"), name="Age"), mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_sex, name="bSex" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex"), mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZ, name="meanG" ), mxAlgebra( expression= meanG + (Age %x% bAge) + (Sex %x% bSex), name="meanDZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labThDZ, name="thinDZ" ), mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ), mxAlgebra( expression= inc %*% thinDZ, name="threDZ" ), mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVaC, lbound=lbVaC, labels=labCvDZ, name="covDZ" ), mxData( observed=dzData, type="raw" ), mxExpectationNormal( covariance="covDZ", means="meanDZ", dimnames=selVars, thresholds="threDZ", threshnames=ordVars ), mxFitFunctionML() ) # Combine Groups Conf <- mxCI (c ('MZ.covMZ','DZ.covDZ') )#'MZ.expCovMZ[3,1]','MZ.expCovMZ[4,2]' SatModel <- mxModel( "Sat", modelMZ, modelDZ,mxFitFunctionMultigroup(c('MZ','DZ')), Conf ) # ------------------------------------------------------------------------------------------------------------------------------- # 1) RUN Saturated Model #SatFit <- mxRun(SatModel, intervals=TRUE) #set.seed(180316) #number means nothing,don't need to change SatFit <- mxTryHardOrdinal(SatModel, intervals=T,extraTries = 15) #,bestInitsOutput=T,silent=F,modify start values (SatSumm <- summary(SatFit,verbose = T)) #warnings() #library(numDeriv) mxCheckIdentification(SatFit, details=T) #, details=TRUE # Generate SatModel Output # round(SatFit$output$estimate,4) # round(SatFit$MZ.expMeanMZ$result,4) round(SatSumm$CI[,c(2,1,3)] , 3) #in the order: Estimate, lower bound, upper bound# # ----------------------x----------------------------------------------------------------------------- # 2) Run SubModels: # Test Significance of Covariates (result:no) dropb <- omxSetParameters( SatModel, label=c(betaLabs_sex,betaLabs_age), free=FALSE, values=0,name="dropb" ) dropbFit <- mxTryHardOrdinal(dropb, intervals=F,extraTries = 20) mxCompare(SatFit, dropbFit ) ### no ## # Constrain expected means (BMI) to be equal across twin order mean1Model <- mxModel( SatModel, name="mean1" ) mean1Model <- omxSetParameters( mean1Model, label=c(labMeMZ[1],labMeMZ[3]), free=frMV[1], values=svMe[1], newlabels=labMeMZ[1] ) mean1Model <- omxSetParameters( mean1Model, label=c(labMeDZ[1],labMeDZ[3]), free=frMV[1], values=svMe[1], newlabels=labMeDZ[1] ) mean1Fit <- mxTryHardOrdinal(mean1Model, intervals=F,extraTries = 15) mxCompare(SatFit, subs <- list(dropbFit,mean1Fit ) ) (mean1Summ <- summary(mean1Fit)) ### yes ## # Constrain expected variance (BMI) to be equal across twin order variance1Model <- mxModel( mean1Model, name="variance1" ) variance1Model <- omxSetParameters( variance1Model, label=c(labVaMZ[1],labVaMZ[3]), free=frMV[1], values=svVa[1],newlabels=labVaMZ[1] ) variance1Model <- omxSetParameters( variance1Model, label=c(labVaDZ[1],labVaDZ[3]), free=frMV[1], values=svVa[1],newlabels=labVaDZ[1] ) variance1Fit <- mxTryHardOrdinal(variance1Model, intervals=F,extraTries = 15) mxCompare(SatFit, subs <- list(dropbFit,mean1Fit,variance1Fit ) ) (variance1Summ <- summary(variance1Fit)) ### yes ## # Constrain expected Thresholds (coronary) to be equal across twin order thresh1Model <- mxModel( variance1Model, name="thresh1" ) thresh1Model <- omxSetParameters( thresh1Model, label=c(labThMZ[1],labThMZ[2]), free=TRUE, values=svTh, newlabels=labThMZ[1] ) thresh1Model <- omxSetParameters( thresh1Model, label=c(labThDZ[1],labThDZ[2]), free=TRUE, values=svTh, newlabels=labThDZ[1] ) thresh1Fit <- mxTryHardOrdinal(thresh1Model, intervals=F,extraTries = 15) mxCompare(SatFit, subs <- list(dropbFit,mean1Fit,variance1Fit ,thresh1Fit ) ) (thresh1Summ <- summary(thresh1Fit)) ### yes ## # Constrain expected means (BMI) to be equal across zygosity meanz1Model <- mxModel( thresh1Model, name="meanz1" ) meanz1Model <- omxSetParameters( meanz1Model, label=c(labMeMZ[1],labMeDZ[1]), free=frMV[1], values=svMe[1], newlabels=labMeMZ[1] ) meanz1Fit <- mxTryHardOrdinal(meanz1Model, intervals=F,extraTries = 15) mxCompare(SatFit, subs <- list(dropbFit,mean1Fit,variance1Fit ,thresh1Fit,meanz1Fit ) ) (meanz1Summ <- summary(meanz1Fit)) ### no ## # Constrain expected variance (BMI) to be equal across zygosity variancez1Model <- mxModel( thresh1Model, name="variancez1" ) variancez1Model <- omxSetParameters( variancez1Model, label=c(labVaMZ[1],labVaDZ[1]), free=frMV[1], values=svVa[1],newlabels=labVaMZ[1] ) variancez1Fit <- mxTryHardOrdinal(variancez1Model, intervals=F,extraTries = 15) mxCompare(SatFit, subs <- list(dropbFit,mean1Fit,variance1Fit ,thresh1Fit,meanz1Fit,variancez1Fit ) ) (variancez1Summ <- summary(variancez1Fit)) ### yes ## # Constrain expected Thresholds (coronary) to be equal across zygosity threshz1Model <- mxModel( variancez1Model, name="threshz1" ) threshz1Model <- omxSetParameters( threshz1Model, label=c(labThMZ[1],labThDZ[1]), free=TRUE, values=svTh, newlabels=labThMZ[1] ) threshz1Fit <- mxTryHardOrdinal(threshz1Model, intervals=F,extraTries = 15) mxCompare(SatFit, subs <- list(dropbFit,mean1Fit,variance1Fit ,thresh1Fit,meanz1Fit,variancez1Fit,threshz1Fit ) ) (thresh1Summ <- summary(threshz1Fit)) ### yes ## ########################################## # tea2g1 smkcur1| tea2g2 smkcur2 # tea2g1 VP1T1 | # smkcur1 CVT1* VP2T1 | #----------------------------------------- # tea2g2 WP1* CTT2T1 | VP1T2 # smkcur2 CTT1T2* WP2* | CVT2 VP2T2 ########################################## #---------------------------------------- # 强制另方差相等输出系??? variance2Model <- mxModel(threshz1Model, name="variance2" ) variance2Model <- omxSetParameters( variance2Model, label=c(labCvMZ[2],labCvMZ[9]), free=frMV[1],newlabels='MZCV' ) variance2Model <- omxSetParameters( variance2Model, label=c(labCvDZ[2],labCvDZ[9]), free=frMV[1],newlabels='DZCV' ) variance2Model <- omxSetParameters( variance2Model, label=c(labCvMZ[4],labCvMZ[6]), free=frMV[1],newlabels='MZCTCT' ) variance2Model <- omxSetParameters( variance2Model, label=c(labCvDZ[4],labCvDZ[6]), free=frMV[1],newlabels='DZCTCT' ) variance2Fit <- mxTryHardOrdinal(variance2Model, intervals=T,extraTries = 15) mxCompare(SatFit, subs <- list(dropbFit,variance2Fit ) ) summary(variance2Fit,verbose = T) tecor <- round(Summ1$CI[c(2,1,3)] , 3 ) for( i in 1:length(tecor[,1])){ tecor[i,1] <- paste(tecor[i,1],"(",tecor[i,2],",",tecor[i,3],")",sep = "") } tecor_vars <- data.frame(type=c("BMI_CHD","-"), zy=c("MZ","DZ"), BMI=c(tecor[1,1],tecor[2,1]), CHD=c(tecor[3,1],tecor[4,1]), CV=c(tecor[5,1],tecor[6,1]), CTCT=c(tecor[7,1],tecor[8,1])) #sink('result/tecor_biv_4.doc',append = T) tecor_vars summary(Fit,verbose = T) #sink() ################################################################################ ### ACE-model ################################################################################ ### State the Quantitative genetic model lbPa <- .0001 # start value for lower bounds lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements laThMZ <- c("thZ_hyp1","thZ_hyp1") laThDZ <- c("thZ_hyp1","thZ_hyp1") labMeMZ <- paste("meanMZ",vars,sep="") labMeDZ <- paste("meanDZ",vars,sep="") labMeZ <- paste("meanZ",vars,sep="") ## Labeling aLabs <- c("aM","aC","aU") cLabs <- c("cM","cC","cU") eLabs <- c("eM","eC","eU") svPa <- .5 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .7 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix ### Some parameters included in all "submodels": baseACE <- mxModel('Base', mxMatrix(name = "a", type = "Lower", nrow = nv, ncol = nv, free=T, labels = aLabs, values=svPaD, lbound=lbPaD), mxMatrix(name = "c", type = "Lower", nrow = nv, ncol = nv, free=T, labels = cLabs, values=svPaD, lbound=lbPaD), mxMatrix(name = "e", type = "Lower", nrow = nv, ncol = nv, free=T, labels = eLabs, values=svPeD, lbound=lbPaD), mxAlgebra(name = "A", expression = a %*% t(a)), mxAlgebra(name = "C", expression = c %*% t(c)), mxAlgebra(name = "E", expression = e %*% t(e)), mxAlgebra( expression=A+C+E, name="V" ), # Constraint on variance of Binary variables mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv1" ), mxConstraint(expression=diag2vec(V) == Unv1, name="Var1"), # Algebra to compute standardized variance components mxAlgebra( expression=A/V, name="h2"), mxAlgebra( expression=C/V, name="c2"), mxAlgebra( expression=E/V, name="e2"), # Algebra to compute total variances and standard deviations (tea2g_ggonal only) mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"), mxAlgebra( expression=solve(sqrt(I*V))%&%V, name ="Rph"), mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="Ra" ), #cov2cor() mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="Rc" ), mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="Re" ), # the proportion of ace contributing to the overlap # IF any r is negative,I agree that it makes more sense to report it as the contribution to the phenotypic correlation rather than # the proportion. So I will not divide by the phenotypic correlation! mxAlgebra( expression=sqrt(h2[1,1])*sqrt(h2[2,2])*Ra[1,2]/Rph[1,2],name = "pa"), mxAlgebra( expression=sqrt(c2[1,1])*sqrt(c2[2,2])*Rc[1,2]/Rph[1,2],name = "pc"), mxAlgebra( expression=sqrt(e2[1,1])*sqrt(e2[2,2])*Re[1,2]/Rph[1,2],name = "pe") ) modelMZ <- mxModel( "MZ", mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_age, name="bAge" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.age1","data.age2"), name="Age"), mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_sex, name="bSex" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex"), mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeMZ, name="meanG" ), mxAlgebra( expression= meanG + (Age %x% bAge) + (Sex %x% bSex), name="meanMZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=laThMZ, name="thinMZ" ), mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ), mxAlgebra( expression= inc %*% thinMZ, name="threMZ" ), mxAlgebra( expression= rbind( cbind(Base.A+Base.C+Base.E , Base.A+Base.C), cbind(Base.A+Base.C , Base.A+Base.C+Base.E)), name="expCovMZ" ), mxData( observed=mzData, type="raw" ), mxExpectationNormal( covariance="expCovMZ", means="meanMZ", dimnames=selVars, thresholds="thinMZ" , threshnames=ordVars), mxFitFunctionML() ) modelDZ <- mxModel( "DZ", mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_age, name="bAge" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.age1","data.age2"), name="Age"), mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, labels=betaLabs_sex, name="bSex" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.sex1","data.sex2"), name="Sex"), mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZ, name="meanG" ), mxAlgebra( expression= meanG + (Age %x% bAge) + (Sex %x% bSex), name="meanDZ" ), mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=laThDZ, name="thinDZ" ), mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" ), mxAlgebra( expression= inc %*% thinDZ, name="threDZ" ), mxAlgebra( expression= rbind( cbind(Base.A+Base.C+Base.E , 0.5%x%Base.A+Base.C), cbind(0.5%x%Base.A+Base.C , Base.A+Base.C+Base.E)), name="expCovDZ" ), mxData( observed=dzData, type="raw" ), mxExpectationNormal( covariance="expCovDZ", means="meanDZ", dimnames=selVars, thresholds="thinDZ" , threshnames=ordVars), mxFitFunctionML() ) Conf1 <- mxCI (c ('Base.h2[1,1]','Base.c2[1,1]','Base.e2[1,1]','Base.h2[2,2]','Base.c2[2,2]','Base.e2[2,2]','Base.Rph[2,1]', 'Base.Re[2,1]',"Base.pe",'Base.Ra[2,1]',"Base.pa")) Confc <- mxCI (c ('Base.Rc[2,1]',"Base.pc") ) AceModel <- mxModel( "ACE", baseACE,modelMZ, modelDZ,Conf1, mxFitFunctionMultigroup(c('MZ.fitfunction','DZ.fitfunction')))# Confa,Confc, # ------------------------------------------------------------------------------ # 4) RUN AceModel AceCI <- mxModel(AceModel,Confc) set.seed(180316) AceFit <- mxTryHardOrdinal(AceCI, intervals=T,extraTries = 15,exhaustive=T) (AceSumm <- summary(AceFit)) mxCompare(SatFit,AceFit) summary(AceFit,verbose = T)