Login

News

threshold ACE model with a covariate

Posted on
No user picture. userzht Joined: 02/19/2011
Forums
Categorical Outcomes
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(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")
)

Replied on Mon, 09/26/2011 - 13:46
Picture of user. Ryne Joined: 07/31/2009

What error are you getting? Without some data or indication of what error you're getting, its tough to help, but I'll hazard a guess.

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
Replied on Tue, 10/11/2011 - 08:56
No user picture. userzht Joined: 02/19/2011

Thanks for your reply. I tried to revise my script which I attached here. Could you please turn over it?

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
Replied on Tue, 10/11/2011 - 10:52
Picture of user. Ryne Joined: 07/31/2009

Can you tell me what you changed, other than adding some new variables at the top of the code? I can't tell if you changed the scaling of your variables or not.

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
Replied on Fri, 10/14/2011 - 23:07
No user picture. userzht Joined: 02/19/2011

Sorry for the delay.

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
Replied on Wed, 10/19/2011 - 08:28
Picture of user. neale Joined: 07/31/2009

Your example file has only zygosity = 1 which would cause the model to be under identified and all sorts of things could then go wrong. If in practice your data do have DZs, then age might reduce h2 if: i) age is correlated with the phenotype; and ii) age is more highly correlated in MZ pairs than in DZ pairs. Usually if the age correlation is the same, then c^2 is reduced by regressing age out of the phenotype.
Log in or register to post comments
  Login

News