I was wondering if there was a reason I can't seem to get confidence intervals for an ordinal ACE model with 2 thresholds.
I made sure to specify the mxCI for my standardized variance components call in my model and to put the "intervals=T" in my mxRun line. In my model fit $output$confidenceIntervals exists, but does not have any values.
Is this normal? If not, is there a way to create confidence intervals in the ordinal case? If it's just a silly code error, my code is below.
Any help would be appreciated!
univACEOrdModel <- mxModel("univACEOrd",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="a" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="c11", name="c" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", 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="sd"),
mxAlgebra( expression=cbind(A/VP,C/VP,E/VP),name="stndVCs"),
# Calculate 95% CIs here
mxAlgebra(A+C+E,name="VP"),
Yes, it's repetitive, but I was desperate and the above was exactly how it ran in my continuous univariate ACE model
mxCI(c("stndVCs")),
# Constraint on variance of ordinal variables
mxConstraint(V == I, name="Var1"),
# Matrix & Algebra for expected means vector
mxMatrix( type="Zero", nrow=1, ncol=nv, name="M" ),
mxAlgebra( expression= cbind(M,M), name="expMean" ),
mxMatrix( type="Full", nrow=2, ncol=nv, free=TRUE, values=c(0.8,1.2), label=c("thre1","thre2"), name="T" ),
mxAlgebra( expression= cbind(T,T), dimnames=list(c('th1','th2'),selVars), name="expThre" ),
# 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, note use of 0.5, converted to 1*1 matrix
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" )
),
mxModel("DZ",
mxData( observed=dzData, type="raw" ),
mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars, thresholds="ACE.expThre" )
),
mxAlgebra( expression=MZ.objective + DZ.objective, name="min2sumll" ),
mxAlgebraObjective("min2sumll")
)
univACEOrdFit <- mxRun(univACEOrdModel,intervals=T)
It was only recently in the OpenMx 1.1 beta release that we started support confidence intervals specified in the submodels. You can try downloading the beta and see if the script works. Another alternative would be to specify
mxCI("ACE.stndVCs")
in the container model.Thanks so much for the quick reply.
That's the problem with being lazy and adapting other people's scripts. I removed the unnecessary submodel and everything is working perfectly now.
Kelly
It is better to spell out TRUE and FALSE as one day you will use a variable T (or F) and wonder why your previously working example stops working
Along the same lines, you should avoid naming you MxMatrix or MxAlgebra object "T" or "F" or "c". The model will execute, but should you try to use mxEval() to examine components of your model, then the model names will override the R function or variable names. For example,
will no longer work because mxEval() beleives that 'c' is a MxMatrix object, not the c() function in R.