Hi,

I have a dataset of twins on 4 time-points. I can fit the Multivariate cholesky decompositon model with ACE components. It considers unstrcutured covariance matrices for A, C and E component. I was wondering if its possible to restrict it to a special structure? Like, say compound symmetric or AR1.

Thanks

Yes, that's possible, but it would be difficult to do under the Cholesky parameterization. It could be done with Cholesky, but you'd need to use MxConstraints on the off-diagonal elements of the

A,C, andEcorrelation matrices.I recommend instead using a correlated-factors paramterization. Under that parameterization, the

A,C, andEcovariance matrices would each be equal to a matrix of loadings times a factor correlation matrix times the transpose of the matrix of loadings. There would be a loadings matrix and a factor correlation matrix for each ofA,C, andE. Each phenotype would load only onto its ownA(or whatever) factor; the loadings matrix could just be a diagonal matrix. The factor correlation matrix would have the structure you want. For compound symmetry, the factor correlation matrix could be an MxMatrix of`type="Stand"`

and the same parameter label for all the off-diagonals. For AR1, I would suggest making a 1x1 matrix for the first-order correlation, and then define the factor correlation matrix as an MxAlgebra assembled, via`cbind()`

and`rbind()`

, from the first-order correlation raised to the appropriate power per cell (with ones on the diagonal of course). Loehlin (1996) and Neale, Røysamb, & Jacobson (2006) may be useful references.Can you guide me how to addd constraints in Cholesky Model? This is how I am currently specifying.

OK, first make

A,C, andEcorrelation matrices:Now, you need to impose the desired structure via MxConstraints. To impose a compound-symmetric correlation structure on

A:Do likewise for

CandE. You can give the MxConstraints names if you want, via argument`name`

. You might be able to come up with a way to do all this using fewer`mxConstraint()`

statements, perhaps with the`vechs()`

function. At any rate, you need to put these new objects into either the MZ or DZ group model.For AR1, try:

And likewise for

CandE.If you take this route, I don't recommend the on-load default optmizer, CSOLNP. Switch to SLSQP (or perhaps NPSOL if your build of OpenMx supports it), e.g.

`mxOption(NULL,"Default optimizer","SLSQP")`

.That is great!!

Can we also add multiple equalities like,

Acon1 <- mxConstraint(corA[1,2] == corA[2,3]==corC[1,2])

That would be nice, but I'm pretty sure the answer is "no" (though I admit that I've never tried). What happens at runtime is that OpenMx converts the

`expression`

you provided to`mxConstraint()`

into a standard form, which has a matrix-valued function on the LHS of the comparator, and on the RHS, a matrix of the same dimensions containing all zeroes. So, you could use fewer`mxConstraint()`

statements (and thereby make fewer MxConstraint objects) if you did something like this:I used the above code and was looking at results, but I still get an unstructured correlation matrix.

Well, something's clearly not right. Please post your full script.

Edit: Also, what's your full output from

`mxVersion()`

?This is the code.

> mxVersion()

OpenMx version: 2.9.9 [GIT v2.9.9]

R version: R version 3.3.1 (2016-06-21)

Platform: x86_64-w64-mingw32

Default optimiser: CSOLNP

NPSOL-enabled?: Yes

OpenMP-enabled?: No

Like I said, the new objects have to actually go into an MxModel object. For instance,

. I have a few additional remarks about the script as well. First, is there a reason you're not correcting for the main effects of age and sex? Second, I don't see any syntax in your script where the dataset is loaded. Are you just restoring an R workspace at the start of the session? If so, I strongly discourage doing that.

I am getting this error,

do you think it is due to the data? I mean dataset does not suit the used model and hence estimates are getting biased and wrong.

Running CholACE with 34 parameters

