I want to fit a threshold ACE model controlling for age. The following is my script. Can anyone tell me what the problem is?
ACE_Cholesky_Model <- mxModel("ACE_Cholesky",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("a11", "a21", "a22"), name="a" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("c11", "c21", "c22"), name="c" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("e11", "e21", "e22"), 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="rg"),
## Note that the rest of the mxModel statements do not change for bivariate/multivariate case
# Matrix & Algebra for expected means vector
mxMatrix( type="Full", nrow=1, ncol=ntv, free=FALSE, values= 0, name="expMeanO"),
# Declare a matrix for the definition variable regression parameters, called beta
mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=0, label=c("beta"), name="betaAge"),
# Algebra for making the means a function of the definition variables age
mxMatrix( type="Full", nrow=1, ncol=2, free=F, values=c(mzData$age1,mzData$age2),label=c("mzData.age1","mzData.age2"), name="DefVarsMZ"),
mxAlgebra( expression=expMeanO + betaAge %% DefVarsMZ, name="expMeanMZ"),
mxMatrix( type="Full", nrow=1, ncol=2, free=F, values=c(dzData$age1,dzData$age2),label=c("dzData.age1","dzData.age2"), name="DefVarsDZ"),
mxAlgebra( expression=expMeanO + betaAge %% DefVarsDZ, name="expMeanDZ"),
# Matrix & Algebra for expected thresholds vector
mxMatrix(type="Full", nrow=maxthresh, ncol=ntv, free=TRUE, values=c(0.1,-.3,0.1,-.3), lbound=-4,name="cutpoints"),
if you have more than 1 threshold you will need to bound the second threshold to stop the thresholds running backwards
ie for 2 thresholds the bounds might look like this
lbound=c(-4,0.001,-4,0.001,-4,0.001,-4,0.001)
mxMatrix( type="Lower", nrow=maxthresh, ncol=maxthresh, free=FALSE, values=1, name="Increment" ),
mxAlgebra( expression= (Increment %*% cutpoints), dimnames=list(thresh,selVars), name="expThresh"),
# 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.expMeanMZ", dimnames=selVars,
thresholds="ACE.expThresh")
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMeanDZ", dimnames=selVars,
thresholds="ACE.expThresh")
),
mxAlgebra( expression=MZ.objective + DZ.objective, name="modelfit" ),
mxAlgebraObjective("modelfit")
)