You are here

Univariate Heterogeneity model with covariates (including DZO)

3 posts / 0 new
Last post
YiTan's picture
Offline
Joined: 10/08/2014 - 05:28
Univariate Heterogeneity model with covariates (including DZO)

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

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
search the script for the name this message says is missing

I think your modelDZo needs meanGf and meanGm as well as meanGfm included in it (because the meanGfm algebra refers to these).

best, tim

YiTan's picture
Offline
Joined: 10/08/2014 - 05:28
means vector of different length from 'dimnames' argument

Thanks Tim, that was very helpful, I included meanGf and meanGm in modelDZo and I got a new error message:

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