SATURATED MODEL describe(twinData, skew=T) # Select Variables for Analysis vars <- 'MidPeriphery' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis covVars <- c('age', "Sex") nc <- 1 # number of covariates # Select Data for Analysis mzData <- subset(twinData, Zyg==1, c(selVars, covVars)) dzData <- subset(twinData, Zyg==2, c(selVars, covVars)) # Set Starting Values svMe <- -0.2 # start value for means svVa <- 1 # start value for variance svVas <- diag(svVa,ntv,ntv) # assign start values to diagonal of matrix lbVa <- .000000000000000000000001 # start value for lower bounds lbVas <- diag(lbVa,ntv,ntv) # assign lower bounds to diagonal of matrix # Create Labels labMeMZ <- c("mMZ1","mMZ2") # labels for means for MZ twins labMeDZ <- c("mDZ1","mDZ2") # labels for means for DZ twins labCvMZ <- c("vMZ1","cMZ21","vMZ2") # labels for (co)variances for MZ twins labCvDZ <- c("vDZ1","cDZ21","vDZ2") # labels for (co)variances for DZ twins labVaMZ <- c("vMZ1","vMZ2") # labels for variances for MZ twins labVaDZ <- c("vDZ1","vDZ2") # labels for variances for DZ twins # ------------------------------------------------------------------------------ # PREPARE MODEL # Saturated Model # Saturated Model defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="b11", name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Sex"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="b12", name="b2" ) # Create Algebra for expected Mean Matrices meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZ, name="meanMZ" ) meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZ, name="meanDZ" ) expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(b1%*%Age,b1%*%Age)+cbind(b2%*%Sex,b2%*%Sex), name="expMeanMZ" ) expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(b1%*%Age,b1%*%Age)+cbind(b2%*%Sex,b2%*%Sex), name="expMeanDZ" ) # Create Algebra for expected Variance/Covariance Matrices covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" ) covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" ) # Create Algebra for Maximum Likelihood Estimates of Twin Correlations matI <- mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="I" ) corMZ <- mxAlgebra( solve(sqrt(I*covMZ)) %&% covMZ, name="corMZ" ) corDZ <- mxAlgebra( solve(sqrt(I*covDZ)) %&% covDZ, name="corDZ" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="covMZ", means="expMeanMZ", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(pathB1, pathB2) defs <- list(defAge, defSex) modelMZ <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, matI, corMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, matI, corDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Confidence Interval Objects ciCov <- mxCI( c('MZ.covMZ','DZ.covDZ' )) ciMean <- mxCI( c('MZ.expMeanMZ','DZ.expMeanDZ' )) MZZZ <-mxAlgebra(expression =cov2cor(MZ.corMZ), name="MZZ") DZZZ <-mxAlgebra(expression =cov2cor(DZ.corDZ), name="DZZ") ci2 <- mxCI( c('MZZ','DZZ') ) # Build Saturated Model with Confidence Intervals model <- mxModel( "oneSATca", pars, modelMZ, modelDZ, multi, ciCov, ciMean, ci2, MZZZ, DZZZ ) # ------------------------------------------------------------------------------ # RUN MODEL # Run Saturated Model fit <- mxRun( model, intervals=T ) sum <- summary( fit ) round(fit$output$estimate,4) fitGofs(fit) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Test Significance of Covariate modelCov <- mxModel( fit, name="testCov" ) modelCov <- omxSetParameters( modelCov, label=c("b11","b12"), free=FALSE, values=0 ) fitCov <- mxRun( modelCov ) mxCompare( fit, fitCov ) # Constrain expected Means to be equal across twin order modelEMO <- mxModel( fit, name="eqMeansTwin" ) modelEMO <- omxSetParameters( modelEMO, label=labMeMZ, free=TRUE, values=svMe, newlabels='mMZ' ) modelEMO <- omxSetParameters( modelEMO, label=labMeDZ, free=TRUE, values=svMe, newlabels='mDZ' ) fitEMO <- mxRun( modelEMO, intervals=F ) mxCompare( fit, subs <- list(fitCov, fitEMO) ) round(fitEMO$output$estimate,4) fitGofs(fitEMO) # Constrain expected Means and Variances to be equal across twin order modelEMVO <- mxModel( fitEMO, name="eqMVarsTwin" ) modelEMVO <- omxSetParameters( modelEMVO, label=labVaMZ, free=TRUE, values=svVa, newlabels='vMZ' ) modelEMVO <- omxSetParameters( modelEMVO, label=labVaDZ, free=TRUE, values=svVa, newlabels='vDZ' ) fitEMVO <- mxRun( modelEMVO, intervals=F ) mxCompare( fit, subs <- list(fitCov, fitEMO, fitEMVO) ) round(fitEMVO$output$estimate,4) fitGofs(fitEMVO) # Constrain expected Means and Variances to be equal across twin order and zygosity modelEMVZ <- mxModel( fitEMVO, name="eqMVarsZyg" ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("mMZ","mDZ"), free=TRUE, values=svMe, newlabels='mZ' ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("vMZ","vDZ"), free=TRUE, values=svVa, newlabels='vZ' ) fitEMVZ <- mxTryHard( modelEMVZ, intervals=F ) mxCompare( fit, subs <- list(fitCov, fitEMO, fitEMVO, fitEMVZ) ) round(fitEMVZ$output$estimate,4) fitGofs(fitEMVZ) mxEval(cov2cor(DZ.covDZ), fit, T) mxEval(cov2cor(MZ.covMZ), fit, T) ACE MODEL describe(twinData, skew=F) # Load Data # Select Variables for Analysis vars <- 'MidPeriphery' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis # Select Covariates for Analysis covVars <- c('age', "Sex") nc <- 1 # Select Data for Analysis mzData <- subset(twinData, Zyg==1, c(selVars, covVars)) dzData <- subset(twinData, Zyg==2, c(selVars, covVars)) # Set Starting Values svMe <- -3 # start value for means svPa <- .8 # start value for path coefficient svPe <- .8 # start value for path coefficient for e lbPa <- .00000001 # start value for lower bounds # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="b11", name="b1" ) defSex <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Sex"), name="Sex" ) pathB2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, label="b12", name="b2" ) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels="xbmi", name="meanG" ) expMean <- mxAlgebra( expression= meanG + cbind(b1%*%Age,b1%*%Age)+cbind(b2%*%Sex,b2%*%Sex), name="expMean" ) # Create Matrices for Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="a11", lbound=lbPa, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="c11", lbound=lbPa, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="e11", lbound=lbPa, name="e" ) # Create Algebra for 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" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= A+C+E, name="V" ) covMZ <- mxAlgebra( expression= A+C, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(pathB1, pathB2, meanG, pathA, pathC, pathE, covA, covC, covE, covP) defs <- list(defAge, defSex) modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Variance Components rowVC <- rep('VC',nv) colVC <- rep(c('A','C','E','SA','SC','SE'),each=nv) estVC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC)) # Create Confidence Interval Objects ciACE <- mxCI( "VC" ) # Build Model with Confidence Intervals modelACE <- mxModel( "oneACEca", pars, modelMZ, modelDZ, multi, estVC, ciACE ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE Model fitACE <- mxTryHard( modelACE, intervals=T ) sumACE <- summary( fitACE ) # Compare with Saturated Model mxCompare( fit, fitACE ) #lrtSAT(fitACE,-2ll,df) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEsts(fitACE) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Test Significance of Covariate modelCov <- mxModel( fitACE, name="testCov" ) modelCov <- omxSetParameters( modelCov, label="b11", free=FALSE, values=0 ) fitCov <- mxRun( modelCov ) mxCompare( fitACE, fitCov ) # Run AE model modelAE <- mxModel( fitACE, name="oneAEca" ) modelAE <- omxSetParameters( modelAE, labels="c11", free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=T ) mxCompare( fitACE, fitAE ) fitGofs(fitAE) fitEsts(fitAE) # Run CE model modelCE <- mxModel( fitACE, name="oneCEca" ) modelCE <- omxSetParameters( modelCE, labels="a11", free=FALSE, values=0 ) fitCE <- mxRun( modelCE, intervals=T ) mxCompare( fitACE, fitCE ) fitGofs(fitCE) fitEsts(fitCE) # Run E model modelE <- mxModel( fitAE, name="oneEca" ) modelE <- omxSetParameters( modelE, labels="a11", free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) mxCompare( fitAE, fitE ) fitGofs(fitE) fitEsts(fitE)