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:
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
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:
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"))
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.
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.
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. *
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.
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!
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"))
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.
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:
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
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"))
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?
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.
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)
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
Try:
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:
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)
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:
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.
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
I recommend writing your script using the following style:
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. *
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.
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)
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.
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
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)
Great this works now!
FYI: You can foreshorten, and I think speed-up, the algebra by replacing statements like:
with:
Which is supported in OpenMx
http://openmx.psyc.virginia.edu/wiki/quadratic-matrix-multiplication
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?
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
ok its up in trunk/demo/
hit "svn update" and let me know what you think
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)