Hi all,
I tried to modify the heterogeneity model for MZM, DZM, MZF, DZF and DZO groups to include covariates. I got the following error message when I tried to run my code:
Running QualACE Error: Unknown reference 'expMeanGf' detected in the entity 'expMeanGfm' in model 'DZo'
May I know what changes should I make to my code to run it successfully? My code is as follows:
require(OpenMx) require(psych) # ------------------------------------------------------------------------------ # PREPARE DATA # Select Variables for Analysis Vars <- c('yrsTraining','Age','Sex','Var1') nvInit <- 4 nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(Vars,c(rep(1,nvInit),rep(2,nvInit)),sep="") #eg c('ht1','wt1,'ht2','wt2') #Select Data for Analysis mzfData <- subset(twinData2015, ZygosityGrp==1, selVars) dzfData <- subset(twinData2015, ZygosityGrp==2, selVars) mzmData <- subset(twinData2015, ZygosityGrp==3, selVars) dzmData <- subset(twinData2015, ZygosityGrp==4, selVars) dzoData <- subset(twinData2015, ZygosityGrp==5, selVars) colnames(mzfData)=colnames(dzfData)=colnames(mzmData)=colnames(dzmData)=colnames(dzoData)=c('yrsTraining1','Age1','Sex1','t1','yrsTraining2','Age2','Sex2','t2') # Generate Descriptive Statistics colMeans(mzfData,na.rm=TRUE) colMeans(dzfData,na.rm=TRUE) colMeans(mzmData,na.rm=TRUE) colMeans(dzmData,na.rm=TRUE) colMeans(dzoData,na.rm=TRUE) cov(mzfData,use="complete") cov(dzfData,use="complete") cov(mzmData,use="complete") cov(dzmData,use="complete") cov(dzoData,use="complete") mzfData$t1[is.na(mzfData$yrsTraining1)] <- NA mzfData$yrsTraining1[is.na(mzfData$yrsTraining1)] <- -999 mzfData$t2[is.na(mzfData$yrsTraining2)] <- NA mzfData$yrsTraining2[is.na(mzfData$yrsTraining2)] <- -999 mzmData$t1[is.na(mzmData$yrsTraining1)] <- NA mzmData$yrsTraining1[is.na(mzmData$yrsTraining1)] <- -999 mzmData$t2[is.na(mzmData$yrsTraining2)] <- NA mzmData$yrsTraining2[is.na(mzmData$yrsTraining2)] <- -999 dzfData$t1[is.na(dzfData$yrsTraining1)] <- NA dzfData$yrsTraining1[is.na(dzfData$yrsTraining1)] <- -999 dzfData$t2[is.na(dzfData$yrsTraining2)] <- NA dzfData$yrsTraining2[is.na(dzfData$yrsTraining2)] <- -999 dzmData$t1[is.na(dzmData$yrsTraining1)] <- NA dzmData$yrsTraining1[is.na(dzmData$yrsTraining1)] <- -999 dzmData$t2[is.na(dzmData$yrsTraining2)] <- NA dzmData$yrsTraining2[is.na(dzmData$yrsTraining2)] <- -999 # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Matrices declared to store a, c, and e Path Coefficients pathAf <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="af11", name="af" ) pathCf <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="cf11", name="cf" ) pathEf <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="ef11", name="ef" ) pathAm <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="am11", name="am" ) pathCm <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="cm11", name="cm" ) pathEm <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="em11", name="em" ) pathRg <- mxMatrix( "Lower", nrow=1, ncol=1, free=TRUE, values=1, label="rg11", name="rg", ubound=1, lbound=0 ) # Matrices generated to hold A, C, and E computed Variance Components covAf <- mxAlgebra( expression=af %*% t(af), name="Af" ) covCf <- mxAlgebra( expression=cf %*% t(cf), name="Cf" ) covEf <- mxAlgebra( expression=ef %*% t(ef), name="Ef" ) covAm <- mxAlgebra( expression=am %*% t(am), name="Am" ) covCm <- mxAlgebra( expression=cm %*% t(cm), name="Cm" ) covEm <- mxAlgebra( expression=em %*% t(em), name="Em" ) # Algebra to compute total variances and standard deviations (diagonal only) covPf <- mxAlgebra( expression=Af+Cf+Ef, name="Vf" ) covPm <- mxAlgebra( expression=Am+Cm+Em, name="Vm" ) # Algebras generated to hold Parameter Estimates and Derived Variance Components colVarsZf <- rep(c('Af','Cf','Ef','SAf','SCf','SEf'),each=nv) estVarsZf <- mxAlgebra( expression=cbind(Af,Cf,Ef,Af/Vf,Cf/Vf,Ef/Vf), name="VarsZf", dimnames=list(NULL,colVarsZf)) colVarsZm <- rep(c('Am','Cm','Em','SAm','SCm','SEm'),each=nv) estVarsZm <- mxAlgebra( expression=cbind(Am,Cm,Em,Am/Vm,Cm/Vm,Em/Vm), name="VarsZm", dimnames=list(NULL,colVarsZm)) # Matrix for expected Mean intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values= 5, label="mean", name="Mean" ) # Matrix for moderating/interacting variable defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.Sex1","data.Sex2"), name="Sex") # Matrices declared to store linear Coefficients for covariate B_Sex <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= .01, label="betaSex", name="bSex" ) meanSex <- mxAlgebra( bSex%*%Sex, name="SexR") #age defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.Age1","data.Age2"), name="Age") # Matrices declared to store linear Coefficients for covariate B_Age <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= .01, label="betaAge", name="bAge" ) meanAge <- mxAlgebra( bAge%*%Age, name="AgeR") #yrsTraining defYTrain <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.yrsTraining1","data.yrsTraining2"), name="YTrain") # Matrices declared to store linear Coefficients for covariate B_YTrain <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= .01, label="betaYTrain", name="bYTrain" ) meanYTrain <- mxAlgebra( bYTrain%*%YTrain, name="YTrainR") #*************************************************************** #*************************************************************** # # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins # expMean <- mxAlgebra( Mean + YTrainR + AgeR + SexR , name="expMean") meanGf <- mxAlgebra( Mean + YTrainR + AgeR + SexR , name="expMeanGf") meanGm <- mxAlgebra( Mean + YTrainR + AgeR + SexR , name="expMeanGm") meanGfm <- mxAlgebra(expression= cbind(expMeanGf, expMeanGm), name="expMeanGfm") defs <- list( intercept, defYTrain, B_YTrain, meanYTrain, defSex, B_Sex, meanSex, defAge, B_Age, meanAge) covMZf <- mxAlgebra( expression= rbind( cbind(Vf, Af+Cf), cbind(Af+Cf, Vf)), name="expCovMZf" ) covDZf <- mxAlgebra( expression= rbind( cbind(Vf, 0.5%x%Af+Cf), cbind(0.5%x%Af+Cf, Vf)), name="expCovDZf" ) covMZm <- mxAlgebra( expression= rbind( cbind(Vm, Am+Cm), cbind(Am+Cm, Vm)), name="expCovMZm" ) covDZm <- mxAlgebra( expression= rbind( cbind(Vm, 0.5%x%Am+Cm), cbind(0.5%x%Am+Cm, Vm)), name="expCovDZm" ) CVfm <- mxAlgebra( expression= 0.5%*%rg%x%(af%*%t(am))+cf%*%t(cm), name="CVfm" ) CVmf <- mxAlgebra( expression= 0.5%*%rg%x%(am%*%t(af))+cm%*%t(cf), name="CVmf" ) covDZo <- mxAlgebra( expression= rbind( cbind(Vf, CVfm), cbind(CVmf, Vm)), name="expCovDZo" ) # Data objects for Multiple Groups dataMZf <- mxData(subset(mzfData, (!is.na(mzfData$t1)|!is.na(mzfData$t2))), "raw") dataDZf <- mxData(subset(dzfData, (!is.na(dzfData$t1)|!is.na(dzfData$t2))), "raw") dataMZm <- mxData(subset(mzmData, (!is.na(mzmData$t1)|!is.na(mzmData$t2))), "raw") dataDZm <- mxData(subset(dzmData, (!is.na(dzmData$t1)|!is.na(dzmData$t2))), "raw") dataDZo <- mxData(subset(dzoData, (!is.na(dzoData$t1)|!is.na(dzoData$t2))), "raw") # Expectation objects for Multiple Groups expMZf <- mxExpectationNormal( covariance="expCovMZf", means="expMeanGf", dimnames=c('t1','t2') ) expDZf <- mxExpectationNormal( covariance="expCovDZf", means="expMeanGf", dimnames=c('t1','t2') ) expMZm <- mxExpectationNormal( covariance="expCovMZm", means="expMeanGm", dimnames=c('t1','t2') ) expDZm <- mxExpectationNormal( covariance="expCovDZm", means="expMeanGm", dimnames=c('t1','t2') ) expDZo <- mxExpectationNormal( covariance="expCovDZo", means="expMeanGfm", dimnames=c('t1','t2') ) funML <- mxFitFunctionML() # Combine Groups parsZf <- list( pathAf, pathCf, pathEf, covAf, covCf, covEf, covPf, estVarsZf ) parsZm <- list( pathAm, pathCm, pathEm, covAm, covCm, covEm, covPm, estVarsZm ) parsZfm <- list( pathRg, CVfm, CVmf) modelMZf <- mxModel( parsZf, defs, meanGf, covMZf, dataMZf, expMZf, funML, name="MZf" ) modelDZf <- mxModel( parsZf, defs, meanGf, covDZf, dataDZf, expDZf, funML, name="DZf" ) modelMZm <- mxModel( parsZm, defs, meanGm, covMZm, dataMZm, expMZm, funML, name="MZm" ) modelDZm <- mxModel( parsZm, defs, meanGm, covDZm, dataDZm, expDZm, funML, name="DZm" ) modelDZo <- mxModel( parsZf, parsZm, parsZfm, meanGfm, covDZo, dataDZo, expDZo, funML, name="DZo" ) fitML <- mxFitFunctionMultigroup(c("MZf.fitfunction","DZf.fitfunction","MZm.fitfunction","DZm.fitfunction","DZo.fitfunction")) QualAceModel <- mxModel( "QualACE", parsZf, parsZm, modelMZf, modelDZf, modelMZm, modelDZm, modelDZo, fitML, mxCI(c("VarsZf[1,4:6]","VarsZm[1,4:6]")) ) # ------------------------------------------------------------------------------ # RUN MODEL # Run Qualitative Sex Differences ACE model QualAceFit <- mxRun(QualAceModel,intervals=T) QualAceSum <- summary(QualAceFit) QualAceSum$pa round(QualAceFit@output$estimate,4) round(cbind(QualAceFit$VarsZf@result,QualAceFit$VarsZm@result),4) mxCompare(QualAceFit) # ------------------------------------------------------------------------------ # FIT SUBMODELS # Run Quantitative non-scalar Sex Differences ACE model QuanAceModel <- mxModel( QualAceFit, name="QuanACE" ) QuanAceModel <- omxSetParameters( QuanAceModel, labels="rg11", free=FALSE, values=1 ) QuanAceFit <- mxRun(QuanAceModel,intervals=T) round(QuanAceFit@output$estimate,4) round(cbind(QuanAceFit$VarsZf@result,QuanAceFit$VarsZm@result),4) # Run Homogeneity ACE model HomoAceModel <- mxModel( QuanAceFit, name="HomoAce" ) HomoAceModel <- omxSetParameters( HomoAceModel, labels=c("af11","am11"), free=TRUE, values=.6, newlabels='a11' ) HomoAceModel <- omxSetParameters( HomoAceModel, labels=c("cf11","cm11"), free=TRUE, values=.6, newlabels='c11' ) HomoAceModel <- omxSetParameters( HomoAceModel, labels=c("ef11","em11"), free=TRUE, values=.6, newlabels='e11' ) HomoAceFit <- mxRun(HomoAceModel,intervals=T) round(HomoAceFit@output$estimate,4) round(cbind(HomoAceFit$VarsZf@result,HomoAceFit$VarsZm@result),4) # ------------------------------------------------------------------------------ # Print Comparative Fit Statistics AceNested <- list(QuanAceFit,HomoAceFit) mxCompare(QualAceFit,AceNested) round(QualAceFit@output$estimate,4) round(QuanAceFit@output$estimate,4) round(HomoAceFit@output$estimate,4) round(rbind(cbind(QualAceFit$VarsZf@result,QualAceFit$VarsZm@result), cbind(QuanAceFit$VarsZf@result,QuanAceFit$VarsZm@result), cbind(HomoAceFit$VarsZf@result,HomoAceFit$VarsZm@result)),4)
Thanks and best regards,
Yi Ting
I think your modelDZo needs meanGf and meanGm as well as meanGfm included in it (because the meanGfm algebra refers to these).
best, tim
Thanks Tim, that was very helpful, I included meanGf and meanGm in modelDZo and I got a new error message:
And so I changed my meanGfm from:
to
Is that conceptually correct?
Thanks and best regards,
Yi Ting