Login

News

Correlated Genetic Factors in Independent Pathways Model

Posted on
No user picture. jarrode28 Joined: 03/01/2010
Forums
Behavioral Genetics Models
Hi,

I'm attempting to run a BG script of an independent pathways model on 9 dichotomous (symptom) manifest variables. That is, the model includes biometric factors specific to each manifest variable and common to all manifest variables. I've been able to get this model to run without any problems.

Based on some theoretical and empirical work, I'm also attempting to run a similar model with two correlated A factors; one accounts for the genetic variance in three manifest variables, with the other accounts for the genetic variance in the remaining six variables. This is where I'm running into some problems.

I've specified the two A factors as follows:

Ac1Free <- c(rep(T,6),rep(F,3))
Ac2Free <- c(rep(F,6),rep(T,3))
Ac1Values <- c(rep(1,6),rep(0,3))
Ac2Values <- c(rep(0,6),rep(1,3))
pathAc1 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=Ac1Free, values=Ac1Values, labels=labFull("ac1",nv,nf), name="ac1" )
pathAc2 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=Ac2Free, values=Ac2Values, labels=labFull("ac2",nv,nf), name="ac2" )

And I've attempted to specify that their correlated as follows:

pathRac <- mxMatrix (type="Diag", nrow=1, ncol=1, free=TRUE, values=0.5, label="rac1", lbound=0, ubound=1, name="rac")
covAc <- mxAlgebra( expression=ac1 %*% rac %*% t(ac2) + ac2 %*% rac %*% t(ac1), name="Ac" )

This model runs, but I'm fairly certain I'm doing something wrong because my degrees of freedom and estimated parameters are unchanged from the first model.

Any input on where I might be going wrong would be greatly appreciated.

Thanks,
Jarrod

Replied on Fri, 11/21/2014 - 16:36
Picture of user. RobK Joined: 04/19/2011

The recommended way to alter parameters in an MxModel object is to use omxSetParameters(). What that function does is return a new MxModel object that has your intended changes made to its parameters. In your case, I think what you'd want to do around line 200 is something like

Ind2A1C1EeModel <- omxSetParameters(IndAceModel, labels=c(pathAc1@labels[7:9],pathAc2@labels[7:9]),
free=c(F,F,F,T,T,T), values=c(0,0,0,0.6,0.6,0.6), name="Ind2Ace")

It looks like you've realized that OpenMx looks only at the contents of an MxModel at runtime, and not at the values of anything in the global environment, so therefore redefining the pathAc1 etc. objects in R's workspace won't free/fix the parameters. It looks like you also attempted to alter the contents of the MxModel, such as with Ind2A1C1EeModel$MZ$ac1 <- pathAc1 on line 223. Though not recommended, that might work in some cases. In any event, omxSetParameters() is the recommended way to go, because it changes all instances of a given parameter (identified by its label) in an MxModel.

Incidentally, I could tell that you had the same free parameters in both models by using omxGetParameters(IndAceModel) and omxGetParameters(Ind2A1C1EeModel).

Log in or register to post comments
Replied on Mon, 11/24/2014 - 00:00
No user picture. jarrode28 Joined: 03/01/2010

This was helpful and did the trick--thanks! I have a follow-up question, which may fall either fall under issues I'm having with model fitting or OpenMx, if you don't mind taking a look at my attached output.

Standard errors are not provided for any of the loadings on my common A factors, or the correlation between my common A factors. There aren't any issues (e.g., Code Red, Green, Blue, etc.) that appear when running the model, but I did have to try several starting values for the genetic correlation between these factors before the model would fit without any problems.

Log in or register to post comments
Replied on Mon, 11/24/2014 - 13:12
Picture of user. RobK Joined: 04/19/2011

For standard errors, you're seeing a mixture of real numbers and NaN. That's almost certainly because the final Hessian matrix, evaluated at the solution, is not positive-definite. You can verify that for yourself via eigen(Ind2A1C1EFit0200@output$calculatedHessian,T,T)$values; if any of the returned values are negative, then you know the Hessian (the matrix of second partial derivatives) isn't positive-definite. The most common reason for a non-positive-definite Hessian is that the optimizer has not actually found a minimum of the objective function. I suppose it could also mean that the solution it found is on a boundary for one or more parameters. That appears to be true in your case, since several of the specific-factor loadings are hitting their lower bounds of 1e-05.

I don't think there's much you can do besides try different start values and see if you can obtain a solution with the same or smaller -2 log likelihood, but with real values for all SEs. If you're willing, you could update to the beta of OpenMx version 2.0, which includes a function, mxTryHard(), that might be useful.

Edit: It looks like your observed variables are all binary 0/1 variables? You should probably treat them as ordinal variables instead.

Log in or register to post comments
Replied on Tue, 01/27/2015 - 13:45
No user picture. jarrode28 Joined: 03/01/2010

A quick follow-up regarding the models we're running here. Is this the appropriate way to specify that I want to estimate the genetic correlation between in an independent pathways model (with 9 manifest variables, 6 loading onto one factor, 3 loading onto the other)?

Ac1Free <- c(rep(T,6),rep(F,3))
Ac2Free <- c(rep(F,6),rep(T,3))
Ac1Values <- c(rep(0.3,6),rep(0,3))
Ac2Values <- c(rep(0,6),rep(0.3,3))

pathAc1 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=Ac1Free, values=Ac1Values, labels=labFull("ac1",nv,nf), name="ac1" )
pathAc2 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=Ac2Free, values=Ac2Values, labels=labFull("ac2",nv,nf), name="ac2" )

pathAs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=labDiag("as",nv), lbound=.00001, name="as" )

covAc <- mxAlgebra( expression=ac1 %*% rac %*% t(ac2) + ac2 %*% rac %*% t(ac1) , name="Ac" )

pathRac <- mxMatrix (type="Diag", nrow=1, ncol=1, free=TRUE, values=0.6, label="rac1", lbound=-.9999, ubound=.9999, name="rac")

Thanks!

Log in or register to post comments
  Login

News