Hi everyone,

I am trying to do a joint analysis between 2 binary variables and 1 continuous variable I am having some problems to fix the variance for the 2 binary variables equal 1.

I think I should change something here

matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )

matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc" )

var1 <- mxConstraint( expression=diag2vec(Oc%&%V)==Unv1, name="Var1" )

NVO=2

NV=3

However, I have done several changes and it doesn’t work since in some tries the matrix are not conformable and in others it doesn’t fix the variance to 1. How could I fix the variance for the two binary variables to 1?

Thank you so much in advance.

Juan J Madrid-Valero

Without the context of the full script, I can't really answer your question.

Thanks Rob for your prompt response here you can find the full script:

At least in the script you posted,

`Oc%&%V`

is going to evaluate to a 1x1 matrix, but`matUnv`

is 2x1, and the two sides of the comparator in an MxConstraint have to be of the same dimensions. A simple way to fix the variance of the two binary variables to 1 is to replace those three lines with:You could have also just used two MxConstraint objects (one for each variance).

Since your analysis involves an MxConstraint, I recommend using SLSQP as the optimizer.

You could also avoid using an MxConstraint entirely if you change the parameterization of the model to direct variance components, and identify the scale of the ordinal variables by fixing their nonshared-environmental variances to 1.

Thank you so much Rob,

I didn’t know why the matrix was ( expression=diag2vec(Oc%&%V)) 1x1.

Also, I have run the script and the results are pretty feasible but I got this message “In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, : The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)”

I have tried tons of starting values but I have not been able to remove the “problem”

By the way I have used the optimizer that you told me and the function MxtryhardOrdinal with extra tries.

Juan J Madrid-Valero

Here the full message:

Retry limit reached

Solution found

Start values from best fit:

0.0873764505129859,-0.00302648728148331,0.00387006292591793,0.446376724057658,0.765174696854374,1.24795735809318,2.27769791884809,-0.336483766998107,-0.275906439822281,0.416526067571617,0.471557647293468,0.00245202168140257,0.202289271880198,0.352135547711807,0.0718736069163078,0.00290227641677533,0.00372693059844336,0.00118080806141676,3.1768820509539,-0.106534386756017,-0.147029899250398,0.76021560181588,0.40208133001671,0.716267708855051

Warning message:

In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, :

The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

`Oc`

is 1x3. Operator`%&%`

is the quadratic-product operator (see the documentation for`mxAlgebra()`

). Thus,`Oc%&%V`

is equivalent to`Oc%*%V%*%t(Oc)`

, which is (1x3) times (3x3) times (3x1), which evaluates to a 1x1 product.Which version of OpenMx are you using? As of v2.8.3, status Red should occur less frequently in ordinal analyses than it used to. Nonetheless, status Red is sometimes unavoidable in ordinal analyses, in that it can still occur even when the optimizer has found a local minimum. You could try using CSOLNP instead of SLSQP. If that still doesn't help, then it might be worth reparameterizing to eliminate the MxConstraint, which I alluded to in another post.

Thanks Rob,

I am using R version 3.4.3 and I have tried with the CSOLNP optimizer but it doesn't work either. I am going to try to change the parameterization of the model to direct variance components. I am not sure how I should do it. Could you please give me some guidance or where I could find information about that?

Juan J Madrid-Valero

What about your OpenMx version?

Let's see...one way, involving minimal changes to your script, would be as follows. First, remove the line that creates

`pathE`

. Then, define`covE`

as:Then, remove the lines that define

`matUnv`

,`matOc`

, and`var1`

. Finally, remove`pathE`

,`matUnv`

,`matOc`

, and`var1`

from the definition of`pars`

. What these changes do is identify the latent scale of the ordinal variables by fixing theirnonshared-environmental variance, rather than theirtotal variance, to 1.0. As the MxConstraint is no longer required, use the CSOLNP optimizer.While I'm thinking about it--consider removing the lower bounds for

`pathA`

and`pathE`

. You can put them back in if the optimizer keeps reaching really crazy solutions.Sorry! I am using OpenMx version: 2.8.3 [GIT v2.8.3] ; R version: R version 3.4.3 (2017-11-30).

I have made all the changes and now It works and I don't have a red mx status but I have this message:

Retry limit reached

Solution found

Running final fit, for Hessian and/or standard errors and/or confidence intervals

Warning messages generated from final fit for final fit for Hessian/SEs/CIs

