Login

News

Trying to run a bivariate saturated model

Posted on
No user picture. JuliaJoplin Joined: 10/12/2014
Forums
Behavioral Genetics Models
Hello! This is a follow-up post to some guidance I got on this board awhile ago that was incredibly helpful. I am trying to run a bivariate saturated model that controls for age and sex. The two variables were assessed at different time points, so one variable needs to be adjusted for age1 and sex, the other variable for age2 and sex.

I haven't found a script that does this, so I attempted to adapt an existing one. The entire purpose of running the saturated model is to 1) get cross-twin cross-trait correlations and 2) check assumptions to be sure everything can be constrained. That's all I want to do! My script for running the bivariate Cholesky seems to work fine.

I assume that this analysis would be fairly common given that assumption checking is part of the process and CTCT correlations are often reported. But, I'm just having a ton of problems with actually running it. Lots of Mx RED error codes, and then when I get to the final steps (constraining means, variances, and covariances, and then constraining CTCT), I am getting error code 10 about implausible start values. Sometimes when I use mxTryHard this goes away, other times it does not and there are errors, and still other times it takes so long that the system crashes.

Granted I'm still a novice, but I have never ever had so many errors (or error code 10), so I'm just wondering if I'm missing something obvious that I've done wrong or if there's a way to make my code less computationally intensive, as I have to repeat this analysis many times.

Any help would be much appreciated.

Replied on Wed, 07/27/2016 - 11:04
Picture of user. AdminRobK Joined: 01/24/2014

From inspecting your script, I don't see anything obviously wrong. The only concern I have is that you're parameterizing the saturated covariance matrices in terms of symmetric matrices. There's nothing wrong with that, but it can make your MxModel harder to optimize, because there's nothing keeping the expected covariance matrices inside the parameter space (i.e., positive-definite). An alternative would be to parameterize each covariance matrix in terms of a lower-triangular matrix times its transpose (i.e., a Cholesky parameterization), which will keep the resulting matrix positive-definite. I suspect this is the reason why mxTryHard() isn't really helping, and why you're getting status code 10: mxTryHard() is able to randomly generate start values such that one of the two model-expected covariance matrices is non-positive-definite; those would indeed be infeasible starting values.

I admit, though, that using the Cholesky parameterization will make it more complicated to impose the equality constraints required for the submodels you fit after the saturated model.

It might help if you posted the output from your script. If you do so, use argument verbose=TRUE with summary().

Edit: You create object lbVas on line 50 of your script. The comment on that line says "lower bound for variances", but you never use lbVas, and don't otherwise enforce positivity of the diagonal elements of the covariance matrices. Try defining lbVas differently:

lbVas <- matrix(NA,ntv,ntv)
diag(lbVas) <- 0.0001

Then, use it when defining the model-expected covariance matrices:

expCovMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=laVaMZ, name="expCovMZ", lbound=lbVas )
expCovDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, labels=laVaDZ, name="expCovDZ", lbound=lbVas )

Log in or register to post comments
Replied on Wed, 07/27/2016 - 12:32
Picture of user. AdminRobK Joined: 01/24/2014

As I look at your script again, I notice that there's something wrong with the model for the means:

BAGE <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .6, label=c("bage_caster", "bage_ogo"), name="betaage")
BSEX <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .6, label=c("bsex_caster", "bsex_ogo"), name="betasex")
BAGE2 <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .6, label=c("bage2_caster", "bage2_ogo"), name="betaage2")

OBSAGE <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="Age" )
OBSAGE2 <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="Age2" )
OBSEX <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="Sex" )

meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=.6, labels=laMeMZ, name="meanMZ" )
meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=.6, labels=laMeDZ, name="meanDZ" )

expMeanMZ <- mxAlgebra( expression= meanMZ + betasex %x%Sex + betaage %x%Age + betaage2 %x%Age2, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + betasex %x%Sex + betaage %x%Age + betaage2 %x%Age2, name="expMeanDZ" )

As it stands, both phenotypes are being regressed onto both age1 and age2. Furthermore, the Kronecker products don't make sense. The Kronecker product betaage %x% Age is the equivalent of

