rm(list=ls()) library(OpenMx) library(psych); library(polycor) source("miFunctions.R") # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA # Load Data data(twinData) dim(twinData) describe(twinData[,1:12], skew=F) twinData[,'ht1'] <- twinData[,'ht1']*10 twinData[,'ht2'] <- twinData[,'ht2']*10 twinData[,'wt1'] <- twinData[,'wt1']/10 twinData[,'wt2'] <- twinData[,'wt2']/10 # Select Variables for Analysis vars <- c('ht','wt') # list of variables names 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 mzData <- subset(twinData, zyg==1, selVars) dzData <- subset(twinData, zyg==3, selVars) m <- 2*2 l <- 3*4 p <- m+l matA <- mxMatrix(type = "Full", nrow = p, ncol = p, byrow = TRUE, values = c(0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0, 0,0,0,0,1,1,0,0,1,1,0,0,1,1,0,0, 0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0, 0,0,0,0,0,0,1,1,0,0,1,1,0,0,1,1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), free = c(FALSE,FALSE,FALSE,FALSE,TRUE,FALSE,FALSE,FALSE,TRUE,FALSE,FALSE,FALSE,TRUE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,TRUE,TRUE,FALSE,FALSE,TRUE,TRUE,FALSE,FALSE,TRUE,TRUE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,TRUE,FALSE,FALSE,FALSE,TRUE,FALSE,FALSE,FALSE,TRUE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,TRUE,TRUE,FALSE,FALSE,TRUE,TRUE,FALSE,FALSE,TRUE,TRUE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE, FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE), labels = c(NA,NA,NA,NA,"axx",NA,NA,NA,"cxx",NA,NA,NA,"exx",NA,NA,NA, NA,NA,NA,NA,"ayx","ayy",NA,NA,"cyx","cyy",NA,NA,"eyx","eyy",NA,NA, NA,NA,NA,NA,NA,NA,"axx",NA,NA,NA,"cxx",NA,NA,NA,"exx",NA, NA,NA,NA,NA,NA,NA,"ayx","ayy",NA,NA,"cyx","cyy",NA,NA,"eyx","eyy", NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA), name = "A") matSM <- mxMatrix(type = "Full", nrow = p, ncol = p, byrow = TRUE, values = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0, 0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1), free = FALSE, name = "SM", lbound = 0.0001) matSD <- mxMatrix(type = "Full", nrow = p, ncol = p, byrow = TRUE, values = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,1,0,0.5,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0.5,0,0,0,0,0,0,0,0, 0,0,0,0,0.5,0,1,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0.5,0,1,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1), free = FALSE, name = "SD", lbound = 0.0001) matF <- mxMatrix(type = "Full", nrow = m, ncol = p, byrow = TRUE, free = FALSE, values = c(1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0), name= "F") matI <- mxMatrix(type = "Iden", nrow = p, ncol = p, free = FALSE, name = "I") meanG <- mxMatrix( type="Full", nrow=1, ncol=m, free=TRUE, values=1, labels=c("meanx","meany","meanx","meany"), name="meanG" ) # Create expected covariance matrix as a function of the A, S and F matrices (following the RAM-Formula) expCovMZ <- mxAlgebra(expression= F%*%(solve(I-A))%*%SM%*%(t(solve(I-A)))%*%t(F), name="expCovMZ" ) expCovDZ <- mxAlgebra(expression= F%*%(solve(I-A))%*%SD%*%(t(solve(I-A)))%*%t(F), 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="meanG", dimnames=selVars ) expDZ <- mxExpectationNormal(covariance="expCovDZ", means="meanG", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups parameters <- list(meanG, matA, matF, matI) modelMZ <- mxModel(parameters, matSM, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel(parameters, matSD, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Build Model with Confidence Intervals modelACE <- mxModel( "oneACEc", parameters, modelMZ, modelDZ, multi) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run ACE Model set.seed(1) fitACE <- mxTryHard(modelACE, extraTries = 10, exhaustive = TRUE) sumACE <- summary( fitACE ) sumACE mxCheckIdentification(fitACE)