Start values from best fit:

0.086752181736907,-0.00393385954682701,0.00457737378618429,0.481260848285712,0.526771126020844,1.19685352917708,2.29257177590332,-0.392479607148631,-0.315671611477633,0.782950339036923,0.501766203589714,-0.266923949247171,10.0657799256565,-0.474366335537433,-0.565466421429395,0.490967007947968

the resuts are the same that with the mx status red.

Thank you so much!!!

Do you think that these results are reliable?

by the way is it possible to know the influences shared by the three phenotypes (all together), not only two by two phenotypes?

and also I would like to know what is the easiest way to compare two correlations in the saturated model

I have tried creating these objects:

CICORMZF <-mxAlgebra(expression =cov2cor(MZf.covMZf), name="CORMZF")

CICORDZF <-mxAlgebra(expression =cov2cor(DZf.covDZf), name="CORDZF")

CICORMZM <-mxAlgebra(expression =cov2cor(MZm.covMZm), name="CORMZM")

CICORDZM <-mxAlgebra(expression =cov2cor(DZm.covDZm), name="CORDZM")

CICORDZO <-mxAlgebra(expression =cov2cor(DZo.covDZo), name="CORDZO")

rAcompare <- mxAlgebra(cbind(CORDZF[1,2] - CORDZM[1,2], name="rAcompare")

and adding them to the model but it doesn't work.

It is a 5group saturated model and I would like to compare the correlations between males and females.

I have thought to equal the correlations in the model using omxsetparameter but in the model I don't have the correlations for the five groups only in the submodels so I am not sure how to proceed.

Thanks!!!

How do you mean?

Yes, that can be done. Look into the Common-Pathway and Independent-Pathway models. Both of them do what you're describing, but in different ways.

One thing wrong I notice immediately is that there's an unbalanced parenthesis in

. I don't think you need the column-bind, so try changing it to

. Otherwise, that block on syntax looks OK to me, but it's hard to say without seeing it in the context of the full script, and without any specifics beyond "it doesn't work." Since you're referencing the covariance matrices by "modelname-dot-objectname", I assume you're placing these algebras in the "container" MxModel, right?

Sorry! When I said if these results are reliable I wanted to say if with this message:

"Retry limit reached

Solution found

Running final fit, for Hessian and/or standard errors and/or confidence intervals

Warning messages generated from final fit for final fit for Hessian/SEs/CIs

Start values from best fit:

0.086752181736907,-0.00393385954682701,0.00457737378618429,0.481260848285712,0...."

I could be sure that the soltuion is "good" I mean since retry limit have been reached I don't know if I should try other starting values or something else (although there is a solution)

Here the 5 group script I am working with:

I would like to add to the model the correlations for each group (MZM,MZF,DZM,DZF,DZO) I have this correlations in the submodels but not in the model. I mean when I run summary(fit) this values are only in the Confident invervals.

Therefore I can't use omxsetparemateres to equal the correlations and test the model (for example corMZF equal to corMZM)

I think I should create five matrix for the correlations but I don't know how to do this. Also I tried adding CICORMZF, CICORMZM... to the model and submodels but I only have the correlations in the submodels.

Aren't you working with two binary phenotypes? If so, the saturated model still needs to have thresholds, and constraints imposed to identify the scale of the latent liability.

Yes sorry, I am working with two binary variables and one continous. This script is for the continous variable. I am testing if there are sex differences (in the continous variable). It is a 5 group saturated univariate model but I don't know how to add the correlations for the 5 groups to the model (actually I have them in the submodels)

What was the warning generated during the final fit? Did you get status Red or not? If you're really concerned, you could increase the number of extra tries. You could also try my Nelder-Mead optimizer, but I doubt it would work any better than CSOLNP.

Ah, I see. The reason for why that doesn't work is unrelated to whether the correlations are in the submodel or in the container (so, you don't need to create new objects for the correlation matrices; the ones you have in the submodels already will do). The reason for why that doesn't work is that the correlations are not statistical parameters of the MxModel object. Parameters are elements of MxMatrices (or values of MxPaths), but the correlations are elements of MxAlgebras. That is, they are not parameters themselves, but are instead

functionsof parameters (informally called "algebras").`omxSetParameters()`

and the 'free parameters' table in`summary()`

