Hi, everyone.
I am a new user of OpenMx. I learned a lot from this forum. But I still often bump into difficulties. Here is one.
I write the following commands to model bivariate (binary variables: smoke and drink) ACE Cholesky model.
qlss <- read.table("qlss.csv", header=TRUE, sep=',')
nvar <- 2 #number of varibles
tnvar <- 4 #number of variblesmax family size
nlower <-nvar(nvar+1)/2 #number of free elements in a lower matrix nvar*nvar
maxcat <- 1
Vars <- c('smoke','drink')
selVars <- paste(Vars,c(rep(1,nvar),rep(2,nvar)),sep="")
this yeilds a list of the 8 varibles for analysis
smoke1 drink1 smoke2 drink2
selVars are binary variables
mzData <- subset(qlss, zyg==1, selVars)
dzData <- subset(qlss, zyg==2, selVars)
tell mx the data is ordinal
---------------------------------------------------------------------
mzData$smoke1=mxFactor(mzData$smoke1, levels= c(0:maxcat))
mzData$smoke2=mxFactor(mzData$smoke2, levels= c(0:maxcat))
dzData$smoke1=mxFactor(dzData$smoke1, levels= c(0:maxcat))
dzData$smoke2=mxFactor(dzData$smoke2, levels= c(0:maxcat))
mzData$drink1=mxFactor(mzData$drink1, levels= c(0:maxcat))
mzData$drink2=mxFactor(mzData$drink2, levels= c(0:maxcat))
dzData$drink1=mxFactor(dzData$drink1, levels= c(0:maxcat))
dzData$drink2=mxFactor(dzData$drink2, levels= c(0:maxcat))
Fit Multivariate ACE Model with RawData and Matrices Input
ACE_Cholesky_Model <- mxModel("ACE_Cholesky",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="a" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE,values=.4, name="c" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="e" ),
# Matrices A, C, and E compute variance components
mxAlgebra( a %% t(a), name="A" ),
mxAlgebra( c %% t(c), name="C" ),
mxAlgebra( e %% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nvar, ncol=nvar, name="I"),
mxAlgebra( solve(sqrt(IV)), name="iSD"),
mxAlgebra( solve(sqrt(IA)%%A%%solve(sqrt(IA))), name="Acorrelation"),
mxAlgebra( solve(sqrt(IC)%%C%%solve(sqrt(IC))), name="Ccorrelation"),
mxAlgebra( solve(sqrt(IE)%%E%%solve(sqrt(IE))), name="Ecorrelation"),
## Note that the rest of the mxModel statements do not change for bi/multivariate case
# Matrix & Algebra for expected means vector
mxMatrix( type="Full", nrow=1, ncol=nvar, free=TRUE, values= 0, name="Mean" ),
mxAlgebra( cbind(Mean,Mean), name="expMean"),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra( rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" ),
# Algebra for expected variance/covariance matrix in DZ
mxAlgebra( rbind ( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
)
,
mxModel("MZ",
mxData( observed=mzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars )
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars )
),
mxAlgebra( MZ.objective + DZ.objective, name="modelfit" ),
mxAlgebraObjective("modelfit"),
mxCI(c("ACE.a","ACE.c","ACE.e"))
)
ACE_Cholesky_Fit <- mxRun(ACE_Cholesky_Model)
ACE_Cholesky_Summ <- summary(ACE_Cholesky_Fit)
ACE_Cholesky_Summ
It worked. But when I tried to drop 'c' or 'c21' & 'c22', in the following commands: (For example,)
AE_Cholesky_Model <- mxModel(ACE_Cholesky_Fit, name="AE_Cholesky",
mxModel(ACE_Cholesky_Fit$ACE,
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=c(T,F,F), values=0, name="c" ) # drop c21 c22 at 0
))
AE_Cholesky_Fit <- mxRun(AE_Cholesky_Model)
the 'error' presented:
The algebra 'ACE.Ccorrelation' in model 'AE_Cholesky' generated the error message: (another sentence was in Chinese) Lapack例行程序dgesv: 系统正好是奇异的. It seems that there were no solution.
I don't understand. What should I do? Any help will be appreciated!
And another problem is that: after the following command
summary(ACE_Cholesky_Fit)
, there were no CI information, even the 'Std. Error' were presented. And, sometimes, 'Std. Error' would be 'Na'. Anyone could tell me what happened?
Printer-friendly version
I think the 'error' I think is that a matrix is not positive definite.
Just jiggle the starting values in your a, c, e, matrices
so instead of
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="a" ),
do something like
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=T, values=c(.41,.42,.39,.38), name="a" ),
this will also fix the flow-on error from summary
Thanks for reply.
I changed into:
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=T, values=c(.41,.42,.39), name="a" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=T, values=c(.41,.42,.39), name="c" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=T, values=c(.41,.42,.39), name="e" ),
but there was still the same error, and no CIs.
Can you generate some simulated data that will replicate the error? We can help you debug the issue if we have some data that reproduces the problem. Here is some help on generating simulated data: http://openmx.psyc.virginia.edu/wiki/generating-simulated-data
Is fakeData() a R function? I can't use it in R 2.12.1. And it seems not an OpenMx function.
FakeData is attached at the bottom of the page about generating fake data. It's an R function you can download and run, or copy-and-paste into your own script.
If you want to have OpenMx calculate the confidence intervals you specified, you'll need to make sure you add
intervals=TRUEto yourmxRun()command. You also won't get standard errors or confidence interval estimates if you don't get a converged solution--that is, if OpenMx returns and error or a code -1 fit. We do that because these take time to calculate (especially CIs), and are not meaningful if you're not at a converged solution.Like Ryne said, the "no solution" error comes because of this line:
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=c(T,F,F), values=0, name="c" ) # drop c21 c22 at 0Here, you explicitly set one of the rows of c to all zeros. I*C therefore has a bottom row of all zeros, is not of full rank, and thus is not invertible. When you try to invert the resulting non-invertable matrix, the error gets thrown.
As a temporary solution, just comment out the CCorrelation algebra and see if that gets rid of the error. If you need the CCorrelation result, you'll want to calculate it differently.
I deleted the following commands in ACE model:
mxAlgebra( solve(sqrt(IA)%%A%%solve(sqrt(IA))), name="Acorrelation"),
mxAlgebra( solve(sqrt(IC)%%C%%solve(sqrt(IC))), name="Ccorrelation"),
mxAlgebra( solve(sqrt(IE)%%E%%solve(sqrt(IE))), name="Ecorrelation"),
then it worked out. Thank you very much. hiahia~~ And CI problems were fixed too.
I'm going to guess that your problem has to do with the inversion in the Ccorrelation algebra. When you rebuild the c matrix for the AE model, you're creating a matrix of all zeros. You do a few operations involving it, yielding sqrt(C %%t(C)) %% C %% sqrt(C%% t(C)) once you remove all of the identity matrices that are in there for some reason. When C is all zero, that algebra is a nvar x nvar zero matrix. When you try to invert that matrix using solve(), you should get an error, as a zero matrix isn't positive definite and can't be inverted. If you have any zeros along the diagonal of C, you'll hit the same problem. I think you're going to have to respecify your AE model, either without a C matrix or with another correction to alleviate this problem.
As a note for future submissions, either attach code in a separate file, surround code blocks with an html tag like "code" or "pre" (which would look like < code > your code < /code > without the spaces) or just put a bunch of extra line breaks. The comments you made between your ACE and AE models helped me make a diagnosis, but it took me a few minutes to figure out this was two separate models with comments in the middle and not one giant model.
should I just drop
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=F,values=0, name="c" ),
or/and
mxAlgebra(c%*%t(c),name="C"),
? I think either is wrong.
Sorry, I am not very good at making comments.
Did you mean this:
Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="a" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=F,values=0, name="c" ),
mxMatrix( type="Lower", nrow=nvar, ncol=nvar, free=TRUE, values=.4, name="e" ),
it ran out this:
error: unknown expected covariance name 'ACE.expCovMZ' detected in the objective function of model 'MZ'
Hi, I could obtain a, c, and e by the command:
mxCI(c("ACE.a","ACE.c","ACE.e"))
, but when I write this:
mxCI(c("ACE.a %*% ACE.iSD")), or
mxCI(c("ACE.A/ACE.V"))
, the error appeared:
"Unknown reference to 'ACE.a*ACE.iSD' (or 'ACE.A/ACE.V')detected in a confidence interval specification in model 'AE_Cholesky' in runHelper (model, frontendStart, intervals, silent, suppressWarnings, unsafe, checkpoint, useSocket, onlyFrontend)"
. Is it possible to obtain CIs by this:
ACEpathMatrices <- c("ACE.a","ACE.c","ACE.e","ACE.iSD","ACE.iSD %% ACE.a","ACE.iSD %% ACE.c","ACE.iSD %*% ACE.e")
ACEpathLabels <- c("pathEst_a","pathEst_c","pathEst_e","sd","stPathEst_a","stPathEst_c","stPathEst_e")
formatOutputMatrices(ACE_Cholesky_Fit,ACEpathMatrices,ACEpathLabels,Vars,4)
? Thank you in advance.
First add mxAlgebra() expressions that compute the desired expression, and use mxCI() on the names that you assign to the new algebras. For example, add the mxAlgebra() inside the construction of the "ACE" model:
And then you can add mxCI("ACE.thing1") to the parent model "ACE_Cholesky".
Thank you for the explanation. I ran this:
mxAlgebra(a %*% iSD, name="sta"),
and obtained standardized a. But the sta[2,1] is different from the a21 which was computed from this:
ACEpathMatrices <- c("ACE.iSD %% ACE.a","ACE.iSD %% ACE.c","ACE.iSD %*% ACE.e")
ACEpathLabels <- c("stPathEst_a","stPathEst_c","stPathEst_e")
formatOutputMatrices(ACE_Cholesky_Fit,ACEpathMatrices,ACEpathLabels,4)
. The results is .39 (.34-.44) for the upper one, and .45 for the lower one. What is the matter?
mxAlgebra(ACE.a %*% ACE.iSD, name="sta"), sta[2,1] # = .39 (.34-.44) ACEpathMatrices = c("ACE.iSD %*% ACE.a") a21 # = .45 for the lower oneI see.
There is another question. When I drop parameters, one submodel (after mxRun()) showed this error:
The job for model 'AE_Cholesky' exited abnormally with the error message: Objective function returned an infinite value.
How can I deal with it? (Other similar type of submodels which dropped different parameters did run.)