You are here

Constrain the total variance as 1?

17 posts / 0 new
Last post
Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
Constrain the total variance as 1?

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~

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
A+C+E=1, standardized raw data, umxACE

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 the umxACE 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.

Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
Why do not constrain

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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
HTML; model identification
I'm sorry for that the codes I posted cannot been fully displayed.

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.

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.

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.

If I constrain the variance to be 1 in a univariate model, will I get the wrong results?

Yes, unless you standardized the phenotype ahead of time (or the phenotypic variance happens to be exactly 1.0).

Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
RED Erorr 5?

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?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
why standardizing?

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?

Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
optimizer NPSOL vs. SLSQP

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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
optimizers & data
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).

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.

It seems that the optimizer SLSQP is better than the NPSOL, right? So there is two optimizer in OpenMx in total?

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.

I standardied the phenotype to acquire the normal distribution of data.

Z-transformation won't change the shape of your variables' distributions. It's a location-scale transformation, i.e. linear.

Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
appropriate optimizer

Thanks. What's the difference between the three optimizer? For different kind of data or condition, such like continuous and non-continuous variables?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
"best" optimizer depends on optimization problem

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() and mxComputeGradientDescent(). Note that both of those help pages have been updated since the release of version 2.6.9.

Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
Ok

Ok. Many thanks for your help to answer my questions~ :D

Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
Covariate regression

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".

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
covariates
Besides, if there are some covariates, such like gender, age and IQ, etc. I can not included them all in the genetic model.

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.

Hence I directly regress them out from the raw data of phenotype. I have no idea if this method is right

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.

Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
Thanks

I see. Many thanks, I'll try the script first.

Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
About the Covariates in Common pathway model

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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
rewrite

There are numerous problems with that syntax. Try this instead:

# -------Covariates-------------------------------------------------------------
 
# Matrix for sex "definition variable":
defSex <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
labels=c("data.Sex_1","data.Sex_2"), name="Sex")
 
# Each of the three phenotypes will need its own regression onto sex:
B_Sex <- mxMatrix( type="Full", nrow=1, ncol=3, free=TRUE,
values= .1, label=c("betaSex_1","betaSex_2","betaSex_3"), name="bSex" )
 
#Note the Kronecker product operator; SexR will be a 1x6 matrix:
meanSex <- mxAlgebra( Sex %x% bSex, name="SexR")
 
#age
defAge <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE,
labels=c("data.Age_1","data.Age_2"), name="Age")
 
# Each of the three phenotypes will need its own regression onto age:
B_Age <- mxMatrix( type="Full", nrow=1, ncol=3, free=TRUE,
values= .1, labels=c("betaAge_1","betaAge_2","betaAge_3") name="bAge" )
 
meanAge <- mxAlgebra( Age %x% bAge, name="AgeR")
# -----------------------------------------------------------------------------
 
# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
# Do you actually want each regression intercept to be different for twin #1 and twin #2?  
# I doubt you do, but that's what your syntax was doing before...
meanMZG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=0,
labels=rep(meanMZLabs[1:3],2), name="expMeanMZG" ) # ntv = 6, there are 3 phenotypes in this model
meanDZG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=0,
labels=rep(meanDZLabs[1:3],2), 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
Liz's picture
Liz
Offline
Joined: 04/21/2016 - 00:27
Thanks

Many thanks, Rob. It's ok now. I think I added "" to the lables of meanMZG and meanDZG that caused the non-compartibale arrays.