output only deal with parameters, not algebras. In contrast, confidence intervals can be calculated for algebras. MxConstraints can also be imposed on algebras.If you want to fit models in which the correlations are constrained across groups in some way, you basically have two options. One option is to reparameterize your MxModels so that the correlations are in fact parameters rather than algebras. For advice on doing that, see this other thread. The other option, which would involve less work, is to use MxConstraints. For instance, if you want to constrain the MZF and MZM correlations to be equal to one another, make an object like this,

and put it in the container MxModel. I recommend using SLSQP when MxConstraints are involved.

You could instead request confidence intervals on the differences in the correlations. For instance,

, and put them into the container MxModel.

Hi Rob! thank you so much again!

1-I think the warnings are for the CI since I have NA in some of them. I don't have mx status red now.

The message is "Retry limit reached

Solution found

Running final fit, for Hessian and/or standard errors and/or confidence intervals

Warning messages generated from final fit for final fit for Hessian/SEs/CIs

Start values from best fit:

0.086752181736907,-0.00393385954682701,0.00457737378618429,0.48126084828...."

I have tried changing the optimizer and 100 extra tries but I have the same resuls and always retry limit is reached.

In this model we have fixed the nonshared-environmental variances to 1 but why not for A and C?

2-Where can I find information about Common-Pathway and Independent-Pathway models? I need to use another script right?

3- With regards to the correations for the five group univariate model it works perfect now! thank you!!!

If you run

`summary()`

with argument`verbose=TRUE`

, you can see diagnostics for the confidence intervals. Also,`mxTryHardOrdinal()`

runs with default argument`exhaustive=TRUE`

, so by default it will always use up all of its extra tries.You could always look around in the materials from the 2016 Boulder Workshop, e.g. here. See also:

Posthuma, D. (2009). Multivariate genetic analysis. In Kim, Y-K. (Ed.),

Handbook of Behavior Genetics(pp. 47-59). New York: Springer Science+Business Media, LLC.Neale, M.C., & Cardon, L.R. (1992).

Methodology for genetic studies of twins and families. Dordrecht, The Netherlands: Kluwer Academic PublishersHi Rob,

Thanks for your response. The diagnosis is "alpha level not reached infeasible start". I have tried with the optimizer "SLSQP" and removing lower bounds on the ACE path coefficients (as I saw in another of your post). But It does not work. It took more than 20 hours and I had to stop the analysis.

What should I do?

Thanks!!

I'd like to see the verbose output from

`summary()`

.Hi Rob,

Here you can find the output (with CSOLNP optimizer):

Hi, Juan. Sorry for the delayed reply. I was at the Boulder Workshop last week.

I am not really sure that CSOLNP found a minimum of the fitfunction. The information matrix isn't positive-definite, and the gradient is asymmetric around every free-parameter estimate. How many extra tries did you run it with? Nonetheless, with ordinal data, I guess it's still possible that's close to a solution. Could you post your current script, so I know exactly what you're running?

