You are here

mxConstraint

18 posts / 0 new
Last post
tbates's picture
Offline
Joined: 07/31/2009 - 14:25
mxConstraint

a place to discuss mxConstraint

mkeller's picture
Offline
Joined: 08/04/2009 - 16:02
Hi all, I'm trying to write a

Hi all,

I'm trying to write a nuclear twin family (NTF) design model in openMx but am running into immediate problems using mxConstraint. The NTF requires several nonlinear constraints, such that, e.g., a parameter will be on both sides of the equation. For example, the phenotypic variance is itself partly a function of the phenotypic variance (through the familial environment). I came up with what I thought was a clever way of doing this in openMx but openMx doesn't like it. Here's what I did:

Fit NTF Model with RawData and Matrices Input

ntf <- mxModel(model="NucTwinFam",
# Matrices
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="Env", name="e"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="FamilialVar", name="f"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="FamilialPath", name="m"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="AddGen", name="a"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="Sib", name="s"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="Dominance", name="d"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="AMCopath", name="mu"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="VarPhen", name="Vp1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="VarF", name="x1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="CovPhenGen", name="delta1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="VarAddGen", name="q1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.4, label="CovFA", name="w1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=2, label="two", name="Two"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="one", name="One"),

#mxAlgebra section
mxAlgebra(expression= Vp1 %% mu %% t(Vp1), name="CvSps"),
mxAlgebra(expression= a %% q1 %% t(a) + f %% x1 %% t(f) + Two %x% a %% w1 %% t(f) + e %% t(e) + d %% t(d) + s %% t(s), name="Vp2"),
mxAlgebra(expression= Two %x% m %
% Vp1 %% t(m) + Two %x% m %% Vp1 %% mu %% t(Vp1) %% t(m), name="x2"),
mxAlgebra(expression= q1 %
% a + w1 %% f, name="delta2"),
mxAlgebra(expression= One + delta1 %
% mu %% t(delta1), name="q2"),
mxAlgebra(expression= delta1 %
% m + delta1 %% mu %% Vp1 %*% t(m), name="w2"),

#constraints
mxConstraint(Vp1,"=",Vp2),
mxConstraint(x1,"=",x2),
mxConstraint(delta1,"=",delta2),
mxConstraint(q1,"=",q2),
mxConstraint(w1,"=",w2))

Here's what I got from R:

Error in typeof(alg1) : object 'Vp1' not found

