Does it mean 'no solution'?
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?
jiggle start values to avoid error
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
Log in or register to post comments
In reply to jiggle start values to avoid error by tbates
changed, but all the same
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.
Log in or register to post comments
generate simulated data
Log in or register to post comments
In reply to generate simulated data by mspiegel
fakeData()
Log in or register to post comments
In reply to fakeData() by userzht
fakeData(). Also, clarification of the error
If you want to have OpenMx calculate the confidence intervals you specified, you'll need to make sure you add
intervals=TRUE
to 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 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.
Log in or register to post comments
In reply to fakeData(). Also, clarification of the error by tbrick
worked out ;)
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.
Log in or register to post comments
I'm going to guess that your
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.
Log in or register to post comments
In reply to I'm going to guess that your by Ryne
I am not sure whether I understand.
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.
Log in or register to post comments
In reply to I'm going to guess that your by Ryne
I am not sure I understand.
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'
Log in or register to post comments
How can I get the CIs of standardized path estimate a, c, and e?
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.
Log in or register to post comments
In reply to How can I get the CIs of standardized path estimate a, c, and e? by userzht
First add mxAlgebra()
And then you can add mxCI("ACE.thing1") to the parent model "ACE_Cholesky".
Log in or register to post comments
In reply to First add mxAlgebra() by mspiegel
different standardized a21?
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?
Log in or register to post comments
In reply to different standardized a21? by userzht
a %*% b != b %*% a
Log in or register to post comments
In reply to a %*% b != b %*% a by tbates
So it is. ^_^
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.)
Log in or register to post comments