Hi!
I'm new to twin data analysis and would like to get help on how to adjust for covariates in sex-lim ACE model. I tried to modify the script twinHet5AceCon.R by adding Age(continuous variable) and Region(binary variable)as covariates in the mean structure of the genetic model but I keep getting an error: Unknown reference 'meanGf' detected in the entity 'expMeanGf' in model 'MZf'.
I guess I may did something wrong somewhere near "expMeanGfmCo <- mxAlgebra( expression= meanGfm +Betaall %*% meanGfmCo , name="expMeanGfm" ) ".
The script is below. Any suggestions would be very much appreciated.
Thanks,
Sophie
————————————————————————————————————————-——————--
Load Data
TwinData<-read.csv(************)
describe(TwinData)
Select Variables for Analysis
Vars <- 'w'
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- c('lnBMI','lnBMI2')
Select Data for Analysis
mzfData <- subset(TwinData, zyg_sex==2, c(lnBMI,lnBMI2,AGE,AGE2,region,region2))
dzfData <- subset(TwinData, zyg_sex==4, c(lnBMI,lnBMI2,AGE,AGE2,region,region2))
mzmData <- subset(TwinData, zyg_sex==1, c(lnBMI,lnBMI2,AGE,AGE2,region,region2))
dzmData <- subset(TwinData, zyg_sex==3, c(lnBMI,lnBMI2,AGE,AGE2,region,region2))
dzoData <- subset(TwinData, zyg_sex==5, c(lnBMI,lnBMI2,AGE,AGE2,region,region2))
svMe <- 20 # start value for means
svPa <- .6 # start value for path
----------------------------------------------------------------
PREPARE MODEL
General non-scalar ACE Model
Matrices declared to store a, c, and e Path Coefficients
pathAf <- mxMatrix( "Lower", nrow=1, ncol=1, free=TRUE,
values=.6, label="af11", name="af" )
pathCf <- mxMatrix( "Lower", nrow=1, ncol=1, free=TRUE,
values=.6, label="cf11", name="cf" )
pathEf <- mxMatrix( "Lower", nrow=1, ncol=1, free=TRUE,
values=.6, label="ef11", name="ef" )
pathAm <- mxMatrix( "Lower", nrow=1, ncol=1, free=TRUE,
values=.6, label="am11", name="am" )
pathCm <- mxMatrix( "Lower", nrow=1, ncol=1, free=TRUE,
values=.6, label="cm11", name="cm" )
pathEm <- mxMatrix( "Lower", nrow=1, ncol=1, free=TRUE,
values=.6, label="em11", name="em" )
pathRg <- mxMatrix( "Lower", nrow=1, ncol=1, free=TRUE,
values=0.5, label="rg11", name="rg", ubound=1, lbound=0 )
Matrices generated to hold A, C, and E computed Variance Components
covAf <- mxAlgebra( af %% t(af), name="Af" )
covCf <- mxAlgebra( cf %% t(cf), name="Cf" )
covEf <- mxAlgebra( ef %*% t(ef), name="Ef" )
covAm <- mxAlgebra( am %% t(am), name="Am" )
covCm <- mxAlgebra( cm %% t(cm), name="Cm" )
covEm <- mxAlgebra( em %*% t(em), name="Em" )
Algebra to compute total variances and standard deviations
(diagonal only)
covPf <- mxAlgebra( Af+Cf+Ef, name="Vf" )
covPm <- mxAlgebra( Am+Cm+Em, name="Vm" )
Algebras generated to hold Parameter Estimates and Derived
Variance Components
colVarsZf <- c('Af','Cf','Ef','SAf','SCf','SEf')
estVarsZf <- mxAlgebra( cbind(Af,Cf,Ef,Af/Vf,Cf/Vf,Ef/Vf),
name="VarsZf", dimnames=list(NULL,colVarsZf))
colVarsZm <- c('Am','Cm','Em','SAm','SCm','SEm')
estVarsZm <- mxAlgebra( cbind(Am,Cm,Em,Am/Vm,Cm/Vm,Em/Vm),
name="VarsZm", dimnames=list(NULL,colVarsZm))
Declare a matrix for the definition variable regression parameters, called beta
Betaall<- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= 0, label=c("betaAge","betaRegion"), name="beta")
Algebra for expected Mean and Variance/Covariance Matrices in
MZ & DZ twins
meanGf <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE,
values=svMe, labels="xbmi", name="meanf" )
meanGfCo <- mxMatrix( type="Full", nrow=2, ncol=2, free=F,
label=c("data.AGE","data.region","data.AGE2","data.region2"), name="meanfCo" )
expMeanGfCo <- mxAlgebra( expression= meanGf+Betaall %*% meanGfCo,
name="expMeanGf" )
meanGm <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE,
values=svMe, labels="xbmi", name="meanm" )
meanGmCo <- mxMatrix( type="Full", nrow=2, ncol=2, free=F,
label=c("data.AGE","data.region","data.AGE2","data.region2"),name="meanmCo" )
expMeanGmCo <- mxAlgebra( expression= meanGm +Betaall %*% meanGmCo ,
name="expMeanGm" )
meanGfm <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE,
values=svMe, labels="xbmi", name="meanfm" )
meanGfmCo <- mxMatrix( type="Full", nrow=2, ncol=2, free=F,
label=c("data.AGE","data.region","data.AGE2","data.region2"), name="meanfmCo" )
expMeanGfmCo <- mxAlgebra( expression= meanGfm +Betaall %*% meanGfmCo ,
name="expMeanGfm" )
covMZf <- mxAlgebra(rbind( cbind(Af+Cf+Ef , Af+Cf),
cbind(Af+Cf , Af+Cf+Ef)),
name="expCovMZf" )
covDZf <- mxAlgebra(rbind( cbind(Af+Cf+Ef , 0.5%x%Af+Cf),
cbind(0.5%x%Af+Cf , Af+Cf+Ef)),
name="expCovDZf" )
covMZm <- mxAlgebra(rbind( cbind(Am+Cm+Em , Am+Cm),
cbind(Am+Cm , Am+Cm+Em)),
name="expCovMZm" )
covDZm <- mxAlgebra(rbind( cbind(Am+Cm+Em , 0.5%x%Am+Cm),
cbind(0.5%x%Am+Cm , Am+Cm+Em)),
name="expCovDZm" )
covDZo <- mxAlgebra(rbind( cbind(Af+Cf+Ef, 0.5%%rg%x%(af%%t(am))+cf%%t(cm)),
cbind(0.5%%rg%x%(am%%t(af))+cm%%t(cf) , Am+Cm+Em)),
name="expCovDZo" )
Data objects for Multiple Groups
dataMZf <- mxData( observed=mzfData, type="raw" )
dataDZf <- mxData( observed=dzfData, type="raw" )
dataMZm <- mxData( observed=mzmData, type="raw" )
dataDZm <- mxData( observed=dzmData, type="raw" )
dataDZo <- mxData( observed=dzoData, type="raw" )
Objective objects for Multiple Groups
objMZf <- mxFIMLObjective( covariance="expCovMZf",
means="expMeanGfCo", dimnames=selVars )
objDZf <- mxFIMLObjective( covariance="expCovDZf",
means="expMeanGfCo", dimnames=selVars )
objMZm <- mxFIMLObjective( covariance="expCovMZm",
means="expMeanGmCo", dimnames=selVars )
objDZm <- mxFIMLObjective( covariance="expCovDZm",
means="expMeanGmCo", dimnames=selVars )
objDZo <- mxFIMLObjective( covariance="expCovDZo",
means="expMeanGfmCo", dimnames=selVars )
Combine Groups
parsZf <- list( pathAf, pathCf, pathEf, covAf, covCf, covEf,
covPf, estVarsZf)
parsZm <- list( pathAm, pathCm, pathEm, covAm, covCm, covEm,
covPm, estVarsZm)
modelMZf <- mxModel( parsZf, expMeanGfCo, meanGfCo, covMZf, dataMZf,
objMZf, name="MZf" )
modelDZf <- mxModel( parsZf, expMeanGfCo, meanGfCo, covDZf, dataDZf,
objDZf, name="DZf" )
modelMZm <- mxModel( parsZm, expMeanGmCo, meanGmCo, covMZm, dataMZm,
objMZm, name="MZm" )
modelDZm <- mxModel( parsZm, expMeanGmCo, meanGmCo, covDZm, dataDZm,
objDZm, name="DZm" )
modelDZo <- mxModel( parsZf, pathRg, parsZm, expMeanGfmCo, meanGfmCo,
covDZo, dataDZo, objDZo, name="DZo" )
minus2ll <- mxAlgebra( MZf.objective+ DZf.objective+ MZm.objective+
DZm.objective+ DZo.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
conf <- mxCI(c('QualACE.VarsZf', 'QualACE.VarsZm'))
QualAceModel <- mxModel( "QualACE", parsZf, parsZm, modelMZf,
modelDZf, modelMZm, modelDZm, modelDZo, minus2ll, obj, conf)
ProTip: When including lots of syntax, either append it to your post as an attachment, or enclose it in the text of your post with the appropriate HTML "code" and "\code" tags.
To answer your question, try changing lines 74-95 of your syntax to the following:
MxAlgebras only look at the
name
s of MxMatrices and other MxAlgebras included in the MxModel; they don't "know" anything about the R symbols for the MxMatrix and MxAlgebra objects sitting in the global environment (a.k.a. your workspace).