With the release of OpenMx 1.1, I started playing around with the new joint ordinal-continuous feature. We had just done a bivariate ACE model for two phenotypes (one ordinal, 3 levels, and one continuous, which we divided into 10 levels) using Classic Mx.
At first glance the results look good. The script seems to run without any problems, and the resulting variance component estimates are comparable to what we see in the univariates - there is some C effect on the ordinal variables that wasn't there in the univariate analysis, but this appears to be due to the C effect on the continuous variable.
What I'm concerned about is the movement that we see in the thresholds for the ordinal phenotype. In some cases the movement is rather small, nothing that I wouldn't attribute to simple differences between univariate and multivariate models. There are other cases, however, where the thresholds more quite a lot. For example, in one model, the thresholds are estimated at 2.59 for T1 and 4.64 for T2 while in the univariate model these are estimated at 1.04 for T1 and 1.81 for T2.
Since this is a new script for me, I'm not sure how things are working. Are differences like this to be expected with this script, or do we need to reconsider our start values and boundaries?
Also, I noticed (as I ran this same bivariate ACE model for 5 ordinal variables against the same continuous variable) that sometimes it would give me NaN standard errors for my parameters and sometimes it wouldn't. Is that also something expected in the joint script?
My joint ordinal-continuous bivariate ACE script is below (probably a bit kludgy):
bivACEModel <- mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholASVnv, lbound=-4, ubound=4,labels=AFac, name="a" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholCSVnv, lbound=-4, ubound=4,labels=CFac, name="c" ),
mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholESVnv, lbound=-4, ubound=4,labels=EFac, name="e" ),
# Matrices A, C, and E compute variance components
mxAlgebra( expression=a %% t(a), name="A" ),
mxAlgebra( expression=c %% t(c), name="C" ),
mxAlgebra( expression=e %% t(e), name="E" ),
# Algebra to compute total variances and standard deviations (diagonal only)
mxAlgebra( expression=A+C+E, name="V" ),
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( expression=solve(sqrt(IV)), name="iSD"),
# Confidence intervals portion for covariance matrices
mxAlgebra( expression=A/V,name="stndVCA"),
mxAlgebra( expression=C/V,name="stndVCC"),
mxAlgebra( expression=E/V,name="stndVCE"),
mxCI(c("stndVCA","stndVCC","stndVCE")),
# Confidence intervals for Phenotypic correlation pieces
mxAlgebra((solve(sqrt(IV)) %% V %% solve(sqrt(IV)))[1,2],name="Vcorr"),
mxAlgebra((solve(sqrt(IA)) %% A %% solve(sqrt(IA)))[1,2],name="Acorr"),
mxAlgebra((solve(sqrt(IC)) %% C %% solve(sqrt(IC)))[1,2],name="Ccorr"),
mxAlgebra((solve(sqrt(IE)) %% E %% solve(sqrt(IE)))[1,2],name="Ecorr"),
mxCI(c("Vcorr","Acorr","Ccorr","Ecorr")),
## 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=nv, free=c(TRUE,FALSE), values=meanSVnv, labels=mean, name="Mean" ),
mxAlgebra( expression= cbind(Mean,Mean), name="expMean"),
mxMatrix(type="Full", nrow=2, ncol=nv, free=c(TRUE,TRUE), values=c(1,1.8),lbound=-4, ubound=6, name="expThre", labels=c("th1","th2"),dimnames=list(c('th1','th2'),selVars[c(2,4)])),
# Algebra for expected variance/covariance matrix in MZ
mxAlgebra( expression= rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" ),
# Algebra for expected variance/covariance matrix in DZ
mxAlgebra( expression= 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, thresholds="ACE.expThre", threshnames=selVars[c(2,4)] )
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThre", threshnames=selVars[c(2,4)] )
),
mxAlgebra( expression=MZ.objective + DZ.objective, name="neg2sumll" ),
mxAlgebraObjective("neg2sumll")
)
bivACEFit <- mxRun(bivACEModel,intervals=TRUE)
bivACESumm <- summary(bivACEFit)
bivACESumm