Any help would be much appreciated. Best...
Matt

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
Try: mxConstraint("Vp1","=","

Try:

mxConstraint("Vp1","=","Vp2"),
mxConstraint("x1","=","x2"),
mxConstraint("delta1","=","delta2"),
mxConstraint("q1","=","q2"),
mxConstraint("w1","=","w2")
mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
As a follow up, the root

As a follow up, the root cause of this error is that algebraic expressions can be written without quotes (only used in mxAlgebra and mxEval), but the name of a named entity must always appear as a string. We're doing magic trickery to parse algebraic expressions. To maintain consistency we could require algebraic expressions to be written as strings. But most users would find this option to be less convenient than the current system.

Note that R provides a similar cryptic error for the following expresssion:

> c(sum(mean(colSums(Vp1))))
Error in inherits(x, "data.frame") : object 'Vp1' not found
mkeller's picture
Offline
Joined: 08/04/2009 - 16:02
Hi Michael and all, Thanks

Hi Michael and all,

Thanks for the help! That got me further, almost to the end. I'm now trying to run the full NTF script and get the following error from R:

> NTFModelFit <- mxRun(NTFModel)
Running NTF
Error: Unknown reference: 'MZNTF.4'

I've looked in the mxModel MZNTF and in the mxModel MZNTF calls but couldn't find any object '4' or anything that could become a '4'.

This may or may not have anything to do with mxConstraint, but I decided to post here because this is the same script as before and because I suspect it does have to do with constraints. Here's the entire script along with simulated data. It all 'works' until the mxRun stage. I'm sure it is (once again) something that isn't too difficult. Any help would be appreciated:

simulate null data

mzf <- matrix(rnorm(400),ncol=4)
dzf <- matrix(rnorm(400),ncol=4)

Fit NTF Model with RawData and Matrices Input

ntf <- mxModel(model="NucTwFam",
# Matrices
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.7, label="Env", name="e"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="FamilialVar", name="f"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.1, label="FamilialPath", name="m"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="AddGen", name="a"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.3, label="Sib", name="s"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2, label="Dominance", name="d"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.1, label="AMCopath", name="mu"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=1.2, label="VarPhen", name="Vp1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.05, label="VarF", name="x1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.25, label="CovPhenGen", name="delta1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.35, label="VarAddGen", name="q1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.15, label="CovFA", name="w1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=2, label="two", name="Two"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="one", name="One"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=.5, label="half", name="Half"),
#mxAlgebra section - nonlinear constraints
mxAlgebra(expression= a %% q1 %% t(a) + f %% x1 %% t(f) + Two %x% a %% w1 %% t(f) + e %% t(e) + d %% t(d) + s %% t(s), name="Vp2"),
mxAlgebra(expression= Two %x% m %
% Vp1 %% t(m) + Two %x% m %% Vp1 %% mu %% t(Vp1) %% t(m), name="x2"),
mxAlgebra(expression= q1 %
% a + w1 %% f, name="delta2"),
mxAlgebra(expression= One + delta1 %
% mu %% t(delta1), name="q2"),
mxAlgebra(expression= delta1 %
% m + delta1 %% mu %% Vp1 %% t(m), name="w2"),
#constraints - equating nonlinear constraints and parameters
mxConstraint('Vp1',"=",'Vp2',name='VpCon'),
mxConstraint('x1',"=",'x2',name='xCon'),
mxConstraint('delta1',"=",'delta2',name='deltaCon'),
mxConstraint('q1',"=",'q2',name='qCon'),
mxConstraint('w1',"=",'w2',name='wCon'),
#mxAlgebra section - relative covariances
mxAlgebra(expression= a %
% q1 %% t(a) + f %% x1 %% t(f) + Two %x% a %% w1 %% t(f) + d %% t(d) + s %% t(s), name="CvMz"),
mxAlgebra(expression= a %
% (q1-Half) %% t(a) + Half%x%Half%x% d%%t(d) + f %% x1 %% t(f) + Two %x% a %% w1 %% t(f) + s %% t(s), name="CvDz"),
mxAlgebra(expression= Half %x% a %
% (q1 %% a + w1 %% f) + Half %x% a %% (q1 %% a + w1 %% f) %% mu %% Vp1 + m %% Vp1 + m %% Vp1 %% mu %% t(Vp1), name="ParChild"),
mxAlgebra(expression= Vp1 %
% mu %*% t(Vp1), name="CvSps"))

mzModel <- mxModel(ntf, name = "MZNTF",
mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(mzf)), name="expMeanMz"),
# Algebra for expected variance/covariance matrix in MZF
mxAlgebra(expression=matrix(
c(Vp1, CvMz, ParChild, ParChild,
CvMz, Vp1, ParChild, ParChild,
ParChild, ParChild, Vp1, CvSps,
ParChild, ParChild, CvSps, Vp1), nrow=4,byrow=TRUE), dimnames=list(colnames(mzf),colnames(mzf)),name="expCovMz"),
mxData(observed=mzf, type="raw"),
mxFIMLObjective(covariance="expCovMz",means="expMeanMz"))

