Does it mean 'no solution'?

14 replies [Last post]
userzht's picture
Offline
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?

userzht's picture
Offline
Joined: 02/19/2011
How can I get the CIs of standardized path estimate a, c, and e?

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.

mspiegel's picture
Offline
Joined: 07/31/2009
First add mxAlgebra()

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".

userzht's picture
Offline
Joined: 02/19/2011
different standardized a21?

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?

tbates's picture
Offline
Joined: 07/31/2009
a %*% b != b %*% a

mxAlgebra(ACE.a %*% ACE.iSD, name="sta"),
sta[2,1] # <i>= .39 (.34-.44)</a>
 
ACEpathMatrices = c("ACE.iSD %*% ACE.a")
a21      # <i>= .45 for the lower one </>

userzht's picture
Offline
Joined: 02/19/2011
So it is. ^_^

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

Ryne's picture
Online
Joined: 07/31/2009
I'm going to guess that your

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.

userzht's picture
Offline
Joined: 02/19/2011
I am not sure I understand.

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'

userzht's picture
Offline
Joined: 02/19/2011
I am not sure whether I understand.

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.

mspiegel's picture
Offline
Joined: 07/31/2009
generate simulated data

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

userzht's picture
Offline
Joined: 02/19/2011
fakeData()

Is fakeData() a R function? I can't use it in R 2.12.1. And it seems not an OpenMx function.

tbrick's picture
Offline
Joined: 07/31/2009
fakeData(). Also, clarification of the error

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.

userzht's picture
Offline
Joined: 02/19/2011
worked out ;)

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.

tbates's picture
Offline
Joined: 07/31/2009
jiggle start values to avoid error

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

userzht's picture
Offline
Joined: 02/19/2011
changed, but all the same

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.