Login

News

Equality constraints in saturated model

Posted on
No user picture. IvanVoronin Joined: 08/18/2013
Forums
Behavioral Genetics Models
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

Replied on Tue, 07/17/2018 - 11:44
Picture of user. AdminRobK Joined: 01/24/2014

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'))
)
}

Log in or register to post comments
Replied on Tue, 07/17/2018 - 12:10
Picture of user. AdminRobK Joined: 01/24/2014

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.

Log in or register to post comments
Replied on Fri, 05/31/2019 - 16:00
Picture of user. AdminRobK Joined: 01/24/2014

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
190315--CSOLNP_hangs3.R
Log in or register to post comments
  Login

News