Login

News

OpenMx sensitivity to start values in a,c,e?

Posted on
Picture of user. tbates Joined: 07/31/2009
Forums
Behavioral Genetics Models
running a tri-bivariate ACE model, and the results can be all over the show (with no convergence warnings), depending on the start values: for instance C wandering from .998 to .054 for one variable.
Any thoughts?
Here are the results, and script fragment

Results starting at a, c, and e at .6

               A1     A2    A3    C1    C2    C3    E1    E2    E3
var1     0.284     NA    NA 0.282    NA    NA 0.434    NA    NA
var2    -0.292  0.027    NA 0.998 0.434    NA 0.294 0.539    NA
var3     0.224 -0.028 0.206 0.366 0.458 0.196 0.411 0.570 0.598
[1] 7550.711

Now, starting at a=.6, c=.3, e=.4

              A1    A2    A3     C1     C2    C3    E1    E2    E3
var1     0.560    NA    NA  0.006     NA    NA 0.434    NA    NA
var2     0.652 0.343    NA  0.054  0.118    NA 0.294 0.539    NA
var3     0.636 0.477 0.277 -0.047 -0.047 0.124 0.411 0.570 0.598
[1] 7550.711

relevant fragment

share = mxModel("share", 
	mxMatrix("Lower", nrow=nVar, ncol=nVar, free=TRUE, values=.6, name="a", byrow=TRUE), # Additive genetic path coefficient
	mxMatrix("Lower", nrow=nVar, ncol=nVar, free=TRUE, values=.6, name="c", byrow=TRUE), # Common environmental path coefficient
	mxMatrix("Lower", nrow=nVar, ncol=nVar, free=TRUE, values=.6, name="e", byrow=TRUE), # Unique environmental path coefficient
	# Multiply by each path coefficient by its inverse to get variance component
  mxAlgebra(a %*% t(a), name="A"), # additive genetic variance
  mxAlgebra(c %*% t(c), name="C"), # common environmental variance
  mxAlgebra(e %*% t(e), name="E") # unique environmental variance
)

mzModel = mxModel(share, name="MZ", 
	expMZMeans,
	# MZ ACE algebra
	mxAlgebra(rbind (cbind(A+C+E , A+C),
                   cbind(A+C   , A+C+E)), dimnames = list(selVars, selVars), name="expCov"),
	mxData(mzData, type="raw"),
	mxFIMLObjective("expCov", "expMean")
)

dzModel = mxModel(share, name="DZ", 
  expDZMeans,
  mxAlgebra(rbind (cbind(A+C+E   , .5%x%A+C),
                   cbind(.5%x%A+C,    A+C+E)  ), dimnames = list(selVars, selVars), name="expCov"),
	mxData(dzData, type="raw"),
	mxFIMLObjective("expCov", "expMean")
)

#Build and Run ACE model
twinACEModel = mxModel("twinACE",mzModel, dzModel, mxAlgebra(MZ.objective + DZ.objective, name="twin"), mxAlgebraObjective("twin") )
twinACEFit = mxRun(twinACEModel)
Replied on Thu, 09/03/2009 - 15:09
Picture of user. tbates Joined: 07/31/2009

Interestingly, rewriting the script so that it uses three-models so that the MZ and DZ sub-models both do their algebra off a single set of ace and ACE matrices and algebras contained in the third shared model eliminates this instability, with the same solution reliably emerging from diverse starts.
Log in or register to post comments
Replied on Fri, 09/04/2009 - 08:32
Picture of user. Steve Joined: 07/30/2009

Hmmm.... That is interesting. We're trying to do the consolidation of the algebra in a smart way. It must be eliminating either some multiplication or a place where your first model had two parameters plus a constraint. I'd need to work through the model step by step and see where the two sets of convariance algebra differ.
Log in or register to post comments
Replied on Fri, 09/11/2009 - 23:08
Picture of user. neale Joined: 07/31/2009

The output fragment you provide indicates strongly that the model is not identified. Different parameter estimates emerge, but the -2lnL is the same. I expect that somehow the information from one group is not being used in the model-fitting.
Log in or register to post comments
Replied on Sat, 09/12/2009 - 19:50
Picture of user. tbates Joined: 07/31/2009

my error was to assume that there was a global space of matrix names across the supermodel and all contained groups so matrix names acted like labels equating all matrices with the same name to be the same.

The fix here was to use the name of the mz group in the algebra of the dz group, linking the two algebras to the same underlying matrices.

mzModel = mxModel(name="mz", 
	mxMatrix("Lower", nVar, nVar, TRUE, .6, name="a", byrow=TRUE), # Additive genetic path coefficient
	mxMatrix("Lower", nVar, nVar, TRUE, .6, name="c", byrow=TRUE), # Common environmental path coefficient
	mxMatrix("Lower", nVar, nVar, TRUE, .6, name="e", byrow=TRUE), # Unique environmental path coefficient
	# Multiply by each path coefficient by its inverse to get variance component
	mxAlgebra(a %*% t(a), name="A"), # additive genetic variance
	mxAlgebra(c %*% t(c), name="C"), # common environmental variance
	mxAlgebra(e %*% t(e), name="E"), # unique environmental variance

	expMZMeans,
	# MZ ACE algebra
	mxAlgebra(rbind (cbind(A+C+E , A+C),
	cbind(A+C   , A+C+E)), dimnames = list(selVars, selVars), name="expCov"),
	mxData(mzData, type="raw"),
	mxFIMLObjective("expCov", "expMeanMZ")
)


dzModel = mxModel(name="dz", 
  expDZMeans,
  mxAlgebra(rbind (cbind(mz.A+mz.C+mz.E   , .5%x%mz.A+mz.C),
                   cbind(.5%x%mz.A+mz.C,    mz.A+mz.C+mz.E)  ), dimnames = list(selVars, selVars), name="expCov"),
	mxData(dzData, type="raw"),
	mxFIMLObjective("expCov", "expMeanDZ")
)

#Build and Run ACE model
model = mxModel("twinACE",mzModel, dzModel, mxAlgebra(mz.objective + dz.objective, name="twin"), mxAlgebraObjective("twin") )
fit = mxRun(model)

It will important to clearly explain the scope of labels and names to users, and that information from different groups has to be imported by using the object's full name, including the originating group.

Easy enough once you know :-)

Log in or register to post comments
  Login

News