dzModel <- mxModel(ntf, name = "DZNTF",
mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(mzf)), name="expMeanDz"),
# Algebra for expected variance/covariance matrix in MZF
mxAlgebra(expression=matrix(
c(Vp1, CvDz, ParChild, ParChild,
CvDz, Vp1, ParChild, ParChild,
ParChild, ParChild, Vp1, CvSps,
ParChild, ParChild, CvSps, Vp1), nrow=4,byrow=TRUE), dimnames=list(colnames(mzf),colnames(mzf)),name="expCovDz"),
mxData(observed=dzf, type="raw"),
mxFIMLObjective(covariance="expCovDz",means="expMeanDz"))

NTFModel <- mxModel(model="NTF", mzModel, dzModel,
mxAlgebra(expression=MZNTF.objective + DZNTF.objective, name="ntffit"), #MZ.objective is the automatic name for the -2LL of mzModel
mxAlgebraObjective("ntffit"))

NTFModelFit <- mxRun(NTFModel)

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
The error is informing you

The error is informing you that have a matrix or an algebra named "4" inside the model named MZNTF. The culprit is the following declaration:

mxAlgebra(expression=matrix(
c(Vp1, CvMz, ParChild, ParChild,
CvMz, Vp1, ParChild, ParChild,
ParChild, ParChild, Vp1, CvSps,
ParChild, ParChild, CvSps, Vp1), nrow=4,byrow=TRUE)

It looks like you've confusing the mxMatrix() and mxAlgebra() notation. mxAlgebra takes 2 arguments, and mxModel takes nrow, ncol, byrow, etc. And at the moment, you can't use 1 x 1 matrices or algebras inside the declarations of mxMatrix() expressions. See ?mxAlgebra for a list of the matrix operations that are supported.

mkeller's picture
Offline
Joined: 08/04/2009 - 16:02
Thanks as always Michael.

Thanks as always Michael. It's helpful to at least know where to start looking (btw, is there any way to figure out where such an error is coming from? I've used traceback() with no success).

I should note, however, that the nrow= argument I used there is an argument to the base R matrix() function, not the mxAlgebra() or mxMatrix() functions. I thought that you could put a matrix into the expression argument of mxAlgebra().

I'll redo the matrix using rbind and cbind I guess, but it seems like you should be able to put a matrix into the expression argument of mxAlgebra, given that I've seen other scripts use rbind(cbind()) to basically make a matrix on the fly and put that into expresion= for mxAlgebra.

Best,
Matt

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
I recommend writing your

I recommend writing your script using the following style:

foo <- mxMatrix(....)
bar <- mxMatrix(....)
baz <- mxMatrix(....)
quux <- mxAlgebra(....)
quuux <- mxAlgebra(...)
...
...
model <- mxModel('model name', foo, bar, baz, quux, quuux)

If there is an error in one of the calls to mxMatrix() or mxAlgebra(), then the error should be easier to find using traceback(). Even if the error doesn't appear until mxModel() or mxRun(), you should be able to inspect each matrix or algebra by just evaluating "foo" or "print(foo)" and you can double-check that it's storing what you expect it to store.

WRT your comments, we have restricted mxAlgebra() expressions only to functions or operations that operate on matrices. Such a restriction greatly simplifies the back-end processing. This is why the matrix() function is not allowed, because the arguments to matrix() are scalars. Similarly cbind() and rbind() are allowed, because the arguments to these functions are matrices. *As of version 0.1.4 you can write numeric constants in mxAlgebra() statements. These constants are converted into 1x1 matrices. *

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
Matt, I'm checking in a patch

Matt, I'm checking in a patch that gives some more useful information for that error message. I ran it on your script and got: Error: Unknown reference 'TRUE' detected in the entity 'expCovMz' in model 'MZNTF'. That's more helpful. The underlying cause is the matrix() expression in the mxAlgebra() declaration.

mkeller's picture
Offline
Joined: 08/04/2009 - 16:02
Hi all, I'm still running

Hi all,

I'm still running into problems with my attempt at a nuclear twin family model using OpenMx. The script has a lot of rather complicated mxConstraints, so I think this is the issue here. Here is the error:

> NTFModelFit <- mxRun(NTFModel)
Running NTF

