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?