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)