Does it mean 'no solution'?

Posted on
No user picture. userzht Joined: 02/19/2011
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 varibles*max 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(I*V)), name="iSD"),
mxAlgebra( solve(sqrt(I*A)%*%A%*%solve(sqrt(I*A))), name="Acorrelation"),
mxAlgebra( solve(sqrt(I*C)%*%C%*%solve(sqrt(I*C))), name="Ccorrelation"),
mxAlgebra( solve(sqrt(I*E)%*%E%*%solve(sqrt(I*E))), 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?

Replied on Wed, 03/30/2011 - 06:59
Picture of user. tbates Joined: 07/31/2009

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

Replied on Wed, 03/30/2011 - 08:48
No user picture. userzht Joined: 02/19/2011

In reply to by tbates

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.

Replied on Wed, 03/30/2011 - 09:38
Picture of user. mspiegel Joined: 07/31/2009

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
Replied on Wed, 03/30/2011 - 10:26
Picture of user. tbrick Joined: 07/31/2009

In reply to by userzht

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=TRUE to your mxRun() 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 0

Here, 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.

Replied on Wed, 03/30/2011 - 20:10
No user picture. userzht Joined: 02/19/2011

In reply to by tbrick

I deleted the following commands in ACE model:

mxAlgebra( solve(sqrt(I*A)%*%A%*%solve(sqrt(I*A))), name="Acorrelation"),
mxAlgebra( solve(sqrt(I*C)%*%C%*%solve(sqrt(I*C))), name="Ccorrelation"),
mxAlgebra( solve(sqrt(I*E)%*%E%*%solve(sqrt(I*E))), name="Ecorrelation"),

then it worked out. Thank you very much. hiahia~~ And CI problems were fixed too.

Replied on Wed, 03/30/2011 - 09:43
Picture of user. Ryne Joined: 07/31/2009

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.

Replied on Wed, 03/30/2011 - 10:19
No user picture. userzht Joined: 02/19/2011

In reply to by Ryne

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'

Replied on Wed, 03/30/2011 - 22:16
No user picture. userzht Joined: 02/19/2011

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.

Replied on Wed, 03/30/2011 - 23:20
Picture of user. mspiegel Joined: 07/31/2009

In reply to by userzht

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:

mxAlgebra(a %*% iSD, name = "thing1")

And then you can add mxCI("ACE.thing1") to the parent model "ACE_Cholesky".

Replied on Thu, 03/31/2011 - 01:17
No user picture. userzht Joined: 02/19/2011

In reply to by mspiegel

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?

Replied on Thu, 03/31/2011 - 06:03
No user picture. userzht Joined: 02/19/2011

In reply to by tbates

I 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.)