Hi,
In the univariate ACE modeling, I constrained the total variance of A, C and E as 1. As suggested by Rob, I put the following lines in the scripts:
unit <- mxMatrix( type="Unit", nrow=1, ncol=1, name="Unit") identifyingConstraint <- mxConstraint(VarP == Unit, name="ic")
I found that the results have been standardized that means the a2 is equils to the a2/(a2+c2+e2). However, the results is different from that if the total variance was not constrained as 1. Furthermore, if I did not constrain the total variance as 1, the scripts always echo red error 5 though the results seems correct. Hence I'm wondering that if to constrain the total variance as 1 is necessary and why? Many thanks~
In a univariate ACE model, there is no need to constrain the variance from A + C + E to be one, other than to standardize the result. If you are working with a standardized variable, then the results will be the same (i.e., the DV is already standardized, and so too, therefore, the model is as well. That would account for the lack of change.
Usually standardization on is handled not by a constraint, but by adding some algebras (either in the model or via mxEval() on the model output) to create Astd Cstd and Estd. That is computationally easier.
I am not sure what the line
unit identifyingConstraint
is doing: that's not valid R/OpenMx code.Not sure if your implementation of the constraint is correct. You might want to explore using the library
umx
, specifically theumxACE
function to get started. It's a well-documented one-line solution to implementing a uni or multivariate Cholesky and also works with umxSummary and plot() gives tabular and graphical output, including standardized, if you add std=TRUE.Many thanks, Tim. I'm sorry for that the codes I posted cannot been fully displayed. I'm wondering that why there is no need to contstrain the variance from A + C + E to be one in a univariate model, but to constrain the variance to be 1 in multivariate Cholesky model. If I constrain the variance to be 1 in a univariate model, will I get the wrong results? I found that if I do not constrain the variance to be 1 in a univariate model, it will echo "5: The Hessian at the solution does not appear to be convex (Mx status RED)". However, if I constrain the variance this error will disappear.
I edited your post so that the code is displayed with syntax highlighting. The main assignment operator in R,
<-
, gets parsed as a broken HTML tags unless it's inside a block of text formatted with the "code" tag. See here for HTML formatting tips.In general, constraining the variance to 1 is unnecessary in both cases. However, in the case of a common-pathway model--which was the context in which I suggested you add the MxConstraint in your previous thread-- there is a need for some kind of constraint to identify the scale of the latent phenotypic factor.
Yes, unless you standardized the phenotype ahead of time (or the phenotypic variance happens to be exactly 1.0).
I see. Many thanks, Rob. About the phenotype standardization, now my method is Z-transformation ((X-Mean)/SD). If the terminal of R reported "5: The Hessian at the solution does not appear to be convex (Mx status RED)" (though the results seems correct), will it affect the acquired results?
It certainly could. At the very least, it doesn't appear the optimizer has found a local minimum of the fitfunction. Which optimizer are you using? More importantly, why are you standardizing the phenotype in the first place?
Hi, Rob. Thanks for your kind help. About the optimizer, I used the default optimizer NPSOL. I have just change the optimizer from the NPSOL to the SLSQP. Although it still echo errors, there is not error 5 (Mx status RED) any more but only error 1 (Mx status GREEN). It seems that the optimizer SLSQP is better than the NPSOL, right? So there is two optimizer in OpenMx in total?
About the standardization of phenotype, for my sample size is limited (dozens of MZ and dozens of DZ) and the values of phenotype ranges between 10 and -10, hence I standardied the phenotype to acquire the normal distribution of data. The results of standardized data was better than that of raw data, at least in my data.
I'm confused here. The on-load default optimizer is SLSQP in every recent release of OpenMx. Furthermore, SLSQP doesn't report a status GREEN; only NPSOL does that.
If you installed OpenMx from our website (which it appears you did), there will be three general-purpose gradient-descent optimizers: SLSQP, NPSOL, and CSOLNP.
Z-transformation won't change the shape of your variables' distributions. It's a location-scale transformation, i.e. linear.
Thanks. What's the difference between the three optimizer? For different kind of data or condition, such like continuous and non-continuous variables?
In the scheme of things, the three gradient-descent optimizers are quite similar in most respects. CSOLNP is the weakest of the three, but it is under active development, and can sometimes be especially useful for certain analyses of ordinal data. Between NPSOL and SLSQP, which one is preferable is heavily context-dependent. There are "rules of thumb" for which seems typically do better in certain cases, but none in which I'm confident enough to endorse in writing.
You can find some details concerning those three optimizers in the help page for
mxOption()
andmxComputeGradientDescent()
. Note that both of those help pages have been updated since the release of version 2.6.9.Ok. Many thanks for your help to answer my questions~ :D
Besides, if there are some covariates, such like gender, age and IQ, etc. I can not included them all in the genetic model. Hence I directly regress them out from the raw data of phenotype. I have no idea if this method is right, especially for the dichotomic variable "gender".
Sure you can. It's certainly possible to have more than one definition variable in the model for the phenotypic mean. It's also possible to treat the covariates as endogenous, random regressors, and incorporate them into the model for the covariance. For that approach, you'll want to either (a) use
mxPath()
specification, (b) use helper functions from the 'umx' package, or (c) follow this article.It IS suboptimal. See the discussion in this other thread. Modeling the covariates as definition variables or as endogenous random regressors both have stronger statistical-theoretical justification, and if there are no missing scores on your covariates, should both give the same answer. The random-regressor approach is preferable if there are missing scores on your covariates.
I see. Many thanks, I'll try the script first.
Hi, Rob. I have included the Age and Gender in the Common pathway model as you suggested, attached please find the added codes about covariates:
-------Covariates-------------------------------------------------------------
Matrix for moderating/interacting variable
defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
labels=c("data.Sex_1","data.Sex_2"), name="Sex")
Matrices declared to store linear Coefficients for covariate
B_Sex <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values= .1, label=c("betaex_1","betaSex_2"), name="bSex" )
meanSex <- mxAlgebra( bSex%*%Sex, name="SexR")
age
defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
labels=c("data.Age_1","data.Age_2"), name="Age")
Matrices declared to store linear Coefficients for covariate
B_Age <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE,
values= .1, label="betaAge", name="bAge" )
meanAge <- mxAlgebra( bAge%*%Age, name="AgeR")
-----------------------------------------------------------------------------
Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanMZG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=0,
labels="meanMZLabs", name="expMeanMZG" ) # ntv = 6, there are 3 phenotypes in this model
meanDZG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=0,
labels="meanDZLabs", name="expMeanDZG" )
defs <- list( defSex, B_Sex, meanSex, defAge, B_Age, meanAge)
meanMZ <- mxAlgebra( expression = expMeanMZG + SexR + AgeR, name = "expMeanMZ") #with covariates
meanDZ <- mxAlgebra( expression = expMeanDZG + SexR + AgeR, name = "expMeanDZ") #with covariates
However, it reported that: The following error occurred while evaluating the subexpression 'MZ.expMeanMZG + MZ.AgeR' during the evaluation of 'MZ.expMeanMZ' in model 'ComACE' : non-conformable arrays. There are three phenotypes and two latent factors in this common pathway model. Did I miss something? BTW, these same codes can work in a univariate ACE model estimation.
There are numerous problems with that syntax. Try this instead:
Many thanks, Rob. It's ok now. I think I added "" to the lables of meanMZG and meanDZG that caused the non-compartibale arrays.