cbind(bage_caster*data.ageT1, bage_caster*data.ageT2, bage_ogo*data.ageT1, bage_ogo*data.ageT2)

which I'm pretty sure isn't what you want--it doesn't match the ordering of your variables.

Which phenotype is supposed to be regressed onto age1, and which is supposed to be regressed onto age2?

Log in or register to post comments
Replied on Wed, 07/27/2016 - 12:49
No user picture. JuliaJoplin Joined: 10/12/2014

Yup, in the original script, both phenotypes are being regressed onto both age1 and age2, which isn't what I ultimately want.

Caster is supposed to be regressed onto age 1, and Ogo is supposed to be regressed onto age 2.

I couldn't figure out how to build that into the original script, so then I had that line that says:

twinSat2Model <- omxSetParameters( twinSatModel, labels=c("bage_ogo", "bage2_caster"), free=FALSE, values=0 )
twinSat2Fit <- mxRun(twinSat2Model)

But it sounds like something still may be wrong?

Log in or register to post comments
Replied on Wed, 07/27/2016 - 13:41
Picture of user. AdminRobK Joined: 01/24/2014


Caster is supposed to be regressed onto age 1, and Ogo is supposed to be regressed onto age 2.

I couldn't figure out how to build that into the original script, so then I had that line that says:

twinSat2Model <- omxSetParameters( twinSatModel, labels=c("bage_ogo", "bage2_caster"), free=FALSE, values=0 )
twinSat2Fit <- mxRun(twinSat2Model)

Ah, I see now. If I were you, I wouldn't bother with fitting a model that doesn't have those two extraneous regression coefficients fixed to zero. You can fix them to zero ab initio with:

BAGE <- mxMatrix( type="Full", nrow=1, ncol=2, free=c(TRUE,FALSE), values= c(.6,0), label=c("bage_caster", "bage_ogo"), name="betaage")
BSEX <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .6, label=c("bsex_caster", "bsex_ogo"), name="betasex")
BAGE2 <- mxMatrix( type="Full", nrow=1, ncol=2, free=c(FALSE,TRUE), values= c(0,.6), label=c("bage2_caster", "bage2_ogo"), name="betaage2")

But it sounds like something still may be wrong?

Yes. Reverse the order of the factors in the Kronecker products:

expMeanMZ <- mxAlgebra( expression= meanMZ + Sex %x%betasex + Age %x% betaage + Age2 %x% betaage2, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + Sex %x%betasex + Age %x% betaage + Age2 %x% betaage2, name="expMeanDZ" )

Log in or register to post comments
Replied on Wed, 07/27/2016 - 22:10
No user picture. JuliaJoplin Joined: 10/12/2014

Thank you! I made these modifications, including the missed lbVas issue you mentioned above. I think I ran into trouble with the ordering of the factors because I adapted this from a univariate Saturated model for just one of the variables that was set up like this:

selVars <- c("happy1", "happy2")

expMeanMZ <- mxAlgebra( expression= meanMZ + betasex %x%Sex + betaage %x%Age, name="expMeanMZ" )

Anyway, I have tried to rerun the script, and it works okay till I get to the very last two models -- eqMVarsCTCTZygModel and eqMVarsZygModel. eqMVarsZygModel I can get to work with mxTryHard, but the other one does not seem to be working. As you mentioned, here's the output!

I haven't paid very close attention to the start values, so that could be a place to start assuming this reparameterization is going to be pretty tricky. (I had another script that I believe did do it in this way but it didn't adjust for covariates.) I know I do need to do this to be sure the assumptions are met.

Log in or register to post comments
Replied on Wed, 07/27/2016 - 11:32
Picture of user. AdminNeale Joined: 03/01/2013

Are the variances and means from which you're starting optimization good for the variables you're analyzing? It's usually best to start the variances a bit larger than the observed, and to try to get the means about right.

I agree with Rob that reparameterizing as a Cholesky would be more robust, but perhaps you would get away with the symmetric specification as long as the starting values are decent.

Log in or register to post comments
  Login

News