*** caught bus error ***
address 0xc, cause 'non-existent physical address'

R doesn't crash after this but hangs indefinitely (eating about 50% CPU), and will not return control to the user, so I can't try traceback() or other methods of debugging. Thanks in advance!

Once again, here is the self-contained script:

simulate data

mzf <- matrix(rnorm(400),ncol=4)
dimnames(mzf) <- list(NULL,c("Tw1","Tw2","Fa","Mo"))
dzf <- matrix(rnorm(400),ncol=4)
dimnames(dzf) <- list(NULL,c("Tw1","Tw2","Fa","Mo"))

Fit NTF Model with RawData and Matrices Input

ntf <- mxModel(model="NucTwFam",
# Matrices
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.7, label="Env", name="e"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="FamilialVar", name="f"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.1, label="FamilialPath", name="m"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="AddGen", name="a"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.3, label="Sib", name="s"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2, label="Dominance", name="d"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.1, label="AMCopath", name="mu"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=1.2, label="VarPhen", name="Vp1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.05, label="VarF", name="x1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.25, label="CovPhenGen", name="delta1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.35, label="VarAddGen", name="q1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.15, label="CovFA", name="w1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=2, label="two", name="Two"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="one", name="One"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=.5, label="half", name="Half"),
#mxAlgebra section - nonlinear constraints
mxAlgebra(expression= a %% q1 %% t(a) + f %% x1 %% t(f) + Two %x% a %% w1 %% t(f) + e %% t(e) + d %% t(d) + s %% t(s), name="Vp2"),
mxAlgebra(expression= Two %x% m %
% Vp1 %% t(m) + Two %x% m %% Vp1 %% mu %% t(Vp1) %% t(m), name="x2"),
mxAlgebra(expression= q1 %
% a + w1 %% f, name="delta2"),
mxAlgebra(expression= One + delta1 %
% mu %% t(delta1), name="q2"),
mxAlgebra(expression= delta1 %
% m + delta1 %% mu %% Vp1 %% t(m), name="w2"),
#constraints - equating nonlinear constraints and parameters
mxConstraint('Vp1',"=",'Vp2',name='VpCon'),
mxConstraint('x1',"=",'x2',name='xCon'),
mxConstraint('delta1',"=",'delta2',name='deltaCon'),
mxConstraint('q1',"=",'q2',name='qCon'),
mxConstraint('w1',"=",'w2',name='wCon'),
#mxAlgebra section - relative covariances
mxAlgebra(expression= a %
% q1 %% t(a) + f %% x1 %% t(f) + Two %x% a %% w1 %% t(f) + d %% t(d) + s %% t(s), name="CvMz"),
mxAlgebra(expression= a %
% (q1-Half) %% t(a) + Half%x%Half%x% d%%t(d) + f %% x1 %% t(f) + Two %x% a %% w1 %% t(f) + s %% t(s), name="CvDz"),
mxAlgebra(expression= Half %x% a %
% (q1 %% a + w1 %% f) + Half %x% a %% (q1 %% a + w1 %% f) %% mu %% Vp1 + m %% Vp1 + m %% Vp1 %% mu %% t(Vp1), name="ParChild"),
mxAlgebra(expression= Vp1 %
% mu %*% t(Vp1), name="CvSps"))

mzModel <- mxModel(ntf, name = "MZNTF",
mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(mzf)), name="expMeanMz"),
# Algebra for expected variance/covariance matrix in MZF
mxAlgebra(expression=rbind(
cbind(Vp1, CvMz, ParChild, ParChild),
cbind(CvMz, Vp1, ParChild, ParChild),
cbind(ParChild, ParChild, Vp1, CvSps),
cbind(ParChild, ParChild, CvSps, Vp1)), dimnames=list(colnames(mzf),colnames(mzf)),name="expCovMz"),
mxData(observed=mzf, type="raw"),
mxFIMLObjective(covariance="expCovMz",means="expMeanMz"))

