I am using a bivariate cholesky model for combined continuous and bivariate variables (bw) with Age effects on thresholds.
I am stuck with the error message in running model 'CholAce':A definition variable has been declared in model 'ACE' that does not contain a data set.
The code:
--------------------------------------------------------------------
CholAceModel <- mxModel("CholAce",
mxModel("ACE",
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,values=svPa, labels=labLower("a",nv),
lbound=lbPa, name="a"),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=svPa, labels=labLower("c",nv), lbound=lbPa, name="c" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=svPa, labels=labLower("e",nv), lbound=lbPa, name="e"),
# Matrices A, C, and E compute variance components
mxAlgebra( expression=a %% t(a), name="A" ),
mxAlgebra( expression=c %% t(c), name="C" ),
mxAlgebra( expression=e %% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( expression=A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( solve(sqrt(IV)), name="iSD"),
mxAlgebra(ACE.iSD %*% ACE.a,name="sta"),
mxAlgebra(ACE.iSD %*% ACE.c,name="stc"),
mxAlgebra(ACE.iSD %*% ACE.e,name="ste"),
mxAlgebra(ACE.A/ACE.V,name="h2"),
mxAlgebra(ACE.C/ACE.V,name="c2"),
mxAlgebra(ACE.E/ACE.V,name="e2"),
mxAlgebra(solve(sqrt(ACE.I*ACE.A)) %*% ACE.A %*% solve(sqrt(ACE.I*ACE.A)),name="ra"),
mxAlgebra(solve(sqrt(ACE.I*ACE.C)) %*% ACE.C %*% solve(sqrt(ACE.I*ACE.C)),name="rc"),
mxAlgebra(solve(sqrt(ACE.I*ACE.E)) %*% ACE.E %*% solve(sqrt(ACE.I*ACE.E)),name="re"),
# Define definition variables to hold the Covariates
mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age1"), name="Age1"),
mxMatrix( type="Full", nrow=1, ncol=1, free=F, labels=c("data.age2"), name="Age2"),
mxMatrix( type="Zero", nrow=1, ncol=ntv, name="M" ),
mxMatrix( type="Full", nrow=maxthresh, ncol=nv, free=T, values=.2, labels=LabCov, name="BageTH" ),
# Matrix & Algebra for expected thresholds vector
mxMatrix( type="Full", nrow=maxthresh, ncol=nvo, free=TRUE, values=c(-.1,-.1), lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThmz, name="ThMZ"),
mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Low" ),
mxAlgebra( expression= cbind(Low%%ThMZ+ BageTH%x%Age1,Low%%ThMZ+ BageTH%x%Age2 ), name="ExpThresMZ"),
mxMatrix( type="Full", nrow=maxthresh, ncol=nvo, free=TRUE, values=c(-.2,-.1), lbound=c(-3,rep(.001,ninc)), ubound=3, labels=LabThdz, name="ThDZ"),
mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Low" ),
mxAlgebra( expression= cbind(Low%%ThDZ+ BageTH%x%Age1,Low%%ThDZ+ BageTH%x%Age2 ), name="ExpThresDZ"),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra( expression= rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" ),
# Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
),
mxModel("MZ",
mxData( observed=mzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.M", dimnames=selVars,thresholds="ACE.ExpThresMZ", threshnames=c("bw1","bw2"))
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.M", dimnames=selVars,thresholds="ACE.ExpThresDZ", threshnames=c("bw1","bw2"))
),
mxAlgebra( expression=MZ.objective + DZ.objective, name="modelfit" ),
mxAlgebraObjective("modelfit"),
mxCI(c("ACE.sta","ACE.stc","ACE.ste","ACE.h2","ACE.c2","ACE.e2","ACE.ra","ACE.rc","ACE.re"))
)
CholAceModel_Fit <- mxRun(CholAceModel,intervals=TRUE)
CholAceModel_Summ <- summary(CholAceModel_Fit)
CholAceModel_Summ
Any help appreciated! Thanks.