library(OpenMx); library(psych) # source("http://ibg.colorado.edu/cdrom2016/maes/UnivariateAnalysis/miFunctions.R") data(twinData) # Do or do not drop Covariates with NA testNA = TRUE if(!testNA){ twinData = twinData[-which(is.na(twinData$age)),] } # Select Variables for Analysis vars = 'bmi' # list of variables names nv = 1 # number of variables ntv = nv*2 # number of total variables selVars = paste0(vars, c(rep(1,nv), rep(2,nv))) covVars = 'age' # Select Data for Analysis mzData = subset(twinData, zyg==1, c(selVars, covVars)) dzData = subset(twinData, zyg==3, c(selVars, covVars)) # Set Starting Values svMe = 20 # start value for means svPa = .4 # start value for path coefficient svPe = .8 # start value for path coefficient for e lbPa = .0001 # start value for lower bounds defAge = mxMatrix(type="Full" , nrow= 1, ncol = 1, free=FALSE, labels=c("data.age"), name="Age") pathB = mxMatrix(type="Full" , nrow= 1, ncol = 1, free=TRUE, values=.01, label="b11", name="b") meanG = mxMatrix(type="Full" , nrow= 1, ncol=ntv, free=TRUE, values=svMe, labels="xbmi", name="meanG") 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") covA = mxAlgebra(name="A", a %*% t(a)) covC = mxAlgebra(name="C", c %*% t(c)) covE = mxAlgebra(name="E", e %*% t(e)) covP = mxAlgebra(name="V", A+C+E) # Create Model Objects for Multiple Groups pars = list(pathB, meanG, pathA, pathC, pathE, covA, covC, covE, covP) modelMZ = mxModel("MZ", pars, defAge, mxAlgebra(name="expMean", meanG + cbind(b%*%Age,b%*%Age)), mxAlgebra(A+C, name="cMZ"), mxAlgebra(rbind(cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ"), mxData(mzData, type="raw"), mxExpectationNormal("expCovMZ", means="expMean", dimnames=selVars), mxFitFunctionML() ) modelDZ = mxModel("DZ", pars, defAge, mxAlgebra(name = "expMean", meanG + cbind(b%*%Age, b%*%Age)), mxAlgebra(0.5%x%A+ C, name="cDZ"), mxAlgebra(rbind(cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ"), mxData(dzData, type="raw"), mxExpectationNormal("expCovDZ", means="expMean", dimnames=selVars), mxFitFunctionML() ) ACE = mxModel("ACE", pars, modelMZ, modelDZ, mxFitFunctionMultigroup(c("MZ","DZ")) ) ACE = mxRun(ACE); summary(ACE) # =============================================================== # = Runs with NAs in def vars and gives start values as results = # =============================================================== # Running ACE with 5 parameters # Warning message: # In model 'ACE' Optimizer returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) # Summary of ACE # # The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED) # # free parameters: # name matrix row col Estimate lbound ubound # 1 b11 b 1 1 0.01 # 2 xbmi meanG 1 1 20.00 # 3 a11 a 1 1 0.40 1e-04 # 4 c11 c 1 1 0.40 1e-04 # 5 e11 e 1 1 0.80 1e-04 # # Model Statistics: # | Parameters | Degrees of Freedom | Fit (-2lnL units) # Model: 5 1772 7.686145e+14 # Saturated: NA NA NA # Independence: NA NA NA # Number of observations/statistics: 920/1777 # # # ** Information matrix is not positive definite (not at a candidate optimum). # Be suspicious of these results. At minimum, do not trust the standard errors. # # Information Criteria: # | df Penalty | Parameters Penalty | Sample-Size Adjusted # AIC: 7.686145e+14 7.686145e+14 NA # BIC: 7.686145e+14 7.686145e+14 7.686145e+14 # CFI: NA # TLI: 1 (also known as NNFI) # RMSEA: 0 [95% CI (NA, NA)] # Prob(RMSEA <= 0.05): NA # To get additional fit indices, see help(mxRefModels) # timestamp: 2017-08-29 21:20:45 # Wall clock time (HH:MM:SS.hh): 00:00:00.97 # optimizer: CSOLNP # OpenMx version number: 2.7.12.47 # Need help? See help(mxSummary)