> require(OpenMx) > require(psych) > source('GenEpiHelperFunctions.R') > mxOption(NULL,"Default optimizer","NPSOL") > vars <- 'mpos' # list of variables names > modVars <- c('age') # moderator variable > nv <- 1 # number of phenotype variables per twin > ntv <- nv*2 # total number of phenotype variables > na <- 1 # number of A factors (per twin) > nc <- 1 # number of C factors (per twin) > ne <- 1 # number of E factors (per twin) > selVars <- paste(vars,c(rep(0,nv),rep(1,nv)),sep="") > mzData <- recipmod[recipmod$zygosity==1, c(selVars,modVars)] > dzData <- recipmod[recipmod$zygosity==2, c(selVars,modVars)] > # Matrices to store a, c, and e Path Coefficients > pathA <- mxMatrix( type="Full", nrow=nv, ncol=na, free=TRUE, values=.1, + label="a11", name="a" ) > pathC <- mxMatrix( type="Full", nrow=nv, ncol=nc, free=TRUE, values=.1, + label="c11", name="c" ) > pathE <- mxMatrix( type="Full", nrow=nv, ncol=ne, free=TRUE, values=.1, + label="e11", name="e" ) > > # Matrices to store the moderating a, c, and e Path Coefficients > modPathA <- mxMatrix( type='Full', nrow=nv, ncol=na, free=TRUE, values=.1, + label="aM11", name="aM" ) > modPathC <- mxMatrix( type='Full', nrow=nv, ncol=nc, free=TRUE, values=.1, + label="cM11", name="cM" ) > modPathE <- mxMatrix( type='Full', nrow=nv, ncol=ne, free=TRUE, values=.1, + label="eM11", name="eM" ) > > # Matrices for the expected means vector > intercept <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=.1, + label="interc", name="int" ) > PathM <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0, + label=c("m11"), name="m" ) > > # Matrices to compute the non-moderated A, C, and E Variance Components > covA <- mxAlgebra( expression=a %*% t(a), name="nA" ) > covC <- mxAlgebra( expression=c %*% t(c), name="nC" ) > covE <- mxAlgebra( expression=e %*% t(e), name="nE" ) > > # Algebra to compute the total variance per twin > covP <- mxAlgebra( expression=nA+nC+nE, name="nV" ) > > # Matrix for the moderator variable > mod <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), + name="D") > > # Matrices to compute the moderated A, C, and E variance components > covAmod <- mxAlgebra( expression=(a+ D%*%aM) %*% t(a+ D%*%aM), name="A" ) > covCmod <- mxAlgebra( expression=(c+ D%*%cM) %*% t(c+ D%*%cM), name="C" ) > covEmod <- mxAlgebra( expression=(e+ D%*%eM) %*% t(e+ D%*%eM), name="E" ) > > # Algebra to compute the total variance per twin > covPmod <- mxAlgebra( expression=A+C+E, name="V" ) > > # Algebra for the expected mean vector and expected variance/covariance matrices > covMZ <- mxAlgebra( expression= rbind ( cbind(A+C+E , A+C), + cbind(A+C, A+C+E)), name="expCovMZ" ) > covDZ <- mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C), + cbind(0.5%x%A+C, A+C+E)), name="expCovDZ" ) > wmod <- mxAlgebra( expression= m %*% D, name="DR" ) > meanG <- mxAlgebra( expression= cbind((int + DR),(int + DR)), name="expMeanG" ) > > # Data objects for Multiple Groups > dataMZ <- mxData( observed=mzData, type="raw" ) > dataDZ <- mxData( observed=dzData, type="raw" ) > > # Objective objects for Multiple Groups > objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanG", + dimnames=selVars ) > objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanG", + dimnames=selVars ) > fitMZ <- mxFitFunctionML() > fitDZ <- mxFitFunctionML() > > # Combine Groups > parsZ <- list( pathA, pathC, pathE, modPathA, modPathC, modPathE, PathM, + intercept) > modelMZ <- mxModel( parsZ, meanG, covAmod, covCmod, covEmod, covPmod, mod, + wmod, covMZ, dataMZ, objMZ, fitMZ, name="MZ" ) > modelDZ <- mxModel( parsZ, meanG, covAmod, covCmod, covEmod, covPmod, mod, + wmod, covDZ, dataDZ, objDZ, fitDZ, name="DZ" ) > minus2ll <- mxAlgebra( expression=MZ.objective+ DZ.objective, name="m2LL" ) > obj <- mxFitFunctionAlgebra( "m2LL" ) > ModAceModel <- mxModel( "ModACE", parsZ, modelMZ, modelDZ, minus2ll, obj ) > # Run Moderation ACE model > ModAceFit <- mxRun(ModAceModel) Running ModACE with 7 parameters Warning message: In model 'ModACE' 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) > ModAceSumm <- summary(ModAceFit) > ModAceSumm Summary of ModACE 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 Std.Error A 1 a11 a 1 1 4.098720e-01 0.048912110 2 c11 c 1 1 -1.935381e-05 0.077327481 3 e11 e 1 1 -7.790336e-01 0.016419892 4 aM11 aM 1 1 -5.948916e-02 0.005364976 5 cM11 cM 1 1 2.630094e-06 0.010116090 6 eM11 eM 1 1 6.845359e-02 0.001445987 7 interc int 1 1 1.989973e+00 0.006811057 observed statistics: 1622 estimated parameters: 7 degrees of freedom: 1615 fit value ( -2lnL units ): 102.7978 number of observations: 840 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: -3127.202 116.7978 NA BIC: -10771.646 149.9316 127.7019 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: 2016-06-01 16:23:16 Wall clock time (HH:MM:SS.hh): 00:00:01.24 optimizer: NPSOL OpenMx version number: 2.5.2 Need help? See help(mxSummary) > round(ModAceFit@output$estimate,4) a11 c11 e11 aM11 cM11 eM11 interc 0.4099 0.0000 -0.7790 -0.0595 0.0000 0.0685 1.9900 > > # Generate Output using Helper Functions > tableFitStatistics(ModAceFit) Name ep -2LL df AIC Model 1 : ModACE 7 102.8 1615 -3127.2 > > # Generate Table of Parameter Estimates using mxEval > pathEstimatesACE <- round(mxEval(cbind(a,c,e,aM,cM,eM), ModAceFit),4) > rownames(pathEstimatesACE) <- 'pathEstimates' > colnames(pathEstimatesACE) <- c('a','c','e','aM','cM','eM') > pathEstimatesACE a c e aM cM eM pathEstimates 0.4099 0 -0.779 -0.0595 0 0.0685 > varComponentsACE <- round(mxEval(cbind(A/V,C/V,E/V), ModAceFit@submodels$MZ),4) > #varComponentsACE <- round(mxEval(cbind(A,C,E), ModAceFit@submodels$MZ),4) > rownames(varComponentsACE) <- 'varComponents' > colnames(varComponentsACE) <- c('%a^2','%c^2','%e^2') > varComponentsACE %a^2 %c^2 %e^2 varComponents 0.0202 0 0.9798 > # Matrix to store the coefficient representing the main effect of Age on mean > PathM <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=0, + label=c("m11"), name="m" ) > > parsZ <- list( pathA, pathC, pathE, modPathA, modPathC, modPathE, PathM, + intercept) > modelMZ <- mxModel( parsZ, meanG, covAmod, covCmod, covEmod, covPmod, mod, + wmod, covMZ, dataMZ, objMZ, fitMZ, name="MZ" ) > modelDZ <- mxModel( parsZ, meanG, covAmod, covCmod, covEmod, covPmod, mod, + wmod, covDZ, dataDZ, objDZ, fitDZ, name="DZ" ) > minus2ll <- mxAlgebra( expression=MZ.objective+ DZ.objective, name="m2LL" ) > obj <- mxFitFunctionAlgebra( "m2LL" ) > ModAceMVModel <- mxModel( "ModACEMV", parsZ, modelMZ, modelDZ, minus2ll, obj ) > > ModAceMVFit <- mxRun(ModAceMVModel) Running ModACEMV with 8 parameters > ModAceMVSumm <- summary(ModAceMVFit) > ModAceMVSumm Summary of ModACEMV free parameters: name matrix row col Estimate Std.Error A 1 a11 a 1 1 0.029346277 0.159965906 2 c11 c 1 1 0.253653884 0.089363406 3 e11 e 1 1 -0.173837282 0.033940808 4 aM11 aM 1 1 -0.017110704 0.016842960 5 cM11 cM 1 1 -0.026102016 0.012496254 6 eM11 eM 1 1 0.001815360 0.003997925 7 m11 m 1 1 -0.004037045 0.004108584 8 interc int 1 1 2.023624488 0.033166406 observed statistics: 1622 estimated parameters: 8 degrees of freedom: 1614 fit value ( -2lnL units ): -666.4603 number of observations: 840 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: -3894.46 -650.4603 NA BIC: -11534.17 -612.5930 -637.9984 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: 2016-06-01 16:23:42 Wall clock time (HH:MM:SS.hh): 00:00:01.50 optimizer: NPSOL OpenMx version number: 2.5.2 Need help? See help(mxSummary) > tableFitStatistics(ModAceMVFit,ModAceFit) Name ep -2LL df AIC diffLL diffdf p Model 1 : ModACEMV 8 -666.46 1614 -3894.46 - - - Model 2 : ModACE 7 102.8 1615 -3127.2 769.26 1 0 > round(ModAceMVFit@output$estimate,4) a11 c11 e11 aM11 cM11 eM11 m11 interc 0.0293 0.2537 -0.1738 -0.0171 -0.0261 0.0018 -0.0040 2.0236