You are here

Home › Forums › OpenSEM Forums › Behavioral Genetics Models

Equality constraints in saturated model

4 posts / 0 new
Log in or register to post comments
Last post
Tue, 07/17/2018 - 11:10
#1
IvanVoronin's picture
IvanVoronin
Offline
Joined: 08/18/2013 - 15:33
Equality constraints in saturated model
AttachmentSize
Binary Data openmx_select.R3.1 KB

I'm writing a function that would define saturated model based on the twin model that had already been fit. For the multivariate twin models I have to constrain the within-twin variances and co-variances to be equal across twins and across zygosity groups. I define the expected co-variance matrices using Cholesky decomposition, so I cannot introduce the constraints by labeling the matrices. At the same time, the reference page for the mxConstraint says that the 'Use of mxConstraint() should be avoided where it is possible to achieve the constraint by equating free parameters by label or position in an MxMatrix or MxAlgebra object'.

So the question is: what is the better way to introduce the constrains given that the constrained elements are computed, not defined as matrices?

Also: I tried to use mxConstraint and omxSelectRowsAndCols to introduce the constraints, and it seems that omxSelectRowsAndCols looks for the selector matrix in the global environment, not in the enclosing environment (the minimal working example attached).

Any response appreciated!
Ivan

Top
Tue, 07/17/2018 - 12:44
AdminRobK's picture
AdminRobK
Offline
Joined: 01/24/2014 - 12:15
I reproduce the behavior you

I reproduce the behavior you report. The part that worried me most was the freezing at runtime. But, that appears to be optimizer-specific. It only happens to me with CSOLNP, the on-load default optimizer. SLSQP and NPSOL don't freeze. In fact, it's possible CSOLNP isn't even freezing, but just taking a really long time.

As for the environment where omxSelectRowsAndCols() looks to find the selector matrix: the reason OpenMx can't find 'selecttw1' on the first mxRun() attempt is that 'selecttw1' does not appear in the MxModel's namespace, nor is there an R object with symbol selecttw1 in the environment at runtime. You could modify your function to create an MxMatrix for both of the selector matrices, and pass the names of the MxMatrices to omxSelectRowsAndCols(). For instance:

def_sat_model <- function(data, zyg, vars) {
    # zyg - zygosity variable (character string)
    # vars - phenotypic variables (character string or vector)
    nv <- length(vars)
    selvars <- paste0(vars, rep(1:2, each = nv))
 
    # Split MZ and DZ data
    mzdata <- data[data[, zyg] == 1, selvars]
    dzdata <- data[data[, zyg] == 2, selvars]
 
    # Starting values - These are for this particular data
    startChol <- matrix(c(0.8, 2.7, 0.4, 1.3,
                                                0.0, 0.3, 0.0, 0.0,
                                                0.0, 0.0, 0.7, 2.4,
                                                0.0, 0.0, 0.0, 0.3),
                                            ncol = 4, byrow = TRUE)
 
    # Selector matrix for the omxSelectRowsAndCols
    #selecttw1 <- matrix(c(rep(1, nv), rep(0, nv)), ncol = 1)
    #selecttw2 <- matrix(c(rep(0, nv), rep(1, nv)), ncol = 1)
 
    mxModel(name = 'Saturated',
                    mxMatrix(type = 'Lower', nrow = 2 * nv, ncol = 2 * nv,
                                     free = TRUE, values = t(startChol), name = 'CholMZ'),
                    mxMatrix(type = 'Lower', nrow = 2 * nv, ncol = 2 * nv,
                                     free = TRUE, values = t(startChol), name = 'CholDZ'),
                    mxAlgebra(CholMZ %*% t(CholMZ), name = 'expCovMZ'),
                    mxAlgebra(CholDZ %*% t(CholDZ), name = 'expCovDZ'),
                    mxMatrix(type="Full",nrow=4,ncol=1,free=F,values=c(rep(1, nv), rep(0, nv)),name="selecttw1"),
                    mxMatrix(type="Full",nrow=4,ncol=1,free=F,values=c(rep(0, nv), rep(1, nv)),name="selecttw2"),
 
                    mxConstraint(omxSelectRowsAndCols(expCovMZ, selecttw1) ==
                                                omxSelectRowsAndCols(expCovMZ, selecttw2)),
                    mxConstraint(omxSelectRowsAndCols(expCovDZ, selecttw1) ==
                                                omxSelectRowsAndCols(expCovDZ, selecttw2)),
                    mxConstraint(omxSelectRowsAndCols(expCovDZ, selecttw1) ==
                                                omxSelectRowsAndCols(expCovMZ, selecttw1)),
 
                    mxMatrix(type = 'Full', nrow = 1, ncol = nv,
                                     free = TRUE, values = c(21, 21), name = 'Mean'),
                    mxAlgebra(cbind(Mean, Mean), name = 'expMeans'),
 
                    mxModel(name = 'MZ',
                                    mxData(type = 'raw', observed = mzdata),
                                    mxExpectationNormal('Saturated.expCovMZ', 'Saturated.expMeans',
                                                                            dimnames = selvars),
                                    mxFitFunctionML()),
                    mxModel(name = 'DZ',
                                    mxData(type = 'raw', observed = dzdata),
                                    mxExpectationNormal('Saturated.expCovDZ', 'Saturated.expMeans',
                                                                            dimnames = selvars),
                                    mxFitFunctionML()),
                    mxFitFunctionMultigroup(c('MZ', 'DZ'))
    )
}
Top
Tue, 07/17/2018 - 13:10
AdminRobK's picture
AdminRobK
Offline
Joined: 01/24/2014 - 12:15
To answer your main question
So the question is: what is the better way to introduce the constrains given that the constrained elements are computed, not defined as matrices?

I'd have to think about it, but I'm not sure there are any alternatives to MxConstraints if you're dead-set on using the Cholesky parameterization. If I were you, I'd give up on the Cholesky parameterization. It's not mandatory; OpenMx has safeguards against model-expected covariance matrices going non-positive-definite.

Top
Fri, 05/31/2019 - 17:00
AdminRobK's picture
AdminRobK
Offline
Joined: 01/24/2014 - 12:15
follow-up

I just wanted to follow-up on this thread. As of OpenMx version 2.13, CSOLNP still runs as slow as molasses on the model in the OP's script, but it at least ephemerally prints progress reporting to the console, and more importantly, it is now user-interruptible.

We now know what causes CSOLNP so much trouble: there are redundant equality MxConstraints in the model, due to the symmetry of the covariance matrices. If we get rid of the redundancies, CSOLNP gets the solution on one try. See the attached script.

My mxVersion():
OpenMx version: 2.13.2.3 [GIT v2.13.2-3-g3d8b302d-dirty]
R version: R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu
Default optimizer: CSOLNP
NPSOL-enabled?: Yes
OpenMP-enabled?: Yes

File attachments: 
Binary Data 190315--CSOLNP_hangs3.R
Top

Copyright © 2007-2024 The OpenMx Project