dzModel <- mxModel(ntf, name = "DZNTF",
mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(mzf)), name="expMeanDz"),
# Algebra for expected variance/covariance matrix in MZF
mxAlgebra(expression=rbind(
cbind(Vp1, CvMz, ParChild, ParChild),
cbind(CvMz, Vp1, ParChild, ParChild),
cbind(ParChild, ParChild, Vp1, CvSps),
cbind(ParChild, ParChild, CvSps, Vp1)), dimnames=list(colnames(mzf),colnames(mzf)),name="expCovDz"),
mxData(observed=dzf, type="raw"),
mxFIMLObjective(covariance="expCovDz",means="expMeanDz"))

NTFModel <- mxModel(model="NTF", mzModel, dzModel,
mxAlgebra(expression=MZNTF.objective + DZNTF.objective, name="ntffit"), #MZ.objective is the automatic name for the -2LL of mzModel
mxAlgebraObjective("ntffit"))

NTFModelFit <- mxRun(NTFModel)

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
Hi Matt, I checked in a fix

Hi Matt,

I checked in a fix to the repository that (I think) fixed the bug. I'm still opening a ticket because there's some hinky code in the constraint assignments in the backend. Are you comfortable with "svn" and "make"? Because you can grab the latest version, if you don't want to wait until the next binary release.

mkeller's picture
Offline
Joined: 08/04/2009 - 16:02
hi Michael, I know what svn

hi Michael,

I know what svn and make do but have never used them. I'm willing to try if you send the code to do it (no worries if this is a hassle). Best,
Matt

mkeller's picture
Offline
Joined: 08/04/2009 - 16:02
Success at last! Thanks for

Success at last! Thanks for working with me Michael.

I downloaded the latest version of OpenMx and used 'svn' and 'make' to create the package. That took care of the "*** caught bus error *** address 0xc, cause 'non-existent physical address' " error message. Next I ran the script below using simulated data (from GeneEvolve). The OpenMx results agree well with the simulated results, as can be seen here:

> res.mat
Var.E Var.F Var.A Var.S Var.D Cor.Sps Var.Phen
OpenMx-Estimated 0.310 0 0.421 0.259 0.049 0.189 1.065
Simulated 0.305 0 0.432 0.200 0.099 0.202 1.036

In case anyone's interested, here is the final NTF script (the part up top is just used for simulating the NTF data). Note that there might be a minor error somewhere in the algebra as the Var(D) is slightly underestimated and Var(S) is slightly overestimated:

NTF design in OpenMx by Matt Keller & Sarah Medland

require(OpenMx)

Run Gene Evolve

setwd("/temp/OMX")
source("GE-73.R")

nms <- c("Famid","Tw1","Tw2","Fa","Mo","bro1","bro2","sis1","sis2","sp.tw1","sp.tw2",
"son1.tw1","son2.tw1","dau1.tw1","dau2.tw1", "son1.tw2","son2.tw2","dau1.tw2","dau2.tw2",
"Tw1.age","Tw2.age","Fa.age","Mo.age","bro1.age","bro2.age","sis1.age","sis2.age","sp.tw1.age",
"sp.tw2.age", "son1.tw1.age","son2.tw1.age","dau1.tw1.age","dau2.tw1.age", "son1.tw2.age",
"son2.tw2.age","dau1.tw2.age","dau2.tw2.age")
mzf <- as.matrix(read.table("MZF",header=FALSE,col.names=nms))
mzf <- mzf[,2:5]
mzf[mzf==-999] <- NA
dzf <- as.matrix(read.table("DZF",header=FALSE,col.names=nms))
dzf <- dzf[,2:5]
dzf[dzf==-999] <- NA

Fit NTF Model with RawData and Matrices Input

