You are here

Adding covariates to sex-lim ACE model

2 posts / 0 new
Last post
liuqq1989's picture
Offline
Joined: 07/09/2015 - 12:37
Adding covariates to sex-lim ACE model

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)

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
OpenMx looks at MxAlgebra names, not R symbols

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:

# 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= meanf+Betaall %*% meanfCo,
     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= meanm +Betaall %*% meanmCo ,
     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= meanfm +Betaall %*% meanfmCo ,
     name="expMeanGfm" )

MxAlgebras only look at the names 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).