Hi! Don't worry! I wanted to go but.... I couldn't :(

Here the script. I think I also tried with 101 extraTries and It was the same result.

Thank you so much

Do you already know for sure that you're going to report CIs from the ACE model? If not, you should first try to be sure you're getting to a decent maximum-likelihood solution for all four models, select the one with smallest AIC, and only then worry about confidence intervals. Try

`mxTryHardOrdinal()`

without`intervals=T`

; stick with CSOLNP for now. You can subsequently use`omxRunCI()`

on a fitted MxModel object to calculate CIs without redoing the primary optimization to find the MLE.I don't see that the primary optimization should be too difficult. Granted, you're working with ordinal data, but your ordinal variables are binary with neither outcome being rare, and you don't have a whole lot of missing data.

The CI details table is so big that it can be hard to read. You could simplify it if you request CIs for a vector of MxAlgebra elements, e.g.,

and just skip the diagonal elements of the correlation matrices.

It would be a good idea to make your

`mxTryHardOrdinal()`

results reproducible. The best way to do so is probably to add a line likebefore each call to

`mxTryHardOrdinal()`

. Also, if you pass argument`bestInitsOutput=TRUE`

to`mxTryHardOrdinal()`

, you can see the start values that the program used to find its best solution, and then modify your script to actually use those start values.Hi Rob,

Thank you so much, now the model estimates all the CI for the best fitted model. I have adjust the starting values and I have modified the line mxCI( c("VC[1,1]…. As you told me.

I had to use the optimizer “SLSQP” is it ok?

Just one more question I am not sure why we have to fix to 1 only nonshared-environmental variance and not for the other components?

Thank you so much again!!!

If it worked, it seems OK to me!

The way OpenMx models ordinal variables is by assuming that the ordinal variable represents a latent, normally-distributed continuous variable that has been discretized into ordered categories by "cutting" at the thresholds (a binary variable will only have 2 categories and 1 threshold). Since the underlying continuous variable isn't observed, the model isn't identified without some additional constraint. With a binary variable, either the threshold or the underlying mean has to be fixed to some value; typically the variance is also fixed, to 1.0. However, in an ACE twin model, fixing the variance (i.e., the total phenotypic variance) to 1.0 requires an MxConstraint. I was concerned the MxConstraint was complicating the optimization of your model, so I suggested an alternative constraint, which is to fix the nonshared-environmental variance to 1.0. When the model is identified in that fashion, the values of the additive-genetic and shared-environmental variance components are interpretable as the factors by which those variance components are proportionately larger than (if >1.0) or smaller than (if <1.0) the nonshared-environmental variance.

Thank you so much!!!

Hi Rob,

Just reading this great thread which has helped me a lot in setting up my model and I am wondering about omxRunCI. This sounds like a great function to calculate the CI's without refitting the whole model. However, I cannot find it. I have OpenMx version 2.6.9 and when I type ?omxRunCI, it says that no documentation is found. Could you please advise?

Thank you,

Julia

Sorry, my bad! Found out that this function is only available with version 2.9.6. Just updated OpenMx, but can't figure out how to use omxRunCI(). The CI's are specified as follows.

And when I pass AceFit into omRunCI, I get the following:

Could you please help me with this one?

Actually, it's available in v2.7.18 or newer.

I am willing to bet a lot of money that the MxModel and other OpenMx objects were created with an older version of the package--something older than v2.7.x--and were saved to disk in an .RData file (possibly as part of R autosaving your workspace at the end of a session) and subsequently reloaded after you updated to v2.9.6. The MxConstraint object class was changed between version 2.6.x and version 2.7.4 by adding new slots to it, in order to allow users to provide analytic constraint Jacobians. You'll need to reconstruct the MxConstraint objects that are in your MxModel, and put them into the MxModel in some way. Probably the best thing to do is just re-run your whole script, since it's very possible that other object classes have had their signatures altered since v2.6.x. Alternately, you could reconstruct the MxConstraints only, and directly overwrite the 'constraints' slot of the model (or of the submodels, if that's where the MxConstraints go).

If you could run

`traceback()`

right after seeing the error message and post the output, that would be helpful, since it's possible some code changes could be made so that future OpenMx releases don't suffer from this lack of backward compatibility.Yes, the models were indeed built under the older 2.6.9 version of OpenMx. Actually now when I updated to 2.9.6 and tried to load workspace with my fitted models, I can't even run the summary of the models. Getting the following:

I can extract estimated parameters though. Is there a way to get the summary of the model fit or will I have to rerun all the models (it takes a couple of hours for each)...

I think you're going to have to re-run them. Or downgrade back to v2.6.9.

You could copy the optimal parameter vector from the old models and use that as starting values. When you re-run the models, it should finish much faster.

Thanks a lot! Will have to do that!

Good idea, Joshua. I have a feeling, though, that the confidence-interval optimization will take the lion's share of the computational effort.

Hello, I want to run a G*E model (2 binary variables ) using the script from “JuanJMV”. Unfortunately, there seems to be some errors in my model

I just met the following Warning message when I add one of the Constraint lines

Warning message:

In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, :

The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

## Constraint on variance of Binary variables

matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )

var_constraint <- mxConstraint( expression=diag2vec(V0)==Unv1, name="Var1" )

or

matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )

var_constraint <- mxConstraint(expression=rbind(V0[2,2],V0[1,1]) == Unv1, name="Var1")

Another errors occurred when I remove the Constraint lines and define covE as:

myVarE <- mxMatrix( type="Symm",nrow=nv, ncol=nv, free=c(F,T,F), values=c(1, 0.1, 1),

labels=c("e_1_1", "e_2_1", "e_2_2"), lbound=c(0.0001,rep(NA,2)),name="E0")

** Information matrix is not positive definite (not at a candidate optimum).

Be suspicious of these results. At minimum, do not trust the standard errors.

Any suggestions from you?

Thanks!