ntf <- mxModel(model="NucTwFam",
# Matrices
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.7, label="Env", name="e"),
mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="FamilialVar", name="f"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.1, label="FamilialPath", name="m"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="AddGen", name="a"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.3, label="Sib", name="s"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2, label="Dominance", name="d"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.1, label="AMCopath", name="mu"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=1.2, label="VarPhen", name="Vp1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.2, label="VarF", name="x1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.25, label="CovPhenGen", name="delta1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=1, label="VarAddGen", name="q1"),
mxMatrix(type="Full", nrow=1, ncol=1, free=TRUE, values=.15, label="CovFA", name="w1"),
#mxAlgebra section - nonlinear constraints
mxAlgebra(expression= a %% q1 %% t(a) + f %% x1 %% t(f) + 2 %x% a %% w1 %% t(f) + e %% t(e) + d %% t(d) + s %% t(s), name="Vp2"),
mxAlgebra(expression= 2 %x% m %
% Vp1 %% t(m) + 2 %x% m %% Vp1 %% mu %% t(Vp1) %% t(m), name="x2"),
mxAlgebra(expression= q1 %
% a + w1 %% f, name="delta2"),
mxAlgebra(expression= 1 + delta1 %
% mu %% t(delta1), name="q2"),
mxAlgebra(expression= delta1 %
% m + delta1 %% mu %% Vp1 %% t(m), name="w2"),
#constraints - equating nonlinear constraints and parameters
mxConstraint('Vp1',"=",'Vp2',name='VpCon'),
mxConstraint('x1',"=",'x2',name='xCon'),
mxConstraint('delta1',"=",'delta2',name='deltaCon'),
mxConstraint('q1',"=",'q2',name='qCon'),
mxConstraint('w1',"=",'w2',name='wCon'),
#mxAlgebra section - relative covariances
mxAlgebra(expression= a %
% q1 %% t(a) + f %% x1 %% t(f) + 2 %x% a %% w1 %% t(f) + d %% t(d) + s %% t(s), name="CvMz"),
mxAlgebra(expression= a %
% (q1-.5) %% t(a) + .25 %x% d %% t(d) + f %% x1 %% t(f) + 2 %x% a %% w1 %% t(f) + s %% t(s), name="CvDz"),
mxAlgebra(expression= .5 %x% a %
% (q1 %% a + w1 %% f) + .5 %x% a %% (q1 %% a + w1 %% f) %% mu %% Vp1 + m %% Vp1 + m %% Vp1 %% mu %% t(Vp1), name="ParChild"),
mxAlgebra(expression= Vp1 %
% mu %*% t(Vp1), name="CvSps"))

mzModel <- mxModel(ntf, name = "MZNTF",
mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(mzf)), name="expMeanMz"),
# Algebra for expected variance/covariance matrix in MZF
mxAlgebra(expression=rbind(
cbind(Vp1, CvMz, ParChild, ParChild),
cbind(CvMz, Vp1, ParChild, ParChild),
cbind(ParChild, ParChild, Vp1, CvSps),
cbind(ParChild, ParChild, CvSps, Vp1)), dimnames=list(colnames(mzf),colnames(mzf)),name="expCovMz"),
mxData(observed=mzf, type="raw"),
mxFIMLObjective(covariance="expCovMz",means="expMeanMz"))

dzModel <- mxModel(ntf, name = "DZNTF",
mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values= .25, label="mean", dimnames=list(NULL, colnames(dzf)), name="expMeanDz"),
# Algebra for expected variance/covariance matrix in MZF
mxAlgebra(expression=rbind(
cbind(Vp1, CvDz, ParChild, ParChild),
cbind(CvDz, Vp1, ParChild, ParChild),
cbind(ParChild, ParChild, Vp1, CvSps),
cbind(ParChild, ParChild, CvSps, Vp1)), dimnames=list(colnames(dzf),colnames(dzf)),name="expCovDz"),
mxData(observed=dzf, type="raw"),
mxFIMLObjective(covariance="expCovDz",means="expMeanDz"))

