Univariate Heterogeneity model with covariates (including DZO)
YiTan
Joined: 10/08/2014
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
search the script for the name this message says is missing
best, tim
Log in or register to post comments
In reply to search the script for the name this message says is missing by tbates
means vector of different length from 'dimnames' argument
Error: The expected means vector associated with the expectation function in model 'DZo' is not of the same length as the 'dimnames' argument provided by the expectation function. The 'dimnames' argument is of length 2 and the expected means vector has 4 columns.
And so I changed my meanGfm from:
meanGfm <- mxAlgebra(expression= cbind(expMeanGf, expMeanGm), name="expMeanGfm")
to
meanGfm <- mxAlgebra(expression= cbind(expMeanGf[,1], expMeanGm[,1]), name="expMeanGfm")
Is that conceptually correct?
Thanks and best regards,
Yi Ting
Log in or register to post comments