# threshold ACE model with a covariate

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(I*V)), 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")

)

## Oh, it is a bivariate

Log in or register to post comments

## What error are you getting?

I *think* that you've failed to properly identify the scale of the categorical variables. You end up defining the expected covariance in terms of A, C and E, which are all completely free. You have to identify the scale of the latent continuous variable underlying each categorical variable, which is most often done by constraining both the mean (typically zero, as you have done) and the variance (either total or residual/E) to a non-zero constant.

I'll add that the lower matrix trick for thresholds (deviations %*% lower triangular = observed thresholds) is (a) of debatable usefulness in OpenMx, as shown by the dev-team debates on this forum, and (b) only useful if you set lower bounds of zero for all but the first row of free parameters. That algebra is there to make sure that the thresholds are strictly increasing, which is easiest to do if you use the above algebra and say that all deviations must be positive. For example, if free parameters a, b and c represent the three free parameters estimated in a thresholds deviation matrix for a four category ordinal variable, the observed threshold between the first and second categories is 'a', between the second and third is 'a+b' and between the third and fourth is 'a+b+c'. You have to constrain b and c to be positive for the algebra to add anything.

For future reference, please submit code as attachments (R or text files). It makes it easier for users and developers to run your code (prevents formatting issues, allows for the inclusion of line breaks and white space, lets us refer to line numbers, etc).

Log in or register to post comments

In reply to What error are you getting? by Ryne

## the approved script

I also wanted to ask whether the covariates (e.g., age) could be fitted in saturated model and common/independent pathway model. If they can, how? Is the method the same as that used in the Cholesky model?

Log in or register to post comments

In reply to the approved script by userzht

## Can you tell me what you

What error are you getting? If you can't describe the error, then can you submit some data (either real, simulated or from the fakeData function on our wiki) so we can run the code? We need to know more than "it's not working" to diagnose the problem.

You can include covariates as definition variables in just about any model. Whatever method you use for Cholesky models should be applicable to other models. What parameters in your model to you think vary with or depend on age?

Log in or register to post comments

In reply to Can you tell me what you by Ryne

## actually, it works after I changed the nrow and ncol

The reason for adding age in the model is that I thought age may have an effect on the genetic correlation between smoke and drink. The script works after I changed the 'nrow' and 'ncol' when defining the covariate-age, specifically, 'betaAge', 'DefVarsMZ', and 'DefVarsDZ'. But the h2, c2, and e2 changed a lot comparing with the model without a covariate, for example, h2 from 0.50 to 0.10. So please diagnosis the script for me. Thank you very much!

Log in or register to post comments

In reply to actually, it works after I changed the nrow and ncol by userzht

## No DZ?

Log in or register to post comments