Error: The job for model 'CholACE' exited abnormally with the error message: fit is not finite (The continuous part of the model implied covariance (loc2) is not positive definite in data 'MZ.data' row 28. Detail:

covariance = matrix(c( # 8x8

0.480003, 0.320338, 0.324131, 1.9695e+019, 0.320001, 0.160337, 0.164131, 1.9695e+019

, 0.320338, 3.22184e+026, 3.25516e+031, 4.14316e+016, 0.160337, 3.22184e+026, 3.25516e+031, 4.14316e+016

, 0.324131, 3.25516e+031, 1.84611e+038, 1.41805e+018, 0.164131, 3.25516e+031, 1.84611e+038, 1.41805e+018

, 1.9695e+019, 4.14316e+016, 1.41805e+018, 2.42432e+039, 1.9695e+019, 4.14316e+016, 1.41805e+018, 2.42432e+039

, 0.320001, 0.160337, 0.164131, 1.9695e+019, 0.480003, 0.320338, 0.324131, 1.9695e+019

, 0.160337, 3.22184e+026, 3.25516e+031, 4.14316e+016, 0.320338, 3.22184e+026, 3.25516e+031, 4.14316e+016

, 0.164131, 3.25516e+031, 1.84611e+038, 1.41805e+018, 0.324131, 3.25516e+031, 1.84611e+038, 1.41805e+018

, 1.9695e+019, 4.14316e+016, 1.41805e+018, 2.424

In addition: Warning message:

In model 'CholACE' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard()

Possibly, but some of those elements of the model-expected covariance matrix just look ridiculously wrong. I think it's more likely that the system of nonlinear constraints your model imposes is making optimization difficult, and/or the start values in your script are poor. A few tips:

I'll remind you that you're using the "Cholesky parameterization plus MxConstraints" approach to fit these structured-correlation-matrix models, which goes against my recommendation (to instead use a correlated-factor parameterization) in the first place. You are not the first user to encounter difficulties using such an approach.

I am getting fine result on a simpler data. I guess for the other one, the starting values need to be more accurate.

Regarding this,

Acon1 <- mxConstraint(corA[1,2] == corA[2,3]==corC[1,2])

Do you think it is possible in OpenMx to allow such restriction in any way, i.e, correlation of A component is the same as that of C? like if I do

Acon1 <- mxConstraint(corA[1,2] == corA[2,3]

Ccon1 <- mxConstraint(corC[1,2] == corA[2,3]

Would it work?

I tried once, it didnt work.

That is indeed possible to do. Assuming you put the new MxConstraints into an MxModel, my guess for why "it didn't work" is that the optimizer is failing to find a solution that satisfies the constraints.

I used this code,

It does gives me A, C and E unstructured. But if I keep only Acon1,..,Acon5 in pars, it gives me compound symmetric correlation structure of A.

Hmm. That's odd. Try again, retaining

`Acon1`

thru`Acon5`

, and including two new MxConstraints:I tried what you said. Only retained con1 through Acon5 and added last two constraints. It is weird it did not work.

Could you post some post-

`mxRun()`

output, such as the estimated values of the 3 correlation matrices, and the`summary()`

output? I might be able to diagnose what's going wrong.My advice to adopt the correlated-factors parameterization still stands.

This is when I specify compound symmetry for all 3,

A

[,1] [,2] [,3] [,4]

[1,] 44.04476 18.84978 24.12863 20.55763

[2,] 18.84978 44.75724 25.76197 20.95504

[3,] 24.12863 25.76197 49.13861 23.88618

[4,] 20.55763 20.95504 23.88618 42.12631

> C

[,1] [,2] [,3] [,4]

[1,] 5.9019408291 4.1313542437 0.0002429391 6.765129377

[2,] 4.1313542437 3.0425636732 0.0002088668 4.935798747

[3,] 0.0002429391 0.0002088668 1.9565944167 0.000469937

[4,] 6.7651293767 4.9357987469 0.0004699370 8.020701197

> E

[,1] [,2] [,3] [,4]

[1,] 13.707064 7.331347 6.261522 5.814909

[2,] 7.331347 12.460888 6.820699 6.240138

[3,] 6.261522 6.820699 10.797324 5.555669

[4,] 5.814909 6.240138 5.555669 11.915593

This is when I specify CS for only A,

[1,] 42.59026 22.69322 23.48385 21.82458

[2,] 22.69322 45.94477 24.39115 22.66777

[3,] 23.48385 24.39115 49.20198 23.45751

[4,] 21.82458 22.66777 23.45751 42.49477

>

> C

[,1] [,2] [,3] [,4]

[1,] 7.8527242 1.2478234 0.6538513 5.8108463

[2,] 1.2478234 1.7890718 0.1040252 3.2639430

[3,] 0.6538513 0.1040252 1.3143163 0.4841337

[4,] 5.8108463 3.2639430 0.4841337 7.7436756

> E

[,1] [,2] [,3] [,4]

[1,] 13.683333 7.053824 6.238655 5.705981

[2,] 7.053824 12.407678 6.930252 6.152708

[3,] 6.238655 6.930252 10.851366 5.556251

[4,] 5.705981 6.152708 5.556251 11.896003

OK, but what do you get from

`mxEval(MZ.corA, CholAceFit, TRUE)`

, and likewise forCandE? Remember, the constraints are imposed on the correlation matrices, not the covariance matrices. And what about the`summary()`

output?