NTFModel <- mxModel(model="NTF", mzModel, dzModel,
mxAlgebra(expression=MZNTF.objective + DZNTF.objective, name="ntffit"), #MZ.objective is the automatic name for the -2LL of mzModel
mxAlgebraObjective("ntffit"))

Run MX

NTFModelFit <- mxRun(NTFModel)

Look at results

summary(NTFModelFit)

res <- NTFModelFit@output$estimate
round(res,3)

compare to simulation

res.mat <- rbind(round(c(res[1:5]^2,res[6:7]),3),round(ALL$track.changes[c('var.U','var.F','var.A','var.S','var.D','cor.spouses','var.cur.phenotype'),'data.t1'],3))
dimnames(res.mat) <- list(c('OpenMx-Estimated','Simulated'),c('Var.E','Var.F','Var.A','Var.S','Var.D','Cor.Sps','Var.Phen'))

look at implied Covariances

round(NTFModelFit@output$algebras$MZNTF.expCovMz,3)
round(NTFModelFit@output$algebras$DZNTF.expCovDz,3)

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
Great this works now! FYI:

Great this works now!

FYI: You can foreshorten, and I think speed-up, the algebra by replacing statements like:

a %*% b %*% t(a) 

with:

a %&% b 

Which is supported in OpenMx
http://openmx.psyc.virginia.edu/wiki/quadratic-matrix-multiplication

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
hi mat, Do you mind if I

hi mat,
Do you mind if I upload the ntf script to the trunk/demo folder?

Then we can all work on it, and add the necessary old-mx verification etc.

I've played with simplifying the algebra. It seems to me that it is not going to work at the minute, as the MZ and DZ groups are using private copies of the matrices and constraints etc (i.e., the constraints are not being satisifed with the same values in both groups.

I'll upload a version re-written with three groups (an alternative would be to refer to the MZ group matrices and algebra in the DZ group by saying MZ.x instead of just x, etc.)

I think that is why some of the numbers might be a bit off.

Also is the code to create a simulate the family data include-able?

mkeller's picture
Offline
Joined: 08/04/2009 - 16:02
Hi Tim, No, go nuts! Feel

Hi Tim,

No, go nuts! Feel free to change the script. I'm not sure that it's possible that the constraints could arrive at different values in the MZ and DZ groups - the model would be seriously underidentified, and if so, I wouldn't think it could arrive by chance at the right answer. But keep me informed. Maybe I'm wrong.

Re the simulation script, I'd be glad to send it to you, but it's awfully involved. It's an MCMC simulation, downloadable at my website (matthewckeller.com). Let me know if you need help using the script.

Best,

Matt

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
ok its up in trunk/demo/ hit

ok its up in trunk/demo/

hit "svn update" and let me know what you think

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
Hi mat, in terms of getting

Hi mat,
in terms of getting this out of failing into passing, it will be good to start with something we know works. To that end...

Following this comment

mxMatrix("Full",1,1, FALSE, 0, label="FamilialPath", name="m"), # fix m=0 if you want Vf=0

If find that if I fix 'm' at 0, the script doesn't run: not positive definite

Running a script like the one you pasted into this thread, but accessing the data file in the repository, and leaving m free so that it runs, I get results (last line) appended to the table included in your script:

> res.mat

Var.E Var.F Var.A Var.S Var.D Cor.Sps Var.Phen

OpenMx 0.443 0.100 0.294 0.149 0 0.209 0.986

Old.mx 0.442 0.099 0.294 0.149 0 0.206 0.986

Simulated 0.433 0.100 0.304 0.197 0 0.200 1.032

new data file 0.294 0.002 0.507 0.151 0.029 0.212 0.986

That last line is the output I am getting. Do the data you sent have a significant D or were you running a script setting D to zero??

try
mate trunk/models/failing/NTF_design.R

( or open trunk/models/failing/NTF_design.R if you don't use TextMate)